Skip to content

GitLab

Commit 671f5c99 authored by Plaszczynski Stephane's avatar Plaszczynski Stephane
Browse files

fix the chunck offset+clean fstream: text OK unless t0+ts<100

parent 217869ae
Hide whitespace changes
Inline Side-by-side
Showing with 102 additions and 62 deletions
cmt/requirements-gcc
View file @ 671f5c99
......@@ -75,6 +75,7 @@ macro_append cppflags ' -DRELPATH=\"$(CAMELROOT)/lik\" '
#CAMEL applications
# pour relinker si changement: toutes les application du package
macro CAMEL_stamps $(PACKAGE_ROOT)/$(CMTCONFIG)/MinuitFit.stamp
macro CAMEL_stamps $(PACKAGE_ROOT)/$(CMTCONFIG)/MCMC.stamp
macro application_suffix ""
#groupe exec
......@@ -97,7 +98,8 @@ macro_append Minimize_dependencies " MinuitFit "
macro_append ScanParam_dependencies " MinuitFit "
macro_append writeChi2_dependencies " MinuitFit "
macro_append writeSpectra_dependencies " MinuitFit "
macro_append mcmc_dependencies " MinuitFit MCMC"
macro_append mcmc_dependencies " MinuitFit"
macro_append mcmc_dependencies " MCMC"
macro_append testKlass_dependencies " MinuitFit"
macro_append clik_example_CC_dependencies " MinuitFit "
macro_append test_clik_dependencies " MinuitFit "
......
cmt/requirements-icc
View file @ 671f5c99
......@@ -60,6 +60,7 @@ macro_append cppflags ' -DRELPATH=\"$(CAMELROOT)/lik\" '
#CAMEL applications
# pour relinker si changement: toutes les application du package
macro CAMEL_stamps $(PACKAGE_ROOT)/$(CMTCONFIG)/MinuitFit.stamp
macro CAMEL_stamps $(PACKAGE_ROOT)/$(CMTCONFIG)/MCMC.stamp
macro application_suffix ""
#groupe exec
......@@ -83,7 +84,8 @@ macro_append Minimize_dependencies " MinuitFit "
macro_append ScanParam_dependencies " MinuitFit "
macro_append writeChi2_dependencies " MinuitFit "
macro_append writeSpectra_dependencies " MinuitFit "
macro_append mcmc_dependencies " MinuitFit MCMC"
macro_append mcmc_dependencies " MinuitFit "
macro_append mcmc_dependencies " MCMC"
macro_append testKlass_dependencies " MinuitFit"
macro_append clik_example_CC_dependencies " MinuitFit "
macro_append test_clik_dependencies " MinuitFit "
......
src/camel/MCMC/MCMC_adaptive.cc
View file @ 671f5c99
......@@ -13,6 +13,7 @@ MCMC_adaptive::MCMC_adaptive(int dim_in,int length_in, unsigned int seedU,bool v
HepVector col(dim_in);
vector<HepVector > matrix(length_in,col);
_mychain=matrix;
cout << "chain with " << _mychain[0].num_row() << " variables of size="<< _mychain.size() << endl;
}
......@@ -80,33 +81,36 @@ void MCMC_adaptive::fill_chain(pdf* target_func, pdf* proposal,int init_step, co
{
post_planck* pl=dynamic_cast<post_planck *>(target_func);
if(pl!=0) cout << "planck likelihood" << endl;
ofstream chain_running(chain_run_name);
///open file to write chain while running
fstream chain_running;
chain_running.open (chain_run_name, fstream::in | fstream::out | fstream::app);
chain_running.precision(8);
//dump header
chain_running << "chi2" << "\t";
for(int l=0;l<_par_name.size();l++) { chain_running << _par_name[l] << "\t";}
chain_running << endl;
chain_running << endl;
chain_running.close();
//output in scientific format
chain_running << scientific;
//chain_running.precision(8);
_add_flag=false;
if(pl && (pl->_likelihood)->_parser.derivedVars().size()!=0)
{
_add_flag=true;
cout << "there will be derived par chain" << endl;
}
if (_verbose) cout << "there will be derived par chain" << endl;
}
//dynamic cast
pdf_gaussian* g=dynamic_cast<pdf_gaussian *>(proposal);
if(g!=0) cout << "gaussian proposal" << endl;
if(g==0) cout << "NON gaussian proposal???" << endl;
_mychain.at(0)=_starting_point; //to be set before calling fill
target_func->evaluate(_starting_point);
cout << "target= " << target_func->_func_value << endl;
if (_verbose) cout << "target= " << target_func->_func_value << endl;
if(_add_flag){
_add_chain.push_back((pl->_likelihood)->_derived_par);
}
......@@ -122,7 +126,7 @@ void MCMC_adaptive::fill_chain(pdf* target_func, pdf* proposal,int init_step, co
vector<double> previous_derived;//for extra par
vector<double> candidate_derived;//for extra par
cout << "I made the first evaluation" << endl;
if (_verbose) cout << "I made the first evaluation" << endl;
//these are only useful for updating covariance
......@@ -181,8 +185,12 @@ void MCMC_adaptive::fill_chain(pdf* target_func, pdf* proposal,int init_step, co
}
cout << "this was my first shoot" << endl;
if (_verbose)
cout << "this was my first shoot" << endl;
int nextind(0);
///////////////////////////////////////main loop////////////////////////
for(int i=2; i<_length; i++)
{
if((_verbose==true)&&(i%1000==0))
......@@ -204,12 +212,12 @@ void MCMC_adaptive::fill_chain(pdf* target_func, pdf* proposal,int init_step, co
//mean to be computed since here
_mean=previous;
compute_mean(i-1,_mean);
cout << "mean="<< _mean << endl;
if (_verbose) cout << "mean="<< _mean << endl;
//first extimation of covariance to be redone
//cout << "matrix_0="<<g->Get_cov_matrix() << endl;
extract_cov(0,i-1,false);
g->Set_cov_matrix(_extracted_cov);
cout << "matrix_t0="<<g->Get_cov_matrix() << endl;
if (_verbose) cout << "matrix_t0="<<g->Get_cov_matrix() << endl;
flag_shoot=proposal->shoot(jump, true);
if(_dim==1) g->Set_scale(0.1);
else g->Set_scale(2.4*2.4/_dim);
......@@ -252,8 +260,7 @@ void MCMC_adaptive::fill_chain(pdf* target_func, pdf* proposal,int init_step, co
if((_verbose==true)||(i%1000==0))
{
cout<<"acc rate="<<ar<<endl;
}
// g->Set_scale(g->Get_scale()*(1-gamma_i)+ (ar-0.25)*gamma_i);
}
g->Set_scale(g->Get_scale()+ (ar-0.25)*gamma_i);
double scale_min=0.001;
if(scale<scale_min) {g->Set_scale(scale_min);}
......@@ -314,33 +321,38 @@ void MCMC_adaptive::fill_chain(pdf* target_func, pdf* proposal,int init_step, co
//for writing on the run
int gate=_bunchSize;
if(i%gate==0)
if(i%gate==0)
{
if (_verbose) cout << "dump at index " << i << " from " << i-gate << " to " << i-1 << endl;
chain_running.open (chain_run_name, fstream::app);
for(int k=i-gate; k<i; k++)
{
chain_running.open (chain_run_name, fstream::in | fstream::out | fstream::app);
for(int k=i-gate; k<i; k++)
{
//chain_running.setf( std::ios::scientific);
chain_running << _lkh_values[k] << "\t";
//chain_running.precision(8);
// chain_running.setf( std::ios::fixed,std:: ios::floatfield );
for(int j=0; j<_dim; j++) chain_running << _mychain.at(k)[j] << "\t";
if(_add_flag) for(int l=0; l< (_add_chain.at(0).size()); l++) chain_running << _add_chain.at(k)[l] << "\t";
chain_running << endl;
}
chain_running.close();
if(_verbose==true)
{
cout << "output written in running chain" <<endl;
if(_add_flag) cout << "extra pars written in added to chain" <<endl;
}
}//end of for for writing
chain_running << _lkh_values[k] << "\t";
for(int j=0; j<_dim; j++) chain_running << _mychain.at(k)[j] << "\t";
if(_add_flag) for(int l=0; l< (_add_chain.at(0).size()); l++) chain_running << _add_chain.at(k)[l] << "\t";
chain_running << endl;
}
chain_running.close();
nextind=i;
}//end of for for writing
}
//if missing data to dump
if (nextind != _mychain.size()){
if (_verbose) cout << "complete chunck" << endl;
chain_running.open (chain_run_name, fstream::app);
for(int k=nextind; k<_mychain.size(); k++)
{
chain_running << _lkh_values[k] << "\t";
for(int j=0; j<_dim; j++) chain_running << _mychain.at(k)[j] << "\t";
if(_add_flag) for(int l=0; l< (_add_chain.at(0).size()); l++) chain_running << _add_chain.at(k)[l] << "\t";
chain_running << endl;
}
chain_running.close();
}
}
src/camel/MCMC/MCMC_sampler.cc
View file @ 671f5c99
......@@ -422,10 +422,11 @@ void MCMC_sampler::set_par(Parser parser)
}
for(int l=0;l<parser.derivedVars().size();l++) _par_name.push_back(parser.derivedVars()[l]);
/*
for(int l=0;l<_par_name.size();l++) cout << _par_name[l] <<"\t";
cout << endl;
cout << _min_par << _max_par << _starting_point << _par_err <<endl;
*/
}
void
MCMC_sampler::write_samples_fits(const char* fn) {
......@@ -440,17 +441,34 @@ MCMC_sampler::write_samples_fits(const char* fn) {
std::vector<fitscolumn> cols;
//define columns
//first is chi2
cols.push_back(fitscolumn("chi2", "no",1,TDOUBLE));
for (size_t i=0;i< _par_name.size();i++){
cols.push_back(fitscolumn(_par_name[i], "?",1,TDOUBLE));
cols.push_back(fitscolumn(_par_name[i], "?",1,TDOUBLE));
}
fout.insert_bintab (cols);
//convert to arr and write staring at index 1
//write chi2
arr<double> tab(_lkh_values.size());
for (size_t j=0;j<tab.size();j++) tab[j]= _lkh_values[j];
fout.write_column(1,tab);
//convert to arr and write starting at index 2
//need to transpoise into columns first
const int N=_mychain.size();
cout <<"chain length="<< N << endl;
vector< arr<double> > out(_dim);
for (size_t ivar=0;ivar<_dim;ivar++){
out[ivar].alloc(N);
for (size_t k=0;k<N;k++) out[ivar][k]=_mychain[k][ivar];
}
for (size_t i=0;i< _par_name.size();i++){
HepVector& vec=_value_chain[i];
arr<double> tab(vec.num_row());
for (size_t j=0;j<tab.size();j++) tab[j]=vec[i];
fout.write_column(i+1,tab);
cout << "writing : parname=" << _par_name[i] << " size=" << out[i].size() << endl;
fout.write_column(i+2,out[i]);
}
fout.close();
......
src/camel/MCMC/exec/mainMCMC.cc
View file @ 671f5c99
......@@ -41,8 +41,13 @@ int main(int argc, char *argv[])
Parser parser(argv[1]);
char *chain_run_name=argv[2];
const char* chain_run_name=argv[2];
const string fileName(chain_run_name);
if (file_present(fileName)) {
std::cout << "output file " + fileName + " already exists: removing it"<<endl;
remove_file(fileName);
}
//MCMC specific parameters
const string algo=parser.params.find<string>("algo","ada");
......@@ -113,12 +118,13 @@ int main(int argc, char *argv[])
chain->SetBunchSize(bunchSize);
//directly from Parser
cout << "using Parser for setting par" << endl;
chain->set_par(parser);
cout << chain->min_par() << endl;
cout << chain->max_par() << endl;
cout << chain->starting_point() << endl;
if (verbose){
cout << chain->min_par() << endl;
cout << chain->max_par() << endl;
cout << chain->starting_point() << endl;
}
//or by hand using set/get functions
......@@ -126,7 +132,7 @@ int main(int argc, char *argv[])
for(int i=0;i<dim;i++)
{
par_err[i]=(chain->par_err())[i];
cout << par_err[i] <<endl;
if (verbose) cout << par_err[i] <<endl;
}
//chose the proposal-------------------------------------------------
......@@ -158,11 +164,11 @@ int main(int argc, char *argv[])
proposal->calc_cov_matrix( cov, par_err);
}
cout << "cov proposal= " << proposal->Get_cov_matrix() << endl;
if (verbose) cout << "cov proposal= " << proposal->Get_cov_matrix() << endl;
//scale
proposal->Set_scale(scale);
cout <<"scale=" << proposal->Get_scale() << endl;
if (verbose) cout <<"scale=" << proposal->Get_scale() << endl;
vector<double> scale_vs_length;
proposal->init_generator(seed);
......@@ -175,14 +181,14 @@ int main(int argc, char *argv[])
{
pdf_gaussian* move= new pdf_gaussian(dim,mean_pro,identity);
move->calc_cov_matrix(identity, par_err);
cout << "cov domove= " <<move->Get_cov_matrix() << endl;
if (verbose) cout << "cov domove= " <<move->Get_cov_matrix() << endl;
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);
cout << "start=" << start << endl;
if (verbose) cout << "start=" << start << endl;
chain->Set_starting_point(chain->starting_point()+start);
cout << "start_new=" << chain->starting_point() << endl;
if (verbose) cout << "start_new=" << chain->starting_point() << endl;
}
......@@ -217,7 +223,7 @@ int main(int argc, char *argv[])
//acceptance rate simple
double rate=chain->acc_rate(0,length);
cout << "acc_rate=" << rate << endl ;
if (verbose) cout << "acc_rate=" << rate << endl ;
//---autocorr----and acc rate-------------
......
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