Commit 67b60f67 authored by Plaszczynski Stephane's avatar Plaszczynski Stephane
Browse files

begin cleaning of MCMC part

parent 9274cda9
......@@ -65,42 +65,28 @@ class MCMC_sampler{
HepVector _mean;
double _t0; //can be use to set a prephase of metropolis
double _ts; //can be use to set init of adapting scale
//private methods
//virtual
//virtual bool choose_next(const HepVector& x, const HepVector& y,double u,pdf* target_func=NULL,pdf* proposal_func=NULL)=0;
//virtual double calc_alpha(const HepVector& x,const HepVector& y,pdf* target_func=NULL,pdf* proposal_func=NULL)=0;
//public:
//default constructor
MCMC_sampler();
//costruttore
//constructors
MCMC_sampler();
MCMC_sampler(int dim_in, int length_in,bool verbose);
//distruttore
// ~MCMC_sampler();
//public methods
//public methods
virtual void fill_chain(pdf* target_func,pdf* proposal_func,int init_step,
const char* chain_run_name)=0;
double acc_rate(int length_0, int length);
void rearrange_chain();
void write_chain(char* str);
double acc_rate(int length_0, int length);
void rearrange_chain();
void write_chain(char* str);
void extract_cov(int burn_in,int stop, bool wr_on_file);
void GR_test_conv(int num_chain,pdf* target_func,pdf* proposal_func);
void estim_autocor(int min_lag,int max_lag, int length);
void compute_mean(int stop,HepVector mean);
void write_rearranged_chain(char* str); //rec, value, lkh
// void set_par(const char* str, int dim); // set parameters value reading from a file (to be rewrited)
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 GR_test_conv(int num_chain,pdf* target_func,pdf* proposal_func);
void estim_autocor(int min_lag,int max_lag, int length);
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);
//double target(vector<double>* a); //qui andranno aggiunti i parametri di target
//visualize parameters
int Get_dim();
......
......@@ -34,41 +34,30 @@ int main(int argc, char *argv[])
planck_assert(argc!=3,"Usage: mcmc parfile(in) samplesfile(out) ");
planck_assert(file_present(argv[1]),"missing par file");
// if(argc!=3)
// {
// cout << "wrong usage" <<endl;
// return EXIT_FAILURE;
// }
Parser parser(argv[1]);
char *chain_run_name=argv[2];
int dim=parser.params.find<int>("dim",1);
const string algo=parser.params.find<string>("algo","metro");
const int length=parser.params.find<int>("length",1);
//MCMC specific parameters
const string algo=parser.params.find<string>("algo","ada");
const int length=parser.params.find<int>("length",100);
const string target_cor_path=parser.params.find<string>("target_cov","");
// const bool use_ext_cor=parser.params.find<bool>("use_ext_cor",false);
// const bool use_ext_cov=parser.params.find<bool>("use_ext_cov",false);
const string proposal_cor=parser.params.find<string>("proposal_cor","");
const string proposal_cov=parser.params.find<string>("proposal_cov","");
const int t0=parser.params.find<int>("t0",length);
const int ts=parser.params.find<int>("ts",0);
double scale=parser.params.find<double>("scale",2.38*2.38/dim);
const int t0=parser.params.find<int>("t0",2000);
const int ts=parser.params.find<int>("ts",10000);
double scale=parser.params.find<double>("scale",0.001);
const bool do_move=parser.params.find<bool>("do_move",false);
bool verbose=parser.params.find<bool>("verbose",false);
const char* target_cor=target_cor_path.c_str();
const unsigned long bunchSize=parser.params.find<unsigned long>("bunchSize",1000);
const int seed=parser.params.find<int>("seed",length);
cout <<"dim=" <<dim <<endl; //proposal dimension
char *chain_run_name=argv[2];
bool use_ext_cor=false;
bool use_ext_cov=false;
if( proposal_cor != "") use_ext_cor = true;
if( proposal_cov != "") use_ext_cov = true;
const bool use_ext_cor=(proposal_cor.size()!=string::npos);
bool use_ext_cov=(proposal_cov.size()!=string::npos);
//---------------------------------------------------------
......@@ -77,10 +66,9 @@ int main(int argc, char *argv[])
chi2->setVerbose(false);
planck_lkh* lkh=new planck_lkh(chi2,parser);
dim=lkh->Get_par_num();
const int dim=lkh->Get_par_num();
cout <<"dim check=" <<dim <<endl;
cout <<"dim=" <<dim <<endl;
//------------------------------------------------------------------
//choosing the algorithm to use
......@@ -119,16 +107,6 @@ int main(int argc, char *argv[])
//chose the proposal-------------------------------------------------
HepVector mean(dim,0);
HepVector mean_pro(dim,0);
// ---------------------------------------
//prior stupidi per toy model
/*
HepVector priors_limit(dim,1);
chain->_min_par=-500*priors_limit;
chain->_max_par=500*priors_limit;
cout << chain->_min_par << endl;
cout << chain->_max_par << endl;
//chain->_starting_point=mean_pro;
*/
//-----------------------------------------
for(int i=0;i<dim;i++) {mean[i]=chain->_starting_point[i]; }
......@@ -137,67 +115,25 @@ int main(int argc, char *argv[])
//-----define the lkh----------------------------------------------
//lkh exemple
//pdf_log_gauss* lkh= new pdf_log_gauss(dim,mean,cov);
// lkh->set_cov(target_cor);
//pdf_gaussian* lkh= new pdf_gaussian(dim,mean,cov);
//cout <<lkh ->Get_cov_matrix() << endl;
pdf_gaussian* proposal= new pdf_gaussian(dim,mean_pro,identity);
if(use_ext_cov=true)
if(use_ext_cov)
{
cout << "update proposal cov matrix.." << endl;
proposal->set_cov(proposal_cov.c_str());
}
else {
if(use_ext_cor==true)
if(use_ext_cor)
{
cout << "update proposal cor matrix.." << endl;
proposal->set_cov(proposal_cor.c_str());
use_ext_cov = true;
use_ext_cov= true;
}
cov = proposal->Get_cov_matrix();
proposal->calc_cov_matrix( cov, par_err);
}
//------we can set the cov to the error in .par
//vector<double> par_err=(parser.upar).Errors();
//for(int i=0;i<dim;i++) { cout << "par err" << par_err[i] << endl;}
//setting proposal
// cov=proposal->Get_cov_matrix();
//proposal->calc_cov_matrix(cov, par_err);
// */
cout << "cov proposal= " << proposal->Get_cov_matrix() << endl;
// cout << "cov target= " <<lkh->Get_cov_matrix() << endl;
//in case of toy distribution you may want to directly shoot from the target
/*
//-----------------------------------------------------------------------------
pdf_gaussian* toy_target= new pdf_gaussian(dim,mean_pro,proposal->Get_cov_matrix());
toy_target->init_generator(time(NULL));
toy_target->Set_scale(1.0);
vector<double> toy_chain;
HepVector toy(dim);
bool do_shoot=toy_target->shoot(toy, true);
for(int i=0;i<length;i++)
{
do_shoot=toy_target->shoot(toy, false);
for(int j=0;j<dim;j++) {toy_chain.push_back(toy[j]+mean[j]);}
}
chain->write_if_useful(toy_chain,"toy_chain.txt");
//-----------------------------------------------------------------------------
*/
//scale
proposal->Set_scale(scale);
......@@ -221,7 +157,6 @@ int main(int argc, char *argv[])
sleep(1);
move->init_generator(time(NULL));
move->Set_scale(scale);
bool yes_move=move->shoot(start, true);
cout << "start=" << start << endl;
chain->Set_starting_point(chain->_starting_point+start);
cout << "start_new=" << chain->_starting_point << endl;
......@@ -241,7 +176,6 @@ int main(int argc, char *argv[])
prior* flat_p=new prior(prior_list);
//pdf* target=new posterior(lkh,flat_p);
pdf* target=new post_planck(lkh,flat_p);
......@@ -262,20 +196,7 @@ int main(int argc, char *argv[])
cout << "acc_rate=" << rate << endl ;
//rearrange chain gives me two vector.
//One with the value and the other with the recurrence
//chain->rearrange_chain();
//chain->write_if_useful(flat_p->_oor_vec_min,"chain_oor_min.txt");
//chain->write_if_useful(flat_p->_oor_vec_max,"chain_oor_max.txt");
// chain->write_chain("chain.txt");
//chain->extract_cov(length/2, true);
//---autocorr----and acc rate-------------
vector<double > ar_vs_length;
int max_lag;
int min_lag;
......@@ -284,8 +205,6 @@ int main(int argc, char *argv[])
{
max_lag=50;
min_lag=max_lag;
//chain->estim_autocor(min_lag,max_lag,step);
//autocorr_vs_length.push_back(chain->_autocor_time);
if((step<t0) && algo=="ada"){ar_vs_length.push_back(chain->acc_rate(0,step)); }
else if(step>=t0 && algo=="ada" )
{ar_vs_length.push_back(chain->acc_rate(step-100,step)); }
......@@ -305,27 +224,14 @@ int main(int argc, char *argv[])
ofstream f("corr.txt");
if(!f) {
cout << "Error in file creation!";
}
f << chain->_mean << endl //mean
<< proposal->Get_cov_matrix() << endl; //cov matrix
f.close();
cout<<"output written in:corr.txt"<<endl;;
cout << "Error in file creation!";
}
f << chain->_mean << endl //mean
<< proposal->Get_cov_matrix() << endl; //cov matrix
f.close();
cout<<"output written in:corr.txt"<<endl;;
}
//cout << "autocorr="<<chain->_autocor_time << endl;
//-----------------------------------------------
//GR test
//int num_chain=4;
//chain->GR_test_conv(num_chain, target,proposal);
//char* name="rearranged_chain.txt";
//chain->write_rearranged_chain(name);
delete chi2;
......
......@@ -62,9 +62,6 @@ void dumpMin(ostream& o,const FunctionMinimum& min1){
o << "Minimum is " << (min1.IsValid() ? "" : "NOT") << " valid" <<endl;
o << "**********************************************************" << endl;
o << "User covariance " << endl;
o << min1.UserCovariance() << endl;
o <<"Minimum: "<<min1 <<std::endl;
}
......
......@@ -105,9 +105,9 @@ fi
jobsub="${QSUB_CMD} -o $PWD -N $(basename $dirout) -t 1-$Ntir camelrun"
echo "about to run : $jobsub"
echo "OK? [o/n]"
echo "OK? [y/n]"
read answer
if [ $answer != 'o' ] ; then
if [ $answer != 'y' ] ; then
echo "exiting"
exit 1
fi
......
......@@ -134,9 +134,9 @@ fi
jobsub="${QSUB_CMD} -o $PWD -N $(basename $dirout) -t 1-$Ntir camelrun"
echo "about to run : $jobsub"
echo "OK? [o/n]"
echo "OK? [y/n]"
read answer
if [ $answer != 'o' ] ; then
if [ $answer != 'y' ] ; then
echo "exiting"
exit 1
fi
......
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