CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
TwoPointCorrelation2D.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 
37 #include "TwoPointCorrelation2D.h"
38 #include "Data2D_extra.h"
39 
40 using namespace std;
41 
42 using namespace cbl;
43 using namespace catalogue;
44 using namespace chainmesh;
45 using namespace data;
46 using namespace pairs;
47 using namespace measure;
48 using namespace twopt;
49 
50 
51 // ===================================================================================================
52 
53 
54 void cbl::measure::twopt::TwoPointCorrelation2D::write_pairs (const std::shared_ptr<pairs::Pair> PP, const std::string dir, const std::string file) const
55 {
56  string MK = "mkdir -p "+dir; if (system (MK.c_str())) {}
57 
58  string file_out = dir+file;
59  ofstream fout(file_out.c_str()); checkIO(fout, file_out);
60 
61 
62  // ----- standard info: scales at the bin centre + number of pairs -----
63 
64  if (PP->pairInfo()==PairInfo::_standard_)
65  for (int i=0; i<PP->nbins_D1(); i++)
66  for (int j=0; j<PP->nbins_D2(); j++)
67  fout << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << i
68  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << j
69  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP->scale_D1(i)
70  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP->scale_D2(j)
71  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP->PP2D(i, j)
72  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP->PP2D_weighted(i, j) << endl;
73 
74 
75  // ----- standard + extra info -----
76 
77  else if (PP->pairInfo()==PairInfo::_extra_)
78  for (int i=0; i<PP->nbins_D1(); i++)
79  for (int j=0; j<PP->nbins_D2(); j++)
80  fout << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << i
81  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << j
82  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP->scale_D1(i)
83  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP->scale_D2(j)
84  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP->PP2D(i, j)
85  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP->PP2D_weighted(i, j)
86  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP->scale_D1_mean(i, j)
87  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP->scale_D1_sigma(i, j)
88  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP->scale_D2_mean(i, j)
89  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP->scale_D2_sigma(i, j)
90  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP->z_mean(i, j)
91  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP->z_sigma(i, j) << endl;
92 
93  else
94  ErrorCBL("no such pairInfo!", "write_pairs", "TwoPointCorrelation2D.cpp");
95 
96  fout.clear(); fout.close(); coutCBL << "I wrote the file " << file_out << endl;
97 
98 }
99 
100 
101 // ============================================================================
102 
103 
104 void cbl::measure::twopt::TwoPointCorrelation2D::read_pairs (std::shared_ptr<pairs::Pair> PP, const std::vector<std::string> dir, const std::string file) const
105 {
106  if (dir.size()==0)
107  ErrorCBL("dir.size()=0!", "read_pairs", "TwoPointCorrelation2D.cpp");
108 
109  string line;
110  int i, j;
111  double pairs, weighted_pairs;
112 
113 
114  // ----- standard info: scales at the bin centre + number of pairs -----
115 
116  if (PP->pairInfo()==PairInfo::_standard_)
117 
118  for (size_t dd=0; dd<dir.size(); dd++) {
119 
120  string file_in = dir[dd]+file;
121  ifstream fin(file_in.c_str()); checkIO(fin, file_in);
122 
123  while (getline(fin, line)) {
124 
125  stringstream ss(line);
126  vector<double> num;
127  double val;
128  while (ss >> val) num.emplace_back(val);
129 
130  if (num.size()!=6)
131  ErrorCBL("the number of lines in the input pair file: "+file_in+" must be 6!", "read_pairs", "TwoPointCorrelation2D.cpp");
132 
133  i = int(num[0]);
134  j = int(num[1]);
135  pairs = num[4];
136  weighted_pairs = num[5];
137 
138  PP->add_data2D(i, j, {pairs, weighted_pairs});
139 
140  }
141 
142  fin.clear(); fin.close(); coutCBL << "I read the file " << file_in << endl;
143  }
144 
145 
146  // ----- standard + extra info -----
147 
148  else if (PP->pairInfo()==PairInfo::_extra_)
149 
150  for (size_t dd=0; dd<dir.size(); dd++) {
151 
152  string file_in = dir[dd]+file;
153  ifstream fin(file_in.c_str()); checkIO(fin, file_in);
154 
155  double scale_D1_mean, scale_D1_sigma, scale_D2_mean, scale_D2_sigma, z_mean, z_sigma;
156 
157  while (getline(fin, line)) {
158 
159  stringstream ss(line);
160  vector<double> num;
161  double val;
162  while (ss >> val) num.emplace_back(val);
163 
164  if (num.size()!=12)
165  ErrorCBL("the number of lines in the input pair file: "+file_in+" must be 12!", "read_pairs", "TwoPointCorrelation2D.cpp");
166 
167  i = int(num[0]);
168  j = int(num[1]);
169  pairs = num[4];
170  weighted_pairs = num[5];
171  scale_D1_mean = num[6];
172  scale_D1_sigma = num[7];
173  scale_D2_mean = num[8];
174  scale_D2_sigma = num[9];
175  z_mean = num[10];
176  z_sigma = num[11];
177 
178  PP->add_data2D(i, j, {pairs, weighted_pairs, scale_D1_mean, pow(scale_D1_sigma, 2)*weighted_pairs, scale_D2_mean, pow(scale_D2_sigma, 2)*weighted_pairs, z_mean, pow(z_sigma, 2)*weighted_pairs});
179 
180  }
181 
182  fin.clear(); fin.close(); coutCBL << "I read the file " << file_in << endl;
183 
184  }
185 
186  else
187  ErrorCBL("no such pairInfo!", "read_pairs", "TwoPointCorrelation2D.cpp");
188 }
189 
190 
191 // ============================================================================
192 
193 
194 void cbl::measure::twopt::TwoPointCorrelation2D::write_pairs (const std::vector<std::shared_ptr<pairs::Pair>> PP, const std::string dir, const std::string file) const
195 {
196  size_t nRegions = m_data->region_list().size();
197 
198  bool cross = (PP.size()==nRegions*nRegions) ? true : false;
199 
200  string MK = "mkdir -p "+dir; if (system (MK.c_str())) {}
201 
202  string file_out = dir+file;
203  ofstream fout(file_out.c_str()); checkIO(fout, file_out);
204 
205 
206  // ----- standard info: scales at the bin centre + number of pairs -----
207 
208  if (PP[0]->pairInfo()==PairInfo::_standard_)
209 
210  for (size_t i=0; i<nRegions; i++)
211  for (size_t j=(cross) ? 0 : i; j<nRegions; j++) {
212  int index = (cross) ? i*nRegions+j : i*nRegions+j-(i-1)*i/2-i;
213  for (int r1=0; r1<PP[index]->nbins_D1(); r1++)
214  for (int r2=0; r2<PP[index]->nbins_D2(); r2++)
215  if (PP[index]->PP2D(r1, r2)>0)
216  fout << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << i
217  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << j
218  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << r1
219  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << r2
220  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP[index]->scale_D1(r1)
221  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP[index]->scale_D2(r2)
222  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP[index]->PP2D(r1, r2)
223  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP[index]->PP2D_weighted(r1, r2) << endl;
224  }
225 
226 
227  // ----- standard + extra info -----
228 
229  else if (PP[0]->pairInfo()==PairInfo::_extra_)
230 
231  for (size_t i=0; i<nRegions; i++)
232  for (size_t j=(cross) ? 0 : i; j<nRegions; j++) {
233  int index = (cross) ? i*nRegions+j : i*nRegions+j-(i-1)*i/2-i;
234  for (int r1=0; r1<PP[index]->nbins_D1(); r1++)
235  for (int r2=0; r2<PP[index]->nbins_D2(); r2++)
236  if (PP[index]->PP2D(r1, r2)>0)
237  fout << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << i
238  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << j
239  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << r1
240  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << r2
241  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP[index]->scale_D1(r1)
242  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP[index]->scale_D2(r2)
243  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP[index]->PP2D(r1, r2)
244  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP[index]->PP2D_weighted(r1, r2)
245  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP[index]->scale_D1_mean(r1, r2)
246  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP[index]->scale_D1_sigma(r1, r2)
247  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP[index]->scale_D2_mean(r1, r2)
248  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP[index]->scale_D2_sigma(r1, r2)
249  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP[index]->z_mean(r1, r2)
250  << " " << setiosflags(ios::fixed) << setprecision(5) << setw(10) << right << PP[index]->z_sigma(r1, r2) << endl;
251  }
252 
253  else
254  ErrorCBL("no such pairInfo!", "write_pairs", "TwoPointCorrelation2D.cpp");
255 
256  fout.clear(); fout.close();
257 
258 }
259 
260 
261 // ============================================================================
262 
263 
264 void cbl::measure::twopt::TwoPointCorrelation2D::read_pairs (std::vector<std::shared_ptr<pairs::Pair>> PP, const std::vector<std::string> dir, const std::string file) const
265 {
266  size_t nRegions = m_data->region_list().size();
267 
268  bool cross = (PP.size() == nRegions*nRegions) ? true : false;
269 
270  int i, j, r1, r2, index;
271  double rad1, rad2, pairs, weighted_pairs;
272 
273 
274  // ----- standard info: scales at the bin centre + number of pairs -----
275 
276  if (PP[0]->pairInfo()==PairInfo::_standard_)
277 
278  for (size_t dd=0; dd<dir.size(); dd++) {
279 
280  string ff = dir[dd]+file;
281  coutCBL << "I'm reading the pair file: " << ff << endl;
282  ifstream fin(ff.c_str()); checkIO(fin, ff);
283 
284  while (fin >> i >> j >> r1 >> r2 >> rad1 >> rad2 >> pairs >> weighted_pairs) {
285  index = (cross) ? i*nRegions+j : i*nRegions-(i-1)*i/2+j-i;
286  PP[index]->add_data2D(r1, r2, {pairs, weighted_pairs});
287  }
288 
289  fin.clear(); fin.close(); coutCBL << "I read the file " << ff << endl;
290  }
291 
292 
293  // ----- standard + extra info -----
294 
295  else if (PP[0]->pairInfo()==PairInfo::_extra_)
296 
297  for (size_t dd=0; dd<dir.size(); dd++) {
298 
299  string ff = dir[dd]+file;
300  coutCBL << "I'm reading the pair file: " << ff << endl;
301  ifstream fin(ff.c_str()); checkIO(fin, ff);
302 
303  double scale_D1_mean, scale_D1_sigma, scale_D2_mean, scale_D2_sigma, z_mean, z_sigma;
304 
305  while (fin >> i >> j >> r1 >> r2 >> rad1 >> rad2 >> pairs >> weighted_pairs >> scale_D1_mean >> scale_D1_sigma >> scale_D2_mean >> scale_D2_sigma >> z_mean >> z_sigma) {
306  index = (cross) ? i*nRegions+j : i*nRegions-(i-1)*i/2+j-i;
307  PP[index]->add_data2D(r1, r2, {pairs, weighted_pairs, scale_D1_mean, pow(scale_D1_sigma, 2)*weighted_pairs, scale_D2_mean, pow(scale_D2_sigma, 2)*weighted_pairs, z_mean, pow(z_sigma, 2)*weighted_pairs});
308  }
309 
310  fin.clear(); fin.close(); coutCBL << "I read the file " << ff << endl;
311  }
312 
313  else
314  ErrorCBL("no such pairInfo!", "read_pairs", "TwoPointCorrelation2D.cpp");
315 
316 }
317 
318 
319 // ============================================================================
320 
321 
322 std::shared_ptr<data::Data> cbl::measure::twopt::TwoPointCorrelation2D::data_with_extra_info (const std::shared_ptr<pairs::Pair> dd, const std::vector<double> scale_D1, const std::vector<double> scale_D2, const std::vector<std::vector<double>> xi, const std::vector<std::vector<double>> error) const
323 {
324  vector<vector<double>> extra(6);
325 
326  for (int i=0; i<dd->nbins_D1(); ++i)
327  for (int j=0; j<dd->nbins_D2(); ++j) {
328  extra[0].push_back(dd->scale_D1_mean(i, j));
329  extra[1].push_back(dd->scale_D1_sigma(i, j));
330  extra[2].push_back(dd->scale_D2_mean(i, j));
331  extra[3].push_back(dd->scale_D2_sigma(i, j));
332  extra[4].push_back(dd->z_mean(i, j));
333  extra[5].push_back(dd->z_sigma(i, j));
334  }
335 
336  return move(unique_ptr<data::Data2D_extra>(new data::Data2D_extra(scale_D1, scale_D2, xi, error, extra)));
337 }
338 
339 
340 // ============================================================================
341 
342 
343 std::shared_ptr<Data> cbl::measure::twopt::TwoPointCorrelation2D::correlation_NaturalEstimator (const std::shared_ptr<pairs::Pair> dd, const std::shared_ptr<pairs::Pair> rr, const int nData, const double nData_weighted, const int nRandom, const double nRandom_weighted)
344 {
345  vector<double> scale_D1, scale_D2;
346  vector<vector<double>> xi, error;
347  scale_D1.resize(m_dd->nbins_D1());
348  scale_D2.resize(m_dd->nbins_D2());
349 
350  xi.resize(m_dd->nbins_D1(), vector<double>(m_dd->nbins_D2(), 0));
351  error.resize(m_dd->nbins_D1(), vector<double>(m_dd->nbins_D2(), 0));
352 
353  // number of objects in the data catalogue
354  int nD = (nData>0) ? nData : m_data->nObjects();
355 
356  // weighted number of objects in the data catalogue
357  double nDw = (nData_weighted>0) ? nData_weighted : m_data->weightedN();
358 
359  // number of objects in the random catalogue
360  int nR = (nRandom>0) ? nRandom : m_random->nObjects();
361 
362  // weighted number of objects in the random catalogue
363  double nRw = (nRandom_weighted>0) ? nRandom_weighted : m_random->weightedN();
364 
365  // inverse of the total number of data-data pairs
366  double nDDi = 1./(nDw*(nDw-1.)*0.5);
367 
368  // inverse of the total number of random-random pairs
369  double nRRi = 1./(nRw*(nRw-1.)*0.5);
370 
371  for (int i=0; i<dd->nbins_D1(); i++) {
372  scale_D1[i] = dd->scale_D1(i);
373  for (int j=0; j<dd->nbins_D2(); j++) {
374  scale_D2[j] = dd->scale_D2(j);
375 
376  xi[i][j] = -1.;
377  error[i][j] = 1000.;
378 
379  if (dd->PP2D_weighted(i, j)>0) {
380 
381  if (rr->PP2D_weighted(i, j)<1.e-30)
382  ErrorCBL("there are no random objects in the bin "+conv(i, par::fINT)+","+conv(j, par::fINT)+"; please, either increase the total number of random objects or enlarge the bin size! (dd="+conv(dd->PP2D_weighted(i, j), par::fDP3)+", rr="+conv(rr->PP2D_weighted(i, j), par::fDP3)+")!", "correlation_NaturalEstimator", "TwoPointCorrelation2D.cpp");
383 
384  // normalised number of data-data weighted pairs
385  double DD_norm = dd->PP2D_weighted(i, j)*nDDi;
386 
387  // normalised number of random-random weighted pairs
388  double RR_norm = rr->PP2D_weighted(i, j)*nRRi;
389 
390  // natural estimator
391  xi[i][j] = max(-1., DD_norm/RR_norm-1.);
392 
393  // Poisson error
394  error[i][j]= PoissonError(Estimator::_natural_, dd->PP2D(i, j), rr->PP2D(i, j), 0, nD, nR);
395  }
396  }
397  }
398 
399  return (!m_compute_extra_info) ? move(unique_ptr<Data2D>(new Data2D(scale_D1, scale_D2, xi, error))) : data_with_extra_info(dd, scale_D1, scale_D2, xi, error);
400 }
401 
402 
403 // ============================================================================
404 
405 
406 std::shared_ptr<Data> cbl::measure::twopt::TwoPointCorrelation2D::correlation_LandySzalayEstimator (const std::shared_ptr<pairs::Pair> dd, const std::shared_ptr<pairs::Pair> rr, const std::shared_ptr<pairs::Pair> dr, const int nData, const double nData_weighted, const int nRandom, const double nRandom_weighted)
407 {
408  vector<double> scale_D1, scale_D2;
409  vector<vector<double>> xi, error;
410  scale_D1.resize(m_dd->nbins_D1());
411  scale_D2.resize(m_dd->nbins_D2());
412 
413  xi.resize(m_dd->nbins_D1(), vector<double>(m_dd->nbins_D2(), 0));
414  error.resize(m_dd->nbins_D1(), vector<double>(m_dd->nbins_D2(), 0));
415 
416  // number of objects in the data catalogue
417  int nD = (nData>0) ? nData : m_data->nObjects();
418 
419  // weighted number of objects in the data catalogue
420  double nDw = (nData_weighted>0) ? nData_weighted : m_data->weightedN();
421 
422  // number of objects in the random catalogue
423  int nR = (nRandom>0) ? nRandom : m_random->nObjects();
424 
425  // weighted number of objects in the random catalogue
426  double nRw = (nRandom_weighted>0) ? nRandom_weighted : m_random->weightedN();
427 
428  // inverse of the total number of data-data pairs
429  double nDDi = 1./(nDw*(nDw-1.)*0.5);
430 
431  // inverse of the total number of random-random pairs
432  double nRRi = 1./(nRw*m_random_dilution_fraction*(nRw*m_random_dilution_fraction-1.)*0.5);
433 
434  // inverse of the total number of data-random pairs
435  double nDRi = 1./(nDw*nRw);
436 
437  for (int i=0; i<dd->nbins_D1(); i++) {
438  scale_D1[i] = dd->scale_D1(i);
439  for (int j=0; j<dd->nbins_D2(); j++) {
440  scale_D2[j] = dd->scale_D2(j);
441 
442  xi[i][j] = -1.;
443  error[i][j] = 1000.;
444 
445  if (dd->PP2D_weighted(i, j)>0) {
446 
447  if (rr->PP2D_weighted(i, j)<1.e-30)
448  ErrorCBL("there are no random objects in the bin "+conv(i, par::fINT)+","+conv(j, par::fINT)+"; please, either increase the total number of random objects or enlarge the bin size! (dd="+conv(dd->PP2D_weighted(i, j), par::fDP3)+", rr="+conv(rr->PP2D_weighted(i, j), par::fDP3)+")!", "correlation_LandySzalayEstimator", "TwoPointCorrelation2D.cpp");
449 
450  // normalised number of data-data weighted pairs
451  double DD_norm = dd->PP2D_weighted(i, j)*nDDi;
452 
453  // normalised number of random-random weighted pairs
454  double RR_norm = rr->PP2D_weighted(i, j)*nRRi;
455 
456  // normalised number of data-random weighted pairs
457  double DR_norm = dr->PP2D_weighted(i, j)*nDRi;
458 
459  // Landy & Szalay estimator
460  xi[i][j] = max(-1., (DD_norm-2.*DR_norm)/RR_norm+1.);
461 
462  // Poisson error
463  error[i][j]= PoissonError(Estimator::_LandySzalay_, dd->PP2D(i, j), rr->PP2D(i, j), dr->PP2D(i, j), nD, nR);
464 
465  }
466  }
467  }
468 
469  return (!m_compute_extra_info) ? move(unique_ptr<Data2D>(new Data2D(scale_D1, scale_D2, xi, error))) : data_with_extra_info(dd, scale_D1, scale_D2, xi, error);
470 }
471 
472 
473 // ============================================================================
474 
475 
476 std::vector<std::shared_ptr<Data>> cbl::measure::twopt::TwoPointCorrelation2D::XiJackknife (const std::vector<std::shared_ptr<pairs::Pair>> dd, const std::vector<std::shared_ptr<pairs::Pair>> rr)
477 {
478  vector<long> region_list = m_data->region_list();
479  size_t nRegions = region_list.size();
480 
481  vector<shared_ptr<Data>> data;
482 
483  for (size_t i=0; i<nRegions; i++) {
484 
485  coutCBL << "analysing region: " << i << " of " << nRegions << "\r"; cout.flush();
486 
487  auto dd_SS = Pair::Create(m_dd->pairType(), m_dd->pairInfo(), m_dd->sMin_D1(), m_dd->sMax_D1(), m_dd->nbins_D1(), m_dd->shift_D1(), m_dd->sMin_D2(), m_dd->sMax_D2(), m_dd->nbins_D2(), m_dd->shift_D2(), m_dd->angularUnits(), m_dd->angularWeight());
488  auto rr_SS = Pair::Create(m_rr->pairType(), m_rr->pairInfo(), m_rr->sMin_D1(), m_rr->sMax_D1(), m_rr->nbins_D1(), m_rr->shift_D1(), m_rr->sMin_D2(), m_rr->sMax_D2(), m_rr->nbins_D2(), m_rr->shift_D2(), m_rr->angularUnits(), m_rr->angularWeight());
489 
490  vector<int> w(nRegions, 1);
491  w[i] = 0;
492 
493  for (size_t j=0; j<nRegions; j++)
494  for (size_t k=j; k<nRegions; k++) {
495  int index = j*nRegions-(j-1)*j/2+k-j;
496  double ww = w[j]*w[k];
497  if (ww>0)
498  for (int bin1=0; bin1<dd_SS->nbins_D1(); bin1++)
499  for (int bin2=0; bin2<dd_SS->nbins_D2(); bin2++) {
500  dd_SS->add_data2D(bin1, bin2, dd[index]);
501  rr_SS->add_data2D(bin1, bin2, rr[index]);
502  }
503  }
504 
505 
506  int nData_SS = m_data->nObjects_condition(Var::_Region_, region_list[i], region_list[i]+1, 1);
507  double nData_SS_weighted = m_data->weightedN_condition(Var::_Region_, region_list[i], region_list[i]+1, 1);
508  int nRandom_SS = m_random->nObjects_condition(Var::_Region_, region_list[i], region_list[i]+1, 1);
509  double nRandom_SS_weighted = m_random->weightedN_condition(Var::_Region_, region_list[i], region_list[i]+1, 1);
510 
511  data.push_back(move(correlation_NaturalEstimator(dd_SS, rr_SS, nData_SS, nData_SS_weighted, nRandom_SS, nRandom_SS_weighted)));
512  }
513 
514  return data;
515 }
516 
517 
518 // ============================================================================
519 
520 
521 std::vector<std::shared_ptr<Data>> cbl::measure::twopt::TwoPointCorrelation2D::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)
522 {
523  vector<long> region_list = m_data->region_list();
524  size_t nRegions = region_list.size();
525 
526  vector<shared_ptr<Data>> data;
527 
528  for (size_t i=0; i<nRegions; i++) {
529 
530  coutCBL << "analysing region: " << i << " of " << nRegions << "\r"; cout.flush();
531 
532  auto dd_SS = Pair::Create(m_dd->pairType(), m_dd->pairInfo(), m_dd->sMin_D1(), m_dd->sMax_D1(), m_dd->nbins_D1(), m_dd->shift_D1(), m_dd->sMin_D2(), m_dd->sMax_D2(), m_dd->nbins_D2(), m_dd->shift_D2(), m_dd->angularUnits(), m_dd->angularWeight());
533  auto rr_SS = Pair::Create(m_rr->pairType(), m_rr->pairInfo(), m_rr->sMin_D1(), m_rr->sMax_D1(), m_rr->nbins_D1(), m_rr->shift_D1(), m_rr->sMin_D2(), m_rr->sMax_D2(), m_rr->nbins_D2(), m_rr->shift_D2(), m_rr->angularUnits(), m_rr->angularWeight());
534  auto dr_SS = Pair::Create(m_dr->pairType(), m_dr->pairInfo(), m_dr->sMin_D1(), m_dr->sMax_D1(), m_dr->nbins_D1(), m_dr->shift_D1(), m_dr->sMin_D2(), m_dr->sMax_D2(), m_dr->nbins_D2(), m_dr->shift_D2(), m_dr->angularUnits(), m_dr->angularWeight());
535 
536  vector<int> w(nRegions, 1);
537  w[i] = 0;
538 
539  for (size_t j=0; j<nRegions; j++){
540 
541  if(w[j]>0){
542  for (size_t k=j; k<nRegions; k++) { // auto pairs
543  if(w[k]>0){
544  int index = j*nRegions+k-(j-1)*j/2-j;
545  for (int bin1=0; bin1<dd_SS->nbins_D1(); bin1++)
546  for (int bin2=0; bin2<dd_SS->nbins_D2(); bin2++) {
547  dd_SS->add_data2D(bin1, bin2, dd[index]);
548  rr_SS->add_data2D(bin1, bin2, rr[index]);
549  }
550  }
551  }
552 
553  for (size_t k=0; k<nRegions; k++) {// cross pairs
554  if (w[k]>0){
555  int index = j*nRegions+k;
556  for (int bin1=0; bin1<dd_SS->nbins_D1(); bin1++)
557  for (int bin2=0; bin2<dd_SS->nbins_D2(); bin2++) {
558  dr_SS->add_data2D(bin1, bin2, dr[index]);
559  }
560  }
561  }
562  }
563  }
564 
565  double nData_SS = m_data->nObjects_condition(Var::_Region_, region_list[i], region_list[i]+1, 1);
566  double nData_SS_weighted = m_data->weightedN_condition(Var::_Region_, region_list[i], region_list[i]+1, 1);
567  double nRandom_SS = m_random->nObjects_condition(Var::_Region_, region_list[i], region_list[i]+1, 1);
568  double nRandom_SS_weighted = m_random->weightedN_condition(Var::_Region_, region_list[i], region_list[i]+1, 1);
569 
570  data.push_back(move(correlation_LandySzalayEstimator(dd_SS, rr_SS, dr_SS, nData_SS, nData_SS_weighted, nRandom_SS, nRandom_SS_weighted)));
571  }
572 
573  return data;
574 }
575 
576 
577 // ============================================================================
578 
579 
580 std::vector<std::shared_ptr<Data>> cbl::measure::twopt::TwoPointCorrelation2D::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)
581 {
582  vector<long> region_list = m_data->region_list();
583  size_t nRegions = region_list.size();
584 
585  vector<shared_ptr<Data>> data;
586  vector<double> nData_reg, nData_reg_weighted, nRandom_reg, nRandom_reg_weighted;
587 
588  for (size_t i=0; i<nRegions; i++) {
589  nData_reg.push_back(m_data->nObjects_condition(Var::_Region_, region_list[i], region_list[i]+1, 0));
590  nData_reg_weighted.push_back(m_data->weightedN_condition(Var::_Region_, region_list[i], region_list[i]+1, 0));
591  nRandom_reg.push_back(m_random->nObjects_condition(Var::_Region_, region_list[i], region_list[i]+1, 0));
592  nRandom_reg_weighted.push_back(m_random->weightedN_condition(Var::_Region_, region_list[i], region_list[i]+1, 0));
593  }
594 
595  random::UniformRandomNumbers_Int ran(0., nRegions-1, seed);
596 
597  int val = 3; // see Norberg et al. 2009
598 
599  int nbins_D1 = m_dd->nbins_D1();
600  int nbins_D2 = m_dd->nbins_D2();
601 
602  for (int i=0; i<nMocks; i++) {
603 
604  coutCBL << "analysing mock: " << i << " of " << nMocks << "\r"; cout.flush();
605 
606  auto dd_SS = Pair::Create(m_dd->pairType(), m_dd->pairInfo(), m_dd->sMin_D1(), m_dd->sMax_D1(), m_dd->nbins_D1(), m_dd->shift_D1(), m_dd->sMin_D2(), m_dd->sMax_D2(), m_dd->nbins_D2(), m_dd->shift_D2(), m_dd->angularUnits(), m_dd->angularWeight());
607  auto rr_SS = Pair::Create(m_rr->pairType(), m_rr->pairInfo(), m_rr->sMin_D1(), m_rr->sMax_D1(), m_rr->nbins_D1(), m_rr->shift_D1(), m_rr->sMin_D2(), m_rr->sMax_D2(), m_rr->nbins_D2(), m_rr->shift_D2(), m_rr->angularUnits(), m_rr->angularWeight());
608 
609  double nData_SS = 0., nData_SS_weighted = 0., nRandom_SS = 0., nRandom_SS_weighted = 0.;
610 
611  vector<int> w(nRegions, 0);
612  for (size_t n=0; n<val*nRegions; n++)
613  w[ran()] ++;
614 
615  for (size_t j=0; j<nRegions; j++) {
616 
617  if (w[j]>0) {
618 
619  nData_SS += w[j]*nData_reg[j];
620  nData_SS_weighted += w[j]*nData_reg_weighted[j];
621  nRandom_SS += w[j]*nRandom_reg[j];
622  nRandom_SS_weighted += w[j]*nRandom_reg_weighted[j];
623 
624  for (size_t k=j; k<nRegions; k++) {
625  if(w[k]>0){
626  int index = j*nRegions-(j-1)*j/2+k-j;
627  double ww = w[j]*w[k];
628  for (int bin1=0; bin1<nbins_D1; bin1++)
629  for (int bin2=0; bin2<nbins_D2; bin2++) {
630  dd_SS->add_data2D(bin1, bin2, dd[index], ww);
631  rr_SS->add_data2D(bin1, bin2, rr[index], ww);
632  }
633  }
634  }
635  }
636  }
637 
638  data.push_back(move(correlation_NaturalEstimator(dd_SS, rr_SS, nData_SS, nData_SS_weighted, nRandom_SS, nRandom_SS_weighted)));
639  }
640 
641 
642  return data;
643 }
644 
645 
646 // ============================================================================
647 
648 
649 std::vector<std::shared_ptr<Data>> cbl::measure::twopt::TwoPointCorrelation2D::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)
650 {
651  vector<long> region_list = m_data->region_list();
652  size_t nRegions = region_list.size();
653 
654  vector<shared_ptr<Data>> data;
655  vector<double> nData_reg, nData_reg_weighted, nRandom_reg, nRandom_reg_weighted;
656 
657  for (size_t i=0; i<nRegions; i++) {
658  nData_reg.push_back(m_data->nObjects_condition(Var::_Region_, region_list[i], region_list[i]+1, 0));
659  nData_reg_weighted.push_back(m_data->weightedN_condition(Var::_Region_, region_list[i], region_list[i]+1, 0));
660  nRandom_reg.push_back(m_random->nObjects_condition(Var::_Region_, region_list[i], region_list[i]+1, 0));
661  nRandom_reg_weighted.push_back(m_random->weightedN_condition(Var::_Region_, region_list[i], region_list[i]+1, 0));
662  }
663 
664  random::UniformRandomNumbers_Int ran(0., nRegions-1, seed);
665 
666  int val = 3; // see Norberg et al. 2009
667 
668  int nbins_D1 = m_dd->nbins_D1();
669  int nbins_D2 = m_dd->nbins_D2();
670 
671  for (int i=0; i<nMocks; i++) {
672 
673  coutCBL << "analysing mock: " << i << " of " << nMocks << "\r"; cout.flush();
674 
675  auto dd_SS = Pair::Create(m_dd->pairType(), m_dd->pairInfo(), m_dd->sMin_D1(), m_dd->sMax_D1(), m_dd->nbins_D1(), m_dd->shift_D1(), m_dd->sMin_D2(), m_dd->sMax_D2(), m_dd->nbins_D2(), m_dd->shift_D2(), m_dd->angularUnits(), m_dd->angularWeight());
676  auto rr_SS = Pair::Create(m_rr->pairType(), m_rr->pairInfo(), m_rr->sMin_D1(), m_rr->sMax_D1(), m_rr->nbins_D1(), m_rr->shift_D1(), m_rr->sMin_D2(), m_rr->sMax_D2(), m_rr->nbins_D2(), m_rr->shift_D2(), m_rr->angularUnits(), m_rr->angularWeight());
677  auto dr_SS = Pair::Create(m_rr->pairType(), m_rr->pairInfo(), m_rr->sMin_D1(), m_rr->sMax_D1(), m_rr->nbins_D1(), m_rr->shift_D1(), m_rr->sMin_D2(), m_rr->sMax_D2(), m_rr->nbins_D2(), m_rr->shift_D2(), m_rr->angularUnits(), m_rr->angularWeight());
678 
679  double nData_SS = 0., nData_SS_weighted = 0., nRandom_SS = 0., nRandom_SS_weighted = 0.;
680 
681  vector<int> w(nRegions, 0);
682  for (size_t n=0; n<val*nRegions; n++)
683  w[ran()] ++;
684 
685  for (size_t j=0; j<nRegions; j++) {
686 
687  if (w[j]>0) {
688 
689  nData_SS += w[j]*nData_reg[j];
690  nData_SS_weighted += w[j]*nData_reg_weighted[j];
691  nRandom_SS += w[j]*nRandom_reg[j];
692  nRandom_SS_weighted += w[j]*nRandom_reg_weighted[j];
693 
694  for (size_t k=j; k<nRegions; k++) {
695  if(w[k]>0){
696  int index = j*nRegions-(j-1)*j/2+k-j;
697  double ww = w[j]*w[k];
698  for (int bin1=0; bin1<nbins_D1; bin1++)
699  for (int bin2=0; bin2<nbins_D2; bin2++) {
700  dd_SS->add_data2D(bin1, bin2, dd[index], ww);
701  rr_SS->add_data2D(bin1, bin2, rr[index], ww);
702  }
703  }
704  }
705  }
706 
707  for (size_t k=0; k<nRegions; k++) {
708  if (w[k]>0) {
709  int index = j*nRegions+k;
710  double ww = w[j]*w[k];
711  for (int bin1=0; bin1<nbins_D1; bin1++)
712  for (int bin2=0; bin2<nbins_D2; bin2++){
713  dr_SS->add_data2D(bin1, bin2, dr[index], ww);
714  }
715  }
716 
717  }
718  }
719  data.push_back(move(correlation_LandySzalayEstimator(dd_SS, rr_SS, dr_SS, nData_SS, nData_SS_weighted, nRandom_SS, nRandom_SS_weighted)));
720  }
721 
722  return data;
723 }
724 
725 
Data2D_extra.h
The class Data2D_extra.
coutCBL
#define coutCBL
CBL print message.
Definition: Kernel.h:734
TwoPointCorrelation2D.h
The class TwoPointCorrelation2D.
cbl::data::Data2D_extra
The class Data2D_extra.
Definition: Data2D_extra.h:52
cbl::data::Data2D
The class Data2D.
Definition: Data2D.h:51
cbl::measure::twopt::TwoPointCorrelation2D::XiBootstrap
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)
Definition: TwoPointCorrelation2D.cpp:580
cbl::measure::twopt::TwoPointCorrelation2D::correlation_LandySzalayEstimator
std::shared_ptr< data::Data > correlation_LandySzalayEstimator(const std::shared_ptr< pairs::Pair > dd, const std::shared_ptr< pairs::Pair > rr, const std::shared_ptr< pairs::Pair > dr, const int nData=0, const double nData_weighted=0., const int nRandom=0, const double nRandom_weighted=0.) override
get a dataset containing the two-point correlation function measured with the Landy-Szalay estimator,...
Definition: TwoPointCorrelation2D.cpp:406
cbl::measure::twopt::TwoPointCorrelation2D::XiJackknife
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)
Definition: TwoPointCorrelation2D.cpp:476
cbl::measure::twopt::TwoPointCorrelation2D::read_pairs
void read_pairs(std::shared_ptr< pairs::Pair > PP, const std::vector< std::string > dir, const std::string file) const override
read the number of pairs
Definition: TwoPointCorrelation2D.cpp:104
cbl::measure::twopt::TwoPointCorrelation2D::data_with_extra_info
std::shared_ptr< data::Data > data_with_extra_info(const std::shared_ptr< pairs::Pair > dd, const std::vector< double > scale_D1, const std::vector< double > scale_D2, const std::vector< std::vector< double >> xi, const std::vector< std::vector< double >> error) const
return a data object with extra info
Definition: TwoPointCorrelation2D.cpp:322
cbl::measure::twopt::TwoPointCorrelation2D::write_pairs
void write_pairs(const std::shared_ptr< pairs::Pair > PP, const std::string dir, const std::string file) const override
write the number of pairs
Definition: TwoPointCorrelation2D.cpp:54
cbl::measure::twopt::TwoPointCorrelation2D::correlation_NaturalEstimator
std::shared_ptr< data::Data > correlation_NaturalEstimator(const std::shared_ptr< pairs::Pair > dd, const std::shared_ptr< pairs::Pair > rr, const int nData=0, const double nData_weighted=0., const int nRandom=0, const double nRandom_weighted=0.) override
get a dataset containing the two-point correlation function measured with the natural estimator,...
Definition: TwoPointCorrelation2D.cpp:343
cbl::random::UniformRandomNumbers_Int
The class UniformRandomNumbers_Int.
Definition: RandomNumbers.h:296
cbl::par::fDP3
static const char fDP3[]
conversion symbol for: double -> std::string
Definition: Constants.h:136
cbl::par::fINT
static const char fINT[]
conversion symbol for: int -> std::string
Definition: Constants.h:121
cbl
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38
cbl::conv
std::string conv(const T val, const char *fact)
convert a number to a std::string
Definition: Kernel.h:903
cbl::checkIO
void checkIO(const std::ifstream &fin, const std::string file="NULL")
check if an input file can be opened
Definition: Kernel.cpp:160
cbl::ErrorCBL
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