mainMCMC.cc 7.9 KB
Newer Older
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
#include "MCMC/MCMC_sampler.hh"
#include "MCMC/MCMC_adaptive.hh"
#include "MCMC/MCMC_metropolis.hh"
#include "MCMC/pdf_gaussian.hh"
#include "MCMC/pdf_log_gauss.hh"
#include "MCMC/pdf_chi_square.hh"
#include "MCMC/pdf.hh"
#include "MCMC/prior.hh"
#include "MCMC/posterior.hh"
#include "MCMC/post_planck.hh"
#include "MCMC/prior_flat.hh"
#include "MCMC/prior_beta.hh"
#include "MCMC/planck_lkh.hh"
#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include "CLHEP/Matrix/Matrix.h"
#include <ctime>
Xavier Garrido's avatar
Xavier Garrido committed
18
//planck lkh
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
19 20 21 22 23 24 25 26 27
#include "Chi2Factory.hh"

#include "cxxsupport/paramfile.h"
#include "Chi2Combiner.hh"
#include "Parser.hh"

#include <iostream>
#include <fstream>
#include <vector>
28

Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
29
//#include <ctime>
Xavier Garrido's avatar
Xavier Garrido committed
30
#include <sys/time.h>
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
31
#include <cstdlib>
32

Xavier Garrido's avatar
Xavier Garrido committed
33
using namespace std;
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
34 35 36

