Commit a3a92f17 authored by Plaszczynski Stephane's avatar Plaszczynski Stephane
Browse files

add some optimized class pk computation untill gets into master branch

parent 28cd9156
......@@ -81,6 +81,10 @@ library MCMC -no_share -s=$(CAMELROOT)/src/camel/MCMC *.cc
macro CAMEL_linkopts "-L$(CAMELROOT)/${CMTCONFIG} -lMinuitFit -lMinuit -lMCMC -lCLHEP -lAbsRand -lUtil "
macro_append cppflags ' -DRELPATH=\"$(CAMELROOT)/lik\" '
#class extras
library class_extra -no_share $(CAMELROOT)/src/class_extra/*.c
macro_append CAMEL_linkopts " -lclass_extra "
#CAMEL applications
# pour relinker si changement: toutes les application du package
......
......@@ -720,21 +720,10 @@ double ClassEngine::get_DMod(double z){
}
#ifdef FASTPK
//fast access
/*
bool
ClassEngine::get_Pkvec(const std::vector<double>& knodes,
std::vector<double>& pks){
planck_assert(sp.ln_tau_size==1,name()+":get_Pkvec cannot proceed since several z values required-> remove z_max_pk/z_pk or use get_Pkvec(k,z,pk)" );
return (spectra_fast_pk_at_kvec(
&ba,
&sp,
const_cast<double*>(&knodes[0]),
knodes.size(),
&pks[0])==_TRUE_);
}
*/
#include "class_extra/spectra_extra.h"
bool
ClassEngine::get_Pkvec(const std::vector<double>& knodes, const std::vector<double>& zvec, std::vector<double>& pks)
......
#include "spectra.h"
/* extra stuf */
int spectra_fast_pk_at_kvec_and_zvec(
struct background * pba,
struct spectra * psp,
double * kvec,
int kvec_size,
double * zvec,
int zvec_size,
double * pk_tot_out /* (must be already allocated with kvec_size*zvec_size) */
) {
/** Summary: */
/** - define local variables */
int index_md;
int index_k, index_knode, index_z;
double *spline, *pk_at_k, *ln_kvec, *ln_pk_table;
double ln_pk_interp, ln_k_interp;
double h, a, b;
index_md = psp->index_md_scalars;
class_test(psp->ic_size[index_md] != 1,
psp->error_message,
"This function has only been coded for pure adiabatic ICs, sorry.");
/** Compute spline over ln(k) */
class_alloc(ln_kvec, sizeof(double)*kvec_size,psp->error_message);
class_alloc(ln_pk_table, sizeof(double)*psp->ln_k_size*zvec_size,psp->error_message);
class_alloc(spline, sizeof(double)*psp->ln_k_size*zvec_size,psp->error_message);
class_alloc(pk_at_k, sizeof(double)*psp->ln_tau_size,psp->error_message);
/** Construct table of log(pk) on the computed k nodes but requested redshifts: */
for (index_z=0; index_z<zvec_size; index_z++){
class_call(spectra_pk_at_z(pba,psp, logarithmic,zvec[index_z],ln_pk_table+index_z*psp->ln_k_size,NULL),
psp->error_message,
psp->error_message);
}
class_call(array_spline_table_columns2(psp->ln_k,
psp->ln_k_size,
ln_pk_table,
zvec_size,
spline,
_SPLINE_NATURAL_,
psp->error_message),
psp->error_message,
psp->error_message);
/** Construct ln(kvec): */
for (index_k=0; index_k<kvec_size; index_k++){
ln_kvec[index_k] = log(kvec[index_k]);
}
/** I will assume that the k vector is sorted in ascending order.
Case k<kmin: */
for(index_k = 0; index_k<kvec_size; index_k++){
if (ln_kvec[index_k] >= psp->ln_k[0])
break;
for (index_z = 0; index_z < zvec_size; index_z++) {
/** If needed, add some extrapolation here */
pk_tot_out[index_z*kvec_size+index_k] = 0.;
}
/** Implement some extrapolation perhaps */
}
/** Case kmin<=k<=kmax. Do not loop through kvec, but loop through the
interpolation nodes. */
for (index_knode=0; index_knode < (psp->ln_k_size-1); index_knode++){
/** Loop through k's that fall in this interval */
//printf("index _k is %d, do we have %g < %g <%g?\n",index_k, psp->ln_k[index_knode],ln_kvec[index_k],psp->ln_k[index_knode+1]);
while ((index_k < kvec_size) && (ln_kvec[index_k] <= psp->ln_k[index_knode+1])){
/** Perform interpolation */
h = psp->ln_k[index_knode+1]-psp->ln_k[index_knode];
b = (ln_kvec[index_k] - psp->ln_k[index_knode])/h;
a = 1.-b;
for (index_z = 0; index_z < zvec_size; index_z++) {
pk_tot_out[index_z*kvec_size+index_k] =
exp(
a * ln_pk_table[index_z*psp->ln_k_size + index_knode]
+ b * ln_pk_table[index_z*psp->ln_k_size + index_knode+1]
+ ((a*a*a-a) * spline[index_z*psp->ln_k_size + index_knode]
+(b*b*b-b) * spline[index_z*psp->ln_k_size + index_knode+1])*h*h/6.0
);
}
index_k++;
}
}
/** case k>kmax */
while (index_k<kvec_size){
for (index_z = 0; index_z < zvec_size; index_z++) {
/** If needed, add some extrapolation here */
pk_tot_out[index_z*kvec_size+index_k] = 0.;
}
index_k++;
}
free(ln_kvec);
free(ln_pk_table);
free(spline);
free(pk_at_k);
return _SUCCESS_;
}
int spectra_fast_pk_nl_at_kvec_and_zvec(
struct background * pba,
struct spectra * psp,
double * kvec,
int kvec_size,
double * zvec,
int zvec_size,
double * pk_tot_out /* (must be already allocated with kvec_size*zvec_size) */
) {
/** Summary: */
/** - define local variables */
int index_md;
int index_k, index_knode, index_z;
double *spline, *pk_at_k, *ln_kvec, *ln_pk_table;
double ln_pk_interp, ln_k_interp;
double h, a, b;
index_md = psp->index_md_scalars;
class_test(psp->ic_size[index_md] != 1,
psp->error_message,
"This function has only been coded for pure adiabatic ICs, sorry.");
/** Compute spline over ln(k) */
class_alloc(ln_kvec, sizeof(double)*kvec_size,psp->error_message);
class_alloc(ln_pk_table, sizeof(double)*psp->ln_k_size*zvec_size,psp->error_message);
class_alloc(spline, sizeof(double)*psp->ln_k_size*zvec_size,psp->error_message);
class_alloc(pk_at_k, sizeof(double)*psp->ln_tau_size,psp->error_message);
/** Construct table of log(pk) on the computed k nodes but requested redshifts: */
for (index_z=0; index_z<zvec_size; index_z++){
class_call(spectra_pk_nl_at_z(pba,psp, logarithmic,zvec[index_z],ln_pk_table+index_z*psp->ln_k_size),
psp->error_message,
psp->error_message);
}
class_call(array_spline_table_columns2(psp->ln_k,
psp->ln_k_size,
ln_pk_table,
zvec_size,
spline,
_SPLINE_NATURAL_,
psp->error_message),
psp->error_message,
psp->error_message);
/** Construct ln(kvec): */
for (index_k=0; index_k<kvec_size; index_k++){
ln_kvec[index_k] = log(kvec[index_k]);
}
/** I will assume that the k vector is sorted in ascending order.
Case k<kmin: */
for(index_k = 0; index_k<kvec_size; index_k++){
if (ln_kvec[index_k] >= psp->ln_k[0])
break;
for (index_z = 0; index_z < zvec_size; index_z++) {
/** If needed, add some extrapolation here */
pk_tot_out[index_z*kvec_size+index_k] = 0.;
}
/** Implement some extrapolation perhaps */
}
/** Case kmin<=k<=kmax. Do not loop through kvec, but loop through the
interpolation nodes. */
for (index_knode=0; index_knode < (psp->ln_k_size-1); index_knode++){
/** Loop through k's that fall in this interval */
//printf("index _k is %d, do we have %g < %g <%g?\n",index_k, psp->ln_k[index_knode],ln_kvec[index_k],psp->ln_k[index_knode+1]);
while ((index_k < kvec_size) && (ln_kvec[index_k] <= psp->ln_k[index_knode+1])){
/** Perform interpolation */
h = psp->ln_k[index_knode+1]-psp->ln_k[index_knode];
b = (ln_kvec[index_k] - psp->ln_k[index_knode])/h;
a = 1.-b;
for (index_z = 0; index_z < zvec_size; index_z++) {
pk_tot_out[index_z*kvec_size+index_k] =
exp(
a * ln_pk_table[index_z*psp->ln_k_size + index_knode]
+ b * ln_pk_table[index_z*psp->ln_k_size + index_knode+1]
+ ((a*a*a-a) * spline[index_z*psp->ln_k_size + index_knode]
+(b*b*b-b) * spline[index_z*psp->ln_k_size + index_knode+1])*h*h/6.0
);
}
index_k++;
}
}
/** case k>kmax */
while (index_k<kvec_size){
for (index_z = 0; index_z < zvec_size; index_z++) {
/** If needed, add some extrapolation here */
pk_tot_out[index_z*kvec_size+index_k] = 0.;
}
index_k++;
}
free(ln_kvec);
free(ln_pk_table);
free(spline);
free(pk_at_k);
return _SUCCESS_;
}
/** @file spectra.h Documented includes for spectra module */
#ifndef __SPECTRA__EXTRA
#define __SPECTRA__EXTRA
#include "transfer.h"
#include "spectra.h"
#ifdef __cplusplus
extern "C" {
#endif
/* extra stuf */
int spectra_fast_pk_at_kvec_and_zvec(
struct background * pba,
struct spectra * psp,
double * kvec,
int kvec_size,
double * zvec,
int zvec_size,
double * pk_tot_out);
int spectra_fast_pk_nl_at_kvec_and_zvec(
struct background * pba,
struct spectra * psp,
double * kvec,
int kvec_size,
double * zvec,
int zvec_size,
double * pk_tot_out);
#ifdef __cplusplus
}
#endif
#endif
/* @endcond */
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