CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
ThreePointCorrelation_comoving_reduced.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 
41 #include "Data1D.h"
42 
43 using namespace std;
44 
45 using namespace cbl;
46 using namespace catalogue;
47 using namespace measure;
48 using namespace triplets;
49 using namespace threept;
50 
51 
52 // ============================================================================
53 
54 
55 void cbl::measure::threept::ThreePointCorrelation_comoving_reduced::measure (const std::string dir_output_triplets, const std::string dir_output_2pt, const std::vector<std::string> dir_input_triplets, const bool count_ddd, const bool count_rrr, const bool count_ddr, const bool count_drr, const bool tcount, const double fact, const int seed)
56 {
57  (void)seed;
58 
59  // ----------- compute the connected three-point correlation function -----------
60 
61  ThreePointCorrelation_comoving_connected::measure(dir_output_triplets, dir_input_triplets, count_ddd, count_rrr, count_ddr, count_drr, tcount, fact);
62  m_scale = ThreePointCorrelation_comoving_connected::m_scale;
63 
64 
65  // ----------- compute the two-point correlation function -----------
66 
67  Catalogue data = *m_data;
68  Catalogue random = *m_random;
69  double rMin = m_ddd->r12();
70  double rMax = m_ddd->r12()+m_ddd->r13()+m_ddd->r13_binSize();
71  double binSize = 0.05;
72  double shift = 0.5;
73  twopt::TwoPointCorrelation1D_monopole TwoP {data, random, BinType::_logarithmic_, rMin, rMax, binSize, shift};
74  TwoP.measure(ErrorType::_Poisson_, dir_output_triplets, {}, par::defaultString, 0, 1, 1, 1, tcount, cbl::measure::twopt::Estimator::_LandySzalay_, fact);
75 
76  vector<double> log_r(TwoP.dd()->nbins()), log_xi(TwoP.dd()->nbins());
77  for (int i=0; i<TwoP.dd()->nbins(); i++) {
78  log_r[i] = log10(TwoP.xx()[i]);
79  log_xi[i] = log10(1.+TwoP.xi1D()[i]);
80  }
81 
82  vector<double> values_interp(m_ddd->nbins()+2, 0.), xi_real_lin(m_ddd->nbins()+2, 0.);
83 
84  values_interp[0] = log10(m_ddd->r12());
85  values_interp[1] = log10(m_ddd->r13());
86 
87  vector<double> theta(m_ddd->nbins(),0.);
88 
89  for (int i=0; i<m_ddd->nbins(); i++) {
90  theta[i]=(i+0.5)*m_ddd->binSize();
91  double tmp_value = sqrt(m_ddd->r12()*m_ddd->r12()+m_ddd->r13()*m_ddd->r13()-2.*m_ddd->r12()*m_ddd->r13()*cos(theta[i]));
92  //double tmp_value = (m_ddd->side_s()+((i+0.5)*m_ddd->binSize()));
93  values_interp[i+2] = log10(tmp_value);
94  }
95 
96  string file_2pt = dir_output_2pt+"2ptCorrelation_3pt.dat";
97  ofstream fout(file_2pt.c_str()); checkIO(fout, file_2pt);
98 
99  for (size_t i=0; i<values_interp.size(); i++) {
100  xi_real_lin[i] = pow(10., interpolated(values_interp[i], log_r, log_xi, "Linear"))-1.;
101  //coutCBL << pow(10.,values_interp[i]) << " --- " << xi_real_lin[i] << endl;
102  fout << pow(10.,values_interp[i]) << " " << xi_real_lin[i] << endl;
103  }
104 
105  fout.clear(); fout.close(); coutCBL <<"I wrote the file "<<file_2pt<<endl;
106 
107 
108  // ----------- compute the reduced three-point correlation function -----------
109 
110  m_QQ.resize(m_ddd->nbins()); m_error.resize(m_ddd->nbins());
111 
112  for (int i=0; i<m_ddd->nbins(); i++) {
113  m_QQ[i] = m_zeta[i]/((xi_real_lin[0]*xi_real_lin[1])+(xi_real_lin[0]*xi_real_lin[i+2])+(xi_real_lin[1]*xi_real_lin[i+2]));
114  m_error[i] = 0.001; // work in progress...
115  }
116 
117  m_dataset = move(unique_ptr<data::Data1D>(new data::Data1D(theta, m_QQ, m_error)));
118 }
119 
120 
121 // ============================================================================
122 
123 
124 void cbl::measure::threept::ThreePointCorrelation_comoving_reduced::measure (const std::vector<std::vector<double>> weight, const bool doJK, const std::string dir_output_triplets, const std::string dir_output_2pt, const std::vector<std::string> dir_input_triplets, const bool count_ddd, const bool count_rrr, const bool count_ddr, const bool count_drr, const bool tcount, const double fact, const int seed)
125 {
126  (void)seed;
127 
128  // ----------- compute the connected three-point correlation function -----------
129 
130  ThreePointCorrelation_comoving_connected::measure(weight, doJK, dir_output_triplets, dir_input_triplets, count_ddd, count_rrr, count_ddr, count_drr, tcount, fact);
131  m_scale = ThreePointCorrelation_comoving_connected::m_scale;
132 
133 
134  // ----------- compute the two-point correlation function -----------
135 
136  Catalogue data = *m_data;
137  Catalogue random = *m_random;
138  double rMin = m_ddd->r12();
139  double rMax = m_ddd->r12()+m_ddd->r13()+m_ddd->r13_binSize();
140  double binSize = 0.05;
141  double shift = 0.5;
142  twopt::TwoPointCorrelation1D_monopole TwoP {data, random, BinType::_logarithmic_, rMin, rMax, binSize, shift};
143  TwoP.measure(ErrorType::_Poisson_,dir_output_triplets, {}, par::defaultString, 0, 1, 1, 1, tcount, cbl::measure::twopt::Estimator::_LandySzalay_, fact);
144 
145  vector<double> log_r(TwoP.dd()->nbins()), log_xi(TwoP.dd()->nbins());
146  for (int i=0; i<TwoP.dd()->nbins(); i++) {
147  log_r[i] = log10(TwoP.xx()[i]);
148  log_xi[i] = log10(1.+TwoP.xi1D()[i]);
149  }
150 
151  vector<double> values_interp(m_ddd->nbins()+2, 0.), xi_real_lin(m_ddd->nbins()+2, 0.);
152 
153  values_interp[0] = log10(m_ddd->r12());
154  values_interp[1] = log10(m_ddd->r13());
155 
156  vector<double> theta(m_ddd->nbins(),0.);
157 
158  for (int i=0; i<m_ddd->nbins(); i++) {
159  theta[i]=(i+0.5)*m_ddd->binSize();
160  double tmp_value = sqrt(m_ddd->r12()*m_ddd->r12()+m_ddd->r13()*m_ddd->r13()-2.*m_ddd->r12()*m_ddd->r13()*cos(theta[i]));
161  //double tmp_value = (m_ddd->side_s()+((i+0.5)*m_ddd->binSize()));
162  values_interp[i+2] = log10(tmp_value);
163  }
164 
165  string file_2pt = dir_output_2pt+"2ptCorrelation_3pt.dat";
166  ofstream fout(file_2pt.c_str()); checkIO(fout, file_2pt);
167 
168  for (size_t i=0; i<values_interp.size(); i++) {
169  xi_real_lin[i] = pow(10., interpolated(values_interp[i], log_r, log_xi, "Linear"))-1.;
170  //coutCBL << pow(10.,values_interp[i]) << " --- " << xi_real_lin[i] << endl;
171  fout << pow(10.,values_interp[i]) << " " << xi_real_lin[i] << endl;
172  }
173 
174  fout.clear(); fout.close(); coutCBL <<"I wrote the file "<<file_2pt<<endl;
175 
176 
177  // ----------- compute the reduced three-point correlation function -----------
178 
179  m_QQ.resize(m_ddd->nbins()); m_error.resize(m_ddd->nbins());
180 
181  for (int i=0; i<m_ddd->nbins(); i++) {
182  m_QQ[i] = m_zeta[i]/((xi_real_lin[0]*xi_real_lin[1])+(xi_real_lin[0]*xi_real_lin[i+2])+(xi_real_lin[1]*xi_real_lin[i+2]));
183  m_error[i] = 0.001; // work in progress...
184  }
185 
187 
188  vector<vector<double>> resampling_threept(weight.size(), vector<double>(m_ddd->nbins(), 0));
189  vector<long> region_list = m_data->region_list();
190 
191  vector<double> nData_reg_weighted, nRandom_reg_weighted;
192 
193  const int nRegions = m_data->nRegions();
194 
195  for (int i=0; i<nRegions; i++) {
196  nData_reg_weighted.push_back(m_data->weightedN_condition(Var::_Region_, region_list[i], region_list[i]+1, 0));
197  nRandom_reg_weighted.push_back(m_random->weightedN_condition(Var::_Region_, region_list[i], region_list[i]+1, 0));
198  }
199 
200  for (size_t i=0; i<weight.size(); i++) {
201 
202  double nData = 0;
203  double nRan = 0;
204 
205  for (size_t j=0; j<weight[i].size(); j++) {
206  nData += weight[i][j]*nData_reg_weighted[j];
207  nRan += weight[i][j]*nRandom_reg_weighted[j];
208  }
209 
210  double norm1 = (double(nData)*double(nData-1)*double(nData-2))/6.;
211  double norm2 = (double(nData)*double(nData-1)*double(nRan))*0.5;
212  double norm3 = (double(nData)*double(nRan)*double(nRan-1))*0.5;
213  double norm4 = (double(nRan)*double(nRan-1)*double(nRan-2))/6.;
214 
215  for (int j=0; j<m_ddd->nbins(); j++)
216  if (m_ddd_regions[i]->TT1D(j)>0 && m_rrr_regions[i]->TT1D(j)>0) {
217  resampling_threept[i][j] = ((m_ddd_regions[i]->TT1D(j)/norm1)/(m_rrr_regions[i]->TT1D(j)/norm4))-3.*((m_ddr_regions[i]->TT1D(j)/norm2)/(m_rrr_regions[i]->TT1D(j)/norm4))+3.*(((m_drr_regions[i]->TT1D(j)/norm3)/(m_rrr_regions[i]->TT1D(j)/norm4)))-1.;
218  resampling_threept[i][j] /= ((xi_real_lin[0]*xi_real_lin[1])+(xi_real_lin[0]*xi_real_lin[j+2])+(xi_real_lin[1]*xi_real_lin[j+2]));
219  }
220  }
221 
222  vector<vector<double>> cov_mat;
223  cbl::covariance_matrix(resampling_threept, cov_mat, doJK);
224  m_dataset = move(unique_ptr<data::Data1D>(new data::Data1D(m_scale, m_zeta, cov_mat)));
225 
226  m_error = m_dataset->error();
227 }
228 
229 
230 // ============================================================================================
231 
232 
233 
234 void cbl::measure::threept::ThreePointCorrelation_comoving_reduced::measure (const cbl::measure::ErrorType errorType, const std::string dir_output_triplets, const std::string dir_output_2pt, const std::vector<std::string> dir_input_triplets, const int nResamplings, const bool count_ddd, const bool count_rrr, const bool count_ddr, const bool count_drr, const bool tcount, const double fact, const int seed)
235 {
236 
237  switch (errorType) {
238 
240  {
241  measure(dir_output_triplets, dir_output_2pt, dir_input_triplets, count_ddd, count_rrr, count_ddr, count_drr, tcount, fact);
242  break;
243  }
244 
246  {
247  const int nRegions = m_data->nRegions();
248 
249  vector<vector<double>> weight(nRegions, vector<double>(nRegions, 1));
250  for (int i=0; i<nRegions; i++)
251  weight[i][i] = 0;
252 
253  measure(weight, true, dir_output_triplets, dir_output_2pt, dir_input_triplets, count_ddd, count_rrr, count_ddr, count_drr, tcount, fact);
254  break;
255  }
256 
258  {
259  const int nRegions = m_data->nRegions();
260 
261  random::UniformRandomNumbers_Int ran(0., nRegions-1, seed);
262 
263  int val = 3; // see Norberg et al. 2009
264 
265  vector<vector<double>> weight(nResamplings, vector<double>(nRegions, 0));
266  for (int i=0; i<nResamplings; i++)
267  for (int j=0; j<val*nRegions; j++)
268  weight[i][ran()] ++;
269 
270  measure(weight, false, dir_output_triplets, dir_output_2pt, dir_input_triplets, count_ddd, count_rrr, count_ddr, count_drr, tcount, fact);
271  break;
272  }
273 
274  default:
275  ErrorCBL("no such kind of error type!", "measure", "hreePointCorrelation_comoving_reduced.cpp");
276  }
277 
278 }
279 
280 
281 // ============================================================================
282 
283 
284 void cbl::measure::threept::ThreePointCorrelation_comoving_reduced::write (const std::string dir, const std::string file, bool connected) const
285 {
286  checkDim(m_scale, m_ddd->nbins(), "scale");
287 
288  string file_out = dir+file;
289  ofstream fout(file_out.c_str()); checkIO(fout, file_out);
290 
291  if (!connected) {
292  fout << "# scale Q error(work in progress)" << endl;
293  for (size_t i=0; i<m_scale.size(); i++)
294  fout << setiosflags(ios::fixed) << setprecision(4) << setw(10) << right << m_scale[i]
295  << " " << setiosflags(ios::fixed) << setprecision(4) << setw(10) << right << m_QQ[i]
296  << " " << setiosflags(ios::fixed) << setprecision(4) << setw(10) << right << m_error[i] << endl;
297  }
298 
299  else {
300  fout << "# scale z error(work in progrss) Q error(work in progress)" << endl;
301  for (size_t i=0; i<m_scale.size(); i++)
302  fout << setiosflags(ios::fixed) << setprecision(4) << setw(10) << right << m_scale[i]
303  << " " << setiosflags(ios::fixed) << setprecision(4) << setw(10) << right << m_zeta[i]
304  << " " << setiosflags(ios::fixed) << setprecision(4) << setw(10) << right << ThreePointCorrelation_comoving_connected::m_error[i]
305  << " " << setiosflags(ios::fixed) << setprecision(4) << setw(10) << right << m_QQ[i]
306  << " " << setiosflags(ios::fixed) << setprecision(4) << setw(10) << right << m_error[i] << endl;
307  }
308 
309  fout.close(); coutCBL << endl << "I wrote the file: " << file_out << endl << endl;
310 }
311 
312 
313 // ============================================================================
314 
315 
316 void cbl::measure::threept::ThreePointCorrelation_comoving_reduced::write_covariance (const std::string dir, const std::string file) const
317 {
318  m_dataset->write_covariance(dir, file);
319 }
The class Data1D.
#define coutCBL
CBL print message.
Definition: Kernel.h:734
The class ThreePointCorrelation_comoving_reduced.
The class TwoPointCorrelation1D_monopole.
The class Catalogue.
Definition: Catalogue.h:654
The class Data1D.
Definition: Data1D.h:51
void measure(const std::string dir_output_triplets, const std::string dir_output_2pt, const std::vector< std::string > dir_input_triplets={}, const bool count_ddd=true, const bool count_rrr=true, const bool count_ddr=true, const bool count_drr=true, const bool tcount=false, const double fact=0.1, const int seed=3213) override
method to measure the three-point correlation function
void write_covariance(const std::string dir, const std::string file) const override
write the measured three-point correlation covariance
void write(const std::string dir, const std::string file, const bool connected) const override
write the monopole of the two-point correlation function
The class UniformRandomNumbers_Int.
static const std::string defaultString
default std::string value
Definition: Constants.h:336
@ _LandySzalay_
Landy&Szalay estimator.
ErrorType
the two-point correlation function error type
Definition: Measure.h:59
@ _Bootstrap_
Bootstrap resampling.
@ _None_
No error computed.
@ _Jackknife_
Jackknife resampling.
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38
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
void checkIO(const std::ifstream &fin, const std::string file="NULL")
check if an input file can be opened
Definition: Kernel.cpp:160
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
double interpolated(const double _xx, const std::vector< double > xx, const std::vector< double > yy, const std::string type)
1D interpolation
Definition: Func.cpp:445
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