diff --git a/augur/analyze.py b/augur/analyze.py index ba3e0aa..5c3ea7b 100644 --- a/augur/analyze.py +++ b/augur/analyze.py @@ -207,52 +207,58 @@ def get_fisher_bias(self): if self.Fij is None: self.get_fisher_matrix() - _calculate_biased_cls = True + if self.bi is not None: + return self.bi - # Try to read the biased data vector - if 'biased_dv' in self.config['fisher_bias']: - _sys_path = self.config['fisher_bias']['biased_dv'] - if (len(_sys_path) < 1) or (os.path.exists(_sys_path) is False): - _calculate_biased_cls = True - else: - import astropy.table - if ('.dat' in _sys_path) or ('.txt' in _sys_path): - _format = 'ascii' - elif ('.fits' in _sys_path): - _format = 'fits' - else: - _format = None - biased_cls = astropy.table.Table.read(_sys_path, format=_format) - if len(biased_cls) != len(self.lk.get_data_vector()): - raise ValueError('The length of the provided Cls should be equal \ - to the length of the data-vector') - _calculate_biased_cls = False - self.biased_cls = biased_cls['dv_sys'] - self.data_fid + else: + _calculate_biased_cls = True - # If there's no biased data vector, calculate it - if _calculate_biased_cls: - _x_here = [] - _labels_here = [] - if 'bias_params' in self.config['fisher_bias'].keys(): - _pars_here = self.pars_fid.copy() - _sys_here = self.req_params.copy() - for key, value in self.config['fisher_bias']['bias_params'].items(): - if key in _pars_here.keys(): - _pars_here[key] = value - _x_here.append(value) - _labels_here.append(key) - elif key in _sys_here.keys(): - _sys_here[key] = value - _x_here.append(value) - _labels_here.append(key) + # Try to read the biased data vector + if 'biased_dv' in self.config['fisher_bias']: + _sys_path = self.config['fisher_bias']['biased_dv'] + if (len(_sys_path) < 1) or (os.path.exists(_sys_path) is False): + _calculate_biased_cls = True + else: + import astropy.table + if ('.dat' in _sys_path) or ('.txt' in _sys_path): + _format = 'ascii' + elif ('.fits' in _sys_path): + _format = 'fits' else: - raise ValueError(f'The requested parameter `{key}` is not recognized. \ - Please make sure that it is part of your model.') - else: - raise ValueError('bias_params is required if no biased_dv file is passed') + _format = None + biased_cls = astropy.table.Table.read(_sys_path, format=_format) + if len(biased_cls) != len(self.lk.get_data_vector()): + raise ValueError('The length of the provided Cls should be equal \ + to the length of the data-vector') + _calculate_biased_cls = False + self.biased_cls = biased_cls['dv_sys'] - self.data_fid + + # If there's no biased data vector, calculate it + if _calculate_biased_cls: + _x_here = [] + _labels_here = [] + if 'bias_params' in self.config['fisher_bias'].keys(): + _pars_here = self.pars_fid.copy() + _sys_here = self.req_params.copy() + for key, value in self.config['fisher_bias']['bias_params'].items(): + if key in _pars_here.keys(): + _pars_here[key] = value + _x_here.append(value) + _labels_here.append(key) + elif key in _sys_here.keys(): + _sys_here[key] = value + _x_here.append(value) + _labels_here.append(key) + else: + raise ValueError(f'The requested parameter `{key}` is not recognized. \ + Please make sure that it is part of your model.') + else: + raise ValueError('bias_params is required if no biased_dv file is passed') - self.biased_cls = self.f(_x_here, _labels_here, _pars_here, _sys_here) - self.data_fid + self.biased_cls = self.f(_x_here, _labels_here, _pars_here, _sys_here) \ + - self.data_fid - Bj = np.einsum('l, lm, jm', self.biased_cls, self.lk.inv_cov, self.derivatives) - bi = np.einsum('ij, j', np.linalg.inv(self.Fij), Bj) - self.bi = bi + Bj = np.einsum('l, lm, jm', self.biased_cls, self.lk.inv_cov, self.derivatives) + bi = np.einsum('ij, j', np.linalg.inv(self.Fij), Bj) + self.bi = bi + return self.bi diff --git a/augur/utils/cov_utils.py b/augur/utils/cov_utils.py index dad48c0..7ff99d8 100644 --- a/augur/utils/cov_utils.py +++ b/augur/utils/cov_utils.py @@ -89,15 +89,17 @@ def get_gaus_cov(S, lk, cosmo, fsky, config): cov_all = np.zeros((len(S.data), len(S.data))) # Loop over statistic in the likelihood (assuming 3x2pt so far) for i, myst1 in enumerate(lk.statistics): + myst1 = myst1.statistic tr1 = myst1.source0.tracers[0].ccl_tracer # Pulling out the tracers tr2 = myst1.source1.tracers[0].ccl_tracer - ell12 = myst1._ell_or_theta + ell12 = myst1.ells # Loop over upper-triangle and fill lower-triangle by symmetry for j in range(i, len(lk.statistics)): myst2 = lk.statistics[j] + myst2 = myst2.statistic tr3 = myst2.source0.tracers[0].ccl_tracer tr4 = myst2.source1.tracers[0].ccl_tracer - ell34 = myst2._ell_or_theta + ell34 = myst2.ells # Assuming that everything has the same ell-edges and we are just changing the length # TODO update this for a more general case if len(ell34) < len(ell12): diff --git a/examples/config_test.yml b/examples/config_test.yml index fc9d670..d7e501a 100644 --- a/examples/config_test.yml +++ b/examples/config_test.yml @@ -70,18 +70,18 @@ statistics: fiducial_sacc_path : test_sacc.sacc cov_options: # Several options implemented -- Gaussian internal - # cov_type : 'gaus_internal' - # fsky : 0.3 + #cov_type : 'gaus_internal' + #fsky : 0.3 # ------ # Or you can also get it from a file - cov_type : 'SRD' - SRD_cov_path : './data/Y1_3x2_SRD_cov.npy' + #cov_type : 'SRD' + #SRD_cov_path : './data/Y1_3x2_SRD_cov.npy' # Or using TJPCov - #cov_type : 'tjpcov' - # IA : 0.0 - # fsky: 0.3 - # binning_info : - # ell_edges : np.geomspace(20, 15000, 21, endpoint=True).astype(np.int32) + cov_type : 'tjpcov' + IA : 0.0 + fsky: 0.3 + binning_info : + ell_edges : np.geomspace(20, 15000, 21, endpoint=True).astype(np.int32) fisher: var_pars: ['Omega_c', 'sigma8', 'n_s', 'w0', 'wa', 'Omega_b', 'h', 'lens0_bias', 'lens1_bias',