Fitting Data and Making FOMs

Here I will go through an example of how to fit data to a transfer function. Later this TF will be used to make a full system model.

First imports:

%load_ext autoreload
%autoreload 2
import numpy as np
import matplotlib.pyplot as plt
import control
import scipy.signal as signal
from scipy import interpolate
import scipy.linalg
from scipy.stats import gmean
import gwinc
import yaml
from wield.iirrational.v2 import data2filter
from wield.control.AAA import tfAAA
from wield.utilities.file_io import load, save
from wield.control.ss_bare.ss import BareStateSpace
from scipy.signal import cheby2
from wield.control import SISO
from buzz import ssutil
from physunits import m, s, kg, Hz
from scipy.constants import c, G

/home/mcculler/local/projects_sync/wield-project-lab/wield-iirrational/src/wield/iirrational/_version.py:35: UserWarning: git base version 0.1 is different than the stored version 2.9.5
  warnings.warn(

Next we will define a function to grab the BNS waveform

# from https://github.com/wield/wield-control/blob/main/src/wield/control/AAA/test/test_AAA_present.py
def gen_waveform(**params):
    """Generate frequency-domain inspiral waveform

    Returns a tuple of (freq, h_plus^tilde, h_cross^tilde).

    The waveform is generated with the lalsimulation
    SimInspiralChooseFDWaveform() function.  Keyword arguments are
    used to update the default waveform parameters (see DEFAULT_PARAMS
    macro).  The mass parameters ('m1' and 'm2') should be specified
    in solar masses and the 'distance' parameter should be specified
    in parsecs**.  Waveform approximants may be given as string names
    (see `lalsimulation` documentation for more info).

    For example, to generate a 20/20 Msolar BBH waveform:

    >>> hp,hc = waveform.gen_waveform('m1'=20, 'm2'=20)

    **NOTE: The requirement that masses be specified in solar masses
    and distances in parsecs is different than that of the underlying
    lalsimulation method which expects mass and distance parameters to
    be in SI units.

    """
    import lalsimulation
    from inspiral_range import waveform
    from inspiral_range import const

    iparams = dict(waveform.DEFAULT_PARAMS)
    iparams.update(**params)

    # convert to SI units
    iparams["distance"] *= const.PC_SI
    iparams["m1"] *= const.MSUN_SI
    iparams["m2"] *= const.MSUN_SI
    iparams["approximant"] = lalsimulation.SimInspiralGetApproximantFromString(
        iparams["approximant"]
    )

    m = iparams["m1"] + iparams["m2"]

    # calculate delta F based on frequency of inner-most stable
    # circular orbit ("fisco")
    fisco = (const.c ** 3) / (const.G * (6 ** 1.5) * 2 * np.pi * m)
    df = 2 ** (np.max([np.floor(np.log(fisco / 4096) / np.log(2)), -6]))

    # FIXME: are these limits reasonable?
    if iparams["deltaF"] is None:
        iparams["deltaF"] = df
    # iparams['f_min'] = 0.1
    # iparams['f_max'] = 10000

    hp, hc = lalsimulation.SimInspiralChooseFDWaveform(**iparams)
    
    freq = hp.f0 + np.arange(len(hp.data.data)) * hp.deltaF

    print("f0", hp.f0)
    select = abs(hp.data.data) > 0

    return freq[select], hp.data.data[select], hc.data.data[select]

Next we will define the parameters of our BNS model and plot its PSD

params = dict()
params["m1"] = 1.4
params["m2"] = 1.4
# rest of params from first BNS detection see https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.119.161101
params["f_min"] = 0.01
params["f_max"] = 1000000
params["deltaF"] = 0.01
params["distance"] = 1e6
F_Hz, hp, hc = gen_waveform(**params) # assuming that this is giving an ASD I.M. Jan 2024

BNS_intrp = interpolate.interp1d(F_Hz, abs(hp)**2) ## we add interpolation to put all of our data with the same FQ spacing

print("F_Hz max", max(F_Hz))

plt.loglog(F_Hz, abs(hp)**2, label = "hp")
plt.loglog(F_Hz, abs(hc)**2, linestyle = ":", label = "hc")
plt.title(str(params["m1"]) +"/"+str(params["m1"])+ "Solarmass BNS Spec")
plt.xlabel("Frequency [Hz]")
plt.ylabel("Power spectral density [strain^2/Hz^2]")
plt.legend()
plt.show()

/opt/conda/user_conda/mcculler/py312/lib/python3.12/site-packages/lalsimulation/lalsimulation.py:8: UserWarning: Wswiglal-redir-stdio:

SWIGLAL standard output/error redirection is enabled in IPython.
This may lead to performance penalties. To disable locally, use:

with lal.no_swig_redirect_standard_output_error():
    ...

To disable globally, use:

lal.swig_redirect_standard_output_error(False)

Note however that this will likely lead to error messages from
LAL functions being either misdirected or lost when called from
Jupyter notebooks.

To suppress this warning, use:

import warnings
warnings.filterwarnings("ignore", "Wswiglal-redir-stdio")
import lal

  import lal

f0 0.0
F_Hz max 14501.800000000001

# # compare the waveform with part of equation 3 from https://arxiv.org/abs/1709.08079
m1 = (params["m1"] * 1.988416e30) * kg
m2 = (params["m2"] * 1.988416e30) * kg
chirp_mass = (m1*m2)**(3/5)/(m1+m2)**(1/5)
print('chirpmass: \t', chirp_mass, ' kg\n\t\t', chirp_mass/1.988416e30, ' solar masses')


# N = kg* m/(s**2)
#c = c #* m/s
# G = G * N * m**2 / (kg**2)
# F_Hz = F_Hz * Hz
# eq3 = ((5/24)**0.5 * c * ( G * chirp_mass * (c**3))**(5/6)* 1/(np.pi**(-2/3)))**2 * F_Hz**(-7/3)
# m2Mpc = 1/3.086e22
# eq3_Mpc = eq3

# print("eq3[0]", eq3[0])
# print("eq3_Mpc[0]", eq3_Mpc[0])

# print("BNS hp[0]", abs(hp[0])**2)

# plt.loglog(F_Hz, eq3_Mpc, label="eq3")
# plt.loglog(F_Hz, abs(hp)**2, label="BNS Strain from LAL Simulation")
# plt.title(str(params["m1"]) +"/"+str(params["m1"])+ "Solarmass BNS Spec")
# plt.xlabel("Frequency [Hz]")
# plt.ylabel("Power spectral density [strain^2/Hz^2]")
# plt.legend()
# plt.show()

chirpmass: 	 2.423423336413827e+30  kg
		 1.2187707886145691  solar masses

Getting the LIGO Sensitivity Curve

Now I must multiply this by the GWINC curve to get our final FOM

F_bug = np.logspace(np.log10(min(F_Hz)), np.log10(max(F_Hz)-1), 10000)
budget = gwinc.load_budget('aLIGO', freq=F_bug)
trace = budget.run(freq=F_bug)

plt.title("ALIGO BUDGET")
plt.loglog(F_bug, trace.psd)
plt.xlabel("Frequency [Hz]")
plt.ylabel("Power Spectral Density [strain/hz]")
# plt.xlim([6, max(F_Hz)])
# plt.ylim([1e-20, 1e-17])

Text(0, 0.5, 'Power Spectral Density [strain/hz]')

prefactor = (2*(5/96)**0.5 * c * ( G * chirp_mass * (c**3))**(5/6) * 1/(np.pi**(-2/3))* 2 / 2.26)**2 
I7 = np.trapz(F_bug**(-7/3)/trace.psd, F_bug)
eq3 = prefactor * I7**0.5
print("eq3: ", eq3/3.086e22, ' Mpc')

eq3:  4.867082672729032e+92  Mpc

/tmp/ipykernel_1398697/2399436929.py:2: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`.
  I7 = np.trapz(F_bug**(-7/3)/trace.psd, F_bug)

from inspiral_range import inspiral_range as ir
from inspiral_range import waveform
from inspiral_range import const

DETECTION_SNR = 8.0
def ian_sensemon_range(freq, m1=1.4, m2=1.4, horizon=False, integrate=False, detection_snr=DETECTION_SNR):
    """Detector inspiral range from closed form expression

    Masses `m1` and `m2` should be specified in solar masses (default:
    m1=m2=1.4).  If the `horizon` keyword is specified the "horizon"
    range will be returned, which differs from the angle-averaged
    range by ~2.26.

    @returns distance in Mpc as a float

    """
    if horizon:
        theta = 4
    else:
        theta = 1.77
    theta /= 1e6 * const.PC_SI
    M_chirp = waveform.M_chirp(m1, m2) * const.MSUN_SI

    i73 = freq ** (-7/3)
    val = theta / detection_snr \
        * waveform.habs_nsp_prefactor(M_chirp) \
        * np.sqrt(i73) / 2

    return val

irdata = ian_sensemon_range(F_bug, integrate=False, horizon=True)**2

range_ir = ir.sensemon_range(F_bug, psd=trace.psd, m1=1.4, m2=1.4, horizon=False, integrate=True, detection_snr=DETECTION_SNR)
print("range: ", range_ir)

#np.savetxt('FBNS_Current_range.txt', [range])

BNS_PSD = BNS_intrp(F_bug)
IR_interp = interpolate.interp1d(F_bug, irdata)

midpoint = len(F_bug)//2
print('relative difference: ', BNS_PSD[midpoint]/irdata[midpoint])


plt.loglog(F_bug, BNS_PSD)
plt.loglog(F_bug, irdata)
plt.scatter(F_bug[midpoint], BNS_PSD[midpoint])
plt.scatter(F_bug[midpoint], irdata[midpoint])
plt.xlabel("Frequency [Hz]")
plt.ylabel("Power spectral density");

range:  194.97198453267936
relative difference:  15.834426264720078

Choose which inspiral waveform PSD to use

#BNS_PSD = BNS_intrp(F_bug) # us the custom inspiral equatin from LAL
BNS_PSD = IR_interp(F_bug) # use the inspiral range aproximation and calibration 

Now I will put the BNS waveform in the Fq space of the GWINC curve

print("BNS_PSD", BNS_PSD)
print("irdata", irdata)
plt.loglog(F_bug, BNS_PSD)
plt.xlabel("Frequency [Hz]")
plt.ylabel("Power spectral density")

BNS_PSD [2.45847163e-35 2.45034592e-35 2.44224707e-35 ... 1.03982333e-49
 1.03638652e-49 1.03296107e-49]
irdata [2.45847163e-35 2.45034592e-35 2.44224707e-35 ... 1.03982333e-49
 1.03638652e-49 1.03296107e-49]

Text(0, 0.5, 'Power spectral density')

Lastly for this step we will combine the twu curves to make a FOM. we divide the BNS PSD by the LIGO Noise spectrum without the controls noise. From the paper: $$\mathrm{SNR}^2 = 4\int_{0}^{\infty} \frac{\left |h_\mathrm{gw}(f)\right |^2}{S_h(f)} df$$

where $S_h(f)$ is the noise spectrum of the detector. This can be written as:

$$ S_h(f) = S_{\mathrm{det}}(f) + |C(f)|^2 S_{\mathrm{p}}(f)$$

where $S_{\mathrm{det}}(f)$ is the detector noise spectrum without controls noise, $S_{\mathrm{p}}(f)$ is the control noise spectrum, and $C(f)$ is the coupling function from the controls noise in the control loop to DARM. By using the taylor series of A/(B+x) to linearize the integral. i.e.:

$$\frac{A}{B+x}=\frac{A}{B} - \frac{Ax}{B^2} + \frac{Ax^2}{B^3} - \frac{Ax^3}{B^4} + \frac{Ax^4}{B^5} - \frac{Ax^5}{B^6} + O(x^6)$$

taking the first two terms of the taylor series and applying it to the SNR equation we get:

$$\mathrm{SNR}^2 = 4\int_{0}^{\infty} \frac{\left | h_{\mathrm{signal}}(f)\right |^2}{S_{\mathrm{det}}(f)} \ df - 4\int_{0}^{\infty} \frac{\left |C\right |^2\left | h_{\mathrm{signal}}(f)\right |^2 }{S_{\mathrm{det}}(f)^2} S_{\mathrm{p}}(f) \ df$$

To get range we need to use the $\mathrm{SNR}$ not the $\mathrm{SNR}^2$ so we take the square root of the above equation which has a different linearization. as $A \rightarrow \infty$

$$\sqrt{A+B} = \sqrt{A} \sqrt{1+\frac{B}{A}} = \sqrt{A} (1+\frac{B}{2 \sqrt{A}})$$

Div_PSD = BNS_PSD / (trace.psd)**2 # divide the BNS PSD by the aLIGO noise psd squared. this will give a PSD of LIGO's sensitivity to a BNS signal
#Div_PSD = Div_PSD * 4000**2 # to put it in units of displacement
F_Div = F_bug
plt.loglog(F_Div, Div_PSD)
plt.xlabel("Frequency [Hz]")
plt.ylabel("Power spectral density")
plt.title("PSD of FOM for BNS Signal")

Text(0.5, 1.0, 'PSD of FOM for BNS Signal')

Fitting the Data to a Transfer Function

First we will Define a function to reduce the number of points to fit to

# edited from https://git.ligo.org/wield/wield-ligo-mcculler/-/blob/main/src/wield/LIGO/mcculler/filter_cavity/test_FC1_SUSPOINT_fits.py
def reduce(F_Hz, PSD):
    F_groups = np.logspace(-2, 4, 50) # originally 200 length
    F_groups = np.geomspace(min(F_Hz), max(F_Hz), 50) # originally 200 length
    
    idx_groups = np.searchsorted(F_Hz, F_groups)
    idx_pairs = list(zip(idx_groups[:-1], idx_groups[1:]))

    lPSDs_min = []
    lPSDs_max = []
    lPSDs_med = []
    lPSDs_mean = []
    lF_Hz = []
    for idx1, idx2 in idx_pairs:
        if idx1 == idx2:
            continue
        lPSDs_max.append(
            np.nanmax(PSD[idx1:idx2])
        )
        lPSDs_med.append(
            np.nanmedian(PSD[idx1:idx2])
        )
        lPSDs_min.append(
            np.nanmin(PSD[idx1:idx2])
        )
        lPSDs_mean.append(
            np.nanmean(PSD[idx1:idx2])
        )
        lF_Hz.append(np.mean(F_Hz[idx1:idx2]))

    return np.asarray(lF_Hz), np.asarray(lPSDs_min), np.asarray(lPSDs_med), np.asarray(lPSDs_max), np.asarray(lPSDs_mean)

Now we will reduce the LIGO sensitivity curve and its reduced form

F_FOM_dwn, PSD_min, PSD_med, PSD_max, PSD_mean = reduce(F_bug, Div_PSD)
FOM_PSD = PSD_mean
# The current form of FOM_PSD is the ASD
FOM_ASD = FOM_PSD**0.5

plt.loglog(F_bug, Div_PSD**0.5, label = 'FOM ASD')
plt.loglog(F_FOM_dwn, FOM_ASD, label = 'Downsampled FOM ASD')
plt.legend()

<matplotlib.legend.Legend at 0x7fad19ccbd40>

Now we will define some functions to work with these transfer functions

def FBNSsimpSS(gain=1, return_name=False, lo_ord=4):
    if return_name:
        name = "FOM BNS Simple"
        return name
    F_Fq_lo = 10 * 2 * np.pi
    F_Fq_hi = 120 * 2 * np.pi
    hp_order = lo_ord
    lp_order = 2
    F_p = []
    F_z = []
    for ord in range(hp_order):
        F_p.append(-F_Fq_lo + 0*1j*F_Fq_lo)
        F_p.append(-F_Fq_lo - 0*1j*F_Fq_lo)
        F_z.append(0)
        F_z.append(0)

    for ord in range(lp_order):
        F_p.append(-F_Fq_hi)
    F_k = 1
    c_zpk = cheby2(8, 170, 1.25*2*np.pi, btype='high', analog=True, output='zpk')
    F_z.extend(c_zpk[0])
    F_p.extend(c_zpk[1])
    F_k *= c_zpk[2]

    ian_lfq = -0.7
    ian_lfq2 = 2
    ian_lfq3 = 0.5
    F_z.append(0)
    F_p.append(ian_lfq)
    F_z.append(0)
    F_p.append(ian_lfq)

    F_z.append(ian_lfq3)
    F_z.append(ian_lfq3)
    F_p.append(ian_lfq2)
    F_p.append(ian_lfq2)
    
    F_p = F_p
    F_z = F_z

    #F_mod = ws_tf2ss(F_z, F_p, F_k)
    #F_mod = ssutil.normalize_gain(F_mod) * gain 
    #F_iod = {"FBNS.in": 0, "FBNS.out": 0}
    F_mod = SISO.zpk(F_z, F_p, F_k, angular=True, fiducial_rtol=1e-5, fiducial_atol=1e-10)
    F_mod = F_mod * F_mod
    F_mod = F_mod.asSS
    F_mod = F_mod / abs(F_mod.asSS.Linf_norm()[0]) 


    F_mod = F_mod * gain
    F_mod.balance_and_truncate()

    #F = ssutil.wieldSS(F_mod, F_iod)
    
    return F_mod.mimo("FBNS.in", "FBNS.out")

plt.loglog(F_bug, Div_PSD**0.5, label = 'FOM ASD')
plt.loglog(F_FOM_dwn, FOM_ASD, label = 'Downsampled FOM ASD')
mag, phase, omega = control.bode(FBNSsimpSS(gain=4 * 100**2).mod, Hz=False, plot=False, omega=F_FOM_dwn*2*np.pi, label = "BNS Simple")
#plt.loglog(omega/(2*np.pi), mag, label="BNS Simple")
plt.legend()

/home/mcculler/local/projects_sync/wield-project-lab/wield-control/src/wield/control/SISO/zpk.py:164: NumericalWarning: StateSpace is large (>50 states), using reduced response fiducial auditing heuristics. TODO to make this smarter
  warnings.warn(f"StateSpace is large (>{self.N_MAX_FID} states), using reduced response fiducial auditing heuristics. TODO to make this smarter", util.NumericalWarning)
/opt/conda/user_conda/mcculler/py312/lib/python3.12/site-packages/control/freqplot.py:435: FutureWarning: bode_plot() return value of mag, phase, omega is deprecated; use frequency_response()
  warnings.warn(

<matplotlib.legend.Legend at 0x7fad1bb86ba0>

from scipy.signal import cheby2
from wield.control import SISO

def FBNSsimpSS(gain=1, return_name=False, lo_ord=4, return_scale=False):
    if return_name:
        name = "FOM BNS Simple"
        return name
    
    F_Fq_lo = 40 * 2 * np.pi
    # aggressive version
    # F_Fq_lo = 30 * 2 * np.pi
    F_Fq_hi = 300 * 2 * np.pi
    hp_order = 1 #lo_ord
    lp_order = 2
    F_p = []
    F_z = []
    for ord in range(hp_order):
        F_p.append(-F_Fq_lo + 0*1j*F_Fq_lo)
        F_p.append(-F_Fq_lo - 0*1j*F_Fq_lo)
        F_z.append(0)
        F_z.append(0)
    for ord in range(lp_order):
        F_p.append(-F_Fq_hi)
    F_k = gain * 1.5e19
    c_zpk = cheby2(8, 160, 4*2*np.pi, btype='high', analog=True, output='zpk')
    # aggressive version
    # c_zpk = cheby2(8, 120, 4*2*np.pi, btype='high', analog=True, output='zpk')
    F_z.extend(c_zpk[0])
    F_p.extend(c_zpk[1])
    F_k *= c_zpk[2]
    
    F = SISO.zpk(F_z, F_p, F_k)
    return F*F
sys_FBNS_simp = FBNSsimpSS()

zpk_dict = {
            'z' : [str(zero) for zero in np.asarray(tuple(sys_FBNS_simp.z)).tolist()],
            'p' : [str(pole) for pole in np.asarray(tuple(sys_FBNS_simp.p)).tolist()],
            'k' : float(sys_FBNS_simp.k)
            }

save('FBNSsimp_FOM.yml', zpk_dict)

with open('BNS_FOM_handfit.yml', 'r') as file:
    filt_dict = yaml.safe_load(file)
p_hand = np.array(filt_dict['p'], dtype=np.complex128)
z_hand = np.array(filt_dict['z'], dtype=np.complex128)
k_hand = np.float64(filt_dict['k'])

hand_zpk = SISO.zpk(z_hand, p_hand, k_hand)
hand_xfr = hand_zpk.fresponse(f=F_FOM_dwn)
BNSsimp_xfr = sys_FBNS_simp.fresponse(f=F_FOM_dwn)

plt.loglog(F_bug, Div_PSD**0.5, label = 'FOM ASD')
plt.loglog(F_FOM_dwn, FOM_ASD, label = 'Downsampled FOM ASD')
plt.loglog(F_FOM_dwn, hand_xfr.mag, label = 'Handfit ASD')
plt.loglog(F_FOM_dwn, BNSsimp_xfr.mag, label = 'BNSSimp ASD')
plt.legend()

/home/mcculler/local/projects_sync/wield-project-lab/wield-control/src/wield/control/SISO/zpk.py:164: NumericalWarning: StateSpace is large (>50 states), using reduced response fiducial auditing heuristics. TODO to make this smarter
  warnings.warn(f"StateSpace is large (>{self.N_MAX_FID} states), using reduced response fiducial auditing heuristics. TODO to make this smarter", util.NumericalWarning)

<matplotlib.legend.Legend at 0x7fad07144c20>

# Divide the FOM ASD by the hand fit to get the reduced ASD
reduced_ASD = FOM_ASD/hand_xfr.mag

plt.loglog(F_FOM_dwn, reduced_ASD, label = 'reduced ASD')
plt.legend()

<matplotlib.legend.Legend at 0x7fad19c22480>

Now we form a guess of the TF

#results = tfAAA(F_Hz=F_FOM_dwn, xfer=FOM_ASD)
results = tfAAA(F_Hz=F_FOM_dwn, xfer=reduced_ASD)

Now we play with the poles and zeros to get a good fit

#print("poles", results.poles)
#print("zeros", results.zeros)
#print("gain", results.gain)

all_z = results.zeros
all_p = results.poles
k = results.gain

for i in range(0, len(all_z)):
    if all_z[i] >= 0:
        all_z[i] = -all_z[i]
        
for i in range(0, len(all_p)):
    if all_p[i] >= 0:
        all_p[i] = -all_p[i]

AAA_zpk = SISO.zpk(all_z, all_p, k, angular=False, fiducial_rtol=1e-5, fiducial_atol=1e-10)
recalibrated_zpk =  hand_zpk * AAA_zpk 
# recalibrated_zpk =  AAA_zpk 
print(recalibrated_zpk.zeros)

2π·[-0.3958406743523140, -0.3958406743523140, -0.3958406743523140,
 -0.3958406743523140, -0.3958406743523140, -0.3958406743523140,
 -0.3958406743523140, -0.3958406743523140, -0.3958406743523140,
 -0.3958406743523140, -0.3958406743523140, -0.3958406743523140,
 -0.3958406743523140, -0.3958406743523140, -0.3958406743523140,
 -0.3958406743523140, -0.3958406743523140, -0.3958406743523140,
 -5.305505757497158e-10, -0.005431630862136328, -0.005431630295139218,
 -0.04216327861663196, -0.04216327868321049, -3.049643520855844,
 -3.049643519812004, -2707.242418662691, -2707.242418662458,
 -2.957619742415325e-12 ± 92.58406532411048j,
 -14.98934980857282     ± 21.76848865621670j,
 -14.98934980875151     ± 21.76848865620463j,
 -2.069960897503006e-10 ± 10.00487816878042j,
 -1.178440698489668e-11 ± 2.829913488939144j,
 -1.029654952340069e-11 ± 1.332579977525743j,
 -0.01709138866381974   ± 0.2654485123713641j,
 -0.01709138868276358   ± 0.2654485123375496j,
 -6.152212595327936e-12 ± 0.3571372397086687j,
 -0.4720203451684382    ± 2.629597995989603j,
 -0.4720203452849284    ± 2.629597996216710j]

plt.loglog(F_FOM_dwn, reduced_ASD, label = 'Reduced ASD')
plt.loglog(F_bug, AAA_zpk.fresponse(f=F_bug).mag, label = 'AAA Fit', linestyle = '--')
plt.legend()

<matplotlib.legend.Legend at 0x7fad196fca70>

# Get the ZPK and the minreal ZPK for all poles and zeros
B_res_all = control.zpk(all_z, all_p, k)
B_res_all_minreal = control.minreal(B_res_all)
minreal_all_z, minreal_all_p, minreal_all_k = signal.tf2zpk(B_res_all_minreal.num[0][0], B_res_all_minreal.den[0][0])

# Get only the real poles and zeros
z = [lz for lz in all_z if lz.real < 0]
p = [lp for lp in all_p if lp.real < 0]

# Get the ZPK for the real poles and zeros and the minreal ZPK for only the real poles and zeros
B_res = control.zpk(z, p, k)
B_res_minreal = control.minreal(B_res)
minreal_z, minreal_p, minreal_k = signal.tf2zpk(B_res_minreal.num[0][0], B_res_minreal.den[0][0])

0 states have been removed from the model
0 states have been removed from the model

/opt/conda/user_conda/mcculler/py312/lib/python3.12/site-packages/control/xferfcn.py:1087: ComplexWarning: Casting complex values to real discards the imaginary part
  den[j, :maxindex+1] = poly(poles[j])
/opt/conda/user_conda/mcculler/py312/lib/python3.12/site-packages/control/xferfcn.py:1117: ComplexWarning: Casting complex values to real discards the imaginary part
  num[i, j, maxindex+1-len(numpoly):maxindex+1] = numpoly

Now we will plot our guess and our othe tfs

_, TF_all = scipy.signal.freqs_zpk(all_z, all_p, k, worN = F_Hz)
_, TF_minreal = scipy.signal.freqs_zpk(minreal_all_z, minreal_all_p, minreal_all_k, worN = F_Hz)

recalibrated_zpk_xfr = recalibrated_zpk.fresponse(f=F_Hz)

plt.loglog(F_bug, Div_PSD**0.5, label = 'Full FOM ASD')
plt.loglog(F_FOM_dwn, hand_xfr.mag, label = 'Handfit ASD')
plt.loglog(F_Hz, recalibrated_zpk_xfr.mag, label = 'Recalibrated AAA')
plt.loglog(F_FOM_dwn, FOM_ASD, label = 'Downsampled FOM ASD')
#plt.loglog(F_Hz, abs(TF_all), label = 'ALL ZPK Fit Bary. (order {})'.format(results.order))
#plt.loglog(F_Hz, abs(TF_minreal), label = 'Minreal ALL ZPK Fit Bary.', linestyle = '--')
plt.legend()

<matplotlib.legend.Legend at 0x7fad1898d460>

Now we fit for the poles and zeros

snr = np.ones_like(F_FOM_dwn)
snr[:-1] = (F_FOM_dwn[1:] - F_FOM_dwn[:-1])**-0.001
from wield.iirrational import fitters_ZPK

all_z = (np.array(recalibrated_zpk.zeros))
all_p = (np.array(recalibrated_zpk.poles))
k = float(recalibrated_zpk.k)

z_iir = tuple(np.array(all_z) / (np.pi * 2))
p_iir = tuple(np.array(all_p) / (np.pi * 2))
print(z_iir)
print(p_iir)

iir_results = data2filter(
    F_Hz=F_FOM_dwn,
    xfer=FOM_ASD,
    #xfer=reduced_ASD,
    mode='reduce',
    zeros=z_iir, 
    poles=p_iir,
    gain=k * (2*np.pi)**(len(z_iir) - len(p_iir)),
    SNR_phase_rel=0,
    SNR=snr,
    #relative_degree=-4,
    # resavg_RthreshOrdDn=1.01,
    baseline_only=True,
    #coding_map=fitters_ZPK.codings_s.coding_maps.RI
    # trust_SNR = True,
)
print("C", iir_results.fit_aid._fitters[0].fitter.poles.fullplane)
print(tuple(np.array(all_p) / (np.pi * 2)))

(np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-0.395840674352314), np.float64(-5.305505757497158e-10), np.float64(-0.005431630862136328), np.float64(-0.005431630295139218), np.float64(-0.04216327861663196), np.float64(-0.042163278683210494), np.float64(-3.0496435208558443), np.float64(-3.049643519812004), np.float64(-2707.2424186626913), np.float64(-2707.242418662458), np.complex128(-2.957619742415325e-12+92.58406532411048j), np.complex128(-2.957619742415325e-12-92.58406532411048j), np.complex128(-14.989349808572818+21.768488656216697j), np.complex128(-14.989349808572818-21.768488656216697j), np.complex128(-14.989349808751511+21.768488656204628j), np.complex128(-14.989349808751511-21.768488656204628j), np.complex128(-2.0699608975030063e-10+10.004878168780424j), np.complex128(-2.0699608975030063e-10-10.004878168780424j), np.complex128(-1.1784406984896683e-11+2.8299134889391437j), np.complex128(-1.1784406984896683e-11-2.8299134889391437j), np.complex128(-1.0296549523400693e-11+1.3325799775257434j), np.complex128(-1.0296549523400693e-11-1.3325799775257434j), np.complex128(-0.017091388663819738+0.2654485123713641j), np.complex128(-0.017091388663819738-0.2654485123713641j), np.complex128(-0.017091388682763577+0.2654485123375496j), np.complex128(-0.017091388682763577-0.2654485123375496j), np.complex128(-6.152212595327936e-12+0.3571372397086687j), np.complex128(-6.152212595327936e-12-0.3571372397086687j), np.complex128(-0.4720203451684382+2.629597995989603j), np.complex128(-0.4720203451684382-2.629597995989603j), np.complex128(-0.47202034528492837+2.6295979962167104j), np.complex128(-0.47202034528492837-2.6295979962167104j))
(np.float64(-200.0), np.float64(-200.0), np.float64(-1999.9999999999998), np.float64(-1.1743521290869434e-10), np.float64(-0.011159044867910854), np.float64(-0.011159044751801689), np.float64(-1302.7841836435434), np.float64(-1302.7841836432744), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-1.276683601197908e-11+91.63193908076074j), np.complex128(-1.276683601197908e-11-91.63193908076074j), np.complex128(-18.709468807770158+25.224614971813303j), np.complex128(-18.709468807770158-25.224614971813303j), np.complex128(-18.709468807798835+25.224614971785446j), np.complex128(-18.709468807798835-25.224614971785446j), np.complex128(-3.5224484988923e-11+10.490717255985867j), np.complex128(-3.5224484988923e-11-10.490717255985867j), np.complex128(-7.755934678337948e-12+2.81092371100195j), np.complex128(-7.755934678337948e-12-2.81092371100195j), np.complex128(-3.1553850493935663e-12+1.2246948964676j), np.complex128(-3.1553850493935663e-12-1.2246948964676j), np.complex128(-8.200117750311727e-13+0.29012748368774854j), np.complex128(-8.200117750311727e-13-0.29012748368774854j), np.complex128(-0.03924083547004946+0.07982846579002265j), np.complex128(-0.03924083547004946-0.07982846579002265j), np.complex128(-0.0392408354694464+0.07982846579497013j), np.complex128(-0.0392408354694464-0.07982846579497013j), np.complex128(-4.26914787276891e-11+0.5783185011865986j), np.complex128(-4.26914787276891e-11-0.5783185011865986j), np.complex128(-4.995095343193837e-11+0.6196917699694076j), np.complex128(-4.995095343193837e-11-0.6196917699694076j), np.complex128(-0.7658207288552369+6.119579611317334j), np.complex128(-0.7658207288552369-6.119579611317334j), np.complex128(-0.7658207289076394+6.119579611415145j), np.complex128(-0.7658207289076394-6.119579611415145j))
TEE_LOGFILE None
A [-2.00000000e+02+0.00000000e+00j -2.00000000e+02+0.00000000e+00j
 -2.00000000e+03+0.00000000e+00j -1.17435213e-10+0.00000000e+00j
 -1.11590449e-02+0.00000000e+00j -1.11590448e-02+0.00000000e+00j
 -1.30278418e+03+0.00000000e+00j -1.30278418e+03+0.00000000e+00j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -1.27668360e-11+9.16319391e+01j -1.27668360e-11-9.16319391e+01j
 -1.87094688e+01+2.52246150e+01j -1.87094688e+01-2.52246150e+01j
 -1.87094688e+01+2.52246150e+01j -1.87094688e+01-2.52246150e+01j
 -3.52244850e-11+1.04907173e+01j -3.52244850e-11-1.04907173e+01j
 -7.75593468e-12+2.81092371e+00j -7.75593468e-12-2.81092371e+00j
 -3.15538505e-12+1.22469490e+00j -3.15538505e-12-1.22469490e+00j
 -8.20011775e-13+2.90127484e-01j -8.20011775e-13-2.90127484e-01j
 -3.92408355e-02+7.98284658e-02j -3.92408355e-02-7.98284658e-02j
 -3.92408355e-02+7.98284658e-02j -3.92408355e-02-7.98284658e-02j
 -4.26914787e-11+5.78318501e-01j -4.26914787e-11-5.78318501e-01j
 -4.99509534e-11+6.19691770e-01j -4.99509534e-11-6.19691770e-01j
 -7.65820729e-01+6.11957961e+00j -7.65820729e-01-6.11957961e+00j
 -7.65820729e-01+6.11957961e+00j -7.65820729e-01-6.11957961e+00j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j]
B [-2.00000000e+03+0.00000000e+00j -3.90376294e-03+0.00000000e+00j
 -1.30278418e+03+2.64289979e-05j -1.30278418e+03-2.64289979e-05j
 -1.11590448e-02+1.16415322e-10j -1.11590448e-02-1.16415322e-10j
 -2.00000000e+02+3.81469727e-06j -2.00000000e+02-3.81469727e-06j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -7.65820729e-01+6.11957961e+00j -7.65820729e-01-6.11957961e+00j
 -7.65820729e-01+6.11957961e+00j -7.65820729e-01-6.11957961e+00j
 -4.64082010e-04+6.19691770e-01j -4.64082010e-04-6.19691770e-01j
 -4.33097913e-04+5.78318501e-01j -4.33097913e-04-5.78318501e-01j
 -3.92408355e-02+7.98284658e-02j -3.92408355e-02-7.98284658e-02j
 -3.92408355e-02+7.98284658e-02j -3.92408355e-02-7.98284658e-02j
 -9.04325228e-03+2.90127484e-01j -9.04325228e-03-2.90127484e-01j
 -8.84127110e-02+1.22469490e+00j -8.84127110e-02-1.22469490e+00j
 -2.88099878e-03+2.81092371e+00j -2.88099878e-03-2.81092371e+00j
 -3.07907738e+00+1.04907173e+01j -3.07907738e+00-1.04907173e+01j
 -1.87094688e+01+2.52246150e+01j -1.87094688e+01-2.52246150e+01j
 -1.87094688e+01+2.52246150e+01j -1.87094688e+01-2.52246150e+01j
 -2.69399238e+01+9.16319391e+01j -2.69399238e+01-9.16319391e+01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j]

/home/mcculler/local/projects_sync/wield-project-lab/wield-iirrational/src/wield/iirrational/fitters_ZPK/codings_s/cplx_sos_NL.py:176: RuntimeWarning: divide by zero encountered in scalar divide
  V_D_c1 = -pD_c1 * c1 / disc * V_D
/home/mcculler/local/projects_sync/wield-project-lab/wield-iirrational/src/wield/iirrational/fitters_ZPK/codings_s/cplx_sos_NL.py:176: RuntimeWarning: invalid value encountered in scalar multiply
  V_D_c1 = -pD_c1 * c1 / disc * V_D
/home/mcculler/local/projects_sync/wield-project-lab/wield-iirrational/src/wield/iirrational/fitters_ZPK/codings_s/cplx_sos_NL.py:177: RuntimeWarning: divide by zero encountered in scalar divide
  V_D_c2 = -pD_c2 / (2 * disc) * V_D
/home/mcculler/local/projects_sync/wield-project-lab/wield-iirrational/src/wield/iirrational/fitters_ZPK/codings_s/cplx_sos_NL.py:177: RuntimeWarning: invalid value encountered in scalar multiply
  V_D_c2 = -pD_c2 / (2 * disc) * V_D
3W   0.56  Fitter_checkpoint improvement succeed, None

------------:Q-ranked order reduction:
4P   1.13    order reduced annealing

3W   3.35    Fitter_checkpoint improvement succeed, None
3W   7.24  Fitter_checkpoint improvement succeed, None

5P   7.25  zero flipping, maxzp 44, residuals=-1.04e-01, -1.04e-01, reldeg=-3
5P   7.80  zero flipped, maxzp 44, residuals=-1.02e-01, reldeg=-3
5P   8.60  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   9.30  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   9.77  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  10.69  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  11.89  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  12.52  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  13.09  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  13.91  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  14.55  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  14.94  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  15.45  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  16.18  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  16.60  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  16.79  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  16.99  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  17.22  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  17.35  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  17.54  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  18.09  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  18.35  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  18.68  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  18.98  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  19.32  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  19.65  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  20.06  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  20.33  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  20.76  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  21.65  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P  21.86  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
2A  88.66  Baseline fit residuals: -1.01e-01, at order 44
BASELINE:  44
------------:investigations:
2I  88.68    max(z, p)       ChiSq.
                   order    avg. res.    med. res.    max. res.
             -----------  -----------  -----------  -----------
                      42    -0.132534    0.0284283      1.05963
                      52    -0.136876    0.011756       4.5478
C [-2.00000000e+03+0.00000000e+00j -3.90376294e-03+0.00000000e+00j
 -1.30278418e+03+2.64289979e-05j -1.30278418e+03-2.64289979e-05j
 -1.11590448e-02+1.16415322e-10j -1.11590448e-02-1.16415322e-10j
 -2.00000000e+02+3.81469727e-06j -2.00000000e+02-3.81469727e-06j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j
 -7.65820729e-01+6.11957961e+00j -7.65820729e-01-6.11957961e+00j
 -7.65820729e-01+6.11957961e+00j -7.65820729e-01-6.11957961e+00j
 -4.64082010e-04+6.19691770e-01j -4.64082010e-04-6.19691770e-01j
 -4.33097913e-04+5.78318501e-01j -4.33097913e-04-5.78318501e-01j
 -3.92408355e-02+7.98284658e-02j -3.92408355e-02-7.98284658e-02j
 -3.92408355e-02+7.98284658e-02j -3.92408355e-02-7.98284658e-02j
 -9.04325228e-03+2.90127484e-01j -9.04325228e-03-2.90127484e-01j
 -8.84127110e-02+1.22469490e+00j -8.84127110e-02-1.22469490e+00j
 -2.88099878e-03+2.81092371e+00j -2.88099878e-03-2.81092371e+00j
 -3.07907738e+00+1.04907173e+01j -3.07907738e+00-1.04907173e+01j
 -1.87094688e+01+2.52246150e+01j -1.87094688e+01-2.52246150e+01j
 -1.87094688e+01+2.52246150e+01j -1.87094688e+01-2.52246150e+01j
 -2.69399238e+01+9.16319391e+01j -2.69399238e+01-9.16319391e+01j
 -8.48230016e+00+8.48230016e-01j -8.48230016e+00-8.48230016e-01j]
(np.float64(-200.0), np.float64(-200.0), np.float64(-1999.9999999999998), np.float64(-1.1743521290869434e-10), np.float64(-0.011159044867910854), np.float64(-0.011159044751801689), np.float64(-1302.7841836435434), np.float64(-1302.7841836432744), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-8.482300164692441+0.8482300164692441j), np.complex128(-8.482300164692441-0.8482300164692441j), np.complex128(-1.276683601197908e-11+91.63193908076074j), np.complex128(-1.276683601197908e-11-91.63193908076074j), np.complex128(-18.709468807770158+25.224614971813303j), np.complex128(-18.709468807770158-25.224614971813303j), np.complex128(-18.709468807798835+25.224614971785446j), np.complex128(-18.709468807798835-25.224614971785446j), np.complex128(-3.5224484988923e-11+10.490717255985867j), np.complex128(-3.5224484988923e-11-10.490717255985867j), np.complex128(-7.755934678337948e-12+2.81092371100195j), np.complex128(-7.755934678337948e-12-2.81092371100195j), np.complex128(-3.1553850493935663e-12+1.2246948964676j), np.complex128(-3.1553850493935663e-12-1.2246948964676j), np.complex128(-8.200117750311727e-13+0.29012748368774854j), np.complex128(-8.200117750311727e-13-0.29012748368774854j), np.complex128(-0.03924083547004946+0.07982846579002265j), np.complex128(-0.03924083547004946-0.07982846579002265j), np.complex128(-0.0392408354694464+0.07982846579497013j), np.complex128(-0.0392408354694464-0.07982846579497013j), np.complex128(-4.26914787276891e-11+0.5783185011865986j), np.complex128(-4.26914787276891e-11-0.5783185011865986j), np.complex128(-4.995095343193837e-11+0.6196917699694076j), np.complex128(-4.995095343193837e-11-0.6196917699694076j), np.complex128(-0.7658207288552369+6.119579611317334j), np.complex128(-0.7658207288552369-6.119579611317334j), np.complex128(-0.7658207289076394+6.119579611415145j), np.complex128(-0.7658207289076394-6.119579611415145j))

3W  88.65  Fitter_checkpoint improvement succeed, None

iir_results.fit_aid._fitters

[Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad19543fe0>,
     log_idx = 0,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad19b17530>,
     log_idx = 1,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad1879d070>,
     log_idx = 12,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad1879db20>,
     log_idx = 14,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad187db170>,
     log_idx = 36,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad187da210>,
     log_idx = 38,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18cab770>,
     log_idx = 62,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18700740>,
     log_idx = 73,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18703620>,
     log_idx = 75,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad186412e0>,
     log_idx = 100,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18760cb0>,
     log_idx = 101,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18643e00>,
     log_idx = 126,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad189d49b0>,
     log_idx = 127,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18641700>,
     log_idx = 152,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad187d9970>,
     log_idx = 153,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18643d10>,
     log_idx = 178,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18ae9df0>,
     log_idx = 179,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18640bc0>,
     log_idx = 204,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18641cd0>,
     log_idx = 205,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18641e50>,
     log_idx = 230,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad187da180>,
     log_idx = 231,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad186a1040>,
     log_idx = 256,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad201cba40>,
     log_idx = 257,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad186dc260>,
     log_idx = 282,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad188cdb50>,
     log_idx = 283,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad186a02c0>,
     log_idx = 308,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad186a0680>,
     log_idx = 309,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad186df350>,
     log_idx = 334,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad186de660>,
     log_idx = 335,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18672330>,
     log_idx = 360,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18671250>,
     log_idx = 361,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18641d60>,
     log_idx = 386,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad186717c0>,
     log_idx = 387,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18672030>,
     log_idx = 412,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18642fc0>,
     log_idx = 413,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad186416d0>,
     log_idx = 438,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad186701d0>,
     log_idx = 439,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18641d00>,
     log_idx = 464,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad188cf950>,
     log_idx = 465,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18640ec0>,
     log_idx = 490,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad1855c3b0>,
     log_idx = 491,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad187031d0>,
     log_idx = 516,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18643d70>,
     log_idx = 517,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18640470>,
     log_idx = 542,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad1855c470>,
     log_idx = 543,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad185ca570>,
     log_idx = 568,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad185c92b0>,
     log_idx = 569,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad1b3321e0>,
     log_idx = 594,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad19da6600>,
     log_idx = 595,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad19eba180>,
     log_idx = 620,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad19dce9c0>,
     log_idx = 621,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad188ce0f0>,
     log_idx = 646,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad19ebabd0>,
     log_idx = 647,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad19540fb0>,
     log_idx = 672,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad19c97ad0>,
     log_idx = 673,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18c5e810>,
     log_idx = 698,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18c5da60>,
     log_idx = 699,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad19eb9e80>,
     log_idx = 724,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad188cd790>,
     log_idx = 725,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad19dceff0>,
     log_idx = 750,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad19ebab10>,
     log_idx = 751,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad19c52090>,
     log_idx = 776,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad188cfb90>,
     log_idx = 777,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad19daf680>,
     log_idx = 802,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad196ac2f0>,
     log_idx = 803,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad1964e870>,
     log_idx = 828,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad19dae660>,
     log_idx = 829,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad1964d610>,
     log_idx = 854,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad1972dc70>,
     log_idx = 855,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad18b08cb0>,
     log_idx = 988,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7fad442f8680>,
     log_idx = 990,
     valid = True,
     )]

And plot the result

#_, TF_all = scipy.signal.freqs_zpk(all_z, all_p, k, worN = F_Hz)
#_, TF_minreal = scipy.signal.freqs_zpk(minreal_all_z, minreal_all_p, minreal_all_k, worN = F_Hz)


plt.loglog(F_bug, Div_PSD**0.5, label = "Full FOM ASD")
plt.loglog(F_FOM_dwn, FOM_ASD, label = "Downsampled FOM ASD")
plt.loglog(F_FOM_dwn, hand_xfr.mag, label = 'Handfit ASD')
plt.loglog(F_Hz, recalibrated_zpk_xfr.mag, label = 'Recalibrated AAA')
plt.loglog(F_FOM_dwn, reduced_ASD, label = "reduced ASD")
#plt.loglog(F_Hz, abs(TF_all), label = 'ALL ZPK Fit Bary. (order {})'.format(results.order))
#plt.loglog(F_Hz, abs(TF_minreal), label = 'Minreal ALL ZPK Fit Bary.', linestyle = '--')
plt.loglog(F_Hz, np.abs(iir_results.fitter.xfer_eval(F_Hz)), label = "data2filter Fit FOM", ls='--', linewidth=4)
plt.grid()
plt.legend()

<matplotlib.legend.Legend at 0x7fad18762480>

iir_results.fitter.xfer_eval(F_Hz)

array([-4.98039622e-01-1.74495439e+00j,  1.12787318e+00-3.58889136e+00j,
        3.36453718e+00-4.56701571e+00j, ...,
       -7.24288676e+18-3.58464631e+19j, -7.24286678e+18-3.58463869e+19j,
       -7.24284680e+18-3.58463107e+19j])

# Python code to merge dict using a single 
# expression
def Merge(dict1, dict2):
	res = {**dict1, **dict2}
	return res
	
# Driver code
dict1 = {'a': 10, 'b': 8}
dict2 = {'d': 6, 'c': 4}
dict3 = Merge(dict1, dict2)
print(dict3)
print(dict1)
print(dict2)

{'a': 10, 'b': 8, 'd': 6, 'c': 4}
{'a': 10, 'b': 8}
{'d': 6, 'c': 4}

# Hand Fit the DARM**2 ZPK
shift_sq = 0.9
z_hand_sq = np.array([-0.07]*18) * 2 * np.pi * shift_sq
mp_sq = -0.75 -0.075j
mul_sq = 9
p_hand_sq = np.array([mp_sq]*mul_sq + [np.real(mp_sq)-np.imag(mp_sq)*1j]*mul_sq)*2 * 2 * np.pi * shift_sq
p_hand_sq = np.append(p_hand_sq, [-200]*2 +[-2000]*1 + [-10000]*0)
# k_hand_sq = 9e39
k_hand_sq = 7e32

#plt.loglog(F_bug, Div_PSD**0.5, label = "Full FOM ASD")
plt.loglog(F_FOM_dwn, FOM_ASD, label = "Downsampled FOM ASD")
#plt.loglog(F_Hz, abs(TF_all), label = 'ALL ZPK Fit Bary. (order {})'.format(results.order))
plt.loglog(F_Hz, np.abs(iir_results.fitter.xfer_eval(F_Hz)), label = "data2filter Fit FOM")
fq_plot = np.logspace(-2, 5, 10000)


hand_sys = SISO.zpk(z_hand_sq, p_hand_sq, k_hand_sq, convention='IIRrational')
fr = hand_sys.fresponse(f=fq_plot)

plt.loglog(-np.real(p_hand_sq), np.array([1e27]*len(p_hand_sq)), color='black', linestyle='', marker='o')
plt.loglog(-np.real(z_hand_sq), np.array([1e27]*len(z_hand_sq)), color='r', linestyle='', marker='x')
plt.loglog(*fr.fplot_mag, label = 'Hand fit to DARM**2 ZPK')
plt.xlabel("Frequency [Hz]")
plt.ylabel("ASD")

plt.grid()
plt.legend()

#plt.savefig('FOM_BNS_fit.pdf', bbox_inches='tight')

/home/mcculler/local/projects_sync/wield-project-lab/wield-control/src/wield/control/SISO/zpk.py:164: NumericalWarning: StateSpace is large (>50 states), using reduced response fiducial auditing heuristics. TODO to make this smarter
  warnings.warn(f"StateSpace is large (>{self.N_MAX_FID} states), using reduced response fiducial auditing heuristics. TODO to make this smarter", util.NumericalWarning)

<matplotlib.legend.Legend at 0x7fad183ba480>

Now we will save our result so it can be loaded by the make full system example file

out_p = iir_results.fitter.poles.fullplane
out_z = iir_results.fitter.zeros.fullplane
out_k = iir_results.fitter.gain

# out_z = z_hand_sq
# out_p = p_hand_sq
# out_k = k_hand_sq

pre_save_sys = SISO.zpk(out_z, out_p, out_k, convention='IIRrational')
from wield.utilities.mpl import mplfigB
plotting_omega = np.logspace(-2, 4, 1000)
plt.close()
axB = mplfigB(Nrows = 2)
fr = pre_save_sys.fresponse(f=plotting_omega)
axB.ax0.loglog(*fr.fplot_mag, label='test')
axB.ax1.semilogx(*fr.fplot_deg225, label='test')


sys_save = SISO.zpk(out_z, out_p, out_k, convention='IIRrational')
zpk_dict = {
            'z' : [str(zero) for zero in np.asarray(tuple(sys_save.z)).tolist()],
            'p' : [str(pole) for pole in np.asarray(tuple(sys_save.p)).tolist()],
            'k' : float(sys_save.k)
            }

save('BNS_FOM.yml', zpk_dict)

Now it is saved. To open it use the following code

with open('BNS_FOM.yml', 'r') as file:
    filt_dict = yaml.safe_load(file)
p = np.array(filt_dict['p'], dtype=np.complex128)
z = np.array(filt_dict['z'], dtype=np.complex128)
k = np.float64(filt_dict['k'])
print(type(k))
print(len(z))
print(len(p))

print(k)
print(z)
print(p)

<class 'numpy.float64'>
39
42
-2.6027954830162637e+34
[-2.06499544e+01 +0.j         -2.31460520e+00 +0.j
 -2.36675714e+04 +0.j         -1.29081156e+01 +0.j
 -1.85249029e+04 +0.j         -4.94022086e+00 +0.j
 -1.29183303e+00 +0.j         -4.94022086e+00 +0.j
 -1.29183303e+00 +0.j         -4.94022086e+00 +0.j
 -1.29183303e+00 +0.j         -4.94022086e+00 +0.j
 -1.29183303e+00 +0.j         -4.94022086e+00 +0.j
 -1.29183303e+00 +0.j         -4.94022086e+00 +0.j
 -1.29183303e+00 +0.j         -4.94022086e+00 +0.j
 -1.29183303e+00 +0.j         -4.94022086e+00 +0.j
 -1.29183303e+00 +0.j         -9.89052486e-01 +0.j
 -3.74963778e-02 +0.j         -7.10782084e-01 +1.50105556j
 -7.10782084e-01 -1.50105556j -3.98864393e+00+13.37090632j
 -3.98864393e+00-13.37090632j -3.98864401e+00+13.37090626j
 -3.98864401e+00-13.37090626j -5.87659068e-01 +1.28119215j
 -5.87659068e-01 -1.28119215j -5.87664308e-01 +1.28119005j
 -5.87664308e-01 -1.28119005j -4.05717328e-02 +0.04292108j
 -4.05717328e-02 -0.04292108j -1.19928970e-01+17.71003385j
 -1.19928970e-01-17.71003385j -7.03816332e-02 +2.25773904j
 -7.03816332e-02 -2.25773904j]
[-5.07240137e+03+0.00000000e+00j -8.37254095e-02+0.00000000e+00j
 -1.35956576e-01+0.00000000e+00j -5.63690459e-02+0.00000000e+00j
 -8.76206516e+01+0.00000000e+00j -3.15322191e+01+0.00000000e+00j
 -8.76206516e+01+0.00000000e+00j -3.15322191e+01+0.00000000e+00j
 -8.76206516e+01+0.00000000e+00j -3.15322191e+01+0.00000000e+00j
 -8.76206516e+01+0.00000000e+00j -3.15322191e+01+0.00000000e+00j
 -8.76206516e+01+0.00000000e+00j -3.15322191e+01+0.00000000e+00j
 -8.76206516e+01+0.00000000e+00j -3.15322191e+01+0.00000000e+00j
 -8.76206516e+01+0.00000000e+00j -3.15322191e+01+0.00000000e+00j
 -8.76206516e+01+0.00000000e+00j -3.15322191e+01+0.00000000e+00j
 -8.76206516e+01+0.00000000e+00j -3.15322191e+01+0.00000000e+00j
 -5.73782567e+00+4.07488474e+01j -5.73782567e+00-4.07488474e+01j
 -5.73782567e+00+4.07488474e+01j -5.73782567e+00-4.07488474e+01j
 -9.17494552e-01+4.67729536e+00j -9.17494552e-01-4.67729536e+00j
 -9.79111001e-03+3.76175886e+00j -9.79111001e-03-3.76175886e+00j
 -1.89534220e-01+6.05891549e-01j -1.89534220e-01-6.05891549e-01j
 -1.82100115e+00+1.68934241e+01j -1.82100115e+00-1.68934241e+01j
 -1.26126746e+00+1.65072231e+00j -1.26126746e+00-1.65072231e+00j
 -1.89534220e-01+6.05891549e-01j -1.89534220e-01-6.05891549e-01j
 -6.24505932e+02+8.48909812e+02j -6.24505932e+02-8.48909812e+02j
 -7.85926170e+03+1.19790252e+04j -7.85926170e+03-1.19790252e+04j]

"""
This File has the sole purpose of creating the ASC DHARD Y Model that is then exported to the buzz repository.
"""

import os
import contextlib

import numpy as np
import matplotlib.pyplot as plt
import control

import pytest
import scipy
from scipy.signal import cheby2
from scipy import signal

from wield.bunch import Bunch
#from wield.control import MIMO
from wield.control import SISO
from wield.control.ss_bare.ss import BareStateSpace as RawStateSpace
from wield.control.ss_bare.ssprint import print_dense_nonzero
from wield.utilities.file_io import load, load_ls, save
from wield.utilities.mpl import mplfigB

from wield.pytest.fixtures import ( 
    tpath_join,
    plot,
    dprint,
    tpath,
    tpath_preclear,
    fpath_join,
    test_trigger,
)

from icecream import ic
from buzz import ssutil

import scipy.optimize
from wield.control import SISO


def bode(
    sys,
    axB=None,
    F_Hz=None,
    omega_limits=None,
    include_zp = True,
    label=None,
    **kwargs
):
    """
    sys should be a wield.control.SISO object
    """
    if axB is None:
        axB = mplfigB(Nrows = 2)

    if include_zp:
        z, p = sys._zp
        z = z[z.imag > 0]
        p = p[p.imag > 0]
        F_include = np.sort(np.concatenate(
            [
                z.imag,
                p.imag,
                z.imag - abs(z.real) - 1e-6,
                z.imag + abs(z.real) + 1e-6,
                p.imag - abs(p.real) - 1e-6,
                p.imag + abs(p.real) + 1e-6,
            ]
        )) / (2 * np.pi)
        F_include = F_include[F_include > 0]

        if omega_limits is None:
            omega_limits = (
                (F_include[0] + 1e-6) / 3,
                F_include[-1] * 3
            )

    if F_Hz is None:
        F_Hz = np.geomspace(omega_limits[0], omega_limits[1], 1000)

    if include_zp:
        print("F_include", F_include)
        F_Hz = np.sort(np.concatenate([F_Hz, F_include]))

    fr = sys.fresponse(f=F_Hz)
    axB.ax0.loglog(*fr.fplot_mag, label=label, **kwargs)
    if hasattr(axB, 'ax1'):
        axB.ax1.semilogx(*fr.fplot_deg225, label=label, **kwargs)
    return axB

FBNS_mod = SISO.zpk(z,p,k, convention='IIRrational')#.asSS
# FBNS_iod = {'F1.in':0, 'F1.out':0}
# FBNS = ssutil.makeSys(FBNS_mod, FBNS_iod)
# omega_lim = [min(F_Hz), max(F_Hz)]

#control.bode(FBNS.mod, Hz=True, wrap_phase=True, omega_limits=omega_lim)

from wield.utilities.mpl import mplfigB
plotting_omega = np.logspace(-2, 4, 1000)
plt.close()
axB = mplfigB(Nrows = 2)
fr = FBNS_mod.fresponse(f=plotting_omega)
axB.ax0.loglog(*fr.fplot_mag, label='test')
axB.ax1.semilogx(*fr.fplot_deg225, label='test')

# tol_list = range(10, 16)
# for tol in tol_list:
#     FBNS_tmp = ssutil.balance_sys_gain(FBNS, nr=tol)
#     nstates_tmp = FBNS.A.shape[0]
#     K_zpk_tmp = SISO.SISOStateSpace(BareStateSpace(FBNS_tmp.A, FBNS_tmp.B, FBNS_tmp.C, FBNS_tmp.D, None)).asZPK
#     control.bode(FBNS.mod, Hz=True, wrap_phase=True, omega_limits=omega_lim, label='sts: {}, tol: {}'.format(nstates_tmp, tol))

plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.savefig('postsave.png', dpi=300, bbox_inches='tight')
plt.show()
plt.close()

/home/mcculler/local/projects_sync/wield-project-lab/wield-control/src/wield/control/SISO/zpk.py:164: NumericalWarning: StateSpace is large (>50 states), using reduced response fiducial auditing heuristics. TODO to make this smarter
  warnings.warn(f"StateSpace is large (>{self.N_MAX_FID} states), using reduced response fiducial auditing heuristics. TODO to make this smarter", util.NumericalWarning)