#!/usr/bin/env python
# coding: utf-8
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
import gwinc
from inspiral_range import inspiral_range as ir
# 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 import (
tjoin,
fjoin,
dprint,
)
from icecream import ic
from buzz import ssutil
import scipy.optimize
from wield.control import SISO
from wield.utilities.file_io import matlab_io, load, save
from buzz import ssutil
from buzz import iodutil
from buzz import buzzutil
from buzz import ssmodels
from buzz.ssutil import wieldSS
from buzz import compute
from buzz import plotting
from icecream import ic
import control
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from wield.control import SISO
from matplotlib.colors import LinearSegmentedColormap
import matplotlib.cm as cm
from wield.control.SISO.zpk_d2c_c2d import d2c_zpk
from . import readFilter
from . import FilteringUtils
# First the State-space model is imported from a file.
from .BH1_solvers import BH1solverset
[docs]
def get_systemB(folder = 'ExampleModels'):
fname = fjoin(folder, "ADY_Sanex.mat")
sys = ssutil.loadSys(fname)
orig_plant_fname = fjoin(folder, "ADY_plant.mat")
orig_plant = ssutil.loadSys(orig_plant_fname)
return Bunch(locals())
#Add scaling for the individual FOMS
total_FOM_scale = 1
FBNS_scale = 1e4 * total_FOM_scale
FFlat_scale = 1 * total_FOM_scale
ct2rad = 5.2e-11 # Calibration from https://alog.ligo-wa.caltech.edu/aLOG/index.php?callRep=69551
strain2m = 1/3995 # Representing the 4km arm length. Strain is h = dL/L where L = 3995 m. The coupling is "I" is in rad to displacement so we need to convert it to strain for the BNS FOM which is in strain
SNR_intg_factor = 1 #2 #from factor of 4 in SNR squared
#divide by sqrt 8 and 2.264
#add avg antenna factor
[docs]
@pytest.mark.parametrize('folder', ['ExampleModels_rescale', 'ExampleModels', 'ExampleModels_working', 'ExampleModels_newFOM', 'ExampleModels_newFOM_F3'])
def T_calc_H2(folder):
sysB = get_systemB(folder = folder)
FBNS_renorm = np.loadtxt(fjoin(folder, 'FBNS_scale.txt')) # Saved by the normalization in the make function
if folder == 'ExampleModels_newFOM_F3':
FOM_out = ["F3.out", "F2.out"]
else:
FOM_out = ["FBNS.out", "F2.out"]
wn_in = ["S.in", "O.in"]
control_in = ["U.in"]
meas_out = ["T.out"]
sysB.sys = ssutil.scale_io(sysB.sys, FOM_out[0], FBNS_scale)
sysB.sys = ssutil.scale_io(sysB.sys, FOM_out[1], FFlat_scale)
# Scale the DARM output
if 'DARM.out' in sysB.sys.outputs.keys():
print('scaling DARM')
DARM_scale_factor = strain2m * SNR_intg_factor
sysB.sys = ssutil.scale_io(sysB.sys, 'DARM.out', DARM_scale_factor)
params = {}
if folder == 'ExampleModels':
params["F1_gain"] = np.geomspace(1e-0, 1e6, 30) * 1/FBNS_scale
elif folder == 'ExampleModels_newFOM_F3':
params["F1_gain"] = np.geomspace(1e-0, 1e6, 30) * 1/FBNS_scale
else:
params["F1_gain"] = np.geomspace(1e2, 1e6, 15) * 1/FBNS_scale
#params["F1_gain"] = np.geomspace(1e2, 5e7, 30) * 1/FBNS_scale
results_H2_bare = compute.calcOpt(
sysB.sys,
control_in,
wn_in,
meas_out,
FOM_out,
params,
solver="LQG",
plant_orig=sysB.orig_plant,
)
print("Length of H2 results: ", len(results_H2_bare))
current_range = np.loadtxt(fjoin(folder, 'FBNS_Current_range.txt')) # Saved by the normalization in the make function
plotting_omega = np.logspace(-3, 4, 1000) * 2 * np.pi
results_H2 = compute.calcFullResults(
results_H2_bare,
colormap="jet",
plot_omega=plotting_omega,
RMS_cal=1,
RMS_cal_F1 = 1/FBNS_scale * 1/FBNS_renorm * strain2m * SNR_intg_factor,
RMS_cal_F2 = 1/FFlat_scale * ct2rad,
current_range = current_range,
)
ofile = "results_H2_ADY.pkl"
tfile = tjoin(ofile)
buzzutil.save_results(results_H2, tfile)
ffile = fjoin(folder, ofile)
buzzutil.save_results(results_H2, ffile)
# copy to the data directory if it is missing
dfile = fjoin(folder, ofile)
# should not copy into rootdir, especially for parametrized test
# should this be for ffile?
# if not os.path.exists(ofile):
# import shutil
# shutil.copy(tfile, ofile)
print("Now running plot test")
T_plot_H2(folder=folder, dfile=tfile)
[docs]
@pytest.mark.parametrize('folder', ['ExampleModels_rescale', 'ExampleModels', 'ExampleModels_working', 'ExampleModels_newFOM', 'ExampleModels_newFOM_F3'])
def T_plot_H2(folder, dfile = None):
"""
Test that runs the plotting code for H2. This test can be run from another test
"""
if folder is None:
folder = "ExampleModels"
sysB = get_systemB(folder=folder)
FBNS_renorm = np.loadtxt(fjoin(folder, 'FBNS_scale.txt')) # Saved by the normalization in the make function
if dfile is None:
ofile = "results_H2_ADY.pkl"
dfile = fjoin(folder, ofile)
print('dfile: ', dfile)
H2_results = buzzutil.load_results(dfile)
print('len H2_results: ', len(H2_results))
# add a gamma value to the results for plotting
for point in H2_results:
if point["gamma"] == None:
point["gamma"] = np.inf
# Remove the labels from the results
for datadict in H2_results:
datadict["label1"] = None
datadict["label2"] = None
# Adding the current controller to the plot
module_gain = -40 # this is the gain of the filter module
filter_list = np.array([1, 3, 4, 5, 6]) # this is a list of the filters in the module that are on
zpksys, filtinfo = readFilter.readFilterSys_Scipy(
fjoin("H1ASC.txt"), "ASC_DHARD_Y", np.array(filter_list)
)
z, p, k = FilteringUtils.d2c(
(zpksys.zeros, zpksys.poles, zpksys.gain), fs=1 / zpksys.dt
)
k = k.real * module_gain
K_mod = SISO.zpk(z, p, k, angular=True, fiducial_rtol=1e-5, fiducial_atol=1e-10).asSS
K_iod = {"C.in": 0, "C.out": 0}
K = wieldSS(K_mod, K_iod)
axB = ssutil.bode(ssutil.asSISO(K)*sysB.orig_plant.siso("P.out", "P.in"), label="K_hand * orig_plant Controller")
axC = ssutil.bode(K, label="K_hand")
K = (ssutil.asSISO(K) * sysB.sys["S.out", "SA.in"]).mimo(
"C.out", "C.in"
) # Add in the part of P that was moved to env. noise block
K = ssutil.balance_sys_gain(K)
axB = ssutil.bode(-1*ssutil.asSISO(K)*H2_results[0]['plant'].siso("P.out", "P.in"), axB =axB, label="K_hand * S_anex * P")
axB = ssutil.bode(H2_results[1]['K'].siso("C.out", "C.in")*H2_results[1]['plant'].siso("P.out", "P.in"), axB =axB, label="K_H2[-15] * P")
axB.save(tjoin("bode_KP.pdf"))
axB.save(tjoin("bode_KP.png"))
axC = ssutil.bode(H2_results[1]['K_usable'], axB =axC, label="K_useable")
axB = ssutil.bode(ssutil.asSISO(H2_results[1]['K'])*sysB.sys.siso("S.out", "SA.in")**(-1), axB =axC, label="K / S_anex")
axC.save(tjoin("bode_K.pdf"))
axC.save(tjoin("bode_K.png"))
curdict = H2_results[0].copy()
curdict["Ac"] = K.A
curdict["Bc"] = K.B
curdict["Cc"] = K.C
curdict["Dc"] = K.D
curdict["K"] = K
curdict["label"] = "Hand Tuned Controller"
extraplots = compute.calcFullResults(
[curdict],
color="dimgrey",
label="Hand Tuned Controller",
plot_omega=curdict["plot_omega"],
)
extraplots[0]["F1_gain"] = sum(buzzutil.listparams(H2_results, "F1_gain"))
extraplots[0]["linecolor"] = "Black" #'DodgerBlue'
extraplots[0]["linestyle"] = "--"
extraplots[0]["solver"] = "Hand Tuned"
extraplots[0]["rms_marker"] = "^"
f_Hz_range = extraplots[0]['FOM1_wn_ASD_Hz']
budget_range = gwinc.load_budget('aLIGO', freq=f_Hz_range)
trace_range = budget_range.run(freq=f_Hz_range)
plt.loglog(f_Hz_range, extraplots[0]['FOM1_wn_ASD']*4000, label='white noise to darm out with hand tuned controller')
plt.loglog(f_Hz_range, trace_range.psd**0.5*4000, label='aLIGO')
plt.legend()
plt.xlabel('Frequency [Hz]')
plt.ylabel('ASD [m/rtHz]')
plt.xlim([6, 6e2])
plt.ylim([4e-21, 3e-17])
plt.savefig(tjoin('DARM_spec.png'), bbox_inches='tight')
plt.savefig(tjoin('DARM_spec.pdf'), bbox_inches='tight')
plt.close()
H2_results_E = H2_results + extraplots
every_nth = 1 # removes every nth element
H2_results_E_short = H2_results[every_nth - 1 :: every_nth] + extraplots
setname_H2 = "ADY_H2"
fig_rms_H2 = plotting.plotRMS(
H2_results_E_short,
setname=setname_H2,
outlineparam="pm",
text="pm",
f1line=False,
h2line=False,
plotrange=True,
test=True,
)
fig_rms_H2[0].savefig(tjoin("rms_H2.pdf"), bbox_inches="tight")
fig_rms_H2[0].savefig(tjoin("rms_H2.png"), bbox_inches="tight")
fig_rms_H2_range_PSD = plotting.plotRMS(
H2_results_E_short,
setname=setname_H2,
outlineparam="pm",
text="pm",
f1line=False,
h2line=False,
plotrange='PSD',
test=True,
)
fig_rms_H2_range_PSD[0].savefig(tjoin("rms_H2_range_PSD.pdf"), bbox_inches="tight")
fig_rms_H2_range_PSD[0].savefig(tjoin("rms_H2_range_PSD.png"), bbox_inches="tight")
ylim = [1e-4, 1e8]
fig_ol_H2, ax = plotting.plotLoop(H2_results_E_short, "OL", setname=setname_H2, ylim=ylim, UG_line=True, test=True)
fig_ol_H2.savefig(tjoin("ol_H2.pdf"), bbox_inches="tight")
fig_ol_H2.savefig(tjoin("ol_H2.png"), bbox_inches="tight")
fig_cl_H2, ax = plotting.plotLoop(H2_results_E_short, "CL", setname=setname_H2, UG_line=True, test=True)
fig_cl_H2.savefig(tjoin("cl_H2.pdf"), bbox_inches="tight")
fig_cl_H2.savefig(tjoin("cl_H2.png"), bbox_inches="tight")
noise_plot_fig, noise_plot_axs = plotting.plot_CLnoises(H2_results_E_short[0], test=True)
noise_plot_fig.savefig(tjoin("noise_H2.pdf"), bbox_inches="tight")
noise_plot_fig.savefig(tjoin("noise_H2.png"), bbox_inches="tight")
fig_darm_H2, ax = plotting.plotDARM(H2_results_E_short, setname=setname_H2, test=True)
fig_darm_H2.savefig(tjoin("All_DARM_Spec.pdf"), bbox_inches="tight")
fig_darm_H2.savefig(tjoin("All_DARM_Spec.png"), bbox_inches="tight")
# setname_vega = 'ASC_vega'
# vega_chart = plotting.vega_plot(H2_results_E, setname=setname_vega, h2line=False, size=500, plotrange=True, test=True)
# vega_chart.save(tjoin("vega_results_plot.html"))
if 'FBNS.DARM.out' in extraplots[0]['CL_all_io'].iod:
DARM_out_name = 'FBNS.DARM.out'
darm_in_mod = True
elif 'F3.DARM.out' in extraplots[0]['CL_all_io'].iod:
DARM_out_name = 'F3.DARM.out'
darm_in_mod = True
else:
DARM_out_name = None
if DARM_out_name:
DARM_plot_omega = np.logspace(-2, 3, 1000) * 2 * np.pi
DARM_noise_plot_fig, DARM_noise_plot_axs = plotting.plot_CLnoises(extraplots[0], extra_fom=DARM_out_name, extra_only=True, plot_omega=DARM_plot_omega, test=True)
DARM_noise_plot_fig.savefig(tjoin("DARM_noise_H2.pdf"), bbox_inches="tight")
DARM_noise_plot_fig.savefig(tjoin("DARM_noise_H2.png"), bbox_inches="tight")
else:
print('DARM not in outputs. Not plotting DARM noise')
print(extraplots[0]['CL_all_io'].outputs)
[docs]
@pytest.mark.parametrize('folder', ['ExampleModels_rescale', 'ExampleModels', 'ExampleModels_working', 'ExampleModels_newFOM', 'ExampleModels_newFOM_F3'])
def T_calc_HiB(folder):
sysB = get_systemB(folder = folder)
FBNS_renorm = np.loadtxt(fjoin(folder, 'FBNS_scale.txt')) # Saved by the normalization in the make function
if folder == 'ExampleModels_newFOM_F3':
FOM_out = ["F3.out", "F2.out"]
else:
FOM_out = ["FBNS.out", "F2.out"]
wn_in = ["S.in", "O.in"]
Zinf = ["Zinf.in"]
control_in = ["U.in"]
meas_out = ["T.out"]
sysB.sys = ssutil.scale_io(sysB.sys, FOM_out[0], FBNS_scale)
sysB.sys = ssutil.scale_io(sysB.sys, FOM_out[1], FFlat_scale)
# Scale the DARM output
if 'DARM.out' in sysB.sys.outputs.keys():
print('scaling DARM')
DARM_scale_factor = strain2m * SNR_intg_factor
sysB.sys = ssutil.scale_io(sysB.sys, 'DARM.out', DARM_scale_factor)
params = {}
if folder == 'ExampleModels':
params['F1_gain'] = np.geomspace(1e1, 1e5, 5) * 1/FBNS_scale # was 1e1 3e3
elif folder == 'ExampleModels_newFOM_F3':
params['F1_gain'] = np.geomspace(1e-6, 1e-3, 5) * 1/FBNS_scale
else:
params["F1_gain"] = np.geomspace(4e1, 1e5, 6) * 1/FBNS_scale
#params["F1_gain"] = np.array([0.09146101038546522])
# params["igsq_capture"] = [
# 0, # H2 constraint
# 1.8e-2, # 1 deg
# 1e-1, # 6 deg
# 0.20, # 11.5 deg
# 0.3, # 17 deg
# 0.5, # 30 deg
# 0.7653668647301797, # 45 deg
# 0.85, # 50.0 deg
# 0.90, # 53.5 deg
# 0.95, # 56.5 deg
# # above this it gets very hairy
# # these are 10 solutions
# ]
# params["igsq_capture"] = np.geomspace(1e-6, 0.99, 40)
params["igsq_capture"] = np.concatenate([np.geomspace(1e-6, 0.1, 6), np.geomspace(.1, 0.95, 6)])
plotting_omega = np.logspace(-3, 4, 1000) * 2 * np.pi
results_Hib_bare = compute.calcOpt(
sysB.sys,
control_in,
wn_in,
meas_out,
FOM_out,
params,
Zinf=Zinf,
solver="HB",
plant_orig=sysB.orig_plant,
BH_solver = BH1solverset,
debug_mode=True,
)
current_range = np.loadtxt(fjoin(folder, 'FBNS_Current_range.txt')) # Saved by the normalization in the make function
print("Calculating Full Results Now:")
results_Hib = compute.calcFullResults(
results_Hib_bare,
colormap="RdYlGn",
plot_omega=plotting_omega,
color_param="igsq",
label=False,
RMS_cal=1,
RMS_cal_F1 = 1/FBNS_scale * 1/FBNS_renorm * strain2m * SNR_intg_factor,
RMS_cal_F2 = 1/FFlat_scale * ct2rad,
current_range = current_range,
)
ofile = "results_Hib_ADY.pkl"
tfile = tjoin(ofile)
buzzutil.save_results(results_Hib, tfile)
ffile = fjoin(folder, ofile)
buzzutil.save_results(results_Hib, ffile)
# TODO, should also check if the data is identical and warn that it should be copied
# copy to the data directory if it is missing
dfile = fjoin(folder, ofile)
# should not save into the root folder for parametrized test - Lee
# if not os.path.exists(ofile):
# import shutil
# shutil.copy(tfile, ofile)
print("Now running plot test")
T_plot_HiB(folder=folder, dfile=tfile)
[docs]
@pytest.mark.parametrize('folder', ['ExampleModels_rescale', 'ExampleModels', 'ExampleModels_working', 'ExampleModels_newFOM', 'ExampleModels_newFOM_F3'])
def T_plot_HiB(folder, dfile = None):
if folder is None:
folder = "ExampleModels"
sysB = get_systemB(folder=folder)
FBNS_renorm = np.loadtxt(fjoin(folder, 'FBNS_scale.txt')) # Saved by the normalization in the make function
H2_results = buzzutil.load_results(fjoin(folder, "results_H2_ADY.pkl"))
print('save fldr: ', fjoin(folder, "results_H2_ADY.pkl"))
if dfile is None:
ofile = "results_Hib_ADY.pkl"
dfile = fjoin(folder, ofile)
results_Hib = buzzutil.load_results(dfile)
# print('Gamma:', buzzutil.listparams(HinfBounded_results, 'gamma'))
# print('igsq:', np.unique(buzzutil.listparams(HinfBounded_results, 'igsq')))
# print('F1_gain', np.unique(buzzutil.listparams(HinfBounded_results, 'F1_gain')))
# print('F1_gain', buzzutil.listparams(HinfBounded_results, 'F1_gain'))
HinfBounded_results = results_Hib
# H2_results = results_H2
# f1_gains = np.unique(buzzutil.listparams(HinfBounded_results, 'F1_gain'))
# HinfBounded_results = buzzutil.filtparam(HinfBounded_results, "F1_gain", f1_gains[3], compare="eq")
for point in HinfBounded_results:
if point["gamma"] == None:
point["gamma"] = np.inf
# Remove labels
for datadict in HinfBounded_results:
datadict["label1"] = None
datadict["label2"] = None
# Adding the current controller to the plot
module_gain = -40 # this is the gain of the filter module
filter_list = np.array([1, 3, 4, 5, 6]) # this is a list of the filters in the module that are on
zpksys, filtinfo = readFilter.readFilterSys_Scipy(
fjoin("H1ASC.txt"), "ASC_DHARD_Y", np.array(filter_list)
)
z, p, k = FilteringUtils.d2c(
(zpksys.zeros, zpksys.poles, zpksys.gain), fs=1 / zpksys.dt
)
k = k.real * module_gain
K_mod = SISO.zpk(z, p, k, angular=True, fiducial_rtol=1e-5, fiducial_atol=1e-10).asSS
K_iod = {"C.in": 0, "C.out": 0}
K = wieldSS(K_mod, K_iod)
K = (K.siso("C.out", "C.in") * sysB.sys["S.out", "SA.in"]).mimo(
"C.out", "C.in"
) # Add in the part of P that was moved to env. noise block
K = ssutil.balance_sys_gain(K)
curdict = HinfBounded_results[0].copy()
curdict["Ac"] = K.A
curdict["Bc"] = K.B
curdict["Cc"] = K.C
curdict["Dc"] = K.D
curdict["K"] = K
curdict["label"] = "Hand Tuned Controller"
extraplots = compute.calcFullResults(
[curdict], color="dimgrey",
label="Hand Tuned Controller",
plot_omega=curdict["plot_omega"],
)
extraplots[0]["F1_gain"] = sum(buzzutil.listparams(HinfBounded_results, "F1_gain"))
extraplots[0]["linecolor"] = "DodgerBlue"
extraplots[0]["linestyle"] = "--"
extraplots[0]["solver"] = "Hand Tuned"
extraplots[0]["rms_marker"] = "^"
HinfBounded_results_E = HinfBounded_results + extraplots
H2_results_E = H2_results + extraplots
every_nth = 2 # removes every nth element
H2_results_E_short = H2_results[every_nth - 1 :: every_nth] + extraplots
xlim= [0.2, 0.5]
ylim = [3e-14, 9e-14]
# Adding the current controller to the plot
setname_Hib = "ADY_Hib"
if len(HinfBounded_results_E + H2_results)>100:
plot_text=False
else:
plot_text='pm'
print('Plotting RMS Plot')
fig_rms_Hib = plotting.plotRMS(
HinfBounded_results_E + H2_results,
setname=setname_Hib,
outlineparam='pm',
text=plot_text,
f1line=True,
h2line=True,
plotrange=False,
#xlim=xlim,
#ylim=ylim
test=True,
)
fig_rms_Hib[0].savefig(tjoin("rms_HiB.pdf"), bbox_inches="tight")
fig_rms_Hib[0].savefig(tjoin("rms_HiB.png"), bbox_inches="tight")
print('Plotting RMS Plot with lost range')
fig_rms_lost_range_Hib = plotting.plotRMS(
HinfBounded_results_E + H2_results,
setname=setname_Hib,
outlineparam='pm',
text=plot_text,
f1line=True,
h2line=True,
plotrange='diff',
test=True,
#xlim=xlim,
#ylim=ylim
)
fig_rms_lost_range_Hib[0].savefig(tjoin("rms_HiB_lost_range.pdf"), bbox_inches="tight")
fig_rms_lost_range_Hib[0].savefig(tjoin("rms_HiB_lost_range.png"), bbox_inches="tight")
f1gains_Hib = np.unique(buzzutil.listparams(HinfBounded_results, "F1_gain"))
try:
HinfBounded_results_single_f1gain = buzzutil.filtparam(HinfBounded_results, "F1_gain", f1gains_Hib[-4], compare="eq")
except:
HinfBounded_results_single_f1gain = buzzutil.filtparam(HinfBounded_results, "F1_gain", f1gains_Hib[0], compare="eq")
HinfBounded_results_single_f1gain = (
compute.calcFullResults(
HinfBounded_results_single_f1gain,
colormap="RdYlGn",
plot_omega=HinfBounded_results[0]["plot_omega"],
color_param="igsq",
label=False,
)
+ extraplots
)
print('Plotting RMS Plot With Range')
fig_rms_range_Hib = plotting.plotRMS(
HinfBounded_results + H2_results,
setname=setname_Hib,
outlineparam='pm',
text=plot_text,
f1line=True,
h2line=True,
plotrange=True,
test=True,
#xlim=xlim,
#ylim=ylim
)
fig_rms_range_Hib[0].savefig(tjoin("rms_HiB_range.pdf"), bbox_inches="tight")
fig_rms_range_Hib[0].savefig(tjoin("rms_HiB_range.png"), bbox_inches="tight")
xlim_PSD=[4.725, 4.805]
print('Plotting RMS Plot With Range from PSD')
fig_rms_range_PSD_Hib = plotting.plotRMS(
HinfBounded_results + H2_results,
setname=setname_Hib,
outlineparam='pm',
text=plot_text,
f1line=True,
h2line=True,
plotrange='PSD',
xlim=xlim_PSD,
test=True,
#ylim=ylim
)
fig_rms_range_PSD_Hib[0].savefig(tjoin("rms_HiB_range_from_PSD.pdf"), bbox_inches="tight")
fig_rms_range_PSD_Hib[0].savefig(tjoin("rms_HiB_range_from_PSD.png"), bbox_inches="tight")
print('Plotting Heatmap')
fig_rmsheatmap = plotting.plotHeatmap(HinfBounded_results + H2_results, test=True)
fig_rmsheatmap[0].savefig(tjoin("heatmap.pdf"), bbox_inches="tight")
fig_rmsheatmap[0].savefig(tjoin("heatmap.png"), bbox_inches="tight")
print('Plotting RMS Plot with Heatmap')
fig_rmsheatmap_Hib = plotting.plotRMS(
HinfBounded_results_E + H2_results,
setname=setname_Hib,
outlineparam='pm',
text=False,
f1line=True,
h2line=True,
plotrange=True,
heatmap=True,
test=True,
)
fig_rmsheatmap_Hib[0].savefig(tjoin("rms_heatmap_HiB.pdf"), bbox_inches="tight")
fig_rmsheatmap_Hib[0].savefig(tjoin("rms_heatmap_HiB.png"), bbox_inches="tight")
print('Plotting RMS Range Plot single result')
fig_single_Hib = plotting.plotRMS(
HinfBounded_results_single_f1gain + H2_results,
setname=setname_Hib,
outlineparam='pm',
text=False,
f1line=True,
h2line=True,
plotrange=True,
heatmap=False,
test=True,
)
fig_single_Hib[0].savefig(tjoin("rms_HiB_singlef1gain_range.pdf"), bbox_inches="tight")
fig_single_Hib[0].savefig(tjoin("rms_HiB_singlef1gain_range.png"), bbox_inches="tight")
"""HinfBounded_results_E = HinfBounded_results + extraplots
H2_results_E = H2_results + extraplots
every_nth = 2 # removes every nth element
H2_results_E_short = H2_results[every_nth-1::every_nth]+ extraplots"""
print('Plotting Open Loop')
ylim = [3e-4, 1e5]
fig_ol_Hib, ax_ol_Hib = plotting.plotLoop(
HinfBounded_results_single_f1gain, "OL", setname=setname_Hib, ylim=ylim, UG_line=True
)
fig_ol_Hib.savefig(tjoin("ol_HiB.pdf"), bbox_inches="tight")
fig_ol_Hib.savefig(tjoin("ol_HiB.png"), bbox_inches="tight")
print('Plotting Closed Loop')
xlim = [1e-1, 1e4]
fig_cl_Hib, ax_cl_Hib = plotting.plotLoop(
HinfBounded_results_single_f1gain, "CL", setname=setname_Hib, UG_line=True, xlim=xlim
)
fig_cl_Hib.savefig(tjoin("cl_HiB.pdf"), bbox_inches="tight")
fig_cl_Hib.savefig(tjoin("cl_HiB.png"), bbox_inches="tight")
print('Plotting gm pm plot')
ylim = None # [0, 12]
plt.clf()
fig_gm_pm, ax_gm_pm = plotting.plot_gm_pm(
HinfBounded_results_E, setname="", f1line=True, text=None, ylim=[0.5,2.5], xlim=[2.5,20]
)
fig_gm_pm.savefig(tjoin("gm_HiB.pdf"), bbox_inches="tight")
fig_gm_pm.savefig(tjoin("gm_HiB.png"), bbox_inches="tight")
# print('Plotting vegaplot')
# setname_vega = 'ASC_vega'
# vega_chart = plotting.vega_plot(HinfBounded_results_E + H2_results, setname=setname_vega, h2line=False, size=500, plotrange=False)
# vega_chart.save(tjoin("vega_results_plot.html"))
DARM_out_name = 'DARM.out'
if DARM_out_name in extraplots[0]['CL_all_io'].outputs.keys():
print('Plotting DARM noise')
DARM_plot_omega = np.logspace(-2, 3, 1000) * 2 * np.pi
DARM_noise_plot_fig, DARM_noise_plot_axs = plotting.plot_CLnoises(extraplots[0], extra_fom=DARM_out_name, extra_only=True, plot_omega=DARM_plot_omega, test=True)
DARM_noise_plot_fig.savefig(tjoin("DARM_noise_H2.pdf"), bbox_inches="tight")
DARM_noise_plot_fig.savefig(tjoin("DARM_noise_H2.png"), bbox_inches="tight")