Skip to content
GitLab
Commit fbf06621 authored by Goran Jelic-Cizmek's avatar Goran Jelic-Cizmek
Browse files

implemented covariance for 2 populations

parent ae2d6df5
Hide whitespace changes
Inline Side-by-side
ccoffe.pxd
View file @ fbf06621
......@@ -198,9 +198,21 @@ cdef extern from "common.h":
int has_cuba
double *density
double *density1
size_t density_len
size_t density1_len
double *density2
size_t density2_len
int covariance_pop1
int covariance_pop2
int covariance_pop3
int covariance_pop4
double *pixelsize
......
coffe.pyx
View file @ fbf06621
......@@ -894,16 +894,16 @@ cdef class Coffe:
@property
def number_density(self):
def number_density1(self):
"""
Number density of tracers (in Mpc^3/h^3) at z_mean.
Number density of first tracers (in Mpc^3/h^3) at z_mean.
"""
return np.array(
[self._parameters.density[i] for i in range(self._parameters.density_len)]
[self._parameters.density1[i] for i in range(self._parameters.density1_len)]
)
@number_density.setter
def number_density(self, value):
@number_density1.setter
def number_density1(self, value):
try:
_ = iter(value)
except TypeError as err:
......@@ -917,18 +917,81 @@ cdef class Coffe:
if not all(_>= 0 and _ < 10 for _ in temp):
raise ValueError
if self._parameters.density_len:
free(self._parameters.density)
if self._parameters.density1_len:
free(self._parameters.density1)
self._parameters.density_len = len(temp)
self._parameters.density = <double *> malloc(sizeof(double) * len(temp))
for i in range(self._parameters.density_len):
self._parameters.density[i] = temp[i]
self._parameters.density1_len = len(temp)
self._parameters.density1 = <double *> malloc(sizeof(double) * len(temp))
for i in range(self._parameters.density1_len):
self._parameters.density1[i] = temp[i]
self._free_corrfunc()
self._free_multipoles()
self._free_covariance_multipoles()
@property
def number_density2(self):
"""
Number density of second tracers (in Mpc^3/h^3) at z_mean.
"""
return np.array(
[self._parameters.density2[i] for i in range(self._parameters.density2_len)]
)
@number_density2.setter
def number_density2(self, value):
try:
_ = iter(value)
except TypeError as err:
raise TypeError(f'The value {value} is not iterable') from err
try:
temp = [float(_) for _ in value]
except TypeError as err:
raise TypeError('Cannot convert all values to floats') from err
if not all(_>= 0 and _ < 10 for _ in temp):
raise ValueError
if self._parameters.density2_len:
free(self._parameters.density2)
self._parameters.density2_len = len(temp)
self._parameters.density2 = <double *> malloc(sizeof(double) * len(temp))
for i in range(self._parameters.density2_len):
self._parameters.density2[i] = temp[i]
self._free_corrfunc()
self._free_multipoles()
self._free_covariance_multipoles()
@property
def covariance_populations(self):
"""
Returns the populations used for the covariance matrix.
"""
return np.array([
self._parameters.covariance_pop1,
self._parameters.covariance_pop2,
self._parameters.covariance_pop3,
self._parameters.covariance_pop4,
])
@covariance_populations.setter
def covariance_populations(self, value):
try:
_ = iter(value)
except TypeError as err:
raise TypeError(f'The value {value} is not iterable') from err
try:
temp = [int(_) for _ in value]
except TypeError as err:
raise TypeError('Cannot convert all values to integers') from err
self._parameters.covariance_pop1, self._parameters.covariance_pop2, self._parameters.covariance_pop3, self._parameters.covariance_pop4 = value
@property
def pixelsize(self):
"""
......@@ -1680,7 +1743,7 @@ cdef class Coffe:
if not all(
len(self.z_mean) == len(_) \
for _ in (self.deltaz, self.number_density, self.pixelsize, self.fsky)
for _ in (self.deltaz, self.number_density1, self.number_density2, self.pixelsize, self.fsky)
):
raise ValueError('Mismatching lengths for covariance parameters')
......
src/common.c
View file @ fbf06621
......@@ -27,6 +27,7 @@
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_sf_coupling.h>
#ifdef HAVE_CLASS
#include "class.h"
......@@ -770,9 +771,13 @@ int coffe_parameters_free(
free(par->fsky);
par->fsky_len = 0;
if (par->density_len)
free(par->density);
par->density_len = 0;
if (par->density1_len)
free(par->density1);
par->density1_len = 0;
if (par->density2_len)
free(par->density2);
par->density2_len = 0;
if (par->pixelsize_len)
free(par->pixelsize);
......@@ -1369,3 +1374,59 @@ int coffe_free_class_struct(
return EXIT_SUCCESS;
}
/**
computes the integral over 4 Legendre polynomials of degrees n, m, a, and b
**/
double coffe_legendre_integral(
int n,
int m,
int a,
int b
)
{
/* shortcut: if n + m + a + b is odd, the integral is identically zero */
if ((n + m + a + b) % 2 != 0)
return 0;
const int lower_limit = (abs(n - m) > abs(a - b) ? abs(n - m) : abs(a - b));
const int upper_limit = (n + m < a + b ? n + m : a + b);
double result = 0;
for (int l = lower_limit; l <= upper_limit; ++l){
result +=
(2. * l + 1)
*pow(gsl_sf_coupling_3j(2 * n, 2 * m, 2 * l, 0, 0, 0), 2)
*pow(gsl_sf_coupling_3j(2 * a, 2 * b, 2 * l, 0, 0, 0), 2);
}
return 2 * result;
}
/**
calculates the Kronecker delta
**/
int coffe_kronecker_delta(
int i,
int j
)
{
return (i == j ? 1 : 0);
}
/**
calculates (-1)^m
**/
int coffe_sign(
int m
)
{
return (abs(m) % 2 == 0 ? 1 : -1);
}
src/common.h
View file @ fbf06621
......@@ -431,14 +431,18 @@ typedef struct coffe_parameters_t
size_t fsky_len;
double *density;
double *density1, *density2;
size_t density_len;
size_t density1_len, density2_len;
double *pixelsize;
size_t pixelsize_len;
/* the covariance is a 4-point function, hence we can in general
specify 4 different populations */
int covariance_pop1, covariance_pop2, covariance_pop3, covariance_pop4;
int covariance_integration_method;
int covariance_integration_bins;
......@@ -802,4 +806,21 @@ int coffe_free_fit_coefficients_array(
);
double coffe_legendre_integral(
int n,
int m,
int a,
int b
);
int coffe_kronecker_delta(
int i,
int j
);
int coffe_sign(int m);
#endif
src/covariance.c
View file @ fbf06621
......@@ -32,6 +32,31 @@
#endif
/**
computes the (flat-sky) coefficients for the covariance of the multipoles.
They are just combinations of the galaxy bias and the growth rate.
**/
double covariance_coefficient(
double bias1,
double bias2,
double growth_rate,
int ell
)
{
switch (ell){
case 0:
return bias1 * bias2 + (bias1 + bias2) * growth_rate / 3 + pow(growth_rate, 2) / 5.;
case 2:
return 2 * (bias1 + bias2) * growth_rate / 3. + 4 * pow(growth_rate, 2) / 7.;
case 4:
return 8 * pow(growth_rate, 2) / 35.;
default:
return 0;
}
}
/**
find the element in the array which corresponds to some values of the parameters (z_mean, l1, l2, sep1, sep2)
TODO what if we can't find anything?
......@@ -132,17 +157,6 @@ static inline double covariance_volume_no_4pi(
}
/**
calculates i^(l1 - l2) when l1 - l2 is even
**/
static int covariance_complex(int l1, int l2)
{
if ((l1 - l2) % 4 == 0) return 1;
else if ((l1 - l2) % 2 == 0) return -1;
else return 0;
}
/**
integrand of the windowed spherical Bessel function
**/
......@@ -723,10 +737,10 @@ int coffe_covariance_init(
}}}}}
double *volume =
(double *)coffe_malloc(sizeof(double)*par->density_len);
(double *)coffe_malloc(sizeof(double)*par->density1_len);
size_t index = 0;
double z_mean = 0;
for (size_t k = 0; k < par->density_len; ++k){
for (size_t k = 0; k < par->density1_len; ++k){
if (par->output_type == COVARIANCE_MULTIPOLES){
z_mean = par->z_mean[k];
......@@ -763,119 +777,291 @@ int coffe_covariance_init(
*/
}
/* TODO implement covariance for two populations */
double galaxy_bias1 = 0, growth_rate = 0;
double galaxy_bias1 = 0;
double galaxy_bias2 = 0;
double growth_rate = 0;
/* set bias to nonzero only if there is a density contribution */
if (par->correlation_contrib.den)
if (par->correlation_contrib.den){
galaxy_bias1 = coffe_interp_spline(&par->galaxy_bias1, z_mean);
galaxy_bias2 = coffe_interp_spline(&par->galaxy_bias2, z_mean);
}
/* set growth rate to nonzero only if there's an rsd contribution */
if (par->correlation_contrib.rsd)
growth_rate = coffe_interp_spline(&bg->f, z_mean);
/* b^2 + 2/3 b growth_rate + growth_rate^2/5 */
const double c0 = pow(galaxy_bias1, 2) + 2 * galaxy_bias1 * growth_rate / 3. + pow(growth_rate, 2) / 5.;
/* 4/3 b growth_rate + 4/7 growth_rate^2 */
const double c2 = 4 * galaxy_bias1 * growth_rate / 3. + 4 * pow(growth_rate, 2) / 7.;
/* 8/35 growth_rate^2 */
const double c4 = 8 * pow(growth_rate, 2) / 35.;
const double c0bar =
c0*c0 + c2*c2/5. + c4*c4/9.;
const double c2bar =
2*c2*(7*c0 + c2)/7. + 4*c2*c4/7. + 100*c4*c4/693.;
const double c4bar =
18*c2*c2/35. + 2*c0*c4 + 40*c2*c4/77. + 162*c4*c4/1001.;
const double c6bar =
10*c4*(9*c2 + 2*c4)/99.;
const double c8bar =
490*c4*c4/1287.;
const double alpha0_11 = covariance_coefficient(
galaxy_bias1,
galaxy_bias1,
growth_rate,
0
);
const double D1z = coffe_interp_spline(&bg->D1, z_mean);
const double alpha0_22 = covariance_coefficient(
galaxy_bias2,
galaxy_bias2,
growth_rate,
0
);
const double coeff_array[] = {c0, c2, c4};
const double coeffbar_array[] = {c0bar, c2bar, c4bar, c6bar, c8bar};
const double alpha0_cross = covariance_coefficient(
galaxy_bias1,
galaxy_bias2,
growth_rate,
0
);
const double alpha2_11 = covariance_coefficient(
galaxy_bias1,
galaxy_bias1,
growth_rate,
2
);
for (size_t i = 0; i<par->multipole_values_len; ++i){
const double alpha2_22 = covariance_coefficient(
galaxy_bias2,
galaxy_bias2,
growth_rate,
2
);
const double alpha2_cross = covariance_coefficient(
galaxy_bias1,
galaxy_bias2,
growth_rate,
2
);
const double alpha4 = covariance_coefficient(
galaxy_bias1,
galaxy_bias2,
growth_rate,
4
);
for (size_t m = 0; m < par->sep_len; ++m){
for (size_t j = 0; j<par->multipole_values_len; ++j){
const double D1z = coffe_interp_spline(&bg->D1, z_mean);
for (size_t i = 0; i < par->multipole_values_len; ++i){
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);
const int l1 = par->multipole_values[i];
const int l2 = par->multipole_values[j];
double coeff_pop1_pop3[] = {0, 0, alpha4};
if (par->covariance_pop1 == 1 && par->covariance_pop3 == 1){
coeff_pop1_pop3[0] = alpha0_11;
coeff_pop1_pop3[1] = alpha2_11;
}
else if (
(par->covariance_pop1 == 1 && par->covariance_pop3 == 2)
||
(par->covariance_pop1 == 2 && par->covariance_pop3 == 1)
){
coeff_pop1_pop3[0] = alpha0_cross;
coeff_pop1_pop3[1] = alpha2_cross;
}
else{
coeff_pop1_pop3[0] = alpha0_22;
coeff_pop1_pop3[1] = alpha2_22;
}
double coeff_pop2_pop4[] = {0, 0, alpha4};
if (par->covariance_pop2 == 1 && par->covariance_pop4 == 1){
coeff_pop2_pop4[0] = alpha0_11;
coeff_pop2_pop4[1] = alpha2_11;
}
else if (
(par->covariance_pop2 == 1 && par->covariance_pop4 == 2)
||
(par->covariance_pop2 == 2 && par->covariance_pop4 == 1)
){
coeff_pop2_pop4[0] = alpha0_cross;
coeff_pop2_pop4[1] = alpha2_cross;
}
else{
coeff_pop2_pop4[0] = alpha0_22;
coeff_pop2_pop4[1] = alpha2_22;
}
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 coeff_pop1_pop4[] = {0, 0, alpha4};
if (par->covariance_pop1 == 1 && par->covariance_pop4 == 1){
coeff_pop1_pop4[0] = alpha0_11;
coeff_pop1_pop4[1] = alpha2_11;
}
else if (
(par->covariance_pop1 == 1 && par->covariance_pop4 == 2)
||
(par->covariance_pop1 == 2 && par->covariance_pop4 == 1)
){
coeff_pop1_pop4[0] = alpha0_cross;
coeff_pop1_pop4[1] = alpha2_cross;
}
else{
coeff_pop1_pop4[0] = alpha0_22;
coeff_pop1_pop4[1] = alpha2_22;
}
/* Kronecker delta (l1, l2) */
double deltal1l2 =
(par->multipole_values[i] == par->multipole_values[j]) ? 1 : 0;
double coeff_pop2_pop3[] = {0, 0, alpha4};
if (par->covariance_pop2 == 1 && par->covariance_pop4 == 1){
coeff_pop2_pop3[0] = alpha0_11;
coeff_pop2_pop3[1] = alpha2_11;
}
else if (
(par->covariance_pop2 == 1 && par->covariance_pop4 == 2)
||
(par->covariance_pop2 == 2 && par->covariance_pop4 == 1)
){
coeff_pop2_pop3[0] = alpha0_cross;
coeff_pop2_pop3[1] = alpha2_cross;
}
else{
coeff_pop2_pop3[0] = alpha0_22;
coeff_pop2_pop3[1] = alpha2_22;
}
double deltaij = (m == n) ? 1 : 0;
/* flat-sky covariance */
const double result_mp_or_ramp =
covariance_complex(
par->multipole_values[i],
par->multipole_values[j]
double coeff_cosmic_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_pop1_pop3[cnt1]
*coeff_pop2_pop4[cnt2]
+
coffe_sign(l2)
*coeff_pop1_pop4[cnt1]
*coeff_pop2_pop3[cnt2]
)
*(
(2 * par->multipole_values[i] + 1)
*deltal1l2
*deltaij
/ 2. / M_PI
/* trigraph */
*(
par->covariance_window
?
1.
/ covariance_volume_no_4pi(
separations[m] - par->pixelsize[k] / 2.,
separations[m] + par->pixelsize[k] / 2.
)
:
1.
/ separations[n]
/ separations[m]
/ par->pixelsize[k]
*coffe_legendre_integral(l1, l2, 2 * cnt1, 2 * cnt2);
}}
coeff_cosmic_cosmic_sum /= 4.;
double coeff_cosmic_poisson_sum = 0;
for (size_t cnt = 0; cnt < COFFE_ARRAY_SIZE(coeff_pop1_pop3); ++cnt){
coeff_cosmic_poisson_sum +=
(
coeff_pop1_pop3[cnt]
*coffe_kronecker_delta(
par->covariance_pop2,
par->covariance_pop4
)
/par->density2[k]
+
coeff_pop2_pop4[cnt]
*coffe_kronecker_delta(
par->covariance_pop1,
par->covariance_pop3
)
/par->density1[k]
+
coeff_pop1_pop4[cnt]
*coffe_kronecker_delta(
par->covariance_pop2,
par->covariance_pop3
)
/par->density[k]
/par->density[k]
/par->density2[k]
+
/* trigraph */
(par->pk_type ? 1 : D1z * D1z)
/* trigraph */
*(
par->pk_type
?
integral_pk[k][i*par->multipole_values_len + j][par->sep_len*n + m]
:
integral_pk[0][i*par->multipole_values_len + j][par->sep_len*n + m]
coeff_pop2_pop3[cnt]
*coffe_kronecker_delta(
par->covariance_pop1,
par->covariance_pop4
)
/par->density1[k]
)
*coffe_legendre_integral(l1, l2, 2 * cnt, 0);
}
coeff_cosmic_poisson_sum /= 8.;
/* Kronecker delta (l1, l2) */
const int deltal1l2 = coffe_kronecker_delta(l1, l2);
/* Kronecker delta of the separations */
const int deltaij = coffe_kronecker_delta(m, n);
const double cov_poisson_poisson =
(2 * l1 + 1)
/ 4. / M_PI
*deltal1l2