clean CLASS Pk aspects- add usefull exec

cmt/requirements-fastpk

@@ -94,7 +94,7 @@ macro application_suffix ""
 #groupe exec
 application writeChi2 -group=exec -s=$(CAMELROOT)/src/camel exec/writeChi2.cc
 application writeSpectra -group=exec -s=$(CAMELROOT)/src/camel exec/writeSpectra.cc
-application writeSpectraPk -group=exec -s=$(CAMELROOT)/src/camel exec/writeSpectraPk.cc
+application writePk -group=exec -s=$(CAMELROOT)/src/camel exec/writePk.cc
 application Minimize -group=exec -s=$(CAMELROOT)/src/camel exec/Minimize.cc
 application Profile -group=exec -s=$(CAMELROOT)/src/camel exec/Profile.cc
 application ScanParam -group=exec -s=$(CAMELROOT)/src/camel exec/ScanParam.cc

src/camel/Class/ClassEngine.cc

@@ -728,14 +728,16 @@ double  ClassEngine::get_DMod(double z){
 bool
 ClassEngine::get_Pkvec(const std::vector<double>& knodes, const std::vector<double>& zvec, std::vector<double>& pks)
 {
-  return ( _nonlin ? (spectra_fast_pk_nl_at_kvec_and_zvec(&ba,&sp,const_cast<double*>(&knodes[0]),knodes.size(),const_cast<double*>(&zvec[0]), zvec.size(),&pks[0])==_TRUE_) : (spectra_fast_pk_at_kvec_and_zvec(&ba,&sp,const_cast<double*>(&knodes[0]),knodes.size(),const_cast<double*>(&zvec[0]), zvec.size(),&pks[0])==_TRUE_) );
+  return ( _nonlin ? 
+(spectra_fast_pk_nl_at_kvec_and_zvec(&ba,&sp,const_cast<double*>(&knodes[0]),knodes.size(),const_cast<double*>(&zvec[0]), zvec.size(),&pks[0])==_TRUE_) 
+: (spectra_fast_pk_at_kvec_and_zvec(&ba,&sp,const_cast<double*>(&knodes[0]),knodes.size(),const_cast<double*>(&zvec[0]), zvec.size(),&pks[0])==_TRUE_) );
 
 }
 
 #endif
 
 bool
-ClassEngine::get_PkNodes(std::vector<double>& knodes,std::vector<double>& pknodes){
+ClassEngine::get_PklinNodes(std::vector<double>& knodes,std::vector<double>& pknodes){
 
   knodes.resize(sp.ln_k_size);
   pknodes.resize(sp.ln_k_size);
@@ -751,3 +753,21 @@ ClassEngine::get_PkNodes(std::vector<double>& knodes,std::vector<double>& pknode
   return _SUCCESS_;
 
 }
+
+bool
+ClassEngine::get_PkNLNodes(std::vector<double>& knodes,std::vector<double>& pknodes){
+
+  knodes.resize(sp.ln_k_size);
+  pknodes.resize(sp.ln_k_size);
+
+  int index_mode=0;
+  int index_ic1_ic2=0;
+  int index_eta = sp.ln_tau_size-1;
+  
+  for (int index_k=0; index_k < sp.ln_k_size; index_k++){
+    knodes[index_k]=std::exp(sp.ln_k[index_k]);
+    pknodes[index_k]=std::exp(sp.ln_pk_nl[(index_eta * sp.ln_k_size + index_k) * sp.ic_ic_size[index_mode] + index_ic1_ic2]);
+  }
+  return _SUCCESS_;
+
+}

src/camel/Class/ClassEngine.hh

@@ -110,11 +110,18 @@ public:
 
   //services
   virtual bool has_Background() const {return true;}
-  virtual bool has_CMB_Cl() const {return (_output.find("tCl")!=std::string::npos);}
-  virtual bool has_CMB_ClPol() const {return (_output.find("pCl")!=std::string::npos);}
-  virtual bool has_CMB_Lensing() const {return (_output.find("lCl")!=std::string::npos);}
-  virtual bool has_Pklin() const {return (_output.find("mPk")!=std::string::npos);}
-  virtual bool has_PkNL() const {return _nonlin;}
+//   //virtual bool has_CMB_Cl() const {return (_output.find("tCl")!=std::string::npos);}
+//   virtual bool has_CMB_ClPol() const {return (_output.find("pCl")!=std::string::npos);}
+//   virtual bool has_CMB_Lensing() const {return (_output.find("lCl")!=std::string::npos);}
+//   virtual bool has_Pklin() const {return (_output.find("mPk")!=std::string::npos);}
+//   virtual bool has_PkNL() const {return _nonlin;}
+
+  //class based
+  inline bool has_CMB_Cl() const {return (pt.has_cl_cmb_temperature==_TRUE_);}
+  inline bool has_CMB_ClPol() const {return (pt.has_cl_cmb_polarization==_TRUE_);}
+  inline bool has_CMB_Lensing() const {return (pt.has_cl_lensing_potential==_TRUE_);}
+  inline bool has_Pklin() const {return (pt.has_pk_matter==_TRUE_);}
+  inline bool has_PkNL() const {return (pt.has_nl_corrections_based_on_delta_m==_TRUE_);}
 
 
   //get value at l ( 2<l<lmax): in units = (micro-K)^2
@@ -157,26 +164,25 @@ public:
   //PK related stuuff
   double get_Pklin(double k, double z);
   double get_PkNL(double k, double z);
-  
-  //vectorized access (if macro defined)
-#ifdef FASTPK  
 
-  //fixed k and z values: output
-  bool get_Pkvec(const std::vector<double>& knodes, const std::vector<double>& z,std::vector<double>& pks);  
+  //nodes used by class to perfrom spline interpolation (at z=0)
+  //knodes/pknodes resized
+  bool get_PklinNodes(std::vector<double>& knodes,std::vector<double>& Pknodes);
+  bool get_PkNLNodes(std::vector<double>& knodes,std::vector<double>& Pknodes);
+  // the followingdepends on switch (prefer previous methods)
+  inline bool get_PkNodes(std::vector<double>& knodes,std::vector<double>& Pknodes) 
+  {return (has_PkNL()? get_PkNLNodes(knodes,Pknodes) : get_PklinNodes(knodes,Pknodes));} 
 
-   inline bool get_Pkvec(const std::vector<double>& kvec, 
-				const double& z,
-				std::vector<double>& pks){
+#ifdef FASTPK  
+  //vectorized access (if macro defined): with switch according to has_NL 
+  // outout vector pks must be properly sized (see Engine.hh)
+  bool get_Pkvec(const std::vector<double>& kin, const std::vector<double>& z,std::vector<double>& pks);  
+  inline bool get_Pkvec(const std::vector<double>& kvec,const double& z, std::vector<double>& pks){
     std::vector<double> zvec(1,z);
     return  this->get_Pkvec(kvec,zvec,pks);
   }
-
-
 #endif
 
-  //nodes used by class to perfrom spline interpolation: at z=0
-  //automatically resized
-  bool get_PkNodes(std::vector<double>& knodes,std::vector<double>& Pknodes);
 
 
   //combile distance

src/camel/Engine.hh

@@ -96,7 +96,7 @@ public:
   //compute Pk for fixed k and z values
   //output Pks should be sized as zvec.size()*kvec.size
   //default implementations calls P(k,z) function -> override for better performances
-  //most general form
+  //most general form:
   inline virtual bool get_Pkvec(const std::vector<double>& kvec, 
 			 const std::vector<double>& zvec,
 			 std::vector<double>& pks){
@@ -105,7 +105,8 @@ public:
 	pks[j*zvec.size()+i]=get_Pk(kvec[i],zvec[j]);
     return true;
   }
-	  
+
+  //defult is just a size 1 vector call
   inline virtual bool get_Pkvec(const std::vector<double>& kvec, 
 				const double& z,
 				std::vector<double>& pks){

src/camel/exec/writePk.cc

+#include "Parser.hh"
+#include "Class/MnClassEngine.hh"
+#include"cxxsupport/cxxutils.h"
+#include "Timer.hh"
+//#include"cxxsupport/arr.h"
+#include"cxxsupport/fitshandle.h"
+#include <iostream>
+#include<fstream>
+#include<cassert>
+#include<numeric>
+#include<string>
+#include<cmath>
+using namespace std;
+
+
+/*
+Utlities to check the CLASS Pk spectrum: fixed nodes location + inetrpolated values
+
+author : S. Plaszczynski, Aug.2017
+*/
+
+
+
+void logspace(double kmin,double kmax, int npks,vector<double>&v){
+
+  double lmax=std::log(kmax);
+  double lmin=std::log(kmin);
+  double step=(lmax-lmin)/(npks-1);
+
+  v.resize(npks);
+  for (int i=0;i<npks;i++){
+    double lval=lmin+step*i;
+    v[i]=std::exp(lval);
+  }
+}
+
+
+void linspace(double kmin,double kmax, int npks,vector<double>&v){
+
+  double step=(kmax-kmin)/(npks-1);
+  v.resize(npks);
+  for (int i=0;i<npks;i++){
+    double kval=kmin+step*i;
+    v[i]=kval;
+  }
+}
+
+
+
+
+
+int main(int argc,char** argv){
+
+  planck_assert(argc==3 || argc==4,"Usage: writePk parfile(in) file.fits(out) z_pk(in)");
+
+  assert_present(argv[1]);
+
+
+  float z0 = 0.;
+  if (argc==4) z0=atof(argv[3]);
+
+  //decodage arguments
+  Parser parser(argv[1]);
+
+  //choose engine
+  string engine= parser.params.find<string>("engine","CLASS");
+  if (engine=="pico") {
+    throw  Message_error("pico cannot be used for Pk !");
+  }
+
+  ClassParams classparms(parser.classparms);
+  classparms.add("output","mPk"); //polar +lens+clphi
+  
+  classparms.add("z_pk",z0); //polar +lens+clphi
+  
+
+  //no lensing
+  try{
+    size_t i=classparms.findIndex("lensing");
+    cout << "setting lensing to false" << endl;
+    classparms.updateVal(i,str(false));
+  }
+  catch (Message_error&){
+    classparms.add("lensing",false);
+  }
+  
+  //define kmax and kvalues either from file either given by k_end
+
+    double kmax=-1;
+    vector<double>  kval;
+    string kvfile=parser.params.find<string>("kval_file","");
+    bool h_rescale=parser.params.find<bool>("h_rescale",false);
+    //if file
+    if (kvfile.size()!=0) {
+      assert_present(kvfile);
+      cout <<"using k values from file=" <<kvfile <<endl;
+      //decode 1st column
+      ifstream f(kvfile.c_str());
+      string line;
+      double kk;
+      while (getline(f,line)) {
+	//skip comments
+	if (line.find("#",0) !=std:: string::npos) continue;
+	std::istringstream str(line);
+	if (str>>kk) {
+	  kval.push_back(kk);
+	}
+      }
+      kmax=kval.back();
+    }
+    else{
+      kmax = parser.params.find<double>("kmax",1);
+      const int npks = parser.params.find<int>("n_pks",1000);
+      const double kmin = parser.params.find<double>("kmin",8.e-4);
+      bool lin=parser.params.find<bool>("linspace",false);
+      lin? linspace(kmin,kmax,npks,kval) : logspace(kmin,kmax,npks,kval) ; 
+      
+    }
+    //modify class param accordingly
+    string label=(h_rescale? "P_k_max_h/Mpc" : "P_k_max_1/Mpc");
+    cout << label << " set to " << kmax <<endl;
+
+    try{
+      size_t i=classparms.findIndex(label);
+      classparms.updateVal(i,dataToString(kmax));
+    }
+    catch (Message_error&){
+      classparms.add(label,kmax);
+    }
+
+
+
+
+  //precision file
+  string pre="";
+  if (parser.params.param_present("precisionFile")){
+    pre=Parser::CheckPath(parser.params.find<string>("precisionFile",""));
+  }
+  
+  
+  cout << "CLASS \t--> precision parameters taken from file=" <<  pre << endl;
+  MnClassEngine* e=new MnClassEngine(parser.vars(),
+		      classparms,
+		      pre);
+ 
+
+  cout << "Computing P(k) at z=" << z0 << endl;
+  ClassEngine* klass=e->getClass();
+
+  //output
+  const std::string fileName(argv[2]);
+  cout << "writing FITS file =" << fileName << endl;
+  if (file_present(fileName)) {
+    //std::cout << "the file " + fileName + " already exists: removing"<<endl;
+    remove_file(fileName);
+  }
+  fitshandle fout(fileName,fitshandle::CREATE, READWRITE);
+
+  bool do_nl =e->has_PkNL();
+  
+  //nodes
+  {
+    std::vector<fitscolumn> pkcols;
+    pkcols.push_back(fitscolumn("knode", "1/Mpc",1,TDOUBLE));
+    pkcols.push_back(fitscolumn("pknode", "(Mpc)^3" ,1,TDOUBLE));
+    if (do_nl) pkcols.push_back(fitscolumn("pknodenl", "(Mpc)^3" ,1,TDOUBLE));
+    fout.insert_bintab(pkcols);
+    
+    vector<double> knode, pknode,pknodeNL;
+    klass->get_PkNodes(knode,pknode);
+     if (do_nl) klass->get_PkNLNodes(knode,pknodeNL);
+
+    for (size_t i=0;i<knode.size();i++){
+      fout.write_column(1,knode[i],i);
+      fout.write_column(2,pknode[i],i);
+       if (do_nl) fout.write_column(3,pknodeNL[i],i);
+    }
+
+
+
+  }
+
+  //interp values
+
+   //test fast Pk
+  /*
+  Timer timer;
+  cout <<"Pk timer test..." << endl;
+  vector<double> pkvec(kval.size());
+  for (size_t i=0;i<100;i++)
+    e->get_Pkvec(kval,z0,pkvec);
+  cout << "fast=" << timer.partial() <<endl;
+
+  for (size_t i=0;i<100;i++){
+    for (int k =0; k<kval.size() ; k++){
+      double pk=e->get_Pklin(kval[k], 0);
+    }
+  }
+  cout << "std=" << timer.partial() <<endl;
+  */
+
+
+  {
+    std::vector<fitscolumn> pkcols;
+    pkcols.push_back(fitscolumn("k", "1/Mpc",1,TDOUBLE));
+    pkcols.push_back(fitscolumn("pklin", "(Mpc)^3" ,1,TDOUBLE));
+    if (do_nl) pkcols.push_back(fitscolumn("pknl", "(Mpc)^3" ,1,TDOUBLE));
+
+    fout.insert_bintab(pkcols);
+    
+    //optional rescale
+    double h=e->get_H0()/100;
+    //warning h_rescaling
+    if (h_rescale)
+      std::transform(kval.begin(),kval.end(),kval.begin(),bind1st(multiplies<double>(),h));
+    cout <<"computing Pk from " << kval[0] << " to " << kval.back() << endl;
+
+    for (size_t i=0;i<kval.size();i++){
+      fout.write_column(1,kval[i],i);
+      fout.write_column(2,klass->get_Pklin(kval[i],z0),i);
+      if (do_nl) fout.write_column(3,klass->get_PkNL(kval[i],z0),i);
+    }
+    cout <<"Growth factor D(z=" << z0 << ")="<< klass->get_growthD(z0) << endl;
+    fout.add_key("D",klass->get_growthD(z0),"growth factor");
+    fout.add_key("f",klass->get_f(z0),"growth rate");
+
+  }
+  
+  fout.close();
+
+
+
+
+  delete e;
+
+  return 0;
+}