CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
VoidCatalogue.cpp
Go to the documentation of this file.
1 /********************************************************************
2  * Copyright (C) 2022 by Federico Marulli, Simone Sartori *
3  * Sofia Contarini, Carlo Cannarozzo and Tommaso Ronconi *
4  * federico.marulli3@unibo.it simone.sartori5@studio.unibo.it *
5  * sofia.contarini3@unibo.it carlo.cannarozzo@studio.unibo.it *
6  * tommaso.ronconi@studio.unibo.it *
7  * *
8  * This program is free software; you can redistribute it and/or *
9  * modify it under the terms of the GNU General Public License as *
10  * published by the Free Software Foundation; either version 2 of *
11  * the License, or (at your option) any later version. *
12  * *
13  * This program is distributed in the hope that it will be useful, *
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
16  * GNU General Public License for more details. *
17  * *
18  * You should have received a copy of the GNU General Public *
19  * License along with this program; if not, write to the Free *
20  * Software Foundation, Inc., *
21  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
22  ********************************************************************/
23 
39 #include "Func.h"
40 #include "Catalogue.h"
41 #include "ChainMesh_Catalogue.h"
42 #include "Object.h"
43 #include "Catalogue.h"
44 #include "Data1D.h"
45 #include "Posterior.h"
46 #include "TwoPointCorrelation.h"
47 #include "TwoPointCorrelation1D_monopole.h"
48 #include "CatalogueChainMesh.h"
49 #include <omp.h>
50 #include <chrono>
51 
52 using namespace std;
53 using namespace cbl;
54 using namespace catalogue;
55 
56 cbl::catalogue::Catalogue::Catalogue (const VoidAlgorithm algorithm, Catalogue tracer_catalogue, Catalogue random_catalogue, const std::string dir_output, const std::string output, const double cellsize, const int n_rec, const double step_size, const double threshold, const std::vector<bool> print)
57 {
58  // -------------------------------------------- //
59  // ---------------- First Step ---------------- //
60  // -------------------------------------------- //
61 
62  cout << endl; coutCBL << "> > > > > >>>>>>>>>> " << par::col_blue << "BACK-IN-TIME VOID FINDER" << par::col_default << " <<<<<<<<<< < < < < <" << endl << endl;
63 
64  const double start_time = omp_get_wtime();
65  shared_ptr<Catalogue> tracer_cat, random_cat;
66  vector<Catalogue> displacement_catalogue(n_rec);
67  vector<vector<unsigned int>> randoms;
68  cosmology::Cosmology cosm;
69  CatalogueChainMesh ChainMesh_tracers;
70  CatalogueChainMesh ChainMesh_randoms;
71 
72  if (algorithm==VoidAlgorithm::_LaZeVo_) {
73 
74  coutCBL << par::col_green << " LaZeVo reconstruction" << par::col_default << endl << endl;
75 
76  Catalogue random_catalogue;
77 
78  int rndd = (int)chrono::steady_clock::now().time_since_epoch().count()%1000000;
79 
80  if (random_catalogue.nObjects()==0) {
81  random_catalogue = Catalogue(RandomType::_createRandom_box_, tracer_catalogue, 1., 10, cosm, false, 10., {}, {}, {}, 10, rndd);
82  cout << endl;
83  }
84  else if (random_catalogue.nObjects() == tracer_catalogue.nObjects()) {
85  coutCBL << "* * * Reading random catalogue of " << random_catalogue.nObjects() << " objects provides by user * * *" << endl << endl;
86  random_catalogue = random_catalogue;
87  }
88  else {
89  coutCBL << "* * * Reading random catalogue of " << random_catalogue.nObjects() << " objects provides by user * * *" << endl;
90  cout << endl;
91  int diff = random_catalogue.nObjects() - tracer_catalogue.nObjects();
92  if (random_catalogue.nObjects() > tracer_catalogue.nObjects()) coutCBL << "WARNING! Random catalogue has " <<
93  diff << " more objects than the tracers catalogue. The number will be equalized by the function catalogue::equalize_random_LC" << endl << endl;
94  if (random_catalogue.nObjects() < tracer_catalogue.nObjects()) coutCBL << "WARNING! Random catalogue has " <<
95  -diff << " less objects than the tracers catalogue. The number will be equalized by the function catalogue::equalize_random_LC" << endl << endl;
96  random_catalogue.equalize_random_box(tracer_catalogue, rndd);
97  }
98 
99  random_cat = make_shared<cbl::catalogue::Catalogue> (cbl::catalogue::Catalogue(move(random_catalogue)));
100  tracer_cat = make_shared<cbl::catalogue::Catalogue> (cbl::catalogue::Catalogue(move(tracer_catalogue)));
101 
102  unsigned int num_objects = tracer_cat->nObjects();
103  std::vector<unsigned int> index_tracer_cat(num_objects);
104  std::iota (std::begin(index_tracer_cat), std::end(index_tracer_cat), 0);
105 
106  coutCBL << "* * * ChainMesh setting * * *" << flush;
107  vector<Var> variables = {Var::_X_, Var::_Y_, Var::_Z_};
108  ChainMesh_tracers = CatalogueChainMesh(variables, cellsize, tracer_cat);
109  ChainMesh_randoms = CatalogueChainMesh(variables, cellsize, random_cat);
110  cout << endl << endl;
111 
112  unsigned int N_near_obj = 3*3*4*par::pi;
113  coutCBL << "* * * Searching close particles * * * " << flush;
114  vector<vector<unsigned int>> near_part = ChainMesh_tracers.N_nearest_objects_cat(N_near_obj);
115  randoms.resize(n_rec, vector<unsigned int>(num_objects));
116  for(int rec=0; rec<n_rec; rec++) {
117  std::iota (std::begin(randoms[rec]), std::end(randoms[rec]), 0);
118  auto time = (long int)chrono::steady_clock::now().time_since_epoch().count()%1000000;
119  std::shuffle(randoms[rec].begin(), randoms[rec].end(), default_random_engine(time));
120  }
121 
122  cout << endl << endl;
123  coutCBL << "* * * Setting starting configuration * * *" << flush;
124 
125  double dist = 4*tracer_cat->mps();
126  double minX = ChainMesh_tracers.lim()[0][0]+dist/10, maxX = ChainMesh_tracers.lim()[0][1]-dist/10;
127  double minY = ChainMesh_tracers.lim()[1][0]+dist/10, maxY = ChainMesh_tracers.lim()[1][1]-dist/10;
128  double minZ = ChainMesh_tracers.lim()[2][0]+dist/10, maxZ = ChainMesh_tracers.lim()[2][1]-dist/10;
129 
130  for (int rec=0; rec<n_rec; rec++) {
131 
132  auto timeS = ((long int)chrono::steady_clock::now().time_since_epoch().count())%1000000;
133  random::UniformRandomNumbers randX(minX, maxX, timeS);
134  random::UniformRandomNumbers randY(minY, maxY, timeS+1);
135  random::UniformRandomNumbers randZ(minZ, maxZ, timeS+2);
136 
137  CatalogueChainMesh ChM_tracer_copy = ChainMesh_tracers;
138  CatalogueChainMesh ChM_random_copy = ChainMesh_randoms;
139  unsigned int removed = 0;
140 
141  while (removed < num_objects) {
142  vector<double> pos = {randX(), randY(), randZ()};
143  vector<unsigned int> close=ChM_tracer_copy.Close_objects(pos, dist);
144  unsigned int tr_to_rmv = std::min(100, (int)close.size());
145  if (tr_to_rmv > 0) {
146  vector<unsigned int> close_random = ChM_random_copy.N_nearest_objects(pos, close.size());
147  std::random_device dev;
148  std::mt19937 rng(dev());
149  removed += tr_to_rmv;
150  while (tr_to_rmv > 0) {
151  std::uniform_int_distribution<std::mt19937::result_type> dist(0, tr_to_rmv-1);
152  int rnd1 = dist(rng), rnd2 = dist(rng);
153  randoms[rec][close[rnd1]] = close_random[rnd2];
154  ChM_tracer_copy.deletePart(close[rnd1]);
155  ChM_random_copy.deletePart(close_random[rnd2]);
156  close.erase(close.begin()+rnd1);
157  close_random.erase(close_random.begin()+rnd2);
158  tr_to_rmv--;
159  }
160  }
161  }
162  }
163 
164  cout << endl << endl;
165 
166  coutCBL << "* * * Performing the iterations * * *" << endl << endl;
167 
168  for (int rec=0; rec<n_rec; rec++) {
169  auto index_tracer_cat_copy = index_tracer_cat;
170  double ratio = 1.;
171  unsigned int n_iter = 0;
172 
173  while (ratio > threshold) {
174 
175  auto time = (long int)chrono::steady_clock::now().time_since_epoch().count()%1000000;
176  std::shuffle(index_tracer_cat_copy.begin(), index_tracer_cat_copy.end(), default_random_engine(time));
177  std::random_device dev;
178  std::mt19937 rng(dev());
179  std::uniform_int_distribution<std::mt19937::result_type> dist(1,N_near_obj-1);
180  vector<bool> index_bool(num_objects, false), used(num_objects, false);
181 
182 #pragma omp parallel num_threads(omp_get_max_threads())
183  {
184 #pragma omp for schedule(static)
185  for (size_t i=0; i<num_objects; i++) {
186  vector<unsigned int> H(4), R_def(4), R(4), R_copy(4);
187  unsigned int rand1 = dist(rng), rand2 = dist(rng), rand3 = dist(rng);
188  while (rand1 == rand2 || rand1==rand3 || rand2 == rand3) {
189  rand2 = dist(rng);
190  rand3 = dist(rng);
191  }
192  if (used[near_part[index_tracer_cat_copy[i]][0]] == false && used[near_part[index_tracer_cat_copy[i]][rand1]] == false &&
193  used[near_part[index_tracer_cat_copy[i]][rand2]] == false && used[near_part[index_tracer_cat_copy[i]][rand3]] == false) {
194  used[near_part[index_tracer_cat_copy[i]][0]] = true;
195  used[near_part[index_tracer_cat_copy[i]][rand1]] = true;
196  used[near_part[index_tracer_cat_copy[i]][rand2]] = true;
197  used[near_part[index_tracer_cat_copy[i]][rand3]] = true;
198  H[0]=near_part[index_tracer_cat_copy[i]][0];
199  H[1]=near_part[index_tracer_cat_copy[i]][rand1];
200  H[2]=near_part[index_tracer_cat_copy[i]][rand2];
201  H[3]=near_part[index_tracer_cat_copy[i]][rand3];
202  R[0]=randoms[rec][H[0]];
203  R[1]=randoms[rec][H[1]];
204  R[2]=randoms[rec][H[2]];
205  R[3]=randoms[rec][H[3]];
206 
207  R_copy = R;
208  double dist_min = 1000000.;
209  double dist;
210  do {
211  dist = 0.;
212  for (size_t j=0; j<R.size(); j++)
213  dist += cbl::Euclidean_distance(tracer_cat->xx(H[j]), random_cat->xx(R[j]), tracer_cat->yy(H[j]), random_cat->yy(R[j]), tracer_cat->zz(H[j]), random_cat->zz(R[j]));
214 
215  if (dist < dist_min) {
216  dist_min = dist;
217  R_def = R;
218  }
219  } while(std::next_permutation(R.begin(), R.end()));
220 
221  if(R_def != R_copy) index_bool[i] = true;
222  for (size_t j=0; j<4; j++) randoms[rec][H[j]] = R_def[j];
223  used[H[0]] = false;
224  used[H[1]] = false;
225  used[H[2]] = false;
226  used[H[3]] = false;
227  }
228  else continue;
229  }
230  }
231 
232  unsigned int changed_couples = count(index_bool.begin(), index_bool.end(), true);
233  ratio = (double)changed_couples/num_objects;
234  if (n_iter%1==0) coutCBL << "Realisation " << rec+1 << " of " << n_rec << " | " << "Iteration N°: " << n_iter << " - " << "current ratio: " << fixed << setprecision(6) << ratio << " - " <<
235  "final threshold: " << threshold << "\r"; cout.flush();
236  n_iter++;
237  }
238 
240 
241  vector<shared_ptr<cbl::catalogue::Object>> displ_catalogue_object;
242  for (size_t ii = 0; ii<randoms[rec].size(); ii++) {
243  cbl::comovingCoordinates displacement_values = {random_cat->xx(randoms[rec][ii])-tracer_cat->xx(ii), random_cat->yy(randoms[rec][ii])-tracer_cat->yy(ii), random_cat->zz(randoms[rec][ii])-tracer_cat->zz(ii)};
244  auto tracer_displ = make_shared<cbl::catalogue::Halo>(displacement_values);
245  displ_catalogue_object.push_back(tracer_displ);
246  }
247 
248  displacement_catalogue[rec].add_objects(displ_catalogue_object);
249 
250  if (rec == n_rec-1 && print[0]) {
251  string name = "../output/Displacement";
252  name.append("_");
253  name.append(output);
254  ofstream fout_displ(name.c_str());
255  cout << endl << endl;
256 
257  coutCBL << "Writing the displacement field in " << name << " ..." << endl;
258 
259  for (unsigned int i=0; i<num_objects; i++) {
260  fout_displ.precision(7);
261  fout_displ << tracer_cat->xx(i) << " " << tracer_cat->yy(i) << " " << tracer_cat->zz(i) << " "
262  << random_cat->xx(randoms[rec][i]) << " " << random_cat->yy(randoms[rec][i]) << " " << random_cat->zz(randoms[rec][i]) << " "
263  << displacement_catalogue[rec].xx(i) << " " << displacement_catalogue[rec].yy(i) << " " << displacement_catalogue[rec].zz(i) << endl;
264  }
265  fout_displ.clear(); fout_displ.close();
266  }
267  }
268  } // End of LaZeVo method
269  else if (algorithm==VoidAlgorithm::_Exact_) {
270 
271  coutCBL << par::col_green << " Exact reconstruction" << par::col_default << endl << endl;
272 
273  randoms.resize(n_rec, vector<unsigned int>(tracer_catalogue.nObjects()));
274  for(int rec=0; rec<n_rec; rec++)
275  std::iota (std::begin(randoms[rec]), std::end(randoms[rec]), 0);
276 
277  if (random_catalogue.nObjects()!=tracer_catalogue.nObjects()) {
278  ErrorCBL ("Random catalogue must have the same dimension of the tracers one!", "Catalogue", "VoidCatalogue.cpp");
279  }
280  random_cat = make_shared<cbl::catalogue::Catalogue> (cbl::catalogue::Catalogue(move(random_catalogue)));
281  tracer_cat = make_shared<cbl::catalogue::Catalogue> (cbl::catalogue::Catalogue(move(tracer_catalogue)));
282  vector<Var> variables = {Var::_X_, Var::_Y_, Var::_Z_};
283  cout << endl;
284  coutCBL << "* * * ChainMesh setting * * *" << endl;
285  ChainMesh_tracers = CatalogueChainMesh(variables, cellsize, tracer_cat);
286 
287  for (int rec=0; rec<n_rec; rec++) {
288  vector<shared_ptr<cbl::catalogue::Object>> displ_catalogue_object;
289  for (size_t ii = 0; ii<randoms[rec].size(); ii++) {
290  cbl::comovingCoordinates displacement_values = {random_cat->xx(randoms[rec][ii])-tracer_cat->xx(ii), random_cat->yy(randoms[rec][ii])-tracer_cat->yy(ii), random_cat->zz(randoms[rec][ii])-tracer_cat->zz(ii)};
291  auto tracer_displ = make_shared<cbl::catalogue::Halo>(displacement_values);
292  displ_catalogue_object.push_back(tracer_displ);
293  }
294  displacement_catalogue[rec].add_objects(displ_catalogue_object);
295 
296  if (rec == n_rec-1 && print[0]) {
297  string name = dir_output + "Displacement";
298  name.append("_");
299  name.append(output);
300  ofstream fout_displ(name.c_str());
301  cout << endl;
302 
303  coutCBL << "Writing the displacement field in " << name << " ..." << endl;
304  for (size_t i=0; i<tracer_cat->nObjects(); i++) {
305  fout_displ.precision(7);
306  fout_displ << tracer_cat->xx(i) << " " << tracer_cat->yy(i) << " " << tracer_cat->zz(i) << " "
307  << random_cat->xx(randoms[rec][i]) << " " << random_cat->yy(randoms[rec][i]) << " " << random_cat->zz(randoms[rec][i]) << " "
308  << displacement_catalogue[rec].xx(i) << " " << displacement_catalogue[rec].yy(i) << " " << displacement_catalogue[rec].zz(i) << endl;
309  }
310  fout_displ.clear(); fout_displ.close();
311  }
312  }
313  }
314  else ErrorCBL ("algorithm type is not correct!", "Catalogue", "VoidCatalogue.cpp");
315 
316  // --------------------------------------------- //
317  // ---------------- Second Step ---------------- //
318  // --------------------------------------------- //
319 
320  vector<shared_ptr<Catalogue>> displ_cat(n_rec);
321  for (int i=0; i<n_rec; i++) displ_cat[i] = make_shared<Catalogue>(Catalogue(move(displacement_catalogue[i])));
322 
323  double density = tracer_cat->numdensity();
324  double step = step_size*tracer_cat->mps();
325 
326  vector<double> min_X = {random_cat->Min(Var::_X_)-0.5*step, tracer_cat->Min(Var::_X_)-0.5*step};
327  vector<double> min_Y = {random_cat->Min(Var::_Y_)-0.5*step, tracer_cat->Min(Var::_Y_)-0.5*step};
328  vector<double> min_Z = {random_cat->Min(Var::_Z_)-0.5*step, tracer_cat->Min(Var::_Z_)-0.5*step};
329  vector<double> min = {*min_element(min_X.begin(),min_X.end()), *min_element(min_Y.begin(),min_Y.end()), *min_element(min_Z.begin(),min_Z.end())};
330  vector<double> max_X = {random_cat->Max(Var::_X_)+0.5*step, tracer_cat->Max(Var::_X_)+0.5*step};
331  vector<double> max_Y = {random_cat->Max(Var::_Y_)+0.5*step, tracer_cat->Max(Var::_Y_)+0.5*step};
332  vector<double> max_Z = {random_cat->Max(Var::_Z_)+0.5*step, tracer_cat->Max(Var::_Z_)+0.5*step};
333  vector<double> max = {*max_element(max_X.begin(),max_X.end()), *max_element(max_Y.begin(),max_Y.end()), *max_element(max_Z.begin(),max_Z.end())};
334 
335  vector<double> cells(3);
336  for(unsigned int&& i=0; i<3; i++) cells[i] = (max[i]-min[i])/step;
337  unsigned int nCells = *max_element(cells.begin(), cells.end());
338  auto min_copy = min;
339  for(unsigned int&& i=0; i<3; i++) {
340  min[i] = min_copy[i] - (step*nCells-(max[i]-min_copy[i]))/2;
341  max[i] = max[i] + (step*nCells-(max[i]-min_copy[i]))/2;
342  }
343 
344  vector<vector<vector<double>>> Divergence(nCells, vector<vector<double>>(nCells, vector<double>(nCells, 0.)));
345 
346  cout << endl;
347  coutCBL << "* * * Estimation of the Divergence field * * *" << endl;
348  cout << endl;
349 
350  for (int rec=0; rec<n_rec; rec++) {
351  for (size_t i=0; i<randoms[rec].size(); i++) {
352  vector<vector<double>> coord_part(2);
353  coord_part[0] = {random_cat->xx(randoms[rec][i]), random_cat->yy(randoms[rec][i]), random_cat->zz(randoms[rec][i])};
354  coord_part[1] = {tracer_cat->xx(i), tracer_cat->yy(i), tracer_cat->zz(i)};
355  vector<double> displ = {displ_cat[rec]->xx(i), displ_cat[rec]->yy(i), displ_cat[rec]->zz(i)};
356  vector<vector<unsigned int>> inds(2, vector<unsigned int>(3));
357  for (unsigned int&& j=0; j<3; j++) {
358  inds[0][j] = (coord_part[0][j]-min[j])/step;
359  inds[1][j] = (coord_part[1][j]-min[j])/step;
360  }
361  if(inds[0]!=inds[1]) {
362  for (unsigned int part=0; part<2; part++) {
363  for (unsigned int j=0; j<2; j++) {
364  for (unsigned int&& k=0; k<3; k++) {
365  double t = (coord_part[0][k]-((inds[part][k]+j)*step+min[k]))/displ[k];
366  if (t>=0. && t<=1.) {
367  vector<double> inters_coord = {coord_part[0][0]-(coord_part[0][0]-coord_part[1][0])*t,
368  coord_part[0][1]-(coord_part[0][1]-coord_part[1][1])*t,
369  coord_part[0][2]-(coord_part[0][2]-coord_part[1][2])*t};
370  if ((k==0) ? (inters_coord[1]>(inds[part][1]*step+min[1]) && inters_coord[1] <= ((inds[part][1]+1)*step+min[1]) &&
371  inters_coord[2]>(inds[part][2]*step+min[2]) && inters_coord[2] <= ((inds[part][2]+1)*step+min[2])) :
372  ((k==1) ? (inters_coord[0]>(inds[part][0]*step+min[0]) && inters_coord[0] <= ((inds[part][0]+1)*step+min[0]) &&
373  inters_coord[2]>(inds[part][2]*step+min[2]) && inters_coord[2] <= ((inds[part][2]+1)*step+min[2])) :
374  (inters_coord[0]>(inds[part][0]*step+min[0]) && inters_coord[0] <= ((inds[part][0]+1)*step+min[0]) &&
375  inters_coord[1]>(inds[part][1]*step+min[1]) && inters_coord[1] <= ((inds[part][1]+1)*step+min[1])))) {
376  if(part==0) Divergence[inds[part][0]][inds[part][1]][inds[part][2]] -= abs(displ[k])/step;
377  if(part==1) Divergence[inds[part][0]][inds[part][1]][inds[part][2]] += abs(displ[k])/step;
378  }
379  }
380  else continue;
381  }
382  }
383  }
384  }
385  }
386  }
387 
388  coutCBL << "* * * Smoothing * * *" << endl << endl;
389 
390  vector<vector<vector<double>>> Divergence_copy = Divergence;
391 
392  for (unsigned int i=0; i<nCells; i++) {
393  for (unsigned int j=0; j<nCells; j++) {
394  for (unsigned int k=0; k<nCells; k++) {
395  Divergence_copy[i][k][j] = Divergence_copy[i][k][j]/n_rec;
396  }
397  }
398  }
399 
400  for (unsigned int i=0; i<nCells; i++) {
401  for (unsigned int j=0; j<nCells; j++) {
402  for (unsigned int k=0; k<nCells; k++) {
403  double sum=0.;
404  double sum_div=0.;
405  for (int ii=std::max(0,(int)i-1); ii<std::min((int)nCells,(int)i+2); ii++) {
406  for (int jj=std::max(0,(int)j-1); jj<std::min((int)nCells,(int)j+2); jj++) {
407  for (int kk=std::max(0,(int)k-1); kk<std::min((int)nCells,(int)k+2); kk++) {
408  double ex = exp(-pow(Euclidean_distance(min[0]+step*(i+0.5), min[0]+step*(ii+0.5), min[1]+step*(j+0.5), min[1]+step*(jj+0.5), min[2]+step*(k+0.5), min[2]+step*(kk+0.5)),2)/(2*step*step));
409  sum += ex;
410  sum_div += Divergence_copy[ii][jj][kk]*ex;
411  }
412  }
413  }
414  Divergence[i][j][k] = sum_div/sum;
415  }
416  }
417  }
418 
419  coutCBL << "* * * Identification of voids * * *" << endl << endl;
420 
421  if (print[1]) {
422  string name;
423  name = dir_output + "Divergence";
424  name.append("_");
425  name.append(output);
426  ofstream fout_divergence(name.c_str());
427 
428  fout_divergence.precision(7);
429 
430  coutCBL << "Writing the divergence field in " << name << " ..." << endl;
431  cout << endl;
432 
433  for (unsigned int&& i=0; i<nCells; i++)
434  for (unsigned int&& j=0; j<nCells; j++)
435  for (unsigned int&& k=0; k<nCells; k++) {
436  fout_divergence << min[0]+step*(i+0.5) << " " << min[1]+step*(j+0.5) << " " << min[2]+step*(k+0.5) << " " << Divergence[i][j][k] << endl;
437  }
438  fout_divergence.clear(); fout_divergence.close();
439  }
440 
441  vector<vector<double>> coord_locmin(6);
442  unsigned int cont = 0;
443 
444  for (unsigned int&& i=0; i<nCells; i++) {
445  for (unsigned int&& j=0; j<nCells; j++) {
446  for (unsigned int&& k=0; k<nCells; k++) {
447  if (Divergence[i][j][k] < 0. && i>0 && j>0 && k>0 && i<nCells-1 && j<nCells-1 && k<nCells-1) {
448  bool control = false;
449  double sum = 0.;
450  vector<double> num(3, 0.), den(3, 0.);
451  vector<int> delta_cells(3, 0);
452  cont++;
453  for (unsigned int&& ii=i-1; ii<i+2; ii++) {
454  delta_cells[0] = i-ii;
455  for (unsigned int&& jj=j-1; jj<j+2; jj++) {
456  delta_cells[1] = j-jj;
457  for (unsigned int&& kk=k-1; kk<k+2; kk++) {
458  delta_cells[2] = k-kk;
459  if (!(ii==i && jj==j && kk==k)) {
460  sum += Divergence[ii][jj][kk];
461  for (unsigned int&& t=0; t<3; t++) {
462  num[t] += Divergence[ii][jj][kk]*(step/Euclidean_distance((double)i, (double)ii, (double)j, (double)jj, (double)k, (double)kk))*delta_cells[t];
463  den[t] += abs(Divergence[ii][jj][kk]);
464  }
465  if (Divergence[ii][jj][kk]<=Divergence[i][j][k] || Divergence[ii][jj][kk] > 0.) {
466  control = true;
467  break;
468  }
469  }
470  }
471  if (control==true) break;
472  }
473  if (control==true) break;
474  }
475  if (control==false) {
476  coord_locmin[0].push_back(min[0]+step*(i+0.5)+num[0]/den[0]);
477  coord_locmin[1].push_back(min[1]+step*(j+0.5)+num[1]/den[1]);
478  coord_locmin[2].push_back(min[2]+step*(k+0.5)+num[2]/den[2]);
479  coord_locmin[3].push_back(i);
480  coord_locmin[4].push_back(j);
481  coord_locmin[5].push_back(k);
482  }
483  }
484  }
485  }
486  }
487 
488  coutCBL << "Number of cells with negative divergence: " << cont << endl;
489 
490  if (coord_locmin[0].size()==0)
491  ErrorCBL("No local minima found!", "Catalogue", "VoidCatalogue.cpp");
492 
493  coutCBL << "Number of voids: " << coord_locmin[0].size() << endl << endl;
494 
495  // -------------------------------------------- //
496  // ---------------- Third Step ---------------- //
497  // -------------------------------------------- //
498 
499  coutCBL << "* * * Voids rescaling * * *" << endl << endl;
500 
501  vector<double> radius_void(coord_locmin[0].size());
502  double mps = tracer_cat->mps();
503 
504  for (size_t i=0; i<coord_locmin[0].size(); i++) {
505  vector<double> pos = {coord_locmin[0][i], coord_locmin[1][i], coord_locmin[2][i]};
506  vector<unsigned int> close = ChainMesh_tracers.Close_objects(pos, 5*mps);
507  vector<double> distances(close.size());
508  for (size_t j=0; j<close.size(); j++) distances[j] = Euclidean_distance(tracer_cat->xx(close[j]), pos[0], tracer_cat->yy(close[j]), pos[1],
509  tracer_cat->zz(close[j]), pos[2]);
510  std::sort(close.begin(), close.end());
511  vector<vector<double>> data(2, vector<double>(distances.size())); //distances, density contrast
512 
513  for (size_t i=0; i<distances.size()-1; i++) {
514  data[0][i]=(distances[i]+distances[i+1])/2;
515  data[1][i]=(i+1)/(volume_sphere(data[0][i])*density);
516  }
517 
518  data[0][distances.size()-1] = 2*data[0][distances.size()-2] - data[0][distances.size()-3];
519  data[1][distances.size()-1] = distances.size()/(volume_sphere(data[0][distances.size()-1])*density);
520 
521  size_t N = 0;
522  while (!(data[1][N] < threshold && data[1][N+1] > threshold) && N < distances.size()-1) N++;
523  if (N == distances.size()-1) {
524  int index = distances.size()-1;
525  double max = data[1][distances.size()-1];
526  for (size_t k=0; k<data[0].size()-1; k++) {
527  if (data[1][k]<data[1][k+1] && data[1][k+1]<0.5 && data[1][k+1]>max) {
528  index = k;
529  max = data[1][k+1];
530  }
531  }
532  radius_void[i] = interpolated(0.5, {0., data[1][index]}, {0., data[0][index]}, "Linear");
533  }
534  else
535  radius_void[i] = interpolated(0.5, {data[1][N], data[1][N+1]}, {data[0][N], data[0][N+1]}, "Linear");
536  }
537 
538  string name;
539  name = dir_output+"Voids_";
540  name.append(output);
541  ofstream fout_voids(name.c_str());
542  for(int&& i=0; i<(int)radius_void.size(); ++i) fout_voids << coord_locmin[0][i] << " " << coord_locmin[1][i] << " " << coord_locmin[2][i] << " " << radius_void[i] << endl;
543  fout_voids.clear(); fout_voids.close();
544 
545 
546  // ---------------------------------------------------------- //
547  // ---------------- Computing execution time ---------------- //
548  // ---------------------------------------------------------- //
549 
550  double seconds = omp_get_wtime() - start_time;
551 
552  cout << endl;
553  coutCBL << "Time spent to compute: " << seconds << " seconds " << endl ;
554  coutCBL << "Time spent to compute: " << seconds/60. << " minutes " << endl ;
555  coutCBL << "Time spent to compute: " << seconds/3600. << " hours " << endl ;
556  cout << endl;
557 }
558 
560 
561 cbl::catalogue::Catalogue::Catalogue (const VoidAlgorithm algorithm, Catalogue tracer_catalogue, Catalogue random_catalogue, const std::string dir_output, const std::string output, const double cellsize, const cbl::cosmology::Cosmology cosm, const std::vector<double> RA_range, const std::vector<double> DEC_range, const int n_rec, const double step_size, const double threshold)
562 {
563  cout << endl; coutCBL << par::col_green << "BitVF for Light cones" << par::col_default << endl << endl;
564  ErrorCBL ("Still not avaiable!", "Catalogue", "VoidCatalogue.cpp");
565  cout << dir_output << output << cellsize << n_rec << step_size << threshold << cosm.Omega_baryon() << endl;
566  Print(RA_range);
567  Print(DEC_range);
568  if (algorithm == VoidAlgorithm::_LaZeVo_) cout << "aaa" << endl;
569  if (tracer_catalogue.nObjects() == 0) cout << "aaa" << endl;
570  if (random_catalogue.nObjects() == 0) cout << "aaa" << endl;
571 
572  /*
573  // -------------------------------------------- //
574  // ---------------- First Step ---------------- //
575  // -------------------------------------------- //
576 
577  const double start_time = omp_get_wtime();
578  shared_ptr<Catalogue> tracer_cat, random_cat;
579  vector<Catalogue> displacement_catalogue(n_rec);
580  vector<vector<unsigned int>> randoms;
581  vector<double> cells_edges;
582  std::setprecision(5);
583  vector<vector<double>> prop = tracer_catalogue.compute_catalogueProperties_lightCone(cosm, RA_range, DEC_range, 25);
584 
585 
586  auto poly = *[](const vector<double> x, vector<double> &parameter) -> const vector<double> {
587 
588  vector<double> model(x.size());
589  for (size_t i=0; i<x.size(); i++)
590  model[i] = parameter[0]*x[i]*x[i]*x[i] + parameter[1]*x[i]*x[i] + parameter[2]*x[i] + parameter[3];
591 
592  return model;
593  };
594 
595  const function<vector<double>(vector<double>,vector<double>&, double)> mps_func = [&poly](const vector<double> x, std::vector<double> &parameter, double stepsize) {
596  vector<double> model(x.size());
597  for (size_t i=0; i<x.size(); i++)
598  model[i] = stepsize*poly({x[i]}, parameter)[0];
599 
600  return model;
601  };
602 
603  if (algorithm==VoidAlgorithm::_LaZeVo_) {
604 
605  Catalogue random_catalogue;
606 
607  long int time = (long int)chrono::steady_clock::now().time_since_epoch().count();
608  //time = (time%1000000);
609  time = 2343;
610  int rndd = time;
611 
612  double mps_trcat = tracer_catalogue.mps();
613 
614  vector<vector<double>> data_trcat(3);
615  data_trcat[0] = prop[0];
616  for (size_t i=0; i<data_trcat[0].size(); i++) data_trcat[0][i] = cosm.D_C(data_trcat[0][i]);
617  data_trcat[1]=prop[4];
618  data_trcat[2]=prop[6];
619 
620  string name = "../output/Mps_";
621  name.append(output);
622  ofstream fout_displ(name.c_str());
623  cout << endl;
624  for (size_t i=0; i<data_trcat[0].size(); i++)
625  {
626  fout_displ.precision(5);
627  fout_displ << data_trcat[0][i] << " " << data_trcat[1][i] << " " << data_trcat[2][i] << endl;
628  }
629  fout_displ.clear(); fout_displ.close();
630 
631  const cbl::data::Data1D data(data_trcat[0], data_trcat[1], data_trcat[2]);
632  shared_ptr<cbl::data::Data> ptr_data = make_shared<cbl::data::Data1D>(data);
633 
634  const int nparameters = 4;
635  const vector<string> parNames = {"A", "B", "C", "D"};
636  double valA = 0., valB = 0., valC = 0., valD = data_trcat[1][0];
637 
638  vector<cbl::statistics::ParameterType> parType(nparameters, cbl::statistics::ParameterType::_Base_);
639  auto ptr_modelInput = make_shared<cbl::cosmology::Cosmology>(cosm);
640 
641  const statistics::Model1D model(poly, nparameters, parType, parNames, ptr_modelInput);
642  auto ptr_model = make_shared<cbl::statistics::Model1D>(model);
643 
644  cbl::statistics::Likelihood likelihood(ptr_data, ptr_model, cbl::statistics::LikelihoodType::_Gaussian_Error_);
645 
646  const vector<double> starting_p = {valA, valB, valC, valD};
647  double minA = -1., maxA = 1.;
648  double minB = -1., maxB = 1.;
649  double minC = -1., maxC = 1.;
650  double minD = 0., maxD = data_trcat[1][0];
651  const vector<vector<double>> limits = { {minA, maxA}, {minB, maxB}, {minC, maxC}, {minD, maxD} };
652 
653  likelihood.maximize(starting_p, limits);
654 
655  const double start=cosm.D_C(cbl::Min(tracer_catalogue.var(Var::_Redshift_)))-data_trcat[1][0]/2;
656  const double end=cosm.D_C(cbl::Max(tracer_catalogue.var(Var::_Redshift_)))+data_trcat[1][data_trcat[1].size()-1]/2;
657  double front=start;
658  double toll=data_trcat[1][0]*0.001;
659  cells_edges.push_back(start);
660 
661  while (front < end) {
662  vector<double> par = likelihood.parameters()->bestfit_value();
663  par[2] = par[2]-1;
664  par[3] = par[3]+front;
665  while (mps_func({front-(data_trcat[0][1]-data_trcat[0][0])/10.}, par, step_size)[0]*
666  mps_func({front+(data_trcat[0][1]-data_trcat[0][0])/10.}, par, step_size)[0] > 0) front+=(data_trcat[0][1]-data_trcat[0][0])/10.;
667  double left = front-(data_trcat[0][1]-data_trcat[0][0])/10., right = front+(data_trcat[0][1]-data_trcat[0][0])/10.;
668  while ( !(mps_func({front}, par, step_size)[0] > -toll && mps_func({front}, par, step_size)[0] < toll) ) {
669  if (left<0.) {
670  if (mps_func({front}, par, step_size)[0]<0.) left = front;
671  else right = front;
672  front = (right+left)/2;
673  }
674  else {
675  if (mps_func({front}, par, step_size)[0]<0.) right = front;
676  else left = front;
677  front = (right+left)/2;
678  }
679  }
680  cells_edges.push_back(front);
681  }
682  if (random_catalogue.nObjects()==0)
683  {
684  random_catalogue = Catalogue(RandomType::_createRandom_homogeneous_LC_, tracer_catalogue, 1., cosm, RA_range, DEC_range, 100, rndd);
685  coutCBL << "... random catalogue created!" << endl;
686  }
687  else if (random_catalogue.nObjects() == tracer_catalogue.nObjects())
688  {
689  coutCBL << "Reading random catalogue of " << random_catalogue.nObjects() << " objects provides by user..." << endl;
690  random_catalogue = random_catalogue;
691  coutCBL << "...Done!" << endl;
692  }
693  else {
694  coutCBL << "Reading random catalogue of " << random_catalogue.nObjects() << " objects provides by user..." << endl;
695  cout << endl;
696  int diff = random_catalogue.nObjects() - tracer_catalogue.nObjects();
697  if (random_catalogue.nObjects() > tracer_catalogue.nObjects()) coutCBL << "WARNING! Random catalogue has " <<
698  diff << " more objects than the tracers catalogue. The number will be equalized by the function catalogue::equalize_random_LC." << endl;
699  if (random_catalogue.nObjects() < tracer_catalogue.nObjects()) coutCBL << "WARNING! Random catalogue has " <<
700  -diff << " less objects than the tracers catalogue. The number will be equalized by the function catalogue::equalize_random_LC." << endl;
701  random_catalogue.equalize_random_lightCone(tracer_catalogue, cosm, RA_range, DEC_range, rndd);
702  coutCBL << "...Done!" << endl;
703  }
704 
705  // shared pointers to the tracer and the random catalogues
706  random_cat = make_shared<cbl::catalogue::Catalogue> (cbl::catalogue::Catalogue(move(random_catalogue)));
707 tracer_cat = make_shared<cbl::catalogue::Catalogue> (cbl::catalogue::Catalogue(move(tracer_catalogue)));
708 
709 // maximum for the number of objects in the tracer catalogue (unsigned int)
710 unsigned int num_objects = ((tracer_cat->nObjects())%2==0) ? tracer_cat->nObjects() : (tracer_cat->nObjects())-1;
711 
712 if ((tracer_cat->nObjects())%2 != 0) {
713  random_cat -> remove_object(num_objects);
714  tracer_cat -> remove_object(num_objects);
715  }
716 
717 // vector for the indices of tracers
718 std::vector<unsigned int> index_tracer_cat(num_objects);
719 std::iota (std::begin(index_tracer_cat), std::end(index_tracer_cat), 0);
720 cout << "num obj: " << num_objects << endl;
721 cout << endl;
722 coutCBL << "Creating ChainMesh... " << flush;
723 // chainmesh setting
724 vector<Var> variables = {Var::_X_, Var::_Y_, Var::_Z_};
725 CatalogueChainMesh ChainMesh_tracers = CatalogueChainMesh(variables, cellsize, tracer_cat);
726 CatalogueChainMesh ChainMesh_randoms = CatalogueChainMesh(variables, cellsize, random_cat);
727 unsigned int N_near_obj = 3*3*4*par::pi;
728 vector<vector<unsigned int>> near_part = ChainMesh_tracers.N_nearest_objects_cat(N_near_obj); //the selected halo (that will have distance = 0) and the nearest 3))
729 randoms.resize(n_rec, vector<unsigned int>(num_objects)); //tracers fixed
730 for(int rec=0; rec<n_rec; rec++) {
731  std::iota (std::begin(randoms[rec]), std::end(randoms[rec]), 0);
732  auto time2 = (long int)chrono::steady_clock::now().time_since_epoch().count();
733  time2 = (time2%1000000);
734  std::shuffle(randoms[rec].begin(), randoms[rec].end(), default_random_engine(time2));
735  }
736 
737 cout << " Done! " << endl;
738 cout << endl;
739 coutCBL << "Setting starting configuration... " << endl;
740 
741 double dist = 4*mps_trcat;
742 double minX = ChainMesh_tracers.lim()[0][0]+dist/10, maxX = ChainMesh_tracers.lim()[0][1]-dist/10;
743 double minY = ChainMesh_tracers.lim()[1][0]+dist/10, maxY = ChainMesh_tracers.lim()[1][1]-dist/10;
744 double minZ = ChainMesh_tracers.lim()[2][0]+dist/10, maxZ = ChainMesh_tracers.lim()[2][1]-dist/10;
745 
746 // FACCIO STAMPARE DISPL1 E CONFRONTO CON DISPL0
747 // RISOLVERE QUESTIONE COPIA CHAIN MESH E POI TESTARE (SEMBRA VADA MALE E LENTA) E RISOLVERE PROBLEMA IDENTIFICAZIONE VUOTI
748 for (int rec=0; rec<n_rec; rec++) {
749 
750  auto timeS = ((long int)chrono::steady_clock::now().time_since_epoch().count())%1000000;
751  random::UniformRandomNumbers randX(minX, maxX, timeS);
752  random::UniformRandomNumbers randY(minY, maxY, timeS+1);
753  random::UniformRandomNumbers randZ(minZ, maxZ, timeS+2);
754 
755  CatalogueChainMesh ChM_tracer_copy = ChainMesh_tracers;
756  CatalogueChainMesh ChM_random_copy = ChainMesh_randoms;
757  cout << ChM_tracer_copy.part(7).size() << " " << ChainMesh_tracers.part(7).size() << endl;
758  cout << ChM_tracer_copy.cellsize() << endl;
759  unsigned int removed = 0;
760 
761  while (removed < num_objects) {
762  vector<double> pos = {randX(), randY(), randZ()};
763  vector<unsigned int> close=ChM_tracer_copy.Close_objects(pos, dist);
764  unsigned int tr_to_rmv = std::min(100, (int)close.size());
765  if (tr_to_rmv > 0) {
766  vector<unsigned int> close_random = ChM_random_copy.N_nearest_objects(pos, close.size());
767  std::random_device dev;
768  std::mt19937 rng(dev());
769  removed += tr_to_rmv;
770  while (tr_to_rmv > 0) {
771  std::uniform_int_distribution<std::mt19937::result_type> dist(0, tr_to_rmv-1);
772  int rnd1 = dist(rng), rnd2 = dist(rng);
773  randoms[rec][close[rnd1]] = close_random[rnd2];
774  ChM_tracer_copy.deletePart(close[rnd1]);
775  ChM_random_copy.deletePart(close_random[rnd2]);
776  close.erase(close.begin()+rnd1);
777  close_random.erase(close_random.begin()+rnd2);
778  tr_to_rmv--;
779  }
780  }
781  }
782  }
783 
784 
785 coutCBL << "Performing the iterations... " << endl;
786 cout << endl;
787 
788 for (int rec=0; rec<n_rec; rec++) {
789  auto index_tracer_cat_copy = index_tracer_cat;
790  double ratio = 1;
791  unsigned int n_iter = 0;
792 
793  while (ratio > threshold) {
794 
795  index_tracer_cat_copy = index_tracer_cat;
796  auto time2 = (long int)chrono::steady_clock::now().time_since_epoch().count();
797  time2 = (time2%1000000);
798  std::shuffle(index_tracer_cat_copy.begin(), index_tracer_cat_copy.end(), default_random_engine(time2));
799  std::random_device dev;
800  std::mt19937 rng(dev());
801  std::uniform_int_distribution<std::mt19937::result_type> dist(1,N_near_obj-1);
802  vector<bool> index_bool(num_objects, false), used(num_objects, false);
803  //vector<vector<unsigned int>> H(tracer_cat_len_4,vector<unsigned int>(4)), R_def(tracer_cat_len_4,vector<unsigned int>(4));
804 
805 #pragma omp parallel num_threads(omp_get_max_threads())
806  {
807 #pragma omp for schedule(static)
808  for (unsigned int i=0; i<num_objects; i++) {
809  vector<unsigned int> H(4), R_def(4), R(4), R_copy(4);
810  unsigned int rand1 = dist(rng);
811  unsigned int rand2 = rand1;
812  unsigned int rand3 = rand1;
813  while (rand1 == rand2 || rand1==rand3 || rand2 == rand3) {
814  rand2 = dist(rng);
815  rand3 = dist(rng);
816  }
817  //vector<unsigned int> R(4);
818  if (used[near_part[index_tracer_cat_copy[i]][0]] == false && used[near_part[index_tracer_cat_copy[i]][rand1]] == false &&
819  used[near_part[index_tracer_cat_copy[i]][rand2]] == false && used[near_part[index_tracer_cat_copy[i]][rand3]] == false) {
820  // selezione particelle coinvolte
821  used[near_part[index_tracer_cat_copy[i]][0]] = true;
822  used[near_part[index_tracer_cat_copy[i]][rand1]] = true;
823  used[near_part[index_tracer_cat_copy[i]][rand2]] = true;
824  used[near_part[index_tracer_cat_copy[i]][rand3]] = true;
825  H[0]=near_part[index_tracer_cat_copy[i]][0];
826  H[1]=near_part[index_tracer_cat_copy[i]][rand1];
827  H[2]=near_part[index_tracer_cat_copy[i]][rand2];
828  H[3]=near_part[index_tracer_cat_copy[i]][rand3];
829  R[0]=randoms[rec][H[0]];
830  R[1]=randoms[rec][H[1]];
831  R[2]=randoms[rec][H[2]];
832  R[3]=randoms[rec][H[3]];
833 
834  R_copy = R;
835  // calcolo configurazione minima
836  double dist_min = 1000000.;
837  double dist;
838  do {
839  dist = 0.;
840  for (size_t j=0; j<R.size(); j++)
841  dist += cbl::Euclidean_distance(tracer_cat->xx(H[j]), random_cat->xx(R[j]), tracer_cat->yy(H[j]), random_cat->yy(R[j]), tracer_cat->zz(H[j]), random_cat->zz(R[j]));
842 
843  if (dist < dist_min) {
844  dist_min = dist;
845  R_def = R;
846  }
847  } while(std::next_permutation(R.begin(), R.end()));
848 
849  //aggiornamento vettori tracers e randoms
850  if(R_def != R_copy) index_bool[i] = true;
851 
852  for (size_t j=0; j<4; j++) randoms[rec][H[j]] = R_def[j];
853 
854  used[H[0]] = false;
855  used[H[1]] = false;
856  used[H[2]] = false;
857  used[H[3]] = false;
858  }
859  else continue;
860  }
861  }
862 
863  unsigned int changed_couples = count(index_bool.begin(), index_bool.end(), true);
864 
865  ratio = (double)changed_couples/num_objects;
866  if (n_iter%100==0) cout << setprecision(5) << "ratio: " << ratio << " " << n_iter << endl;
867  n_iter++;
868  }
869 
870  double distanza = 0.;
871  for (unsigned int i=0; i<num_objects; i++) {
872  distanza += cbl::Euclidean_distance(tracer_cat->xx(i), random_cat->xx(randoms[rec][i]),
873  tracer_cat->yy(i), random_cat->yy(randoms[rec][i]), tracer_cat->zz(i), random_cat->zz(randoms[rec][i]));
874  }
875  cout << "distanza finale: " << distanza/num_objects << endl;
876 
877  vector<shared_ptr<cbl::catalogue::Object>> displ_catalogue_object;
878 
879  for (size_t ii = 0; ii<randoms[rec].size(); ii++)
880  {
881  cbl::comovingCoordinates displacement_values = {random_cat->xx(randoms[rec][ii])-tracer_cat->xx(ii), random_cat->yy(randoms[rec][ii])-tracer_cat->yy(ii), random_cat->zz(randoms[rec][ii])-tracer_cat->zz(ii)};
882  auto tracer_displ = make_shared<cbl::catalogue::Halo>(displacement_values);
883  displ_catalogue_object.push_back(tracer_displ);
884  }
885 
887  displacement_catalogue[rec].add_objects(displ_catalogue_object);
888 
889  if (rec == 0) {
890  stringstream ss;
891  ss << rec;
892  string str = ss.str();
893  // string output = ss.str();
894  string name = "../output/Displacement";
895  name.append("_");
896  name.append(str);
897  name.append("_");
898  // output = "coord_LCDM.dat";
899  name.append(output);
900  ofstream fout_displ(name.c_str());
901  cout << endl;
902 
903  coutCBL << "Writing the displacement ... " << endl;
904  for (unsigned int i=0; i<num_objects; i++)
905  {
906  fout_displ.precision(7);
907  fout_displ << tracer_cat->xx(i) << " " << tracer_cat->yy(i) << " " << tracer_cat->zz(i) << " "
908  << random_cat->xx(randoms[rec][i]) << " " << random_cat->yy(randoms[rec][i]) << " " << random_cat->zz(randoms[rec][i]) << " "
909  << displacement_catalogue[rec].xx(i) << " " << displacement_catalogue[rec].yy(i) << " " << displacement_catalogue[rec].zz(i) << endl;
910  } // modified
911 
912  fout_displ.clear(); fout_displ.close();
913  }
914  }
915 
916 
917 } // End of LaZeVo method
918  else if (algorithm==VoidAlgorithm::_Exact_)
919  {
920 
921  } // End of Exact method
922  else ErrorCBL ("algorithm type is not correct!", "Catalogue", "VoidCatalogue.cpp");
923 
924 // --------------------------------------------- //
925 // ---------------- Second Step ---------------- //
926 // --------------------------------------------- //
927 
928 vector<shared_ptr<Catalogue>> displ_cat(n_rec);
929 for (int i=0; i<n_rec; i++) displ_cat[i] = make_shared<Catalogue>(Catalogue(move(displacement_catalogue[i])));
930 
931 double THETA_min = -DEC_range[1]+par::pi/2, THETA_max = -DEC_range[0]+par::pi/2;
932 double PHI_min = RA_range[0], PHI_max = RA_range[1];
933 
934 cout << endl;
935 coutCBL << "RA range: [" << PHI_min << ", " << PHI_max << "]" << endl;
936 coutCBL << "Dec range: [" << DEC_range[0] << ", " << DEC_range[1] << "]" << endl;
937 cout << endl;
938 
939 double delta_THETA = THETA_max - THETA_min, delta_PHI = RA_range[1]-RA_range[0];
940 vector<double> volume = prop[2];
941 vector<double> density = prop[3];
942 vector<double> step(cells_edges.size()-1), z_cells_edges(cells_edges.size()), z_cells_centres(step.size()), cells_centre(step.size()), THETA_step(step.size()),
943  real_delta_THETA(step.size()), real_THETA_min(step.size());
944 vector<unsigned int> THETA_ncells(step.size());
945 
946 for (size_t i=0; i<step.size(); i++) step[i]=cells_edges[i+1]-cells_edges[i];
947 for (size_t i=0; i<cells_edges.size(); i++) z_cells_edges[i] = cosm.Redshift(cells_edges[i]);
948 for (size_t i=0; i<step.size(); i++) cells_centre[i]=(cells_edges[i+1]+cells_edges[i])/2;
949 for (size_t i=0; i<step.size(); i++) z_cells_centres[i]=cosm.Redshift(cells_centre[i]);
950 for (size_t i=0; i<step.size(); i++) {
951  THETA_step[i] = step[i]/cells_centre[i];
952  if ((((int)(delta_THETA/(step[i]/cells_centre[i]))+1)*THETA_step[i]) > par::pi) {
953  step[i]=par::pi*cells_centre[i]/((int)(delta_THETA/(step[i]/cells_centre[i]))+1);
954  THETA_step[i] = step[i]/cells_centre[i];
955  }
956  }
957 for (size_t i=0; i<step.size(); i++) THETA_ncells[i]=delta_THETA/THETA_step[i]+1;
958 for (size_t i=0; i<step.size(); i++) real_delta_THETA[i]=THETA_ncells[i]*THETA_step[i];
959 for (size_t i=0; i<step.size(); i++) real_THETA_min[i]= (THETA_min - (real_delta_THETA[i]-delta_THETA)/2 < 0.) ? 0. : THETA_min - (real_delta_THETA[i]-delta_THETA)/2;
960 
961 vector<vector<vector<double>>> Divergence(step.size());
962 vector<vector<double>> real_delta_PHI(step.size()), real_PHI_min(step.size()), PHI_step(step.size());
963 
964 for (size_t i=0; i<step.size(); i++) {
965  Divergence[i].resize(THETA_ncells[i]);
966  PHI_step[i].resize(THETA_ncells[i]);
967  real_PHI_min[i].resize(THETA_ncells[i]);
968  real_delta_PHI[i].resize(THETA_ncells[i]);
969  for (unsigned int j=0; j<THETA_ncells[i]; j++) {
970  double angle = real_THETA_min[i] + (j+0.5)*THETA_step[i];
971  double temp_PHI_step = THETA_step[i]/sin(angle);
972  PHI_step[i][j] = temp_PHI_step;
973  if ((((int)(delta_PHI/PHI_step[i][j])+1)*PHI_step[i][j]) > 2*par::pi) PHI_step[i][j] = 2*par::pi/((int)(delta_PHI/PHI_step[i][j])+1);
974  int temp_nCells = delta_PHI/PHI_step[i][j] + 1;
975  real_delta_PHI[i][j] = temp_nCells*PHI_step[i][j];
976  real_PHI_min[i][j] = PHI_min - (real_delta_PHI[i][j]-delta_PHI)/2;
977  Divergence[i][j].resize(temp_nCells);
978  }
979  }
980 
981 vector<double> cell_vol(step.size()), sph_area(step.size()+1), udrl_area(step.size());
982 for (size_t i=0; i<step.size(); i++) cell_vol[i] = (-cos(par::pi/2+THETA_step[i]/2)+cos(par::pi/2-THETA_step[i]/2))*THETA_step[i]/(4*par::pi)*
983  (volume_sphere(cells_edges[i+1])-volume_sphere(cells_edges[i]));
984 for (size_t i=0; i<=step.size()+1; i++) sph_area[i] = (-cos(par::pi/2+THETA_step[i]/2)+cos(par::pi/2-THETA_step[i]/2))*THETA_step[i]*cells_edges[i]*cells_edges[i];
985 for (size_t i=0; i<step.size(); i++) udrl_area[i] = THETA_step[i]*(cells_edges[i+1]*cells_edges[i+1]-cells_edges[i]*cells_edges[i])/2;
986 
987 unsigned int nCells = 0;
988 for (size_t i=0; i<step.size(); i++)
989  for (size_t j=0; j<Divergence[i].size(); j++) nCells+=Divergence[i][j].size();
990 
991 cout << endl;
992 coutCBL << "Compute the Divergence field on " << nCells << " cells ..." << endl;
993 cout << endl;
994 
995 for(int rec=0; rec<n_rec; rec++) {
996  for (size_t i=0; i<randoms[rec].size(); i++) {
997  vector<vector<double>> coord_part(2);
998  coord_part[0] = {random_cat->redshift(randoms[rec][i]), -(random_cat->dec(randoms[rec][i]))+par::pi/2,
999  (DEC_range[0] < 0. && random_cat->ra(randoms[rec][i]) > par::pi) ? random_cat->ra(randoms[rec][i])-2*par::pi : random_cat->ra(randoms[rec][i]),
1000  random_cat->xx(randoms[rec][i]), random_cat->yy(randoms[rec][i]), random_cat->zz(randoms[rec][i])};
1001  coord_part[1] = {tracer_cat->redshift(i), -(tracer_cat->dec(i))+par::pi/2,
1002  (DEC_range[0] < 0. && tracer_cat->ra(i) > par::pi) ? tracer_cat->ra(i)-2*par::pi : tracer_cat->ra(i),
1003  tracer_cat->xx(i), tracer_cat->yy(i), tracer_cat->zz(i)};
1004  vector<double> displ = {displ_cat[rec]->xx(randoms[rec][i]), displ_cat[rec]->yy(randoms[rec][i]), displ_cat[rec]->zz(randoms[rec][i])};
1005  vector<vector<unsigned int>> inds(2, vector<unsigned int>(3));
1006  // selezione cella
1007  for (int p=0; p<2; p++) {
1008  int j = 0;
1009  while (coord_part[p][0] > z_cells_edges[j+1]) j++;
1010  inds[p][0] = j;
1011  inds[p][1] = (coord_part[p][1]-real_THETA_min[j])/THETA_step[j];
1012  inds[p][2] = (coord_part[p][2]-real_PHI_min[j][inds[p][1]])/PHI_step[j][inds[p][1]];
1013  }
1014 
1015  if(inds[0]!=inds[1]) {
1016  for (int p=0; p<2; p++) {
1017  // possibili intersezioni
1018  // intersezione facce sferiche
1019  for (int sph=0; sph<2; sph++) {
1020  vector<double> coord_intersect;
1021  double a, b, c;
1022  a = (coord_part[1][3]-coord_part[0][3])*(coord_part[1][3]-coord_part[0][3]) +
1023  (coord_part[1][4]-coord_part[0][4])*(coord_part[1][4]-coord_part[0][4]) +
1024  (coord_part[1][5]-coord_part[0][5])*(coord_part[1][5]-coord_part[0][5]);
1025  b = 2*coord_part[0][3]*(coord_part[1][3]-coord_part[0][3]) + 2*coord_part[0][4]*(coord_part[1][4]-coord_part[0][4]) +
1026  2*coord_part[0][5]*(coord_part[1][5]-coord_part[0][5]);
1027  c = coord_part[0][3]*coord_part[0][3] + coord_part[0][4]*coord_part[0][4] + coord_part[0][5]*coord_part[0][5] -
1028  cells_edges[inds[p][0]+sph]*cells_edges[inds[p][0]+sph];
1029  double del = b*b - 4*a*c;
1030  if (del>0) {
1031  double u1, u2;
1032  u1 = (-b+sqrt(del))/(2*a);
1033  u2 = (-b-sqrt(del))/(2*a);
1034  if( u1<=1. && u1>=0. ) coord_intersect = {coord_part[0][3]+(coord_part[1][3]-coord_part[0][3])*u1,
1035  coord_part[0][4]+(coord_part[1][4]-coord_part[0][4])*u1,
1036  coord_part[0][5]+(coord_part[1][5]-coord_part[0][5])*u1};
1037  else if ( u2<=1. && u2>=0. ) coord_intersect = {coord_part[0][3]+(coord_part[1][3]-coord_part[0][3])*u2,
1038  coord_part[0][4]+(coord_part[1][4]-coord_part[0][4])*u2,
1039  coord_part[0][5]+(coord_part[1][5]-coord_part[0][5])*u2};
1040  else continue;
1041  //calcolo ora il versore alla superficie
1042  double phi, theta;
1043  phi = (RA_range[0] > 0. && atan(coord_intersect[1]/coord_intersect[0]) < 0.) ?
1044  atan(coord_intersect[1]/coord_intersect[0]) + 2*par::pi : atan(coord_intersect[1]/coord_intersect[0]);
1045  theta = acos(coord_intersect[2]/cells_edges[inds[p][0]+sph]);
1046  if (theta>inds[p][1]*THETA_step[inds[p][0]]+real_THETA_min[inds[p][0]] && theta <= (inds[p][1]+1)*THETA_step[inds[p][0]]+real_THETA_min[inds[p][0]] &&
1047  phi>inds[p][2]*PHI_step[inds[p][0]][inds[p][1]]+real_PHI_min[inds[p][0]][inds[p][1]] && phi<=(inds[p][2]+1)*PHI_step[inds[p][0]][inds[p][1]]+real_PHI_min[inds[p][0]][inds[p][1]]) {
1048  double normaliz = sqrt(coord_intersect[0]*coord_intersect[0]+coord_intersect[1]*coord_intersect[1]+coord_intersect[2]*coord_intersect[2]);
1049  vector<double> norm = {coord_intersect[0]/normaliz, coord_intersect[1]/normaliz, coord_intersect[2]/normaliz};
1050  // calcolo il prodotto scalare tra normale e displ
1051  if (sph == 0) {
1052  norm[0] = -norm[0];
1053  norm[1] = -norm[1];
1054  norm[2] = -norm[2];
1055  }
1056  double scal = norm[0]*displ[0] + norm[1]*displ[1] + norm[2]*displ[2];
1057  //computo DIVERGENZA
1058  Divergence[inds[p][0]][inds[p][1]][inds[p][2]] += scal*sph_area[inds[p][0]+sph]/cell_vol[inds[p][0]];
1059  }
1060  }
1061  }
1062 
1063  // calcolo divergenza sulle 2 facce della cella sopra e sotto (intersezione cono/retta)
1064  for (int k=0; k<2; k++) {
1065  double a, b, c, theta = THETA_step[inds[p][0]]*(inds[p][1]+k)+real_THETA_min[inds[p][0]];
1066  double l = coord_part[1][3]-coord_part[0][3], m = coord_part[1][4]-coord_part[0][4], n = coord_part[1][5]-coord_part[0][5];
1067 
1068  a = (l*l + m*m - n*n*tan(theta)*tan(theta));
1069  b = 2*(l*coord_part[0][3] + m*coord_part[0][4] - n*coord_part[0][5]*tan(theta)*tan(theta));
1070  c = coord_part[0][3]*coord_part[0][3] + coord_part[0][4]*coord_part[0][4] - coord_part[0][5]*coord_part[0][5]*tan(theta)*tan(theta);
1071 
1072  double t1, t2, delta;
1073  delta = b*b - 4*a*c;
1074 
1075  if (delta > 0) {
1076  t1 = (-b+sqrt(delta))/(2*a);
1077  t2 = (-b-sqrt(delta))/(2*a);
1078  }
1079  else continue;
1080  vector<double> coord_intersect;
1081  if( t1<=1. && t1>=0. ) coord_intersect = {coord_part[0][3]+(coord_part[1][3]-coord_part[0][3])*t1,
1082  coord_part[0][4]+(coord_part[1][4]-coord_part[0][4])*t1,
1083  coord_part[0][5]+(coord_part[1][5]-coord_part[0][5])*t1};
1084  else if ( t2<=1. && t2>=0. ) coord_intersect = {coord_part[0][3]+(coord_part[1][3]-coord_part[0][3])*t2,
1085  coord_part[0][4]+(coord_part[1][4]-coord_part[0][4])*t2,
1086  coord_part[0][5]+(coord_part[1][5]-coord_part[0][5])*t2};
1087  else continue;
1088 
1089  double r = sqrt(coord_intersect[0]*coord_intersect[0]+coord_intersect[1]*coord_intersect[1]+coord_intersect[2]*coord_intersect[2]);
1090  double phi = (RA_range[0] > 0. && atan(coord_intersect[1]/coord_intersect[0]) < 0.) ? atan(coord_intersect[1]/coord_intersect[0]) + 2*par::pi : atan(coord_intersect[1]/coord_intersect[0]);
1091  if (r>cells_edges[inds[p][0]] && r<cells_edges[inds[p][0]+1] &&
1092  phi>inds[p][2]*PHI_step[inds[p][0]][inds[p][1]]+real_PHI_min[inds[p][0]][inds[p][1]] && phi<=(inds[p][2]+1)*PHI_step[inds[p][0]][inds[p][1]]+real_PHI_min[inds[p][0]][inds[p][1]]) {
1093  vector<double> norm;
1094  if (k==0) norm = {2*coord_intersect[0]*cos(theta)*cos(theta), 2*coord_intersect[1]*cos(theta)*cos(theta), -2*coord_intersect[0]*sin(theta)*sin(theta)};
1095  if (k==1) norm = {-2*coord_intersect[0]*cos(theta)*cos(theta), -2*coord_intersect[1]*cos(theta)*cos(theta), 2*coord_intersect[0]*sin(theta)*sin(theta)};
1096  double normaliz = sqrt(norm[0]*norm[0]+norm[1]*norm[1]+norm[2]*norm[2]);
1097  norm = {norm[0]/normaliz, norm[1]/normaliz, norm[2]/normaliz};
1098  double scal = norm[0]*displ[0] + norm[1]*displ[1] + norm[2]*displ[2];
1099  Divergence[inds[p][0]][inds[p][1]][inds[p][2]] += scal*udrl_area[inds[p][0]+k]/cell_vol[inds[p][0]];
1100  }
1101  }
1102  vector<double> temp_PHI(2);
1103  vector<vector<double>> cent_coord(2);
1104  for(int fac=0; fac<2; fac++) temp_PHI[fac] = PHI_step[inds[p][0]][inds[p][1]]*(inds[p][2]+fac) + real_PHI_min[inds[p][0]][inds[p][1]];
1105 
1106  cent_coord[0] = {par::pi/2, temp_PHI[0]+3*par::pi/2};
1107  cent_coord[1] = {par::pi/2, temp_PHI[1]+par::pi/2};
1108 
1109  for (int n=0; n<2; n++) {
1110  //eq piano per tre punti (di cui uno (0,0,0)) e coeff versore normale
1111  double a,b,c;
1112 
1113  a = sin(cent_coord[n][0])*cos(cent_coord[n][1]);
1114  b = sin(cent_coord[n][0])*sin(cent_coord[n][1]);
1115  c = cos(cent_coord[n][0]);
1116 
1117  // intersezione
1118  double l = coord_part[1][3]-coord_part[0][3], m = coord_part[1][4]-coord_part[0][4], n_ = coord_part[1][5]-coord_part[0][5];
1119 
1120  vector<double> coord_intersect(3);
1121  double t = (a*coord_part[1][3] + b*coord_part[1][4] + c*coord_part[1][5])/(a*l + b*m + c*n_);
1122  if (t>0 && t<1) coord_intersect = {coord_part[1][3]-(coord_part[1][3]-coord_part[0][3])*t,
1123  coord_part[1][4]-(coord_part[1][4]-coord_part[0][4])*t,
1124  coord_part[1][5]-(coord_part[1][5]-coord_part[0][5])*t};
1125  else continue;
1126 
1127  double r, theta;
1128  r = sqrt(coord_intersect[0]*coord_intersect[0] + coord_intersect[1]*coord_intersect[1] + coord_intersect[2]*coord_intersect[2]);
1129  theta = acos(coord_intersect[2]/r);
1130 
1131  //verifica intersezione in cella. ordine: dx, up, sx, down
1132 
1133  if (r >= cells_edges[inds[p][0]] && r < cells_edges[inds[p][0]+1] && theta > THETA_step[inds[p][0]]*inds[p][1]+real_THETA_min[inds[p][0]]
1134  && theta < THETA_step[inds[p][0]]*(inds[p][1]+1)+real_THETA_min[inds[p][0]]) {
1135  // calcolo il prodotto scalare tra normale e displ
1136  double scal = a*displ[0] + b*displ[1] + c*displ[2];
1137  //computo DIVERGENZA
1138  Divergence[inds[p][0]][inds[p][1]][inds[p][2]] += scal*udrl_area[inds[p][0]]/cell_vol[inds[p][0]];
1139  }
1140  }
1141  }
1142  }
1143  }
1144  }
1145 
1146 Catalogue Divergence_cell;
1147 
1148 int n_cell = 0;
1149 for (size_t i=0; i<Divergence.size(); i++) {
1150  for (size_t j=0; j<Divergence[i].size(); j++) {
1151  for (size_t k=0; k<Divergence[i][j].size(); k++) {
1152  double radius = ((cells_edges[i+1]+cells_edges[i])/2);
1153  double theta = (j+0.5)*THETA_step[i] + real_THETA_min[i];
1154  double phi = (k+0.5)*PHI_step[i][j] + real_PHI_min[i][j];
1155  vector<double> cell_values = {radius*sin(theta)*cos(phi), radius*sin(theta)*sin(phi), radius*cos(theta)};
1156  shared_ptr<cbl::catalogue::Object> cell_obj = make_shared<cbl::catalogue::Halo>(cell_values[0], cell_values[1], cell_values[2],
1157  phi, theta, z_cells_centres[i], Divergence[i][j][k]/n_rec);
1158  Divergence_cell.add_object(cell_obj);
1159  n_cell++;
1160  }
1161  }
1162  }
1163 
1164 shared_ptr<Catalogue> div_cat = make_shared<cbl::catalogue::Catalogue> (cbl::catalogue::Catalogue(move(Divergence_cell)));
1165 vector<Var> variables = {Var::_X_, Var::_Y_, Var::_Z_};
1166 CatalogueChainMesh ChainMesh_div = CatalogueChainMesh(variables, cellsize, div_cat);
1167 
1168 coutCBL << "Smoothing..." << endl;
1169 
1170 n_cell = 0;
1171 
1172 for (size_t i=0; i<Divergence.size(); i++) {
1173  for (size_t j=0; j<Divergence[i].size(); j++) {
1174  for (size_t k=0; k<Divergence[i][j].size(); k++) {
1175  if ((i>0 && i<Divergence.size()-1) && (j>0 && j<Divergence[i].size()-1) && (k>0 && k<Divergence[i][j].size()-1)) {
1176  vector<unsigned int> close_cells = ChainMesh_div.Close_objects({div_cat->xx(n_cell), div_cat->yy(n_cell), div_cat->zz(n_cell)}, 2*step[i]);
1177  double sum=0.;
1178  double sum_div=0.;
1179  for(auto&& t : close_cells) {
1180  double ex = exp(-pow(Euclidean_distance(div_cat->xx(n_cell),div_cat->xx(t), div_cat->yy(n_cell), div_cat->yy(t),
1181  div_cat->zz(n_cell), div_cat->zz(t)),2)/(0.5*step[i]*step[i]));
1182  sum += ex;
1183  sum_div += div_cat->var(t, Var::_Weight_)*ex;
1184  }
1185  div_cat->set_var(n_cell, Var::_Weight_, sum_div/sum);
1186  }
1187  n_cell++;
1188  }
1189  }
1190  }
1191 
1192 string name;
1193 name = dir_output + "Divergence";
1194 name.append("_");
1195 name.append(output);
1196 ofstream fout_divergence(name.c_str());
1197 
1198 fout_divergence.precision(7);
1199 vector<vector<double>> coord_locmin(9);
1200 n_cell = 0;
1201 int cont = 0;
1202 
1203 for (size_t i=0; i<Divergence.size(); i++) {
1204  for (size_t j=0; j<Divergence[i].size(); j++) {
1205  for (size_t k=0; k<Divergence[i][j].size(); k++) {
1206  fout_divergence << div_cat->xx(n_cell) << " " << div_cat->yy(n_cell) << " " << div_cat->zz(n_cell) << " " <<
1207  div_cat->ra(n_cell) << " " << div_cat->dec(n_cell) << " " << div_cat->redshift(n_cell) << " " <<
1208  div_cat->weight(n_cell) << endl;
1209 
1210  if (div_cat->weight(n_cell)<0. && (i>0 && i<Divergence.size()-1) && (j>0 && j<Divergence[i].size()-1) && (k>0 && k<Divergence[i][j].size()-1)) {
1211  vector<unsigned int> close_cells = ChainMesh_div.Close_objects({div_cat->xx(n_cell), div_cat->yy(n_cell), div_cat->zz(n_cell)}, 1.5*step[i]);
1212  bool control = false;
1213  double sum = 0., den = 0.;
1214  vector<double> num(3, 0.);
1215  for(auto&& t : close_cells) {
1216  if ((int)t != n_cell) {
1217  if (div_cat->weight(t) < div_cat->weight(n_cell)) {
1218  control = true;
1219  break;
1220  }
1221  sum += div_cat->weight(t);
1222  num[0] += div_cat->weight(t)*(step[i]/Euclidean_distance(div_cat->xx(n_cell),div_cat->xx(t), div_cat->yy(n_cell), div_cat->yy(t),
1223  div_cat->zz(n_cell), div_cat->zz(t))*(div_cat->xx(n_cell)-div_cat->xx(t)));
1224  num[1] += div_cat->weight(t)*(step[i]/Euclidean_distance(div_cat->xx(n_cell),div_cat->xx(t), div_cat->yy(n_cell), div_cat->yy(t),
1225  div_cat->zz(n_cell), div_cat->zz(t))*(div_cat->yy(n_cell)-div_cat->yy(t)));
1226  num[2] += div_cat->weight(t)*(step[i]/Euclidean_distance(div_cat->xx(n_cell),div_cat->xx(t), div_cat->yy(n_cell), div_cat->yy(t),
1227  div_cat->zz(n_cell), div_cat->zz(t))*(div_cat->zz(n_cell)-div_cat->zz(t)));
1228  den += abs(div_cat->weight(t));
1229  }
1230  }
1231  if (control==false) {
1232  coord_locmin[0].push_back(div_cat->xx(n_cell)+num[0]/den);
1233  coord_locmin[1].push_back(div_cat->yy(n_cell)+num[1]/den);
1234  coord_locmin[2].push_back(div_cat->zz(n_cell)+num[2]/den);
1235  double r, theta, phi;
1236  r = sqrt(coord_locmin[0][cont]*coord_locmin[0][cont]+coord_locmin[1][cont]*coord_locmin[1][cont]+coord_locmin[2][cont]*coord_locmin[2][cont]);
1237  phi = atan2(coord_locmin[1][cont], coord_locmin[0][cont]);
1238  theta = acos(coord_locmin[2][cont]/r);
1239  coord_locmin[3].push_back(phi);
1240  coord_locmin[4].push_back(theta);
1241  coord_locmin[5].push_back(cosm.Redshift(r));
1242  coord_locmin[6].push_back(i);
1243  coord_locmin[7].push_back(j);
1244  coord_locmin[8].push_back(k);
1245  cont++;
1246  }
1247  }
1248  n_cell++;
1249  }
1250  }
1251  }
1252 
1253 fout_divergence.clear(); fout_divergence.close();
1254 
1255 unsigned int num_voids = coord_locmin[0].size();
1256 //int number = x_void.size();
1257 if (num_voids==0)
1258  ErrorCBL("No local minima found!", "Catalogue", "VoidCatalogue.cpp");
1259 
1260 coutCBL << "Number of local minima: " << num_voids << endl << endl;
1261 
1262 // -------------------------------------------- //
1263 // ---------------- Third Step ---------------- //
1264 // -------------------------------------------- //
1265 
1266 // Identification of voids //
1267 
1268 
1269 coutCBL << "* * * Identification of voids * * *" << endl << endl;
1270 
1271 vector<double> radius_void(num_voids, 0.);
1272 for (unsigned int&& i=0; i<num_voids; i++) {
1273  double sum = -1., prev_sum;
1274  double incr = 0.;
1275  vector<double> temp_vect = {cosm.D_C(z_cells_edges[z_cells_edges.size()-1]-z_cells_edges[coord_locmin[6][num_voids]]),
1276  cosm.D_C(z_cells_edges[coord_locmin[6][num_voids]]-z_cells_edges[0]),
1277  (coord_locmin[4][num_voids]-THETA_min)*cells_centre[coord_locmin[6][num_voids]], (THETA_max-coord_locmin[4][num_voids])*cells_centre[coord_locmin[6][num_voids]],
1278  (coord_locmin[3][num_voids]-PHI_min)*cells_centre[coord_locmin[6][num_voids]]*sin(coord_locmin[3][num_voids]),
1279  (PHI_max-coord_locmin[3][num_voids])*cells_centre[coord_locmin[6][num_voids]]*sin(coord_locmin[3][num_voids])};
1280  double dist_lim = *min_element(temp_vect.begin(), temp_vect.end());
1281  while (sum<=0.) {
1282  prev_sum = sum;
1283  sum = 0.;
1284  incr += 0.5;
1285  vector<unsigned int> close_cells = ChainMesh_div.Close_objects({coord_locmin[0][num_voids], coord_locmin[1][num_voids],
1286  coord_locmin[2][num_voids]}, incr*step[coord_locmin[6][num_voids]]);
1287 
1288  for (auto&& j : close_cells) sum += div_cat->weight(j);
1289  if (sum>0.) radius_void[i] = interpolated(0., {prev_sum, sum},
1290  {(incr-0.5)*step[coord_locmin[6][num_voids]], incr*step[coord_locmin[6][num_voids]]}, "Linear");
1291  if (sum <0. && incr*step[coord_locmin[6][num_voids]]>=dist_lim) radius_void[i] = dist_lim;
1292  }
1293  }
1294 
1295 name = dir_output+"Voids_";
1296 name.append(output);
1297 ofstream fout_voids(name.c_str());
1298 
1299 for(int&& i=0; i<(int)radius_void.size(); ++i) fout_voids << coord_locmin[0][i] << " " << coord_locmin[1][i] << " " << coord_locmin[2][i] << " " <<
1300  coord_locmin[3][i] << " " << coord_locmin[4][i] << " " << coord_locmin[5][i] << " " << radius_void[i] << endl;
1301 
1302 fout_voids.clear(); fout_voids.close();
1303 
1304 cout << endl;
1305 coutCBL << "Number of voids: " << (int)radius_void.size() << endl << endl;
1306 
1307 
1308 // ---------------------------------------------------------- //
1309 // ---------------- Computing execution time ---------------- //
1310 // ---------------------------------------------------------- //
1311 
1312 //float seconds = float( clock () - begin_time ) / CLOCKS_PER_SEC;
1313 double seconds = omp_get_wtime() - start_time;
1314 
1315 cout << endl;
1316 coutCBL << "Time spent to compute: " << seconds << " seconds " << endl ;
1317 coutCBL << "Time spent to compute: " << seconds/60. << " minutes " << endl ;
1318 coutCBL << "Time spent to compute: " << seconds/3600. << " hours " << endl ;
1319 cout << endl;
1320 cout << dir_output << endl;
1321 */
1322 }
1323 
1324 
1325 
1326 //********************************** VERSION 1.0 **********************************//
1327 
1328 void cbl::catalogue::Catalogue::clean_void_catalogue (const bool initial_radius, const std::vector<double> delta_r, const double threshold, bool rescale, const std::shared_ptr<Catalogue> tracers_catalogue, chainmesh::ChainMesh3D ChM, const double ratio, const bool checkoverlap, const Var ol_criterion)
1329 {
1330  const double start_time = omp_get_wtime();
1331 
1332  // ---------------------------------------------------- //
1333  // ---------------- Cleaning Procedure ---------------- //
1334  // ---------------------------------------------------- //
1335 
1336  cout << endl;
1337  coutCBL << par::col_blue << "> > > > > >>>>>>>>>> CLEANING PROCEDURE STARTED <<<<<<<<< < < < < <" << par::col_default << endl << endl;
1338  coutCBL << "Voids in the initial Catalogue: " << nObjects() << endl << endl;
1339  if (nObjects()==0) ErrorCBL("Empty void catalogue!", "Catalogue", "VoidCatalogue.cpp");
1340 
1341  coutCBL << par::col_green << " --- Removing spurious voids --- " << par::col_default << endl << endl;
1342  coutCBL << "Removed voids: " << endl;
1343 
1344  double spurious_time = omp_get_wtime();
1345  if (initial_radius) {
1346  vector<bool> remove(nObjects(), false);
1347  for (size_t i=0; i<nObjects(); i++)
1348  if (radius(i) < delta_r[0] || delta_r[1] < radius(i)) remove[i] = true;
1349 
1350  remove_objects(remove);
1351  cout << "\t r_min-r_max criterion: " << count(remove.begin(), remove.end(), true) << endl;
1352  }
1353 
1354  compute_centralDensity(tracers_catalogue, ChM, threshold);
1355  coutCBL << "Voids in the Catalogue: " << nObjects() << endl;
1356  if (nObjects()==0) ErrorCBL("Empty void catalogue!", "Catalogue", "VoidCatalogue.cpp");
1357  coutCBL << "Time spent by the spurious voids-checking procedure: " << fixed << setprecision(3) << omp_get_wtime()-spurious_time << " seconds" << endl;
1358 
1359 
1360  // ---------------------------------------------------- //
1361  // ----------------- Radius Rescaling ----------------- //
1362  // ---------------------------------------------------- //
1363 
1364  if (rescale) {
1365 
1366  double density = tracers_catalogue->numdensity();
1367 
1368  cout << endl;
1369  coutCBL << par::col_green << " --- Rescaling radii --- \n" << par::col_default << endl;
1370 
1371  double rescaling_time = omp_get_wtime();
1372 
1373  vector<bool> remove(nObjects(), false), bad_rescale(nObjects(), false);
1374 
1375  double min_X = tracers_catalogue->Min(Var::_X_), max_X = tracers_catalogue->Max(Var::_X_);
1376  double min_Y = tracers_catalogue->Min(Var::_Y_), max_Y = tracers_catalogue->Max(Var::_Y_);
1377  double min_Z = tracers_catalogue->Min(Var::_Z_), max_Z = tracers_catalogue->Max(Var::_Z_);
1378 
1379 #pragma omp parallel num_threads(omp_get_max_threads())
1380  {
1381 #pragma omp for ordered schedule(dynamic)
1382  for (size_t j=0; j<nObjects(); j++) {
1383  vector<double> temp_val = {3.*radius(j), xx(j)-min_X, max_X-xx(j), yy(j)-min_Y, max_Y-yy(j), zz(j)-min_Z, max_Z-zz(j), delta_r[1]};
1384  double value = *min_element(temp_val.begin(), temp_val.end());
1385 #pragma omp ordered
1386  ChM.get_searching_region(value);
1387  vector<long> close = ChM.close_objects(coordinate(j));
1388  vector<double> radii;
1389  for (auto&& k : close) {
1390  double distance = Catalogue::distance(j, tracers_catalogue->catalogue_object(k));
1391  if (distance < value) radii.emplace_back(distance);
1392  }
1393  if (radii.size() < 3) {
1394  remove[j] = true;
1395  bad_rescale[j] = true;
1396  continue;
1397  }
1398 
1399  std::sort(radii.begin(), radii.end());
1400  vector<vector<double>> data(2, vector<double>(radii.size())); //distances, density contrast
1401  for (size_t i=0; i<radii.size()-1; i++) {
1402  data[0][i] = (radii[i]+radii[i+1])*0.5;
1403  data[1][i] = (i+1)/(volume_sphere(data[0][i])*density);
1404  }
1405  data[0][radii.size()-1] = 2*data[0][radii.size()-2] - data[0][radii.size()-3];
1406  data[1][radii.size()-1] = radii.size()/(volume_sphere(data[0][radii.size()-1])*density);
1407 
1408  unsigned int N = std::round(radii.size()*0.25);
1409  bool expand = false;
1410 
1411  if (N>1) {
1412  while (!(data[1][N-2] < threshold && data[1][N-1] > threshold))
1413  if (N>2)
1414  N--;
1415  else {
1416  expand = true;
1417  break;
1418  }
1419  if (expand) {
1420  N = std::round(radii.size()*0.25);
1421  while (!(data[1][N-2] < threshold && data[1][N-1] > threshold))
1422  if (N<radii.size()-2) N++;
1423  else break;
1424  }
1425  }
1426 
1427  double new_radius;
1428  if (N<=2 or N>=radii.size()-2) {
1429  remove[j] = true;
1430  bad_rescale[j] = true;
1431  continue;
1432  }
1433  else
1434  new_radius = interpolated(threshold, {data[1][N-1], data[1][N]}, {data[0][N-1], data[0][N]}, "Linear");
1435 
1436  if (new_radius > 0) set_var(j, Var::_Radius_, new_radius);
1437  else {
1438  remove[j] = true;
1439  bad_rescale[j] = true;
1440  continue;
1441  }
1442 
1443  // -----------------------
1444 
1445  N = std::round(radii.size()*0.25);
1446  expand = false;
1447 
1448  if (N>1) {
1449  while (!(data[1][N-2] < 0.5 && data[1][N-1] > 0.5))
1450  if (N>2)
1451  N--;
1452  else {
1453  expand = true;
1454  break;
1455  }
1456  if (expand) {
1457  N = std::round(radii.size()*0.25);
1458  while (!(data[1][N-2] < 0.5 && data[1][N-1] > 0.5))
1459  if (N<radii.size()-2)
1460  N++;
1461  else
1462  break;
1463  }
1464  }
1465 
1466  double new_generic;
1467  if (N<=2 or N>=radii.size()-2) {
1468  remove[j] = true;
1469  bad_rescale[j] = true;
1470  continue;
1471  }
1472  else
1473  new_generic = interpolated(threshold, {data[1][N-1], data[1][N]}, {data[0][N-1], data[0][N]}, "Linear");
1474 
1475  if (new_generic > 0) set_var(j, Var::_Generic_, new_generic);
1476  else {
1477  remove[j] = true;
1478  bad_rescale[j] = true;
1479  continue;
1480  }
1481  }
1482  }
1483 
1484  remove_objects(remove);
1485  coutCBL << "Removed voids:" << endl;
1486  cout << "\t Bad rescaled: " << count(bad_rescale.begin(), bad_rescale.end(), true) << endl;
1487 
1488  vector<bool> remove_outofrange(nObjects(), false);
1489  for (size_t ii = 0; ii<nObjects(); ii++)
1490  if (radius(ii) < delta_r[0] || delta_r[1] < radius(ii))
1491  remove_outofrange[ii] = true;
1492 
1493  remove_objects(remove_outofrange);
1494  cout << "\t Out of range ["+conv(delta_r[0],par::fDP2)+","+conv(delta_r[1],par::fDP2)+"]: " << count(remove_outofrange.begin(), remove_outofrange.end(), true) << endl;
1495 
1496  compute_densityContrast(tracers_catalogue, ChM, ratio);
1497 
1498  coutCBL << "Time spent by the rescaling procedure: " << omp_get_wtime()-rescaling_time << " seconds" << endl;
1499 
1500  }
1501 
1502  coutCBL << "Voids in the Catalogue: " << nObjects() << endl;
1503  if (nObjects()==0) ErrorCBL("Empty void catalogue!", "Catalogue", "VoidCatalogue.cpp");
1504 
1505  // ---------------------------------------------------- //
1506  // ------------------ Overlap Check ------------------- //
1507  // ---------------------------------------------------- //
1508 
1509  if (checkoverlap) {
1510  cout << endl;
1511  coutCBL << par::col_green << " --- Checking for overlapping voids --- \n" << par::col_default << endl;
1512 
1513  double ol_time = omp_get_wtime();
1514  vector<bool> remove(nObjects(), false);
1515  vector<double> criteriumOrder(nObjects());
1516  string criterium;
1517 
1518  if (ol_criterion == Var::_CentralDensity_) {
1519  criteriumOrder = var(Var::_CentralDensity_);
1520  criterium = "central density";
1521  }
1522  else if (ol_criterion == Var::_DensityContrast_) {
1523  criteriumOrder = var(Var::_DensityContrast_);
1524  criterium = "density contrast";
1525  }
1526 
1527  else ErrorCBL("allowed overlap criteria are '_CentralDensity_' or '_DensityContrast_' .", "Catalogue", "VoidCatalogue.cpp");
1528 
1529  std::vector<int> indices(nObjects(), 0);
1530  std::iota(indices.begin(), indices.end(), 0);
1531  std::sort(indices.begin(), indices.end(), [&](int A, int B)
1532  -> bool {return (ol_criterion == Var::_CentralDensity_) ? criteriumOrder[A] < criteriumOrder[B] : criteriumOrder[A] > criteriumOrder[B];});
1533 
1534  Order(indices);
1535  coutCBL << "Catalogue ordered according to the " << criterium << endl << endl;
1536  coutCBL << "* * * Generating ChainMesh for void centres * * *" << endl;
1537 
1538  chainmesh::ChainMesh3D ChM_voids(Min(Var::_Radius_), var(Var::_X_), var(Var::_Y_), var(Var::_Z_), 2*Max(Var::_Radius_));
1539 
1540  for (size_t i=0; i<nObjects(); i++) {
1541  if (!remove[i]) {
1542  vector<long> close = ChM_voids.close_objects(coordinate(i));
1543  std::sort(close.begin(), close.end());
1544  for (auto && j : close) {
1545  if (!remove[j]) {
1546  double distance = Catalogue::distance(i, catalogue_object(j));
1547  if (distance < radius(i)+radius(j) && (long)i!=j) {
1548  if (ol_criterion == Var::_CentralDensity_) {
1549  if (centralDensity(i) < centralDensity(j)) {
1550  remove[j] = true;
1551  }
1552  else if (centralDensity(i) > centralDensity(j)) {
1553  remove[i] = true;
1554  break;
1555  }
1556  else if (centralDensity(i) == centralDensity(j)) {
1557  if (densityContrast(i) < densityContrast(j)) {
1558  remove[i] = true;
1559  break;
1560  }
1561  else
1562  remove[j] = true;
1563  }
1564  }
1565  else if (ol_criterion == Var::_DensityContrast_) {
1566  if (densityContrast(i) < densityContrast(j)) {
1567  remove[i] = true;
1568  break;
1569  }
1570  else if (densityContrast(i) > densityContrast(j))
1571  remove[j] = true;
1572  else if (densityContrast(i) == densityContrast(j)) {
1573  if (centralDensity(i) < centralDensity(j))
1574  remove[j] = true;
1575  else {
1576  remove[i] = true;
1577  break;
1578  }
1579  }
1580  }
1581  }
1582  }
1583  }
1584  }
1585  }
1586 
1587  remove_objects(remove);
1588 
1589  cout << endl;
1590  coutCBL << "Voids removed to avoid overlap: " << count(remove.begin(), remove.end(), true) << endl;
1591  coutCBL << "Time spent by the overlap-checking procedure: " << omp_get_wtime()-ol_time << " seconds" << endl;
1592  }
1593 
1594  cout << endl;
1595  coutCBL << "Voids in the final Catalogue: " << nObjects() << endl;
1596  if (nObjects()==0) ErrorCBL("Empty void catalogue!", "Catalogue", "VoidCatalogue.cpp");
1597 
1598  cout << endl;
1599  coutCBL << "Total time spent: " << omp_get_wtime()-start_time << " seconds \n" << endl;
1600 
1601 }
1602 
1603 //********************************** VERSION 2.0 **********************************//
1604 
1605 void cbl::catalogue::Catalogue::clean_void_catalogue (const std::vector<double> par_numdensity, const bool initial_radius, const std::vector<double> delta_r, const double threshold, bool rescale, const std::shared_ptr<Catalogue> tracers_catalogue, chainmesh::ChainMesh3D ChM, const double ratio, const bool checkoverlap, const Var ol_criterion)
1606 {
1607  const double start_time = omp_get_wtime();
1608 
1609  // ---------------------------------------------------- //
1610  // ---------------- Cleaning Procedure ---------------- //
1611  // ---------------------------------------------------- //
1612 
1613  cout << endl;
1614  coutCBL << par::col_blue << "> > > > > >>>>>>>>>> CLEANING PROCEDURE STARTED <<<<<<<<< < < < < <" << par::col_default << endl << endl;
1615  coutCBL << "Voids in the initial Catalogue: " << nObjects() << endl << endl;
1616  if (nObjects()==0) ErrorCBL("Empty void catalogue!", "Catalogue", "VoidCatalogue.cpp");
1617 
1618  coutCBL << par::col_green << " --- Removing spurious voids --- " << par::col_default << endl << endl;
1619  coutCBL << "Removed voids: " << endl;
1620 
1621  if (initial_radius) {
1622  vector<bool> remove(nObjects(), false);
1623  for (size_t i=0; i<nObjects(); i++)
1624  if (radius(i) < delta_r[0] || delta_r[1] < radius(i)) remove[i] = true;
1625 
1626  remove_objects(remove);
1627  cout << "\t r_min-r_max criterion: " << count(remove.begin(), remove.end(), true) << endl;
1628  }
1629 
1630  compute_centralDensity(tracers_catalogue, ChM, par_numdensity, threshold);
1631  coutCBL << "Voids in the Catalogue: " << nObjects() << endl;
1632  if (nObjects()==0) ErrorCBL("Empty void catalogue!", "Catalogue", "VoidCatalogue.cpp");
1633 
1634  // ---------------------------------------------------- //
1635  // ----------------- Radius Rescaling ----------------- //
1636  // ---------------------------------------------------- //
1637 
1638  if (rescale) {
1639 
1640  cout << endl;
1641  coutCBL << par::col_green << " --- Rescaling radii --- \n" << par::col_default << endl;
1642 
1643  double rescaling_time = omp_get_wtime();
1644 
1645  vector<bool> remove(nObjects(), false), bad_rescale(nObjects(), false);
1646 
1647 #pragma omp parallel num_threads(omp_get_max_threads())
1648  {
1649 #pragma omp for ordered schedule(dynamic)
1650  for (size_t j=0; j<nObjects(); j++) {
1651  double value = 2*radius(j);
1652 #pragma omp ordered
1653  ChM.get_searching_region(value);
1654  vector<long> close = ChM.close_objects(coordinate(j));
1655  vector<double> radii;
1656  for (auto&& k : close) {
1657  double distance = Catalogue::distance(j, tracers_catalogue->catalogue_object(k));
1658  if (distance < value) radii.emplace_back(distance);
1659  }
1660  if (radii.size() < 3) {
1661  remove[j] = true;
1662  bad_rescale[j] = true;
1663  continue;
1664  }
1665 
1666  double zz = redshift(j);
1667  double density = 0.;
1668  for (size_t N=par_numdensity.size(); N-->0;) density += par_numdensity[par_numdensity.size()-1-N]*pow(zz,N);
1669 
1670  std::sort(radii.begin(), radii.end());
1671 
1672  vector<vector<double>> data(2, vector<double>(radii.size())); //distances, density contrast
1673  for (size_t i=0; i<radii.size()-1; i++) {
1674  data[0][i] = (radii[i]+radii[i+1])*0.5;
1675  data[1][i] = (i+1)/(volume_sphere(data[0][i])*density);
1676  }
1677  data[0][radii.size()-1] = 2*data[0][radii.size()-2] - data[0][radii.size()-3];
1678  data[1][radii.size()-1] = radii.size()/(volume_sphere(data[0][radii.size()-1])*density);
1679 
1680  unsigned int N = std::round(radii.size()*0.25);
1681  bool expand = false;
1682 
1683  if (N>1) {
1684  while (!(data[1][N-2] < threshold && data[1][N-1] > threshold))
1685  if (N>2) N--;
1686  else {
1687  expand = true;
1688  break;
1689  }
1690  if (expand) {
1691  N = std::round(radii.size()*0.25);
1692  while (!(data[1][N-2] < threshold && data[1][N-1] > threshold))
1693  if (N<radii.size()-2) N++;
1694  else break;
1695  }
1696  }
1697 
1698  double new_radius;
1699  if (N<=2 or N>=radii.size()-2) {
1700  remove[j] = true;
1701  bad_rescale[j] = true;
1702  continue;
1703  }
1704  else
1705  new_radius = interpolated(threshold, {data[1][N-1], data[1][N]}, {data[0][N-1], data[0][N]}, "Linear");
1706 
1707  if (new_radius > 0) set_var(j, Var::_Radius_, new_radius);
1708  else {
1709  remove[j] = true;
1710  bad_rescale[j] = true;
1711  continue;
1712  }
1713 
1714  // -----------------------
1715 
1716  N = std::round(radii.size()*0.25);
1717  expand = false;
1718 
1719  if (N>1) {
1720  while (!(data[1][N-2] < 0.5 && data[1][N-1] > 0.5))
1721  if (N>2)
1722  N--;
1723  else {
1724  expand = true;
1725  break;
1726  }
1727  if (expand) {
1728  N = std::round(radii.size()*0.25);
1729  while (!(data[1][N-2] < 0.5 && data[1][N-1] > 0.5))
1730  if (N<radii.size()-2)
1731  N++;
1732  else
1733  break;
1734  }
1735  }
1736 
1737  double new_generic;
1738  if (N<=2 or N>=radii.size()-2) {
1739  remove[j] = true;
1740  bad_rescale[j] = true;
1741  continue;
1742  }
1743  else
1744  new_generic = interpolated(threshold, {data[1][N-1], data[1][N]}, {data[0][N-1], data[0][N]}, "Linear");
1745 
1746  if (new_generic > 0) set_var(j, Var::_Generic_, new_generic);
1747  else {
1748  remove[j] = true;
1749  bad_rescale[j] = true;
1750  continue;
1751  }
1752  }
1753  }
1754 
1755  remove_objects(remove);
1756  coutCBL << "Removed voids:" << endl;
1757  cout << "\t Bad rescaled: " << count(bad_rescale.begin(), bad_rescale.end(), true) << endl;
1758 
1759  vector<bool> remove_outofrange(nObjects(), false);
1760  for (size_t ii = 0; ii<nObjects(); ii++)
1761  if (radius(ii) < delta_r[0] || delta_r[1] < radius(ii))
1762  remove_outofrange[ii] = true;
1763 
1764  remove_objects(remove_outofrange);
1765  cout << "\t Out of range ["+conv(delta_r[0],par::fDP2)+","+conv(delta_r[1],par::fDP2)+"]: " << count(remove_outofrange.begin(), remove_outofrange.end(), true) << endl;
1766 
1767  compute_densityContrast(tracers_catalogue, ChM, ratio);
1768 
1769  coutCBL << "Time spent by the rescaling procedure: " << omp_get_wtime()-rescaling_time << " seconds \n" << endl;
1770 
1771  }
1772 
1773  coutCBL << "Voids in the Catalogue: " << nObjects() << endl;
1774  if (nObjects()==0) ErrorCBL("Empty void catalogue!", "Catalogue", "VoidCatalogue.cpp");
1775 
1776  // ---------------------------------------------------- //
1777  // ------------------ Overlap Check ------------------- //
1778  // ---------------------------------------------------- //
1779 
1780  if (checkoverlap) {
1781  cout << endl;
1782  coutCBL << par::col_green << " --- Checking for overlapping voids --- \n" << par::col_default << endl;
1783 
1784  double ol_time = omp_get_wtime();
1785  vector<bool> remove(nObjects(), false);
1786 
1787  vector<double> criteriumOrder(nObjects());
1788  string criterium;
1789 
1790  if (ol_criterion == Var::_CentralDensity_) {
1791  criteriumOrder = var(Var::_CentralDensity_);
1792  criterium = "central density";
1793  }
1794  else if (ol_criterion == Var::_DensityContrast_) {
1795  criteriumOrder = var(Var::_DensityContrast_);
1796  criterium = "density contrast";
1797  }
1798 
1799  else ErrorCBL("allowed overlap criteria are '_CentralDensity_' or '_DensityContrast_' .", "Catalogue", "VoidCatalogue.cpp");
1800 
1801  std::vector<int> indices(nObjects(), 0);
1802  std::iota(indices.begin(), indices.end(), 0);
1803  std::sort(indices.begin(), indices.end(), [&](int A, int B)
1804  -> bool {return (ol_criterion == Var::_CentralDensity_) ? criteriumOrder[A] < criteriumOrder[B] : criteriumOrder[A] > criteriumOrder[B];});
1805 
1806  Order(indices);
1807  coutCBL << "Catalogue ordered according to the " << criterium << endl << endl;
1808  coutCBL << "* * * Generating ChainMesh for void centres * * *" << endl;
1809 
1810  chainmesh::ChainMesh3D ChM_voids(Min(Var::_Radius_), var(Var::_X_), var(Var::_Y_), var(Var::_Z_), 2*Max(Var::_Radius_));
1811 
1812  for (size_t i=0; i<nObjects(); i++) {
1813  if (!remove[i]) {
1814  vector<long> close = ChM_voids.close_objects(coordinate(i));
1815  std::sort(close.begin(), close.end());
1816  for (auto &&j : close) {
1817  if (!remove[j]) {
1818  double distance = Catalogue::distance(i, catalogue_object(j));
1819  if (distance<radius(i)+radius(j) && (long)i!=j) {
1820  if (ol_criterion == Var::_CentralDensity_) {
1821  if (centralDensity(i) < centralDensity(j)) {
1822  remove[j] = true;
1823  }
1824  else if (centralDensity(i) > centralDensity(j)) {
1825  remove[i] = true;
1826  break;
1827  }
1828  else if (centralDensity(i) == centralDensity(j)) {
1829  if (densityContrast(i) < densityContrast(j)) {
1830  remove[i] = true;
1831  break;
1832  }
1833  else
1834  remove[j] = true;
1835  }
1836  }
1837  else if (ol_criterion == Var::_DensityContrast_) {
1838  if (densityContrast(i) < densityContrast(j)) {
1839  remove[i] = true;
1840  break;
1841  }
1842  else if (densityContrast(i) > densityContrast(j))
1843  remove[j] = true;
1844  else if (densityContrast(i) == densityContrast(j)) {
1845  if (centralDensity(i) < centralDensity(j))
1846  remove[j] = true;
1847  else {
1848  remove[i] = true;
1849  break;
1850  }
1851  }
1852  }
1853  }
1854  }
1855  }
1856  }
1857  }
1858 
1859  remove_objects(remove);
1860 
1861  cout << endl;
1862  coutCBL << "Voids removed to avoid overlap: " << count(remove.begin(), remove.end(), true) << endl;
1863  coutCBL << "Time spent by the overlap-checking procedure: " << omp_get_wtime()-ol_time << " seconds" << endl;
1864  }
1865 
1866  cout << endl;
1867  coutCBL << "Voids in the final Catalogue: " << nObjects() << endl;
1868  if (nObjects()==0) ErrorCBL("Empty void catalogue!", "Catalogue", "VoidCatalogue.cpp");
1869 
1870  cout << endl;
1871  coutCBL << "Total time spent: " << omp_get_wtime()-start_time << " seconds \n" << endl;
1872 
1873 }
1874 
1875 
1876 //********************************** VERSION 3.0 **********************************//
1877 
1878 void cbl::catalogue::Catalogue::clean_void_catalogue (const std::vector<std::vector<double>> data_numdensity, const std::string method_interpolation, const bool initial_radius, const std::vector<double> delta_r, const double threshold, bool rescale, const std::shared_ptr<Catalogue> tracers_catalogue, chainmesh::ChainMesh3D ChM, const double ratio, const bool checkoverlap, const Var ol_criterion)
1879 {
1880 
1881  const double start_time = omp_get_wtime();
1882 
1883  // ---------------------------------------------------- //
1884  // ---------------- Cleaning Procedure ---------------- //
1885  // ---------------------------------------------------- //
1886 
1887  cout << endl;
1888  coutCBL << par::col_blue << "> > > > > >>>>>>>>>> CLEANING PROCEDURE STARTED <<<<<<<<< < < < < <" << par::col_default << endl << endl;
1889  coutCBL << "Voids in the initial Catalogue: " << nObjects() << endl << endl;
1890  if (nObjects()==0) ErrorCBL("Empty void catalogue!", "Catalogue", "VoidCatalogue.cpp");
1891 
1892  coutCBL << par::col_green << " --- Removing spurious voids --- " << par::col_default << endl << endl;
1893  coutCBL << "Removed voids: " << endl;
1894 
1895  if (initial_radius) {
1896  vector<bool> remove(nObjects(), false);
1897  for (size_t i=0; i<nObjects(); i++)
1898  if (radius(i) < delta_r[0] || delta_r[1] < radius(i)) remove[i] = true;
1899 
1900  remove_objects(remove);
1901  cout << "\t r_min-r_max criterion: " << count(remove.begin(), remove.end(), true) << endl;
1902  }
1903 
1904  compute_centralDensity(tracers_catalogue, ChM, data_numdensity, method_interpolation, threshold);
1905  coutCBL << "Voids in the Catalogue: " << nObjects() << endl;
1906  if (nObjects()==0) ErrorCBL("Empty void catalogue!", "Catalogue", "VoidCatalogue.cpp");
1907 
1908  // ---------------------------------------------------- //
1909  // ----------------- Radius Rescaling ----------------- //
1910  // ---------------------------------------------------- //
1911 
1912  if (rescale) {
1913 
1914  cout << endl;
1915  coutCBL << par::col_green << " --- Rescaling radii --- \n" << par::col_default << endl;
1916 
1917  double rescaling_time = omp_get_wtime();
1918 
1919  vector<bool> remove(nObjects(), false), bad_rescale(nObjects(), false);
1920 
1921 #pragma omp parallel num_threads(omp_get_max_threads())
1922  {
1923 #pragma omp for ordered schedule(dynamic)
1924  for (size_t j=0; j<nObjects(); j++) {
1925  double value = 2*radius(j);
1926 #pragma omp ordered
1927  ChM.get_searching_region(value);
1928  vector<long> close = ChM.close_objects(coordinate(j));
1929  vector<double> radii;
1930  for (auto&& k : close) {
1931  double distance = Catalogue::distance(j, tracers_catalogue->catalogue_object(k));
1932  if (distance < value) radii.emplace_back(distance);
1933  }
1934  if (radii.size() < 3) {
1935  remove[j] = true;
1936  bad_rescale[j] = true;
1937  continue;
1938  }
1939 
1940  double zz = redshift(j);
1941  double density = interpolated(zz, data_numdensity[0], data_numdensity[1], method_interpolation);
1942 
1943  std::sort(radii.begin(), radii.end());
1944 
1945  vector<vector<double>> data(2, vector<double>(radii.size())); //distances, density contrast
1946  for (size_t i=0; i<radii.size()-1; i++) {
1947  data[0][i] = (radii[i]+radii[i+1])*0.5;
1948  data[1][i] = (i+1)/(volume_sphere(data[0][i])*density);
1949  }
1950  data[0][radii.size()-1] = 2*data[0][radii.size()-2] - data[0][radii.size()-3];
1951  data[1][radii.size()-1] = radii.size()/(volume_sphere(data[0][radii.size()-1])*density);
1952 
1953  unsigned int N = std::round(radii.size()*0.25);
1954  bool expand = false;
1955 
1956  if (N>1) {
1957  while (!(data[1][N-2] < threshold && data[1][N-1] > threshold))
1958  if (N>2) N--;
1959  else {
1960  expand = true;
1961  break;
1962  }
1963  if (expand) {
1964  N = std::round(radii.size()*0.25);
1965  while (!(data[1][N-2] < threshold && data[1][N-1] > threshold))
1966  if (N<radii.size()-2) N++;
1967  else break;
1968  }
1969  }
1970 
1971  double new_radius;
1972  if (N<=2 or N>=radii.size()-2) {
1973  remove[j] = true;
1974  bad_rescale[j] = true;
1975  continue;
1976  }
1977  else
1978  new_radius = interpolated(threshold, {data[1][N-1], data[1][N]}, {data[0][N-1], data[0][N]}, "Linear");
1979 
1980  if (new_radius > 0) set_var(j, Var::_Radius_, new_radius);
1981  else {
1982  remove[j] = true;
1983  bad_rescale[j] = true;
1984  continue;
1985  }
1986 
1987  // -----------------------
1988 
1989  N = std::round(radii.size()*0.25);
1990  expand = false;
1991 
1992  if (N>1) {
1993  while (!(data[1][N-2] < 0.5 && data[1][N-1] > 0.5))
1994  if (N>2)
1995  N--;
1996  else {
1997  expand = true;
1998  break;
1999  }
2000  if (expand) {
2001  N = std::round(radii.size()*0.25);
2002  while (!(data[1][N-2] < 0.5 && data[1][N-1] > 0.5))
2003  if (N<radii.size()-2)
2004  N++;
2005  else
2006  break;
2007  }
2008  }
2009 
2010  double new_generic;
2011  if (N<=2 or N>=radii.size()-2) {
2012  remove[j] = true;
2013  bad_rescale[j] = true;
2014  continue;
2015  }
2016  else
2017  new_generic = interpolated(threshold, {data[1][N-1], data[1][N]}, {data[0][N-1], data[0][N]}, "Linear");
2018 
2019  if (new_generic > 0) set_var(j, Var::_Generic_, new_generic);
2020  else {
2021  remove[j] = true;
2022  bad_rescale[j] = true;
2023  continue;
2024  }
2025  }
2026  }
2027 
2028  remove_objects(remove);
2029  coutCBL << "Removed voids:" << endl;
2030  cout << "\t Bad rescaled: " << count(bad_rescale.begin(), bad_rescale.end(), true) << endl;
2031 
2032  vector<bool> remove_outofrange(nObjects(), false);
2033  for (size_t ii = 0; ii<nObjects(); ii++)
2034  if (radius(ii) < delta_r[0] || delta_r[1] < radius(ii))
2035  remove_outofrange[ii] = true;
2036 
2037  remove_objects(remove_outofrange);
2038  cout << "\t Out of range ["+conv(delta_r[0],par::fDP2)+","+conv(delta_r[1],par::fDP2)+"]: " << count(remove_outofrange.begin(), remove_outofrange.end(), true) << endl;
2039 
2040  compute_densityContrast(tracers_catalogue, ChM, ratio);
2041 
2042  coutCBL << "Time spent by the rescaling procedure: " << omp_get_wtime()-rescaling_time << " seconds \n" << endl;
2043 
2044  }
2045 
2046  coutCBL << "Voids in the Catalogue: " << nObjects() << endl;
2047  if (nObjects()==0) ErrorCBL("Empty void catalogue!", "Catalogue", "VoidCatalogue.cpp");
2048 
2049  // ---------------------------------------------------- //
2050  // ------------------ Overlap Check ------------------- //
2051  // ---------------------------------------------------- //
2052 
2053  if (checkoverlap) {
2054  cout << endl;
2055  coutCBL << par::col_green << " --- Checking for overlapping voids --- \n" << par::col_default << endl;
2056 
2057  double ol_time = omp_get_wtime();
2058  vector<bool> remove(nObjects(), false);
2059 
2060  vector<double> criteriumOrder(nObjects());
2061  string criterium;
2062 
2063  if (ol_criterion == Var::_CentralDensity_) {
2064  criteriumOrder = var(Var::_CentralDensity_);
2065  criterium = "central density";
2066  }
2067  else if (ol_criterion == Var::_DensityContrast_) {
2068  criteriumOrder = var(Var::_DensityContrast_);
2069  criterium = "density contrast";
2070  }
2071 
2072  else ErrorCBL("allowed overlap criteria are '_CentralDensity_' or '_DensityContrast_' .", "Catalogue", "VoidCatalogue.cpp");
2073 
2074  std::vector<int> indices(nObjects(), 0);
2075  std::iota(indices.begin(), indices.end(), 0);
2076  std::sort(indices.begin(), indices.end(), [&](int A, int B)
2077  -> bool {return (ol_criterion == Var::_CentralDensity_) ? criteriumOrder[A] < criteriumOrder[B] : criteriumOrder[A] > criteriumOrder[B];});
2078 
2079  Order(indices);
2080  coutCBL << "Catalogue ordered according to the " << criterium << endl << endl;
2081  coutCBL << "* * * Generating ChainMesh for void centres * * *" << endl;
2082 
2083  chainmesh::ChainMesh3D ChM_voids(Min(Var::_Radius_), var(Var::_X_), var(Var::_Y_), var(Var::_Z_), 2*Max(Var::_Radius_));
2084 
2085  for (size_t i=0; i<nObjects(); i++) {
2086  if (!remove[i]) {
2087  vector<long> close = ChM_voids.close_objects(coordinate(i));
2088  std::sort(close.begin(), close.end());
2089  for (auto &&j : close) {
2090  if (!remove[j]) {
2091  double distance = Catalogue::distance(i, catalogue_object(j));
2092  if (distance < radius(i)+radius(j) && (long)i!=j) {
2093  if (ol_criterion == Var::_CentralDensity_) {
2094  if (centralDensity(i) < centralDensity(j)) {
2095  remove[j] = true;
2096  }
2097  else if (centralDensity(i) > centralDensity(j)) {
2098  remove[i] = true;
2099  break;
2100  }
2101  else if (centralDensity(i) == centralDensity(j)) {
2102  if (densityContrast(i) < densityContrast(j)) {
2103  remove[i] = true;
2104  break;
2105  }
2106  else
2107  remove[j] = true;
2108  }
2109  }
2110  else if (ol_criterion == Var::_DensityContrast_) {
2111  if (densityContrast(i) < densityContrast(j)) {
2112  remove[i] = true;
2113  break;
2114  }
2115  else if (densityContrast(i) > densityContrast(j))
2116  remove[j] = true;
2117  else if (densityContrast(i) == densityContrast(j)) {
2118  if (centralDensity(i) < centralDensity(j))
2119  remove[j] = true;
2120  else {
2121  remove[i] = true;
2122  break;
2123  }
2124  }
2125  }
2126  }
2127  }
2128  }
2129  }
2130  }
2131 
2132  remove_objects(remove);
2133 
2134  cout << endl;
2135  coutCBL << "Voids removed to avoid overlap: " << count(remove.begin(), remove.end(), true) << endl;
2136  coutCBL << "Time spent by the overlap-checking procedure: " << omp_get_wtime()-ol_time << " seconds" << endl;
2137  }
2138 
2139  cout << endl;
2140  coutCBL << "Voids in the final Catalogue: " << nObjects() << endl;
2141  if (nObjects()==0) ErrorCBL("Empty void catalogue!", "Catalogue", "VoidCatalogue.cpp");
2142 
2143  cout << endl;
2144  coutCBL << "Total time spent: " << omp_get_wtime()-start_time << " seconds \n" << endl;
2145 
2146 }
2147 
2148 // ============================================================================
2149 
2150 void cbl::catalogue::Catalogue::compute_centralDensity (const std::shared_ptr<Catalogue> tracers_catalogue, chainmesh::ChainMesh3D ChM, const double threshold)
2151 {
2152  const double density = tracers_catalogue->numdensity();
2153  const double mps = tracers_catalogue->mps();
2154  vector<bool> remove(nObjects(), false);
2155 
2156  for (size_t j=0; j<nObjects(); j++) {
2157  ChM.get_searching_region(2*mps);
2158  int cont=0;
2159  vector<long> close = ChM.close_objects(coordinate(j));
2160  for (auto&& k : close) {
2161  double distance = cbl::catalogue::Catalogue::distance(j, tracers_catalogue->catalogue_object(k));
2162  if (distance <= 2*mps) cont++;
2163  }
2164 
2165  if (cont/(volume_sphere(2*mps)*density) > threshold) remove[j]=true;
2166  else
2167  set_var(j, Var::_CentralDensity_, cont/(volume_sphere(2*mps)*density));
2168  }
2169  remove_objects(remove);
2170  cout << "\t High central density: " << std::count(remove.begin(), remove.end(), true) << endl;
2171 
2172 }
2173 
2174 // ============================================================================
2175 
2176 void cbl::catalogue::Catalogue::compute_centralDensity (const std::shared_ptr<Catalogue> tracers_catalogue, chainmesh::ChainMesh3D ChM, const std::vector<double> par_numdensity, const double threshold)
2177 {
2178  vector<bool> remove(nObjects(), false);
2179  for (size_t j=0; j<nObjects(); j++) {
2180 
2181  double zz = m_object[j]->redshift();
2182  double density = 0.;
2183  for (size_t N=par_numdensity.size(); N-->0;) density += par_numdensity[par_numdensity.size()-1-N]*pow(zz,N);
2184  double mps = pow(density, -1./3.);
2185 
2186  ChM.get_searching_region(2.*mps);
2187  int cont = 0;
2188  vector<long> close = ChM.close_objects(coordinate(j));
2189  for (auto&& k : close) {
2190  double distance = cbl::catalogue::Catalogue::distance(j, tracers_catalogue->catalogue_object(k));
2191  if (distance <= 2.*mps) cont++;
2192  }
2193 
2194  if (cont/(volume_sphere(2*mps)*density) > threshold) remove[j]=true;
2195  else
2196  set_var(j, Var::_CentralDensity_, cont/(volume_sphere(2*mps)*density));
2197  }
2198  remove_objects(remove);
2199  cout << "\t High central density: " << std::count(remove.begin(), remove.end(), true) << endl;
2200 }
2201 
2202 // ============================================================================
2203 
2204 void cbl::catalogue::Catalogue::compute_centralDensity (const std::shared_ptr<Catalogue> tracers_catalogue, chainmesh::ChainMesh3D ChM, const std::vector<std::vector<double>> data_numdensity, const std::string method_interpolation, const double threshold)
2205 {
2206  vector<bool> remove(nObjects(), false);
2207  for (size_t j=0; j<nObjects(); j++) {
2208 
2209  double zz = m_object[j]->redshift();
2210  double density = interpolated(zz, data_numdensity[0], data_numdensity[1], method_interpolation);
2211  double mps = pow(density, -1./3.);
2212 
2213  ChM.get_searching_region(2.*mps);
2214  int cont = 0;
2215  vector<long> close = ChM.close_objects(coordinate(j));
2216  for (auto&& k : close) {
2217  double distance = cbl::catalogue::Catalogue::distance(j, tracers_catalogue->catalogue_object(k));
2218  if (distance <= 2.*mps) cont++;
2219  }
2220 
2221  if (cont/(volume_sphere(2*mps)*density) > threshold) remove[j]=true;
2222  else
2223  set_var(j, Var::_CentralDensity_, cont/(volume_sphere(2*mps)*density));
2224  }
2225  remove_objects(remove);
2226  cout << "\t High central density: " << std::count(remove.begin(), remove.end(), true) << endl;
2227 
2228 }
2229 
2230 // ============================================================================
2231 
2232 
2233 void cbl::catalogue::Catalogue::compute_densityContrast (const shared_ptr<Catalogue> tracers_catalogue, cbl::chainmesh::ChainMesh3D ChM, const double ratio) {
2234 
2235  if (ratio==0. or ratio==1. or ratio>3.)
2236  ErrorCBL ("Please choose an appropriate value for the parameter ratio!", "compute_densityContrast", "VoidCatalogue.cpp");
2237 
2238  vector<bool> remove(nObjects(), false), void_voids(nObjects(), false), cloud_in_void(nObjects(), false);
2239 
2240 #pragma omp parallel num_threads(omp_get_max_threads())
2241  {
2242 #pragma omp for ordered schedule(static)
2243  for (size_t j = 0; j<nObjects(); j++) {
2244 
2245  double outer_radius = std::max(radius(j)*ratio, radius(j));
2246  double inner_radius = std::min(radius(j)*ratio, radius(j));
2247 
2248 #pragma omp ordered
2249  ChM.get_searching_region(outer_radius);
2250  vector<long> close = ChM.close_objects(coordinate(j));
2251 
2252  //compute distances between the void and the surrounding particles
2253  vector<double> distances;
2254  for (auto&& k : close) {
2255  double distance = Catalogue::distance(j, tracers_catalogue->catalogue_object(k));
2256  if (distance <= outer_radius) distances.emplace_back(distance);
2257  }
2258 
2259  std::sort(distances.begin(), distances.end());
2260 
2261  if (distances.size() > 1) {
2262 
2263  //FIND CENTRAL DENSITY
2264 
2265  int NN_in = 0;
2266  int NN_out = 0;
2267  while (distances[NN_in]<inner_radius && NN_in<(int)distances.size()-1) NN_in++; // check-1
2268  while (distances[NN_out]<outer_radius && NN_out<(int)distances.size()-1) NN_out++;
2269 
2270  double delta_in = NN_in/cbl::volume_sphere(inner_radius);
2271  double delta_out = NN_out/cbl::volume_sphere(outer_radius);
2272 
2273  if (NN_in>0) {
2274  double compare = delta_out/delta_in;
2275  if (compare==0.) {
2276  remove[j] = true;
2277  void_voids[j] = true;
2278  }
2279  else if (compare>=1.)
2280  set_var(j, Var::_DensityContrast_, compare);
2281 
2282  else if (compare<1.) {
2283  remove[j] = true;
2284  cloud_in_void[j] = true;
2285  }
2286  }
2287  else set_var(j, Var::_DensityContrast_, 9999.);
2288 
2289  }
2290  }//for
2291  }
2292  remove_objects(remove);
2293 
2294  cout << "\t Cloud-in-void: " << std::count(cloud_in_void.begin(), cloud_in_void.end(), true) << endl;
2295  cout << "\t Empty voids: " << std::count(cloud_in_void.begin(), cloud_in_void.end(), true) << endl;
2296 
2297 }
CatalogueChainMesh.h
The class catalogue::CatalogueChainMesh.
Catalogue.h
The class Catalogue
ChainMesh_Catalogue.h
Implementation of the chain-mesh data structure.
Data1D.h
The class Data1D.
Func.h
Useful generic functions.
coutCBL
#define coutCBL
CBL print message.
Definition: Kernel.h:734
Object.h
The class Object.
Posterior.h
The class Posterior.
TwoPointCorrelation1D_monopole.h
The class TwoPointCorrelation1D_monopole.
TwoPointCorrelation.h
The class TwoPointCorrelation.
cbl::catalogue::Catalogue
The class Catalogue.
Definition: Catalogue.h:654
cbl::catalogue::Catalogue::Catalogue
Catalogue()=default
default constructor
cbl::catalogue::Catalogue::compute_centralDensity
void compute_centralDensity(const std::shared_ptr< Catalogue > tracers_catalogue, chainmesh::ChainMesh3D ChM, const double threshold)
compute the central density of each object in a void catalogue. We define as central density the dens...
Definition: VoidCatalogue.cpp:2150
cbl::catalogue::Catalogue::nObjects
size_t nObjects() const
get the number of objects of the catalogue
Definition: Catalogue.h:2264
cbl::catalogue::Catalogue::compute_densityContrast
void compute_densityContrast(const std::shared_ptr< Catalogue > tracers_catalogue, chainmesh::ChainMesh3D ChM, const double ratio=0.1)
compute density contrast of cosmic voids in catalogue as the ratio between an inner and an outer sphe...
Definition: VoidCatalogue.cpp:2233
cbl::catalogue::Catalogue::equalize_random_box
void equalize_random_box(cbl::catalogue::Catalogue tracer_catalogue, const int seed)
equalize the number of objects in two Box catalogues
Definition: RandomCatalogue.cpp:694
cbl::catalogue::Catalogue::clean_void_catalogue
void clean_void_catalogue(const bool initial_radius=false, const std::vector< double > delta_r={-1, 1000}, const double threshold=0.205, const bool rescale=true, const std::shared_ptr< Catalogue > tracers_catalogue={}, chainmesh::ChainMesh3D ChM={}, const double ratio=1.5, const bool checkoverlap=true, const Var ol_criterion=Var::_CentralDensity_)
function that modifies a void catalogue according to a set of user selected criteria....
Definition: VoidCatalogue.cpp:1328
cbl::catalogue::Catalogue::distance
double distance(const int i, std::shared_ptr< Object > obj) const
get the distrance between the i-th object of the catalogue and another object
Definition: Catalogue.cpp:2103
cbl::chainmesh::ChainMesh3D
The class ChainMesh3D.
Definition: ChainMesh.h:392
cbl::chainmesh::ChainMesh::close_objects
std::vector< long > close_objects(std::vector< double > center, long ii=-1) const
get the indeces of the objects close to an object
Definition: ChainMesh.cpp:298
cbl::chainmesh::ChainMesh::get_searching_region
void get_searching_region(const double r_max, const double r_min=-1.)
set the internal variable m_search_region, the list of cell around a generic center
Definition: ChainMesh.cpp:226
cbl::cosmology::Cosmology
The class Cosmology.
Definition: Cosmology.h:277
cbl::cosmology::Cosmology::Omega_baryon
double Omega_baryon() const
get the private member Cosmology::m_Omega_baryon
Definition: Cosmology.h:1126
cbl::random::UniformRandomNumbers
The class UniformRandomNumbers.
Definition: RandomNumbers.h:254
cbl::par::col_green
static const std::string col_green
green colour (used when printing something on the screen)
Definition: Constants.h:309
cbl::par::col_default
static const std::string col_default
default colour (used when printing something on the screen)
Definition: Constants.h:297
cbl::par::col_blue
static const std::string col_blue
blue colour (used when printing something on the screen)
Definition: Constants.h:312
cbl::par::fDP2
static const char fDP2[]
conversion symbol for: double -> std::string
Definition: Constants.h:133
cbl::par::pi
static const double pi
: the ratio of a circle's circumference to its diameter
Definition: Constants.h:199
cbl::catalogue::Var
Var
the catalogue variables
Definition: Catalogue.h:70
cbl::catalogue::VoidAlgorithm
VoidAlgorithm
the algorithm used to look for Voids
Definition: Catalogue.h:452
cbl
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38
cbl::Min
T Min(const std::vector< T > vect)
minimum element of a std::vector
Definition: Kernel.h:1324
cbl::conv
std::string conv(const T val, const char *fact)
convert a number to a std::string
Definition: Kernel.h:903
cbl::Euclidean_distance
double Euclidean_distance(const double x1, const double x2, const double y1, const double y2, const double z1, const double z2)
the Euclidean distance in 3D relative to the origin (0,0,0), i.e. the Euclidean norm
Definition: Func.cpp:250
cbl::Print
void Print(const T value, const int prec, const int ww, const std::string header="", const std::string end="\n", const bool use_coutCBL=true, std::ostream &stream=std::cout, const std::string colour=cbl::par::col_default)
function to print values with a proper homegenised format
Definition: Kernel.h:1142
cbl::volume_sphere
T volume_sphere(const T RR)
the volume of a sphere of a given radius
Definition: Func.h:1231
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
cbl::Max
T Max(const std::vector< T > vect)
maximum element of a std::vector
Definition: Kernel.h:1336
cbl::interpolated
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