From ff862a89b3dc25323e4bf111e0321a2fdd9d6a35 Mon Sep 17 00:00:00 2001 From: Javier Sanchez Date: Wed, 30 Oct 2024 10:11:32 -0700 Subject: [PATCH 1/5] Addressed #67 and started implementation of #68 --- augur/analyze.py | 26 +++++++++++++++++++++----- 1 file changed, 21 insertions(+), 5 deletions(-) diff --git a/augur/analyze.py b/augur/analyze.py index baa8748..e949d48 100644 --- a/augur/analyze.py +++ b/augur/analyze.py @@ -32,7 +32,7 @@ def __init__(self, config, likelihood=None, tools=None, req_params=None): Output Fisher matrix """ - _config = parse_config(config) # Load full config + config = parse_config(config) # Load full config # Load the likelihood if no likelihood is passed along if likelihood is None: @@ -78,7 +78,7 @@ def __init__(self, config, likelihood=None, tools=None, req_params=None): self.pars_fid = tools.get_ccl_cosmology().__dict__['_params_init_kwargs'] # Load the relevant section of the configuration file - self.config = _config['fisher'] + self.config = config['fisher'] # Initialize pivot point self.x = [] @@ -187,13 +187,29 @@ 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'): + def get_derivatives(self, force=False, method='5pt_stencil', normalize_params=True): + """ + 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. + """ + if normalize_params: + x_piv = self.x # Modify + else: + x_piv = self.x + # 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), - self.x, h=float(self.config['step'])) + x_piv, h=float(self.config['step'])) elif 'numdifftools' in method: import numdifftools as nd if 'numdifftools_kwargs' in self.config.keys(): @@ -203,7 +219,7 @@ def get_derivatives(self, force=False, method='5pt_stencil'): self.derivatives = nd.Gradient(lambda y: self.f(y, self.var_pars, self.pars_fid, self.req_params), step=float(self.config['step']), - **ndkwargs)(self.x).T + **ndkwargs)(x_piv).T else: raise ValueError(f'Selected method: `{method}` is not available. \ Please select 5pt_stencil or numdifftools.') From 79187fe4a64010417d77df97cdccd7745b0b15e3 Mon Sep 17 00:00:00 2001 From: Javier Sanchez Date: Wed, 30 Oct 2024 14:39:41 -0700 Subject: [PATCH 2/5] 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.') From 157905760aa0baefb81e90670e99c449c93972d9 Mon Sep 17 00:00:00 2001 From: Javier Sanchez Date: Wed, 30 Oct 2024 15:31:02 -0700 Subject: [PATCH 3/5] added ability to change step size on the fly --- augur/analyze.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/augur/analyze.py b/augur/analyze.py index c137730..c6c5391 100644 --- a/augur/analyze.py +++ b/augur/analyze.py @@ -212,7 +212,7 @@ 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'): + def get_derivatives(self, force=False, method='5pt_stencil', step=None): """ Auxiliary function to compute numerical derivatives of the helper function `f` @@ -223,9 +223,12 @@ def get_derivatives(self, force=False, method='5pt_stencil'): method : str Method to compute derivatives, currently only `5pt_stencil` or `numdifftools` are allowed. + step : float + Step size for numerical differentiation """ - step = float(self.config['step']) + if step is None: + 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: From a6d9a8a0578e7b0670435d39424d2f4f9824e206 Mon Sep 17 00:00:00 2001 From: Javier Sanchez Date: Thu, 31 Oct 2024 10:52:23 -0700 Subject: [PATCH 4/5] clean up debugging print statements --- augur/analyze.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/augur/analyze.py b/augur/analyze.py index c6c5391..3338f4a 100644 --- a/augur/analyze.py +++ b/augur/analyze.py @@ -166,9 +166,7 @@ def f(self, x, labels, pars_fid, sys_fid): 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() @@ -179,8 +177,6 @@ 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) From d380288934432dd712f0d3ae55b63218f3152f1d Mon Sep 17 00:00:00 2001 From: Javier Sanchez Date: Thu, 31 Oct 2024 16:19:40 -0700 Subject: [PATCH 5/5] fixed bug when using normalized step --- augur/analyze.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/augur/analyze.py b/augur/analyze.py index 3338f4a..0bb254c 100644 --- a/augur/analyze.py +++ b/augur/analyze.py @@ -245,6 +245,8 @@ def get_derivatives(self, force=False, method='5pt_stencil', step=None): else: raise ValueError(f'Selected method: `{method}` is not available. \ Please select 5pt_stencil or numdifftools.') + if self.norm is not None: + self.derivatives /= self.norm[:, None] return self.derivatives else: return self.derivatives