-
Notifications
You must be signed in to change notification settings - Fork 0
/
gsf_pull_surface_slider.py
93 lines (73 loc) · 3.92 KB
/
gsf_pull_surface_slider.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
import logging
import pickle
from itertools import cycle
from pathlib import Path
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.widgets import Button
from apyts.constants import *
from apyts.gsf_utils import *
from test_utils import *
class ResidualSlideShow(SurfaceSlideShow):
def __init__(self):
self.n_surfaces = 0
self.n_rows = 0
self.residuals = []
self.row_names = []
self.n_cols = 4
self.col_names = ["LOC", "PHI", "QOP", "P"]
self.col_ranges = [(-0.1,0.1), (-0.05,0.05), (-0.5,0.5), (-1,1)]
def plot(self, ax, surface, row, col):
rmin = min(self.col_ranges[col][0], min(self.residuals[surface][row][:,col]))
rmax = max(self.col_ranges[col][1], max(self.residuals[surface][row][:,col]))
ax.hist(self.residuals[surface][row][:,col], bins="rice", range=(rmin, rmax), histtype="step")
ax.set_title("{} ({})".format(self.col_names[col], self.row_names[row]))
def add_row(self, sim_data, fit_data_surface_fn, rowname=""):
self.n_rows += 1
self.row_names.append(rowname)
n_surfaces = max([ len(truth) for _, truth, _ in sim_data ])
for surface in range(n_surfaces):
surface_filtered = fit_data_surface_fn(surface)
surface_truth = np.stack([ truth[surface] for _, truth, _ in sim_data if len(truth) > surface ])
surface_filtered = np.append(surface_filtered, abs(1./surface_filtered[:,eBoundQoP]).reshape(-1,1), axis=1)
surface_truth = np.append(surface_truth, abs(1./surface_truth[:,eBoundQoP]).reshape(-1,1), axis=1)
try:
self.residuals[surface].append(surface_filtered - surface_truth)
except:
self.residuals.append([surface_filtered - surface_truth])
self.n_surfaces = max(n_surfaces, self.n_surfaces)
if __name__ == "__main__":
input_path = Path("output/20230117_132335-500particles")
with open(input_path / "sim_result.pickle", "rb") as f:
data = pickle.load(f)
true_pars = data["true_pars"]
sim_results = data["sim_results"]
with open(input_path / "start_parameters.pickle", "rb") as f:
data = pickle.load(f)
start_pars = data["parameters"]
start_covs = data["covs"]
with open(input_path / "gsf_fit_result.pickle", "rb") as f:
gsf_fit_results = pickle.load(f)["fit_result"]
with open(input_path / "kf_fit_result.pickle", "rb") as f:
kf_fit_results = pickle.load(f)["fit_result"]
kf_mask = np.array([ False if data is None else True for data in kf_fit_results ])
kf_fit_results = [ d for d in kf_fit_results if d is not None ]
kf_sim_results = [ d for d, keep in zip(sim_results, kf_mask) if keep ]
kf_final_pars = np.stack([ data[0][0] for data in kf_fit_results ])
gsf_mask = np.array([ False if data is None else True for data in gsf_fit_results ])
gsf_fit_results = [ d for d in gsf_fit_results if d is not None ]
gsf_sim_results = [ d for d, keep in zip(sim_results, gsf_mask) if keep ]
gsf_final_pars = np.stack([ max(cmps, key=lambda c: c[0])[1] for cmps, _, _, _ in gsf_fit_results ])
kf_surface_filtered_fn = lambda i: \
np.stack([ filtered[i][0] for _, _, filtered, _ in kf_fit_results if len(filtered) > i ])
gsf_surface_filtered_maxw_fn = lambda i: \
np.stack([ max(filtered[i], key=lambda c: c[0])[1] for _, _, filtered, _ in gsf_fit_results if len(filtered) > i ])
gsf_surface_filtered_mean_fn = lambda i: \
np.stack([ gaussian_mixture_moments(filtered[i])[0] for _, _, filtered, _ in gsf_fit_results if len(filtered) > i ])
rs = ResidualSlideShow()
rs.add_row(kf_sim_results, kf_surface_filtered_fn, "KF")
rs.add_row(gsf_sim_results, gsf_surface_filtered_maxw_fn, "GSF|maxw")
rs.add_row(gsf_sim_results, gsf_surface_filtered_mean_fn, "GSF|mean")
rs.make_slideshow()
plt.show()