From 79187fe4a64010417d77df97cdccd7745b0b15e3 Mon Sep 17 00:00:00 2001 From: Javier Sanchez Date: Wed, 30 Oct 2024 14:39:41 -0700 Subject: [PATCH] adding ability to use step size proportional to sampling interval --- augur/analyze.py | 65 ++++++++++++++++++++++++++++++++---------------- 1 file changed, 44 insertions(+), 21 deletions(-) diff --git a/augur/analyze.py b/augur/analyze.py index e949d48..c137730 100644 --- a/augur/analyze.py +++ b/augur/analyze.py @@ -5,10 +5,12 @@ from augur.utils.config_io import parse_config from firecrown.parameters import ParamsMap from astropy.table import Table +import warnings class Analyze(object): - def __init__(self, config, likelihood=None, tools=None, req_params=None): + def __init__(self, config, likelihood=None, tools=None, req_params=None, + norm_step=True): """ Run numerical derivatives of a likelihood to obtain a Fisher matrix estimate. @@ -25,6 +27,9 @@ def __init__(self, config, likelihood=None, tools=None, req_params=None): req_params: dict Dictionary containing the required parameters for the likelihood. Required if a likelihood object is passed. + norm_step : bool + If `True` it internally normalizes the step size using the sampling bounds in order + to determine the samples. Returns: -------- @@ -73,7 +78,7 @@ def __init__(self, config, likelihood=None, tools=None, req_params=None): self.tools = tools self.req_params = req_params self.data_fid = self.lk.get_data_vector() - + self.norm_step = norm_step # Get the fiducial cosmological parameters self.pars_fid = tools.get_ccl_cosmology().__dict__['_params_init_kwargs'] @@ -87,22 +92,27 @@ def __init__(self, config, likelihood=None, tools=None, req_params=None): self.Fij = None self.bi = None self.biased_cls = None + self.par_bounds = [] + self.norm = None # Load the parameters to vary # We will allow 2 options -- one where we pass something # a la cosmosis with parameters and minimum, central, and max # we can also allow priors if set(['parameters', 'var_pars']).issubset(set(self.config.keys())): - raise Warning('Both `parameters` and `var_pars` found in Fisher. \ + warnings.warn('Both `parameters` and `var_pars` found in Fisher. \ Ignoring `parameters and using fiducial cosmology \ as pivot.`') + if 'parameters' in self.config.keys(): self.var_pars = list(self.config['parameters'].keys()) for var in self.var_pars: _val = self.config['parameters'][var] if isinstance(_val, list): self.x.append(_val[1]) + self.par_bounds.append([_val[0], _val[-1]]) else: self.x.append(_val) + # The other option is to pass just the parameter names and evaluate around # the fiducial values elif 'var_pars' in self.config.keys(): @@ -117,6 +127,14 @@ def __init__(self, config, likelihood=None, tools=None, req_params=None): in the list of parameters in the likelihood.') # Cast to numpy array (this will be done later anyway) self.x = np.array(self.x) + self.par_bounds = np.array(self.par_bounds) + if (len(self.par_bounds) < 1) & (self.norm_step): + self.norm_step = False + warnings.warn('Parameter bounds not provided -- the step will not be normalized') + # Normalize the pivot point given the sampling region + if self.norm_step: + self.norm = self.par_bounds[:, 1] - self.par_bounds[:, 0] + self.x = (self.x - self.par_bounds[:, 0]) * 1/self.norm def f(self, x, labels, pars_fid, sys_fid): """ @@ -127,14 +145,14 @@ def f(self, x, labels, pars_fid, sys_fid): ----------- x : float, list or np.ndarray - Point at which to compute the theory vector + Point at which to compute the theory vector. labels : list - Names of parameters to vary + Names of parameters to vary. pars_fid : dict Dictionary containing the fiducial ccl cosmological parameters sys_fid: dict Dictionary containing the fiducial `systematic` (required) parameters - for the likelihood + for the likelihood. Returns: -------- f_out : np.ndarray @@ -146,6 +164,11 @@ def f(self, x, labels, pars_fid, sys_fid): else: if isinstance(x, list): x = np.array(x) + # If we normalize the sampling we need to undo the normalization + if self.norm_step: + print(x) + x = self.norm * x + self.par_bounds[:, 0] + print(x) if x.ndim == 1: _pars = pars_fid.copy() _sys_pars = sys_fid.copy() @@ -156,6 +179,8 @@ def f(self, x, labels, pars_fid, sys_fid): _sys_pars.update({labels[i]: x[i]}) else: raise ValueError(f'Parameter name {labels[i]} not recognized!') + print(_pars) + print(_sys_pars) self.tools.reset() self.lk.reset() pmap = ParamsMap(_sys_pars) @@ -187,39 +212,37 @@ def f(self, x, labels, pars_fid, sys_fid): f_out.append(self.lk.compute_theory_vector(self.tools)) return np.array(f_out) - def get_derivatives(self, force=False, method='5pt_stencil', normalize_params=True): + def get_derivatives(self, force=False, method='5pt_stencil'): """ Auxiliary function to compute numerical derivatives of the helper function `f` Parameters: ----------- - force : bool, If `True` force recalculation of the derivatives - method : str, Method to compute derivatives, currently only `5pt_stencil` or - numdifftools are accepted. - normalize_params : bool, If `True` it normalizes the parameters before computing the - derivatives. + force : bool + If `True` force recalculation of the derivatives. + method : str + Method to compute derivatives, currently only `5pt_stencil` or `numdifftools` + are allowed. """ - if normalize_params: - x_piv = self.x # Modify - else: - x_piv = self.x + step = float(self.config['step']) # Compute the derivatives with respect to the parameters in var_pars at x if (self.derivatives is None) or (force): if '5pt_stencil' in method: self.derivatives = five_pt_stencil(lambda y: self.f(y, self.var_pars, self.pars_fid, self.req_params), - x_piv, h=float(self.config['step'])) + self.x, h=step) elif 'numdifftools' in method: import numdifftools as nd if 'numdifftools_kwargs' in self.config.keys(): ndkwargs = self.config['numdifftools_kwargs'] else: ndkwargs = {} - self.derivatives = nd.Gradient(lambda y: self.f(y, self.var_pars, self.pars_fid, - self.req_params), - step=float(self.config['step']), - **ndkwargs)(x_piv).T + jacobian_calc = nd.Jacobian(lambda y: self.f(y, self.var_pars, self.pars_fid, + self.req_params), + step=step, + **ndkwargs) + self.derivatives = jacobian_calc(self.x).T else: raise ValueError(f'Selected method: `{method}` is not available. \ Please select 5pt_stencil or numdifftools.')