Skip to content
Snippets Groups Projects
Commit 4a93744c authored by Danilo Ferreira de Lima's avatar Danilo Ferreira de Lima
Browse files

Adapted offline script for paper submission.

parent 9949e757
No related branches found
No related tags found
1 merge request!14Corrected bugs in the BNN and added many plotting scripts adapted for the paper
Hide whitespace changes
Inline Side-by-side
pes_to_spec/bnn.py
+ 1
− 1
View file @ 4a93744c
...@@ -222,7 +222,7 @@ class BNNModel(RegressorMixin, BaseEstimator): ...@@ -222,7 +222,7 @@ class BNNModel(RegressorMixin, BaseEstimator):
# train # train
self.model.train() self.model.train()
epochs = 500 epochs = 1000
for epoch in range(epochs): for epoch in range(epochs):
meter = {k: AverageMeter(k, ':6.3f') meter = {k: AverageMeter(k, ':6.3f')
for k in ('loss', '-log(lkl)', '-log(prior)', '-log(hyper)', 'sigma', 'w.prec.')} for k in ('loss', '-log(lkl)', '-log(prior)', '-log(hyper)', 'sigma', 'w.prec.')}
......
pes_to_spec/test/offline_analysis.py
+ 103
− 33
View file @ 4a93744c
...@@ -11,6 +11,10 @@ import numpy as np ...@@ -11,6 +11,10 @@ import numpy as np
from extra_data import open_run, by_id, RunDirectory from extra_data import open_run, by_id, RunDirectory
from pes_to_spec.model import Model, matching_ids 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 from itertools import product
import matplotlib import matplotlib
...@@ -71,11 +75,46 @@ def plot_pes(filename: str, pes_raw_int: np.ndarray, first: int, last: int): ...@@ -71,11 +75,46 @@ def plot_pes(filename: str, pes_raw_int: np.ndarray, first: int, last: int):
fig = plt.figure(figsize=(16, 8)) fig = plt.figure(figsize=(16, 8))
gs = GridSpec(1, 1) gs = GridSpec(1, 1)
ax = fig.add_subplot(gs[0, 0]) 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.plot(np.arange(first, last), pes_raw_int, c='b', lw=3)
ax.legend() #ax.legend()
ax.set(title=f"", ax.set(title=f"",
xlabel="ToF index", xlabel="Time-of-flight index",
ylabel="Intensity") 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) fig.savefig(filename)
plt.close(fig) plt.close(fig)
...@@ -113,11 +152,12 @@ def plot_result(filename: str, ...@@ -113,11 +152,12 @@ def plot_result(filename: str,
unc=unc, unc=unc,
beam_intensity=intensity*1e-3*np.ones_like(spec_raw_pe) beam_intensity=intensity*1e-3*np.ones_like(spec_raw_pe)
)) ))
df.to_csv(filename.replace('.png', '.csv')) df.to_csv(filename.replace('.pdf', '.csv'))
pes_data = deepcopy(pes) if pes is not None:
pes_data['bin'] = np.arange(len(pes['channel_1_D'])) pes_data = deepcopy(pes)
df = pd.DataFrame(pes_data) pes_data['bin'] = np.arange(len(pes['channel_1_D']))
df.to_csv(filename.replace('.png', '_pes.csv')) df = pd.DataFrame(pes_data)
df.to_csv(filename.replace('.pdf', '_pes.csv'))
fig = plt.figure(figsize=(12, 8)) fig = plt.figure(figsize=(12, 8))
gs = GridSpec(1, 1) gs = GridSpec(1, 1)
...@@ -140,9 +180,9 @@ def plot_result(filename: str, ...@@ -140,9 +180,9 @@ def plot_result(filename: str,
ax.spines['right'].set_visible(False) ax.spines['right'].set_visible(False)
ax.set( ax.set(
xlabel="Photon energy [eV]", xlabel="Photon energy [eV]",
ylabel="Intensity", ylabel="Intensity [a.u.]",
ylim=(0, 1.3*Y)) 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]) ax2 = plt.axes([0,0,1,1])
# Manually set the position and relative size of the inset axes within ax1 # Manually set the position and relative size of the inset axes within ax1
#ip = InsetPosition(ax, [0.65,0.6,0.35,0.4]) #ip = InsetPosition(ax, [0.65,0.6,0.35,0.4])
...@@ -180,7 +220,7 @@ def main(): ...@@ -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('-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('-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('-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('-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('-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') parser.add_argument('-o', '--offset', type=int, metavar='INT', default=0, help='Train ID offset')
...@@ -253,8 +293,10 @@ def main(): ...@@ -253,8 +293,10 @@ def main():
print(f"Number of test IDs: {len(test_tids)}") print(f"Number of test IDs: {len(test_tids)}")
# read the PES data for each channel # 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}" 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() pes_raw = {ch: run[pes_name, f"digitizers.{ch}.raw.samples"].select_trains(by_id[tids]).ndarray()
for ch in channels} for ch in channels}
pes_raw_t = {ch: run_test[pes_name, f"digitizers.{ch}.raw.samples"].select_trains(by_id[test_tids]).ndarray() 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(): ...@@ -296,14 +338,42 @@ def main():
t = list() t = list()
t_names = 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) 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 # 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(): for k in pes_raw.keys():
pes_raw[k] = pes_raw[k][train_idx] 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: if len(args.model) == 0:
print("Fitting") print("Fitting")
start = time_ns() start = time_ns()
...@@ -344,7 +414,7 @@ def main(): ...@@ -344,7 +414,7 @@ def main():
ylabel="Response [a.u.]", ylabel="Response [a.u.]",
yscale='log', yscale='log',
) )
fig.savefig(os.path.join(args.directory, "impulse.png")) fig.savefig(os.path.join(args.directory, "impulse.pdf"))
plt.close(fig) plt.close(fig)
print(f"Resolution: {model.resolution:.2f} eV") print(f"Resolution: {model.resolution:.2f} eV")
...@@ -358,7 +428,7 @@ def main(): ...@@ -358,7 +428,7 @@ def main():
ylabel="Filter intensity [a.u.]", ylabel="Filter intensity [a.u.]",
yscale='log', 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) plt.close(fig)
fig = plt.figure(figsize=(12, 8)) fig = plt.figure(figsize=(12, 8))
...@@ -370,7 +440,7 @@ def main(): ...@@ -370,7 +440,7 @@ def main():
ylabel="Filter value [a.u.]", ylabel="Filter value [a.u.]",
yscale='log', yscale='log',
) )
fig.savefig(os.path.join(args.directory, "wiener.png")) fig.savefig(os.path.join(args.directory, "wiener.pdf"))
plt.close(fig) plt.close(fig)
print("Check consistency") print("Check consistency")
...@@ -434,7 +504,7 @@ def main(): ...@@ -434,7 +504,7 @@ def main():
ax2.xaxis.label.set_size(SMALL_SIZE) ax2.xaxis.label.set_size(SMALL_SIZE)
ax2.yaxis.label.set_size(SMALL_SIZE) ax2.yaxis.label.set_size(SMALL_SIZE)
ax2.tick_params(axis='both', which='major', labelsize=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) plt.close(fig)
fig = plt.figure(figsize=(12, 8)) fig = plt.figure(figsize=(12, 8))
...@@ -454,7 +524,7 @@ def main(): ...@@ -454,7 +524,7 @@ def main():
# verticalalignment='top', horizontalalignment='right', # verticalalignment='top', horizontalalignment='right',
# transform=ax.transAxes, # transform=ax.transAxes,
# color='black', fontsize=15) # 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) plt.close(fig)
spec_smooth_pca = model.y_model['pca'].transform(spec_smooth) spec_smooth_pca = model.y_model['pca'].transform(spec_smooth)
...@@ -478,7 +548,7 @@ def main(): ...@@ -478,7 +548,7 @@ def main():
# verticalalignment='top', horizontalalignment='right', # verticalalignment='top', horizontalalignment='right',
# transform=ax.transAxes, # transform=ax.transAxes,
# color='black', fontsize=15) # 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) plt.close(fig)
fig = plt.figure(figsize=(12, 8)) fig = plt.figure(figsize=(12, 8))
...@@ -491,7 +561,7 @@ def main(): ...@@ -491,7 +561,7 @@ def main():
xlim=(0, 5), xlim=(0, 5),
ylim=(0, np.mean(xgm_flux_t) + 3*np.std(xgm_flux_t)) 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) plt.close(fig)
res_prepca = np.sum((spec_smooth_pca[:, np.newaxis, :] - spec_pred["expected_pca"])/spec_pred["expected_pca_unc"], axis=1) 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(): ...@@ -508,7 +578,7 @@ def main():
xlim=(-3, 3), xlim=(-3, 3),
) )
ax.legend(frameon=False) 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) plt.close(fig)
fig = plt.figure(figsize=(12, 8)) fig = plt.figure(figsize=(12, 8))
...@@ -528,7 +598,7 @@ def main(): ...@@ -528,7 +598,7 @@ def main():
# transform=ax.transAxes, # transform=ax.transAxes,
# color='black', fontsize=15) # color='black', fontsize=15)
plt.tight_layout() plt.tight_layout()
fig.savefig(os.path.join(args.directory, "intensity.png")) fig.savefig(os.path.join(args.directory, "intensity.pdf"))
plt.close(fig) plt.close(fig)
# rmse # rmse
...@@ -542,7 +612,7 @@ def main(): ...@@ -542,7 +612,7 @@ def main():
xlabel=r"Root-mean-squared error", xlabel=r"Root-mean-squared error",
ylabel="XGM intensity [uJ]", 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) plt.close(fig)
fig = plt.figure(figsize=(12, 8)) fig = plt.figure(figsize=(12, 8))
...@@ -561,7 +631,7 @@ def main(): ...@@ -561,7 +631,7 @@ def main():
# verticalalignment='top', horizontalalignment='right', # verticalalignment='top', horizontalalignment='right',
# transform=ax.transAxes, # transform=ax.transAxes,
# color='black', fontsize=15) # 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) plt.close(fig)
## SPEC integral w.r.t XGM intensity ## SPEC integral w.r.t XGM intensity
...@@ -600,16 +670,16 @@ def main(): ...@@ -600,16 +670,16 @@ def main():
#plt.close(fig) #plt.close(fig)
first, last = model.get_low_resolution_range() first, last = model.get_low_resolution_range()
first = max(0, first+250) first = max(0, first+200)
last = min(last, pes_raw_t["channel_1_D"].shape[1]-1) last = min(last-200, pes_raw_t["channel_1_D"].shape[1]-1)
pes_to_show = 'sum' pes_to_show = "" #'channel_1_D'
# plot # plot
high_int_idx = np.argsort(xgm_flux_t[:,0]) high_int_idx = np.argsort(xgm_flux_t[:,0])
for q in [10, 25, 50, 75, 100]: for q in [10, 25, 50, 75, 100]:
qi = int(len(high_int_idx)*(q/100.0)) qi = int(len(high_int_idx)*(q/100.0))
for idx in high_int_idx[qi-10:qi]: for idx in high_int_idx[qi-10:qi]:
tid = test_tids[idx] 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" {k: item[idx, 0, ...] if k != "pca"
else item[0, ...] else item[0, ...]
for k, item in spec_pred.items()}, for k, item in spec_pred.items()},
...@@ -622,9 +692,9 @@ def main(): ...@@ -622,9 +692,9 @@ def main():
pes_to_show=pes_to_show, pes_to_show=pes_to_show,
pes_bin=np.arange(first, last), pes_bin=np.arange(first, last),
) )
#for ch in channels: for ch in channels:
# plot_pes(os.path.join(args.directory, f"test_pes_{tid}_{ch}.png"), plot_pes(os.path.join(args.directory, f"test_pes_{tid}_{ch}.pdf"),
# pes_raw_t[ch][idx, first:last], first, last) -pes_raw_t[ch][idx, first:last], first, last)
if __name__ == '__main__': if __name__ == '__main__':
main() main()
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Legal Notice | Contacts