CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
TwoPointCorrelation_projected.cpp
Go to the documentation of this file.
1 /********************************************************************
2  * Copyright (C) 2015 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 
39 #include "Data1D_extra.h"
40 
41 using namespace std;
42 
43 using namespace cbl;
44 using namespace catalogue;
45 using namespace chainmesh;
46 using namespace pairs;
47 using namespace measure;
48 using namespace twopt;
49 
50 
51 // ============================================================================
52 
53 
54 std::shared_ptr<data::Data> cbl::measure::twopt::TwoPointCorrelation_projected::data_with_extra_info (const std::vector<double> rp, const std::vector<double> ww, const std::vector<double> error) const
55 {
57 
58  vector<double> weightTOT(dd2D->nbins_D1(), 0.), scale_mean(dd2D->nbins_D1(), 0.), scale_sigma(dd2D->nbins_D1(), 0.), z_mean(dd2D->nbins_D1(), 0.), z_sigma(dd2D->nbins_D1(), 0.);
59 
60  double fact_err, fact_scale, fact_z;
61 
62  for (int i=0; i<dd2D->nbins_D1(); ++i) {
63 
64  for (int j=0; j<dd2D->nbins_D2(); ++j)
65  weightTOT[i] += dd2D->PP2D_weighted(i, j);
66 
67  for (int j=0; j<dd2D->nbins_D2(); ++j) {
68  scale_mean[i] += dd2D->scale_D1_mean(i, j)*dd2D->PP2D_weighted(i, j)/weightTOT[i];
69  z_mean[i] += dd2D->z_mean(i, j)*dd2D->PP2D_weighted(i, j)/weightTOT[i];
70  }
71 
72  scale_sigma[i] = pow(dd2D->scale_D1_sigma(i, 0), 2)*dd2D->PP2D_weighted(i, 0);
73  z_sigma[i] = pow(dd2D->z_sigma(i, 0), 2)*dd2D->PP2D_weighted(i, 0);
74 
75  for (int j=1; j<dd2D->nbins_D2(); ++j) {
76  if (dd2D->PP2D_weighted(i, j)>0) {
77  fact_err = dd2D->PP2D_weighted(i, j)*dd2D->PP2D_weighted(i, j-1)/(dd2D->PP2D_weighted(i, j)+dd2D->PP2D_weighted(i, j-1));
78  fact_scale = pow(dd2D->scale_D1_mean(i, j)-dd2D->scale_D1_mean(i, j-1), 2)*fact_err;
79  fact_z = pow(dd2D->z_mean(i, j)-dd2D->z_mean(i, j-1), 2)*fact_err;
80  scale_sigma[i] += pow(dd2D->scale_D1_sigma(i, j), 2)*dd2D->PP2D_weighted(i, j)+fact_scale;
81  z_sigma[i] += pow(dd2D->z_sigma(i, j),2)*weightTOT[i]+fact_z;
82  }
83  }
84  }
85 
86  vector<vector<double>> extra(4);
87 
88  for (int i=0; i<dd2D->nbins_D1(); ++i) {
89  extra[0].push_back(scale_mean[i]);
90  extra[1].push_back(sqrt(scale_sigma[i]/weightTOT[i]));
91  extra[2].push_back(z_mean[i]);
92  extra[3].push_back(sqrt(z_sigma[i]/weightTOT[i]));
93  }
94 
95  return move(unique_ptr<data::Data1D_extra>(new data::Data1D_extra(rp, ww, error, extra)));
96 }
97 
98 
99 // ============================================================================================
100 
101 
102 std::shared_ptr<data::Data> cbl::measure::twopt::TwoPointCorrelation_projected::Projected (const std::vector<double> rp, const std::vector<double> pi, const std::vector<std::vector<double>> xi, const std::vector<std::vector<double>> error_xi)
103 {
104  vector<double> ww, error;
105 
106  ww.resize(0); ww.resize(rp.size(), 0.);
107  error.resize(0); error.resize(rp.size(), 0.);
108 
109  const double binSize = 1./m_dd->binSize_inv_D2();
110 
111  const int pim = nint((m_piMax_integral-Min(pi))/binSize); // to convert from Mpc/h into the vector index
112 
113  for (size_t i=0; i<rp.size(); i++) {
114  ww[i] = 0.;
115  error[i] = 0.;
116 
117  for (size_t j=0; j<pim*1.000001; j++) {
118  if (j<xi[i].size() and j<error_xi[i].size()) {
119  ww[i] = ww[i]+2.*binSize*xi[i][j];
120  if (ww[i]>-1.) error[i] += pow(2.*binSize*error_xi[i][j], 2); // check!!!!
121  }
122  }
123  }
124 
125  for_each( error.begin(), error.end(), [] (double &vv) { vv = sqrt(vv);} );
126 
127  return (!m_compute_extra_info) ? move(unique_ptr<data::Data1D>(new data::Data1D(rp, ww, error))) : data_with_extra_info(rp, ww, error);
128 
129 }
130 
131 
132 // ============================================================================================
133 
134 
135 void cbl::measure::twopt::TwoPointCorrelation_projected::read (const std::string dir, const std::string file)
136 {
137  m_dataset->read(dir+file);
138 }
139 
140 
141 // ============================================================================================
142 
143 
144 void cbl::measure::twopt::TwoPointCorrelation_projected::write (const std::string dir, const std::string file, const int rank) const
145 {
146  vector<double> xx = m_dataset->xx();
147 
148  checkDim(xx, m_dd->nbins_D1(), "rp");
149 
150  string header = "[1] perpendicular separation at the bin centre # [2] projected two-point correlation function # [3] error";
151  if (m_compute_extra_info) header += " # [4] mean perpendicular separation # [5] standard deviation of the distribution of perpendicular separations # [6] mean redshift # [7] standard deviation of the redshift distribution";
152 
153  m_dataset->write(dir, file, header, 5, rank);
154 }
155 
156 
157 // ============================================================================================
158 
159 
160 void cbl::measure::twopt::TwoPointCorrelation_projected::measure (const ErrorType errorType, const std::string dir_output_pairs, const std::vector<std::string> dir_input_pairs, const std::string dir_output_resample, const int nMocks, const bool count_dd, const bool count_rr, const bool count_dr, const bool tcount, const Estimator estimator, const double fact, const int seed)
161 {
162  switch (errorType) {
163  case (ErrorType::_Poisson_) :
164  measurePoisson(dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
165  break;
166  case (ErrorType::_Jackknife_) :
167  measureJackknife(dir_output_pairs, dir_input_pairs, dir_output_resample, count_dd, count_rr, count_dr, tcount, estimator, fact);
168  break;
169  case (ErrorType::_Bootstrap_) :
170  measureBootstrap(nMocks, dir_output_pairs, dir_input_pairs, dir_output_resample, count_dd, count_rr, count_dr, tcount, estimator, fact, seed);
171  break;
172  default:
173  ErrorCBL("unknown type of error!", "measure", "TwoPointCorrelation_projected.cpp");
174  }
175 }
176 
177 
178 // ============================================================================================
179 
180 
181 void cbl::measure::twopt::TwoPointCorrelation_projected::measurePoisson (const std::string dir_output_pairs, const std::vector<std::string> dir_input_pairs, const bool count_dd, const bool count_rr, const bool count_dr, const bool tcount, const Estimator estimator, const double fact)
182 {
183  // ----------- measure the 2D two-point correlation function, xi(rp,pi) -----------
184 
185  TwoPointCorrelation2D_cartesian::measurePoisson(dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
186 
187 
188  // ----------- integrate the 2D two-point correlation function along the parallel direction -----------
189 
190  m_dataset = Projected(TwoPointCorrelation2D_cartesian::xx(), TwoPointCorrelation2D_cartesian::yy(), TwoPointCorrelation2D_cartesian::xi2D(), TwoPointCorrelation2D_cartesian::error2D());
191 
192 }
193 
194 
195 // ============================================================================================
196 
197 
198 void cbl::measure::twopt::TwoPointCorrelation_projected::measureJackknife (const std::string dir_output_pairs, const std::vector<std::string> dir_input_pairs, const std::string dir_output_resample, const bool count_dd, const bool count_rr, const bool count_dr, const bool tcount, const Estimator estimator, const double fact)
199 {
200  if (dir_output_resample!=par::defaultString && dir_output_resample!="") {
201  string mkdir = "mkdir -p "+dir_output_resample;
202  if (system(mkdir.c_str())) {}
203  }
204 
205  vector<shared_ptr<data::Data>> data;
206  vector<shared_ptr<pairs::Pair>> dd_regions, rr_regions, dr_regions;
207  count_allPairs_region(dd_regions, rr_regions, dr_regions, TwoPType::_2D_Cartesian_, dir_output_pairs,dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
208 
209  auto data_cart = (estimator==Estimator::_natural_) ? correlation_NaturalEstimator(m_dd, m_rr) : correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
210 
211  if (estimator==Estimator::_natural_)
212  data = XiJackknife(dd_regions, rr_regions);
213  else if (estimator==Estimator::_LandySzalay_)
214  data = XiJackknife(dd_regions, rr_regions, dr_regions);
215  else
216  ErrorCBL("the chosen estimator is not implemented!", "measureJackknife", "TwoPointCorrelation_projected.cpp");
217 
218  vector<vector<double>> ww, covariance;
219  for (size_t i=0; i<data.size(); i++) {
220  ww.push_back(data[i]->data());
221  if (dir_output_resample != par::defaultString && dir_output_resample != "") {
222  string file = "xi_projected_Jackkknife_"+conv(i, par::fINT)+".dat";
223  string header = "[1] perpendicular separation at the bin centre # [2] projected two-point correlation function # [3] error";
224  if (m_compute_extra_info) header += " # [4] mean perpendicular separation # [5] standard deviation of the distribution of perpendicular separations # [6] mean redshift # [7] standard deviation of the redshift distribution";
225  data[i]->write(dir_output_resample, file, header, 10, 0);
226  }
227  }
228 
229  covariance_matrix(ww, covariance, true);
230 
231  vector<double> xx_cart = data_cart->xx(), yy_cart = data_cart->yy();
232  vector<vector<double> > dd_cart, error_cart;
233  data_cart->get_data(dd_cart); data_cart->get_error(error_cart);
234 
235  m_dataset = Projected(xx_cart, yy_cart, dd_cart, error_cart);
236  m_dataset->set_covariance(covariance);
237 
238 }
239 
240 
241 // ============================================================================================
242 
243 
244 void cbl::measure::twopt::TwoPointCorrelation_projected::measureBootstrap (const int nMocks, const std::string dir_output_pairs, const std::vector<std::string> dir_input_pairs, const std::string dir_output_resample, const bool count_dd, const bool count_rr, const bool count_dr, const bool tcount, const Estimator estimator, const double fact, const int seed)
245 {
246  if (nMocks<=0)
247  ErrorCBL("the number of mocks must be >0!", "measureBootstrap", "TwoPointCorrelation1D_monopole.cpp");
248 
249  if (dir_output_resample!=par::defaultString && dir_output_resample!="") {
250  string mkdir = "mkdir -p "+dir_output_resample;
251  if (system(mkdir.c_str())) {}
252  }
253 
254  vector<shared_ptr<data::Data>> data;
255  vector<shared_ptr<pairs::Pair>> dd_regions, rr_regions, dr_regions;
256  count_allPairs_region(dd_regions, rr_regions, dr_regions, TwoPType::_2D_Cartesian_, dir_output_pairs,dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
257 
258  auto data_cart = (estimator==Estimator::_natural_) ? correlation_NaturalEstimator(m_dd, m_rr) : correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
259 
260  if (estimator==Estimator::_natural_)
261  data = XiBootstrap(nMocks, dd_regions, rr_regions, seed);
262  else if (estimator==Estimator::_LandySzalay_)
263  data = XiBootstrap(nMocks, dd_regions, rr_regions, dr_regions, seed);
264  else
265  ErrorCBL("the chosen estimator is not implemented!", "measureBootstrap", "TwoPointCorrelation_projected.cpp");
266 
267  vector<vector<double>> ww, covariance;
268  for (size_t i=0; i<data.size(); i++) {
269  ww.push_back(data[i]->data());
270  if (dir_output_resample!=par::defaultString && dir_output_resample!="") {
271  string file = "xi_projected_Bootstrap_"+conv(i, par::fINT)+".dat";
272  string header = "[1] perpendicular separation at the bin centre # [2] projected two-point correlation function # [3] error";
273  if (m_compute_extra_info) header += " # [4] mean perpendicular separation # [5] standard deviation of the distribution of perpendicular separations # [6] mean redshift # [7] standard deviation of the redshift distribution";
274  data[i]->write(dir_output_resample, file, header, 10, 0);
275  }
276  }
277 
278  covariance_matrix(ww, covariance, false);
279 
280  vector<double> xx_cart = data_cart->xx(), yy_cart = data_cart->yy();
281  vector<vector<double> > dd_cart, error_cart;
282  data_cart->get_data(dd_cart); data_cart->get_error(error_cart);
283 
284  m_dataset = Projected(xx_cart, yy_cart, dd_cart, error_cart);
285  m_dataset->set_covariance(covariance);
286 }
287 
288 
289 // ============================================================================================
290 
291 
292 std::vector<std::shared_ptr<data::Data>> cbl::measure::twopt::TwoPointCorrelation_projected::XiJackknife (const std::vector<std::shared_ptr<pairs::Pair>> dd, const std::vector<std::shared_ptr<pairs::Pair>> rr)
293 {
294  vector<shared_ptr<data::Data>> data;
295 
296  auto data2d = TwoPointCorrelation2D_cartesian::XiJackknife(dd, rr);
297 
298  for (size_t i=0; i<data2d.size(); i++) {
299  vector<double> xx_cart = data2d[i]->xx(), yy_cart = data2d[i]->yy();
300  vector<vector<double>> dd_cart, error_cart;
301  data2d[i]->get_data(dd_cart); data2d[i]->get_error(error_cart);
302  data.push_back(move(Projected(xx_cart, yy_cart, dd_cart, error_cart)));
303  }
304 
305  return data;
306 }
307 
308 
309 // ============================================================================================
310 
311 
312 std::vector<std::shared_ptr<data::Data>> cbl::measure::twopt::TwoPointCorrelation_projected::XiJackknife (const std::vector<std::shared_ptr<pairs::Pair>> dd, const std::vector<std::shared_ptr<pairs::Pair>> rr, const std::vector<std::shared_ptr<pairs::Pair>> dr)
313 {
314  vector<shared_ptr<data::Data>> data;
315 
316  auto data2d = TwoPointCorrelation2D_cartesian::XiJackknife(dd, rr, dr);
317 
318  for (size_t i=0; i<data2d.size(); i++) {
319  vector<double> xx_cart = data2d[i]->xx(), yy_cart = data2d[i]->yy();
320  vector<vector<double>> dd_cart, error_cart;
321  data2d[i]->get_data(dd_cart); data2d[i]->get_error(error_cart);
322  data.push_back(move(Projected(xx_cart, yy_cart, dd_cart, error_cart)));
323  }
324 
325  return data;
326 }
327 
328 
329 // ============================================================================================
330 
331 
332 std::vector<std::shared_ptr<data::Data>> cbl::measure::twopt::TwoPointCorrelation_projected::XiBootstrap (const int nMocks, const std::vector<std::shared_ptr<pairs::Pair>> dd, const std::vector<std::shared_ptr<pairs::Pair>> rr, const int seed)
333 {
334  vector<shared_ptr<data::Data>> data;
335 
336  auto data2d = TwoPointCorrelation2D_cartesian::XiBootstrap(nMocks, dd, rr, seed);
337 
338  for (size_t i=0; i<data2d.size(); i++) {
339  vector<double> xx_cart = data2d[i]->xx(), yy_cart = data2d[i]->yy();
340  vector<vector<double>> dd_cart, error_cart;
341  data2d[i]->get_data(dd_cart); data2d[i]->get_error(error_cart);
342  data.push_back(move(Projected(xx_cart, yy_cart, dd_cart, error_cart)));
343  }
344 
345  return data;
346 }
347 
348 
349 // ============================================================================================
350 
351 
352 std::vector<std::shared_ptr<data::Data>> cbl::measure::twopt::TwoPointCorrelation_projected::XiBootstrap (const int nMocks, const std::vector<std::shared_ptr<pairs::Pair>> dd, const std::vector<std::shared_ptr<pairs::Pair>> rr, const std::vector<std::shared_ptr<pairs::Pair>> dr, const int seed)
353 {
354  vector<shared_ptr<data::Data>> data;
355 
356  auto data2d = TwoPointCorrelation2D_cartesian::XiBootstrap(nMocks, dd, rr, dr, seed);
357 
358  for (size_t i=0; i<data2d.size(); i++) {
359  vector<double> xx_cart = data2d[i]->xx(), yy_cart = data2d[i]->yy();
360  vector<vector<double>> dd_cart, error_cart;
361  data2d[i]->get_data(dd_cart); data2d[i]->get_error(error_cart);
362  data.push_back(move(Projected(xx_cart, yy_cart, dd_cart, error_cart)));
363  }
364 
365  return data;
366 }
367 
368 
369 // ============================================================================
370 
371 
372 void cbl::measure::twopt::TwoPointCorrelation_projected::read_covariance (const std::string dir, const std::string file)
373 {
374  m_dataset->set_covariance(dir+file);
375 }
376 
377 
378 // ============================================================================
379 
380 
381 void cbl::measure::twopt::TwoPointCorrelation_projected::write_covariance (const std::string dir, const std::string file) const
382 {
383  m_dataset->write_covariance(dir, file);
384 }
385 
386 
387 // ============================================================================
388 
389 
390 void cbl::measure::twopt::TwoPointCorrelation_projected::compute_covariance (const std::vector<std::shared_ptr<data::Data>> xi, const bool JK)
391 {
392  vector<vector<double>> Xi;
393 
394  for (size_t i=0; i<xi.size(); i++)
395  Xi.push_back(xi[i]->data());
396 
397  vector<vector<double>> cov_mat;
398  cbl::covariance_matrix(Xi, cov_mat, JK);
399 
400  m_dataset->set_covariance(cov_mat);
401 }
402 
403 
404 // ============================================================================
405 
406 
407 void cbl::measure::twopt::TwoPointCorrelation_projected::compute_covariance (const std::vector<std::string> file, const bool JK)
408 {
409  vector<double> rad, mean;
410  vector<vector<double>> cov_mat;
411 
412  cbl::covariance_matrix(file, rad, mean, cov_mat, JK);
413  m_dataset->set_covariance(cov_mat);
414 }
The class Data1D_extra.
The class TwoPointCorrelation_projected.
The class Data1D_extra.
Definition: Data1D_extra.h:52
The class Data1D.
Definition: Data1D.h:51
void write_covariance(const std::string dir, const std::string file) const override
write the measured two-point correlation
void measurePoisson(const std::string dir_output_pairs=par::defaultString, const std::vector< std::string > dir_input_pairs={}, const bool count_dd=true, const bool count_rr=true, const bool count_dr=true, const bool tcount=true, const Estimator estimator=Estimator::_LandySzalay_, const double fact=0.1) override
measure the projected two-point correlation function with Poisson errors
std::shared_ptr< data::Data > Projected(const std::vector< double > rp, const std::vector< double > pi, const std::vector< std::vector< double >> xi, const std::vector< std::vector< double >> error_xi) override
measure projected correlation function
void measure(const ErrorType errorType=ErrorType::_Poisson_, const std::string dir_output_pairs=par::defaultString, const std::vector< std::string > dir_input_pairs={}, const std::string dir_output_resample=par::defaultString, const int nMocks=0, const bool count_dd=true, const bool count_rr=true, const bool count_dr=true, const bool tcount=true, const Estimator estimator=Estimator::_LandySzalay_, const double fact=0.1, const int seed=3213) override
measure the projected two-point correlation function
void measureJackknife(const std::string dir_output_pairs=par::defaultString, const std::vector< std::string > dir_input_pairs={}, const std::string dir_output_resample=par::defaultString, const bool count_dd=true, const bool count_rr=true, const bool count_dr=true, const bool tcount=true, const Estimator estimator=Estimator::_LandySzalay_, const double fact=0.1) override
measure the projected two-point correlation function estimating the covariance with Jackknife resampl...
std::shared_ptr< data::Data > data_with_extra_info(const std::vector< double > rp, const std::vector< double > ww, const std::vector< double > error) const
return a data object with extra info
std::vector< std::shared_ptr< data::Data > > XiBootstrap(const int nMocks, const std::vector< std::shared_ptr< pairs::Pair >> dd, const std::vector< std::shared_ptr< pairs::Pair >> rr, const int seed=3213) override
measure the bootstrap resampling of the two-point correlation function, ξ(r)
void read(const std::string dir, const std::string file) override
read the projected two-point correlation function
void write(const std::string dir=par::defaultString, const std::string file=par::defaultString, const int rank=0) const override
write the projected two-point correlation function
void compute_covariance(const std::vector< std::shared_ptr< data::Data >> xi, const bool JK) override
compute the covariance matrix
void measureBootstrap(const int nMocks, const std::string dir_output_pairs, const std::vector< std::string > dir_input_pairs={}, const std::string dir_output_resample=par::defaultString, const bool count_dd=true, const bool count_rr=true, const bool count_dr=true, const bool tcount=true, const Estimator estimator=Estimator::_LandySzalay_, const double fact=0.1, const int seed=3213) override
measure the projected two-point correlation function estimating the covariance with Bootstrap resampl...
void read_covariance(const std::string dir, const std::string file) override
read the measured covariance matrix
std::vector< std::shared_ptr< data::Data > > XiJackknife(const std::vector< std::shared_ptr< pairs::Pair >> dd, const std::vector< std::shared_ptr< pairs::Pair >> rr) override
measure the jackknife resampling of the two-point correlation function, ξ(r)
std::shared_ptr< pairs::Pair > m_dd
number of data-data pairs
static const char fINT[]
conversion symbol for: int -> std::string
Definition: Constants.h:121
static const std::string defaultString
default std::string value
Definition: Constants.h:336
static const double pi
: the ratio of a circle's circumference to its diameter
Definition: Constants.h:199
Estimator
the two-point correlation estimator
ErrorType
the two-point correlation function error type
Definition: Measure.h:59
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38
T Min(const std::vector< T > vect)
minimum element of a std::vector
Definition: Kernel.h:1324
std::string conv(const T val, const char *fact)
convert a number to a std::string
Definition: Kernel.h:903
void checkDim(const std::vector< T > vect, const int val, const std::string vector, bool equal=true)
check if the dimension of a std::vector is equal/lower than an input value
Definition: Kernel.h:1532
int ErrorCBL(const std::string msg, const std::string functionCBL, const std::string fileCBL, const cbl::glob::ExitCode exitCode=cbl::glob::ExitCode::_error_)
throw an exception: it is used for handling exceptions inside the CosmoBolognaLib
Definition: Kernel.h:780
void covariance_matrix(const std::vector< std::vector< double >> mat, std::vector< std::vector< double >> &cov, const bool JK=false)
compute the covariance matrix from an input dataset
Definition: Func.cpp:761
int nint(const T val)
the nearest integer
Definition: Kernel.h:915