import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib
import matplotlib.patheffects as pe
import matplotlib.ticker as ticker
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.lines import Line2D
from scipy.interpolate import griddata
from matplotlib.ticker import MaxNLocator
import numpy as np
import os
from inspiral_range import inspiral_range as ir
from buzz import compute
from buzz.compute import color_map_color
from buzz import buzzutil
from buzz import ssutil
import gwinc
from wield.bunch import Bunch
# For Dark Plots
#plt.style.use('dark_background')
#import matplotlib as mpl
#mpl.rcParams.update(mpl.rcParamsDefault)
# for full list https://matplotlib.org/stable/tutorials/introductory/customizing.html
style_6p5in = {
"savefig.dpi": 400,
"figure.dpi": 400,
"lines.linewidth": 1.25,
"axes.facecolor": "FFFFFF",
"figure.facecolor": "FFFFFF",
"figure.figsize": [6.5**2/5.72, 3.0**2/2.76],
"mathtext.fontset": "cm",
"font.family": "DejaVu Serif",
"figure.titlesize": "large", # size of the figure title (Figure.suptitle())
"figure.titleweight": "normal", # weight of the figure title
"pdf.fonttype": 3, # Output Type 3 (Type3) or Type 42 (TrueType)
"font.size": 10.0,
"axes.grid": True, # display grid or not
"axes.grid.axis": "both", # which axis the grid should apply to
"axes.grid.which": "major", # gridlines at {major, minor, both} ticks
"axes.titlesize": "large", # fontsize of the axes title
"axes.titlepad": 6.0, # pad between axes and title in points
"axes.labelsize": 10, # fontsize of the x any y labels
"axes.labelpad": 3.0, # space between label and axis
"axes.formatter.useoffset": False, # If True, the tick label formatter
"axes.formatter.offset_threshold": 4, # When useoffset is True, the offset
"legend.fontsize": 8.0, #"medium",
"legend.borderpad": 0.4, # border whitespace
"legend.labelspacing": 0.5, # the vertical space between the legend entries
"legend.handlelength": 2.0, # the length of the legend lines
"legend.handleheight": 0.7, # the height of the legend handle
"legend.handletextpad": 0.8, # the space between the legend line and legend text
"legend.borderaxespad": 0.5, # the border between the axes and legend edge
"legend.columnspacing": 2.0, # column separation
'xtick.major.size' : 6,
'xtick.minor.size' : 3,
'ytick.major.size' : 6,
'ytick.minor.size' : 3,
'xtick.direction': 'in', # direction: {in, out, inout}
'ytick.direction': 'in', # direction: {in, out, inout}
'xtick.minor.bottom': True, # draw x axis bottom minor ticks
'xtick.minor.visible': True, # visibility of minor ticks on x-axis
'grid.alpha': 0.5,
}
style_3p5in = {
"savefig.dpi": 400,
"figure.dpi": 400,
"lines.linewidth": 1.0,
"axes.facecolor": "FFFFFF",
"figure.facecolor": "FFFFFF",
"figure.figsize": [3.5+0.08, 2.5+0.13],
"mathtext.fontset": "cm",
"font.family": "DejaVu Serif",
"figure.titlesize": "medium", # size of the figure title (Figure.suptitle())
"figure.titleweight": "normal", # weight of the figure title
"pdf.fonttype": 3, # Output Type 3 (Type3) or Type 42 (TrueType)
"font.size": 10.0,
"axes.grid": False, # display grid or not
"axes.titlesize": 10, # fontsize of the axes title
"axes.titlepad": 6.0, # pad between axes and title in points
"axes.labelsize": 10, # fontsize of the x any y labels
"axes.labelpad": 0.0, # space between label and axis
"axes.formatter.useoffset": True, # If True, the tick label formatter
"axes.formatter.offset_threshold": 4, # When useoffset is True, the offset
"legend.fontsize": 8,
"legend.borderpad": 0.4, # border whitespace
"legend.labelspacing": 0.5, # the vertical space between the legend entries
"legend.handlelength": 2.0, # the length of the legend lines
"legend.handleheight": 0.7, # the height of the legend handle
"legend.handletextpad": 0.8, # the space between the legend line and legend text
"legend.borderaxespad": 0.5, # the border between the axes and legend edge
"legend.columnspacing": 2.0, # column separation
'xtick.major.size' : 6,
'xtick.minor.size' : 3,
'ytick.major.size' : 6,
'ytick.minor.size' : 3,
'xtick.direction': 'in', # direction: {in, out, inout}
'ytick.direction': 'in', # direction: {in, out, inout}
'xtick.minor.bottom': True, # draw x axis bottom minor ticks
'ytick.minor.left': True, # draw x axis bottom minor ticks
'ytick.major.left': True, # draw x axis bottom minor ticks
'xtick.minor.visible': True, # visibility of minor ticks on x-axis
}
# def color_map_color(value, cmap_name="winter", vmin=0, vmax=1):
# if type(value)==list or type(value)==np.ndarray:
# colors=[]
# for val in value:
# colors.append(color_map_color(val, cmap_name=cmap_name, vmin=min(value), vmax=max(value)))
# return colors
# if value < vmin:
# value = vmin
# if value > vmax:
# value = vmax
# norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax)
# cmap = cm.get_cmap(cmap_name) # PiYG
# cmap.set_bad('lightgrey')
# rgb = cmap(norm(value))[:3] # will return rgba, we take only first 3 so we get rgb
# color = matplotlib.colors.rgb2hex(rgb)
# return color
[docs]
def makePlotFolder(setname=''):
"""Makes a folder to store the plots in if it doesn't already exist
"""
if setname=='':
if not os.path.exists('plots'):
os.makedirs('plots')
print('Created /plots folder')
else:
makePlotFolder()
if not os.path.exists('plots/' + setname):
os.makedirs('plots/' + setname)
print('Created /plots/' + setname + ' folder')
else:
pass
return
[docs]
def addplotcolor(results,color):
for result in results:
result['plotcolor'] = color
return results
[docs]
def plotLoop(results, loop, extraplots=None, setname = '', test=False, xlim=None, ylim=None, UG_line=False, save=True):
if setname != '':
setname_folder = setname + '/'
setname_u = setname + '_'
else:
setname_folder = ''
matplotlib.rcParams.update(style_6p5in) # Add the full width style
matplotlib.rcParams.update({"figure.figsize": [6.5**2/5.09, 3**2/2.74]})
fig_l, axs_l = plt.subplots(2, 1, sharex=True, gridspec_kw={'height_ratios': [2.5,1]})
#axs_l[0].set_title(loop+' Gain')
axs_l[0].set_ylabel('Magnitude (Abs)')
axs_l[0].set_yscale('log')
axs_l[1].set_ylabel('Phase (Deg)')
axs_l[1].set_xlabel('Frequency (Hz)')
axs_l[1].set_xscale('log')
axs_l[0].grid(alpha=0.5)
axs_l[1].grid(alpha=0.5)
x_toplot_name = 'FR_' + loop + '_omega_Hz'
y_mag_toplot_name = 'FR_' + loop + '_mag'
y_ph_toplot_name = 'FR_' + loop + '_ph_deg'
#y_ph_toplot_name = 'FR_' + loop + '_ph_unwrapped_Hz' # For Unwrapped Phase
axs_l[1].set_xlim([min(results[0][x_toplot_name]), max(results[0][x_toplot_name])])
for idx, result in enumerate(results):
if result['label1'] == '':
axs_l[0].plot(result[x_toplot_name], result[y_mag_toplot_name], color=result['linecolor'], linestyle=result['linestyle'])
else:
axs_l[0].plot(result[x_toplot_name], result[y_mag_toplot_name], color=result['linecolor'], linestyle=result['linestyle'], label=result['label1'])
if result['label2'] == '':
axs_l[1].plot(result[x_toplot_name], result[y_ph_toplot_name], color=result['linecolor'], linestyle=result['linestyle'])
else:
axs_l[1].plot(result[x_toplot_name], result[y_ph_toplot_name], color=result['linecolor'], linestyle=result['linestyle'], label=result['label2'])
if UG_line:
for result in results:
axs_l[0].plot(result[x_toplot_name], np.ones(len(result[x_toplot_name])), color='black', linestyle='--', alpha=0.4, label='Unity Gain', linewidth=0.6)
break
plt.subplots_adjust(hspace=0.1)
#axs_l[0].legend(loc='upper right')
leg = axs_l[0].legend()
#axs_l[0].legend(loc='center left', bbox_to_anchor=(1, 0.5))
#axs_l[1].legend(loc='center left', bbox_to_anchor=(1, 0.5))
for line in leg.get_lines():
line.set_linewidth(1.5)
if xlim:
axs_l[0].set_xlim(xlim)
if ylim:
axs_l[0].set_ylim(ylim)
norm = results[0]['colormap_norm']
if results[0]['color_param'] == 'F1_gain':
cbar_ticks = None
loop_cbar_label = 'BNS FOM Weight, $\zeta$'
elif results[0]['color_param'] == 'gamma':
cbar_ticks = np.rint(np.geomspace(norm.vmin, norm.vmax, 7)).astype(int)
loop_cbar_label = '$\mathcal{H}_\infty$ bound, $\gamma$'
elif results[0]['color_param'] == 'igsq':
cbar_ticks = None
loop_cbar_label ='Inverse $\mathcal{H}_\infty$ bound squared $(1/\gamma^2)$'
else:
cbar_ticks = None
loop_cbar_label = results[0]['color_param']
if cbar_ticks is None:
cbar = fig_l.colorbar(cm.ScalarMappable(norm=norm, cmap=results[0]['colormap']), label=loop_cbar_label, ax=axs_l)
else:
cbar = fig_l.colorbar(cm.ScalarMappable(norm=norm, cmap=results[0]['colormap']), ax=axs_l, label=loop_cbar_label, ticks=cbar_ticks)
cbar.ax.set_yticklabels(cbar_ticks)
# Adjust the distance between the colorbar and the plot
cbar.ax.yaxis.set_offset_position('right') # To keep the offset on the right side
cbar.ax.yaxis.set_label_position('right')
cbar.ax.set_position([0.76, 0.11, 0.03, 0.77]) # [left, bottom, width, height]
cbar.ax.yaxis.set_label_coords(4, 0.5) # Adjust the position of the label (optional)
if save:
if not test:
makePlotFolder(setname=setname)
save_fname ='plots/'+ setname_folder + setname_u + loop + '_BODE'
fig_l.savefig(save_fname + '.pdf', bbox_inches='tight', pad_inches=0.05)
fig_l.savefig(save_fname+ '.png', bbox_inches='tight', pad_inches=0.05)
plt.show()
return fig_l, axs_l
[docs]
def plotRMS(
results,
extras=None,
setname='',
outlineparam='gm',
text=False,
test=False,
plotrange=False,
zero_current_range=True,
f1line=False,
h2line=False,
xlim=None,
ylim=None,
ax_labels=None,
title=False,
cb2_label=None,
save=True,
heatmap=False,
onlyCLstable=False,
):
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib
import matplotlib.patheffects as pe
import matplotlib.ticker as ticker
from matplotlib.colors import LinearSegmentedColormap
from matplotlib.lines import Line2D
from scipy.interpolate import griddata
import numpy as np
import os
from buzz import compute
from buzz.compute import color_map_color
from buzz import buzzutil
from buzz import ssutil
from wield.bunch import Bunch
if setname != '':
setname_folder = setname + '/'
setname_u = setname + '_'
else:
setname_folder = ''
results = results.copy()
matplotlib.rcParams.update(style_6p5in) # Add the full width style
if onlyCLstable:
all_ctrl_len = len(results)
results = buzzutil.filtparam(results, 'isCLstable', True, compare='equal')
print('removed {} CL unstable controllers'.format(all_ctrl_len - len(results)))
if h2line:
h2line_results = buzzutil.filtparam(results, 'solver', 'LQG', compare='equal')
results = buzzutil.filtparam(results, 'solver', 'LQG', compare='not equal')
else:
h2line_results = []
fom1_name = results[0]['FOM_out'][0].split('.')[0]
fom2_name = results[0]['FOM_out'][1].split('.')[0]
fig = plt.figure(figsize=(6.5**2/5.77, 3.0**2/2.86))
fig.subplots_adjust(wspace=0.35)
plt.rc('axes', axisbelow=True)
ax = fig.add_axes([0.07, 0.07, 0.7, 0.8])
axs1 = fig.add_axes([0.8, 0.07, 0.02, 0.8])
axs2 = fig.add_axes([0.92, 0.07, 0.02, 0.8])
if title:
ax.set_title(fom1_name + " RMS vs. " + fom2_name + " RMS")
#ax.set_xlabel(fom1_name + " RMS")
if ax_labels:
ax.set_xlabel(ax_labels[0])
ax.set_ylabel(ax_labels[1])
else:
if plotrange:
if plotrange=='diff' or plotrange=='diff_range':
ax.set_xlabel('Range Lost (Mpc)')
else:
ax.set_xlabel('Range (Mpc)')
else:
ax.set_xlabel('Scaled RMS From BNS FOM')
ax.set_ylabel('Scaled RMS From Flat FOM')
ax.grid(alpha=0.5)
fom1_H2_norm = []
fom2_H2_norm = []
FOM1_gain = []
RMS_markers = []
for result in results:
FOM1_gain.append(result['F1_gain'])
RMS_markers.append(result['rms_marker'])
if plotrange:
if plotrange=='diff' or plotrange=='diff_range':
fom1_H2_norm.append(result['abs_diff_range'])
elif plotrange=='diff_psd':
fom1_H2_norm.append(result['abs_diff_range_PSD'])
elif plotrange=='psd' or plotrange=='PSD':
fom1_H2_norm.append(result['range_from_PSD'])
else:
fom1_H2_norm.append(result['range'])
else:
fom1_H2_norm.append(result['F1_H2_Norm'])
fom2_H2_norm.append(result['F2_H2_Norm'])
#color_list = color_map_color(FOM1_gain)
if results[0]['solver'] == 'LQG':
resultcolorname = 'F1_gain'
colorlog = True
else:
resultcolorname = 'igsq'
colorlog = False
outlinecolor_nonabs = buzzutil.listparams(results + h2line_results, outlineparam)
outlinecolor_values = np.abs(outlinecolor_nonabs)
# outlinecolor_tmp = []
# for v in outlinecolor:
# if not (np.isinf(v) or np.isnan(v)):
# outlinecolor_tmp.append(v)
# else:
# outlinecolor_tmp.append(v)#'white')
# outlinecolor = outlinecolor_tmp
#outline_cmap = 'inferno'
if heatmap:
outline_cmap = 'jet_r'
else:
outline_cmap = 'YlGnBu_r'
outlinecolor_list = color_map_color(outlinecolor_values, cmap_name=outline_cmap)
if heatmap:
heatmap_f1_points = []
heatmap_f2_points = []
heatmap_color = []
color_list = buzzutil.listparams(results + h2line_results, 'linecolor')
for idx, result in enumerate(results + h2line_results):
if plotrange:
if plotrange=='diff' or plotrange=='diff_range':
f1_plotpoint = result['abs_diff_range']
elif plotrange=='diff_psd':
f1_plotpoint = result['abs_diff_range_PSD']
elif plotrange=='psd' or plotrange=='PSD':
f1_plotpoint = result['range_from_PSD']
else:
f1_plotpoint = result['range']
else:
f1_plotpoint = result['F1_H2_Norm']
#ax.scatter(f1_plotpoint, result['F2_H2_Norm'], color=color_list[idx], s=100, edgecolors=outlinecolor_list[idx], marker=result['rms_marker'], linewidth=4, label=result['label1'])
if heatmap and result['solver'] != 'Hand Tuned':
heatmap_f1_points.append(f1_plotpoint)
heatmap_f2_points.append(result['F2_H2_Norm'])
heatmap_color.append(outlinecolor_values[idx])
#ax.scatter(f1_plotpoint, result['F2_H2_Norm'], color='black', s=50, marker='o', linewidth=1)
ax.scatter(f1_plotpoint, result['F2_H2_Norm'], color=outlinecolor_list[idx], s=30, edgecolors='black', linewidth=1)
else:
if not (h2line and result['solver']=='LQG'):
ax.scatter(f1_plotpoint, result['F2_H2_Norm'], color=color_list[idx], s=70, edgecolors=outlinecolor_list[idx], marker=result['rms_marker'], linewidth=2.75, label=result['label1'])
if f1line:
for gain in np.unique(buzzutil.listparams(results, 'F1_gain')):
results_line = buzzutil.filtparam(results, 'F1_gain', [gain])
fom1_H2_norm_line = []
fom2_H2_norm_line = []
for result in results_line:
if plotrange:
if plotrange=='diff' or plotrange=='diff_range':
fom1_H2_norm_line.append(result['abs_diff_range'])
elif plotrange=='diff_psd':
fom1_H2_norm_line.append(result['abs_diff_range_PSD'])
elif plotrange=='psd' or plotrange=='PSD':
fom1_H2_norm_line.append(result['range_from_PSD'])
else:
fom1_H2_norm_line.append(result['range'])
else:
fom1_H2_norm_line.append(result['F1_H2_Norm'])
fom2_H2_norm_line.append(result['F2_H2_Norm'])
if gain == np.unique(buzzutil.listparams(results, 'F1_gain'))[-1]:
f1gain_lab = 'Lines of Constant $\zeta$'
else:
f1gain_lab = None
ax.plot(fom1_H2_norm_line, fom2_H2_norm_line, color='black', linewidth=1.4, linestyle='dotted', alpha=0.7, label=f1gain_lab)
if h2line:
h2_line_f1_rms = []
h2_line_f2_rms = []
for result in h2line_results:
if plotrange:
if plotrange=='diff' or plotrange=='diff_range':
h2_line_f1_rms.append(result['abs_diff_range'])
elif plotrange=='diff_psd':
h2_line_f1_rms.append(result['abs_diff_range_PSD'])
elif plotrange=='psd' or plotrange=='PSD':
h2_line_f1_rms.append(result['range_from_PSD'])
else:
h2_line_f1_rms.append(result['range'])
else:
h2_line_f1_rms.append(result['F1_H2_Norm'])
h2_line_f2_rms.append(result['F2_H2_Norm'])
#plot dotted line
#ax.plot(h2_line_f1_rms, h2_line_f2_rms, color='black', linewidth=1.5, linestyle='dashed', alpha=0.2)
ax.plot(h2_line_f1_rms, h2_line_f2_rms, color='red', linestyle='dashed', alpha=0.5, label='H2 Optimal Limit')
if text:
if isinstance(text, str):
if text == 'idx':
text_vals_non_int =np.linspace(0,len(results)-1,len(results))
text_vals = np.array([int(val) for val in text_vals_non_int])
else:
try:
text_vals = np.abs(buzzutil.listparams(results, text))
except:
KeyError('Text parameter not found in results')
else:
text_vals = outlinecolor_nonabs
for i, txt in enumerate(text_vals):
ax.annotate(txt, (fom1_H2_norm[i], fom2_H2_norm[i]), color='black',xytext=(5, 5), textcoords="offset pixels", path_effects=[pe.withStroke(linewidth=1.5, foreground='w')])
if text!=outlineparam:
if text == 'gm':
text_ann = 'Abs Gain Margin (dB) (Text Near Point)'
elif text == 'pm':
text_ann = 'Phase Margin (deg) (Text Near Point)'
elif text == 'idx':
text_ann = 'index of controller in list (Text Near Point)'
else:
text_ann = outlineparam + ' (Text Near Point)'
ax.text(0.014,0.01,text_ann, horizontalalignment='left',verticalalignment='bottom', transform=ax.transAxes)
ax.set_xscale('log')
ax.set_yscale('log')
# # Force the x-axis to use a logarithmic scale
# ax.xaxis.set_major_locator(ticker.LogLocator(base=10.0, subs='auto', numticks=10))
# ax.xaxis.set_minor_locator(ticker.LogLocator(base=10.0, subs='auto', numticks=10))
# ax.xaxis.set_major_formatter(ticker.LogFormatter())
# ax.xaxis.set_minor_formatter(ticker.LogFormatter(minor_thresholds=(2, 0.4)))
# if plotrange:
# # Prevent scientific notation on x-axis
# ax.xaxis.set_major_locator(ticker.MaxNLocator(integer=True))
# ax.xaxis.set_minor_locator(ticker.MaxNLocator(integer=True))
# ax.xaxis.set_major_formatter(ticker.ScalarFormatter())
# ax.xaxis.set_minor_formatter(ticker.ScalarFormatter())
leg = ax.legend()
# change the line width for the legend
for line in leg.get_lines():
line.set_linewidth(1.0)
if xlim or ylim:
if xlim=='H2' or xlim=='h2' or ylim=='H2' or ylim=='h2':
if xlim=='H2' or xlim=='h2':
xlim_h2 = [min(h2_line_f1_rms), max(h2_line_f1_rms)]
ax.set_xlim(xlim_h2)
if ylim=='H2' or ylim=='h2':
ylim_h2 = [min(h2_line_f2_rms), max(h2_line_f2_rms)]
ax.set_ylim(ylim_h2)
else:
if xlim and xlim!='H2' and xlim!='h2':
ax.set_xlim(xlim)
if ylim and ylim!='H2' and ylim!='h2':
ax.set_ylim(ylim)
else:
pass
outlinecolor_vmin=np.nanmin([v for v in outlinecolor_values if v != np.inf and v != -np.inf and v != None])
outlinecolor_vmax=np.nanmax([v for v in outlinecolor_values if v != np.inf and v != -np.inf and v != None])
norm1 = matplotlib.colors.Normalize(vmin=outlinecolor_vmin, vmax=outlinecolor_vmax)
cb1 = matplotlib.colorbar.ColorbarBase(axs1, cmap=outline_cmap, norm=norm1)
if outlineparam == 'gm':
cb1.set_label('Abs Gain Margin (dB) (outline color)')
elif outlineparam == 'pm':
cb1.set_label('Phase Margin (deg) (outline color)')
else:
cb1.set_label(outlineparam + ' (outline color)')
norm2 = buzzutil.listparams(results, 'colormap_norm')[0]
cmapname = buzzutil.listparams(results, 'colormap')[0]
cb2 = matplotlib.colorbar.ColorbarBase(axs2, cmap=cmapname, norm=norm2)
if cb2_label:
cb2.set_label(cb2_label)
elif 'color_param' in results[0]:
cb2_label = buzzutil.listparams(results, 'color_param')[0]
if cb2_label == 'F1_gain':
cb2.set_label('BNS FOM Weight, $\zeta$')
elif cb2_label == 'igsq':
cb2.set_label('Inverse $\mathcal{H}_\infty$ bound squared $(1/\gamma^2)$')
elif cb2_label == 'gamma':
cb2.set_label('$\mathcal{H}_\infty$ bound, $\gamma$')
else:
cb2.set_label(cb2_label)
else:
cb2.set_label('BNS FOM Weight, $\zeta$')
if heatmap:
# Transform data to log scale
log_x = np.log(heatmap_f1_points)
log_y = np.log(heatmap_f2_points)
grid_log_x = np.linspace(log_x.min(), log_x.max(), 500)
grid_log_y = np.linspace(log_y.min(), log_y.max(), 500)
grid_log_x, grid_log_y = np.meshgrid(grid_log_x, grid_log_y)
# Interpolate the data in log scale
grid_v = griddata((log_x, log_y), heatmap_color, (grid_log_x, grid_log_y), method='cubic')
# Filter out points between the H2 line and the origin in the mesh grid
mesh_pts_x, mesh_pts_y = np.exp(grid_log_x), np.exp(grid_log_y)
h2_pts_1 = np.array([h2_line_f1_rms[:-1], h2_line_f2_rms[:-1]]).T
h2_pts_2 = np.array([h2_line_f1_rms[1:], h2_line_f2_rms[1:]]).T
xmax_axis = np.array([np.maximum(h2_line_f1_rms[:-1], h2_line_f1_rms[1:]), np.zeros(len(h2_line_f1_rms) - 1)]).T
ymax_axis = np.array([np.zeros(len(h2_line_f2_rms) - 1), np.maximum(h2_line_f2_rms[:-1], h2_line_f2_rms[1:])]).T
def is_point_in_triangle_vectorized(pts_x, pts_y, v1, v2, v3):
def sign(p1, p2, p3):
return (p1[0] - p3[0]) * (p2[1] - p3[1]) - (p2[0] - p3[0]) * (p1[1] - p3[1])
d1 = sign((pts_x, pts_y), v1, v2)
d2 = sign((pts_x, pts_y), v2, v3)
d3 = sign((pts_x, pts_y), v3, v1)
has_neg = (d1 < 0) | (d2 < 0) | (d3 < 0)
has_pos = (d1 > 0) | (d2 > 0) | (d3 > 0)
return ~(has_neg & has_pos)
mask = np.ones(grid_v.shape, dtype=bool)
for k in range(len(h2_line_f1_rms) - 1):
mask &= ~is_point_in_triangle_vectorized(mesh_pts_x, mesh_pts_y, h2_pts_1[k], h2_pts_2[k], (0, 0))
mask &= ~is_point_in_triangle_vectorized(mesh_pts_x, mesh_pts_y, h2_pts_1[k], xmax_axis[k], (0, 0))
mask &= ~is_point_in_triangle_vectorized(mesh_pts_x, mesh_pts_y, ymax_axis[k], h2_pts_2[k], (0, 0))
grid_v = np.where(mask, grid_v, np.nan)
ax.pcolormesh(np.exp(grid_log_x), np.exp(grid_log_y), grid_v, alpha=0.5, cmap=outline_cmap, norm=norm1)
#fig.set_size_inches(3*1.10, 2*1.10)
if save:
if not test:
makePlotFolder(setname=setname)
save_fname ='plots/'+ setname_folder + setname_u + 'FOM_RMS_Scatter'
fig.savefig(save_fname + '.pdf', bbox_inches='tight', pad_inches=0.05)
fig.savefig(save_fname + '.png', bbox_inches='tight', pad_inches=0.05)
plt.show()
return fig, ax
[docs]
def plotHeatmap(results, param='pm', abs=True, setname='', test=False, save=True):
results = buzzutil.filtparam(results, 'solver', 'Hand Tuned', compare='not equal')
# Extract data from results
x = np.array(buzzutil.listparams(results, 'F1_H2_Norm'))
y = np.array(buzzutil.listparams(results, 'F2_H2_Norm'))
v = np.array(buzzutil.listparams(results, param))
if abs:
v = np.abs(v)
c_vmin=np.nanmin([h for h in v if h != np.inf and h != -np.inf and h != None])
c_vmax=np.nanmax([h for h in v if h != np.inf and h != -np.inf and h != None])
for idx, val in enumerate(v):
if val == np.inf or val == np.nan or val == None:
v[idx] = c_vmax
if val == -np.inf:
v[idx] = c_vmin
norm1 = matplotlib.colors.Normalize(vmin=c_vmin, vmax=c_vmax)
# Transform data to log scale
log_x = np.log(x)
log_y = np.log(y)
grid_log_x = np.geomspace(log_x.min(), log_x.max(), 500)
grid_log_y = np.geomspace(log_y.min(), log_y.max(), 500)
grid_log_x, grid_log_y = np.meshgrid(grid_log_x, grid_log_y)
# Interpolate the data in log scale
grid_v = griddata((log_x, log_y), v, (grid_log_x, grid_log_y), method='cubic')
# Create the plot with fig and ax
fig, ax = plt.subplots(figsize=(8, 6))
cax = ax.pcolormesh(np.exp(grid_log_x), np.exp(grid_log_y), grid_v, shading='auto', cmap='viridis', alpha=0.8, norm=norm1)
fig.colorbar(cax, ax=ax, label='abs(phase margin (deg))')
ax.scatter(x, y, c=v, cmap='viridis', edgecolor='k', norm=norm1)
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_xlabel('FOM_1_H2_Norm')
ax.set_ylabel('FOM_2_H2_Norm')
# Save the plot if required
if save:
if not test:
makePlotFolder(setname=setname)
save_fname = 'plots/' + setname + '_heatmap'
fig.savefig(save_fname + '.pdf', bbox_inches='tight', pad_inches=0.05)
fig.savefig(save_fname + '.png', bbox_inches='tight', pad_inches=0.05)
plt.show()
return fig, ax
[docs]
def plotNoise(result, setname='', test=False):
if not isinstance(result, (dict, Bunch)):
print('result must be a dictionary using first entry')
result = result[0]
matplotlib.rcParams.update(style_6p5in) # Add the full width style
if setname != '':
setname_folder = setname + '/'
setname_u = setname + '_'
else:
setname_folder = ''
#matplotlib.rcParams.update(style_3p5in) # Add the half width style
fig, axs = plt.subplots(2, 1, figsize=(10, 10), sharex=True)
axs[0].set_title('Noise to Measured Output')
axs[0].set_ylabel('Magnitude abs')
axs[0].set_yscale('log')
axs[1].set_ylabel('Phase')
axs[1].set_xlabel('Frequency (Hz)')
#axs[1].set_xlim([min(results[0]['FR_OL_omega']), max(results[0]['FR_OL_omega'])])
axs[1].set_xscale('log')
axs[0].grid(alpha=0.9)
axs[1].grid(alpha=0.9)
for key in result['FR_WN'].keys():
axs[0].plot(result['FR_WN'][key]['omega'], result['FR_WN'][key]['mag'], label=key)
axs[1].plot(result['FR_WN'][key]['omega'], result['FR_WN'][key]['phase'])
axs[0].set_xlim([min(result['FR_WN'][key]['omega']), max(result['FR_WN'][key]['omega'])])
#axs[0].legend(loc='center left', bbox_to_anchor=(1, 0.5))
axs[0].legend()
if not test:
makePlotFolder(setname=setname)
save_fname ='plots/'+ setname_folder + setname_u + 'Noise_BODE'
fig.savefig(save_fname + '.pdf', bbox_inches='tight', pad_inches=0.05)
fig.savefig(save_fname+ '.png', bbox_inches='tight', pad_inches=0.05)
plt.show()
return fig, axs
[docs]
def plotDARM(results, setname='', test=False, plot_omega=None, total_darm=False):
matplotlib.rcParams.update(style_6p5in) # Add the full width style
fig, axs = plt.subplots(1, 1, sharex=True)
base_F_Hz = results[0]['FOM1_wn_ASD_Hz']
minf = np.min(base_F_Hz)
maxf = np.max(base_F_Hz)
for result in results:
if total_darm:
ASD = result['total_DARM_ASD']
F_Hz = result['total_DARM_ASD_Hz']
else:
ASD = result['FOM1_wn_ASD']*4000
F_Hz = result['FOM1_wn_ASD_Hz']
temp_minf = min(F_Hz)
temp_maxf = max(F_Hz)
if temp_minf<minf:
minf = temp_minf
if temp_maxf>maxf:
maxf = temp_maxf
if result['solver'] == 'Hand Tuned':
l_style = '--'
line_label = 'Hand Tuned'
#plot the reverse BNS Range
BNSfo = result['plant'].siso('FBNS.out', 'FBNS.in').fresponse(F_Hz)
BNS_range_inv = ir.sensemon_range(F_Hz, psd=np.ones_like(F_Hz), m1=1.4, m2=1.4, horizon=False, integrate=False, detection_snr=8)**2
det_ASD_from_FOM = ((BNS_range_inv)/BNSfo.mag**2)**0.5
axs.loglog(F_Hz, det_ASD_from_FOM, linestyle='--', color='black', label='FOM optimization line')
else:
l_style = '-'
line_label = None
axs.loglog(F_Hz, ASD, linestyle=l_style, color=result['linecolor'], label=line_label)
budget_range = gwinc.load_budget('aLIGO', freq=base_F_Hz)
trace_range = budget_range.run(freq=base_F_Hz)
axs.loglog(base_F_Hz, trace_range.psd**0.5*4000, label='aLIGO')
axs.set_xlim([minf, maxf])
axs.legend()
axs.set_xlabel('Frequency [Hz]')
axs.set_ylabel('ASD [m/rtHz]')
norm = results[0]['colormap_norm']
if results[0]['color_param'] == 'F1_gain':
cbar_ticks = None
loop_cbar_label = 'BNS FOM Weight, $\zeta$'
elif results[0]['color_param'] == 'gamma':
cbar_ticks = np.rint(np.geomspace(norm.vmin, norm.vmax, 7)).astype(int)
loop_cbar_label = '$\mathcal{H}_\infty$ bound, $\gamma$'
elif results[0]['color_param'] == 'igsq':
cbar_ticks = None
loop_cbar_label ='Inverse $\mathcal{H}_\infty$ bound squared $(1/\gamma^2)$'
else:
cbar_ticks = None
loop_cbar_label = results[0]['color_param']
cbar = fig.colorbar(cm.ScalarMappable(norm=norm, cmap=results[0]['colormap']), label=loop_cbar_label, ax=axs)
# Adjust the distance between the colorbar and the plot
cbar.ax.yaxis.set_offset_position('right') # To keep the offset on the right side
cbar.ax.yaxis.set_label_position('right')
cbar.ax.set_position([0.76, 0.11, 0.03, 0.77]) # [left, bottom, width, height]
cbar.ax.yaxis.set_label_coords(4, 0.5) # Adjust the position of the label (optional)
if setname != '':
setname_folder = setname + '/'
setname_u = setname + '_'
else:
setname_folder = ''
setname_u = setname + '_'
if not test:
makePlotFolder(setname=setname)
save_fname ='plots/'+ setname_folder + setname_u + 'Noise_Plot_' + results[0]['solver']
fig.savefig(save_fname + '.pdf', bbox_inches='tight', pad_inches=0.05)
fig.savefig(save_fname+ '.png', bbox_inches='tight', pad_inches=0.05)
plt.show()
return fig, axs
[docs]
def plot_CLnoises(results, setname='', extra_fom=None, extra_only=False, scale=1, test=False, plot_omega=None):
if not isinstance(results, list):
results = [results]
ex_fom_list = []
if extra_fom:
if not isinstance(extra_fom, list):
extra_fom = [extra_fom]
for ex_fom in extra_fom:
if ex_fom not in results[0]['CL_all_io'].outputs.keys():
print('Extra FOM, ' + ex_fom + ' not valid output. Skiping this FOM')
else:
ex_fom_list.append(ex_fom)
#print(results[0]['CL_all_io'].outputs)
#print(results[0]['wn_in'])
#print(results[0]['FOM_out'])
if setname != '':
setname_folder = setname + '/'
setname_u = setname + '_'
else:
setname_folder = ''
setname_u = setname + '_'
matplotlib.rcParams.update(style_6p5in) # Add the full width style
fig, axs = plt.subplots(1, 1, sharex=True)
if plot_omega is None:
F_Hz = results[0]['FR_CL_omega_Hz']
else:
F_Hz = plot_omega/(2*np.pi)
# for res in results[0]:
# CL = res['CL_all_io']
# CLNs = CL.siso('F3.out', 'S.in')
# xfs = CLNs.fresponse(F_Hz)
# line, = axs.loglog(*xfs.fplot_mag)
# CLNo = CL.siso('F3.out', 'O.in')
# xfo = CLNo.fresponse(F_Hz)
# line, = axs.loglog(*xfo.fplot_mag, color = line.get_color(), alpha = 0.5)
# rmssq1 = np.sum((xfs.mag)**2 * xfs.df) / (2 * np.pi)
# rmssq2 = np.sum((xfo.mag)**2 * xfs.df) / (2 * np.pi)
# rmssq = rmssq1 + rmssq2
# print("RMS CMP: ", rmssq, CLNs.L2_norm() + CLNo.L2_norm() )
for res in results:
if not extra_only:
FOM_list = res['FOM_out'] + ex_fom_list
else:
if len(ex_fom_list) == 0:
KeyError('No extra FOMs to plot')
else:
FOM_list = ex_fom_list
wn_CL_dict = {}
for fom in FOM_list:
wn_CL_dict[fom] = {}
for wn in res['wn_in']:
CL = res['CL_all_io']
CLNs = CL.siso(fom, wn) * scale
xfs = CLNs.fresponse(f=F_Hz)
wn_CL_dict[fom][wn] = xfs
line, = axs.loglog(*xfs.fplot_mag, label=wn + ' to '+ fom)
xfs_total = np.zeros(len(F_Hz))
for wn in res['wn_in']:
xfs_total += wn_CL_dict[fom][wn].mag**2
line, = axs.loglog(F_Hz, xfs_total**0.5, label='sq all noises to '+ fom)
rmssq1 = np.sum((xfs.mag)**2 * xfs.df) / (2 * np.pi)
rmssq = rmssq1
print("RMS: ", rmssq, CLNs.L2_norm() )
axs.set_xlim([min(F_Hz), max(F_Hz)])
axs.set_ylabel('Magnitude (Abs)')
axs.set_xlabel('Frequency (Hz)')
axs.legend()
if not test:
makePlotFolder(setname=setname)
save_fname ='plots/'+ setname_folder + setname_u + 'Noise_Plot_' + res['solver']
fig.savefig(save_fname + '.pdf', bbox_inches='tight', pad_inches=0.05)
fig.savefig(save_fname+ '.png', bbox_inches='tight', pad_inches=0.05)
plt.show()
for res in results:
if res['solver']=="Hand Tuned" and 'F3.out' in res['CL_all_io'].iod and 'F3.DARM.out' in res['CL_all_io'].iod:
if not test:
makePlotFolder(setname=setname)
save_fname ='plots/'+ setname_folder + setname_u + 'Noise_Plot_hand_' + res['solver']
fig.savefig(save_fname + '.pdf', bbox_inches='tight', pad_inches=0.05)
fig.savefig(save_fname+ '.png', bbox_inches='tight', pad_inches=0.05)
oml=np.array([1e1,1e3])*2*np.pi
oml=None
tmp_bode = ssutil.bode(res['CL_all_io'].siso('F3.out', 'S.in'), omega_limits=oml, label='f3 env noise')
tmp_bode = ssutil.bode(res['CL_all_io'].siso('F3.out', 'O.in'), ax=tmp_bode, omega_limits=oml, label='f3 meas noise')
tmp_bode = ssutil.bode(res['CL_all_io'].siso('F3.DARM.out', 'S.in'), ax=tmp_bode, omega_limits=oml, label='Darm env noise')
tmp_bode = ssutil.bode(res['CL_all_io'].siso('F3.DARM.out', 'O.in'), ax=tmp_bode, omega_limits=oml, label='Darm meas noise')
if not test:
tmp_bode.save('plots/'+ setname_folder + setname_u + 'tmp_Bode.pdf')
return fig, axs
"""def plotFOM(result, setname='', test=False, plot_omega=None, from_plant=True):
#matplotlib.rcParams.update(matplotlib.rcParamsDefault)
matplotlib.rcParams.update(style_3p5in)
if type(result) != dict:
if type(result) == list:
print('result must be a dictionary. Given list. Using first entry')
return plotFOM(result[0], setname=setname, test=test, plot_omega=plot_omega)
else:
print('result must be a dictionary. or list of dictionaries where the first result will be used')
if setname != '':
setname_folder = setname + '/'
setname_u = setname + '_'
else:
setname_folder = ''
if not plot_omega:
plot_omega = result['plot_omega']
FR_FOM ={}
for fom in result['FOM_out']:
fom_plot_in = result['control_in'][0]
fom_plot_out = fom
if from_plant:
fom_plot_siso = ssutil.truncate_io(result['plant'], [fom_plot_in], [fom_plot_out])
else:
try:
from T_Make_ASC_Model import FBNSsimpSS, FFlatSS
if fom == 'FBNS.out':
fom_plot_siso_mod = FBNSsimpSS(gain=50**2)
elif fom == 'F3.out':
fom_plot_siso_mod = FFlatSS(gain=1)
else:
KeyError('FOM_out must be FBNS.out or F3.out')
except:
print('Could not import FBNSsimpSS or FFlatSS. Plotting FOM from plant')
return plotFOM(result, setname=setname, test=test, plot_omega=plot_omega, from_plant=True)
if fom == 'FBNS.out':
name = 'LIGO\'s BNS\nSensitivity FOM'
plotcolor ='r'
elif fom == 'F3.out':
name='Broadband FOM'
plotcolor ='dodgerblue'
else:
name = fom.split('.')[0]
fom_plot_mag, fom_plot_ph, fom_plot_omega = control.freqresp(fom_plot_siso.mod, plot_omega)
wndict = {'mag': fom_plot_mag,
'phase': (fom_plot_ph + np.pi) % (2 * np.pi) - np.pi,
'phase_unwrapped': fom_plot_ph,
'omega': fom_plot_omega,
'omega_Hz': fom_plot_omega / (2 * np.pi),
'name': name,
'color': plotcolor,
}
FR_FOM[fom] = wndict
fig, axs = plt.subplots()
axs.set_ylabel('Magnitude (Abs)')
axs.set_xlabel('Frequency (Hz)')
axs.grid()
for key in FR_FOM.keys():
axs.plot(FR_FOM[key]['omega_Hz'], FR_FOM[key]['mag'], label=FR_FOM[key]['name'], color=FR_FOM[key]['color'])
#axs[1].plot(FR_FOM[key]['omega'], FR_FOM[key]['phase'])
axs.set_yscale('log')
axs.set_xscale('log')
axs.set_xlim([min(FR_FOM[key]['omega_Hz']), max(FR_FOM[key]['omega_Hz'])])
axs.legend(loc='lower right')
if not test:
makePlotFolder(setname=setname)
save_fname ='plots/'+ setname_folder + setname_u + 'FOM_Plot'
fig.savefig(save_fname + '.pdf', bbox_inches='tight', pad_inches=0.05)
fig.savefig(save_fname+ '.png', bbox_inches='tight', pad_inches=0.05, dpi=500)
plt.show()
return fig"""
[docs]
def plot_gm_pm(HinfBounded_results, setname='', f1line=False, text=False, xlim=None, ylim=None, test=True, io_scales=None):
# Create a custom colormap for the legend
num_points = 20 # You can adjust this to control the number of colors in the legend point
gradient = np.linspace(0, 1, num_points).reshape(-1, 1)
gradient = np.vstack((gradient, gradient))
extraplots = buzzutil.filtparam(HinfBounded_results, 'solver', 'Hand Tuned', compare='eq')
#HinfBounded_results_E_gamma_clr = compute.calcFullResults(HinfBounded_results, colormap='RdYlGn', plot_omega=np.logspace(-2,3,1000)*2*np.pi, color_param='gamma', label=False)
#plot_results_no_H2_gamma_gm_pm = buzzutil.filtparam(HinfBounded_results_E_gamma_clr, 'solver', 'Hand Tuned', compare='neq')
plot_results_no_H2_gamma_gm_pm = buzzutil.filtparam(HinfBounded_results, 'gamma', np.inf, compare='neq')
plot_results_no_H2_gamma_gm_pm = compute.calcFullResults(plot_results_no_H2_gamma_gm_pm, colormap='RdYlGn', color_param='gamma', label=False, io_scales=io_scales)
fig, ax = plt.subplots()
cmap = plt.get_cmap(plot_results_no_H2_gamma_gm_pm[2]['colormap'])
cmap.set_bad('lightgrey')
cmap_norm = plot_results_no_H2_gamma_gm_pm[2]['colormap_norm']
legend_cmap = LinearSegmentedColormap.from_list('legend_cmap', cmap(gradient[:, 0]))
pm_list = np.abs(buzzutil.listparams(plot_results_no_H2_gamma_gm_pm, 'pm'))
gm_list = np.abs(buzzutil.listparams(plot_results_no_H2_gamma_gm_pm, 'gm'))
clr_list = buzzutil.listparams(plot_results_no_H2_gamma_gm_pm, 'linecolor')
ax.scatter(pm_list, gm_list, c=clr_list, s=30)#, edgecolors='black')
try:
extraplots_gm_pm = extraplots
extraplots_gm_pm[0]['linecolor'] = 'DodgerBlue'
ax.scatter(np.abs(extraplots_gm_pm[0]['pm']), np.abs(extraplots_gm_pm[0]['gm']), c=extraplots_gm_pm[0]['linecolor'], s=80, label=extraplots_gm_pm[0]['label1'], marker=extraplots_gm_pm[0]['rms_marker'], edgecolors='none')
except:
extraplots_gm_pm = None
if f1line:
for gain in np.unique(buzzutil.listparams(HinfBounded_results, 'F1_gain')):
results_line = buzzutil.filtparam(HinfBounded_results, 'F1_gain', [gain])
pm_line = []
gm_line = []
for result in results_line:
pm_line.append(abs(result['pm']))
gm_line.append(abs(result['gm']))
if gain == np.unique(buzzutil.listparams(HinfBounded_results, 'F1_gain'))[-1]:
f1gain_lab = 'Lines of Constant $\zeta$'
else:
f1gain_lab = None
ax.plot(pm_line, gm_line, color='black', linewidth=1.4, linestyle='dotted', alpha=0.35, label=f1gain_lab)
if text:
if isinstance(text, str):
if text == 'idx':
text_vals_non_int =np.linspace(0,len(HinfBounded_results)-1,len(HinfBounded_results))
text_vals = np.array([int(val) for val in text_vals_non_int])
else:
try:
text_vals = np.abs(buzzutil.listparams(HinfBounded_results, text))
except:
KeyError('Text parameter not found in results')
else:
TypeError('Text parameter must be a string or None')
for i, txt in enumerate(text_vals):
ax.annotate(txt, (np.abs(HinfBounded_results[i]['pm']), np.abs(HinfBounded_results[i]['gm'])), color='black',xytext=(5, 5), textcoords="offset pixels", path_effects=[pe.withStroke(linewidth=1.5, foreground='w')])
if text == 'gm':
text_ann = 'Abs Gain Margin (dB) (Text Near Point)'
elif text == 'pm':
text_ann = 'Phase Margin (deg) (Text Near Point)'
elif text == 'idx':
text_ann = 'index of controller in list (Text Near Point)'
else:
text_ann = text + ' (Text Near Point)'
ax.text(1-0.014,0.01,text_ann, horizontalalignment='right',verticalalignment='bottom', transform=ax.transAxes)
# Create a custom legend
legend_handles = []
legend_labels = ["$\mathcal{H}_\infty$ Bounded Controllers", "Hand Tuned Controller"]
if extraplots_gm_pm:
# Define the number of dots in the legend point
num_dots = 5
# Define the colors for the legend dots based on the gradient
legend_dot_colors = legend_cmap(np.linspace(0.99, 0.51, num_dots))
# Create a custom handler to create the legend point with multiple dots
class MultiDotHandler(object):
def legend_artist(self, legend, orig_handle, fontsize, handlebox):
move_dots_right = 3.3
x0, y0 = handlebox.xdescent+move_dots_right, handlebox.ydescent
width, height = handlebox.width, handlebox.height
for i in range(num_dots):
dot_x = x0 + i * width / num_dots
dot_y = y0 + height / 2.0
legend_marker = Line2D([dot_x], [dot_y], marker='o', markersize=8, color=legend_dot_colors[i])
handlebox.add_artist(legend_marker)
# Create a custom handler for the second entry (blue 'x' marker)
class BlueXHandler(object):
def legend_artist(self, legend, orig_handle, fontsize, handlebox):
move_dots_right = 3.3
x0, y0 = handlebox.xdescent, handlebox.ydescent
width, height = handlebox.width, handlebox.height
#x_marker = Line2D([10], [3], markersize=6.5, markerfacecolor=extraplots_gm_pm[0]['linecolor'], marker=extraplots_gm_pm[0]['rms_marker'], markeredgecolor=extraplots_gm_pm[0]['linecolor'], markeredgewidth=0.0)
x_marker = Line2D([10], [3], markersize=8, c=extraplots_gm_pm[0]['linecolor'], marker=extraplots_gm_pm[0]['rms_marker'], markeredgewidth=0.0)
#x_marker = plt.scatter([x0 + width/2], [y0 + height/2], c=extraplots_gm_pm[0]['linecolor'], marker=extraplots_gm_pm[0]['rms_marker'], edgecolors=extraplots_gm_pm[0]['linecolor'])
handlebox.add_artist(x_marker)
# Create a new legend handle for the first scatter plot
first_handle = Line2D([0], [0])
# Add the new handle and label to the existing lists
legend_handles.append(first_handle)
# Create a new legend handle for the second scatter plot
second_handle = Line2D([0], [0])
# Add the new handle and label to the existing lists
legend_handles.append(second_handle)
# Create a handler map for the custom legend handlers
handler_map = {legend_handles[0]: MultiDotHandler(), legend_handles[1]: BlueXHandler()}
# Add the legend to the plot
ax.legend(legend_handles, legend_labels, handler_map=handler_map, loc='upper left')
if xlim:
ax.set_xlim(xlim)
if ylim:
ax.set_ylim(ylim)
#norm1 = matplotlib.colors.Normalize(vmin=np.min(outlinecolor), vmax=np.max(outlinecolor))
#cb1 = matplotlib.colorbar.ColorbarBase(axs1, cmap=outline_cmap, norm=norm1)
#cb1 = matplotlib.colorbar.ColorbarBase(ax2, cmap=cmap, norm=cmap_norm)
norm = plot_results_no_H2_gamma_gm_pm[0]['colormap_norm']
cbar_ticks = np.rint(np.geomspace(norm.vmin, norm.vmax, 7)).astype(int)
cbar_tick_labels = cbar_ticks #np.round(cbar_ticks, 0)
loop_cbar_label = '$\mathcal{H}_\infty$ bound, $\gamma$'
cbar = fig.colorbar(cm.ScalarMappable(norm=norm, cmap=plot_results_no_H2_gamma_gm_pm[0]['colormap']), ax=ax, label=loop_cbar_label, ticks=cbar_ticks)
cbar.ax.set_yticklabels(cbar_tick_labels)
cbar.ax.yaxis.set_offset_position('right') # To keep the offset on the right side
cbar.ax.yaxis.set_label_position('right')
cbar.ax.set_position([0.76, 0.11, 0.03, 0.77]) # [left, bottom, width, height]
cbar.ax.yaxis.set_label_coords(3, 0.5) # Adjust the position of the label (optional)
#cbar = fig.colorbar(cm.ScalarMappable(norm=cmap_norm, cmap=cmap), label='Gamma', ticks=[1.5, 2,3,4,5,6,8,10,13, 30, 100], ax=ax) # Add a colorbar to show the color scale
#cbar.ax.set_yticklabels(['1.5','2','3','4', '5','6','8', '10','13','30', '100'])
ax.set_xlabel('Abs. Val. Phase Margin (Degrees)')
ax.set_ylabel('Abs. Val. Gain Margin (dB)')
setname_folder = setname + '/'
setname_u = setname + '_'
if not test:
save_fname ='plots/'+ setname_folder + setname_u + 'gmpm_Plot'
fig.savefig(save_fname + '.pdf', bbox_inches='tight', pad_inches=0.05)
fig.savefig(save_fname+ '.png', bbox_inches='tight', pad_inches=0.05, dpi=500)
plt.show()
return fig, ax
[docs]
def vega_plot(results, setname='', test=False, h2line=True, size=600, plotrange=False):
import altair as alt
import polars as pl
H2_results = buzzutil.filtparam(results, 'solver', 'LQG', compare='equal')
results_HiB = buzzutil.filtparam(results, 'solver', ['HB','mix','BH', 'HiB'], compare='equal')
def keyfilt(d, **kw):
if isinstance(d, (list, tuple)):
return [keyfilt_(v, **kw) for v in d]
else:
return keyfilt_(d, **kw)
def keyfilt_(d):
d2 = {}
for k, v in d.items():
v = np.asarray(v)
dt = v.dtype
if v.shape != () and np.issubdtype(dt, np.number):
d2[k] = list(v)
elif v.shape != ():
pass
# d2[k] = list(v)
elif dt == object:
pass
# d2[k] = v
else:
d2[k] = v.item()
return d2
#######################################################
# df_H2 = pl.DataFrame([H2_results[10]]).with_row_index()
df_HiB = pl.DataFrame(keyfilt(results_HiB[5:10])).with_row_index()
#df_HiB = df_HiB.to_pandas()
###########################################################
# reduced dataset for debugging
#df_HiB = pl.DataFrame(keyfilt(results_HiB[20:25])).with_row_index()
# full data for plotting
#df_HiB = pl.DataFrame(keyfilt(results_HiB + H2_results)).with_row_index()
df_HiB = pl.DataFrame(keyfilt(results_HiB + H2_results)).with_row_index()
df_HiB = df_HiB.with_columns(pm = abs(pl.col("pm")))
df_ids = df_HiB.select('index')
# df_H2 = pl.DataFrame(keyfilt(results_HiB)).with_row_index()
# df_H2 = df_H2.with_columns(pm = abs(pl.col("pm")))
tooltip_info = [
'index',
'solver',
'F1_gain',
'F2_gain',
'gamma',
'igsq',
'isoptimal',
'F1_H2_Norm',
'F2_H2_Norm',
'isKstable',
'gm',
'pm',
'wcg',
'wcp',
]
alt_tooltip_info = [
'index:O',
'solver:N',
'F1_gain:Q',
'F2_gain:Q',
'gamma:Q',
'igsq:Q',
'isoptimal:N',
'F1_H2_Norm:Q',
'F2_H2_Norm:Q',
'isKstable:N',
'gm:Q',
'pm:Q',
'wcg:Q',
'wcp:Q',
]
#df_HiB
###########################################################
selector = alt.selection_point(fields=['index'])
if plotrange:
F1_plot_param = 'range'
else:
F1_plot_param = 'F1_H2_Norm'
alt_HiB = alt.LookupData(df_HiB, 'index', list(set([F1_plot_param, 'F2_H2_Norm', 'pm']) | set(tooltip_info)))
baseheight = size
magratio = 0.66
bodemagheight = int(baseheight*magratio)
bodephaseheight = baseheight - bodemagheight
basewidth = size
base = alt.Chart(df_ids).add_params(selector)
main_color = dict(
color = alt.condition(
selector,
alt.Color('index:N').legend(None),
alt.value('lightgray')
))
main_color2 = dict(
color = alt.Color('index:N'),
opacity = alt.condition(
selector,
alt.value(1),
alt.value(0),
))
# F1_line_color = dict(
# color=alt.condition(
# alt.selection_single(),
# alt.Color('f1_gain:N').legend(None),
# alt.value('lightgray')
# )
# )
if plotrange:
xtitle = 'Range (Mpc)'
else:
xtitle = "RMS of 8-30 bandpass, control counts"
#autofind domain with margins of 10%
domain_margin = 0.1
domain_dif = max(df_HiB[F1_plot_param]) - min(df_HiB[F1_plot_param])
custom_domain = [min(df_HiB[F1_plot_param])-domain_dif*domain_margin, max(df_HiB[F1_plot_param])+domain_dif*domain_margin]
rms_points = base.transform_lookup(lookup='index', from_=alt_HiB).mark_point(filled=True, size=200).encode(
x=alt.X(F1_plot_param+':Q').scale(type='log').title(xtitle),
y=alt.Y('F2_H2_Norm:Q').scale(type='log').title("RMS of Flat FOM, error counts"),
tooltip=alt_tooltip_info,
**main_color,
).properties(
width=basewidth,
height= baseheight,
)
rms_points_cl = base.transform_lookup(lookup='index', from_=alt_HiB).mark_point(size=float('NaN')).encode(
**main_color,
)
rms_pm = base.transform_lookup(lookup='index', from_=alt_HiB).mark_point(filled=False, size=350).encode(
#x=alt.X(F1_plot_param+':Q').scale(type='log', domain=custom_domain),
x=alt.X(F1_plot_param+':Q').scale(type='log'),
y=alt.Y('F2_H2_Norm:Q').scale(type='log'),
color=alt.condition(
selector,
alt.Color('pm:Q').scale(scheme="redblue", reverse=False),
alt.value('lightgray')
),
tooltip=alt_tooltip_info,
).properties(
width=basewidth,
height= baseheight,
)
# Define the lines connecting points with the same phase margin (pm)
H2_line = base.transform_lookup(lookup='index', from_=alt_HiB).transform_filter(
alt.datum.solver == 'LQG').mark_line(opacity=0.4).encode(
#x=alt.X(F1_plot_param + ':Q').scale(type='log', domain=custom_domain).title(xtitle),
x=alt.X(F1_plot_param + ':Q').scale(type='log').title(xtitle),
y=alt.Y('F2_H2_Norm:Q').scale(type='log'),
color=alt.value('red'),
detail='solver:N', # Group lines by solver
tooltip=alt_tooltip_info,
).properties(
width=basewidth,
height=baseheight,
)
# Define the lines connecting points with the same phase margin (pm)
F1_gain_lines = base.transform_lookup(lookup='index', from_=alt_HiB).transform_filter(
alt.datum.solver == 'HB').mark_line(strokeDash=[5, 5], opacity=0.2).encode(
#x=alt.X(F1_plot_param + ':Q').scale(type='log', domain=custom_domain).title(xtitle),
x=alt.X(F1_plot_param + ':Q').scale(type='log').title(xtitle),
y=alt.Y('F2_H2_Norm:Q').scale(type='log'),
color=alt.value('black'),
detail='F1_gain:Q', # Group lines by f1 gain
tooltip=alt_tooltip_info,
).properties(
width=basewidth,
height=baseheight,
)
# because the dataset is so large, we are using tricks to synthisize subsets using the lookup join
alt_HiB_OL = alt.LookupData(df_HiB, 'index', [F1_plot_param, 'F2_H2_Norm', 'pm', 'FR_OL_omega_Hz', 'FR_OL_mag', 'FR_OL_ph_deg'])
base_OL = base.transform_lookup(lookup='index', from_=alt_HiB_OL).transform_flatten(['FR_OL_omega_Hz', 'FR_OL_mag', 'FR_OL_ph_deg'])
alt_HiB_CL = alt.LookupData(df_HiB, 'index', [F1_plot_param, 'F2_H2_Norm', 'pm', 'FR_CL_omega_Hz', 'FR_CL_mag', 'FR_CL_ph_deg'])
base_CL = base.transform_lookup(lookup='index', from_=alt_HiB_CL).transform_flatten(['FR_CL_omega_Hz', 'FR_CL_mag', 'FR_CL_ph_deg'])
OL_mag = base_OL.mark_line().encode(
x=alt.X('FR_OL_omega_Hz:Q').scale(type='log'),
y=alt.Y('FR_OL_mag:Q').scale(type='log'),
**main_color,
).properties(
width=basewidth,
height=bodemagheight
).transform_filter(selector)
OL_phase = base_OL.mark_line().encode(
x=alt.X('FR_OL_omega_Hz:Q').scale(type='log'),
y=alt.Y('FR_OL_ph_deg:Q'),
**main_color,
).properties(
width=basewidth,
height=bodephaseheight
).transform_filter(selector)
CL_mag = base_CL.mark_line().encode(
x=alt.X('FR_CL_omega_Hz:Q').scale(type='log'),
y=alt.Y('FR_CL_mag:Q').scale(type='log'),
**main_color,
).properties(
width=basewidth,
height=bodemagheight
).transform_filter(selector)
CL_phase = base_CL.mark_line().encode(
x=alt.X('FR_CL_omega_Hz:Q').scale(type='log'),
y=alt.Y('FR_CL_ph_deg:Q'),
**main_color,
).properties(
width=basewidth,
height=bodephaseheight
).transform_filter(selector)
if h2line and len(H2_results)!=0:
h2line = base.mark_line().encode(
x=alt.X(F1_plot_param+':Q').scale(type='log'),
y=alt.Y('F2_H2_Norm:Q').scale(type='log'),
**main_color,
).properties(
width=basewidth,
height=baseheight,
).transform_filter(selector)
# the inclusion of rms_points_cl is a hack to have the line plots use consistent coloring with the rms_points plot. It normally would, but the inclusion of rms_pm which has a different color map screws it up
fullchart = (
(alt.layer( H2_line, F1_gain_lines, rms_pm, rms_points, title = 'RMS') ) | alt.vconcat(rms_points_cl + OL_mag, OL_phase, title='Open Loop Response') | alt.vconcat(CL_mag, CL_phase, title='Closed Loop Response')
).interactive()
if not test:
makePlotFolder(setname=setname)
setname_folder = setname + '/'
setname_u = setname + '_'
save_fname ='plots/'+ setname_folder + setname_u + '_plot'
fullchart.save(save_fname + '.html')
fullchart.save('plots/vega.html')
fullchart.show()
return fullchart
[docs]
def MIMObode(sys, plotting_omega=None, omega_limits=None, UG_line=False, grid=True, **kwargs):
"""
Create a grid of Bode plots for a MIMO system.
Parameters:
sys (control.StateSpace or control.TransferFunction): The MIMO system.
plotting_omega (np.array): The frequency range to plot. If None, it will be calculated automatically.
omega_limits (list): The limits of the frequency range to plot. If None, it will be calculated automatically.
UG_line (bool): If True, a line will be plotted at unity gain.
grid (bool): If True, a grid will be plotted.
**kwargs: Additional keyword arguments to pass to control.bode.
Returns:
fig (matplotlib.figure.Figure): The figure object.
axes (np.array): The axes objects.
"""
from buzz import iodutil
print('MIMObode is a work in progress and may not work as expected.')
print(' for systems with many inputs and outputs, the plot take a while and may be too large to display.')
outputs, inputs = iodutil.listoutputs(sys), iodutil.listinputs(sys)
fig, axes = plt.subplots(len(outputs)*2, len(inputs), figsize=(3*len(inputs), 2*len(outputs)), sharex=True)
fig.suptitle('MIMO Bode Plot')
for i in range(len(outputs)):
for j in range(len(inputs)):
ax_single_bode = Bunch()
ax_single_bode.ax0 = axes[i*2, j]
ax_single_bode.ax1 = axes[i*2+1, j]
if i == len(outputs)-1:
xlabel=True
else:
xlabel=False
if j == 0:
ylabel=True
else:
ylabel=False
_ = ssutil.bode(sys.siso(outputs[i], inputs[j]), ax=ax_single_bode, plotting_omega=plotting_omega, omega_limits=omega_limits, UG_line=UG_line, axB=None, xlabel=xlabel, ylabel=ylabel, grid=grid, **kwargs)
#axes[i, j].set_title(f'Output {i+1} to Input {j+1}')
plt.subplots_adjust(hspace=0.1, wspace=0.1)
plt.show()
return fig, axes