met au propre la pertie lmax- clarifie les services des Engine- change get_pk

cmt/requirements-gcc

@@ -4,7 +4,7 @@ use class HEAD
 
 #compiler options
 #default is gcc here are the C++ compiler options
-macro cppflags " -O3  -pipe -ansi -Woverloaded-virtual -Wpedantic -Wunused -Wmaybe-uninitialized -Wno-deprecated -Wno-div-by-zero "
+macro cppflags " -g  -pipe -ansi -Woverloaded-virtual -Wpedantic -Wunused -Wmaybe-uninitialized -Wno-deprecated -Wno-div-by-zero "
 
 # although it may be strange we do not compile CAMEL with OMP because CLASS is not thread safe
 # no worries: CLASS is compiled with OMP

src/camel/Chi2Factory.cc

@@ -159,9 +159,10 @@ Chi2Factory::gimeChi2(Parser& parser){
     allCl->add(acte);
   }
 
-  cout << "NUMBER OF LIKELIHOODS=" << allCl->nLik() << endl;
-  cout << "TTmax=" << allCl->getTTmax() << endl;
+  cout << "NUMBER OF CMB LIKELIHOODS=" << allCl->nLik() << endl;
   int lmax=allCl->getTTmax();
+  cout << "TTmax=" << lmax << endl;
+
 
   //lensing may affect lmax too
 #ifdef CLIK
@@ -175,18 +176,10 @@ Chi2Factory::gimeChi2(Parser& parser){
     if (!fid && lensing->lmax() > lmax) lmax=lensing->lmax();
   }
 #endif
-
-//overide if you want
-if (parser.params.param_present("lmax")){
-  int Lmax=parser.params.find<int>("lmax");
-  if (lmax != Lmax) 
-    cout << " WARNING your specified LMAX=" << Lmax <<" differs from the guessed TT range:" << lmax <<" use at you own risks!" << endl;
-lmax=Lmax;
-}
+  
 
   cout << "#########################################" << endl;
 
-
   cout << "Engine lmax=" << lmax << endl;
 
   //2. creation ENGINE
@@ -207,17 +200,24 @@ lmax=Lmax;
 //extra class params
   ClassParams classparms(parser.classparms);
 
+
   //CMB
   string output="";
   if (lmax>0) {
     classparms.add("l_max_scalars",lmax);
-    output+="tCl,pCl,lCl";
+    output="tCl,pCl,lCl";
     
   }
   //if P(k) needed
   const bool doPk=parser.params.find<bool>("do_mPk",false);
-  if (doPk) output+=",mPk";
 
+  //if pk asked or no CMB
+  if (doPk || lmax<2) {
+    output+=(output.size()==0? "mPk" : ",mPk");
+  }
+
+  cout << ">>> CLASS output will be: " << output << endl;
+ 
   classparms.add("output",output); //polar +lens+clphi
 
   for (size_t i=0;i<classparms.size();i++) 
@@ -230,7 +230,6 @@ lmax=Lmax;
   }
   cout << "CLASS \t--> precision =" <<  pre << endl;
 
-
   //le calculateur de spectres
   //MnUserVariables upar(parser.upar);
   e=new MnClassEngine(parser.vars(),
@@ -244,9 +243,10 @@ lmax=Lmax;
   //4. observables:
 
   //rajoute l'engine aux Cl
-  allCl->setEngine(e);
-
-  comb->add(allCl);
+  if (lmax>=2 ){
+    allCl->setEngine(e);
+    comb->add(allCl);
+  }
 
   //lensing
 #ifdef CLIK 

src/camel/Class/ClassEngine.cc

@@ -85,7 +85,7 @@ void ClassParams::updateVal(const unsigned& i,const string& newval) {pars[i].sec
 //---------------
 // Constructors --
 //----------------
-ClassEngine::ClassEngine(const ClassParams& pars): cl(0),dofree(true),_hasPk(false){
+ClassEngine::ClassEngine(const ClassParams& pars): cl(0),dofree(true),_nonlin(false){
 
   _lmax=-1; //default
 
@@ -108,13 +108,13 @@ ClassEngine::ClassEngine(const ClassParams& pars): cl(0),dofree(true),_hasPk(fal
       istringstream strstrm(pars.value(i));
       strstrm >> _lmax;
     }
-    if (pars.key(i)=="output") {
-	if (pars.value(i).find("mPk") != std::string::npos) _hasPk=true;
-    }
-  }
-  cout << __FILE__ << " : using lmax=" << _lmax <<endl;
-  assert(_lmax>0);
+    //store class output
+    if (pars.key(i)=="output") _output=pars.value(i);
 
+    //is nonlinear set
+    if (pars.key(i)=="non linear") _nonlin=true;
+
+  }
   //input
   if (input_init(&fc,&pr,&ba,&th,&pt,&tr,&pm,&sp,&nl,&le,&op,_errmsg) == _FAILURE_) 
     throw invalid_argument(_errmsg);
@@ -134,7 +134,7 @@ ClassEngine::ClassEngine(const ClassParams& pars): cl(0),dofree(true),_hasPk(fal
 
   }
 
-ClassEngine::ClassEngine(const ClassParams& pars,const string & precision_file): cl(0),dofree(true),_hasPk(false){
+ClassEngine::ClassEngine(const ClassParams& pars,const string & precision_file): cl(0),dofree(true),_nonlin(false){
 
   cout << "Running CLASS version " << _VERSION_ << endl;
 
@@ -162,17 +162,12 @@ ClassEngine::ClassEngine(const ClassParams& pars,const string & precision_file):
       istringstream strstrm(pars.value(i));
       strstrm >> _lmax;
     }
-    if (pars.key(i)=="output") {
-      if (pars.value(i).find("mPk") != std::string::npos) {
-	_hasPk=true;
-	cout << "found Pk: status" << (int)hasPk() << endl;
-      }
-    }
-  }
-  cout << __FILE__ << " : using lmax=" << _lmax <<endl;
-  assert(_lmax>0);
-  
+    //store class output
+    if (pars.key(i)=="output") _output=pars.value(i);
+    //is nonlinear set?
+    if (pars.key(i)=="non linear") _nonlin=true;
 
+  }
 
   //concatenate both
   if (parser_cat(&fc_input,&fc_precision,&fc,_errmsg) == _FAILURE_) throw invalid_argument(_errmsg);
@@ -550,7 +545,7 @@ double ClassEngine::get_f(double z)
   return f_z;
 }
 
