CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
Forecast.cpp
Go to the documentation of this file.
1 /********************************************************************
2  * Copyright (C) 2010 by Federico Marulli and Alfonso Veropalumbo *
3  * federico.marulli3@unibo.it *
4  * *
5  * This program is free software; you can redistribute it and/or *
6  * modify it under the terms of the GNU General Public License as *
7  * published by the Free Software Foundation; either version 2 of *
8  * the License, or (at your option) any later version. *
9  * *
10  * This program is distributed in the hope that it will be useful, *
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13  * GNU General Public License for more details. *
14  * *
15  * You should have received a copy of the GNU General Public *
16  * License along with this program; if not, write to the Free *
17  * Software Foundation, Inc., *
18  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19  ********************************************************************/
20 
35 #include "Data1D.h"
36 #include "GlobalFunc.h"
37 
38 using namespace std;
39 
40 using namespace cbl;
41 
42 
43 // ============================================================================
44 
45 
46 std::vector<double> cbl::fit_covariance_matrix_2PCF_monopole (const std::vector<double> mean, const std::vector<std::vector<double>> mock_xi0, const bool doJK, const cbl::cosmology::Cosmology cosmology, const double nObjects, const double Volume, const double bias, const double redshift, const double rMin, const double rMax, const int nbins, const cbl::BinType bin_type, const std::string method_Pk, const double sigma_NL, const bool NL)
47 {
48  coutCBL << "Fitting the Covariance Matrix" << endl;
49  const int nmock = mock_xi0.size();
50 
51  double term = (doJK) ? nmock-1 : 1.;
52  vector<double> mock_mean(nbins, 0);
53 
54  if (mean.size()==0) {
55  vector<vector<double>> vv = transpose(mock_xi0);
56  for (int i=0; i<nbins; i++)
57  mock_mean[i] = Average(vv[i]);
58  }
59  else
60  for (int i=0; i<nbins; i++)
61  mock_mean[i] = mean[i];
62 
63  vector<vector<double>> diff(nmock, vector<double>(nbins, 0));
64  for (int i=0; i<nmock; i++)
65  for (int j=0; j<nbins; j++)
66  diff[i][j] = mock_xi0[i][j]-mock_mean[j];
67 
68  cosmology::Cosmology cosmo = cosmology;
69  vector<double> kk = logarithmic_bin_vector(200, 1.e-4, 1.e2);
70  vector<double> Pk0(kk.size()), Pk2(kk.size()), Pk4(kk.size());
71  double linear_growth_rate = cosmo.linear_growth_rate(redshift, 1);
72  double beta = linear_growth_rate/bias;
73 
74  double pk0_fact = bias*bias*(1+2.*beta/3+beta*beta/5);
75  double pk2_fact = 4./3*beta+4./7*beta*beta;
76  double pk4_fact = 8./35*beta+beta;
77 
78  if (sigma_NL==0 or NL==true) {
79  vector<double> pkDM = cosmo.Pk_matter(kk, method_Pk, NL, redshift);
80  for (size_t i=0; i<kk.size(); i++) {
81  Pk0[i] = pkDM[i]*pk0_fact;
82  Pk2[i] = pkDM[i]*pk2_fact;
83  Pk4[i] = pkDM[i]*pk4_fact;
84  }
85  }
86  else {
87  vector<double> PkDM = cosmo.Pk_matter_DeWiggled(method_Pk, "EisensteinHu", kk, redshift, sigma_NL);
88  for (size_t i=0; i<kk.size(); i++) {
89  Pk0[i] = PkDM[i]*pk0_fact;
90  Pk2[i] = PkDM[i]*pk2_fact;
91  Pk4[i] = PkDM[i]*pk4_fact;
92  }
93  }
94 
95  vector<double> nObj = linear_bin_vector(20, 0.2*nObjects, nObjects);
96  vector<double> Vol = linear_bin_vector(20, 0.2*Volume, Volume);
97 
98  auto func = [&] (vector<double> &parameters)
99  {
100  const double _nObjects = parameters[0];
101  const double _Volume = parameters[1];
102  if (_nObjects<0.2*nObjects || _nObjects>nObjects ) return 1.e30;
103  if (_Volume<0.2*Volume || _Volume>Volume ) return 1.e30;
104 
105  vector<double> rr;
106  vector<vector<double>> covariance, icovariance;
107  Covariance_XiMultipoles (rr, covariance, nbins, rMin, rMax, _nObjects, _Volume, kk, {Pk0, Pk2, Pk4}, {0}, bin_type);
108  for (int i=0; i<nbins; i++)
109  for (int j=0; j<nbins; j++)
110  covariance[i][j] /=term;
111 
112  invert_matrix(covariance, icovariance);
113  double chi2 = 0.;
114  for (int i=0; i<nmock; i++)
115  chi2+=v_M_vt(diff[i], icovariance);
116  double LL = nmock*nbins*nbins*log(2*par::pi)+nmock*log(determinant_matrix(covariance))+chi2;
117  return LL;
118  };
119 
120  double min_nObj = 0., min_Vol = 0.;
121  double max_nObj = 50*nObjects, max_Vol = 50*Volume;
122  double min_LL = 1.e30;
123  for (size_t i=0; i<nObj.size(); i++)
124  for (size_t j=0; j<Vol.size(); j++) {
125  vector<double> pars = {nObj[i], Vol[j]};
126  double LL = func(pars);
127  if (LL<min_LL) {
128  min_LL = LL;
129  min_nObj = nObj[i];
130  min_Vol = Vol[j];
131  }
132  }
133 
134  vector<double> start = {nObjects, Volume};
135  vector<vector<double>> limits = { {min_nObj, max_nObj}, {min_Vol, max_Vol}};
136  return wrapper::gsl::GSL_minimize_nD(func, start, limits);
137 }
138 
139 
140 // ============================================================================
141 
142 
143 std::shared_ptr<cbl::data::Data> cbl::generate_mock_2PCF_monopole (const cbl::cosmology::Cosmology cosmology, const double bias, const double nObjects, const double Volume, const double redshift, const double rMin, const double rMax, const int nbins, const cbl::BinType bin_type, const std::string method_Pk, const double sigma_NL, const bool NL)
144 {
145  auto data = generate_mock_2PCF_multipoles(cosmology, bias, nObjects, Volume, redshift, rMin, rMax, nbins, bin_type, method_Pk, sigma_NL, NL);
146 
147  int ndata = data->ndata()/3;
148 
149  vector<double> rad(ndata, 0), xi0(ndata, 0);
150  vector<vector<double>> covariance(ndata, vector<double>(ndata, 0));
151 
152  for (int i=0; i<ndata; i++) {
153  rad[i] = data->xx(i);
154  xi0[i] = data->data(i);
155  for (int j=0; j<ndata; j++)
156  covariance[i][j] = data->covariance(i, j);
157  }
158 
159  return make_shared<data::Data1D> (data::Data1D(rad, xi0, covariance));
160 }
161 
162 // ============================================================================
163 
164 
165 std::shared_ptr<cbl::data::Data> cbl::generate_mock_2PCF_multipoles (const cbl::cosmology::Cosmology cosmology, const double bias, const double nObjects, const double Volume, const double redshift, const double rMin, const double rMax, const int nbins, const cbl::BinType bin_type, const std::string method_Pk, const double sigma_NL, const bool NL)
166 {
167  cosmology::Cosmology cosmo = cosmology;
168  vector<double> kk = logarithmic_bin_vector(200, 1.e-4, 1.e2);
169  vector<double> Pk0(kk.size()), Pk2(kk.size()), Pk4(kk.size());
170  double linear_growth_rate = cosmo.linear_growth_rate(redshift, 1);
171  double beta = linear_growth_rate/bias;
172 
173  double pk0_fact = bias*bias*(1+2.*beta/3+beta*beta/5);
174  double pk2_fact = 4./3*beta+4./7*beta*beta;
175  double pk4_fact = 8./35*beta+beta;
176 
177  if (sigma_NL==0 or NL==true) {
178  vector<double> pkDM = cosmo.Pk_matter(kk, method_Pk, NL, redshift);
179  for (size_t i=0; i<kk.size(); i++) {
180  Pk0[i] = pkDM[i]*pk0_fact;
181  Pk2[i] = pkDM[i]*pk2_fact;
182  Pk4[i] = pkDM[i]*pk4_fact;
183  }
184  }
185  else {
186  vector<double> PkDM = cosmo.Pk_matter_DeWiggled(method_Pk, "EisensteinHu", kk, redshift, sigma_NL);
187  for (size_t i=0; i<kk.size(); i++) {
188  Pk0[i] = PkDM[i]*pk0_fact;
189  Pk2[i] = PkDM[i]*pk2_fact;
190  Pk4[i] = PkDM[i]*pk4_fact;
191  }
192  }
193 
194 
195  vector<double> rr;
196  vector<vector<double>> xil(3), covariance;
197 
198  Covariance_XiMultipoles (rr, covariance, nbins, rMin, rMax, nObjects, Volume, kk, {Pk0, Pk2, Pk4}, {0, 2, 4}, bin_type);
199 
200  xil[0] = wrapper::fftlog::transform_FFTlog(rr, 1, kk, Pk0, 0);
201  xil[1] = wrapper::fftlog::transform_FFTlog(rr, 1, kk, Pk2, 2);
202  xil[2] = wrapper::fftlog::transform_FFTlog(rr, 1, kk, Pk4, 4);
203 
204  vector<double> rad, xi;
205  for (int i=0; i<3; i++) {
206  for (int j=0; j<nbins; j++) {
207  rad.push_back(rr[j]);
208  xi.push_back(xil[i][j]);
209  }
210  }
211 
212  return make_shared<data::Data1D>(data::Data1D(rad, xi, covariance));
213 }
214 
The class Data1D.
Generic functions that use one or more classes of the CosmoBolognaLib.
#define coutCBL
CBL print message.
Definition: Kernel.h:734
The class Cosmology.
Definition: Cosmology.h:277
std::vector< double > Pk_matter_DeWiggled(const std::string linear_method, const std::string nowiggles_method, const std::vector< double > kk, const double redshift, const double sigma_NL, const int order=4, const int nknots=10, const double lambda=0.25, const bool store_output=true, const std::string output_root="test", const bool norm=1, const double prec=1.e-4)
the dark matter power spectrum, de-wiggled (see e.g. Anderson et al 2014)
Definition: PkXi.cpp:2193
std::vector< double > Pk_matter(const std::vector< double > kk, const std::string method_Pk, const bool NL, const double redshift, const bool store_output=true, const std::string output_root="test", const int norm=-1, const double k_min=0.001, const double k_max=100., const double prec=1.e-2, const std::string file_par=par::defaultString, const bool unit1=false)
the dark matter power spectrum
Definition: PkXi.cpp:1331
double linear_growth_rate(const double redshift, const double prec=1.e-4) const
the linear growth rate at a given redshift,
Definition: Cosmology.cpp:662
The class Data1D.
Definition: Data1D.h:51
static const double pi
: the ratio of a circle's circumference to its diameter
Definition: Constants.h:199
double bias(const double Mmin, const double sigmalgM, const double M0, const double M1, const double alpha, const std::shared_ptr< void > inputs)
the mean galaxy bias
std::vector< double > transform_FFTlog(const std::vector< double > yy, const int dir, const std::vector< double > xx, const std::vector< double > fx, const double mu=0, const double q=0, const double kr=1, const int kropt=0)
wrapper of the FFTlog to compute the discrete Fourier sine or cosine transform of a logarithmically...
Definition: FFTlog.cpp:43
std::vector< double > GSL_minimize_nD(FunctionDoubleVector func, const std::vector< double > start, const std::vector< std::vector< double >> ranges, const unsigned int max_iter=1000, const double tol=1.e-6, const double epsilon=0.1)
minimize the provided function using GSL procedure
Definition: GSLwrapper.cpp:497
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38
double Average(const std::vector< double > vect)
the average of a std::vector
Definition: Func.cpp:870
std::vector< T > logarithmic_bin_vector(const size_t nn, const T min, const T max)
fill a std::vector with logarithmically spaced values
Definition: Kernel.h:1621
std::vector< T > linear_bin_vector(const size_t nn, const T min, const T max)
fill a std::vector with linearly spaced values
Definition: Kernel.h:1604
void invert_matrix(const std::vector< std::vector< double >> mat, std::vector< std::vector< double >> &mat_inv, const double prec=1.e-10, const int Nres=-1)
function to invert a matrix
Definition: Func.cpp:673
T v_M_vt(const std::vector< T > vv, const std::vector< std::vector< T >> MM)
return the value of
Definition: Kernel.h:1833
std::shared_ptr< cbl::data::Data > generate_mock_2PCF_multipoles(const cbl::cosmology::Cosmology cosmology, const double bias, const double nObjects, const double Volume, const double redshift, const double rMin, const double rMax, const int nbins, const cbl::BinType bin_type, const std::string method_Pk="CAMB", const double sigma_NL=0., const bool NL=true)
generate mock measurementes of the 2PCF multipoles from gaussian covariance matrix
Definition: Forecast.cpp:165
std::vector< std::vector< T > > transpose(std::vector< std::vector< T >> matrix)
transpose a matrix
Definition: Kernel.h:1729
void Covariance_XiMultipoles(std::vector< double > &rr, std::vector< std::vector< double >> &covariance, const int nbins, const double rMin, const double rMax, const double nObjects, const double Volume, const std::vector< double > kk, const std::vector< std::vector< double >> Pk_multipoles, const std::vector< int > orders, const cbl::BinType bin_type=cbl::BinType::_linear_)
Covariance matrix for two-point correlation multipoles (see Grieb et al. 2016, Eq....
std::vector< double > fit_covariance_matrix_2PCF_monopole(const std::vector< double > mean, const std::vector< std::vector< double >> mock_xi0, const bool doJK, const cbl::cosmology::Cosmology cosmology, const double nObjects, const double Volume, const double bias, const double redshift, const double rMin, const double rMax, const int nbins, const cbl::BinType bin_type, const std::string method_Pk="CAMB", const double sigma_NL=0., const bool NL=true)
fit the input covariance matrix using the gaussian model, varying the number of objects and the volum...
Definition: Forecast.cpp:46
BinType
the binning type
Definition: Kernel.h:505
double determinant_matrix(const std::vector< std::vector< double >> mat)
compute the determinant of a matrix
Definition: Func.cpp:648
std::shared_ptr< cbl::data::Data > generate_mock_2PCF_monopole(const cbl::cosmology::Cosmology cosmology, const double bias, const double nObjects, const double Volume, const double redshift, const double rMin, const double rMax, const int nbins, const cbl::BinType bin_type, const std::string method_Pk="CAMB", const double sigma_NL=0., const bool NL=true)
generate mock measurementes of the 2PCF monopole from gaussian covariance matrix
Definition: Forecast.cpp:143