From f7e9cf77e5e026b496bec32003265902ec56f00a Mon Sep 17 00:00:00 2001 From: Mikael Mieskolainen Date: Sun, 3 Nov 2024 23:30:51 +0000 Subject: [PATCH] peakfit: add 2D profile scan visualization and logging --- configs/peakfit/tune0.yml | 20 +++++++++------- configs/peakfit/tune2.yml | 20 +++++++++------- icefit/icepeak.py | 6 ++--- icefit/peakfit.py | 49 ++++++++++++++++++++++++++++++++------- requirements-aux.txt | 2 +- 5 files changed, 66 insertions(+), 31 deletions(-) diff --git a/configs/peakfit/tune0.yml b/configs/peakfit/tune0.yml index fb98cb10..61a5b821 100644 --- a/configs/peakfit/tune0.yml +++ b/configs/peakfit/tune0.yml @@ -118,15 +118,17 @@ techno: minos: True # If True, Minuit MINOS uncertainties (slower but best), otherwise HESSE migrad_first_step: 0.001 # First step size (relative percentage), set None for migrad heuristic - min_ndof: 1 # Minimum number of d.o.f required - max_chi2: 1200 # Maximum chi2/ndf for a succesful fit - min_count: 5 # Minimum number of total entries (sum of bin content) in the histogram - trials: 10 # Number of random restarts (all are executed and inspected, the best is kept) - rand_sigma: 0.5 # Random perturbation scale for initial values per trial + min_ndof: 1 # Minimum number of d.o.f required + max_chi2: 1200 # Maximum chi2/ndf for a succesful fit + min_count: 5 # Minimum number of total entries (sum of bin content) in the histogram + trials: 10 # Number of random restarts (all are executed and inspected, the best is kept) + rand_sigma: 0.5 # Random perturbation scale for initial values per trial - set_to_nan: True # Set parameters after the fit to NaN if not passing max_chi2 or min_count + set_to_nan: True # Set parameters after the fit to NaN if not passing max_chi2 or min_count - strategy: 2 # Default 1 (0 ~ fast, 1 ~ balanced, 2 ~ accurate) - tol: 0.1 # Default 0.1, https://iminuit.readthedocs.io/en/stable/reference.html#iminuit.Minuit.tol + strategy: 2 # Default 1 (0 ~ fast, 1 ~ balanced, 2 ~ accurate) + tol: 0.1 # Default 0.1, https://iminuit.readthedocs.io/en/stable/reference.html#iminuit.Minuit.tol - cov_eps: 0.0 # Covariance matrix post-regularization (added to diagonal) (set to 0 for none, can bias!) + cov_eps: 0.0 # Covariance matrix post-regularization (added to diagonal) (set to 0 for none, can bias!) + + draw_mnmatrix: False # Draw MINOS profile likelihood scan 2D contours (very slow) diff --git a/configs/peakfit/tune2.yml b/configs/peakfit/tune2.yml index 4715ffa2..4e6c8e34 100644 --- a/configs/peakfit/tune2.yml +++ b/configs/peakfit/tune2.yml @@ -109,15 +109,17 @@ techno: minos: True # If True, Minuit MINOS uncertainties (slower but best), otherwise HESSE migrad_first_step: 0.001 # First step size (relative percentage), set None for migrad heuristic - min_ndof: 1 # Minimum number of d.o.f required - max_chi2: 1200 # Maximum chi2/ndf for a succesful fit - min_count: 5 # Minimum number of total entries (sum of bin content) in the histogram - trials: 10 # Number of random restarts (all are executed and inspected, the best is kept) - rand_sigma: 0.5 # Random perturbation scale for initial values per trial + min_ndof: 1 # Minimum number of d.o.f required + max_chi2: 1200 # Maximum chi2/ndf for a succesful fit + min_count: 5 # Minimum number of total entries (sum of bin content) in the histogram + trials: 10 # Number of random restarts (all are executed and inspected, the best is kept) + rand_sigma: 0.5 # Random perturbation scale for initial values per trial - set_to_nan: True # Set parameters after the fit to NaN if not passing max_chi2 or min_count + set_to_nan: True # Set parameters after the fit to NaN if not passing max_chi2 or min_count - strategy: 2 # Default 1 (0 ~ fast, 1 ~ balanced, 2 ~ accurate) - tol: 0.1 # Default 0.1, https://iminuit.readthedocs.io/en/stable/reference.html#iminuit.Minuit.tol + strategy: 2 # Default 1 (0 ~ fast, 1 ~ balanced, 2 ~ accurate) + tol: 0.1 # Default 0.1, https://iminuit.readthedocs.io/en/stable/reference.html#iminuit.Minuit.tol - cov_eps: 0.0 # Covariance matrix post-regularization (added to diagonal) (set to 0 for none, can bias!) + cov_eps: 0.0 # Covariance matrix post-regularization (added to diagonal) (set to 0 for none, can bias!) + + draw_mnmatrix: False # Draw MINOS profile likelihood scan 2D contours (very slow) diff --git a/icefit/icepeak.py b/icefit/icepeak.py index 3608c9c8..af0df837 100644 --- a/icefit/icepeak.py +++ b/icefit/icepeak.py @@ -704,7 +704,7 @@ def poiss_nll_loss(par): cprint(f"[{param['fit_type']}] chi2 / ndf = {chi2:.2f} / {ndof} = {chi2/ndof:.2f} [{str}]", 'yellow') print('') - return par, cov, var2pos + return par, cov, var2pos, best_m1 def generate_start_values(trial, param, techno): """ @@ -898,7 +898,7 @@ def analyze_1D_fit(hist, param: dict, techno: dict, fitfunc, Returns: output dictionary """ - print(__name__ + f'.analyze_1D_fit:') + cprint(__name__ + f'.analyze_1D_fit:', 'magenta') h = TH1_to_numpy(hist) d = hist_decompose(h, param=param, techno=techno) @@ -1071,7 +1071,7 @@ def analyze_1D_fit(hist, param: dict, techno: dict, fitfunc, # chi2 / ndf chi2_ndf = chi2 / ndof if ndof > 0 else -999 title = f"$\\chi^2 / n_\\mathrm{{dof}} = {chi2:.2f} / {ndof:0.0f} = {chi2_ndf:.2f}$" - print(f'plot title: {title}') + print(f"plot: {title.replace('$', '')} \n") plt.title(title) # --------------------------------------------------------------- diff --git a/icefit/peakfit.py b/icefit/peakfit.py index c069ad50..a1b54c81 100644 --- a/icefit/peakfit.py +++ b/icefit/peakfit.py @@ -205,11 +205,11 @@ def fit_task(f, inputparam, savepath, YEAR, GENTYPE): fitfunc_ = {f'{key}': fitfunc} # Fit and analyze - par,cov,var2pos = icepeak.binned_1D_fit(hist=hist_, param=param, fitfunc=fitfunc_, + par,cov,var2pos,m1 = icepeak.binned_1D_fit(hist=hist_, param=param, fitfunc=fitfunc_, techno=techno, par_fixed=None) output = icepeak.analyze_1D_fit(hist=hist, param=param, fitfunc=fitfunc, techno=techno, par=par, cov=cov, par_fixed=None) - + # Create savepath total_savepath = f'{savepath}/Run{YEAR}/{GENTYPE}/{SYST}' if not os.path.exists(total_savepath): @@ -222,8 +222,23 @@ def fit_task(f, inputparam, savepath, YEAR, GENTYPE): # Save the pull histogram plt.figure(output['fig_pulls']) - plt.savefig(f'{total_savepath}/{tree}_pulls.pdf', bbox_inches='tight') - plt.savefig(f'{total_savepath}/{tree}_pulls.png', bbox_inches='tight', dpi=300) + plt.savefig(f'{total_savepath}/{tree}__pulls.pdf', bbox_inches='tight') + plt.savefig(f'{total_savepath}/{tree}__pulls.png', bbox_inches='tight', dpi=300) + + # Save parameter outputs + with open(f'{total_savepath}/{tree}.log', "w") as file: + print(par, file=file) + print(m1.params, file=file) + print(cov, file=file) + + # Draw MINOS contours + if techno['draw_mnmatrix']: + cprint('Draw MINOS profile likelihood scan 2D contours', 'yellow') + + m1.draw_mnmatrix(figsize=(2.5*len(par), 2.5*len(par))) + plt.figure(plt.gcf()) + plt.savefig(f'{total_savepath}/{tree}__mnmatrix.pdf', bbox_inches='tight') + plt.savefig(f'{total_savepath}/{tree}__mnmatrix.png', bbox_inches='tight', dpi=300) # Save the fit numerical data par_dict, cov_arr = icepeak.iminuit2python(par=par, cov=cov, var2pos=var2pos) @@ -280,11 +295,11 @@ def fit_task_multi(f, inputparam, savepath, YEAR, GENTYPE): # -------------------------------------------------------------------- # Simultaneous fit of both - par,cov,var2pos = icepeak.binned_1D_fit(hist=hist, param=param, fitfunc=fitfunc, + par,cov,var2pos,m1 = icepeak.binned_1D_fit(hist=hist, param=param, fitfunc=fitfunc, techno=techno, par_fixed=par_fixed) # Analyze and plot Pass, Fail seperately - for key in f.keys(): + for idx, key in enumerate(f.keys()): output = icepeak.analyze_1D_fit(hist=hist[key], param=param, fitfunc=fitfunc[key], techno=techno, par=par, cov=cov, par_fixed=par_fixed) @@ -297,15 +312,31 @@ def fit_task_multi(f, inputparam, savepath, YEAR, GENTYPE): if not os.path.exists(total_savepath): os.makedirs(total_savepath) - # Save the fit plot + # Save the fit plot plt.figure(output['fig']) plt.savefig(f'{total_savepath}/{tree}.pdf', bbox_inches='tight') plt.savefig(f'{total_savepath}/{tree}.png', bbox_inches='tight', dpi=300) # Save the pull histogram plt.figure(output['fig_pulls']) - plt.savefig(f'{total_savepath}/{tree}_pulls.pdf', bbox_inches='tight') - plt.savefig(f'{total_savepath}/{tree}_pulls.png', bbox_inches='tight', dpi=300) + plt.savefig(f'{total_savepath}/{tree}__pulls.pdf', bbox_inches='tight') + plt.savefig(f'{total_savepath}/{tree}__pulls.png', bbox_inches='tight', dpi=300) + + # Save parameter outputs + if (idx == 0): + with open(f'{total_savepath}/{tree}.log', "w") as file: + print(par, file=file) + print(m1.params, file=file) + print(cov, file=file) + + # Draw MINOS contours + if techno['draw_mnmatrix'] and (idx == 0): + cprint('Draw MINOS profile likelihood scan 2D contours', 'yellow') + + m1.draw_mnmatrix(figsize=(2.5*len(par), 2.5*len(par))) + plt.figure(plt.gcf()) + plt.savefig(f'{total_savepath}/{tree}__mnmatrix.pdf', bbox_inches='tight') + plt.savefig(f'{total_savepath}/{tree}__mnmatrix.png', bbox_inches='tight', dpi=300) # Save the fit numerical data par_dict, cov_arr = icepeak.iminuit2python(par=par, cov=cov, var2pos=var2pos) diff --git a/requirements-aux.txt b/requirements-aux.txt index 6c9714a6..4bde751f 100644 --- a/requirements-aux.txt +++ b/requirements-aux.txt @@ -13,7 +13,7 @@ jax==0.4.28 jaxlib==0.4.28 # IMINUIT -iminuit==2.25.2 +iminuit==2.30.1 # MYSTIC mystic==0.4.2