Commit e32fdef4 authored by Plaszczynski Stephane's avatar Plaszczynski Stephane
Browse files

continue cleaning

parent 6d89f2ee
#include "MCMC_adaptive.hh"
#include "AbsRand/AbsFlat.hh"
#include "AbsRand/planck_rng.hh"
#include <fstream>
using namespace std;
......
#include "MCMC_metropolis.hh"
#include "pdf.hh"
#include "AbsRand/AbsFlat.hh"
#include "AbsRand/planck_rng.hh"
#include <fstream>
using namespace std;
......
#include "MCMC_sampler.hh"
#include "AbsRand/planck_rng.hh"
#include "AbsRand/GaussW.hh"
#include "pdf.hh"
#include <vector>
#include <stdlib.h>
......
......@@ -8,8 +8,6 @@
#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include <cmath>
#include <cstdio>
#include <vector>
#include <string>
......@@ -48,7 +46,7 @@ class MCMC_sampler{
std::vector<HepVector > _autocorr;
HepVector _autocor_time;
public:
//constructors
MCMC_sampler();
MCMC_sampler(int dim_in, int length_in,bool verbose);
......@@ -69,8 +67,9 @@ class MCMC_sampler{
void write_if_useful(std::vector<double>& vec,char* str);
void write_vecHep(std::vector<HepVector >& vec,char* str);
//visualize parameters
//accessors
int Get_dim();
int Get_length();
int Get_burn_in();
......@@ -79,6 +78,11 @@ class MCMC_sampler{
std::vector<HepVector > Get_value_chain();
std::vector<int> Get_recur_chain();
HepMatrix Get_extracted_cov();
const HepVector& min_par() {return _min_par;}
const HepVector& max_par() {return _max_par;}
const HepVector& starting_point() {return _starting_point;}
const HepVector& par_err() {return _par_err;}
//setting parameters
void Set_dim(int a);
......
......@@ -91,16 +91,16 @@ int main(int argc, char *argv[])
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;
cout << chain->min_par() << endl;
cout << chain->max_par() << endl;
cout << chain->starting_point() << endl;
//or by hand using set/get functions
vector<double> par_err(dim);
for(int i=0;i<dim;i++)
{
par_err[i]=chain->_par_err[i];
par_err[i]=(chain->par_err())[i];
cout << par_err[i] <<endl;
}
......@@ -108,7 +108,7 @@ int main(int argc, char *argv[])
HepVector mean(dim,0);
HepVector mean_pro(dim,0);
//-----------------------------------------
for(int i=0;i<dim;i++) {mean[i]=chain->_starting_point[i]; }
for(int i=0;i<dim;i++) {mean[i]=(chain->starting_point())[i]; }
HepSymMatrix cov(dim,1);
HepSymMatrix identity(dim,1);
......@@ -158,8 +158,8 @@ int main(int argc, char *argv[])
move->init_generator(time(NULL));
move->Set_scale(scale);
cout << "start=" << start << endl;
chain->Set_starting_point(chain->_starting_point+start);
cout << "start_new=" << chain->_starting_point << endl;
chain->Set_starting_point(chain->starting_point()+start);
cout << "start_new=" << chain->starting_point() << endl;
}
......@@ -169,7 +169,7 @@ int main(int argc, char *argv[])
vector<pdf*> prior_list(dim);
for(int i=0; i<dim; i++)
{
pdf* p=new prior_flat(chain->_min_par[i],chain->_max_par[i]);
pdf* p=new prior_flat((chain->min_par())[i],(chain->max_par())[i]);
//pdf* p=new prior_beta(chain->_min_par[i],chain->_max_par[i], 1,1);
prior_list[i]=p;
}
......
#include "pdf.hh"
#include <vector>
#include <stdlib.h>
#include <iostream>
using namespace std;
......
#ifndef __pdf_class_guard__
#define __pdf_class_guard__
#include <math.h>
#include <stdio.h>
#include <vector>
#include "AbsRand.hh"
#include "AbsFlat.hh"
#include "AbsGauss.hh"
#include "Drand48.hh"
#include "planck_rng.hh"
#include "GaussW.hh"
#include "AbsRand/AbsRand.hh"
#include "CLHEP/Matrix/Vector.h"
class pdf{
......
#include "pdf_chi_square.hh"
#include <math.h>
#include <cmath>
#include "CLHEP/Matrix/Vector.h"
#include "AbsRand/planck_rng.hh"
#include "AbsRand/GaussW.hh"
using namespace std;
//default constructor
......
#include "pdf_gaussian.hh"
#include <math.h>
#include <iostream>
#include <fstream>
#include "CLHEP/Matrix/Vector.h"
#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include "AbsRand/planck_rng.hh"
#include "AbsRand/GaussW.hh"
#include <cmath>
#include <cassert>
#include <iostream>
#include <fstream>
using namespace std;
......
......@@ -3,9 +3,13 @@
#include "pdf.hh"
#include "CLHEP/Matrix/Vector.h"
#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include <iostream>
class pdf_gaussian : public pdf
......
#include "pdf_log_gauss.hh"
#include <math.h>
#include <iostream>
#include <fstream>
#include "CLHEP/Matrix/Vector.h"
#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include "AbsRand/planck_rng.hh"
#include "AbsRand/GaussW.hh"
#include <cmath>
#include <iostream>
#include <fstream>
#include<cassert>
using namespace std;
......
#include "pdf_log_gemma.hh"
#include <math.h>
#include <iostream>
#include <fstream>
#include "CLHEP/Matrix/Vector.h"
#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/SymMatrix.h"
using namespace std;
//default constructor
pdf_log_gemma::pdf_log_gemma(int dim_gauss_in,const HepSymMatrix& cov_matrix_in):
_dim_gauss(dim_gauss_in),_cov_matrix(cov_matrix_in), _scale(1), _t_func(0)
{
int ifail=0;
_inv_cov=_cov_matrix;
_inv_cov.invert(ifail);
assert(ifail==0);
}
//visualize parameters
HepVector pdf_log_gemma::Get_mean() {return _mean;}
HepSymMatrix pdf_log_gemma::Get_cov_matrix() {return _cov_matrix;}
int pdf_log_gemma::Get_dim_gauss() {return _dim_gauss;}
int pdf_log_gemma::Get_t_func() {return _t_func;}
HepMatrix pdf_log_gemma::Get_L() {return _L;}
double pdf_log_gemma::Get_scale() {return _scale;}
//setting parameters
void pdf_log_gemma::Set_mean(HepVector a) {_mean=a;}
void pdf_log_gemma::Set_cov_matrix(HepSymMatrix a) {_cov_matrix=a;}
void pdf_log_gemma::Set_dim_gauss(int a) {_dim_gauss=a;}
void pdf_log_gemma::Set_t_func(int a) {_t_func=a;}
void pdf_log_gemma::Set_L(HepMatrix a) {_L=a;}
void pdf_log_gemma::Set_scale(double a) { _scale=a;}
/*!
This is a simple 5 points Newton-Cotes formula that allows numerical integration of the function given by x
and y. This (simple) algorithm requires the x values to be equally spaced and in increasing order. It is very
accurate as long as the function is well sampled (see Numerical Recipes in C for details).
*/
double pdf_log_gemma::integrate_nc5(const vector<double> & x, const vector<double> & y)
{
// This is a five points Newton-Cotes (Bode's formula) integrator
// See NumRec for details
// We assume that the data is regularly gridded
unsigned int npts = x.size();
double h = x[1]-x[0];
vector<unsigned int> ii;
unsigned int nbii = (unsigned int)floor((npts-1.)/4);
unsigned int rest = (npts-1)-nbii*4;
unsigned int nbii2;
if(rest == 1 || rest == 2) nbii2 = nbii-1;
else nbii2 = nbii;
ii.resize(nbii2);
for(unsigned int i=0; i<nbii2; i++) ii[i] = (i+1)*4;
double integral = 0;
for(unsigned int i=0; i<nbii2; i++)
{
integral += 2.*h*(7.*(y[ii[i]-4]+y[ii[i]])+32.*(y[ii[i]-3]+y[ii[i]-1])+12.*y[ii[i]-2])/45.;
}
if( rest+1 == 2 )
{
// decoupage 4-3
unsigned int shift = npts-2;
// 4 points
integral += 3*h*(y[shift-4]+3*y[shift-3]+3*y[shift-2]+y[shift-1])/8.;
// 3 points
integral += h*(y[npts-3]+4*y[npts-2]+y[npts-1])/3.;
}
else if( rest+1 == 3 )
{
// decoupage 4-4
unsigned int shift = npts-3;
// 4 points
integral += 3*h*(y[shift-4]+3*y[shift-3]+3*y[shift-2]+y[shift-1])/8.;
// 3 points
shift = npts;
integral += 3*h*(y[shift-4]+3*y[shift-3]+3*y[shift-2]+y[shift-1])/8.;
}
else if(rest+1 == 4)
{
integral += 3*h*(y[npts-4]+3*y[npts-3]+3*y[npts-2]+y[npts-1])/8.;
}
return integral;
}
//function for T(z)
void pdf_log_gemma::t_z(const HepVector& value)
{
//order of paameters
//omega_m, k, h100, w_0,w_a
//cout << "value=" << value << endl;
HepVector y(_dim_gauss);
HepVector integrated(_dim_gauss);
double factor=3.0*value[1]*(1.0-value[0])/((2.47E-05)/(value[2]*value[2]));
//cout << "factor=" << factor << endl;
//cout << " dim gauss= " << _dim_gauss <<endl;
for(int i=0; i < _dim_gauss ; i++)
{
integrated[i]=0;
vector<double> x(_steps+1);
vector<double> to_integer(_steps+1);
for(int j=0; j <= _steps; j++)
{
x[j]=j*_z[i]/double(_steps);
//cout << " x=" << x[j]<< endl;
to_integer[j]=pow((1.0+x[j]),3*(value[1]+value[3]+value[4]))*pow((1.0+x[j]),-2)*exp(-3.0*value[4]*x[j]/(1.0+x[j]));
//cout << " to_int=" << to_integer[j]<< endl;
}
integrated[i]=integrate_nc5(x,to_integer);
//cout << "integral=" << integrated[i] << endl;
y[i]=pow(1.0-factor*integrated[i],0.25);
//cout << "y=" << y[i] << endl;
}
//fill HepVector _t_guess in the end.
_t_guess=y;
//cout << "_t_guess ready" << endl;
for(int i=0; i < _dim_gauss ; i++)
{
_t_guess[i]=_T_0*(1.0 + _z[i])*y[i];
//cout << _t_guess[i] << ", ";
//cout << y[i] << ", ";
}
}
//low z approx
void pdf_log_gemma::t_lowz(const HepVector& value)
{
//order of paameters
//omega_m, k, h100, w_0,w_a
//cout << "value=" << value << endl;
HepVector y(_dim_gauss);
HepVector integrated(_dim_gauss);
double factor=3.0*value[1]*(1.0-value[0])/((2.47E-05)/(value[2]*value[2]))/4.;
//cout << "factor=" << factor << endl;
//cout << " dim gauss= " << _dim_gauss <<endl;
for(int i=0; i < _dim_gauss ; i++)
{
y[i]=1.0-factor*_z[i];
//cout << "y=" << y[i] << endl;
}
//fill HepVector _t_guess in the end.
_t_guess=y;
//cout << "_t_guess ready" << endl;
for(int i=0; i < _dim_gauss ; i++)
{
_t_guess[i]=_T_0*(1.0 + _z[i])*y[i];
//cout << _t_guess[i] << ", ";
//cout << y[i] << ", ";
}
}
//simple model
void pdf_log_gemma::t_beta(const HepVector& value)
{
//---------------------------------
//if I only need T(z)=T_0*(i+z)^(1-beta)
HepVector z(_dim_gauss);
//fill HepVector _t_guess
_t_guess=z;
for(int i=0; i < _dim_gauss ; i++)
{
_t_guess[i]=_T_0*pow((1.0 + _z[i]),1-value[0]);
//cout << "value= " << value[0] << ",";
}
}
//methods no longer virtual
void pdf_log_gemma::evaluate(const HepVector& value)
{
//every time I call this function _t_guess is updated!
if(_t_func==1) t_beta(value);
else if(_t_func==2) t_lowz(value);
else t_z(value);
//cout << "t_z called!" << endl;
HepVector copy_value=_t_guess-_mean; //N.B here mean is data!
double beta=_inv_cov.similarity(copy_value);
//the multivariate gaussian is not normalised. No need in MCMC
_func_value=beta;
_lkh_value=_func_value; //this is needed for generic case
}
void pdf_log_gemma::init_generator()
{
_normal=new GaussW(new planck_rng(),0,1);
}
void pdf_log_gemma::calc_corr_matrix(HepSymMatrix& corr_matrix)
{
corr_matrix=_cov_matrix;
for(int i=0; i < _dim_gauss ; i++)
{
for(int j=0; j < _dim_gauss ; j++)
{
corr_matrix[i][j]=_cov_matrix[i][j]/(sqrt(_cov_matrix[i][i])*sqrt(_cov_matrix[j][j]));
}
}
//cout << "cov_matrix="<< _cov_matrix << endl;
}
//to shoot from a multivariate gaussian we need Cholesky's decomposition
void pdf_log_gemma::Cholesky_decompose(const HepSymMatrix& _cov_matrix)
{
//transform in corr matrix
HepSymMatrix corr_matrix(_dim_gauss);
calc_corr_matrix(corr_matrix);
HepSymMatrix cov_matrix=corr_matrix;
//cout << "corr_matrix="<< cov_matrix << endl;
HepMatrix prova(_dim_gauss,_dim_gauss);
_L=prova;
_L[0][0]=sqrt(cov_matrix[0][0]);
double ss=0;
double sp=0;
for(int i=1; i < _dim_gauss ; i++)
{
_L[i][0]=cov_matrix[i][0]/_L[0][0];
}
for(int j=1; j < _dim_gauss ; j++)
{
ss=0;
for(int k=0; k < j; k++)
{
ss+=pow(_L[j][k],2);
}
_L[j][j]=sqrt(cov_matrix[j][j]-ss);
for(int i=j+1; i < _dim_gauss; i++)
{
sp=0;
for(int k=0; k < j ; k++)
{
sp+=_L[i][k]*_L[j][k];
}
_L[i][j]=(cov_matrix[i][j]-sp)/_L[j][j];
}
}
//transform back to cov matrix
//cout << "L'="<< _L << endl;
for(int i=0; i < _dim_gauss ; i++)
{
for(int j=0; j < _dim_gauss ; j++)
{
_L[i][j]= _L[i][j]*sqrt(_scale)*sqrt(_cov_matrix[i][i]);
}
}
// cout << "L="<< _L << endl;
}
bool pdf_log_gemma::shoot(HepVector& jump, bool flag)
{
HepVector eta(_dim_gauss);
vector<double> eta_v(_dim_gauss);
_normal -> fill(eta_v);
if(flag==true)
{
Cholesky_decompose(_cov_matrix);
cout << "decompose=" << _L << endl;
}
for(int i=0; i < _dim_gauss ; i++)
{
eta[i]=eta_v[i];
}
jump=_L*eta;
return true;
}
void pdf_log_gemma::set_cov(const char* cov_file)
{
cout << "using:"<< cov_file << endl;
ifstream file_in;
file_in.open(cov_file);
if(!file_in)
{
cout <<"file not found";
}
for(int i=0; i<_dim_gauss; i++)
{
for(int j=0; j<_dim_gauss; j++)
{
file_in >> _cov_matrix[i][j];
}
}
file_in.close();
double det=_cov_matrix.determinant();
cout << "det= " << det << endl;
cout << "cov= " << _cov_matrix << endl;
int ifail=0;
_inv_cov=_cov_matrix;
_inv_cov.invert(ifail);
assert(ifail==0);
}
void pdf_log_gemma::calc_cov_matrix(HepSymMatrix& corr_matrix,vector<double>& cov_vector )
{
_cov_matrix=corr_matrix;
for(int i=0; i < _dim_gauss ; i++)
{
for(int j=0; j < _dim_gauss ; j++)
{
_cov_matrix[i][j]=corr_matrix[i][j]*(cov_vector[i]*cov_vector[j]);
}
}
int ifail=0;
_inv_cov=_cov_matrix;
_inv_cov.invert(ifail);
assert(ifail==0);
cout << "inverse=" << _inv_cov << endl;
}
void pdf_log_gemma::write_cov(const char* str)
{
cout << "writing covariance in file: "<< str << endl;
ofstream f(str);