#!/usr/bin/env python
# coding: utf-8
import os
import numpy as np
import matplotlib.pyplot as plt
import pytest
import scipy
from wield.bunch import Bunch
import gwinc
# from wield.control import MIMO
from wield.control import SISO
from icecream import ic
from buzz import ssutil
from wield.utilities import file_io
from buzz import buzzutil
from buzz.ssutil import wieldSS
from buzz import compute
from buzz import plotting
from icecream import ic
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
from matplotlib.colors import LinearSegmentedColormap
import matplotlib.cm as cm
from wield.pytest import (
tjoin,
fjoin,
dprint,
)
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 = 1 * 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_working_F3',
'ExampleModels_newFOM',
'ExampleModels_newFOM_F3'
]
)
def T_calc_H2(folder):
sysB = get_systemB(folder = folder)
io_scales = file_io.load(fjoin(folder, 'io_scales.yml')) # Saved by the normalization in the make function
if '_F3' in folder:
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"]
# currently NO-OP and could be removed
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))
plotting_omega = np.logspace(-3, 4, 1000) * 2 * np.pi
results_H2 = compute.calcFullResults(
results_H2_bare,
colormap="jet",
plot_omega=plotting_omega,
io_scales = io_scales,
)
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_working_F3',
'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"
io_scales = file_io.load(fjoin(folder, 'io_scales.yml')) # Saved by the normalization in the make function
sysB = get_systemB(folder=folder)
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"],
io_scales = io_scales,
)
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")
print('Plotting RMS Plot with lost range')
fig_rms_lost_range_H2 = plotting.plotRMS(
H2_results_E,
setname=setname_H2,
outlineparam='pm',
text='pm',
f1line=True,
h2line=True,
plotrange='diff',
test=True,
#xlim=xlim,
#ylim=ylim
)
fig_rms_lost_range_H2[0].savefig(tjoin("rms_H2_lost_range.pdf"), bbox_inches="tight")
fig_rms_lost_range_H2[0].savefig(tjoin("rms_H2_lost_range.png"), bbox_inches="tight")
print('Plotting RMS Plot with lost range (PSD computation)')
fig_rms_lost_range_H2_psd = plotting.plotRMS(
H2_results_E,
setname=setname_H2,
outlineparam='pm',
text='pm',
f1line=True,
h2line=True,
plotrange='diff_psd',
test=True,
#xlim=xlim,
#ylim=ylim
)
fig_rms_lost_range_H2_psd[0].savefig(tjoin("rms_H2_lost_range_PSD.pdf"), bbox_inches="tight")
fig_rms_lost_range_H2_psd[0].savefig(tjoin("rms_H2_lost_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=False, 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, scale=1/io_scales.w2_rescale)
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_working_F3',
'ExampleModels_newFOM',
'ExampleModels_newFOM_F3'
]
)
def T_calc_HiB(folder):
sysB = get_systemB(folder = folder)
io_scales = file_io.load(fjoin(folder, 'io_scales.yml')) # Saved by the normalization in the make function
if '_F3' in folder:
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"] = [
1e-4, # gain-100 limit
1.8e-2, # 1 deg, around gain-10
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,
)
print("Calculating Full Results Now:")
results_Hib = compute.calcFullResults(
results_Hib_bare,
colormap="RdYlGn",
plot_omega=plotting_omega,
color_param="igsq",
label=False,
io_scales = io_scales,
)
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_working_F3',
'ExampleModels_newFOM',
'ExampleModels_newFOM_F3'
]
)
def T_plot_HiB(folder, dfile = None):
if folder is None:
folder = "ExampleModels"
sysB = get_systemB(folder=folder)
io_scales = file_io.load(fjoin(folder, 'io_scales.yml')) # 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"],
io_scales=io_scales,
)
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")
print('Plotting RMS Plot with lost range (PSD computation)')
fig_rms_lost_range_Hib_psd = plotting.plotRMS(
HinfBounded_results_E + H2_results,
setname=setname_Hib,
outlineparam='pm',
text=plot_text,
f1line=True,
h2line=True,
plotrange='diff_psd',
test=True,
#xlim=xlim,
#ylim=ylim
)
fig_rms_lost_range_Hib_psd[0].savefig(tjoin("rms_HiB_lost_range_PSD.pdf"), bbox_inches="tight")
fig_rms_lost_range_Hib_psd[0].savefig(tjoin("rms_HiB_lost_range_PSD.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,
io_scales=io_scales,
)
+ 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")
try:
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")
except Exception as e:
print("Heatmap Failed:, ", e)
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, test=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, test=True,
)
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], io_scales=io_scales, test=True,
)
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, test=True)
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, scale=1/io_scales.w2_rescale)
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")