Commit f886d286 authored by Goran Jelic-Cizmek's avatar Goran Jelic-Cizmek
Browse files

Fixed ordering in covariance

This is to be able to add the `covariance_matrix` method, which returns the
block-diagonal matrix that can be used for, say, Fisher forecasts.
parent 5607ad09
......@@ -20,6 +20,7 @@ from typing import Any, Callable, List, Tuple, Union
from dataclasses import dataclass
import os
import numpy as np
from scipy.linalg import block_diag
def _check_parameter(
......@@ -1512,6 +1513,44 @@ cdef class Coffe:
])
def covariance_matrix(self, recompute : bool = False):
"""
Convenience function that returns the covariance of multipoles as a
numpy matrix.
"""
result = [_.value for _ in self.compute_covariance_bulk(recompute=recompute)]
z_size = len(self.z_mean)
rowsize = round(np.sqrt(len(result) // z_size))
return block_diag(
*[np.reshape(
result[i * len(result) // z_size : (i + 1) * len(result) // z_size],
(rowsize, rowsize)
) for i in range(z_size)]
)
def covariance_matrix_inverse(self, recompute : bool = False):
"""
Convenience function that returns the inverse of the covariance of
multipoles as a numpy matrix.
"""
result = [_.value for _ in self.compute_covariance_bulk(recompute=recompute)]
z_size = len(self.z_mean)
rowsize = round(np.sqrt(len(result) // z_size))
return block_diag(
*[np.linalg.inv(np.reshape(
result[i * len(result) // z_size : (i + 1) * len(result) // z_size],
(rowsize, rowsize)
)) for i in range(z_size)]
)
def _background_init(self):
ccoffe.coffe_background_init(
&self._parameters,
......
......@@ -785,29 +785,29 @@ int coffe_covariance_init(
const double coeffbar_array[] = {c0bar, c2bar, c4bar, c6bar, c8bar};
for (size_t i = 0; i<par->multipole_values_len; ++i){
for (size_t j = 0; j<par->multipole_values_len; ++j){
/* the sums c_i wigner3j^2(l1, l2, i) */
double coeff_sum = 0;
for (size_t cnt = 0; cnt < COFFE_ARRAY_SIZE(coeff_array); ++cnt){
coeff_sum +=
coeff_array[cnt]
*pow(gsl_sf_coupling_3j(2*par->multipole_values[i], 2*par->multipole_values[j], 4*cnt, 0, 0, 0), 2);
}
for (size_t m = 0; m < par->sep_len; ++m){
for (size_t j = 0; j<par->multipole_values_len; ++j){
for (size_t n = 0; n < par->sep_len; ++n){
/* the sums c_i wigner3j^2(l1, l2, i) */
double coeff_sum = 0;
for (size_t cnt = 0; cnt < COFFE_ARRAY_SIZE(coeff_array); ++cnt){
coeff_sum +=
coeff_array[cnt]
*pow(gsl_sf_coupling_3j(2*par->multipole_values[i], 2*par->multipole_values[j], 4*cnt, 0, 0, 0), 2);
}
double coeffbar_sum = 0;
for (size_t cnt = 0; cnt < COFFE_ARRAY_SIZE(coeffbar_array); ++cnt){
coeffbar_sum +=
coeffbar_array[cnt]
*pow(gsl_sf_coupling_3j(2*par->multipole_values[i], 2*par->multipole_values[j], 4*cnt, 0, 0, 0), 2);
}
double coeffbar_sum = 0;
for (size_t cnt = 0; cnt < COFFE_ARRAY_SIZE(coeffbar_array); ++cnt){
coeffbar_sum +=
coeffbar_array[cnt]
*pow(gsl_sf_coupling_3j(2*par->multipole_values[i], 2*par->multipole_values[j], 4*cnt, 0, 0, 0), 2);
}
/* Kronecker delta (l1, l2) */
double deltal1l2 =
(par->multipole_values[i] == par->multipole_values[j]) ? 1 : 0;
/* Kronecker delta (l1, l2) */
double deltal1l2 =
(par->multipole_values[i] == par->multipole_values[j]) ? 1 : 0;
for (size_t m = 0; m < par->sep_len; ++m){
for (size_t n = 0; n < par->sep_len; ++n){
double deltaij = (m == n) ? 1 : 0;
/* flat-sky covariance */
const double result_mp_or_ramp =
......@@ -880,8 +880,8 @@ int coffe_covariance_init(
size_t counter = 0;
for (size_t zi = 0; zi < par->z_mean_len; ++zi){
for (size_t l1 = 0; l1 < par->multipole_values_len; ++l1){
for (size_t l2 = 0; l2 < par->multipole_values_len; ++l2){
for (size_t r1 = 0; r1 < par->sep_len; ++r1){
for (size_t l2 = 0; l2 < par->multipole_values_len; ++l2){
for (size_t r2 = 0; r2 < par->sep_len; ++r2){
cov_mp->array[counter].coords.z_mean = par->z_mean[zi];
cov_mp->array[counter].coords.l1 = par->multipole_values[l1];
......
......@@ -149,3 +149,13 @@ class TestCoffe:
assert np.allclose(df.loc[(df.l1 == mp1) & (df.l2 == mp2)].r1.values, x)
assert np.allclose(df.loc[(df.l1 == mp1) & (df.l2 == mp2)].r2.values, y)
assert np.allclose(df.loc[(df.l1 == mp1) & (df.l2 == mp2)].value.values, z, rtol=5e-4)
assert np.allclose(
cosmo.covariance_matrix(),
np.transpose(cosmo.covariance_matrix())
)
assert np.allclose(
cosmo.covariance_matrix_inverse(),
np.linalg.inv(cosmo.covariance_matrix())
)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment