Skip to content

Commit

Permalink
remove deprecated correct_analytical_cov_eigenspectrum_ratio_gp
Browse files Browse the repository at this point in the history
  • Loading branch information
zatkins2 committed Oct 21, 2024
1 parent 703b545 commit 2b38925
Showing 1 changed file with 0 additions and 86 deletions.
86 changes: 0 additions & 86 deletions pspipe_utils/covariance.py
Original file line number Diff line number Diff line change
Expand Up @@ -368,92 +368,6 @@ def skew(cov, dir=1):
return corrected_cov



def correct_analytical_cov_eigenspectrum_ratio_gp(lb, an_full_cov, mc_full_cov,
var_eigenspectrum_ratios_by_block=None,
idx_arrs_by_block=None, return_all=False):
"""Correct an analytical covariance matrix with a monte carlo covariance
matrix. Assumes the following rotated monte carlo matrix is diagonal:
mc_rot = (ana**-.5) @ mc @ (ana**-.5).T
and then fits the diagonal in that basis with a Gaussian process (GP) on
the observed values:
ana_corrected = (ana**.5) @ np.diag(GP(np.diag(mc_rot))) @ (ana**.5).T
The GP is applied to each block diagonal separately.
Parameters
----------
lb : (nbin,) np.ndarray
The locations in ell of the bins.
an_full_cov : (nblock*nbin, nblock*nbin) np.ndarray
Analytic covariance matrix to be corrected.
mc_full_cov : (nblock*nbin, nblock*nbin) np.ndarray
Noisy monte carlo matrix to use to correct the analytic matrix.
var_eigenspectrum_ratios_by_block : (nblock, nbin) np.ndarray, optional
The variance of the diagonal elements of mc_rot, by default None. These
would need to be precomputed in the rotated basis. If None, the
noise level in the diagonal elements of mc_rot is estimated from the
elements themselves by the Gaussian process.
idx_arrs_by_block : (nblock,) list, optional
A list of np.ndarrays, each of which may have between 0 and nbin elements,
that gives for each block the bin indices that should actually be used
in the Gaussian Process, by default None. If the list is None, or if any
element of the list is None, all elements for all blocks (or for the
specific block) are used in the Gaussian process. If an element of the list
is an empty array, then no Gaussian process is run on that block; its
observed values are used without any smoothing applied to them.
return_all : bool, optional
If True, in addition to returning ana_corrected, also return the
diagonal of the mc_rot matrix, the smoothed diagonal, and a list
of the Gaussian process regression objects for each block,
by default False.
Returns
-------
(nblock*nbin, nblock*nbin) np.ndarray, {(nblock*nbin,) np.ndarray,
(nblock*nbin,) np.ndarray, (nblock,) list of sklearn.GaussianProcessRegressor
objects}
ana_corrected as above, and if return_all, np.diag(mc_rot),
GP(np.diag(mc_rot)), and a list of the GPs for each block. If for any
block the idx_arrs_by_block array is empty, the returned GP is None
since no GP was actually used for that block.
"""
sqrt_an_full_cov = utils.eigpow(an_full_cov, 0.5)
inv_sqrt_an_full_cov = np.linalg.inv(sqrt_an_full_cov)
res = inv_sqrt_an_full_cov @ mc_full_cov @ inv_sqrt_an_full_cov.T # res should be close to the identity if an_full_cov is good
res_diag = np.diag(res)

n_spec = len(res_diag) // len(lb)

res_diags = np.split(res_diag, n_spec)

if var_eigenspectrum_ratios_by_block is None:
var_eigenspectrum_ratios_by_block = [None] * n_spec

if idx_arrs_by_block is None:
idx_arrs_by_block = [None] * n_spec

smoothed_res_diags = []
gprs = []
for i in range(n_spec):
smoothed_gp_diag, gpr = smooth_gp_diag(lb, res_diags[i], var_eigenspectrum_ratios_by_block[i],
idx_arr=idx_arrs_by_block[i], return_gpr=True)
smoothed_res_diags.append(smoothed_gp_diag)
gprs.append(gpr)

smoothed_res_diag = np.hstack(smoothed_res_diags)

corrected_cov = sqrt_an_full_cov @ np.diag(smoothed_res_diag) @ sqrt_an_full_cov.T

if return_all:
return corrected_cov, res_diag, smoothed_res_diag, gprs
else:
return corrected_cov


def correct_analytical_cov_block_diag_gp(lb, an_full_cov, mc_full_cov,
var_mc_cov_anaflat=None,
idx_arrs_by_block=None, return_all=False):
Expand Down

0 comments on commit 2b38925

Please sign in to comment.