Commit 5ca0e3c6 authored by Plaszczynski Stephane's avatar Plaszczynski Stephane
Browse files

explicit the 2 seeds in MCMC

parent e2a37d04
......@@ -6,8 +6,8 @@
using namespace std;
MCMC_adaptive::MCMC_adaptive(int dim_in,int length_in, bool verbose):
MCMC_sampler(dim_in,length_in,verbose)
MCMC_adaptive::MCMC_adaptive(int dim_in,int length_in, unsigned int seedU,bool verbose):
MCMC_sampler(dim_in,length_in,seedU,verbose)
{
HepVector col(dim_in);
......@@ -131,7 +131,7 @@ void MCMC_adaptive::fill_chain(pdf* target_func, pdf* proposal,int init_step, co
HepSymMatrix new_cov(_dim);
//creating random generator
AbsFlat *flat=new planck_rng();
AbsFlat *flat=new planck_rng(0,1,_seedU);
vector<double> u(_length-1);
flat ->fill(u);
......
......@@ -27,7 +27,7 @@ protected:
public:
//constructor
MCMC_adaptive(int dim_in, int length_in, bool verbose);
MCMC_adaptive(int dim_in, int length_in, unsigned int seedU,bool verbose);
virtual void fill_chain(pdf* target_func,pdf* proposal,int init_step,
......
......@@ -7,8 +7,8 @@
using namespace std;
MCMC_metropolis::MCMC_metropolis(int dim_in,int length_in, bool verbose):
MCMC_sampler(dim_in,length_in, verbose)
MCMC_metropolis::MCMC_metropolis(int dim_in,int length_in, unsigned int seedU,bool verbose):
MCMC_sampler(dim_in,length_in, seedU,verbose)
{
HepVector col(dim_in,0.0);
......@@ -119,7 +119,7 @@ void MCMC_metropolis::fill_chain(pdf* target_func, pdf* proposal,int init_step,c
vector<double> previous_derived;//for extra par
vector<double> candidate_derived;//for extra par
AbsFlat *flat=new planck_rng();
AbsFlat *flat=new planck_rng(0,1,_seedU);
vector<double> u(_length);
flat ->fill(u);
......
......@@ -14,7 +14,7 @@ class MCMC_metropolis : public MCMC_sampler
public:
//costruttore
MCMC_metropolis(int dim_in, int length_in, bool verbose);
MCMC_metropolis(int dim_in, int length_in, unsigned int seedU, bool verbose);
virtual void fill_chain(pdf* target_func,pdf* proposal,int init_step,const char* chain_run_name);
......
......@@ -18,9 +18,9 @@ MCMC_sampler::MCMC_sampler()
{}
MCMC_sampler::MCMC_sampler(int dim_in,int length_in, bool verbose):
MCMC_sampler::MCMC_sampler(int dim_in,int length_in, unsigned int seed,bool verbose):
_dim(dim_in), _length(length_in), _verbose(verbose),_bunchSize(1000)
_dim(dim_in), _length(length_in), _seedU(seed),_verbose(verbose),_bunchSize(1000)
{
HepVector col(dim_in,0.0);
vector<HepVector > matrix(length_in,col);
......
......@@ -45,11 +45,14 @@ class MCMC_sampler{
HepVector _GR_param;
std::vector<HepVector > _autocorr;
HepVector _autocor_time;
//needs a seed for uniform shoot
unsigned int _seedU;
public:
//constructors
MCMC_sampler();
MCMC_sampler(int dim_in, int length_in,bool verbose);
MCMC_sampler(int dim_in, int length_in,unsigned int seedU,bool verbose);
//public methods
......@@ -95,8 +98,11 @@ public:
void Set_recur_chain(std::vector<int> a);
void Set_starting_point(HepVector a);
void SetBunchSize(const unsigned long b){_bunchSize=b;}
template<typename T> inline void SetUseed(const T& s){
_seedU=static_cast<unsigned int>(std::abs(s));
}
//verbose
void Set_verbose(bool verbose) {_verbose=verbose;}
......
......@@ -25,7 +25,10 @@
#include <iostream>
#include <fstream>
#include <vector>
#include <unistd.h>
#include <ctime>
#include <stdint.h>
using namespace std;
int main(int argc, char *argv[])
......@@ -54,13 +57,11 @@ int main(int argc, char *argv[])
const unsigned long bunchSize=parser.params.find<unsigned long>("bunchSize",1000);
const int seed=parser.params.find<int>("seed",length);
const bool use_ext_cor=(proposal_cor.size()!=0);
bool use_ext_cor=(proposal_cor.size()!=0);
bool use_ext_cov=(proposal_cov.size()!=0);
//---------------------------------------------------------
Chi2Combiner* chi2=Chi2Factory::gimeChi2(parser);
chi2->setVerbose(false);
......@@ -70,18 +71,40 @@ int main(int argc, char *argv[])
cout <<"dim=" <<dim <<endl;
//--------defining seeds-----------------------------------------------
unsigned int seed(0);
if (parser.params.param_present("seed")){
seed=parser.params.find<int>("seed",1234);
}
else {//pick automatiocally
time_t result = time(NULL);
seed=static_cast<uint32_t>(result);
}
cout << "seed for proposal=" << seed << endl;
unsigned int seedU(0);
//need anothe rone for the algorithm itsel
if (parser.params.param_present("seedU")){
seedU=parser.params.find<int>("seedU",1234);
}
else {//pick automatiocally: time too close to previous one so divide by 2
time_t result = time(NULL);
seedU=static_cast<uint32_t>(result/2);
}
cout << "seed for algorithm=" << seedU << endl;
//------------------------------------------------------------------
//choosing the algorithm to use
//for the moment just metropolis or adaptive.Metropolis is the default one
MCMC_sampler* chain;
if(algo=="ada")
{
chain= new MCMC_adaptive(dim,length,verbose);
chain= new MCMC_adaptive(dim,length,seedU,verbose);
cout << "using Adaptive algorithm" << endl;
}
else
{
chain= new MCMC_metropolis(dim,length,verbose);
chain= new MCMC_metropolis(dim,length,seedU,verbose);
cout << "using Metropolis algorithm" << endl;
}
......@@ -140,8 +163,7 @@ int main(int argc, char *argv[])
cout <<"scale=" << proposal->Get_scale() << endl;
vector<double> scale_vs_length;
// proposal->init_generator(time(NULL));
proposal->init_generator(seed+1234);
proposal->init_generator(seed);
......
......@@ -25,37 +25,35 @@ void prior_flat::Set_max(double a) { _max=a;}
//methods no longer virtual
void prior_flat::evaluate(const HepVector& value)
{
if(value[0]<_min)
{
_func_value=0;
_oor_flag_min=true;
}
else if(value[0]>_max)
{
_func_value=0;
_oor_flag_min=false;
}
else
{
_func_value=1;
}
}
void prior_flat::init_generator()
{
_flat =new planck_rng();
}
{
if(value[0]<_min)
{
_func_value=0;
_oor_flag_min=true;
}
else if(value[0]>_max)
{
_func_value=0;
_oor_flag_min=false;
}
else
{
_func_value=1;
}
}
void prior_flat::init_generator(int seed)
{
_flat =new planck_rng(_min,_max,seed);
}
bool prior_flat::shoot(HepVector& jump, bool flag)
{
double x= _flat->next();
jump[0]=x*(_max-_min)+_min;
jump[0]= _flat->next();
return true;
}
......
......@@ -20,7 +20,7 @@ public:
//method no longer virtual
void evaluate(const HepVector& value);
void init_generator();
void init_generator(int seed);
bool shoot(HepVector& jump, bool flag);
//visualize parameters
......
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