Docker-in-Docker (DinD) capabilities of public runners deactivated. More info

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

add optiopnal fast Pk access

parent 7e8fa522
......@@ -718,3 +718,23 @@ double ClassEngine::get_DMod(double z){
return(mu);
}
#ifdef FASTPK
//fast access
void
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)" );
spectra_fast_pk_at_kvec(
&ba,
&sp,
const_cast<double*>(&knodes[0]),
knodes.size(),
&pks[0]);
}
#endif
......@@ -120,7 +120,7 @@ public:
//get value at l ( 2<l<lmax): in units = (micro-K)^2
//don't call if FAILURE returned previously
//throws std::execption if pb
virtual std::string name() const {return "CLASS";}
virtual std::string name() const {return "ClassEngine";}
double getCl(Engine::cltype t,const long &l);
void getCls(const std::vector<unsigned>& lVec, //input
......@@ -155,6 +155,13 @@ public:
double get_Pklin(double k, double z);
double get_PkNL(double k, double z);
#ifdef FASTPK
//should be used with care only for z=0 so no z_max_pk and z_pk should be given
//pks should be sized as knodes
void get_Pkvec(const std::vector<double>& knodes,
std::vector<double>& pks);
#endif
//combile distance
double com_distance(double z);
......@@ -233,6 +240,8 @@ protected:
std::string _output;
double* pvecback; // helper to avoid too much allocations
};
#endif
......
......@@ -45,7 +45,7 @@ public:
bool updateParValues(const std::vector<double>& minuitPars);
virtual std::string name() const {return "Minuit-CLASS";}
virtual std::string name() const {return "MnClassEngine";}
//then read with this one
inline void getCls(const std::vector<unsigned>& lVec, //input
......@@ -71,6 +71,7 @@ public:
inline double z_drag() const {return klass->z_drag();}
inline double rs_drag() const {return klass->rs_drag();}
double get_Dv(double z) {return klass->get_Dv(z);}
......@@ -94,6 +95,10 @@ public:
virtual double get_Pklin(double k, double z) { return klass->get_Pklin(k,z);}
virtual double get_PkNL(double k, double z) { return klass->get_PkNL(k,z);}
inline void get_Pkvec(const std::vector<double>& knodes,
std::vector<double>& pks)
{ return klass->get_Pkvec(knodes,pks);}
void writeCls(std::ostream &o) {klass->writeCls(o);}
......
......@@ -93,12 +93,14 @@ public:
virtual double get_Pklin(double k, double z) { return undef("get_Pklin");}
virtual double get_PkNL(double k, double z) { return undef("get_Pknl");}
//access through k input vector. output vector pks should be sized the same
// P(k,z=0) access through vectors for given k values.
// output vector pks should be sized the same
//default implementation: calls P(k,z) -> override for better performances
inline virtual
void get_Pks(const std::vector<double>& k, std::vector<double>& pks,double z){
for (size_t i=0;i<k.size();i++) pks[i]=get_Pk(k[i],z);
inline virtual void get_Pkvec(const std::vector<double>& knodes,
std::vector<double>& pks){
for (size_t i=0;i<knodes.size();i++) pks[i]=get_Pk(knodes[i],0);
}
......
......@@ -3,6 +3,7 @@
#ifdef PICO
#include "pico/MnPicoEngine.hh"
#endif
#include "Timer.hh"
#include"cxxsupport/arr.h"
#include"cxxsupport/fitshandle.h"
#include <iostream>
......@@ -38,6 +39,7 @@ int main(int argc,char** argv){
float z0 = 0.;
if (argc==5) z0=atof(argv[4]);
//decodage arguments
Parser parser(argv[1]);
int lmax=atoi(argv[2]);
......@@ -56,25 +58,20 @@ int main(int argc,char** argv){
classparms.add("output","tCl,pCl,lCl,mPk"); //polar +lens+clphi
classparms.add("z_pk",z0); //polar +lens+clphi
classparms.add("z_max_pk",z0); //polar +lens+clphi
//classparms.add("z_max_pk",z0); //polar +lens+clphi
//increase kmax if higher value requested
const double kmax = parser.params.find<double>("k_end",1);
for (size_t i=0;i<classparms.size();i++) {
if (classparms.key(i)=="P_k_max_1/Mpc"){
double val;
stringToData(classparms.value(i),val);
if ( val < kmax) {
cout << "Warning: requested k_end=" <<kmax<< " while P_k_max_1/Mpc=" << classparms.value(i) << " increasing parameter" <<endl;
classparms.updateVal(i, dataToString(kmax));
}
}
cout << "CLASS \t--> "<< classparms.key(i) << "\t" << classparms.value(i) <<endl;
try{
size_t i=classparms.findIndex("P_k_max_1/Mpc");
classparms.updateVal(i,dataToString(kmax));
}
catch (Message_error&){
classparms.add("P_k_max_1/Mpc",kmax);
}
for (size_t i=0;i<classparms.size();i++)
cout << "CLASS \t--> "<< classparms.key(i) << "\t" << classparms.value(i) <<endl;
//precision file
string pre="";
......@@ -201,6 +198,7 @@ int main(int argc,char** argv){
pkcols.push_back(fitscolumn("k", "1/Mpc",1,TDOUBLE));
pkcols.push_back(fitscolumn("Pklin", "(1/Mpc)^-3" ,1,TDOUBLE));
if (do_nl) pkcols.push_back(fitscolumn("Pknl", "(1/Mpc)^-3",1,TDOUBLE));
pkcols.push_back(fitscolumn("Pkfast", "(1/Mpc)^-3",1,TDOUBLE));
fout.insert_bintab(pkcols);
// def values
......@@ -228,31 +226,53 @@ int main(int argc,char** argv){
const int npks = parser.params.find<int>("n_pks",1000);
const double kmin = parser.params.find<double>("k_beg",8.e-4);
logspace(kmin,kmax,npks,kval) ;
}
//test fast Pk
Timer timer;
cout <<"Pk timer test..." << endl;
vector<double> pkvec(kval.size());
for (size_t i=0;i<1000;i++)
e->get_Pkvec(kval,pkvec);
cout << "fast=" << timer.partial() <<endl;
for (size_t i=0;i<1000;i++){
for (int k =0; k<kval.size() ; k++){
double pk=e->get_Pklin(kval[k], 0);
}
}
cout << "std=" << timer.partial() <<endl;
vector<double> pklin;
vector<double> pknl ;
for (int k =0; k<kval.size() ; k++){
double myk=kval[k];
pklin.push_back( e->get_Pklin(myk, z0) );
if ( do_nl ) pknl.push_back( e->get_PkNL(myk, z0) );
}
//check fast
vector<double> pkfast(kval.size());
try{
e->get_Pkvec(kval,pkfast);
}
catch (Message_error &e){};
int a(0);
arr<double> apk(kval.size());
for (size_t i=0; i<kval.size() ; i++) apk[i]=kval[i];
fout.write_column(1,apk);
fout.write_column(++a,apk);
for (size_t i=0; i<kval.size() ; i++) apk[i]=pklin[i];
fout.write_column(2,apk);
fout.write_column(++a,apk);
if ( do_nl ){
for (size_t i=0; i<kval.size() ; i++) apk[i]=pknl[i];
fout.write_column(3,apk);
fout.write_column(++a,apk);
}
for (size_t i=0; i<kval.size() ; i++) apk[i]=pkfast[i];
fout.write_column(++a,apk);
fout.close();
}
delete e;
return 0;
}
......@@ -51,7 +51,7 @@ int main(int argc,char** argv){
pars.add("output","mPk"); //pol +clphi-> mPk added for sigma8
//need to compute P(k,z) for sigma8(z)
//pars.add("z_pk","0,0.57");
pars.add("z_pk","0,0.57");
pars.add("z_max_pk",1);
......@@ -66,7 +66,8 @@ int main(int argc,char** argv){
else{
engine=new ClassEngine(pars);
}
cout << "*****************************************" << endl;
//test keyword access
const char* derived_par[]={ "H0","H(0)","z_reio","tau_reio","z_drag","rs_drag","sigma8","Dv(0.57)","Da(0.57)","H(0.57)","sigma8(0.57)", "f(0.57)","growthD(0.57)"};
const int Nder=sizeof(derived_par)/sizeof(char*);
......@@ -99,7 +100,18 @@ int main(int argc,char** argv){
cout << "z=" << zi[i] << " D=" << engine->get_growthD(zi[i]) <<endl;
}
//test pk
struct spectra& sp=engine->sp;
cout << "ln_tau_size=" << sp.ln_tau_size << endl;
double kval[]={1e-4,3e-4,1e-3,2e-3,1e-2,3e-2,0.1};
const vector<double> kvec(kval,kval+sizeof(kval)/sizeof(double));
vector<double> pkvec(kvec.size());
//spectra_fast_pk_at_kvec(&ba,&sp,const_cast<double*>(&kvec[0]),kvec.size(),&pkvec[0]);
engine->get_Pkvec(kvec,pkvec);
for (size_t i=0;i<kvec.size();i++) cout << "k=" << kvec[i] << "\t p(k)=" << pkvec[i] << endl;
//clean
delete engine;
......
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