int main(int argc, char *argv[])
{
Xavier Garrido's avatar
Xavier Garrido committed
37

Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
38
  planck_assert(argc==3,"Usage: mcmc parfile(in) txtfile(out)");
Matthieu Tristram's avatar
Matthieu Tristram committed
39 40

  planck_assert(file_present(argv[1]),"missing par file");
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
41

Xavier Garrido's avatar
Xavier Garrido committed
42 43

  Parser parser(argv[1]);
44
  const char* chain_run_name=argv[2];
45

46 47 48 49 50
  const string fileName(chain_run_name);
  if (file_present(fileName)) {
    std::cout << "output file " + fileName + " already exists: removing it"<<endl;
    remove_file(fileName);
  }
51 52
  string ar_file=string("ar_vs_length_")+fileName;
  string scale_file=string("scale_")+fileName;
53
  string cor_file=string("last_step_")+fileName;
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
54

55 56
  cout << "output files:" <<endl;
  cout << ar_file << " " << scale_file << " " << cor_file << " " << endl;
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
57

58 59
  //MCMC specific parameters
  const string algo=parser.params.find<string>("algo","ada");
60
  const int length=parser.params.find<int>("length",1000);
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
61 62 63
  const string target_cor_path=parser.params.find<string>("target_cov","");
  const string proposal_cor=parser.params.find<string>("proposal_cor","");
  const string proposal_cov=parser.params.find<string>("proposal_cov","");
64
  const int t0=parser.params.find<int>("t0",2000);
65
  const int ts=parser.params.find<int>("ts",10000);
66
  double scale=parser.params.find<double>("scale",-1);
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
67 68 69
  const bool do_move=parser.params.find<bool>("do_move",false);
  bool verbose=parser.params.find<bool>("verbose",false);

70 71
  const unsigned long bunchSize=parser.params.find<unsigned long>("bunchSize",1000);

Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
72

73
  bool use_ext_cor=(proposal_cor.size()!=0);
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
74
  bool use_ext_cov=(proposal_cov.size()!=0);
Xavier Garrido's avatar
Xavier Garrido committed
75

Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
76
  //---------------------------------------------------------
Xavier Garrido's avatar
Xavier Garrido committed
77

Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
78 79 80 81
  Chi2Combiner* chi2=Chi2Factory::gimeChi2(parser);
  chi2->setVerbose(false);

  planck_lkh* lkh=new planck_lkh(chi2,parser);
Xavier Garrido's avatar
Xavier Garrido committed
82
  const int dim=lkh->Get_par_num();
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
83

84
  cout <<"dim=" <<dim <<endl;
Xavier Garrido's avatar
Xavier Garrido committed
85

86 87 88 89 90
  //SP 25/06/19: compute scale if undefined depending on dim
  if (scale<0) {
      scale=2.38*2.38/dim;
      cout << "refining scale to " << scale << endl;
  }
91 92 93
  //--------defining seeds-----------------------------------------------
  unsigned int seed(0);
  if (parser.params.param_present("seed")){
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
94
    seed=parser.params.find<int>("seed",0);
95
  }
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
96 97 98 99
  else {//pick automatiocally: microseconds from current time
    timeval CurrentTime;
    gettimeofday( &CurrentTime, NULL );
    seed=static_cast<unsigned int>(CurrentTime.tv_usec);
100 101 102 103
  }
  cout << "seed for proposal=" << seed << endl;

   unsigned int seedU(0);
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
104
  //need another one for the algorithm itsel
105
  if (parser.params.param_present("seedU")){
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
106
    seedU=parser.params.find<int>("seedU",0);
107
  }
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
108 109 110
  else {//throw a random number to avoid correlations
    std::srand(seed);
    seedU=static_cast<unsigned int>(std::rand()%32768);
111 112 113
  }
  cout << "seed for algorithm=" << seedU << endl;

Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
114 115 116 117 118 119
  //------------------------------------------------------------------
  //choosing the algorithm to use
  //for the moment just metropolis or adaptive.Metropolis is the default one
  MCMC_sampler* chain;
   if(algo=="ada")
   {
120
     chain= new MCMC_adaptive(dim,length,seedU,verbose);
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
121 122 123
     cout << "using Adaptive algorithm" << endl;
   }
   else
Xavier Garrido's avatar
Xavier Garrido committed
124
   {
125
     chain= new MCMC_metropolis(dim,length,seedU,verbose);
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
126 127
     cout << "using Metropolis algorithm" << endl;
   }
Xavier Garrido's avatar
Xavier Garrido committed
128

129
   chain->SetBunchSize(bunchSize);
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
130 131

   //directly from Parser
Xavier Garrido's avatar
Xavier Garrido committed
132 133
   chain->set_par(parser);

134 135 136 137 138
   if (verbose){
     cout <<  chain->min_par() << endl;
     cout <<  chain->max_par() << endl;
     cout << chain->starting_point() << endl;
   }
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
139

Xavier Garrido's avatar
Xavier Garrido committed
140 141
   //or by hand using set/get functions

142 143
   vector<double> par_err(dim);
   for(int i=0;i<dim;i++)
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
144
    {
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
145
      par_err[i]=(chain->par_err())[i];
146
      if (verbose) cout << par_err[i] <<endl;
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
147
    }
Xavier Garrido's avatar
Xavier Garrido committed
148

Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
149 150 151 152
  //chose the proposal-------------------------------------------------
  HepVector  mean(dim,0);
  HepVector  mean_pro(dim,0);
  //-----------------------------------------
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
153
  for(int i=0;i<dim;i++) {mean[i]=(chain->starting_point())[i]; }
Xavier Garrido's avatar
Xavier Garrido committed
154

Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
155 156 157
  HepSymMatrix cov(dim,1);
  HepSymMatrix identity(dim,1);

Xavier Garrido's avatar
Xavier Garrido committed
158

Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
159 160 161
  //-----define the lkh----------------------------------------------
  pdf_gaussian* proposal= new pdf_gaussian(dim,mean_pro,identity);

162
  if(use_ext_cov)
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
163
  {
164
    cout << "Read proposal cov matrix=" << proposal_cov.c_str() << endl;
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
165 166 167
    proposal->set_cov(proposal_cov.c_str());
  }
  else {
168
    if(use_ext_cor)
Xavier Garrido's avatar
Xavier Garrido committed
169
    {
Matthieu Tristram's avatar
Matthieu Tristram committed
170 171
      cout << "Read proposal cor matrix.." << endl;
      cout << proposal_cor.c_str() << endl;
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
172
      proposal->set_cov(proposal_cor.c_str());
173
      use_ext_cov= true;
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
174 175 176 177
    }
    cov = proposal->Get_cov_matrix();
    proposal->calc_cov_matrix( cov, par_err);
  }
Xavier Garrido's avatar
Xavier Garrido committed
178

179
  if (verbose) cout << "cov proposal= " << proposal->Get_cov_matrix() << endl;
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
180 181 182

  //scale
  proposal->Set_scale(scale);
183
  if (verbose) cout <<"scale=" << proposal->Get_scale() << endl;
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
184 185

  vector<double> scale_vs_length;
186
  proposal->init_generator(seed);
Xavier Garrido's avatar
Xavier Garrido committed
187

Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
188 189 190 191 192 193 194


  //if we want to displace the starting point
   if(do_move==true)
     {
       pdf_gaussian* move= new pdf_gaussian(dim,mean_pro,identity);
       move->calc_cov_matrix(identity, par_err);
195
       if (verbose) cout <<  "cov domove= " <<move->Get_cov_matrix() << endl;
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
196 197 198 199
       if(use_ext_cov==true){ move->Set_cov_matrix(proposal->Get_cov_matrix());}
       HepVector start(dim);
       move->init_generator(time(NULL));
       move->Set_scale(scale);
200
       if (verbose) cout << "start=" << start << endl;
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
201
       chain->Set_starting_point(chain->starting_point()+start);
202
       if (verbose) cout << "start_new=" << chain->starting_point() << endl;
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
203
     }
Xavier Garrido's avatar
Xavier Garrido committed
204

Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
205 206 207 208 209 210 211

   chain->Set_burn_in(0);


   vector<pdf*> prior_list(dim);
   for(int i=0; i<dim; i++)
     {
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
212
       pdf* p=new prior_flat((chain->min_par())[i],(chain->max_par())[i]);
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
213 214 215
       //pdf* p=new prior_beta(chain->_min_par[i],chain->_max_par[i], 1,1);
       prior_list[i]=p;
     }
Xavier Garrido's avatar
Xavier Garrido committed
216 217


Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
218 219 220 221
   prior* flat_p=new prior(prior_list);
   pdf* target=new post_planck(lkh,flat_p);


Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
222 223 224 225 226 227
   //set t0/ts for adaptive
   if (algo=="ada") {
     MCMC_adaptive* ada=dynamic_cast<MCMC_adaptive*>(chain);
     ada->Set_t0(t0);
     ada->Set_ts(ts);
   }
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
228 229 230 231 232 233

   //int init_step=atoi(argv[???]); //step to restart from to continue the chain
   chain->fill_chain(target,proposal,0, chain_run_name);

   //proposal->write_cov("cov_proposal.txt");

Xavier Garrido's avatar
Xavier Garrido committed
234

Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
235 236
   //acceptance rate simple
   double rate=chain->acc_rate(0,length);
237
   if (verbose) cout << "acc_rate=" << rate << endl ;
Xavier Garrido's avatar
Xavier Garrido committed
238

Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
239 240 241

   //---autocorr----and acc rate-------------
   vector<double > ar_vs_length;
Xavier Garrido's avatar
Xavier Garrido committed
242 243
   // int max_lag;
   // int min_lag;
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
244 245 246
   vector<HepVector > autocorr_vs_length;
   for(int step=100; step<=length; step+=100)
     {
Xavier Garrido's avatar
Xavier Garrido committed
247 248
       // max_lag=50;
       // min_lag=max_lag;
Matthieu Tristram's avatar
Matthieu Tristram committed
249 250
       if((step<t0) && algo=="ada")
         {ar_vs_length.push_back(chain->acc_rate(0,step)); }
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
251
       else if(step>=t0 && algo=="ada" )
Matthieu Tristram's avatar
Matthieu Tristram committed
252
         {ar_vs_length.push_back(chain->acc_rate(step-100,step)); }
Xavier Garrido's avatar
Xavier Garrido committed
253
       else
Matthieu Tristram's avatar
Matthieu Tristram committed
254
         {ar_vs_length.push_back(chain->acc_rate(0,step));}
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
255 256
     }

257
   chain->write_if_useful(ar_vs_length,ar_file.c_str());
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
258 259 260

   if(algo=="ada")
   {
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
261
     MCMC_adaptive* ada=dynamic_cast<MCMC_adaptive*>(chain);
262
     chain->write_if_useful(proposal->_scale_vs_length,scale_file.c_str());
Xavier Garrido's avatar
Xavier Garrido committed
263

264
     ofstream f(cor_file.c_str());
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
265
     if(!f) {
266 267
       cout << "Error in file creation!";
     }
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
268
     f << ada->mean() << endl //mean
269
       << proposal->Get_cov_matrix() << endl; //cov matrix
Xavier Garrido's avatar
Xavier Garrido committed
270
     f.close();
271
     cout<<"proposal cov written in " << cor_file <<endl;;
Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
272
   }
Xavier Garrido's avatar
Xavier Garrido committed
273

Plaszczynski Stephane's avatar
Plaszczynski Stephane committed
274 275 276
   delete chi2;

}