CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
TwoPointCorrelation1D_angular.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 
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::TwoPointCorrelation1D_angular::set_parameters (const BinType binType, const double thetaMin, const double thetaMax, const int nbins, const double shift, const CoordinateUnits angularUnits, std::function<double(double)> angularWeight, const bool compute_extra_info)
55 {
56  if (!compute_extra_info)
57  m_dd = (binType==BinType::_logarithmic_) ? move(Pair::Create(PairType::_angular_log_, PairInfo::_standard_, thetaMin, thetaMax, nbins, shift, angularUnits, angularWeight))
58  : move(Pair::Create(PairType::_angular_lin_, PairInfo::_standard_, thetaMin, thetaMax, nbins, shift, angularUnits, angularWeight));
59  else
60  m_dd = (binType==BinType::_logarithmic_) ? move(Pair::Create(PairType::_angular_log_, PairInfo::_extra_, thetaMin, thetaMax, nbins, shift, angularUnits, angularWeight))
61  : move(Pair::Create(PairType::_angular_lin_, PairInfo::_extra_, thetaMin, thetaMax, nbins, shift, angularUnits, angularWeight));
62 
63  m_rr = (binType==BinType::_logarithmic_) ? move(Pair::Create(PairType::_angular_log_, PairInfo::_standard_, thetaMin, thetaMax, nbins, shift, angularUnits))
64  : move(Pair::Create(PairType::_angular_lin_, PairInfo::_extra_, thetaMin, thetaMax, nbins, shift, angularUnits));
65 
66  m_dr = (binType==BinType::_logarithmic_) ? move(Pair::Create(PairType::_angular_log_, PairInfo::_standard_, thetaMin, thetaMax, nbins, shift, angularUnits))
67  : move(Pair::Create(PairType::_angular_lin_, PairInfo::_extra_, thetaMin, thetaMax, nbins, shift, angularUnits));
68 }
69 
70 
71 // ============================================================================================
72 
73 
74 void cbl::measure::twopt::TwoPointCorrelation1D_angular::set_parameters (const BinType binType, const double thetaMin, const double thetaMax, const double binSize, const double shift, const CoordinateUnits angularUnits, std::function<double(double)> angularWeight, const bool compute_extra_info)
75 {
76  if (!compute_extra_info)
77  m_dd = (binType==BinType::_logarithmic_) ? move(Pair::Create(PairType::_angular_log_, PairInfo::_standard_, thetaMin, thetaMax, binSize, shift, angularUnits, angularWeight))
78  : move(Pair::Create(PairType::_angular_lin_, PairInfo::_standard_, thetaMin, thetaMax, binSize, shift, angularUnits, angularWeight));
79  else
80  m_dd = (binType==BinType::_logarithmic_) ? move(Pair::Create(PairType::_angular_log_, PairInfo::_extra_, thetaMin, thetaMax, binSize, shift, angularUnits, angularWeight))
81  : move(Pair::Create(PairType::_angular_lin_, PairInfo::_extra_, thetaMin, thetaMax, binSize, shift, angularUnits, angularWeight));
82 
83  m_rr = (binType==BinType::_logarithmic_) ? move(Pair::Create(PairType::_angular_log_, PairInfo::_standard_, thetaMin, thetaMax, binSize, shift, angularUnits))
84  : move(Pair::Create(PairType::_angular_lin_, PairInfo::_standard_, thetaMin, thetaMax, binSize, shift, angularUnits));
85 
86  m_dr = (binType==BinType::_logarithmic_) ? move(Pair::Create(PairType::_angular_log_, PairInfo::_standard_, thetaMin, thetaMax, binSize, shift, angularUnits))
87  : move(Pair::Create(PairType::_angular_lin_, PairInfo::_standard_, thetaMin, thetaMax, binSize, shift, angularUnits));
88 }
89 
90 
91 // ============================================================================================
92 
93 
94 void cbl::measure::twopt::TwoPointCorrelation1D_angular::read (const string dir, const string file)
95 {
96  m_dataset->read(dir+file);
97 }
98 
99 
100 // ============================================================================================
101 
102 
103 void cbl::measure::twopt::TwoPointCorrelation1D_angular::write (const string dir, const string file, const int rank) const
104 {
105  vector<double> xx = m_dataset->xx();
106 
107  checkDim(xx, m_dd->nbins(), "theta");
108 
109  string header = "[1] angular separation at the bin centre # [2] angular two-point correlation function # [3] error";
110  if (m_compute_extra_info) header += " # [4] mean angular separation # [5] standard deviation of the distribution of angular separations # [6] mean redshift # [7] standard deviation of the redshift distribution";
111 
112  m_dataset->write(dir, file, header, 5, rank);
113 }
114 
115 
116 // ============================================================================================
117 
118 
119 void cbl::measure::twopt::TwoPointCorrelation1D_angular::measure (const ErrorType errorType, const string dir_output_pairs, const vector<string> dir_input_pairs, const string dir_output_resample, int nMocks, bool count_dd, const bool count_rr, const bool count_dr, const bool tcount, const Estimator estimator, const double fact, const int seed)
120 {
121  switch (errorType) {
122  case(ErrorType::_Poisson_):
123  measurePoisson(dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
124  break;
125  case(ErrorType::_Jackknife_):
126  measureJackknife(dir_output_pairs, dir_input_pairs, dir_output_resample, count_dd, count_rr, count_dr, tcount, estimator, fact);
127  break;
128  case(ErrorType::_Bootstrap_):
129  measureBootstrap(nMocks, dir_output_pairs, dir_input_pairs, dir_output_resample, count_dd, count_rr, count_dr, tcount, estimator, fact, seed);
130  break;
131  default:
132  ErrorCBL("unknown type of error!", "measure", "TwoPointCorrelation1D_angular.cpp");
133  }
134 }
135 
136 
137 // ============================================================================================
138 
139 
140 void cbl::measure::twopt::TwoPointCorrelation1D_angular::measurePoisson (const string dir_output_pairs, const vector<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)
141 {
142  // ----------- count the data-data, random-random and data-random pairs, or read them from file -----------
143 
144  count_allPairs(m_twoPType, dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
145 
146 
147  // ----------- compute the angular two-point correlation function -----------
148 
149  if (estimator==Estimator::_natural_)
150  m_dataset = correlation_NaturalEstimator(m_dd, m_rr);
151  else if (estimator==Estimator::_LandySzalay_)
152  m_dataset = correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
153  else
154  ErrorCBL("the chosen estimator is not implemented!", "measurePoisson", "TwoPointCorrelation1D_angular.cpp");
155 }
156 
157 
158 // ============================================================================================
159 
160 
161 void cbl::measure::twopt::TwoPointCorrelation1D_angular::measureJackknife (const string dir_output_pairs, const vector<string> dir_input_pairs, const 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)
162 {
163  if (dir_output_resample!=par::defaultString && dir_output_resample!="") {
164  string mkdir = "mkdir -p "+dir_output_resample;
165  if (system(mkdir.c_str())) {}
166  }
167 
168  vector<long> region_list = m_data->region_list();
169  size_t nRegions = region_list.size();
170 
171  vector<vector<double> > xi_SubSamples, covariance;
172 
173  vector<shared_ptr<Pair> > dd_regions, rr_regions, dr_regions;
174 
175  count_allPairs_region(dd_regions, rr_regions, dr_regions, m_twoPType, dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
176 
177  vector<shared_ptr<Data> > data_SS = (estimator==Estimator::_natural_) ? XiJackknife(dd_regions, rr_regions) : XiJackknife(dd_regions, rr_regions, dr_regions);
178 
179  for (size_t i=0; i<nRegions; i++) {
180 
181  if (dir_output_resample!=par::defaultString && dir_output_resample!="") {
182  string file = "xi_Jackknife_"+conv(i, par::fINT)+".dat";
183  string header = "[1] angular separation at the bin centre # [2] angular two-point correlation function # [3] error";
184  if (m_compute_extra_info) header += " # [4] mean separation # [5] standard deviation of the separation distribution # [6] mean redshift # [7] standard deviation of the redshift distribution";
185  data_SS[i]->write(dir_output_resample, file, header, 10, 0);
186  }
187 
188  vector<double> dd; data_SS[i]->get_data(dd);
189 
190  xi_SubSamples.push_back(dd);
191  }
192 
193  covariance_matrix(xi_SubSamples, covariance, 1);
194 
195  if (estimator==Estimator::_natural_)
196  m_dataset = correlation_NaturalEstimator(m_dd, m_rr);
197  else if (estimator==Estimator::_LandySzalay_)
198  m_dataset = correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
199  else
200  ErrorCBL("the chosen estimator is not implemented!", "measureJackknife", "TwoPointCorrelation1D_angular.cpp");
201 
202  m_dataset->set_covariance(covariance);
203 
204 }
205 
206 
207 // ============================================================================================
208 
209 
210 void cbl::measure::twopt::TwoPointCorrelation1D_angular::measureBootstrap (const int nMocks, const string dir_output_pairs, const vector<string> dir_input_pairs, const 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)
211 {
212  if (nMocks<=0)
213  ErrorCBL("the number of mocks must be >0!", "measureBootstrap", "TwoPointCorrelation1D_angular.cpp");
214 
215  if (dir_output_resample!=par::defaultString && dir_output_resample!="") {
216  string mkdir = "mkdir -p "+dir_output_resample;
217  if (system(mkdir.c_str())) {}
218  }
219 
220  vector<vector<double> > xi_SubSamples, covariance;
221 
222  vector<shared_ptr<Pair> > dd_regions, rr_regions, dr_regions;
223 
224  count_allPairs_region(dd_regions, rr_regions, dr_regions, m_twoPType, dir_output_pairs, dir_input_pairs, count_dd, count_rr, count_dr, tcount, estimator, fact);
225 
226  vector<shared_ptr<Data> > data_SS = (estimator==Estimator::_natural_) ? XiBootstrap(nMocks, dd_regions, rr_regions, seed) : XiBootstrap(nMocks, dd_regions, rr_regions, dr_regions, seed);
227 
228  for (int i=0; i<nMocks; i++) {
229 
230  if (dir_output_resample!=par::defaultString && dir_output_resample!="") {
231  string file = "xi_Bootstrap_"+conv(i, par::fINT)+".dat";
232  string header = "[1] angular separation at the bin centre # [2] angular two-point correlation function # [3] error";
233  if (m_compute_extra_info) header += " # [4] mean separation # [5] standard deviation of the separation distribution # [6] mean redshift # [7] standard deviation of the redshift distribution";
234  data_SS[i]->write(dir_output_resample, file, header, 10, 0);
235  }
236 
237  vector<double> dd; data_SS[i]->get_data(dd);
238 
239  xi_SubSamples.push_back(dd);
240  }
241 
242  covariance_matrix(xi_SubSamples, covariance, 0);
243 
244  if (estimator==Estimator::_natural_)
245  m_dataset = correlation_NaturalEstimator(m_dd, m_rr);
246  else if (estimator==Estimator::_LandySzalay_)
247  m_dataset = correlation_LandySzalayEstimator(m_dd, m_rr, m_dr);
248  else
249  ErrorCBL("the chosen estimator is not implemented!", "measureBootstrap", "TwoPointCorrelation1D_angular.cpp");
250 
251  m_dataset->set_covariance(covariance);
252 
253 }
The class TwoPointCorrelation1D_angular.
void measureBootstrap(const int nMocks, 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, const int seed=3213) override
measure the angular two-point correlation function estimating the covariance with Bootstrap resamplin...
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 angular two-point correlation function with Poisson errors
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 angular two-point correlation function estimating the covariance with Jackknife resamplin...
void read(const std::string dir, const std::string file) override
read the angular two-point 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 angular two-point correlation function
void set_parameters(const BinType binType, const double thetaMin, const double thetaMax, const int nbins, const double shift, const CoordinateUnits angularUnits=CoordinateUnits::_radians_, std::function< double(double)> angularWeight=nullptr, const bool compute_extra_info=false)
set the binning parameters
void write(const std::string dir=par::defaultString, const std::string file=par::defaultString, const int rank=0) const override
write the angular two-point correlation function
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
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
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
CoordinateUnits
the coordinate units
Definition: Kernel.h:562
BinType
the binning type
Definition: Kernel.h:505