Compadre 1.5.5
Loading...
Searching...
No Matches
NeighborSearchTest.cpp
Go to the documentation of this file.
1#include <vector>
2#include <map>
3#include <random>
4
5#include <Compadre_Config.h>
6#include <Compadre_GMLS.hpp>
8
9#ifdef COMPADRE_USE_MPI
10#include <mpi.h>
11#endif
12
13template<typename T1, typename T2>
14std::pair<T2,T1> swap_pair(const std::pair<T1,T2> &item)
15{
16 return std::pair<T2,T1>(item.second, item.first);
17}
18
19template<typename T1, typename T2>
20std::multimap<T2,T1> invert_map(const std::map<T1,T2> &original_map)
21{
22 std::multimap<T2,T1> new_map;
23 std::transform(original_map.begin(), original_map.end(), std::inserter(new_map, new_map.begin()), swap_pair<T1,T2>);
24 return new_map;
25}
26
27using namespace Compadre;
28
29int main (int argc, char* args[]) {
30
31// initializes MPI (if available) with command line arguments given
32#ifdef COMPADRE_USE_MPI
33MPI_Init(&argc, &args);
34#endif
35
36// initializes Kokkos with command line arguments given
37Kokkos::initialize(argc, args);
38
39// becomes false if the computed solution not within the failure_threshold of the actual solution
40bool all_passed = true;
41
42{
43 // seed random generator
44 srand(1234);
45 Kokkos::Timer timer;
46
47 int search_type = 0; // 0 - radius, 1 - knn
48 if (argc >= 5) {
49 int arg5toi = atoi(args[4]);
50 if (arg5toi >= 0) {
51 search_type = arg5toi;
52 }
53 }
54 bool do_radius_search = (search_type==0);
55 bool do_knn_search = (search_type==1);
56 printf("do_radius_search: %d\n", do_radius_search);
57 printf("do_knn_search: %d\n", do_knn_search);
58
59 // check if 4 arguments are given from the command line
60 // multiplier times h spacing for search
61 double h_multiplier = 3.0; // dimension 3 by default
62 if (argc >= 4) {
63 double arg4tof = atof(args[3]);
64 if (arg4tof > 0) {
65 h_multiplier = arg4tof;
66 }
67 }
68
69 // check if 3 arguments are given from the command line
70 // set the number of target sites where we will reconstruct the target functionals at
71 int number_target_coords = 200; // 200 target sites by default
72 if (argc >= 3) {
73 int arg3toi = atoi(args[2]);
74 if (arg3toi > 0) {
75 number_target_coords = arg3toi;
76 }
77 }
78
79 // check if 2 arguments are given from the command line
80 // set the number of dimensions for points and the search
81 int dimension = 3; // 3D by default
82 if (argc >= 2) {
83 int arg2toi = atoi(args[1]);
84 if (arg2toi > 0) {
85 dimension = arg2toi;
86 }
87 }
88
89 //// minimum neighbors for unisolvency is the same as the size of the polynomial basis
90 //const int min_neighbors = Compadre::GMLS::getNP(order, dimension);
91
92 // approximate spacing of source sites
93 double h_spacing = 1./std::pow(number_target_coords, 0.50 + 0.2*(dimension==2));
94 int n_neg1_to_1 = 2*(1/h_spacing) + 1; // always odd
95
96 // number of source coordinate sites that will fill a box of [-1,1]x[-1,1]x[-1,1] with a spacing approximately h
97 int number_source_coords = (n_neg1_to_1+1)*(n_neg1_to_1+1);
98 if (dimension==3) number_source_coords *= (n_neg1_to_1+1);
99 printf("target coords: %d\n", number_target_coords);
100 printf("source coords: %d\n", number_source_coords);
101
102 // coordinates of source sites
103 Kokkos::View<double**, Kokkos::DefaultHostExecutionSpace> source_coords("source coordinates",
104 number_source_coords, 3);
105
106 // coordinates of target sites
107 Kokkos::View<double**, Kokkos::DefaultHostExecutionSpace> target_coords("target coordinates", number_target_coords, 3);
108
109 // fill source coordinates with a uniform grid
110 int source_index = 0;
111 double this_coord[3] = {0,0,0};
112 for (int i=-n_neg1_to_1/2; i<n_neg1_to_1/2+1; ++i) {
113 this_coord[0] = i*h_spacing;
114 for (int j=-n_neg1_to_1/2; j<n_neg1_to_1/2+1; ++j) {
115 this_coord[1] = j*h_spacing;
116 for (int k=-n_neg1_to_1/2; k<n_neg1_to_1/2+1; ++k) {
117 this_coord[2] = k*h_spacing;
118 if (dimension==3) {
119 source_coords(source_index,0) = this_coord[0];
120 source_coords(source_index,1) = this_coord[1];
121 source_coords(source_index,2) = this_coord[2];
122 source_index++;
123 }
124 }
125 if (dimension==2) {
126 source_coords(source_index,0) = this_coord[0];
127 source_coords(source_index,1) = this_coord[1];
128 source_coords(source_index,2) = 0;
129 source_index++;
130 }
131 }
132 if (dimension==1) {
133 source_coords(source_index,0) = this_coord[0];
134 source_coords(source_index,1) = 0;
135 source_coords(source_index,2) = 0;
136 source_index++;
137 }
138 }
139
140 // fill target coords somewhere inside of [-0.5,0.5]x[-0.5,0.5]x[-0.5,0.5]
141 for(int i=0; i<number_target_coords; i++){
142
143 // first, we get a uniformly random distributed direction
144 double rand_dir[3] = {0,0,0};
145
146 for (int j=0; j<dimension; ++j) {
147 // rand_dir[j] is in [-0.5, 0.5]
148 rand_dir[j] = ((double)rand() / (double) RAND_MAX) - 0.5;
149 }
150
151 // then we get a uniformly random radius
152 for (int j=0; j<dimension; ++j) {
153 target_coords(i,j) = rand_dir[j];
154 }
155
156 }
157
158 if (do_radius_search) {
159 timer.reset();
160 { // 1D compressed row neighbor lists
161 printf("\n1D compressed row neighbor lists:\n");
162 // Point cloud construction for neighbor search
163 // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
164 auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
165
166 Kokkos::View<int*, Kokkos::DefaultHostExecutionSpace> neighbor_lists("neighbor lists",
167 number_target_coords); // first column is # of neighbors
168 Kokkos::View<int*, Kokkos::DefaultHostExecutionSpace> number_neighbors_list("number of neighbors list",
169 number_target_coords); // first column is # of neighbors
170
171 // each target site has a window size
172 Kokkos::View<double*, Kokkos::DefaultHostExecutionSpace> epsilon("h supports", number_target_coords);
173 Kokkos::View<double*, Kokkos::DefaultHostExecutionSpace> random_vec("random vector", number_target_coords);
174
175 Kokkos::parallel_for("random perturbation of search size", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
176 (0,number_target_coords), [&](const int i) {
177 // value in [-.25, .25]
178 random_vec(i) = 0.5 * ((double)rand() / (double) RAND_MAX) - 0.5;
179 epsilon(i) = (h_multiplier + random_vec(i)) * h_spacing;
180 });
181 Kokkos::fence();
182
183 // query the point cloud to generate the neighbor lists using a radius search
184
185 // start with a dry-run, but without enough room to store the results
186 size_t total_num_neighbors = point_cloud_search.generateCRNeighborListsFromRadiusSearch(true /* dry run */,
187 target_coords, neighbor_lists, number_neighbors_list, epsilon);
188 printf("total num neighbors: %lu\n", total_num_neighbors);
189
190 // resize neighbor lists to be large enough to hold the results
191 Kokkos::resize(neighbor_lists, total_num_neighbors);
192
193 // search again, now that we know that there is enough room to store the results
194 point_cloud_search.generateCRNeighborListsFromRadiusSearch(false /* dry run */,
195 target_coords, neighbor_lists, number_neighbors_list, epsilon);
196
197 auto nla(CreateNeighborLists(neighbor_lists, number_neighbors_list));
198
199 double radius_search_time = timer.seconds();
200 printf("nanoflann search time: %f s\n", radius_search_time);
201
202 // convert point cloud search to vector of maps
203 timer.reset();
204 std::vector<std::map<int, double> > point_cloud_neighbor_list(number_target_coords, std::map<int, double>());
205 Kokkos::parallel_for("point search conversion", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
206 (0,number_target_coords), [&](const int i) {
207 auto &point_cloud_neighbor_list_i = const_cast<std::map<int, double>& >(point_cloud_neighbor_list[i]);
208 for (int j=0; j<nla.getNumberOfNeighborsHost(i); ++j) {
209 point_cloud_neighbor_list_i[nla.getNeighborHost(i,j)] = 1.0;
210 }
211 });
212 double point_cloud_convert_time = timer.seconds();
213 printf("point cloud convert time: %f s\n", point_cloud_convert_time);
214 Kokkos::fence();
215
216 // loop over targets, finding all of their neighbors by brute force
217 timer.reset();
218 std::vector<std::map<int, double> > brute_force_neighbor_list(number_target_coords, std::map<int, double>());
219 Kokkos::parallel_for("brute force radius search", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
220 (0,number_target_coords), [&](const int i) {
221 auto &brute_force_neighbor_list_i = const_cast<std::map<int, double>& >(brute_force_neighbor_list[i]);
222 for (int j=0; j<number_source_coords; ++j) {
223 double dist = 0;
224 for (int k=0; k<dimension; ++k) {
225 dist += (target_coords(i,k)-source_coords(j,k))*(target_coords(i,k)-source_coords(j,k));
226 }
227 dist = std::sqrt(dist);
228
229 if (dist<(h_multiplier + random_vec(i))*h_spacing) {
230 brute_force_neighbor_list_i[j]=dist;
231 }
232 }
233 });
234 double brute_force_search_time = timer.seconds();
235 printf("brute force search time: %f s\n", brute_force_search_time);
236 Kokkos::fence();
237
238 timer.reset();
239 // check that all neighbors in brute force list are in point cloud search list
240 int num_passed = 0;
241 Kokkos::parallel_reduce("check brute force search in point cloud search", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
242 (0,number_target_coords), [&](const int i, int& t_num_passed) {
243 bool all_found = true;
244 for (auto it=brute_force_neighbor_list[i].begin(); it!=brute_force_neighbor_list[i].end(); ++it) {
245 if (point_cloud_neighbor_list[i].count(it->first)!=1) {
246 all_found = false;
247 }
248 }
249 if (all_found) t_num_passed++;
250 }, Kokkos::Sum<int>(num_passed));
251 Kokkos::fence();
252 if (num_passed != number_target_coords) {
253 printf("Brute force neighbor not found in point cloud list\n");
254 all_passed = false;
255 }
256
257 num_passed = 0;
258 // check that all neighbors in point cloud search list are in brute force list
259 Kokkos::parallel_reduce("original in brute force search", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
260 (0,number_target_coords), KOKKOS_LAMBDA(const int i, int& t_num_passed) {
261 bool all_found = true;
262 for (int j=0; j<nla.getNumberOfNeighborsHost(i); ++j) {
263 if (brute_force_neighbor_list[i].count(nla.getNeighborHost(i,j))!=1) {
264 all_found = false;
265 }
266 }
267 if (all_found) t_num_passed++;
268 }, Kokkos::Sum<int>(num_passed));
269 Kokkos::fence();
270 if (num_passed != number_target_coords) {
271 printf("Point cloud neighbor not found in brute force list\n");
272 all_passed = false;
273 }
274
275 double check_time = timer.seconds();
276 printf("check time: %f s\n", check_time);
277 }
278 timer.reset();
279 { // 2D neighbor lists
280 printf("\n2D neighbor lists:\n");
281 // Point cloud construction for neighbor search
282 // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
283 auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
284
285 Kokkos::View<int**, Kokkos::DefaultHostExecutionSpace> neighbor_lists("neighbor lists",
286 number_target_coords, 2); // first column is # of neighbors
287
288 // each target site has a window size
289 Kokkos::View<double*, Kokkos::DefaultHostExecutionSpace> epsilon("h supports", number_target_coords);
290 Kokkos::View<double*, Kokkos::DefaultHostExecutionSpace> random_vec("random vector", number_target_coords);
291
292 Kokkos::parallel_for("random perturbation of search size", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
293 (0,number_target_coords), [&](const int i) {
294 // value in [-.25, .25]
295 random_vec(i) = 0.5 * ((double)rand() / (double) RAND_MAX) - 0.5;
296 epsilon(i) = (h_multiplier + random_vec(i)) * h_spacing;
297 });
298 Kokkos::fence();
299
300 // query the point cloud to generate the neighbor lists using a radius search
301
302 // start with a dry-run, but without enough room to store the results
303 auto max_num_neighbors = point_cloud_search.generate2DNeighborListsFromRadiusSearch(true /* dry run */,
304 target_coords, neighbor_lists, epsilon);
305 printf("max num neighbors: %lu\n", max_num_neighbors);
306
307 // resize neighbor lists to be large enough to hold the results
308 neighbor_lists = Kokkos::View<int**, Kokkos::DefaultHostExecutionSpace>("neighbor lists",
309 number_target_coords, 1+max_num_neighbors); // first column is # of neighbors
310
311 // search again, now that we know that there is enough room to store the results
312 max_num_neighbors = point_cloud_search.generate2DNeighborListsFromRadiusSearch(false /* dry run */,
313 target_coords, neighbor_lists, epsilon);
314
315 double radius_search_time = timer.seconds();
316 printf("nanoflann search time: %f s\n", radius_search_time);
317
318 // convert point cloud search to vector of maps
319 timer.reset();
320 std::vector<std::map<int, double> > point_cloud_neighbor_list(number_target_coords, std::map<int, double>());
321 Kokkos::parallel_for("point search conversion", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
322 (0,number_target_coords), [&](const int i) {
323 auto &point_cloud_neighbor_list_i = const_cast<std::map<int, double>& >(point_cloud_neighbor_list[i]);
324 for (int j=1; j<=neighbor_lists(i,0); ++j) {
325 point_cloud_neighbor_list_i[neighbor_lists(i,j)] = 1.0;
326 }
327 });
328 double point_cloud_convert_time = timer.seconds();
329 printf("point cloud convert time: %f s\n", point_cloud_convert_time);
330 Kokkos::fence();
331
332 // loop over targets, finding all of their neighbors by brute force
333 timer.reset();
334 std::vector<std::map<int, double> > brute_force_neighbor_list(number_target_coords, std::map<int, double>());
335 Kokkos::parallel_for("brute force radius search", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
336 (0,number_target_coords), [&](const int i) {
337 auto &brute_force_neighbor_list_i = const_cast<std::map<int, double>& >(brute_force_neighbor_list[i]);
338 for (int j=0; j<number_source_coords; ++j) {
339 double dist = 0;
340 for (int k=0; k<dimension; ++k) {
341 dist += (target_coords(i,k)-source_coords(j,k))*(target_coords(i,k)-source_coords(j,k));
342 }
343 dist = std::sqrt(dist);
344
345 if (dist<(h_multiplier + random_vec(i))*h_spacing) {
346 brute_force_neighbor_list_i[j]=dist;
347 }
348 }
349 });
350 double brute_force_search_time = timer.seconds();
351 printf("brute force search time: %f s\n", brute_force_search_time);
352 Kokkos::fence();
353
354 timer.reset();
355 // check that all neighbors in brute force list are in point cloud search list
356 int num_passed = 0;
357 Kokkos::parallel_reduce("check brute force search in point cloud search", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
358 (0,number_target_coords), [&](const int i, int& t_num_passed) {
359 bool all_found = true;
360 for (auto it=brute_force_neighbor_list[i].begin(); it!=brute_force_neighbor_list[i].end(); ++it) {
361 if (point_cloud_neighbor_list[i].count(it->first)!=1) {
362 all_found = false;
363 }
364 }
365 if (all_found) t_num_passed++;
366 }, Kokkos::Sum<int>(num_passed));
367 Kokkos::fence();
368 if (num_passed != number_target_coords) {
369 printf("Brute force neighbor not found in point cloud list\n");
370 all_passed = false;
371 }
372
373 num_passed = 0;
374 // check that all neighbors in point cloud search list are in brute force list
375 Kokkos::parallel_reduce("original in brute force search", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
376 (0,number_target_coords), KOKKOS_LAMBDA(const int i, int& t_num_passed) {
377 bool all_found = true;
378 for (int j=1; j<=neighbor_lists(i,0); ++j) {
379 if (brute_force_neighbor_list[i].count(neighbor_lists(i,j))!=1) {
380 all_found = false;
381 }
382 }
383 if (all_found) t_num_passed++;
384 }, Kokkos::Sum<int>(num_passed));
385 Kokkos::fence();
386 if (num_passed != number_target_coords) {
387 printf("Point cloud neighbor not found in brute force list\n");
388 all_passed = false;
389 }
390
391 double check_time = timer.seconds();
392 printf("check time: %f s\n", check_time);
393 }
394 timer.reset();
395 { // convert 2D neighbor lists to 1D compressed row neighbor lists
396 printf("\n2D neighbor lists converted to 1D compressed row neighbor lists:\n");
397 // Point cloud construction for neighbor search
398 // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
399 auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
400
401 Kokkos::View<int**, Kokkos::DefaultHostExecutionSpace> neighbor_lists("neighbor lists",
402 number_target_coords, 2); // first column is # of neighbors
403
404 // each target site has a window size
405 Kokkos::View<double*, Kokkos::DefaultHostExecutionSpace> epsilon("h supports", number_target_coords);
406 Kokkos::View<double*, Kokkos::DefaultHostExecutionSpace> random_vec("random vector", number_target_coords);
407
408 Kokkos::parallel_for("random perturbation of search size", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
409 (0,number_target_coords), [&](const int i) {
410 // value in [-.25, .25]
411 random_vec(i) = 0.5 * ((double)rand() / (double) RAND_MAX) - 0.5;
412 epsilon(i) = (h_multiplier + random_vec(i)) * h_spacing;
413 });
414 Kokkos::fence();
415
416 // query the point cloud to generate the neighbor lists using a radius search
417
418 // start with a dry-run, but without enough room to store the results
419 auto max_num_neighbors = point_cloud_search.generate2DNeighborListsFromRadiusSearch(true /* dry run */,
420 target_coords, neighbor_lists, epsilon);
421 printf("max num neighbors: %lu\n", max_num_neighbors);
422
423 // resize neighbor lists to be large enough to hold the results
424 neighbor_lists = Kokkos::View<int**, Kokkos::DefaultHostExecutionSpace>("neighbor lists",
425 number_target_coords, 1+max_num_neighbors); // first column is # of neighbors
426
427 // search again, now that we know that there is enough room to store the results
428 max_num_neighbors = point_cloud_search.generate2DNeighborListsFromRadiusSearch(false /* dry run */,
429 target_coords, neighbor_lists, epsilon);
430
431 auto nla = Convert2DToCompressedRowNeighborLists(neighbor_lists);
432
433 double radius_search_time = timer.seconds();
434 printf("nanoflann search time: %f s\n", radius_search_time);
435
436 // convert point cloud search to vector of maps
437 timer.reset();
438 std::vector<std::map<int, double> > point_cloud_neighbor_list(number_target_coords, std::map<int, double>());
439 Kokkos::parallel_for("point search conversion", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
440 (0,number_target_coords), [&](const int i) {
441 auto &point_cloud_neighbor_list_i = const_cast<std::map<int, double>& >(point_cloud_neighbor_list[i]);
442 for (int j=0; j<nla.getNumberOfNeighborsHost(i); ++j) {
443 point_cloud_neighbor_list_i[nla.getNeighborHost(i,j)] = 1.0;
444 }
445 });
446 double point_cloud_convert_time = timer.seconds();
447 printf("point cloud convert time: %f s\n", point_cloud_convert_time);
448 Kokkos::fence();
449
450 // loop over targets, finding all of their neighbors by brute force
451 timer.reset();
452 std::vector<std::map<int, double> > brute_force_neighbor_list(number_target_coords, std::map<int, double>());
453 Kokkos::parallel_for("brute force radius search", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
454 (0,number_target_coords), [&](const int i) {
455 auto &brute_force_neighbor_list_i = const_cast<std::map<int, double>& >(brute_force_neighbor_list[i]);
456 for (int j=0; j<number_source_coords; ++j) {
457 double dist = 0;
458 for (int k=0; k<dimension; ++k) {
459 dist += (target_coords(i,k)-source_coords(j,k))*(target_coords(i,k)-source_coords(j,k));
460 }
461 dist = std::sqrt(dist);
462
463 if (dist<(h_multiplier + random_vec(i))*h_spacing) {
464 brute_force_neighbor_list_i[j]=dist;
465 }
466 }
467 });
468 double brute_force_search_time = timer.seconds();
469 printf("brute force search time: %f s\n", brute_force_search_time);
470 Kokkos::fence();
471
472 timer.reset();
473 // check that all neighbors in brute force list are in point cloud search list
474 int num_passed = 0;
475 Kokkos::parallel_reduce("check brute force search in point cloud search", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
476 (0,number_target_coords), [&](const int i, int& t_num_passed) {
477 bool all_found = true;
478 for (auto it=brute_force_neighbor_list[i].begin(); it!=brute_force_neighbor_list[i].end(); ++it) {
479 if (point_cloud_neighbor_list[i].count(it->first)!=1) {
480 all_found = false;
481 }
482 }
483 if (all_found) t_num_passed++;
484 }, Kokkos::Sum<int>(num_passed));
485 Kokkos::fence();
486 if (num_passed != number_target_coords) {
487 printf("Brute force neighbor not found in point cloud list\n");
488 all_passed = false;
489 }
490
491 num_passed = 0;
492 // check that all neighbors in point cloud search list are in brute force list
493 Kokkos::parallel_reduce("original in brute force search", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
494 (0,number_target_coords), KOKKOS_LAMBDA(const int i, int& t_num_passed) {
495 bool all_found = true;
496 for (int j=0; j<nla.getNumberOfNeighborsHost(i); ++j) {
497 if (brute_force_neighbor_list[i].count(nla.getNeighborHost(i,j))!=1) {
498 all_found = false;
499 }
500 }
501 if (all_found) t_num_passed++;
502 }, Kokkos::Sum<int>(num_passed));
503 Kokkos::fence();
504 if (num_passed != number_target_coords) {
505 printf("Point cloud neighbor not found in brute force list\n");
506 all_passed = false;
507 }
508
509 double check_time = timer.seconds();
510 printf("check time: %f s\n", check_time);
511 }
512 }
513 else if (do_knn_search) {
514
515 // do knn search
516 timer.reset();
517 { // 1D compressed row neighbor lists
518 printf("\n1D compressed row neighbor lists:\n");
519 // Point cloud construction for neighbor search
520 // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
521 auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
522
523 const int min_neighbors = Compadre::GMLS::getNP(3, dimension);
524 Kokkos::View<int*, Kokkos::DefaultHostExecutionSpace> neighbor_lists("neighbor lists",
525 number_target_coords); // first column is # of neighbors
526 Kokkos::View<int*, Kokkos::DefaultHostExecutionSpace> number_neighbors_list("number of neighbors list",
527 number_target_coords); // first column is # of neighbors
528
529 // each target site has a window size
530 Kokkos::View<double*, Kokkos::DefaultHostExecutionSpace> epsilon("h supports", number_target_coords);
531
532 // query the point cloud to generate the neighbor lists using a KNN search
533
534 // start with a dry-run, but without enough room to store the results
535 auto total_num_neighbors = point_cloud_search.generateCRNeighborListsFromKNNSearch(true /* dry-run for sizes */,
536 target_coords, neighbor_lists, number_neighbors_list, epsilon, min_neighbors, 1.5 /* cutoff_multiplier */);
537 printf("total num neighbors: %lu\n", total_num_neighbors);
538
539 // resize with room to store results
540 Kokkos::resize(neighbor_lists, total_num_neighbors);
541
542 // real knn search with space to store
543 point_cloud_search.generateCRNeighborListsFromKNNSearch(false /*not dry run*/,
544 target_coords, neighbor_lists, number_neighbors_list, epsilon, min_neighbors, 1.5 /* cutoff_multiplier */);
545
546 auto nla(CreateNeighborLists(neighbor_lists, number_neighbors_list));
547
548 // convert point cloud search to vector of maps
549 timer.reset();
550 std::vector<std::map<int, double> > point_cloud_neighbor_list(number_target_coords, std::map<int, double>());
551 Kokkos::parallel_for("point search conversion", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
552 (0,number_target_coords), [&](const int i) {
553 auto &point_cloud_neighbor_list_i = const_cast<std::map<int, double>& >(point_cloud_neighbor_list[i]);
554 for (int j=0; j<nla.getNumberOfNeighborsHost(i); ++j) {
555 double dist = 0;
556 for (int k=0; k<dimension; ++k) {
557 dist += (target_coords(i,k)-source_coords(nla.getNeighborHost(i,j),k))*(target_coords(i,k)-source_coords(nla.getNeighborHost(i,j),k));
558 }
559 dist = std::sqrt(dist);
560 point_cloud_neighbor_list_i[nla.getNeighborHost(i,j)] = dist;
561 }
562 });
563 double point_cloud_convert_time = timer.seconds();
564 printf("point cloud convert time: %f s\n", point_cloud_convert_time);
565 Kokkos::fence();
566
567 double eps = 1e-14;
568 // need to do a sort by dist, then verify epsilon(i) = 1.5 * knn_distance
569 int epsilon_diff_from_knn_dist_multiplied = 0;
570 Kokkos::parallel_reduce("check kth neighbor in point cloud", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
571 (0,number_target_coords), [&](const int i, int& t_epsilon_diff_from_knn_dist_multiplied) {
572 auto &point_cloud_neighbor_list_i = const_cast<std::map<int, double>& >(point_cloud_neighbor_list[i]);
573 std::multimap<double, int> inverted_map = invert_map(point_cloud_neighbor_list_i);
574 int j = 0;
575 for (auto it=inverted_map.begin(); it!=inverted_map.end(); ++it) {
576 if (j==min_neighbors-1) {
577 if (std::abs(epsilon(i)/1.5 - it->first) > eps) {
578 t_epsilon_diff_from_knn_dist_multiplied++;
579 }
580 break;
581 }
582 j++;
583 }
584 }, Kokkos::Sum<int>(epsilon_diff_from_knn_dist_multiplied));
585 if (epsilon_diff_from_knn_dist_multiplied > 0) {
586 printf("1.5*kth_neighbor_distance differs from calculation in epsilon vector.\n");
587 all_passed = false;
588 }
589
590 // loop over targets, finding knn by brute force
591 // point cloud multiplied kth by 1.5, but we are searching by brute force to ensure that for the distance h,
592 // we find at least as many neighbors as min_neighbors
593 // (not for the 1.5x as many, but rather for 1.0x as many)
594 timer.reset();
595 int total_insufficient_neighbors = 0;
596 std::vector<std::map<int, double> > brute_force_neighbor_list(number_target_coords, std::map<int, double>());
597 Kokkos::parallel_reduce("brute force knn search", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
598 (0,number_target_coords), [&](const int i, int& t_total_insufficient_neighbors) {
599 auto &brute_force_neighbor_list_i = const_cast<std::map<int, double>& >(brute_force_neighbor_list[i]);
600 for (int j=0; j<number_source_coords; ++j) {
601 double dist = 0;
602 for (int k=0; k<dimension; ++k) {
603 dist += (target_coords(i,k)-source_coords(j,k))*(target_coords(i,k)-source_coords(j,k));
604 }
605 dist = std::sqrt(dist);
606
607 // epsilon was set to 1.5 * distance to kth neighbor, but we need distance to kth neighbor
608 if (dist<(epsilon(i)/1.5 + eps)) {
609 brute_force_neighbor_list_i[j]=dist;
610 }
611 }
612 // verify k neighbors found
613 if (brute_force_neighbor_list_i.size() < (size_t)min_neighbors) t_total_insufficient_neighbors++;
614 }, Kokkos::Sum<int>(total_insufficient_neighbors));
615 if (total_insufficient_neighbors > 0) {
616 printf("less neighbors found using kth_neighbor_distance+eps than using knn search.\n");
617 all_passed = false;
618 }
619 double brute_force_search_time = timer.seconds();
620 printf("brute force search time: %f s\n", brute_force_search_time);
621 Kokkos::fence();
622 }
623 timer.reset();
624 { // 2D neighbor lists
625 printf("\n2D neighbor lists:\n");
626 // Point cloud construction for neighbor search
627 // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
628 auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
629
630 const int min_neighbors = Compadre::GMLS::getNP(3, dimension);
631 Kokkos::View<int**, Kokkos::DefaultHostExecutionSpace> neighbor_lists("neighbor lists",
632 number_target_coords, 1+min_neighbors); // first column is # of neighbors
633
634 // each target site has a window size
635 Kokkos::View<double*, Kokkos::DefaultHostExecutionSpace> epsilon("h supports", number_target_coords);
636
637 // query the point cloud to generate the neighbor lists using a KNN search
638 // start with a dry-run, but without enough room to store the results
639 auto max_num_neighbors = point_cloud_search.generate2DNeighborListsFromKNNSearch(true /* dry-run for sizes */,
640 target_coords, neighbor_lists, epsilon, min_neighbors, 1.5 /* cutoff_multiplier */);
641 printf("max num neighbors: %lu\n", max_num_neighbors);
642
643 // resize with room to store results
644 Kokkos::resize(neighbor_lists, neighbor_lists.extent(0), 1+max_num_neighbors);
645
646 // real knn search with space to store
647 max_num_neighbors = point_cloud_search.generate2DNeighborListsFromKNNSearch(false /*not dry run*/,
648 target_coords, neighbor_lists, epsilon, min_neighbors, 1.5 /* cutoff_multiplier */);
649
650 // convert point cloud search to vector of maps
651 timer.reset();
652 std::vector<std::map<int, double> > point_cloud_neighbor_list(number_target_coords, std::map<int, double>());
653 Kokkos::parallel_for("point search conversion", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
654 (0,number_target_coords), [&](const int i) {
655 auto &point_cloud_neighbor_list_i = const_cast<std::map<int, double>& >(point_cloud_neighbor_list[i]);
656 for (int j=1; j<=neighbor_lists(i,0); ++j) {
657 double dist = 0;
658 for (int k=0; k<dimension; ++k) {
659 dist += (target_coords(i,k)-source_coords(neighbor_lists(i,j),k))*(target_coords(i,k)-source_coords(neighbor_lists(i,j),k));
660 }
661 dist = std::sqrt(dist);
662 point_cloud_neighbor_list_i[neighbor_lists(i,j)] = dist;
663 }
664 });
665 double point_cloud_convert_time = timer.seconds();
666 printf("point cloud convert time: %f s\n", point_cloud_convert_time);
667 Kokkos::fence();
668
669 double eps = 1e-14;
670 // need to do a sort by dist, then verify epsilon(i) = 1.5 * knn_distance
671 int epsilon_diff_from_knn_dist_multiplied = 0;
672 Kokkos::parallel_reduce("check kth neighbor in point cloud", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
673 (0,number_target_coords), [&](const int i, int& t_epsilon_diff_from_knn_dist_multiplied) {
674 auto &point_cloud_neighbor_list_i = const_cast<std::map<int, double>& >(point_cloud_neighbor_list[i]);
675 std::multimap<double, int> inverted_map = invert_map(point_cloud_neighbor_list_i);
676 int j = 0;
677 for (auto it=inverted_map.begin(); it!=inverted_map.end(); ++it) {
678 if (j==min_neighbors-1) {
679 if (std::abs(epsilon(i)/1.5 - it->first) > eps) {
680 t_epsilon_diff_from_knn_dist_multiplied++;
681 }
682 break;
683 }
684 j++;
685 }
686 }, Kokkos::Sum<int>(epsilon_diff_from_knn_dist_multiplied));
687 if (epsilon_diff_from_knn_dist_multiplied > 0) {
688 printf("1.5*kth_neighbor_distance differs from calculation in epsilon vector.\n");
689 all_passed = false;
690 }
691
692 // loop over targets, finding knn by brute force
693 // point cloud multiplied kth by 1.5, but we are searching by brute force to ensure that for the distance h,
694 // we find at least as many neighbors as min_neighbors
695 // (not for the 1.5x as many, but rather for 1.0x as many)
696 timer.reset();
697 int total_insufficient_neighbors = 0;
698 std::vector<std::map<int, double> > brute_force_neighbor_list(number_target_coords, std::map<int, double>());
699 Kokkos::parallel_reduce("brute force knn search", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
700 (0,number_target_coords), [&](const int i, int& t_total_insufficient_neighbors) {
701 auto &brute_force_neighbor_list_i = const_cast<std::map<int, double>& >(brute_force_neighbor_list[i]);
702 for (int j=0; j<number_source_coords; ++j) {
703 double dist = 0;
704 for (int k=0; k<dimension; ++k) {
705 dist += (target_coords(i,k)-source_coords(j,k))*(target_coords(i,k)-source_coords(j,k));
706 }
707 dist = std::sqrt(dist);
708
709 // epsilon was set to 1.5 * distance to kth neighbor, but we need distance to kth neighbor
710 if (dist<(epsilon(i)/1.5 + eps)) {
711 brute_force_neighbor_list_i[j]=dist;
712 }
713 }
714 // verify k neighbors found
715 if (brute_force_neighbor_list_i.size() < (size_t)min_neighbors) t_total_insufficient_neighbors++;
716 }, Kokkos::Sum<int>(total_insufficient_neighbors));
717 if (total_insufficient_neighbors > 0) {
718 printf("less neighbors found using kth_neighbor_distance+eps than using knn search.\n");
719 all_passed = false;
720 }
721 double brute_force_search_time = timer.seconds();
722 printf("brute force search time: %f s\n", brute_force_search_time);
723 Kokkos::fence();
724 }
725 timer.reset();
726 { // convert 2D neighbor lists to 1D compressed row neighbor lists
727 printf("\n2D neighbor lists converted to 1D compressed row neighbor lists:\n");
728 // Point cloud construction for neighbor search
729 // CreatePointCloudSearch constructs an object of type PointCloudSearch, but deduces the templates for you
730 auto point_cloud_search(CreatePointCloudSearch(source_coords, dimension));
731
732 const int min_neighbors = Compadre::GMLS::getNP(3, dimension);
733 Kokkos::View<int**, Kokkos::DefaultHostExecutionSpace> neighbor_lists("neighbor lists",
734 number_target_coords, 1+min_neighbors); // first column is # of neighbors
735
736 // each target site has a window size
737 Kokkos::View<double*, Kokkos::DefaultHostExecutionSpace> epsilon("h supports", number_target_coords);
738
739 // query the point cloud to generate the neighbor lists using a KNN search
740 // start with a dry-run, but without enough room to store the results
741 auto max_num_neighbors = point_cloud_search.generate2DNeighborListsFromKNNSearch(true /* dry-run for sizes */,
742 target_coords, neighbor_lists, epsilon, min_neighbors, 1.5 /* cutoff_multiplier */);
743 printf("max num neighbors: %lu\n", max_num_neighbors);
744
745 // resize with room to store results
746 Kokkos::resize(neighbor_lists, neighbor_lists.extent(0), 1+max_num_neighbors);
747
748 // real knn search with space to store
749 max_num_neighbors = point_cloud_search.generate2DNeighborListsFromKNNSearch(false /*not dry run*/,
750 target_coords, neighbor_lists, epsilon, min_neighbors, 1.5 /* cutoff_multiplier */);
751
752 auto nla = Convert2DToCompressedRowNeighborLists(neighbor_lists);
753
754 // convert point cloud search to vector of maps
755 timer.reset();
756 std::vector<std::map<int, double> > point_cloud_neighbor_list(number_target_coords, std::map<int, double>());
757 Kokkos::parallel_for("point search conversion", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
758 (0,number_target_coords), [&](const int i) {
759 auto &point_cloud_neighbor_list_i = const_cast<std::map<int, double>& >(point_cloud_neighbor_list[i]);
760 for (int j=0; j<nla.getNumberOfNeighborsHost(i); ++j) {
761 double dist = 0;
762 for (int k=0; k<dimension; ++k) {
763 dist += (target_coords(i,k)-source_coords(nla.getNeighborHost(i,j),k))*(target_coords(i,k)-source_coords(nla.getNeighborHost(i,j),k));
764 }
765 dist = std::sqrt(dist);
766 point_cloud_neighbor_list_i[nla.getNeighborHost(i,j)] = dist;
767 }
768 });
769 double point_cloud_convert_time = timer.seconds();
770 printf("point cloud convert time: %f s\n", point_cloud_convert_time);
771 Kokkos::fence();
772
773 double eps = 1e-14;
774 // need to do a sort by dist, then verify epsilon(i) = 1.5 * knn_distance
775 int epsilon_diff_from_knn_dist_multiplied = 0;
776 Kokkos::parallel_reduce("check kth neighbor in point cloud", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
777 (0,number_target_coords), [&](const int i, int& t_epsilon_diff_from_knn_dist_multiplied) {
778 auto &point_cloud_neighbor_list_i = const_cast<std::map<int, double>& >(point_cloud_neighbor_list[i]);
779 std::multimap<double, int> inverted_map = invert_map(point_cloud_neighbor_list_i);
780 int j = 0;
781 for (auto it=inverted_map.begin(); it!=inverted_map.end(); ++it) {
782 if (j==min_neighbors-1) {
783 if (std::abs(epsilon(i)/1.5 - it->first) > eps) {
784 t_epsilon_diff_from_knn_dist_multiplied++;
785 }
786 break;
787 }
788 j++;
789 }
790 }, Kokkos::Sum<int>(epsilon_diff_from_knn_dist_multiplied));
791 if (epsilon_diff_from_knn_dist_multiplied > 0) {
792 printf("1.5*kth_neighbor_distance differs from calculation in epsilon vector.\n");
793 all_passed = false;
794 }
795
796 // loop over targets, finding knn by brute force
797 // point cloud multiplied kth by 1.5, but we are searching by brute force to ensure that for the distance h,
798 // we find at least as many neighbors as min_neighbors
799 // (not for the 1.5x as many, but rather for 1.0x as many)
800 timer.reset();
801 int total_insufficient_neighbors = 0;
802 std::vector<std::map<int, double> > brute_force_neighbor_list(number_target_coords, std::map<int, double>());
803 Kokkos::parallel_reduce("brute force knn search", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
804 (0,number_target_coords), [&](const int i, int& t_total_insufficient_neighbors) {
805 auto &brute_force_neighbor_list_i = const_cast<std::map<int, double>& >(brute_force_neighbor_list[i]);
806 for (int j=0; j<number_source_coords; ++j) {
807 double dist = 0;
808 for (int k=0; k<dimension; ++k) {
809 dist += (target_coords(i,k)-source_coords(j,k))*(target_coords(i,k)-source_coords(j,k));
810 }
811 dist = std::sqrt(dist);
812
813 // epsilon was set to 1.5 * distance to kth neighbor, but we need distance to kth neighbor
814 if (dist<(epsilon(i)/1.5 + eps)) {
815 brute_force_neighbor_list_i[j]=dist;
816 }
817 }
818 // verify k neighbors found
819 if (brute_force_neighbor_list_i.size() < (size_t)min_neighbors) t_total_insufficient_neighbors++;
820 }, Kokkos::Sum<int>(total_insufficient_neighbors));
821 if (total_insufficient_neighbors > 0) {
822 printf("less neighbors found using kth_neighbor_distance+eps than using knn search.\n");
823 all_passed = false;
824 }
825 double brute_force_search_time = timer.seconds();
826 printf("brute force search time: %f s\n", brute_force_search_time);
827 Kokkos::fence();
828 }
829 } else {
830 all_passed = false;
831 printf("do_radius_search and do_knn_search both false. Invalid combination.\n");
832 }
833}
834
835// finalize Kokkos and MPI (if available)
836Kokkos::finalize();
837#ifdef COMPADRE_USE_MPI
838MPI_Finalize();
839#endif
840
841// output to user that test passed or failed
842if(all_passed) {
843 fprintf(stdout, "Passed test \n");
844 return 0;
845} else {
846 fprintf(stdout, "Failed test \n");
847 return -1;
848}
849
850} // main
std::pair< T2, T1 > swap_pair(const std::pair< T1, T2 > &item)
std::multimap< T2, T1 > invert_map(const std::map< T1, T2 > &original_map)
int main(int argc, char *args[])
Manifold GMLS Example.
static KOKKOS_INLINE_FUNCTION int getNP(const int m, const int dimension=3, const ReconstructionSpace r_space=ReconstructionSpace::ScalarTaylorPolynomial)
Returns size of the basis for a given polynomial order and dimension General to dimension 1....
PointCloudSearch< view_type > CreatePointCloudSearch(view_type src_view, const local_index_type dimensions=-1, const local_index_type max_leaf=-1)
CreatePointCloudSearch allows for the construction of an object of type PointCloudSearch with templat...
NeighborLists< view_type > CreateNeighborLists(view_type number_of_neighbors_list)
CreateNeighborLists allows for the construction of an object of type NeighborLists with template dedu...
NeighborLists< view_type_1d > Convert2DToCompressedRowNeighborLists(view_type_2d neighbor_lists)
Converts 2D neighbor lists to compressed row neighbor lists.