diff --git a/pes_to_spec/bnn.py b/pes_to_spec/bnn.py
index 887512e093e555bb65bef8442585eb0c84801146..32461c28a8c9313fe8adea470240d397f12e5627 100644
--- a/pes_to_spec/bnn.py
+++ b/pes_to_spec/bnn.py
@@ -222,7 +222,7 @@ class BNNModel(RegressorMixin, BaseEstimator):
# train
self.model.train()
- epochs = 500
+ epochs = 1000
for epoch in range(epochs):
meter = {k: AverageMeter(k, ':6.3f')
for k in ('loss', '-log(lkl)', '-log(prior)', '-log(hyper)', 'sigma', 'w.prec.')}
diff --git a/pes_to_spec/test/offline_analysis.py b/pes_to_spec/test/offline_analysis.py
index a9ba79450dffcacbcc0b93b19eec85337c741174..871670e92e0c8b13091b1fa6cbf5128656a592d2 100755
--- a/pes_to_spec/test/offline_analysis.py
+++ b/pes_to_spec/test/offline_analysis.py
@@ -11,6 +11,10 @@ import numpy as np
from extra_data import open_run, by_id, RunDirectory
from pes_to_spec.model import Model, matching_ids
+# for helper plots
+from pes_to_spec.model import SelectRelevantLowResolution
+from sklearn.decomposition import PCA
+
from itertools import product
import matplotlib
@@ -71,11 +75,46 @@ def plot_pes(filename: str, pes_raw_int: np.ndarray, first: int, last: int):
fig = plt.figure(figsize=(16, 8))
gs = GridSpec(1, 1)
ax = fig.add_subplot(gs[0, 0])
- ax.plot(np.arange(first, last), pes_raw_int, c='b', lw=3, label="Low-resolution measurement")
- ax.legend()
+ ax.plot(np.arange(first, last), pes_raw_int, c='b', lw=3)
+ #ax.legend()
ax.set(title=f"",
- xlabel="ToF index",
- ylabel="Intensity")
+ xlabel="Time-of-flight index",
+ ylabel="Counts [a.u.]")
+ ax.spines['top'].set_visible(False)
+ ax.spines['right'].set_visible(False)
+ fig.savefig(filename)
+ plt.close(fig)
+
+def pca_variance_plot(filename: str, variance_ratio: np.ndarray, n_comp: int):
+ """
+ Plot variance contribution.
+
+ Args:
+ filename: Output file name.
+ variance_ratio: Contribution of each component's variance.
+
+ """
+ fig = plt.figure(figsize=(8, 8))
+ gs = GridSpec(1, 1)
+ ax = fig.add_subplot(gs[0, 0])
+ c = np.cumsum(variance_ratio)
+ ax.bar(1+np.arange(len(variance_ratio)), variance_ratio*100, color='tab:red', alpha=0.3, label="Per component")
+ ax.plot(1+np.arange(len(variance_ratio)), c*100, c='tab:blue', lw=5, label="Cumulative")
+ ax.plot([n_comp, n_comp], [0, c[n_comp]*100], lw=3, ls='--', c='m', label="Components kept")
+ ax.plot([0, n_comp], [c[n_comp]*100, c[n_comp]*100], lw=3, ls='--', c='m')
+ ax.legend(frameon=False)
+ print(f"PCA plot: total n. components: {len(variance_ratio)}")
+ x_max = np.where(c > 0.99)[0][0]
+ print(f"Fraction of variance: {c[n_comp]}")
+ ax.set_yscale('log')
+ ax.set(title=f"",
+ xlabel="Component",
+ ylabel="Variance [%]",
+ xlim=(1, x_max),
+ ylim=(0.1, 100))
+ ax.spines['top'].set_visible(False)
+ ax.spines['right'].set_visible(False)
+ plt.tight_layout()
fig.savefig(filename)
plt.close(fig)
@@ -113,11 +152,12 @@ def plot_result(filename: str,
unc=unc,
beam_intensity=intensity*1e-3*np.ones_like(spec_raw_pe)
))
- df.to_csv(filename.replace('.png', '.csv'))
- pes_data = deepcopy(pes)
- pes_data['bin'] = np.arange(len(pes['channel_1_D']))
- df = pd.DataFrame(pes_data)
- df.to_csv(filename.replace('.png', '_pes.csv'))
+ 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']))
+ df = pd.DataFrame(pes_data)
+ df.to_csv(filename.replace('.pdf', '_pes.csv'))
fig = plt.figure(figsize=(12, 8))
gs = GridSpec(1, 1)
@@ -140,9 +180,9 @@ def plot_result(filename: str,
ax.spines['right'].set_visible(False)
ax.set(
xlabel="Photon energy [eV]",
- ylabel="Intensity",
+ ylabel="Intensity [a.u.]",
ylim=(0, 1.3*Y))
- if pes is not None:
+ if (pes is not None) and (pes_to_show != ""):
ax2 = plt.axes([0,0,1,1])
# Manually set the position and relative size of the inset axes within ax1
#ip = InsetPosition(ax, [0.65,0.6,0.35,0.4])
@@ -180,7 +220,7 @@ def main():
parser.add_argument('-t', '--test-run', type=int, metavar='INT', help='Run to test', default=None)
parser.add_argument('-m', '--model', type=str, metavar='FILENAME', default="", help='Model to load. If given, do not train a model and just do inference with this one.')
parser.add_argument('-d', '--directory', type=str, metavar='DIRECTORY', default=".", help='Where to save the results.')
- parser.add_argument('-S', '--spec', type=str, metavar='NAME', default="SA3_XTD10_SPECT/MDL/SPECTROMETER_SQS_NAVITAR:output", help='SPEC name')
+ parser.add_argument('-S', '--spec', type=str, metavar='NAME', default="SA3_XTD10_SPECT/MDL/SPECTROMETER_SCS_NAVITAR:output", help='SPEC name')
parser.add_argument('-P', '--pes', type=str, metavar='NAME', default="SA3_XTD10_PES/ADC/1:network", help='PES name')
parser.add_argument('-X', '--xgm', type=str, metavar='NAME', default="SA3_XTD10_XGM/XGM/DOOCS:output", help='XGM name')
parser.add_argument('-o', '--offset', type=int, metavar='INT', default=0, help='Train ID offset')
@@ -253,8 +293,10 @@ def main():
print(f"Number of test IDs: {len(test_tids)}")
# read the PES data for each channel
+ #channels = [f"channel_{i}_{l}"
+ # for i, l in product(range(1, 5), ["A", "B", "C", "D"])]
channels = [f"channel_{i}_{l}"
- for i, l in product(range(1, 5), ["A", "B", "C", "D"])]
+ for i, l in product([1,3,4], ["A", "B", "C", "D"])]
pes_raw = {ch: run[pes_name, f"digitizers.{ch}.raw.samples"].select_trains(by_id[tids]).ndarray()
for ch in channels}
pes_raw_t = {ch: run_test[pes_name, f"digitizers.{ch}.raw.samples"].select_trains(by_id[test_tids]).ndarray()
@@ -296,14 +338,42 @@ def main():
t = list()
t_names = list()
- model = Model(model_type=args.model_type)
+ model = Model(channels=channels, model_type=args.model_type)
train_idx = np.isin(tids, train_tids) & (xgm_flux[:,0] > args.xgm_cut)
# we just need this for training and we need to avoid copying it, which blows up the memoray usage
for k in pes_raw.keys():
pes_raw[k] = pes_raw[k][train_idx]
- model.debug_peak_finding(pes_raw, os.path.join(args.directory, "test_peak_finding.png"))
+ # apply only PES selection for plotting
+ x_select = SelectRelevantLowResolution(channels, tof_start=None, delta_tof=300, poly=True)
+ x_select.fit(pes_raw)
+ pes_raw_select = x_select.transform(pes_raw, pulse_energy=xgm_flux[train_idx])
+ ch = channels[0]
+ idx = 0
+ first = x_select.tof_start - x_select.delta_tof
+ last = x_select.tof_start + x_select.delta_tof
+ plot_pes(os.path.join(args.directory, f"pes_example_{ch}.pdf"),
+ -pes_raw[ch][idx, first:last],
+ first=first,
+ last=last)
+
+ # apply PCA for plotting
+ B, P, _ = pes_raw_select.shape
+ pes_raw_select = pes_raw_select.reshape((B*P, -1))
+ pca = PCA(None, whiten=True)
+ pca.fit(pes_raw_select)
+ pca_variance_plot(os.path.join(args.directory, f"pca_pes.pdf"),
+ pca.explained_variance_ratio_,
+ 600)
+
+ 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)
+
+ model.debug_peak_finding(pes_raw, os.path.join(args.directory, "test_peak_finding.pdf"))
if len(args.model) == 0:
print("Fitting")
start = time_ns()
@@ -344,7 +414,7 @@ def main():
ylabel="Response [a.u.]",
yscale='log',
)
- fig.savefig(os.path.join(args.directory, "impulse.png"))
+ fig.savefig(os.path.join(args.directory, "impulse.pdf"))
plt.close(fig)
print(f"Resolution: {model.resolution:.2f} eV")
@@ -358,7 +428,7 @@ def main():
ylabel="Filter intensity [a.u.]",
yscale='log',
)
- fig.savefig(os.path.join(args.directory, "wiener_ft.png"))
+ fig.savefig(os.path.join(args.directory, "wiener_ft.pdf"))
plt.close(fig)
fig = plt.figure(figsize=(12, 8))
@@ -370,7 +440,7 @@ def main():
ylabel="Filter value [a.u.]",
yscale='log',
)
- fig.savefig(os.path.join(args.directory, "wiener.png"))
+ fig.savefig(os.path.join(args.directory, "wiener.pdf"))
plt.close(fig)
print("Check consistency")
@@ -434,7 +504,7 @@ def main():
ax2.xaxis.label.set_size(SMALL_SIZE)
ax2.yaxis.label.set_size(SMALL_SIZE)
ax2.tick_params(axis='both', which='major', labelsize=SMALL_SIZE)
- fig.savefig(os.path.join(args.directory, "intensity_vs_chi2.png"))
+ fig.savefig(os.path.join(args.directory, "intensity_vs_chi2.pdf"))
plt.close(fig)
fig = plt.figure(figsize=(12, 8))
@@ -454,7 +524,7 @@ def main():
# verticalalignment='top', horizontalalignment='right',
# transform=ax.transAxes,
# color='black', fontsize=15)
- fig.savefig(os.path.join(args.directory, "chi2.png"))
+ fig.savefig(os.path.join(args.directory, "chi2.pdf"))
plt.close(fig)
spec_smooth_pca = model.y_model['pca'].transform(spec_smooth)
@@ -478,7 +548,7 @@ def main():
# verticalalignment='top', horizontalalignment='right',
# transform=ax.transAxes,
# color='black', fontsize=15)
- fig.savefig(os.path.join(args.directory, "chi2_prepca.png"))
+ fig.savefig(os.path.join(args.directory, "chi2_prepca.pdf"))
plt.close(fig)
fig = plt.figure(figsize=(12, 8))
@@ -491,7 +561,7 @@ def main():
xlim=(0, 5),
ylim=(0, np.mean(xgm_flux_t) + 3*np.std(xgm_flux_t))
)
- fig.savefig(os.path.join(args.directory, "intensity_vs_chi2_prepca.png"))
+ fig.savefig(os.path.join(args.directory, "intensity_vs_chi2_prepca.pdf"))
plt.close(fig)
res_prepca = np.sum((spec_smooth_pca[:, np.newaxis, :] - spec_pred["expected_pca"])/spec_pred["expected_pca_unc"], axis=1)
@@ -508,7 +578,7 @@ def main():
xlim=(-3, 3),
)
ax.legend(frameon=False)
- fig.savefig(os.path.join(args.directory, "res_prepca.png"))
+ fig.savefig(os.path.join(args.directory, "res_prepca.pdf"))
plt.close(fig)
fig = plt.figure(figsize=(12, 8))
@@ -528,7 +598,7 @@ def main():
# transform=ax.transAxes,
# color='black', fontsize=15)
plt.tight_layout()
- fig.savefig(os.path.join(args.directory, "intensity.png"))
+ fig.savefig(os.path.join(args.directory, "intensity.pdf"))
plt.close(fig)
# rmse
@@ -542,7 +612,7 @@ def main():
xlabel=r"Root-mean-squared error",
ylabel="XGM intensity [uJ]",
)
- fig.savefig(os.path.join(args.directory, "intensity_vs_rmse.png"))
+ fig.savefig(os.path.join(args.directory, "intensity_vs_rmse.pdf"))
plt.close(fig)
fig = plt.figure(figsize=(12, 8))
@@ -561,7 +631,7 @@ def main():
# verticalalignment='top', horizontalalignment='right',
# transform=ax.transAxes,
# color='black', fontsize=15)
- fig.savefig(os.path.join(args.directory, "rmse.png"))
+ fig.savefig(os.path.join(args.directory, "rmse.pdf"))
plt.close(fig)
## SPEC integral w.r.t XGM intensity
@@ -600,16 +670,16 @@ def main():
#plt.close(fig)
first, last = model.get_low_resolution_range()
- first = max(0, first+250)
- last = min(last, pes_raw_t["channel_1_D"].shape[1]-1)
- pes_to_show = 'sum'
+ first = max(0, first+200)
+ last = min(last-200, pes_raw_t["channel_1_D"].shape[1]-1)
+ pes_to_show = "" #'channel_1_D'
# plot
high_int_idx = np.argsort(xgm_flux_t[:,0])
for q in [10, 25, 50, 75, 100]:
qi = int(len(high_int_idx)*(q/100.0))
for idx in high_int_idx[qi-10:qi]:
tid = test_tids[idx]
- plot_result(os.path.join(args.directory, f"test_q{q}_{tid}.png"),
+ plot_result(os.path.join(args.directory, f"test_q{q}_{tid}.pdf"),
{k: item[idx, 0, ...] if k != "pca"
else item[0, ...]
for k, item in spec_pred.items()},
@@ -622,9 +692,9 @@ def main():
pes_to_show=pes_to_show,
pes_bin=np.arange(first, last),
)
- #for ch in channels:
- # plot_pes(os.path.join(args.directory, f"test_pes_{tid}_{ch}.png"),
- # pes_raw_t[ch][idx, first:last], first, last)
+ for ch in channels:
+ plot_pes(os.path.join(args.directory, f"test_pes_{tid}_{ch}.pdf"),
+ -pes_raw_t[ch][idx, first:last], first, last)
if __name__ == '__main__':
main()