-double ClassEngine::get_Pk(double k, double z){
+double ClassEngine::get_Pklin(double k, double z){
   double tau;
   int index;
   double *pvecback;

src/camel/Class/ClassEngine.hh

@@ -108,6 +108,15 @@ public:
   //update from ClassParams class (should be same size than in constructor)
   bool updateParValues(const ClassParams& par);
 
+  //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;}
+
+
   //get value at l ( 2<l<lmax): in units = (micro-K)^2
   //don't call if FAILURE returned previously
   //throws std::execption if pb
@@ -126,7 +135,6 @@ public:
 	      std::vector<double>& clephi);
 
 
-  bool hasPk() const {return _hasPk;}
   //recombination
   inline double z_rec() const {return th.z_rec;}
   inline double rs_rec() const {return th.rs_rec;} 
@@ -143,7 +151,8 @@ public:
   double get_Fz(double z);
   double get_H(double z);
   double get_DMod(double z);
-  double get_Pk(double k, double z);
+
+  double get_Pklin(double k, double z);
   double get_PkNL(double k, double z);
 
   double com_distance(double z);
@@ -215,8 +224,8 @@ private:
 
 protected:
  
-  bool _hasPk;
-  
+  std::string _output;
+  bool _nonlin;
 };
 
 #endif

src/camel/Class/MnClassEngine.hh

@@ -58,9 +58,14 @@ public:
 		  std::vector<double>& clep) {return klass->getLensing(lVec,clpp,cltp,clep);}
 
   //forward
-  bool hasPk() const {return klass->hasPk();}
-  double get_Pk(double k, double z)  { return klass->get_Pk(k,z);}
-  double get_PkNL(double k, double z)  { return klass->get_PkNL(k,z);}
+  //services
+  virtual bool has_Background() const {return klass-> has_Background();}
+  virtual bool has_CMB_Cl() const {return klass->has_CMB_Cl();}
+  virtual bool has_CMB_ClPol() const {return klass->has_CMB_ClPol();}
+  virtual bool has_CMB_Lensing() const {return klass->has_CMB_Lensing();}
+  virtual bool has_Pklin() const {return klass->has_Pklin();}
+  virtual bool has_PkNL() const {return klass->has_PkNL();}
+
 
 
   inline double z_drag() const {return klass->z_drag();}

src/camel/Engine.hh

@@ -32,6 +32,17 @@ public:
   //pure virtual:
   virtual bool updateParValues(const std::vector<double>& cosmopars)=0;
 
+  //which services are provides
+  virtual bool has_Background() const {return false;}
+
+  virtual bool has_CMB_Cl() const {return false;}
+  virtual bool has_CMB_ClPol() const {return false;}
+  virtual bool has_CMB_Lensing() const {return false;}
+
+  virtual bool has_Pklin() const {return false;}
+  virtual bool has_PkNL() const {return false;}
+
+  
   // units = (micro-K)^2
   virtual void getCls(const std::vector<unsigned>& lVec, //input
 		      std::vector<double>& cltt,
@@ -45,8 +56,6 @@ public:
 			  std::vector<double>& cltp,
 			  std::vector<double>& clep)=0;
 
