#!/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 . import plotting_paper as 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_paper'
]
)
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
elif folder == 'ExampleModels_paper':
params["F1_gain"] = np.geomspace(1e2, 1e4, 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_paper'
]
)
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=False,
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_paper'
]
)
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_paper'
]
)
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")