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)
48 coutCBL <<
"Fitting the Covariance Matrix" << endl;
49 const int nmock = mock_xi0.size();
51 double term = (doJK) ? nmock-1 : 1.;
52 vector<double> mock_mean(nbins, 0);
55 vector<vector<double>> vv =
transpose(mock_xi0);
56 for (
int i=0; i<nbins; i++)
60 for (
int i=0; i<nbins; i++)
61 mock_mean[i] = mean[i];
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];
70 vector<double> Pk0(kk.size()), Pk2(kk.size()), Pk4(kk.size());
72 double beta = linear_growth_rate/
bias;
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;
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;
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;
98 auto func = [&] (vector<double> ¶meters)
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;
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;
114 for (
int i=0; i<nmock; i++)
115 chi2+=
v_M_vt(diff[i], icovariance);
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);
134 vector<double> start = {nObjects, Volume};
135 vector<vector<double>> limits = { {min_nObj, max_nObj}, {min_Vol, max_Vol}};
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)
145 auto data =
generate_mock_2PCF_multipoles(cosmology,
bias, nObjects, Volume, redshift, rMin, rMax, nbins, bin_type, method_Pk, sigma_NL, NL);
147 int ndata = data->ndata()/3;
149 vector<double> rad(ndata, 0), xi0(ndata, 0);
150 vector<vector<double>> covariance(ndata, vector<double>(ndata, 0));
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);
159 return make_shared<data::Data1D> (
data::Data1D(rad, xi0, covariance));
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)
169 vector<double> Pk0(kk.size()), Pk2(kk.size()), Pk4(kk.size());
171 double beta = linear_growth_rate/
bias;
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;
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;
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;
196 vector<vector<double>> xil(3), covariance;
198 Covariance_XiMultipoles (rr, covariance, nbins, rMin, rMax, nObjects, Volume, kk, {Pk0, Pk2, Pk4}, {0, 2, 4}, bin_type);
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]);
212 return make_shared<data::Data1D>(
data::Data1D(rad, xi, covariance));
Generic functions that use one or more classes of the CosmoBolognaLib.
#define coutCBL
CBL print message.
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)
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
double linear_growth_rate(const double redshift, const double prec=1.e-4) const
the linear growth rate at a given redshift,
static const double pi
: the ratio of a circle's circumference to its diameter
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...
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
The global namespace of the CosmoBolognaLib
double Average(const std::vector< double > vect)
the average of a std::vector
std::vector< T > logarithmic_bin_vector(const size_t nn, const T min, const T max)
fill a std::vector with logarithmically spaced values
std::vector< T > linear_bin_vector(const size_t nn, const T min, const T max)
fill a std::vector with linearly spaced values
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
T v_M_vt(const std::vector< T > vv, const std::vector< std::vector< T >> MM)
return the value of
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
std::vector< std::vector< T > > transpose(std::vector< std::vector< T >> matrix)
transpose a matrix
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...
double determinant_matrix(const std::vector< std::vector< double >> mat)
compute the determinant of a matrix
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