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

Added flags for the various terms in the covariance

parent dd64bda7
......@@ -214,6 +214,12 @@ cdef extern from "common.h":
int covariance_pop4
int covariance_cosmic
int covariance_mixed
int covariance_poisson
double *pixelsize
size_t pixelsize_len
......
......@@ -996,6 +996,47 @@ cdef class Coffe:
self._free_covariance_multipoles()
@property
def covariance_cosmic(self):
"""
Whether the CV-CV term should be taken into account when computing the
covariance of multipoles
"""
return bool(self._parameters.covariance_cosmic)
@covariance_cosmic.setter
def covariance_cosmic(self, value):
self._parameters.covariance_cosmic = int(bool(value))
self._free_covariance_multipoles()
@property
def covariance_mixed(self):
"""
Whether the mixed (CV-Poisson + Poisson-CV) term should be taken into account when computing the
covariance of multipoles
"""
return bool(self._parameters.covariance_mixed)
@covariance_mixed.setter
def covariance_mixed(self, value):
self._parameters.covariance_mixed = int(bool(value))
self._free_covariance_multipoles()
@property
def covariance_poisson(self):
"""
Whether the Poisson-Poisson term should be taken into account when computing the
covariance of multipoles
"""
return bool(self._parameters.covariance_poisson)
@covariance_poisson.setter
def covariance_poisson(self, value):
self._parameters.covariance_poisson = int(bool(value))
self._free_covariance_multipoles()
@property
def pixelsize(self):
"""
......
......@@ -443,6 +443,9 @@ typedef struct coffe_parameters_t
specify 4 different populations */
int covariance_pop1, covariance_pop2, covariance_pop3, covariance_pop4;
/* which contributions should be considered for the covariance (CV-CV, CV-Poisson, and Poisson-Poisson */
int covariance_cosmic, covariance_mixed, covariance_poisson;
int covariance_integration_method;
int covariance_integration_bins;
......
......@@ -622,6 +622,7 @@ int coffe_covariance_init(
#pragma omp parallel for num_threads(par->nthreads) collapse(2)
for (size_t m = 0; m < par->sep_len; ++m){
for (size_t n = 0; n < par->sep_len; ++n){
if (par->covariance_mixed){
integral_pk[index_redshift][i * par->multipole_values_len + j][par->sep_len * n + m] =
(2 * par->multipole_values[i] + 1)
*(2 * par->multipole_values[j] + 1)
......@@ -639,10 +640,11 @@ int coffe_covariance_init(
k_max
)
/M_PI;
}}
}}}
#pragma omp parallel for num_threads(par->nthreads) collapse(2)
for (size_t m = 0; m < par->sep_len; ++m){
for (size_t n = 0; n < par->sep_len; ++n){
if (par->covariance_cosmic){
integral_pk2[index_redshift][i * par->multipole_values_len + j][par->sep_len * n + m] =
(2 * par->multipole_values[i] + 1)
*(2 * par->multipole_values[j] + 1)
......@@ -661,7 +663,7 @@ int coffe_covariance_init(
)
/2.
/M_PI;
}}
}}}
}}
}
else if (par->covariance_integration_method == 2){
......@@ -921,10 +923,10 @@ int coffe_covariance_init(
coeff_pop2_pop3[1] = alpha2_22;
}
double coeff_cosmic_cosmic_sum = 0;
double coeff_cosmic_sum = 0;
for (size_t cnt1 = 0; cnt1 < COFFE_ARRAY_SIZE(coeff_pop1_pop3); ++cnt1){
for (size_t cnt2 = 0; cnt2 < COFFE_ARRAY_SIZE(coeff_pop1_pop3); ++cnt2){
coeff_cosmic_cosmic_sum += (
coeff_cosmic_sum += (
coeff_pop1_pop3[cnt1]
*coeff_pop2_pop4[cnt2]
+
......@@ -934,11 +936,11 @@ int coffe_covariance_init(
)
*coffe_legendre_integral(l1, l2, 2 * cnt1, 2 * cnt2);
}}
coeff_cosmic_cosmic_sum /= 4.;
coeff_cosmic_sum /= 4.;
double coeff_cosmic_poisson_sum = 0;
double coeff_mixed_sum = 0;
for (size_t cnt = 0; cnt < COFFE_ARRAY_SIZE(coeff_pop1_pop3); ++cnt){
coeff_cosmic_poisson_sum +=
coeff_mixed_sum +=
(
coeff_pop1_pop3[cnt]
*coffe_kronecker_delta(
......@@ -970,7 +972,7 @@ int coffe_covariance_init(
)
*coffe_legendre_integral(l1, l2, 2 * cnt, 0);
}
coeff_cosmic_poisson_sum /= 8.;
coeff_mixed_sum /= 8.;
/* Kronecker delta (l1, l2) */
const int deltal1l2 = coffe_kronecker_delta(l1, l2);
......@@ -978,7 +980,8 @@ int coffe_covariance_init(
/* Kronecker delta of the separations */
const int deltaij = coffe_kronecker_delta(m, n);
const double cov_poisson_poisson =
const double cov_poisson = (
par->covariance_poisson ?
(2 * l1 + 1)
/ 4. / M_PI
*deltal1l2
......@@ -1020,9 +1023,10 @@ int coffe_covariance_init(
)
/par->density1[k]
/par->density2[k]
/volume[k];
/volume[k] : 0);
const double cov_cosmic_poisson =
const double cov_mixed = (
par->covariance_mixed ?
coffe_sign((l2 - l1) / 2)
*coffe_sign(l1 + l2)
/* trigraph */
......@@ -1035,10 +1039,11 @@ int coffe_covariance_init(
:
integral_pk[0][i*par->multipole_values_len + j][par->sep_len*n + m]
)
*coeff_cosmic_poisson_sum
/volume[k];
*coeff_mixed_sum
/volume[k] : 0);
const double cov_cosmic_cosmic =
const double cov_cosmic = (
par->covariance_cosmic ?
coffe_sign((l2 - l1) / 2) *
/* trigraph */
(par->pk_type ? 1 : D1z * D1z * D1z * D1z)
......@@ -1051,16 +1056,16 @@ int coffe_covariance_init(
:
integral_pk2[0][i*par->multipole_values_len + j][par->sep_len*n + m]
)
*coeff_cosmic_cosmic_sum
/volume[k];
*coeff_cosmic_sum
/volume[k] : 0);
/* flat-sky covariance */
const double result_mp_or_ramp =
cov_poisson_poisson
cov_poisson
+
cov_cosmic_poisson
cov_mixed
+
cov_cosmic_cosmic;
cov_cosmic;
if (par->output_type == COVARIANCE_MULTIPOLES){
cov_mp->array[index].value = result_mp_or_ramp;
......
......@@ -818,6 +818,9 @@ int coffe_parse_default_parameters(
par->covariance_pop2 = 1;
par->covariance_pop3 = 1;
par->covariance_pop4 = 1;
par->covariance_cosmic = 1;
par->covariance_mixed = 1;
par->covariance_poisson = 1;
par->have_window = 0;
par->window_size = 0;
......
......@@ -258,6 +258,59 @@ class TestCoffe:
cosmo.compute_multipole(l=4, r=10, z=1.5)
def test_multiple_populations(self):
"""
Tests for 2 populations of galaxies
"""
# no covariance contributions
cosmo = coffe.Coffe(
has_density=True,
has_rsd=True,
number_density1=[1e-5],
number_density2=[1e-5],
z_mean=[1.0],
deltaz=[0.1],
fsky=[0.2],
pixelsize=[50],
sep=np.arange(50, 350, 50),
l=[0, 2, 4],
covariance_cosmic=False,
covariance_mixed=False,
covariance_poisson=False,
)
k, pk = np.transpose(np.loadtxt('PkL_CLASS.dat'))
cosmo.set_power_spectrum_linear(k, pk)
cov = covariance_matrix(cosmo.compute_covariance_bulk())
assert np.allclose(cov, 0)
# only Poisson contribution
cosmo.covariance_poisson = True
k, pk = np.transpose(np.loadtxt('PkL_CLASS.dat'))
cosmo.set_power_spectrum_linear(k, pk)
cov = covariance_matrix(cosmo.compute_covariance_bulk())
# the diagonal should not be zero
assert not np.allclose(np.diag(cov), 0)
# for 4 different poulations, covariance with mixed and Poisson
# contributions only should be zero
cosmo.set_parameters(
covariance_populations=[1, 2, 3, 4],
covariance_poisson=True,
covariance_mixed=True,
covariance_cosmic=False,
)
cov = covariance_matrix(cosmo.compute_covariance_bulk())
assert np.allclose(cov, 0)
class TestRepresentation:
"""
......
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