-  virtual bool hasPk() const {return false;}
-
   //"derived" parameters
   //if not overriden default is to throw an exception
   //be aware for CLASS you may need to turn P(k) computation on ("do_mPk=true" in CAMEL) for sigma8.
@@ -65,14 +74,21 @@ public:
   virtual double get_DMod(double z) { return undef("Dmod(z)");}
   virtual double Omega_m() const { return undef("Omega_m");}
   virtual double YHe() const { return undef("YHe");}
+  virtual double z_rec() const {return undef("zrec");}
+  virtual double rs_rec() const {return undef("rs(rec)");} 
+
   //distance (comobile Mpc)
   virtual double com_distance(double z) { return undef("com_distance");}
 
-  virtual double get_Pk(double k, double z)  { return undef("getpk");}
-  virtual double get_PkNL(double k, double z)  { return undef("getpknl");}
 
-  virtual double z_rec() const {return undef("zrec");}
-  virtual double rs_rec() const {return undef("rs(rec)");} 
+  //P(K): convention if NL defined used it, oterwise lin
+  virtual double get_Pk(double k,double z=0) 
+  {
+    return (has_PkNL() ? get_PkNL(k,z) : get_Pklin(k,z));
+  }
+  //rathe prefer explicit calls
+  virtual double get_Pklin(double k, double z)  { return undef("getpklin");}
+  virtual double get_PkNL(double k, double z)  { return undef("getpknl");}
 
   //access thtought keyword: convention is that kw is name of previous methods without "get_"
   // z values can be specified by appending them in parenthesis, eg. "H(0.57)" for get_H(0.57).
@@ -89,7 +105,7 @@ public:
 protected:
   int _lmax;
 
-  virtual double undef(const char* c) const;
+  virtual double undef(const char* msg) const;
 
 };
 

src/camel/exec/writeChi2.cc

@@ -52,7 +52,7 @@ int main(int argc,char** argv){
     cout <<  extra[i] << "=" << val << endl;
   }
   //sigma8 can only be compute if Pk available
-  if (theChi2->getEngine()->hasPk()){
+  if (theChi2->getEngine()->has_Pklin()){
     cout << "sigma8" << "=" << theChi2->getEngine()->get("sigma8") << endl;
   } else
     {

src/camel/exec/writeSpectraPk.cc

@@ -58,16 +58,7 @@ int main(int argc,char** argv){
     classparms.add("z_pk",z0); //polar +lens+clphi
     classparms.add("z_max_pk",z0); //polar +lens+clphi
 
-    bool do_nl =false;
-    try{
-      size_t test =  classparms.findIndex("non linear");
-      cout << " >>> CLASS nl par found !! "<< test<<endl;
-      do_nl = true;
-    }catch (Message_error& msg){
-      
-      cout << " >>> CLASS nl par NOT found !! " <<endl;      
-      
-    }
+  
 
     for (size_t i=0;i<classparms.size();i++) 
       cout << "CLASS \t--> "<< classparms.key(i) << "\t" << classparms.value(i) <<endl;
@@ -83,6 +74,8 @@ int main(int argc,char** argv){
 			classparms,
 			pre);
 
+    bool do_nl =e->has_PkNL();
+
     //derived
     //expansion
     double H0=e->get_H0();
@@ -230,7 +223,7 @@ int main(int argc,char** argv){
   vector<double>  pknl ;
   for (int k =0; k<kval.size() ; k++){
     double myk=kval[k];
-    pklin.push_back( e->get_Pk(myk,  z0)   );
+    pklin.push_back( e->get_Pklin(myk,  z0)   );
     if ( do_nl ) pknl.push_back(  e->get_PkNL(myk,  z0) );
 
   }

src/camel/pico/MnPicoEngine.hh

@@ -40,6 +40,11 @@ public:
   // destructor
   ~MnPicoEngine();
 
+  //services
+  bool has_CMB_Cl() const {return true;}
+  bool has_CMB_ClPol() const {return true;}
+  bool has_CMB_Lensing() const {return true;}
+
   //this will update the computations for a realization of par 
   //conistent with the MnUserVariables
   bool updateParValues(const std::vector<double>& minuitPars);

src/camel/test/testKlass.cc

@@ -44,9 +44,11 @@ int main(int argc,char** argv){
   pars.add("k_pivot",0.05);
   pars.add("sBBN file",Parser::CheckPath("bbn/sBBN.dat"));
 
-  pars.add("l_max_scalars",l_max_scalars);
-  pars.add("lensing",true); //note boolean
-  pars.add("output","tCl,pCl,lCl,mPk"); //pol +clphi-> mPk added for sigma8
+  //pars.add("l_max_scalars",l_max_scalars);
+  //pars.add("lensing",true); //note boolean
+
+  //pars.add("output","tCl,pCl,lCl,mPk"); //pol +clphi-> mPk added for sigma8
+  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");