Commit 5908e1c3 authored by lesgourg's avatar lesgourg
Browse files

several bug and memory leak fixes

* 2.6:
  changed version number to 2.6.1
  corrected default value of selection_magnification_bias=0
  corrected a factor in the definition of C_l^shearshear
  removed a memory leak in functions raw_cl ans lensed_cl
  removed a memory leak in functions pk and pk_lin of the wrapper
  added to the wrapper new functions f() and sigma()
  corrected bugs affecting only P_NL(k) in presence of isocurvature modes
parents 992b18b4 ffe71e10
......@@ -347,7 +347,7 @@ Cls or the shear/convergence Cls. The relations between them are trivial:
--> deflection d:
Cl^dd = l(l+1) C_l^phiphi
--> convergence kappa and shear gamma: the share the same harmonic power spectrum:
Cl^gamma-gamma = 1/4 * [l(l+1)]^2 C_l^phi-phi
Cl^gamma-gamma = 1/4 * [(l+2)!/(l-2)!] C_l^phi-phi
By defaut, the code will try to compute the following cross-correlation Cls (if
available): temperature-polarisation, temperature-CMB lensing, polarization-CMB
lensing, CMB lensing-density, and density-lensing. Other cross-correlations are
......
......@@ -15,7 +15,7 @@
#ifndef __COMMON__
#define __COMMON__
#define _VERSION_ "v2.6.0"
#define _VERSION_ "v2.6.1"
/* @cond INCLUDE_WITH_DOXYGEN */
#define _TRUE_ 1 /**< integer associated to true statement */
......
......@@ -38,6 +38,7 @@ cdef extern from "class.h":
int index_bg_conf_distance
int index_bg_H
int index_bg_D
int index_bg_f
short long_info
short inter_normal
double T_cmb
......@@ -95,6 +96,7 @@ cdef extern from "class.h":
int store_perturbations
int k_output_values_num
double k_output_values[_MAX_NUMBER_OF_K_FILES_]
double k_max_for_pk
int index_k_output_values[_MAX_NUMBER_OF_K_FILES_]
char scalar_titles[_MAXTITLESTRINGLENGTH_]
char vector_titles[_MAXTITLESTRINGLENGTH_]
......@@ -320,3 +322,11 @@ cdef extern from "class.h":
int nonlinear_k_nl_at_z(void* pba, void* pnl, double z, double* k_nl)
int spectra_firstline_and_ic_suffix(void *ppt, int index_ic, char first_line[_LINE_LENGTH_MAX_], FileName ic_suffix)
int spectra_sigma(
void * pba,
void * ppm,
void * psp,
double R,
double z,
double * sigma)
\ No newline at end of file
......@@ -8,6 +8,8 @@
This module defines a class called Class. It is used with Monte Python to
extract cosmological parameters.
# JL 14.06.2017: TODO: check whether we should free somewhere the allocated fc.filename and titles, data (4 times)
"""
from math import exp,log
import numpy as np
......@@ -465,6 +467,12 @@ cdef class Class:
cl['ell'] = np.arange(lmax+1)
free(rcl)
for index_md in range(self.sp.md_size):
free(cl_md[index_md])
free(cl_md_ic[index_md])
free(cl_md)
free(cl_md_ic)
return cl
def lensed_cl(self, lmax=-1,nofail=False):
......@@ -636,6 +644,12 @@ cdef class Class:
cl['ell'] = np.arange(lmax+1)
free(dcl)
for index_md in range(self.sp.md_size):
free(cl_md[index_md])
free(cl_md_ic[index_md])
free(cl_md)
free(cl_md_ic)
return cl
def z_of_r (self,z_array):
......@@ -712,6 +726,8 @@ cdef class Class:
else:
if spectra_pk_nl_at_k_and_z(&self.ba,&self.pm,&self.sp,k,z,&pk) ==_FAILURE_:
raise CosmoSevereError(self.sp.error_message)
free(pk_ic)
return pk
# Gives the linear pk for a given (k,z)
......@@ -740,6 +756,7 @@ cdef class Class:
if spectra_pk_at_k_and_z(&self.ba,&self.pm,&self.sp,k,z,&pk,pk_ic)==_FAILURE_:
raise CosmoSevereError(self.sp.error_message)
free(pk_ic)
return pk
def get_pk(self, np.ndarray[DTYPE_t,ndim=3] k, np.ndarray[DTYPE_t,ndim=1] z, int k_size, int z_size, int mu_size):
......@@ -764,6 +781,36 @@ cdef class Class:
pk[index_k,index_z,index_mu] = self.pk_lin(k[index_k,index_z,index_mu],z[index_z])
return pk
# Gives sigma(R,z) for a given (R,z)
def sigma(self,double R,double z):
"""
Gives the pk for a given R and z
(R is the radius in units of Mpc, so if R=8/h this will be the usual sigma8(z)
.. note::
there is an additional check to verify whether output contains `mPk`,
and whether k_max > ...
because otherwise a segfault will occur
"""
cdef double sigma
if (self.pt.has_pk_matter == _FALSE_):
raise CosmoSevereError(
"No power spectrum computed. In order to get sigma(R,z) you must add mPk to the list of outputs."
)
if (self.pt.k_max_for_pk < self.ba.h):
raise CosmoSevereError(
"In order to get sigma(R,z) you must set 'P_k_max_h/Mpc' to 1 or bigger, in order to have k_max > 1 h/Mpc."
)
if spectra_sigma(&self.ba,&self.pm,&self.sp,R,z,&sigma)==_FAILURE_:
raise CosmoSevereError(self.sp.error_message)
return sigma
def age(self):
self.compute(["background"])
return self.ba.age
......@@ -855,6 +902,36 @@ cdef class Class:
return D
def scale_independent_growth_factor_f(self, z):
"""
scale_independent_growth_factor_f(z)
Return the scale invariant growth factor f(z)=d ln D / d ln a for CDM perturbations
(exactly, the quantity defined by Class as index_bg_f in the background module)
Parameters
----------
z : float
Desired redshift
"""
cdef double tau
cdef int last_index #junk
cdef double * pvecback
pvecback = <double*> calloc(self.ba.bg_size,sizeof(double))
if background_tau_of_z(&self.ba,z,&tau)==_FAILURE_:
raise CosmoSevereError(self.ba.error_message)
if background_at_tau(&self.ba,tau,self.ba.long_info,self.ba.inter_normal,&last_index,pvecback)==_FAILURE_:
raise CosmoSevereError(self.ba.error_message)
f = pvecback[self.ba.index_bg_f]
free(pvecback)
return f
def Hubble(self, z):
"""
Hubble(z)
......
......@@ -3088,7 +3088,7 @@ int input_default_params(
/** - transfer structure */
ptr->selection_bias[0]=1.;
ptr->selection_magnification_bias[0]=1.;
ptr->selection_magnification_bias[0]=0.;
ptr->lcmb_rescale=1.;
ptr->lcmb_pivot=0.1;
ptr->lcmb_tilt=0.;
......
......@@ -464,7 +464,7 @@ int spectra_pk_at_z(
/** - fourth step: depending on requested mode (linear or logarithmic), apply necessary transformation to the output arrays */
/** - --> (a) linear mode: if only one initial condition, convert output_pk to linear format; if several initial conditions, convert output_ic to linear format, output_tot is already in this format */
/** - --> (a) linear mode: if only one initial condition, convert output_tot to linear format; if several initial conditions, convert output_ic to linear format, output_tot is already in this format */
if (mode == linear) {
......@@ -2551,13 +2551,13 @@ int spectra_pk(
/** - define local variables */
int index_md;
int index_ic1,index_ic2,index_ic1_ic2;
int index_ic1,index_ic2,index_ic1_ic1,index_ic2_ic2,index_ic1_ic2;
int index_k;
int index_tau;
double * primordial_pk; /* array with argument primordial_pk[index_ic_ic] */
double source_ic1;
double source_ic2;
double ln_pk_tot;
double pk_tot=0.,ln_pk_tot=0.;
/** - check the presence of scalar modes */
......@@ -2593,7 +2593,7 @@ int spectra_pk(
ppm->error_message,
psp->error_message);
ln_pk_tot =0;
pk_tot =0;
/* curvature primordial spectrum:
P_R(k) = 1/(2pi^2) k^3 <R R>
......@@ -2619,7 +2619,7 @@ int spectra_pk(
*source_ic1*source_ic1
*exp(primordial_pk[index_ic1_ic2]));
ln_pk_tot += psp->ln_pk[(index_tau * psp->ln_k_size + index_k)* psp->ic_ic_size[index_md] + index_ic1_ic2];
pk_tot += exp(psp->ln_pk[(index_tau * psp->ln_k_size + index_k)* psp->ic_ic_size[index_md] + index_ic1_ic2]);
}
......@@ -2628,6 +2628,8 @@ int spectra_pk(
for (index_ic2 = index_ic1+1; index_ic2 < psp->ic_size[index_md]; index_ic2++) {
index_ic1_ic2 = index_symmetric_matrix(index_ic1,index_ic2,psp->ic_size[index_md]);
index_ic1_ic1 = index_symmetric_matrix(index_ic1,index_ic1,psp->ic_size[index_md]);
index_ic2_ic2 = index_symmetric_matrix(index_ic2,index_ic2,psp->ic_size[index_md]);
if (psp->is_non_zero[index_md][index_ic1_ic2] == _TRUE_) {
......@@ -2642,7 +2644,10 @@ int spectra_pk(
psp->ln_pk[(index_tau * psp->ln_k_size + index_k)* psp->ic_ic_size[index_md] + index_ic1_ic2] =
primordial_pk[index_ic1_ic2]*SIGN(source_ic1)*SIGN(source_ic2);
ln_pk_tot += psp->ln_pk[(index_tau * psp->ln_k_size + index_k)* psp->ic_ic_size[index_md] + index_ic1_ic2];
pk_tot += psp->ln_pk[(index_tau * psp->ln_k_size + index_k)* psp->ic_ic_size[index_md] + index_ic1_ic2]
* sqrt(psp->ln_pk[(index_tau * psp->ln_k_size + index_k)* psp->ic_ic_size[index_md] + index_ic1_ic1]
* psp->ln_pk[(index_tau * psp->ln_k_size + index_k)* psp->ic_ic_size[index_md] + index_ic2_ic2]);
}
else {
......@@ -2651,6 +2656,8 @@ int spectra_pk(
}
}
ln_pk_tot = log(pk_tot);
/* if non-linear corrections required, compute the total non-linear matter power spectrum */
if (pnl->method != nl_none) {
......
Markdown is supported
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