diff --git a/pes_to_spec/test/offline_analysis.py b/pes_to_spec/test/offline_analysis.py
index 871670e92e0c8b13091b1fa6cbf5128656a592d2..cf5e482491e8c5432a7599b8bdabd0c0bad3441a 100755
--- a/pes_to_spec/test/offline_analysis.py
+++ b/pes_to_spec/test/offline_analysis.py
@@ -126,7 +126,8 @@ def plot_result(filename: str,
#spec_raw_int: Optional[np.ndarray]=None,
pes: Optional[np.ndarray]=None,
pes_to_show: Optional[str]="",
- pes_bin: Optional[np.ndarray]=None,
+ first: Optional[int]=None,
+ last: Optional[int]=None,
):
"""
Plot result with uncertainty band.
@@ -139,7 +140,6 @@ def plot_result(filename: str,
spec_raw_int: Original true expected result with shape (features,).
pes: PES spectrum for the inset.
pes_to_show: Name of the channel shown.
- pes_bin: PES bins.
intensity: The XGM intensity in uJ.
"""
@@ -150,22 +150,27 @@ def plot_result(filename: str,
spec=spec_smooth,
prediction=spec_pred["expected"],
unc=unc,
+ unc_pca=unc_pca,
+ unc_stat=unc_stat,
beam_intensity=intensity*1e-3*np.ones_like(spec_raw_pe)
))
df.to_csv(filename.replace('.pdf', '.csv'))
if pes is not None:
pes_data = deepcopy(pes)
pes_data['bin'] = np.arange(len(pes['channel_1_D']))
+ pes_data['first'] = first*np.ones_like(pes_data['bin'])
+ pes_data['last'] = last*np.ones_like(pes_data['bin'])
df = pd.DataFrame(pes_data)
df.to_csv(filename.replace('.pdf', '_pes.csv'))
fig = plt.figure(figsize=(12, 8))
gs = GridSpec(1, 1)
ax = fig.add_subplot(gs[0, 0])
- ax.plot(spec_raw_pe, spec_smooth, c='b', lw=3, label="High-res. measurement (smoothened)")
- ax.plot(spec_raw_pe, spec_pred["expected"], c='r', ls='--', lw=3, label="High-res. prediction")
+ ax.plot(spec_raw_pe, spec_smooth, c='b', lw=3, label="Grating spectrometer")
+ ax.plot(spec_raw_pe, spec_pred["expected"], c='r', ls='--', lw=3, label="Prediction")
#ax.fill_between(spec_raw_pe, spec_pred["expected"] - unc, spec_pred["expected"] + unc, facecolor='green', alpha=0.6, label="68% unc.")
- ax.fill_between(spec_raw_pe, spec_pred["expected"] - unc, spec_pred["expected"] + unc, facecolor='gold', alpha=0.5, label="68% unc.")
+ ax.fill_between(spec_raw_pe, spec_pred["expected"] - unc, spec_pred["expected"] + unc, facecolor='gold', alpha=0.5, label="68% unc. (total)")
+ ax.fill_between(spec_raw_pe, spec_pred["expected"] - unc_pca, spec_pred["expected"] + unc_pca, facecolor='magenta', alpha=0.5, label="68% unc. (PCA only)")
#ax.fill_between(spec_raw_pe, spec_pred["expected"] - unc_stat, spec_pred["expected"] + unc_stat, facecolor='red', alpha=0.6, label="68% unc. (stat.)")
#ax.fill_between(spec_raw_pe, spec_pred["expected"] - unc_pca, spec_pred["expected"] + unc_pca, facecolor='magenta', alpha=0.6, label="68% unc. (syst., PCA)")
#if spec_raw_int is not None:
@@ -188,6 +193,7 @@ def plot_result(filename: str,
#ip = InsetPosition(ax, [0.65,0.6,0.35,0.4])
ip = InsetPosition(ax, [0.72,0.7,0.35,0.4])
ax2.set_axes_locator(ip)
+ pes_bin = np.arange(first, last)
if pes_to_show == "sum":
pes_plot = sum([pes[k][pes_bin] for k in pes.keys()])
pes_label = r"$\sum$ PES channels"
@@ -366,12 +372,20 @@ def main():
pca_variance_plot(os.path.join(args.directory, f"pca_pes.pdf"),
pca.explained_variance_ratio_,
600)
+ df = pd.DataFrame(dict(variance_ratio=pca.explained_variance_ratio_,
+ n_comp=600*np.ones_like(pca.explained_variance_ratio_),
+ ))
+ df.to_csv(os.path.join(args.directory, "pca_pes.csv"))
pca_spec = PCA(None, whiten=True)
pca_spec.fit(spec_raw_int[train_idx])
pca_variance_plot(os.path.join(args.directory, f"pca_spec.pdf"),
pca_spec.explained_variance_ratio_,
20)
+ df = pd.DataFrame(dict(variance_ratio=pca_spec.explained_variance_ratio_,
+ n_comp=20*np.ones_like(pca_spec.explained_variance_ratio_),
+ ))
+ df.to_csv(os.path.join(args.directory, "pca_spec.csv"))
model.debug_peak_finding(pes_raw, os.path.join(args.directory, "test_peak_finding.pdf"))
if len(args.model) == 0:
@@ -486,7 +500,7 @@ def main():
ax.scatter(chi2/ndof, xgm_flux_t[:,0], c='r', s=20)
ax.set(title=f"", #avg(stat unc) = {unc_stat}, avg(pca unc) = {unc_pca}",
xlabel=r"$\chi^2/$ndof",
- ylabel="XGM intensity [uJ]",
+ ylabel="Beam intensity [uJ]",
xlim=(0, 5),
)
ax2 = plt.axes([0,0,1,1])
@@ -498,7 +512,7 @@ def main():
#ax2.scatter(chi2/ndof, np.sum(spec_raw_int, axis=1)*de, c='g', s=30)
ax2.set(title="",
xlabel=r"$\chi^2/$ndof",
- ylabel=f"XGM intensity [uJ]",
+ ylabel=f"Beam intensity [uJ]",
)
ax2.title.set_size(SMALL_SIZE)
ax2.xaxis.label.set_size(SMALL_SIZE)
@@ -529,6 +543,7 @@ def main():
spec_smooth_pca = model.y_model['pca'].transform(spec_smooth)
chi2_prepca = np.sum((spec_smooth_pca[:, np.newaxis, :] - spec_pred["expected_pca"])**2/(spec_pred["expected_pca_unc"]**2), axis=(-1, -2))
+
ndof_prepca = float(spec_smooth_pca.shape[-1])
print(f"Chi2 before PCA: {np.mean(chi2_prepca):.2f}, ndof: {ndof_prepca}, chi2/ndof: {np.mean(chi2_prepca/ndof_prepca):.2f} +/- {np.std(chi2_prepca/ndof_prepca):.2f}")
fig = plt.figure(figsize=(12, 8))
@@ -536,7 +551,7 @@ def main():
ax = fig.add_subplot(gs[0, 0])
sns.histplot(x=chi2_prepca/ndof_prepca, kde=True, linewidth=3, ax=ax)
ax.set(title=f"",
- xlabel=r"$\chi^2/$ndof before undoing PCA",
+ xlabel=r"$\chi^2/$ndof",
ylabel="Counts [a.u.]",
xlim=(0, 5),
)
@@ -556,8 +571,8 @@ def main():
ax = fig.add_subplot(gs[0, 0])
ax.scatter(chi2_prepca/ndof_prepca, xgm_flux_t[:,0], c='r', s=20)
ax.set(title=f"",
- xlabel=r"$\chi^2/$ndof before undoing PCA",
- ylabel="XGM intensity [uJ]",
+ xlabel=r"$\chi^2/$ndof",
+ ylabel="Beam intensity [uJ]",
xlim=(0, 5),
ylim=(0, np.mean(xgm_flux_t) + 3*np.std(xgm_flux_t))
)
@@ -586,7 +601,7 @@ def main():
ax = fig.add_subplot(gs[0, 0])
sns.histplot(x=xgm_flux_t[:,0], kde=True, linewidth=3, ax=ax)
ax.set(title=f"",
- xlabel="XGM intensity [uJ]",
+ xlabel="Beam intensity [uJ]",
ylabel="Counts [a.u.]",
)
#ax.text(0.90, 0.95, fr"$\mu = ${np.mean(xgm_flux_t[:,0]):.2f}",
@@ -610,7 +625,7 @@ def main():
ax = plt.gca()
ax.set(title=f"",
xlabel=r"Root-mean-squared error",
- ylabel="XGM intensity [uJ]",
+ ylabel="Beam intensity [uJ]",
)
fig.savefig(os.path.join(args.directory, "intensity_vs_rmse.pdf"))
plt.close(fig)
@@ -669,6 +684,13 @@ def main():
#fig.savefig(os.path.join(args.directory, "xgm_vs_expected.png"))
#plt.close(fig)
+ df = pd.DataFrame(dict(chi2_prepca=chi2_prepca,
+ ndof=spec_smooth_pca.shape[-1]*np.ones_like(chi2_prepca),
+ xgm_flux_t=xgm_flux_t[:,0],
+ rmse=rmse,
+ ))
+ df.to_csv(os.path.join(args.directory, "quality.csv"))
+
first, last = model.get_low_resolution_range()
first = max(0, first+200)
last = min(last-200, pes_raw_t["channel_1_D"].shape[1]-1)
@@ -690,7 +712,8 @@ def main():
pes={k: -item[idx, :]
for k, item in pes_raw_t.items()},
pes_to_show=pes_to_show,
- pes_bin=np.arange(first, last),
+ first=first,
+ last=last,
)
for ch in channels:
plot_pes(os.path.join(args.directory, f"test_pes_{tid}_{ch}.pdf"),
diff --git a/pes_to_spec/test/prepare_plots.py b/pes_to_spec/test/prepare_plots.py
index a348a88b2afdf339baad6f3a24fb88824b1f4497..50343609101884fe255ef660dafa43dc3f53f64b 100755
--- a/pes_to_spec/test/prepare_plots.py
+++ b/pes_to_spec/test/prepare_plots.py
@@ -6,6 +6,18 @@ import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import seaborn as sns
+SMALL_SIZE = 12
+MEDIUM_SIZE = 18
+BIGGER_SIZE = 22
+
+plt.rc('font', size=BIGGER_SIZE) # controls default text sizes
+plt.rc('axes', titlesize=BIGGER_SIZE) # fontsize of the axes title
+plt.rc('axes', labelsize=BIGGER_SIZE) # fontsize of the x and y labels
+plt.rc('xtick', labelsize=BIGGER_SIZE) # fontsize of the tick labels
+plt.rc('ytick', labelsize=BIGGER_SIZE) # fontsize of the tick labels
+plt.rc('legend', fontsize=MEDIUM_SIZE) # legend fontsize
+plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title
+
def plot_final(df: pd.DataFrame, filename: str):
fig = plt.figure(figsize=(12, 8))
gs = GridSpec(1, 1)
@@ -16,7 +28,7 @@ def plot_final(df: pd.DataFrame, filename: str):
ax.fill_between(df.energy, df.prediction - df.unc_pca, df.prediction + df.unc_pca, facecolor='magenta', alpha=0.5, label="68% unc. (PCA only)")
Y = np.amax(df.spec)
ax.legend(frameon=False, borderaxespad=0, loc='upper left')
- ax.set_title(f"Beam intensity: {df.beam_intensity.iloc[0]*1e-3:.1f} mJ", loc="left")
+ ax.set_title(f"Beam intensity: {df.beam_intensity.iloc[0]:.1f} mJ", loc="left")
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.set(
@@ -39,6 +51,18 @@ def plot_chi2(df: pd.DataFrame, filename: str):
fig.savefig(filename)
plt.close(fig)
+def plot_rmse(df: pd.DataFrame, filename: str):
+ fig = plt.figure(figsize=(12, 8))
+ gs = GridSpec(1, 1)
+ ax = fig.add_subplot(gs[0, 0])
+ sns.histplot(x=df.rmse, kde=True, linewidth=3, ax=ax)
+ ax.set(title=f"",
+ xlabel=r"Root-mean-square-error",
+ ylabel="Counts [a.u.]",
+ )
+ fig.savefig(filename)
+ plt.close(fig)
+
def plot_rmse_intensity(df: pd.DataFrame, filename: str):
fig = plt.figure(figsize=(12, 8))
gs = GridSpec(1, 1)
@@ -113,7 +137,8 @@ def plot_pes(df: pd.DataFrame, channel:str, filename: str):
fig = plt.figure(figsize=(16, 8))
gs = GridSpec(1, 1)
ax = fig.add_subplot(gs[0, 0])
- ax.plot(df.bin[df.range], df.loc[:, channel], c='b', lw=3)
+ first, last = df.loc[:, 'first'].iloc[0], df.loc[:, 'last'].iloc[0]
+ ax.plot(df.loc[(df.bin >= first) & (df.bin < last), "bin"], df.loc[(df.bin >= first) & (df.bin < last), channel], c='b', lw=3)
#ax.legend()
ax.set(title=f"",
xlabel="Time-of-flight index",
@@ -124,19 +149,21 @@ def plot_pes(df: pd.DataFrame, channel:str, filename: str):
plt.close(fig)
if __name__ == '__main__':
+ indir = 'p900331r69t70'
channel = 'channel_4_A'
fname = 'test_q100_1724098413'
- plot_final(pd.read_csv(f'{fname}.csv'), f'{fname}.pdf')
- plot_pes(pd.read_csv(f'{fname}_pes.csv'), channel, f'{fname}_{channel}.pdf')
+ plot_final(pd.read_csv(f'{indir}/{fname}.csv'), f'{fname}.pdf')
+ plot_pes(pd.read_csv(f'{indir}/{fname}_pes.csv'), channel, f'{fname}_{channel}.pdf')
fname = 'test_q100_1724098596'
- plot_final(pd.read_csv(f'{fname}.csv'), f'{fname}.pdf')
- plot_pes(pd.read_csv(f'{fname}_pes.csv'), channel, f'{fname}_{channel}.pdf')
+ plot_final(pd.read_csv(f'{indir}/{fname}.csv'), f'{fname}.pdf')
+ plot_pes(pd.read_csv(f'{indir}/{fname}_pes.csv'), channel, f'{fname}_{channel}.pdf')
- plot_chi2(pd.read_csv(f'quality.csv'), f'chi2_prepca.pdf')
- plot_chi2_intensity(pd.read_csv(f'quality.csv'), f'intensity_vs_chi2_prepca.pdf')
- plot_rmse_intensity(pd.read_csv(f'quality.csv'), f'intensity_vs_rmse.pdf')
+ plot_chi2(pd.read_csv(f'{indir}/quality.csv'), f'chi2_prepca.pdf')
+ plot_chi2_intensity(pd.read_csv(f'{indir}/quality.csv'), f'intensity_vs_chi2_prepca.pdf')
+ plot_rmse(pd.read_csv(f'{indir}/quality.csv'), f'rmse.pdf')
+ plot_rmse_intensity(pd.read_csv(f'{indir}/quality.csv'), f'intensity_vs_rmse.pdf')
- pca_variance_plot(pd.read_csv(f'pca_spec.csv'), f'pca_spec.pdf')
- pca_variance_plot(pd.read_csv(f'pca_pes.csv'), f'pca_pes.pdf')
+ pca_variance_plot(pd.read_csv(f'{indir}/pca_spec.csv'), f'pca_spec.pdf')
+ pca_variance_plot(pd.read_csv(f'{indir}/pca_pes.csv'), f'pca_pes.pdf')