diff --git a/bin/rabbit_fit.py b/bin/rabbit_fit.py index 44893f1..2db7ee6 100755 --- a/bin/rabbit_fit.py +++ b/bin/rabbit_fit.py @@ -182,6 +182,13 @@ def make_parser(): action="store_true", help="compute impacts of frozen (non-profiled) systematics", ) + parser.add_argument( + "--setConstraintMinimum", + default=[], + nargs=2, + action="append", + help="Set the constraint minima of specified parameter to specified value", + ) return parser.parse_args() diff --git a/bin/rabbit_plot_likelihood_scan.py b/bin/rabbit_plot_likelihood_scan.py index a6a4752..a052db1 100755 --- a/bin/rabbit_plot_likelihood_scan.py +++ b/bin/rabbit_plot_likelihood_scan.py @@ -1,9 +1,7 @@ #!/usr/bin/env python3 import argparse -import os -import matplotlib.pyplot as plt import mplhep as hep import numpy as np @@ -16,32 +14,6 @@ logger = None -def writeOutput(fig, outfile, extensions=[], postfix=None, args=None, meta_info=None): - name, _ = os.path.splitext(outfile) - - if postfix: - name += f"_{postfix}" - - for ext in extensions: - if ext[0] != ".": - ext = "." + ext - output = name + ext - print(f"Write output file {output}") - plt.savefig(output) - - output = name.rsplit("/", 1) - output[1] = os.path.splitext(output[1])[0] - if len(output) == 1: - output = (None, *output) - if args is None and meta_info is None: - return - output_tools.write_logfile( - *output, - args=args, - meta_info=meta_info, - ) - - def parseArgs(): parser = argparse.ArgumentParser() parser.add_argument( @@ -73,6 +45,11 @@ def parseArgs(): default="./test", help="Folder path for output", ) + parser.add_argument( + "--eoscp", + action="store_true", + help="Override use of xrdcp and use the mount instead", + ) parser.add_argument( "-p", "--postfix", type=str, help="Postfix for output file name" ) @@ -113,6 +90,13 @@ def parseArgs(): default=None, help="x axis limits", ) + parser.add_argument( + "--ylim", + type=float, + nargs=2, + default=None, + help="y axis limits", + ) parser.add_argument("--titlePos", type=int, default=2, help="title position") parser.add_argument( "--config", @@ -133,6 +117,7 @@ def plot_scan( subtitle=None, titlePos=0, xlim=None, + ylim=None, combine=None, ylabel=r"$-2\,\Delta \log L$", config={}, @@ -148,13 +133,15 @@ def plot_scan( if xlim is None: xlim = (min(x), max(x)) + if ylim is None: + ylim = (min(y), max(y)) fig, ax = plot_tools.figure( x, xlabel, ylabel, xlim=xlim, - ylim=(min(y), max(y)), # logy=args.logy + ylim=ylim, # logy=args.logy ) ax.axhline(y=1, color="gray", linestyle="--", alpha=0.5) @@ -227,6 +214,7 @@ def plot_scan( def main(): args = parseArgs() + outdir = output_tools.make_plot_dir(args.outpath, eoscp=args.eoscp) global logger logger = logging.setup_logger(__file__, args.verbose, args.noColorLogger) @@ -280,19 +268,20 @@ def main(): subtitle=args.subtitle, titlePos=args.titlePos, xlim=args.xlim, + ylim=args.ylim, config=config, combine=(vals, nlls) if args.combine is not None else None, no_hessian=args.noHessian, ) - os.makedirs(args.outpath, exist_ok=True) - outfile = os.path.join(args.outpath, f"nll_scan_{param}") - writeOutput( - fig, - outfile=outfile, - extensions=["png", "pdf"], - meta_info=meta, + + to_join = [f"nll_scan_{param}", args.postfix] + outfile = "_".join(filter(lambda x: x, to_join)) + plot_tools.save_pdf_and_png(outdir, outfile) + output_tools.write_index_and_log( + outdir, + outfile, + analysis_meta_info=meta, args=args, - postfix=args.postfix, ) diff --git a/rabbit/fitter.py b/rabbit/fitter.py index 538eb12..a61c9d0 100644 --- a/rabbit/fitter.py +++ b/rabbit/fitter.py @@ -45,9 +45,7 @@ def __init__(self, xv): def __call__(self, intermediate_result): loss = intermediate_result.fun - logger.debug( - f"Iteration {self.iiter}: loss value {loss}" - ) # ; status {intermediate_result.status}") + logger.debug(f"Iteration {self.iiter}: loss value {loss}") if np.isnan(loss): raise ValueError(f"Loss value is NaN at iteration {self.iiter}") @@ -1979,14 +1977,12 @@ def _compute_lbeta(self, beta, full_nll=False): return None def _compute_nll_components(self, profile=True, full_nll=False): - nexpfullcentral, _, beta = self._compute_yields_with_beta( + nexp, _, beta = self._compute_yields_with_beta( profile=profile, compute_norm=False, full=False, ) - nexp = nexpfullcentral - if self.chisqFit: ln = 0.5 * tf.reduce_sum((nexp - self.nobs) ** 2 / self.varnobs, axis=-1) elif self.covarianceFit: @@ -2111,6 +2107,68 @@ def loss_val_grad_hess_beta(self, profile=True): return val, grad, hess + def fit(self): + def scipy_loss(xval): + self.x.assign(xval) + val, grad = self.loss_val_grad() + return val.__array__(), grad.__array__() + + def scipy_hessp(xval, pval): + self.x.assign(xval) + p = tf.convert_to_tensor(pval) + val, grad, hessp = self.loss_val_grad_hessp(p) + return hessp.__array__() + + def scipy_hess(xval): + self.x.assign(xval) + val, grad, hess = self.loss_val_grad_hess() + if self.diagnostics: + cond_number = tfh.cond_number(hess) + logger.info(f" - Condition number: {cond_number}") + edmval = tfh.edmval(grad, hess) + logger.info(f" - edmval: {edmval}") + return hess.__array__() + + xval = self.x.numpy() + + callback = FitterCallback(self) + + if self.minimizer_method in [ + "trust-krylov", + "trust-ncg", + ]: + info_minimize = dict(hessp=scipy_hessp) + elif self.minimizer_method in [ + "trust-exact", + "dogleg", + ]: + info_minimize = dict(hess=scipy_hess) + else: + info_minimize = dict() + + try: + res = scipy.optimize.minimize( + scipy_loss, + xval, + method=self.minimizer_method, + jac=True, + tol=0.0, + callback=callback, + **info_minimize, + ) + except Exception as ex: + # minimizer could have called the loss or hessp functions with "random" values, so restore the + # state from the end of the last iteration before the exception + xval = callback.xval + logger.debug(ex) + else: + xval = res["x"] + logger.debug(res) + + self.x.assign(xval) + + return callback + def minimize(self): if self.is_linear: logger.info( @@ -2138,65 +2196,7 @@ def minimize(self): callback = None else: - - def scipy_loss(xval): - self.x.assign(xval) - val, grad = self.loss_val_grad() - return val.__array__(), grad.__array__() - - def scipy_hessp(xval, pval): - self.x.assign(xval) - p = tf.convert_to_tensor(pval) - val, grad, hessp = self.loss_val_grad_hessp(p) - return hessp.__array__() - - def scipy_hess(xval): - self.x.assign(xval) - val, grad, hess = self.loss_val_grad_hess() - if self.diagnostics: - cond_number = tfh.cond_number(hess) - logger.info(f" - Condition number: {cond_number}") - edmval = tfh.edmval(grad, hess) - logger.info(f" - edmval: {edmval}") - return hess.__array__() - - xval = self.x.numpy() - - callback = FitterCallback(xval) - - if self.minimizer_method in [ - "trust-krylov", - "trust-ncg", - ]: - info_minimize = dict(hessp=scipy_hessp) - elif self.minimizer_method in [ - "trust-exact", - "dogleg", - ]: - info_minimize = dict(hess=scipy_hess) - else: - info_minimize = dict() - - try: - res = scipy.optimize.minimize( - scipy_loss, - xval, - method=self.minimizer_method, - jac=True, - tol=0.0, - callback=callback, - **info_minimize, - ) - except Exception as ex: - # minimizer could have called the loss or hessp functions with "random" values, so restore the - # state from the end of the last iteration before the exception - xval = callback.xval - logger.debug(ex) - else: - xval = res["x"] - logger.debug(res) - - self.x.assign(xval) + callback = self.fit() return callback @@ -2204,6 +2204,9 @@ def nll_scan(self, param, scan_range, scan_points, use_prefit=False): # make a likelihood scan for a single parameter # assuming the likelihood is minimized + # freeze minimize which mean to not update it in the fit + self.freeze_params(param) + idx = np.where(self.parms.astype(str) == param)[0][0] # store current state of x temporarily @@ -2229,51 +2232,28 @@ def nll_scan(self, param, scan_range, scan_points, use_prefit=False): if i == 0: continue + logger.debug(f"Now at i={i} x={ixval}") self.x.assign(tf.tensor_scatter_nd_update(self.x, [[idx]], [ixval])) - def scipy_loss(xval): - self.x.assign(xval) - val, grad = self.loss_val_grad() - grad = grad.numpy() - grad[idx] = 0 # Zero out gradient for the frozen parameter - return val.numpy(), grad - - def scipy_hessp(xval, pval): - self.x.assign(xval) - pval[idx] = ( - 0 # Ensure the perturbation does not affect frozen parameter - ) - p = tf.convert_to_tensor(pval) - val, grad, hessp = self.loss_val_grad_hessp(p) - hessp = hessp.numpy() - # TODO: worth testing modifying the loss/grad/hess functions to imply 1 - # for the corresponding hessian element instead of 0, - # since this might allow the minimizer to converge more efficiently - hessp[idx] = ( - 0 # Zero out Hessian-vector product at the frozen index - ) - return hessp - - res = scipy.optimize.minimize( - scipy_loss, - self.x, - method="trust-krylov", - jac=True, - hessp=scipy_hessp, - ) - if res["success"]: - dnlls[nscans // 2 + sign * i] = ( - self.reduced_nll().numpy() - nll_best - ) - scan_vals[nscans // 2 + sign * i] = ixval + self.fit() + + dnlls[nscans // 2 + sign * i] = self.reduced_nll().numpy() - nll_best + + scan_vals[nscans // 2 + sign * i] = ixval # reset x to original state self.x.assign(xval) + # let the parameter be free again + self.defreeze_params(param) + return scan_vals, dnlls def nll_scan2D(self, param_tuple, scan_range, scan_points, use_prefit=False): + # freeze minimize which mean to not update it in the fit + self.freeze_params(param_tuple) + idx0 = np.where(self.parms.astype(str) == param_tuple[0])[0][0] idx1 = np.where(self.parms.astype(str) == param_tuple[1])[0][0] @@ -2320,7 +2300,9 @@ def nll_scan2D(self, param_tuple, scan_range, scan_points, use_prefit=False): ix = best_fit - i iy = best_fit + j - # print(f"i={i}, j={j}, r={r} drow={drow}, dcol={dcol} | ix={ix}, iy={iy}") + logger.debug( + f"Now at (ix,iy) = ({ix},{iy}) (x,y)= ({x_scans[ix]},{y_scans[iy]})" + ) self.x.assign( tf.tensor_scatter_nd_update( @@ -2328,41 +2310,15 @@ def nll_scan2D(self, param_tuple, scan_range, scan_points, use_prefit=False): ) ) - def scipy_loss(xval): - self.x.assign(xval) - val, grad = self.loss_val_grad() - grad = grad.numpy() - grad[idx0] = 0 - grad[idx1] = 0 - return val.numpy(), grad - - def scipy_hessp(xval, pval): - self.x.assign(xval) - pval[idx0] = 0 - pval[idx1] = 0 - p = tf.convert_to_tensor(pval) - val, grad, hessp = self.loss_val_grad_hessp(p) - hessp = hessp.numpy() - hessp[idx0] = 0 - hessp[idx1] = 0 - - if np.allclose(hessp, 0, atol=1e-8): - return np.zeros_like(hessp) + self.fit() - return hessp + dnlls[ix, iy] = self.reduced_nll().numpy() - nll_best - res = scipy.optimize.minimize( - scipy_loss, - self.x, - method="trust-krylov", - jac=True, - hessp=scipy_hessp, - ) + self.x.assign(xval) - if res["success"]: - dnlls[ix, iy] = self.reduced_nll().numpy() - nll_best + # let the parameter be free again + self.defreeze_params(param_tuple) - self.x.assign(xval) return x_scans, y_scans, dnlls def contour_scan(self, param, nll_min, q=1, signs=[-1, 1], fun=None):