Commit 838e4498 authored by Plaszczynski Stephane's avatar Plaszczynski Stephane
Browse files

all data likelihoods now derive from Chi2Data

parent c0b2ad9e
......@@ -10,32 +10,23 @@
#include <iostream>
#include <cstring>
#include <string>
#include<fstream>
#include<cassert>
#include <libgen.h> //for basename
//#define SPEED_OF_LIGHT 299792.458
// km/s
// premiere implementation - mesure a un z unique
using namespace std;
//constructor
//should take in input BAO.par and fill the Matrix and vectors with it
BAO2D_chi2::BAO2D_chi2(const string paramfile, Engine* klass ) : Chi2Cosmo(basename(const_cast<char*>(paramfile.c_str()))),engine(klass)
BAO2D_chi2::BAO2D_chi2(const string& fileBAO, const Parameters& par, Engine* klass ):Chi2Data(par,klass)
{
Chi2Data::_name=basename(const_cast<char*>(fileBAO.c_str()));
//need it to fill _data, _cov_matrix, _z;
_is_normalized = 0; // check is still needed
_hinverted = 0;
read_par_file(paramfile);
}
//toy construcor (usefull?)
BAO2D_chi2::BAO2D_chi2(const string fileBAO)
{
//cout <<"I was here" << endl;
read_par_file(fileBAO);
}
......@@ -43,11 +34,10 @@ BAO2D_chi2::BAO2D_chi2(const string fileBAO)
//destructor
BAO2D_chi2::~BAO2D_chi2() {}
//functions
// compute actual likelihood (sais pas ce que c'est que tout ces const...)
double BAO2D_chi2::operator()(const std::vector<double>& par) const
double BAO2D_chi2::chi2_eff(const std::vector<double>& par) const
{
double chi2=0;
......@@ -79,11 +69,7 @@ double BAO2D_chi2::operator()(const std::vector<double>& par) const
cout << "B2D for "<< _Experiment_name <<" Chi2= "<<chi2<<endl;
}
Chi2Cosmo::registerChi2(chi2);
return chi2;
return chi2;
}
void BAO2D_chi2::read_par_file(const string fileName)
......
#ifndef BAO2D_chi2_hh
#define BAO2D_chi2_hh
#include "Chi2Cosmo.hh"
#include "Chi2Data.hh"
#include "Engine.hh"
#include "Parameters.hh"
#include "CLHEP/Matrix/Vector.h"
#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include <iostream>
#include<fstream>
#include<cassert>
#include<vector>
#include<string>
class BAO2D_chi2 : public Chi2Cosmo
class BAO2D_chi2 : public Chi2Data
{
public:
//constructors
//should take in input BAO.par and fill the Matrix and vectors with it
BAO2D_chi2(const std::string fileBAO, Engine * klass);
BAO2D_chi2(const std::string& fileBAO, const Parameters& par, Engine * klass);
//toy constructor
BAO2D_chi2(const std::string fileBAO);
......@@ -30,14 +30,13 @@ public:
// destructor
~BAO2D_chi2();
//to override: must return something comptible with -2log(L)
double chi2_eff(const std::vector<double>& par) const;
//no nuisances so use empty
std::vector<std::string> requires() const {return std::vector<std::string>();}
//operator to override: must return something comptible with -2log(L)
double operator()(const std::vector<double>& par) const;
void read_par_file(const std::string filename);
double Up() const {return 1.;} //what is this for??
//void setEngine(Engine* e){engine=e;} //to check
//maybe this are not useful...
//visualize parameters
HepSymMatrix Get_cov_matrix();
......
......@@ -7,17 +7,21 @@
#include <iostream>
#include <cstring>
#include <string>
#include <iostream>
#include <fstream>
#include <cassert>
#include <libgen.h> //for basename
//#define SPEED_OF_LIGHT 299792.458
using namespace std;
//constructor
//should take in input BAO.par and fill the Matrix and vectors with it
BAO_chi2::BAO_chi2(const string fileBAO, Engine* klass ):Chi2Cosmo(basename(const_cast<char*>(fileBAO.c_str()))),engine(klass)
BAO_chi2::BAO_chi2(const string& fileBAO, const Parameters& par, Engine* klass ):Chi2Data(par,klass)
{
Chi2Data::_name=basename(const_cast<char*>(fileBAO.c_str()));
//need it to fill _data, _cov_matrix, _z;
read_par_file(fileBAO);
......@@ -30,13 +34,6 @@ BAO_chi2::BAO_chi2(const string fileBAO, Engine* klass ):Chi2Cosmo(basename(cons
}
//toy construcor
BAO_chi2::BAO_chi2(const string fileBAO)
{
//cout <<"I was here" << endl;
read_par_file(fileBAO);
}
//destructor
BAO_chi2::~BAO_chi2() {}
......@@ -45,7 +42,7 @@ BAO_chi2::~BAO_chi2() {}
//functions
double BAO_chi2::operator()(const std::vector<double>& par) const
double BAO_chi2::chi2_eff(const std::vector<double>& par) const
{
......@@ -82,8 +79,6 @@ double BAO_chi2::operator()(const std::vector<double>& par) const
//this computes xT*_inv_cov*x
chi2=_inv_cov.similarity(diff);
//cout<< "chi2 bao="<< chi2 <<endl;
Chi2Cosmo::registerChi2(chi2);
return chi2;
}
......
......@@ -4,40 +4,37 @@
/////////////WARNING: THIS IS ACTAULLY A CHI2 (NOT LIKELIHOOD) /////////
#include "Chi2Cosmo.hh"
#include "Chi2Data.hh"
#include "Engine.hh"
#include "Parameters.hh"
#include "CLHEP/Matrix/Vector.h"
#include "CLHEP/Matrix/Matrix.h"
#include "CLHEP/Matrix/SymMatrix.h"
#include <iostream>
#include <fstream>
#include <cassert>
#include <vector>
#include <string>
class BAO_chi2 : public Chi2Cosmo
class BAO_chi2 : public Chi2Data
{
public:
//constructors
//should take in input BAO.par and fill the Matrix and vectors with it
BAO_chi2(const std::string fileBAO, Engine * klass);
BAO_chi2(const std::string& fileBAO, const Parameters& par, Engine * klass);
//toy constructor
BAO_chi2(const std::string fileBAO);
// destructor
~BAO_chi2();
//operator to override: must return something comptible with -2log(L)
double operator()(const std::vector<double>& par) const;
double chi2_eff(const std::vector<double>& par) const;
//no nuisances so used default requires
std::vector<std::string> requires() const {return std::vector<std::string>();}
void read_par_file(const std::string filename);
double Up() const {return 1.;} //what is this for??
//void setEngine(Engine* e){engine=e;} //to check
//maybe this are not useful...
//visualize parameters
......@@ -50,10 +47,7 @@ public:
void Set_data(HepVector a);
void Set_z(HepVector a);
//private:
Engine* engine; //when should class be called???
std::vector<std::string> _data_name;
HepVector _z;
HepVector _data;
......
......@@ -12,7 +12,6 @@
// Collaborating classes --
//-------------------------
#include"Engine.hh"
#include"Timer.hh"
#include "ClLikelihood.hh"
#include"cxxutils.h"
//--------------------
......@@ -35,15 +34,19 @@ using namespace std;
// Constructors --
//----------------
Chi2CMB::Chi2CMB(std::vector<ClLikelihood*>& mylik, Engine* klass,std::vector<string>& _parNames):Chi2Cosmo("CMB-all"),vlik(mylik),engine(klass),parNames(_parNames),iter(0),chi2_prev(100),timer(new Timer()),verbose(false),_lmax(-1)
Chi2CMB::Chi2CMB(std::vector<ClLikelihood*>& mylik, Engine* klass,const Parameters& p):Chi2Data(p,klass),vlik(mylik),verbose(false),_lmax(-1)
{
Chi2Data::_name="CMB_all";
//max des max
_lmax=vlik[0]->getTTmax();
for (size_t i=1;i<vlik.size();i++)
if (vlik[i]->getTTmax()>_lmax)
_lmax=vlik[i]->getTTmax();
for (size_t i=1;i<vlik.size();i++){
if (vlik[i]->getTTmax()>_lmax) _lmax=vlik[i]->getTTmax();
}
buildIndex();
}
Chi2CMB::Chi2CMB(std::vector<string>& _parNames):Chi2Cosmo("CMB-all"),parNames(_parNames),iter(0),chi2_prev(100),timer(new Timer()),verbose(false),_lmax(-1){
Chi2CMB::Chi2CMB(const Parameters& p):Chi2Data(p),verbose(false),_lmax(-1){
Chi2Data::_name="CMB_all";
}
......@@ -58,9 +61,6 @@ Chi2CMB::~Chi2CMB()
std::cout << "===== " << vlik[i]->name() << ": chi2=" << vlik[i]->chi2() << std::endl;
delete vlik[i];
}
//double tot=timer->total();
//cout <<"TIMER TOTAL TIME=" << tot << " s\t=" << tot/60. << " min\t=" << tot/3600 << "h" << endl;
delete timer;
}
//-----------------
......@@ -70,25 +70,30 @@ void
Chi2CMB::add(ClLikelihood* l) {
vlik.push_back(l);
if (l->getTTmax()>_lmax) _lmax= l->getTTmax();
buildIndex();
}
std::vector<std::string>
Chi2CMB::requires() const{
std::vector<std::string> nuiNames;
for (size_t i=0;i<vlik.size();i++){
std::copy (vlik[i]->nuisanceNames().begin(), vlik[i]->nuisanceNames().end(), std::back_inserter(nuiNames));
}
//cout << "requires called. vector is " << endl;
//copy(nuiNames.begin(),nuiNames.end(),ostream_iterator<string>(cout,"\n"));
return nuiNames;
}
double
Chi2CMB::operator()(const std::vector<double>& par) const {
Chi2CMB::chi2_eff(const std::vector<double>& par) const {
double chi2=0;
double chi2_prev=Chi2Data::_chi2;
//horrible trick pour changer la classe
Chi2CMB* self=const_cast<Chi2CMB*>(this);
ostringstream os;
if(verbose) {
os << iter << "\t";
copy(par.begin(),par.end(),ostream_iterator<double>(os,"\t"));
}
const int lmin=2;
vector<unsigned int> lvec(_lmax-lmin+1);
const int lmin=2;
vector<unsigned int> lvec(_lmax-lmin+1);
for (size_t i=0;i<lvec.size();i++){
lvec[i]=lmin+i;
}
......@@ -104,17 +109,9 @@ Chi2CMB::operator()(const std::vector<double>& par) const {
}
}
catch(exception &e){
cout << __FILE__ << ":\tfail getting Cls: " << e.what() << endl;
cerr << __FILE__ << ":\tfail getting Cls: " << e.what() << endl;
chi2=chi2_prev*1.2;
self->chi2_prev=chi2;
self->iter++;
if(verbose)
{
os << _FAILURE_ << "\t" << chi2 << "\t" << timer->partial();
cout << os.str() << endl;
}
return chi2;
}
//boucle vlikelihoods
......@@ -131,27 +128,19 @@ Chi2CMB::operator()(const std::vector<double>& par) const {
vector<string> nuiNames(lik->nuisanceNames());
for (int j=0;j<nuiNames.size();j++){
string n=nuiNames[j];
vector<string>::const_iterator it=find(parNames.begin(),parNames.end(),n);
planck_assert(it!=parNames.end(),"missing nuisance param="+n) ;
nuiPar.push_back(par[it-parNames.begin()]);
map<string,size_t>::const_iterator it = Chi2Data::index.find(n);
planck_assert(it!=index.end(),"Chi2CMB : parameter not found "+n);
nuiPar.push_back(par[it->second]);
}
//sait + si il faut couper lec cl...
dchi2=2*lik->computeLikelihood(lvec,cltt,clte,clee,clbb,nuiPar) ;
//cout << "chi2 " << i << " =" << dchi2 << endl;
chi2+=dchi2;
}
}
if(verbose){
os << fixed << setprecision(3) ;
os << _SUCCESS_ << "\t" << chi2 << "\t" << timer->partial();
cout << os.str() << endl;
}
self->iter++;
self->chi2_prev=chi2;
Chi2Cosmo::registerChi2(chi2);
return chi2;
}
......@@ -16,25 +16,25 @@
#ifndef Chi2CMB_hh
#define Chi2CMB_hh
#include "Chi2Cosmo.hh"
#include "Chi2Data.hh"
//forward declaration
class Engine;
class ClLikelihood;
class Timer;
#include<vector>
#include<string>
class Chi2CMB : public Chi2Cosmo
class Chi2CMB : public Chi2Data
{
public:
//constructors
Chi2CMB(std::vector<ClLikelihood*>& _mylik, Engine* _klass,std::vector<std::string>& _parNames);
//better : use dynamic methods to add likelihoods+engine
Chi2CMB(std::vector<std::string>& _parNames);
Chi2CMB(std::vector<ClLikelihood*>& _mylik, Engine* _klass,const Parameters& params);
//better : use the dynamic add method afterwards
Chi2CMB(const Parameters& params);
// destructor
~Chi2CMB();
......@@ -43,30 +43,26 @@ public:
void add(ClLikelihood* l) ;
void setEngine(Engine* e){engine=e;}
unsigned niter() const {return iter;};
void setVerbose(bool choice) { verbose=choice;}
//operator to override: must return something comptible with -2log(L)
double operator()(const std::vector<double>& par) const;
double chi2_eff(const std::vector<double>& par) const;
double Up() const {return 1.;}
//nuisances from concatenation of all likelihoods
std::vector<std::string> requires() const;
//the lmax of all likelihoods
int getTTmax() const {return _lmax;}
size_t nLik() const {return vlik.size();}
private:
Engine* engine;
std::vector<ClLikelihood*> vlik;
Timer* timer;
unsigned iter;
double chi2_prev;
bool verbose;
std::vector<std::string> parNames;
int _lmax;
//std::vector<std::string> nuiNames;
};
#endif
......
......@@ -64,7 +64,7 @@ ClikLikelihood::ClikLikelihood(char* file)
}
_name=string(basename(file));
dump();
//dump();
}
//--------------
// Destructor --
......
......@@ -28,8 +28,9 @@ using namespace std;
//---------------
// Constructors --
//----------------
Chi2Combiner::Chi2Combiner(Engine* e):Chi2Cosmo("Chi2Combiner"),iter(0),chi2_prev(100),timer(new Timer()),verbose(true),engine(e)
Chi2Combiner::Chi2Combiner(const Parameters& p,Engine* e):Chi2Data(p,e),iter(0),timer(new Timer()),verbose(true)
{
Chi2Data::_name="Chi2Combiner";
}
//--------------
// Destructor --
......@@ -49,11 +50,32 @@ Chi2Combiner::~Chi2Combiner()
//-----------------
// Member functions --
//-----------------
void Chi2Combiner::add(Chi2Data* c) {
obs.push_back(c);
buildIndex();
}
std::vector<std::string>
Chi2Combiner::requires() const{
std::vector<std::string> nuiNames;
for (size_t i=0;i<obs.size();i++){
vector<string> req=obs[i]->requires();
std::copy (req.begin(), req.end(), std::back_inserter(nuiNames));
}
return nuiNames;
}
double
Chi2Combiner::operator()(const std::vector<double>& par) const {
Chi2Combiner::chi2_eff(const std::vector<double>& par) const {
double chi2=0;
double chi2_prev=Chi2Data::_chi2;
//horrible trick pour changer la classe
Chi2Combiner* self=const_cast<Chi2Combiner*>(this);
......@@ -75,7 +97,6 @@ Chi2Combiner::operator()(const std::vector<double>& par) const {
if (!OK) {
//A TUNER!!!!!!!!!!!!!1
chi2=chi2_prev*10.;
self->chi2_prev=chi2;
self->iter++;
if(verbose)
{
......@@ -99,11 +120,7 @@ Chi2Combiner::operator()(const std::vector<double>& par) const {
os << "\t" << chi2 << "\t" << timer->partial();
cout << os.str() << endl;
}
self->iter++;
self->chi2_prev=chi2;
Chi2Cosmo::registerChi2(chi2);
return chi2;
}
......@@ -16,47 +16,47 @@
#ifndef Chi2Combiner_hh
#define Chi2Combiner_hh
#include "Chi2Cosmo.hh"
#include "Chi2Data.hh"
#include "Engine.hh"
#include "Parameters.hh"
#include<vector>
#include<string>
class Timer;
class Chi2Combiner : public Chi2Cosmo
class Chi2Combiner : public Chi2Data
{
public:
//constructors: takes ownership of engine : deletes it.
Chi2Combiner(Engine* e=0);
Chi2Combiner(const Parameters& p,Engine* e=0);
// destructor
~Chi2Combiner();
void add(Chi2Cosmo* c) {obs.push_back(c);}
void add(Chi2Data* c);
unsigned niter() const {return iter;};
void setVerbose(bool choice) { verbose=choice;}
//vector of cosmological parameters + nuisance : returns a Chi2
double operator()(const std::vector<double>& par) const;
double chi2_eff(const std::vector<double>& par) const;
//all nuisance names
virtual std::vector<std::string> requires() const;
//is the computation a success
bool success() const {return OK;}
//accessor to engine
Engine* getEngine() {return engine;}
//for minimization
double Up() const {return 1.;}
private:
Engine* engine;
std::vector<Chi2Cosmo*> obs;
std::vector<Chi2Data*> obs;
Timer* timer;
unsigned iter;
double chi2_prev;
bool verbose;
bool OK;
......
......@@ -16,19 +16,20 @@
// C++
//--------------------
#include<iostream>
#include<iomanip>
using namespace std;
using namespace ROOT::Minuit2;
//--------------------
// C
//----------------
#define OFMT(o,w) o.width(w); o<<left
//---------------
// Constructors --
//----------------
Chi2Data::Chi2Data(const Parameters& upar,Engine* e):user_par(upar),engine(e)
Chi2Data::Chi2Data(const Parameters& upar,Engine* e):user_par(upar),engine(e),_name("unknown")
{
updateIndex();
//updateIndex();
}
//--------------
// Destructor --
......@@ -44,16 +45,28 @@ Chi2Data::~Chi2Data()
double
Chi2Data::operator()(const std::vector<double>& x) const{
double val=chi2(x);
double val=chi2_eff(x);
//register (breaks const)
Chi2Data* me=const_cast<Chi2Data*>(this);
me->_chi2=val;
return val;
}
void
Chi2Data::buildIndex(){
vector<string> nui=requires();
//check nuisances present
for (size_t i=0;i<nui.size