Commit 0a1ca84a authored by Xavier Garrido's avatar Xavier Garrido
Browse files

[WIP] new hillipop with CIB model inside

parent a7f46ffb
......@@ -25,7 +25,7 @@
#include<iterator>
#include<iomanip>
//--------------------
// C
// C
//----------------
#include"common.h" //de CLASS
......@@ -56,17 +56,17 @@ Chi2CMB::~Chi2CMB()
//-----------------
// Member functions --
//-----------------
void
void
Chi2CMB::add(ClLikelihood* l) {
vlik.push_back(l);
vlik.push_back(l);
if (l->getTTmax()>_lmax) _lmax= l->getTTmax();
buildIndex();
}
std::vector<std::string>
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));
......@@ -76,12 +76,12 @@ Chi2CMB::requires() const{
return nuiNames;
}
double
double
Chi2CMB::chi2_eff(const std::vector<double>& par) const {
double chi2=0;
double chi2_prev=Chi2Data::_chi2;
const int lmin=2;
vector<unsigned int> lvec(_lmax-lmin+1);
for (size_t i=0;i<lvec.size();i++){
......@@ -93,7 +93,7 @@ Chi2CMB::chi2_eff(const std::vector<double>& par) const {
vector<double> clbb;
try{
#pragma omp critical
#pragma omp critical
{
engine->getCls(lvec,cltt,clte,clee,clbb);
}
......@@ -106,11 +106,11 @@ Chi2CMB::chi2_eff(const std::vector<double>& par) const {
//boucle vlikelihoods
double dchi2(0);
#pragma omp parallel
#pragma omp parallel
{
int i;
#pragma omp for schedule(static)
#pragma omp for schedule(static)
for (i=0;i<vlik.size();i++){
ClLikelihood* lik=vlik[i];
//need to build nuisance list
......@@ -122,7 +122,7 @@ Chi2CMB::chi2_eff(const std::vector<double>& par) const {
}
//sait + si il faut couper lec cl...
dchi2=2*lik->computeLikelihood(lvec,cltt,clte,clee,clbb,nuiPar) ;
dchi2=2*lik->computeLikelihood(lvec,cltt,clte,clee,clbb,nuiPar,engine) ;
//cout << "chi2 " << i << " =" << dchi2 << endl;
chi2+=dchi2;
......
......@@ -9,7 +9,7 @@
// Stephane Plaszczynski (plaszczy@lal.in2p3.fr)
//
// History (add to end):
// creation: Fri Feb 17 11:09:55 CET 2012
// creation: Fri Feb 17 11:09:55 CET 2012
//
//------------------------------------------------------------------------
......@@ -23,6 +23,8 @@ using std::vector;
using std::string;
#include"cxxutils.h"
class Engine;
class ClLikelihood
{
......@@ -35,7 +37,7 @@ public:
//TT lmax returned (maybe not a good idea)
virtual int getTTmax() const=0;
//_n should contain the name of each nuisance par
//_n should contain the name of each nuisance par
inline virtual const vector<std::string>& nuisanceNames() const{ return _n;};
......@@ -43,23 +45,34 @@ public:
// ouptut = - log(L)
//most general form
virtual double computeLikelihood(const std::vector<unsigned int>& l,
std::vector<double>& cltt,
std::vector<double>& clte,
std::vector<double>& clee,
virtual double computeLikelihood(const std::vector<unsigned int>& l,
std::vector<double>& cltt,
std::vector<double>& clte,
std::vector<double>& clee,
std::vector<double>& clbb,
std::vector<double>& nuisance)=0;
//if no nuisance calls previous with empty list
virtual double computeLikelihood(const std::vector<unsigned int>& l,
std::vector<double>& cltt,
std::vector<double>& clte,
std::vector<double>& clee,
virtual double computeLikelihood(const std::vector<unsigned int>& l,
std::vector<double>& cltt,
std::vector<double>& clte,
std::vector<double>& clee,
std::vector<double>& clbb){
std::vector<double> nui;
return computeLikelihood(l,cltt,clte,clee,clbb,nui);
}
//passing Boltzman solver
virtual double computeLikelihood(const std::vector<unsigned int>& l,
std::vector<double>& cltt,
std::vector<double>& clte,
std::vector<double>& clee,
std::vector<double>& clbb,
std::vector<double>& nuisance,
Engine * engine) {
return computeLikelihood(l,cltt,clte,clee,clbb,nuisance);
}
std::string name() const {return _name;} //lkl name
double chi2() const {return 2*_lnlkl;};
......@@ -82,4 +95,3 @@ protected:
};
#endif
......@@ -9,7 +9,7 @@
// Benjamin Rouille d'Orfeuil (rouille@lal.in2p3.fr)
//
// History (add to end):
// creation: Tue Oct 8 14:40:00 CET 2013
// creation: Tue Oct 8 14:40:00 CET 2013
// test for DirNames in case it is relative : 31/08/2015
//
//------------------------------------------------------------------------
......@@ -35,13 +35,13 @@ public:
// Get maximum multipole for a given mode
inline int getLmax(cltype t) const {return _maxOflmax[t];}
// Get TT maximum multipole. This will be the one use by CLASS.
int getTTmax() const {return getLmax(TT);}
// Tells you which mode are considered
bool hasCl(cltype t) const {return _typeStatus[t] == 1;}
// Returns what is currently used for the calculation.
void dump() const;
......@@ -50,7 +50,7 @@ public:
// Overidden
double computeLikelihood(const vector<unsigned int>& l,std::vector<double>& cltt,std::vector<double>& clte,std::vector<double>& clee,std::vector<double>& clbb,std::vector<double>& nuisance);
// Use base definition
double computeLikelihood(const std::vector<unsigned int>& l,std::vector<double>& cltt,std::vector<double>& clte,std::vector<double>& clee,std::vector<double>& clbb) {
std::vector<double> nuisance;
......@@ -58,20 +58,20 @@ public:
}
// destructor
~Hillipop();
virtual ~Hillipop();
private:
protected:
// Process the parameters file
void Init(const std::string fileName);
// Index - cross-spectra matching
void ProduceList();
// Reads the file that enclosed the multipole extrema of each cross-spectrum
// Reads the file that enclosed the multipole extrema of each cross-spectrum
// for each mode
void ProcessMultipoles(const std::string fileName);
// Reads the file that enclosed the eigen values of the beams
void ProcessBeamsEigenmodes(const std::string pathToBeams);
......@@ -81,7 +81,7 @@ private:
// Reads the cross-spectra
void ProcessXSpectraErrors(const std::string pathToXSpectraErrors);
// Reads the covariance matrix. Only correlation between cross-spectra are
// Reads the covariance matrix. Only correlation between cross-spectra are
// assumed for the moment
void ProcessCovMatrix(const std::string pathToMatrix);
......@@ -185,8 +185,7 @@ private:
// Residuals
vector<vector<vector<double> > > _residual, _weight, _instrumental;
};
#endif
#include "CMB/HillipopCIB.hh"
#include "Class/MnClassEngine.hh"
HillipopCIB::HillipopCIB(const string fileName) : Hillipop(fileName)
{
std::cout << "HillipopCIB::HillipopCIB: Entering..." << std::endl;
}
HillipopCIB::~HillipopCIB()
{
}
double HillipopCIB::computeLikelihood(const vector<unsigned int>& l,
vector<double>& cltt,
vector<double>& clte,
vector<double>& clee,
vector<double>& clbb,
vector<double>& nuisance,
Engine * engine)
{
std::cout << "HillipopCIB::computeLikelihood: Entering..." << std::endl;
// Reshaping the theoretical CMB power spectra
vector<double> ClCMB(5*(_maxOflmax[0]+1),0);
for(unsigned int i = 0; i <= _maxOflmax[0]-l[0]; i++) {
ClCMB[0*(_maxOflmax[0]+1)+l[i]] = cltt[i];
ClCMB[1*(_maxOflmax[0]+1)+l[i]] = clee[i];
ClCMB[2*(_maxOflmax[0]+1)+l[i]] = clbb[i];
ClCMB[3*(_maxOflmax[0]+1)+l[i]] = clte[i];
ClCMB[4*(_maxOflmax[0]+1)+l[i]] = clte[i];
}
//Compute Residuals
computeResiduals(ClCMB, nuisance, false, engine);
double chi2 = 0;
double ndof = 0;
// Begin: chi2
// diagonal
for(unsigned int m = 0; m < 4; m++) {
if( _typeStatus[m] ) {
for(unsigned int f1 = 0; f1 < _nXFreq; f1++) {
if( _XFreqStatus[m][f1] ) {
for(unsigned int l1 = _lminXFreq[m][f1]; l1 <= _lmaxXFreq[m][f1]; l1++) {
ndof++;
for(unsigned int l2 = _lminXFreq[m][f1]; l2 <= _lmaxXFreq[m][f1]; l2++) {
chi2 += _residual[m][f1][l1]*_invCovMat[m][m][f1][f1][l1-_lminXFreq[m][f1]][l2-_lminXFreq[m][f1]]*_residual[m][f1][l2];
}
}
for(unsigned int f2 = f1+1; f2 < _nXFreq; f2++) {
if( _XFreqStatus[m][f2] ) {
for(unsigned int l1 = _lminXFreq[m][f1]; l1 <= _lmaxXFreq[m][f1]; l1++) {
for(unsigned int l2 = _lminXFreq[m][f2]; l2 <= _lmaxXFreq[m][f2]; l2++ ) {
chi2 += 2*_residual[m][f1][l1]*_invCovMat[m][m][f1][f2][l1-_lminXFreq[m][f1]][l2-_lminXFreq[m][f2]]*_residual[m][f2][l2];
}
}
}
}
}
}
}
}
// off-diagonal
for(unsigned int m1 = 0; m1 < 3; m1++) {
if( _typeStatus[m1] ) {
for(unsigned int m2 = m1+1; m2 < 4; m2++) {
if( _typeStatus[m2] ) {
for(unsigned int f1 = 0; f1 < _nXFreq; f1++) {
if( _XFreqStatus[m1][f1] ) {
for(unsigned int f2 = 0; f2 < _nXFreq; f2++) {
if( _XFreqStatus[m2][f2] ) {
for(unsigned int l1 = _lminXFreq[m1][f1]; l1 <= _lmaxXFreq[m1][f1]; l1++) {
for(unsigned int l2 = _lminXFreq[m2][f2]; l2 <= _lmaxXFreq[m2][f2]; l2++) {
chi2 += 2*_residual[m1][f1][l1]*_invCovMat[m1][m2][f1][f2][l1-_lminXFreq[m1][f1]][l2-_lminXFreq[m2][f2]]*_residual[m2][f2][l2];
}
}
}
}
}
}
}
}
}
}
_residual.clear();
_instrumental.clear();
_weight.clear();
// End: chi2
_chi2 = chi2;
_ndof = ndof;
// if( _fileOut != "" ) {
// computeResiduals( ClCMB, nuisance, true);
// WriteOutput(ClCMB,nuisance);
// _residual.clear();
// _instrumental.clear();
// _weight.clear();
// }
_lnlkl= _chi2/2.;
return _chi2/2.;
}
void HillipopCIB::computeResiduals(vector<double> ClCMB, vector<double>& nuisance, bool fullRange, Engine * engine)
{
std::cout << "HillipopCIB::computeResiduals: Entering..." << std::endl;
MnClassEngine * klass = dynamic_cast<MnClassEngine*>(engine);
if (!klass) {
throw std::logic_error("HillipopCIB::computeResiduals: Casting to ClassEngine fails!");
}
std::cout << "HillipopCIB::computeResiduals: Tcmb = " << klass->Tcmb() << std::endl;
unsigned int lminXFreq, lmaxXFreq, lminXSpectra, lmaxXSpectra = 0;
double tmpBeamEigenmodes = 1.;
double foregrounds = 0;
unsigned int modeOffset = 0;
double ApsRadio = 0.0, ApsDusty = 0.0, Asz = 0.0, Acib = 0.0, Aksz = 0.0, Aszxcib = 0.0;
double AdustTT = 0.0, AdustTP = 0.0, AdustPP = 0.0;
_residual.resize(5);
_instrumental.resize(5);
_weight.resize(5);
// Begin: gathering information regarding the parameters
// Calibration
vector<double> cal(_nMap);
for(unsigned int i = 0; i < _nMap; i++) cal[i] = nuisance[i];
//Absolute Calibration
const double Aplanck=nuisance[getIndex("A_planck")];
#ifdef BeamEigenmodes
// First beam eigenmode amplitude
vector<double> beta(_nXSpectra);
for(unsigned int c = 0; c < _nXSpectra; c++) beta[c] = nuisance[_nMap+c];
#endif
// Amplitudes of the point sources
if( _typeStatus[0]) ApsRadio = nuisance[getIndex("Aradio")];
if( _typeStatus[0]) ApsDusty = nuisance[getIndex("Adusty")];
// tSZ amplitude
if( _typeStatus[0]) Asz = nuisance[getIndex("Asz")];
// CIB amplitude
if( _typeStatus[0]) Acib = nuisance[getIndex("Acib")];
// Dust amplitude
if( _typeStatus[0]) AdustTT = nuisance[getIndex("AdustTT")];
if( _typeStatus[1] || _typeStatus[2]) AdustPP = nuisance[getIndex("AdustPP")];
if( _typeStatus[3]) AdustTP = nuisance[getIndex("AdustTP")];
// kSZ amplitude
if( _typeStatus[0]) Aksz = nuisance[getIndex("Aksz")];
// tSZxCIB amplitude
if( _typeStatus[0]) Aszxcib = nuisance[getIndex("Aszxcib")];
#ifdef PolarizationEfficiency
// Polarization efficiency
vector<double> epsilon(_nMap);
for(unsigned int i = 0; i < _nMap; i++) epsilon[i] = nuisance[getIndex("epsilon"+str(i))];
#endif
// End: gathering information regarding the parameters
// Init
for(unsigned int m = 0; m < 5; m++) {
_residual[m].resize(_nXFreq); _weight[m].resize(_nXFreq); _instrumental[m].resize(_nXFreq);
for(unsigned int f = 0; f < _nXFreq; f++) {
_residual[m][f].resize(_maxOflmax[0]+1,0);
_weight[m][f].resize(_maxOflmax[0]+1,0);
_instrumental[m][f].resize(_maxOflmax[0]+1,0);
}
}
lminXFreq = lminXSpectra = 10;
lmaxXFreq = lmaxXSpectra = 2500;
// lminXFreq = 50;
// lmaxXFreq = 2500;
// lminXSpectra = 50;
// lmaxXSpectra = 2500;
// Begin: Residuals
double tmpCal, tmpWeight;
#ifdef PolarizationEfficiency
double tmpPolEffPP, tmpPolEffTP, tmpPolEffPT = 0.;
#endif
for(unsigned int m = 0; m < 5; m++) {
if( _modeStatus[m] ) {
modeOffset = m*(_maxOflmax[0]+1);
for(unsigned int f = 0; f < _nXFreq; f++) {
if( _XFreqStatus[m][f] ) {
if( !fullRange) {
lminXFreq = _lminXFreq[m][f];
lmaxXFreq = _lmaxXFreq[m][f];
}
for(unsigned int l = lminXFreq; l <= lmaxXFreq; l++) {
if( m == 0 ) foregrounds = ApsRadio*_ClpsRadio[f]+ApsDusty*_ClpsDusty[f]+Asz*_ClSZ[f][l]+Acib*_ClCIB[f][l]+AdustTT*_ClDust[f][modeOffset+l]+Aksz*_ClkSZ[l]+Aszxcib*_ClSZxCIB[f][l];
if( m == 1 || m == 2 ) foregrounds = AdustPP*_ClDust[f][modeOffset+l];
if( m == 3 || m == 4 ) foregrounds = AdustTP*_ClDust[f][modeOffset+l];
for(unsigned int c = 0; c < _XFreq2XSpectra[f].size(); c++) {
if( !fullRange) {
lminXSpectra = _lminXSpectra[m][_XFreq2XSpectra[f][c]];
lmaxXSpectra = _lmaxXSpectra[m][_XFreq2XSpectra[f][c]];
}
if( _XSpectraStatus[m][_XFreq2XSpectra[f][c]] && l >= lminXSpectra && l <= lmaxXSpectra ) {
tmpWeight = _ClWeightData[_XFreq2XSpectra[f][c]][modeOffset+l];
#ifdef BeamEigenmodes
tmpBeamEigenmodes = pow(1+beta[_XFreq2XSpectra[f][c]]*_beamEigenmodes[_XFreq2XSpectra[f][c]][l],2);
#endif
tmpCal = (Aplanck*Aplanck) * (1+cal[_XSpectra2Maps[_XFreq2XSpectra[f][c]][0]]+cal[_XSpectra2Maps[_XFreq2XSpectra[f][c]][1]]);
#ifdef PolarizationEfficiency
if( m == 0 ) {
_residual[m][f][l] += tmpWeight*(_ClData[_XFreq2XSpectra[f][c]][modeOffset+l]-tmpCal*tmpBeamEigenmodes*(ClCMB[modeOffset+l]+foregrounds));
_instrumental[m][f][l] += tmpWeight*tmpCal*tmpBeamEigenmodes;
}
if( m == 1 || m == 2 ) {
tmpPolEffPP = (1+epsilon[_XSpectra2Maps[_XFreq2XSpectra[f][c]][0]])*(1+epsilon[_XSpectra2Maps[_XFreq2XSpectra[f][c]][1]]);
_residual[m][f][l] += tmpWeight*(_ClData[_XFreq2XSpectra[f][c]][modeOffset+l]-tmpCal*tmpPolEffPP*tmpBeamEigenmodes*(ClCMB[modeOffset+l]+foregrounds));
_instrumental[m][f][l] += tmpWeight*tmpCal*tmpPolEffPP*tmpBeamEigenmodes;
}
if( m == 3 ) {
tmpPolEffTP = 1+epsilon[_XSpectra2Maps[_XFreq2XSpectra[f][c]][1]];
_residual[m][f][l] += tmpWeight*(_ClData[_XFreq2XSpectra[f][c]][modeOffset+l]-tmpCal*tmpPolEffTP*tmpBeamEigenmodes*(ClCMB[modeOffset+l]+foregrounds));
_instrumental[m][f][l] += tmpWeight*tmpCal*tmpPolEffTP*tmpBeamEigenmodes;
}
if( m == 4 ) {
tmpPolEffPT = 1+epsilon[_XSpectra2Maps[_XFreq2XSpectra[f][c]][0]];
_residual[m][f][l] += tmpWeight*(_ClData[_XFreq2XSpectra[f][c]][modeOffset+l]-tmpCal*tmpPolEffPT*tmpBeamEigenmodes*(ClCMB[modeOffset+l]+foregrounds));
_instrumental[m][f][l] += tmpWeight*tmpCal*tmpPolEffPT*tmpBeamEigenmodes;
}
#else
_residual[m][f][l] += tmpWeight*(_ClData[_XFreq2XSpectra[f][c]][modeOffset+l]-tmpCal*tmpBeamEigenmodes*(ClCMB[modeOffset+l]+foregrounds));
_instrumental[m][f][l] += tmpWeight*tmpCal*tmpBeamEigenmodes;
#endif
_weight[m][f][l] += tmpWeight;
}
}
}
}
}
}
}
for(unsigned int m = 0; m < 4; m++) {
if( _typeStatus[m] ) {
for(unsigned int f = 0; f < _nXFreq; f++) {
if( _XFreqStatus[m][f] ) {
if( !fullRange) {
lminXFreq = _lminXFreq[m][f];
lmaxXFreq = _lmaxXFreq[m][f];
}
for(unsigned int l = lminXFreq; l <= lmaxXFreq; l++) {
if( m < 3) {
_instrumental[m][f][l] /= _weight[m][f][l];
_residual[m][f][l] /= _weight[m][f][l];
}
else {
_instrumental[3][f][l] = (_instrumental[3][f][l]+_instrumental[4][f][l])/(_weight[3][f][l]+_weight[4][f][l]);
_residual[3][f][l] = (_residual[3][f][l]+_residual[4][f][l])/(_weight[3][f][l]+_weight[4][f][l]);
}
}
}
}
}
}
_residual.resize(4);
_instrumental.resize(4);
// End: Residuals
}
#ifndef HillipopCIB_hh
#define HillipopCIB_hh
#include "CMB/Hillipop.hh"
class Engine;
class HillipopCIB : public Hillipop
{
public:
// Constructor
HillipopCIB(const std::string fileName);
// destructor
~HillipopCIB();
// Overidden
double computeLikelihood(const vector<unsigned int>& l,
std::vector<double>& cltt,
std::vector<double>& clte,
std::vector<double>& clee,
std::vector<double>& clbb,
std::vector<double>& nuisance,
Engine * engine);
// Compute residuals
void computeResiduals(vector<double> ClCMB, vector<double>& nuisance, bool fullRange, Engine * engine);
};
#endif
......@@ -30,6 +30,7 @@
#include "CMB/ClikLikelihood.hh"
#include "CMB/HiLLiPOP.hh"
#include "CMB/Hillipop.hh"
#include "CMB/HillipopCIB.hh"
#include "CMB/Lollipop.hh"
#include "CMB/HighEll.hh"
#include "CMB/SPT_high.hh"
......@@ -75,22 +76,22 @@ using namespace std;
//----------------
Chi2Combiner*
Chi2Combiner*
Chi2Factory::gimeChi2(Parser& parser){
//chi2 construit a partir de tous le Cl-Likeilihood
Chi2CMB* allCl= new Chi2CMB(parser.vars());
//CLIK case
#ifdef CLIK
#ifdef CLIK
if (parser.params.param_present("clikfile") || parser.params.param_present("clikfile1")){
const string clikfile=( parser.params.param_present("clikfile")? Parser::CheckPath(parser.params.find<string>("clikfile","")) : Parser::CheckPath(parser.params.find<string>("clikfile1","")) );
planck_assert(file_present(clikfile),"mssing CLIK file");
ClLikelihood* lik=new ClikLikelihood(const_cast<char*>(clikfile.c_str()));
allCl->add(lik);
//autres clik?
int iclik=2;
while(parser.params.param_present(string("clikfile")+dataToString(iclik))){
......@@ -126,13 +127,19 @@ Chi2Factory::gimeChi2(Parser& parser){
allCl->add(lal);
}
// et les high ell
//Hillipop (with PS and CIB model)
if (parser.params.param_present("HillipopCIB")){
const string file=Parser::CheckPath(parser.params.find<string>("HillipopCIB",""));
allCl->add(new HillipopCIB(file));
}
// et les high ell
// SPT low
if (parser.params.param_present("SPT_Low")){
const string file=Parser::CheckPath(parser.params.find<string>("SPT_Low",""));
SPTLow * sptl = new SPTLow(file);
allCl->add(sptl);
}
}
// SPThigh
if (parser.params.param_present("SPT_High")){
const string file=Parser::CheckPath(parser.params.find<string>("SPT_High",""));
......@@ -176,7 +183,7 @@ Chi2Factory::gimeChi2(Parser& parser){
if (!fid && lensing->lmax() > lmax) lmax=lensing->lmax();
}