CosmoBolognaLib
Free Software C++/Python libraries for cosmological calculations
CatalogueChainMesh.cpp
Go to the documentation of this file.
1 /********************************************************************
2  * Copyright (C) 2022 by Federico Marulli, Simone Sartori *
3  * federico.marulli3@unibo.it simone.sartori5@studio.unibo.it *
4  * *
5  * This program is free software; you can redistribute it and/or *
6  * modify it under the terms of the GNU General Public License as *
7  * published by the Free Software Foundation; either version 2 of *
8  * the License, or (at your option) any later version. *
9  * *
10  * This program is distributed in the hope that it will be useful, *
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13  * GNU General Public License for more details. *
14  * *
15  * You should have received a copy of the GNU General Public *
16  * License along with this program; if not, write to the Free *
17  * Software Foundation, Inc., *
18  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19  ********************************************************************/
20 
34 #include "CatalogueChainMesh.h"
35 
36 using namespace std;
37 using namespace cbl;
38 
39 
40 // ============================================================================
41 
42 cbl::catalogue::CatalogueChainMesh::CatalogueChainMesh (std::vector<catalogue::Var> variables, double cellsize, std::shared_ptr<cbl::catalogue::Catalogue> cat, std::shared_ptr<cbl::catalogue::Catalogue> cat2)
43 {
44  m_part_catalogue = cat;
45  m_part_catalogue2 = cat2;
46  unsigned int nDim=variables.size();
47  vector<vector<double>> data(nDim);
48  for (size_t i=0; i<nDim; i++) data[i] = m_part_catalogue->var(variables[i]);
49  vector<vector<double>> data2={};
50  if (m_part_catalogue2!=NULL) {
51  data2.resize(nDim, vector<double>(m_part_catalogue2->var(variables[0])));
52  for (size_t i=0; i<nDim; i++) data2[i] = m_part_catalogue2->var(variables[i]);
53  }
54 
55  vector<unsigned int> num_cells(nDim);
56  vector<double> delta(nDim);
57  double start_cellsize = cellsize;
58  m_lim.resize(nDim);
59  for (size_t i=0; i<nDim; i++) {
60  m_lim[i].resize(2);
61  m_lim[i][0] = *min_element(data[i].begin(), data[i].end());
62  m_lim[i][0] = (m_part_catalogue2!=NULL) ? min(m_lim[i][0], *min_element(data2[i].begin(), data2[i].end()))-0.05*cellsize : m_lim[i][0]-0.05*cellsize;
63  m_lim[i][1] = *max_element(data[i].begin(), data[i].end());
64  m_lim[i][1] = (m_part_catalogue2!=NULL) ? max(m_lim[i][1], *max_element(data2[i].begin(), data2[i].end()))+0.05*cellsize : m_lim[i][1]+0.05*cellsize;
65  delta[i] = m_lim[i][1]-m_lim[i][0];
66  }
67 
68  unsigned int nCells=0;
69  double delta_max = *max_element(delta.begin(),delta.end());
70  cellsize = delta_max;
71  vector<vector<unsigned int>> cells(1);
72 
73  while(cellsize > start_cellsize) {
74  nCells++;
75  cellsize = delta_max/nCells;
76  cells.resize((int)pow(nCells, nDim));
77  }
78 
79  m_dimension = nCells;
80  m_cellsize = cellsize;
81 
82  for (unsigned int i=0; i<data[0].size(); i++) {
83  unsigned int index=0;
84  for (unsigned int j=0; j<nDim; j++) index += ((unsigned int)((data[j][i] - m_lim[j][0])/cellsize))*pow(nCells,nDim-j-1);
85  cells[index].emplace_back(i);
86  }
87 
88 
89  unsigned int dim_nearCells = nCells*sqrt(nDim)+1;
90  // include the objects in the catalogue
91 
92  for (size_t i=0; i<cells.size(); i++) {
93  vector<vector<unsigned int>> nearCells(dim_nearCells);
94  for (size_t j=0; j<cells.size(); j++) {
95  double sum = 0;
96  unsigned int index1 = i, index2 = j;
97  for (unsigned int k=0; k<nDim; k++) {
98  double temp_dist = (double)(index1%nCells)-(double)(index2%nCells);
99  if ((int)temp_dist != 0) temp_dist -= 0.5;
100  sum += temp_dist*temp_dist;
101  index1 = (int)(index1/nCells);
102  index2 = (int)(index2/nCells);
103  }
104  nearCells[(int)sqrt(sum)].emplace_back(j);
105  }
106  check_memory(8.0, true, "cbl::catalogue::Catalogue::Catalogue of ChainMeshCatalogue.cpp. Increase cellsize");
107  add_object(move(Object::Create(i, cells[i], nearCells)));
108  }
109 }
110 
111 // =========================================================================
112 
113 std::vector<unsigned int> cbl::catalogue::CatalogueChainMesh::Closer_object (const std::vector<double> pos)
114 {
115  vector<int> inds(3);
116  for (int i=0; i<3; i++) inds[i] = (int)((pos[i] - m_lim[i][0])/m_cellsize);
117 
118  int center_index = inds[2] + inds[1]*m_dimension + inds[0]*m_dimension*m_dimension;
119 
120  unsigned int n_start_max = 1;
121  unsigned int n_start_min = 0;
122  unsigned int n_start = n_start_max;
123  vector<unsigned int> cObject(2);
124  double distance = (unsigned int)m_cellsize*10000;
125  vector<vector<unsigned int>> nCells = nearCells(center_index);
126  double radius = m_cellsize;
127 
128  int add = 3;
129 
130  while(distance == (unsigned int)m_cellsize*10000) {
131  for(unsigned int i = n_start_min; i<n_start_max; i++) {
132  if (i == nCells.size()) break;
133  for(unsigned int j : nCells[i]) {
134  for (unsigned int k : part(j)) {
135  double new_distance = Euclidean_distance(pos[0], m_part_catalogue->xx(k), pos[1], m_part_catalogue->yy(k), pos[2], m_part_catalogue->zz(k));
136  if(new_distance < distance && new_distance < radius && new_distance > n_start_min*m_cellsize) {
137  distance = new_distance;
138  cObject[0] = k;
139  cObject[1] = j;
140  }
141  }
142  }
143  }
144  n_start_max += 1;
145  n_start_min = (n_start_max <= n_start+add) ? 0 : n_start_min+1;
146  radius += m_cellsize;
147  }
148 
149  return cObject;
150 }
151 
152 // =========================================================================
153 
154 std::vector<unsigned int> cbl::catalogue::CatalogueChainMesh::Close_objects (const std::vector<double> pos, const double rMax, const double rMin)
155 {
156  vector<int> inds(3);
157  for (int i=0; i<3; i++) inds[i] = (int)((pos[i] - m_lim[i][0])/m_cellsize);
158  int center_index = inds[2] + inds[1]*m_dimension + inds[0]*m_dimension*m_dimension;
159 
160  unsigned int n_max = (rMax/m_cellsize)+1;
161  n_max = (m_cellsize*n_max < rMax) ? n_max+2 : n_max+1;
162  unsigned int n_min = max(nint(rMin/m_cellsize)-3, 0);
163 
164  vector<vector<unsigned int>> nCells = nearCells(center_index);
165  vector<unsigned int> objects;
166 
167  for(unsigned int i = n_min; i<n_max; i++) {
168  if (i == nCells.size()) break;
169  for(unsigned int j : nCells[i]) {
170  for (unsigned int k : part(j)) {
171  double new_distance = Euclidean_distance(pos[0], m_part_catalogue->xx(k), pos[1], m_part_catalogue->yy(k), pos[2], m_part_catalogue->zz(k));
172  if(new_distance < rMax && new_distance >= rMin) objects.emplace_back(k);
173  }
174  }
175  }
176 
177  return objects;
178 }
179 
180 // =========================================================================
181 
182 unsigned int cbl::catalogue::CatalogueChainMesh::Count_objects (const std::vector<double> pos, const double rMax, const double rMin)
183 {
184  vector<int> inds(3);
185  for (int i=0; i<3; i++) inds[i] = (int)((pos[i] - m_lim[i][0])/m_cellsize);
186  int center_index = inds[2] + inds[1]*m_dimension + inds[0]*m_dimension*m_dimension;
187  int add = 3;
188  unsigned int n_max = nint(rMax/m_cellsize);
189  n_max = (m_cellsize*n_max < rMax) ? n_max+2 : n_max+1;
190  unsigned int n_min = max(nint(rMin/m_cellsize)-add, 0);
191 
192  vector<vector<unsigned int>> nCells = nearCells(center_index);
193  unsigned int objects = 0;
194 
195  for(unsigned int i = n_min; i<n_max; i++) {
196  if (i == nCells.size()) break;
197  for(unsigned int j : nCells[i]) {
198  for (unsigned int k : part(j)) {
199  double new_distance = Euclidean_distance(pos[0], m_part_catalogue->xx(k), pos[1], m_part_catalogue->yy(k), pos[2], m_part_catalogue->zz(k));
200  if(new_distance < rMax && new_distance >= rMin)
201  objects+=1;
202 
203  }
204  }
205  }
206 
207  return objects;
208 }
209 
210 // =========================================================================
211 
212 std::vector<unsigned int> cbl::catalogue::CatalogueChainMesh::N_nearest_objects (const std::vector<double> pos, const unsigned int N)
213 {
214  vector<int> inds(3);
215  for (int i=0; i<3; i++) inds[i] = (int)((pos[i] - m_lim[i][0])/m_cellsize);
216 
217  int center_index = inds[2] + inds[1]*m_dimension + inds[0]*m_dimension*m_dimension;
218 
219  vector<unsigned int> Object={}, temp_Object={};
220  vector<double> distance_obj={}, temp_distance_obj={};;
221  vector<vector<unsigned int>> nCells = nearCells(center_index);
222  double radius=m_cellsize, new_distance;
223  size_t index=0;
224  vector<bool> mask;
225  while (Object.size() < N) {
226 
227  for (size_t i=0; i<temp_Object.size(); i++) if (temp_distance_obj[i]<radius*(index+1) && mask[i]==false) mask[i] = true;
228 
229  for(unsigned int j : nCells[index]) {
230  for (unsigned int k : part(j)) {
231  new_distance = Euclidean_distance(pos[0], m_part_catalogue->xx(k), pos[1], m_part_catalogue->yy(k), pos[2], m_part_catalogue->zz(k));
232  temp_Object.push_back(k);
233  temp_distance_obj.push_back(new_distance);
234  if(new_distance < radius*(index+1)) mask.emplace_back(true);
235  else mask.emplace_back(false);
236  }
237  }
238 
239  if (index==nCells.size()-1) {
240  Object = temp_Object;
241  distance_obj = temp_distance_obj;
242  break;
243  }
244  else {
245  unsigned int Ntrue=count(mask.begin(), mask.end(), true);
246  if (Ntrue>=N) {
247  for (size_t i=0; i<temp_Object.size(); i++) {
248  if (mask[i]==true) {
249  Object.emplace_back(temp_Object[i]);
250  distance_obj.emplace_back(temp_distance_obj[i]);
251  }
252  }
253  }
254  }
255 
256  index++;
257  }
258 
259  vector<int> ind(Object.size(), 0);
260  vector<unsigned int> output(N);
261  std::iota(ind.begin(), ind.end(), 0);
262  std::sort(ind.begin(), ind.end(), [&](int A, int B) -> bool {return distance_obj[A] < distance_obj[B];});
263  for (size_t i=0; i<N; i++) output[i] = Object[ind[i]];
264 
265  return output;
266 }
267 
268 // =========================================================================
269 
270 std::vector<unsigned int> cbl::catalogue::CatalogueChainMesh::N_nearest_objects (const unsigned int obj, const unsigned int N)
271 {
272  vector<int> inds(3);
273  vector<double> pos = {m_part_catalogue->xx(obj), m_part_catalogue->yy(obj), m_part_catalogue->zz(obj)};
274  for (int i=0; i<3; i++) inds[i] = (int)((pos[i] - m_lim[i][0])/m_cellsize);
275  int center_index = inds[2] + inds[1]*m_dimension + inds[0]*m_dimension*m_dimension;
276 
277  vector<unsigned int> Object={}, temp_Object={};
278  vector<double> distance_obj={}, temp_distance_obj={};;
279  vector<vector<unsigned int>> nCells = nearCells(center_index);
280  double radius=m_cellsize, new_distance;
281  size_t index=0;
282  vector<bool> mask;
283  while (Object.size() < N) {
284 
285  for (size_t i=0; i<temp_Object.size(); i++) if (temp_distance_obj[i]<radius*(index+1) && mask[i]==false) mask[i] = true;
286 
287  for(unsigned int j : nCells[index]) {
288  for (unsigned int k : part(j)) {
289  new_distance = Euclidean_distance(pos[0], m_part_catalogue->xx(k), pos[1], m_part_catalogue->yy(k), pos[2], m_part_catalogue->zz(k));
290  temp_Object.push_back(k);
291  temp_distance_obj.push_back(new_distance);
292  if(new_distance < radius*(index+1)) mask.emplace_back(true);
293  else mask.emplace_back(false);
294  }
295  }
296 
297  if (index==nCells.size()-1) {
298  Object = temp_Object;
299  distance_obj = temp_distance_obj;
300  break;
301  }
302  else {
303  unsigned int Ntrue=count(mask.begin(), mask.end(), true);
304  if (Ntrue>=N) {
305  for (size_t i=0; i<temp_Object.size(); i++) {
306  if (mask[i]==true) {
307  Object.emplace_back(temp_Object[i]);
308  distance_obj.emplace_back(temp_distance_obj[i]);
309  }
310  }
311  }
312  }
313 
314  index++;
315  }
316 
317  vector<int> ind(Object.size(), 0);
318  vector<unsigned int> output(N);
319  std::iota(ind.begin(), ind.end(), 0);
320  std::sort(ind.begin(), ind.end(), [&](int A, int B) -> bool {return distance_obj[A] < distance_obj[B];});
321  for (size_t i=0; i<N; i++) output[i] = Object[ind[i]];
322 
323  return output;
324 }
325 
326 // =========================================================================
327 
328 std::vector<std::vector<unsigned int>> cbl::catalogue::CatalogueChainMesh::N_nearest_objects_cat (const unsigned int N)
329 {
330  vector<vector<unsigned int>> near_part_cat(m_part_catalogue->nObjects());
331 #pragma omp parallel num_threads(omp_get_max_threads())
332  {
333 #pragma omp for schedule(static)
334  for (size_t i=0; i<m_part_catalogue->nObjects(); i++) {
335  near_part_cat[i] = N_nearest_objects(i, N);
336  }
337  }
338 
339  return near_part_cat;
340 }
341 
342 // =========================================================================
343 
344 void cbl::catalogue::CatalogueChainMesh::deletePart (const unsigned int index)
345 {
346  vector<int> inds(3);
347  vector<double> pos = {m_part_catalogue -> xx(index), m_part_catalogue -> yy(index), m_part_catalogue -> zz(index)};
348 
349  for (int i=0; i<3; i++) inds[i] = (int)((pos[i] - m_lim[i][0])/m_cellsize);
350 
351  int center_index = inds[2] + inds[1]*m_dimension + inds[0]*m_dimension*m_dimension;
352  vector<unsigned int> particles = part(center_index);
353  particles.erase(remove(particles.begin(), particles.end(), index), particles.end());
354 
355  set_part(center_index, particles);
356 
357 }
358 
359 // =========================================================================
360 
361 void cbl::catalogue::CatalogueChainMesh::deletePart (const unsigned int i, const unsigned int p)
362 {
363 
364  vector<unsigned int> particles = part(i);
365  particles.erase(remove(particles.begin(), particles.end(), p), particles.end());
366  set_part(i, particles);
367 
368 }
CatalogueChainMesh.h
The class catalogue::CatalogueChainMesh.
cbl
The global namespace of the CosmoBolognaLib
Definition: CAMB.h:38
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::nint
int nint(const T val)
the nearest integer
Definition: Kernel.h:915
cbl::check_memory
int check_memory(const double frac, const bool exit=true, const std::string func="", const int type=1)
check if the memory used by current process is larger than a given fraction of the available memory
Definition: Kernel.cpp:252