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_969275/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 0x7f0410dec380>

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 0x7f04112dbfe0>

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)

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.legend()

<matplotlib.legend.Legend at 0x7f0410abcd70>

# 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 0x7f04107e3aa0>

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 0x7f041050c380>

# 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 0x7f041029c080>

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]
------------:Q-ranked order reduction:
4P   0.17    order reduced annealing

/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.09  Fitter_checkpoint improvement succeed, None
3W   0.69    Fitter_checkpoint improvement succeed, None
3W   1.06  Fitter_checkpoint improvement succeed, None

5P   1.06  zero flipping, maxzp 44, residuals=-1.04e-01, -1.04e-01, reldeg=-3
5P   1.12  zero flipped, maxzp 44, residuals=-1.02e-01, reldeg=-3
5P   1.21  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.28  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.36  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.40  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.44  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.48  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.53  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.58  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.63  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.66  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.69  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.73  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.77  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.81  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.85  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.88  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.91  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.94  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   1.98  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   2.02  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   2.06  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   2.09  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   2.13  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   2.17  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   2.21  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   2.24  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   2.28  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   2.32  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
5P   2.36  zero flipped, maxzp 44, residuals=-1.01e-01, reldeg=-3
2A  10.36  Baseline fit residuals: -1.01e-01, at order 44
BASELINE:  44
------------:investigations:
2I  10.39    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  10.36  Fitter_checkpoint improvement succeed, None

iir_results.fit_aid._fitters

[Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f041074d8b0>,
     log_idx = 0,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f0410dc8560>,
     log_idx = 1,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040fff1040>,
     log_idx = 12,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040fff2000>,
     log_idx = 14,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f0410037860>,
     log_idx = 36,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f04100370b0>,
     log_idx = 38,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f041083d9d0>,
     log_idx = 62,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f0410037920>,
     log_idx = 73,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f0410102b10>,
     log_idx = 75,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040fff2810>,
     log_idx = 100,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f04104b4770>,
     log_idx = 101,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f041008e270>,
     log_idx = 126,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f041083f890>,
     log_idx = 127,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040fff3020>,
     log_idx = 152,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f041022a1b0>,
     log_idx = 153,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f0410848740>,
     log_idx = 178,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040fff3b30>,
     log_idx = 179,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f0410177620>,
     log_idx = 204,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f0410272930>,
     log_idx = 205,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f041008ede0>,
     log_idx = 230,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f041008edb0>,
     log_idx = 231,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f041008d1f0>,
     log_idx = 256,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f0410100aa0>,
     log_idx = 257,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f041008d7c0>,
     log_idx = 282,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f041071aff0>,
     log_idx = 283,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ee0bef0>,
     log_idx = 308,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f0410101a30>,
     log_idx = 309,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f04106baed0>,
     log_idx = 334,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ee0b7a0>,
     log_idx = 335,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ee4e810>,
     log_idx = 360,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040fff2c30>,
     log_idx = 361,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040fff3d70>,
     log_idx = 386,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ee4e360>,
     log_idx = 387,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ee4e780>,
     log_idx = 412,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f041029fdd0>,
     log_idx = 413,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f04102c0170>,
     log_idx = 438,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ee81610>,
     log_idx = 439,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ecf6e40>,
     log_idx = 464,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040fff1a90>,
     log_idx = 465,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f041071b020>,
     log_idx = 490,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040fff08f0>,
     log_idx = 491,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f041053faa0>,
     log_idx = 516,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f041053e150>,
     log_idx = 517,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ecf4950>,
     log_idx = 542,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040eebfb60>,
     log_idx = 543,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ee0a600>,
     log_idx = 568,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ee09d60>,
     log_idx = 569,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f04100374a0>,
     log_idx = 594,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040eebd970>,
     log_idx = 595,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ee0b710>,
     log_idx = 620,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f0410583710>,
     log_idx = 621,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ecf5e50>,
     log_idx = 646,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f0410037b60>,
     log_idx = 647,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ee0a6f0>,
     log_idx = 672,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040eebd7c0>,
     log_idx = 673,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ecf6f30>,
     log_idx = 698,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ecf40e0>,
     log_idx = 699,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f04100ecd70>,
     log_idx = 724,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040eebd5e0>,
     log_idx = 725,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f0410422f30>,
     log_idx = 750,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ecf46b0>,
     log_idx = 751,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ede2810>,
     log_idx = 776,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ede2240>,
     log_idx = 777,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ed4e720>,
     log_idx = 802,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f04101af530>,
     log_idx = 803,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f041008ed20>,
     log_idx = 828,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f041008f410>,
     log_idx = 829,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040ed2aae0>,
     log_idx = 854,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f04104b7cb0>,
     log_idx = 855,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f040edb2000>,
     log_idx = 988,
     valid = True,
     ),
 Bunch(
     checkpoint_idx = -1,
     fitter = <wield.iirrational.fitters_ZPK.MRF.MultiReprFilterS at 0x7f04101af6b0>,
     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 0x7f040ee4c950>

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')

<matplotlib.legend.Legend at 0x7f041029f110>

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)