Skip to content

Commit

Permalink
bug fixes using internal covariance, adding method consistency
Browse files Browse the repository at this point in the history
  • Loading branch information
Javier Sanchez authored and Javier Sanchez committed Jun 13, 2024
1 parent 2c25bcd commit 06bea74
Show file tree
Hide file tree
Showing 3 changed files with 63 additions and 55 deletions.
94 changes: 50 additions & 44 deletions augur/analyze.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
6 changes: 4 additions & 2 deletions augur/utils/cov_utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down
18 changes: 9 additions & 9 deletions examples/config_test.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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',
Expand Down

0 comments on commit 06bea74

Please sign in to comment.