Merge branch 'dev' into 'master'

Dev

See merge request !15

cmt/requirements-lite.generic

@@ -4,7 +4,7 @@ use class *
 
 #compiler options
 #default is gcc here are the C++ compiler options
-macro cppflags " -O2 -pipe -fPIC -Wall "
+macro cppflags " -O2 -pipe -fPIC -Wall -Wno-reorder -Wno-sign-compare -Wno-unknown-pragmas"
 
 #do we need optimzations here?
 #macro_append cppflags "-ffast-math -m64 "

cmt/requirements-lsst.cc

@@ -5,7 +5,7 @@ use class HEAD
 #compiler options
 macro cpp "icpc"
 
-macro cppflags " -O2 -ip -ansi_alias -align -Wbrief -Wdeprecated -Wuninitialized -Wbrief -Wunused-function  -Wpointer-arith -Wreturn-type"
+macro cppflags " -O2 -ip -ansi_alias -align -Wall  -Wno-unknown-pragmas"
 #macro_append cppflags " -g "
 
 #CAMEL includes

src/AbsRand/OofNoise.cc

@@ -75,7 +75,9 @@ OofNoise::OofNoise(AbsGauss* r,double slope,double fmin,double fknee,double fsam
 
 OofNoise::~OofNoise(){
   for (unsigned i=0;i<ff.size();i++) delete ff[i];
-  if (_r) delete _r; _r=NULL;
+  if (_r) 
+      delete _r; 
+  _r=NULL;
 }
 
 double

src/camel/CMB/Lollipop.cc

@@ -86,7 +86,7 @@ void Lollipop::CreateDataSet( string DataFile)
   fp.close();
   
   if( (long)dataset.size() != Lmax-Lmin+1) {
-    printf( "Wrong number of l in datafile (%d)\n", (long)dataset.size());
+    printf( "Wrong number of l in datafile (%ld)\n", (long)dataset.size());
     exit(-1);
   }
 

src/camel/Class/ClassAdapters.hh

@@ -27,6 +27,7 @@ class ClassAdapter
 {
 public:
   virtual std::string className() {return name;};
+  virtual ~ClassAdapter(){};
   virtual std::string toClass(const std::vector<double>& par) const=0;
 
 protected:

src/camel/Class/ClassEngine.cc

@@ -94,7 +94,7 @@ ClassEngine::ClassEngine(const ClassParams& pars): cl(0),dofree(true),pvecback(0
   //prepare fp structure
   size_t n=pars.size();
   //
-  parser_init(&fc,n,"pipo",_errmsg);
+  parser_init(&fc,n,const_cast<char*>("pipo"),_errmsg);
   
   //config
   for (size_t i=0;i<pars.size();i++){
@@ -151,7 +151,7 @@ ClassEngine::ClassEngine(const ClassParams& pars,const string & precision_file):
   fc_input.filename=new char[1];
  //prepare fc par structure
   size_t n=pars.size();
-  parser_init(&fc_input,n,"pipo",_errmsg);
+  parser_init(&fc_input,n,const_cast<char*>("pipo"),_errmsg);
   //config
   for (size_t i=0;i<pars.size();i++){
     strcpy(fc_input.name[i],pars.key(i).c_str());
@@ -568,9 +568,9 @@ double ClassEngine::get_Pklin(double k, double z){
 #ifdef CLASSV_ABOVE_2_7
   double pk_cb=0.;
   double pk_cb_ic=0.;
-  int ret=spectra_pk_at_k_and_z(&ba,&pm,&sp,k,z,&mypk,pk_ic,&pk_cb,&pk_cb_ic);
+  assert(spectra_pk_at_k_and_z(&ba,&pm,&sp,k,z,&mypk,pk_ic,&pk_cb,&pk_cb_ic)==_SUCCESS_);
 #else
-  int ret=spectra_pk_at_k_and_z(&ba,&pm,&sp,k,z,&mypk,pk_ic);
+  assert(spectra_pk_at_k_and_z(&ba,&pm,&sp,k,z,&mypk,pk_ic)==_SUCCESS_);
 #endif
 
   return mypk;
@@ -587,9 +587,9 @@ double ClassEngine::get_PkNL(double k, double z){
 
 #ifdef CLASSV_ABOVE_2_7
   double pk_cb_tot=0;
-  int ret = spectra_pk_nl_at_k_and_z(&ba,&pm,&sp,k,z,&mypk,&pk_cb_tot);
+  assert(spectra_pk_nl_at_k_and_z(&ba,&pm,&sp,k,z,&mypk,&pk_cb_tot)==_SUCCESS_);
 #else
-  int ret = spectra_pk_nl_at_k_and_z(&ba,&pm,&sp,k,z,&mypk);
+  assert(spectra_pk_nl_at_k_and_z(&ba,&pm,&sp,k,z,&mypk)==_SUCCESS_);
 #endif
   return mypk;
 }
@@ -654,7 +654,6 @@ double ClassEngine::get_Dv(double z)
   //call to fill pvecback
   background_at_tau(&ba,tau,ba.long_info,ba.inter_normal, &index, pvecback);
 
-
   double H_z=pvecback[ba.index_bg_H];
   double D_ang=pvecback[ba.index_bg_ang_distance];
 #ifdef DBUG
@@ -663,8 +662,7 @@ double ClassEngine::get_Dv(double z)
 #endif
   double D_v;
 
-  D_v=pow(D_ang*(1+z),2)*z/H_z;
-  D_v=pow(D_v,1./3.);
+  D_v=std::pow(D_ang*(1+z)*D_ang*(1+z)*z/H_z,1./3.);
 #ifdef DBUG
   cout << D_v << endl;
 #endif
@@ -729,17 +727,9 @@ double ClassEngine::get_Da(double z)
   int index;
   //transform redshift in conformal time
   background_tau_of_z(&ba,z,&tau);
-
   //call to fill pvecback
   background_at_tau(&ba,tau,ba.long_info,ba.inter_normal, &index, pvecback);
-
-
-  double H_z=pvecback[ba.index_bg_H];
   double D_ang=pvecback[ba.index_bg_ang_distance];
-#ifdef DBUG
-  cout << "H_z= "<< H_z <<endl;
-  cout << "D_ang= "<< D_ang <<endl;
-#endif
   return D_ang;
 }
 

src/camel/Class/MnClassEngine.hh

@@ -38,7 +38,7 @@ public:
 		const std::string& prefile);
 
   // destructor
-  ~MnClassEngine();
+  virtual ~MnClassEngine();
 
   //this will update the CLASS computations for a realization of par 
   //conistent with the MnUserVariables

src/camel/Lensing/Chi2Lensing.hh

@@ -35,7 +35,7 @@ public:
   Chi2Lensing(LensingLikelihood* ll, const Variables& pars,Engine* _klass);
 
   // destructor
-  ~Chi2Lensing();
+  virtual ~Chi2Lensing();
 
   //operator to override: must return something comptible with -2log(L)
   double chi2_eff(const std::vector<double>& par) const;

src/camel/Lensing/LensingLikelihood.hh

@@ -22,7 +22,7 @@ class LensingLikelihood
 {
 
 public:
-
+  virtual ~LensingLikelihood(){};
   virtual int lmax() const =0;
  //convention: return -log(L)
   virtual double computePPLikelihood(const std::vector<unsigned int>& l,

src/camel/MCMC/MCMC_adaptive.cc

@@ -323,7 +323,8 @@ void MCMC_adaptive::fill_chain(pdf* target_func, pdf* proposal,int init_step, co
 	    int gate=_bunchSize;
 	    if(i%gate==0)
 	      {
-		if (_verbose) cout << "dump at index " << i << " from " << i-gate << " to " << i-1 << endl;
+              //if (_verbose) cout << "dump at index " << i << " from " << i-gate << " to " << i-1 << endl;
+              cout << i << " written" <<endl;
 		chain_running.open (chain_run_name, fstream::app);
 		
 		for(int k=i-gate; k<i; k++)

src/camel/MCMC/MCMC_sampler.cc

@@ -136,7 +136,7 @@ void MCMC_sampler::write_chain(char* str)
 
 }
 
-void MCMC_sampler::write_if_useful(vector<double>& vec,char* str)
+void MCMC_sampler::write_if_useful(vector<double>& vec,const char* str)
 {
 
   cout << "output written in:"<< str << endl;
@@ -149,7 +149,7 @@ void MCMC_sampler::write_if_useful(vector<double>& vec,char* str)
 
 } 
 
-void MCMC_sampler::write_vecHep(vector<HepVector >& vec,char* str)
+void MCMC_sampler::write_vecHep(vector<HepVector >& vec,const char* str)
 {
 
   cout << "output written in:"<< str << endl;

src/camel/MCMC/MCMC_sampler.hh

@@ -67,8 +67,8 @@ public:
   void compute_mean(int stop,HepVector mean);
   void write_rearranged_chain(char* str); //rec, value, lkh
   void set_par(Parser parser); // set parameters value using Parser
-  void write_if_useful(std::vector<double>& vec,char* str);
-  void write_vecHep(std::vector<HepVector >& vec,char* str);
+  void write_if_useful(std::vector<double>& vec,const char* str);
+  void write_vecHep(std::vector<HepVector >& vec,const char* str);
 
   
   //accessors

src/camel/MCMC/exec/mainMCMC.cc

@@ -19,6 +19,8 @@
 #include "Chi2Factory.hh"
 
 #include "cxxsupport/paramfile.h"
+#include "cxxsupport/cxxutils.h"
+
 #include "Chi2Combiner.hh"
 #include "Parser.hh"
 
@@ -48,7 +50,17 @@ int main(int argc, char *argv[])
     std::cout << "output file " + fileName + " already exists: removing it"<<endl;
     remove_file(fileName);
   }
+  string base=BaseName(fileName);
+
+  string ar_file=string("ar_vs_length_")+base;
+  string scale_file=string("scale_")+base;
+  string cor_file=string("last_step_")+base;
 
+  cout << "MCMC output files:" <<endl;
+  cout << ar_file << " " << scale_file << " " << cor_file << " " << endl;
+
+  cout << "output files:" <<endl;
+  cout << ar_file << " " << scale_file << " " << cor_file << " " << endl;
 
   //MCMC specific parameters
   const string algo=parser.params.find<string>("algo","ada");
@@ -77,7 +89,6 @@ int main(int argc, char *argv[])
   const int dim=lkh->Get_par_num();
 
   cout <<"dim=" <<dim <<endl;
-
   //--------defining seeds-----------------------------------------------
   unsigned int seed(0);
   if (parser.params.param_present("seed")){
@@ -151,8 +162,7 @@ int main(int argc, char *argv[])
 
   if(use_ext_cov)
   {
-    cout << "Read proposal cov matrix.." << endl;
-    cout << proposal_cov.c_str() << endl;
+    cout << "Read proposal cov matrix=" << proposal_cov.c_str() << endl;
     proposal->set_cov(proposal_cov.c_str());
   }
   else {
@@ -245,21 +255,21 @@ int main(int argc, char *argv[])
          {ar_vs_length.push_back(chain->acc_rate(0,step));}
      }
 
-   chain->write_if_useful(ar_vs_length,"ar_vs_length.txt");
+   chain->write_if_useful(ar_vs_length,ar_file.c_str());
 
    if(algo=="ada")
    {
      MCMC_adaptive* ada=dynamic_cast<MCMC_adaptive*>(chain);
-     chain->write_if_useful(proposal->_scale_vs_length,"scale_vs_length.txt");
+     chain->write_if_useful(proposal->_scale_vs_length,scale_file.c_str());
 
-     ofstream f("corr.txt");
+     ofstream f(cor_file.c_str());
      if(!f) {
        cout << "Error in file creation!";
      }
      f << ada->mean() << endl //mean
        << proposal->Get_cov_matrix() << endl; //cov matrix
      f.close();
-     cout<<"output written in:corr.txt"<<endl;;
+     cout<<"proposal cov written in " << cor_file <<endl;;
    }
 
    delete chi2;

src/cxxsupport/cxxutils.cc

@@ -294,3 +294,15 @@ string DirName(string source)
     source.erase(std::find(source.rbegin(), source.rend(), '/').base(), source.end());
     return source;
 }
+
+string BaseName(string source)
+{
+    if (source.size() <= 1) //Make sure it's possible to check the last character.
+    {
+        return source;
+    }
+    planck_assert(*(source.rbegin() + 1) != '/',source +" had no BaseName");
+
+    source.erase(source.begin(),std::find(source.rbegin(), source.rend(), '/').base());
+    return source;
+}

src/cxxsupport/cxxutils.h

@@ -264,6 +264,8 @@ inline unsigned int healpix_repcount (int npix)
 
   //acess directory name without trailing slash
   std::string DirName (std::string x);
+  std::string BaseName (std::string x);
+
 
 
 

work/batch/cc/Minimize5.sh

@@ -49,8 +49,7 @@ fi
 
 
 parbase=$(basename $file)
-
-dirout=${parbase%".par"}_"min"
+dirout=$(echo $parbase | awk -F"." '{NF--;printf($0)}')_"min"
 OUTDIR=$PWD/$dirout
 
 
@@ -115,7 +114,7 @@ RANDOM=\${SGE_TASK_ID}
 awk -v seed=\$RANDOM -f genrand.awk parfile_in > parfile
 fi
 
-cp parfile $OUTDIR/camel\${SGE_TASK_ID}.par
+cp parfile $OUTDIR/camel\${SGE_TASK_ID}.cml
 echo "***************************************************************"
   #================================
   # Camel run

work/batch/cc/runMCMC_class.sh

@@ -26,7 +26,8 @@ fi
 
 parFile=$(readlink -e $file)
 parbase=$(basename $parFile)
-dirout=${parbase%".par"}_"MC"
+#rm extension (whatever it is)
+dirout=$(echo $parbase | awk -F"." '{NF--;printf($0)}')_"MC"
 zeroot=$PWD/$dirout
 
 
@@ -59,7 +60,8 @@ if [ $ndim -ne $nlines ] ; then
 echo "covariance matrix not compatible in size with parFile"
 exit 1
 fi
-scale=$(echo "$ndim" | awk '{printf(2.38*2.38/$0)}')
+# now defult behavior of teh application
+#scale=$(echo "$ndim" | awk '{printf(2.38*2.38/$0)}')
 
 
 cat > mc_setup <<EOF
@@ -67,7 +69,6 @@ length=$NSAMPLES
 proposal_cov=$cov
 derived=sigma8 
 do_mPK=true
-scale=$scale
 EOF
 
 echo "**************************************"
@@ -111,9 +112,7 @@ cp $parbase $PWD/camel.par
 
 ./$localexec $parbase $PWD/samples\${SGE_TASK_ID}.txt
 
-cp ar_vs_length.txt $PWD/ar_vs_length\${SGE_TASK_ID}.txt
-cp scale_vs_length.txt $PWD/scale_vs_length\${SGE_TASK_ID}.txt
-cp corr.txt $PWD/corr\${SGE_TASK_ID}.txt
+cp *.txt $PWD/
 
 qstat -j \${JOB_ID} -nenv
 eof

work/ref/hlpTT_PS_H0.par → work/ref/hlpTT_PS_H0.cml

@@ -71,10 +71,8 @@ class     sBBN\ file            bbn/sBBN_2017_marcucci.dat
 
 precisionFile=class_pre/hpjul3.pre
 ###############################################################
+# minimization
 
-
-###############################################################
-#fitter
 nitermax=50000
 set_stra=2
 set_tol=0.00001
@@ -83,3 +81,14 @@ set_tol=0.00001
 remove_cosmo_limits=false
 doHesse=true
 ###############################################################
+#MCMC
+#1st samples are drawn up to t0 from given proposal (scaled by scale)
+#if your matrix is nice scale_th~2.4**2/ndim otherwise low (default=0.001)
+#t0=2000
+#scale=0.3
+
+#then cov is adapted up to ts with scale_th. If your matrix is nice do
+# not destroy things there.
+#ts=0
+
+#then scale will auto-adapt to reach ar~0.25