5 #include <Compadre_Config.h>
9 #ifdef COMPADRE_USE_MPI
13 template<
typename T1,
typename T2>
14 std::pair<T2,T1>
swap_pair(
const std::pair<T1,T2> &item)
16 return std::pair<T2,T1>(item.second, item.first);
19 template<
typename T1,
typename T2>
20 std::multimap<T2,T1>
invert_map(
const std::map<T1,T2> &original_map)
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>);
27 using namespace Compadre;
29 int main (
int argc,
char* args[]) {
32 #ifdef COMPADRE_USE_MPI
33 MPI_Init(&argc, &args);
37 Kokkos::initialize(argc, args);
40 bool all_passed =
true;
49 int arg5toi = atoi(args[4]);
51 search_type = arg5toi;
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);
61 double h_multiplier = 3.0;
63 double arg4tof = atof(args[3]);
65 h_multiplier = arg4tof;
71 int number_target_coords = 200;
73 int arg3toi = atoi(args[2]);
75 number_target_coords = arg3toi;
83 int arg2toi = atoi(args[1]);
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;
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);
103 Kokkos::View<double**, Kokkos::DefaultHostExecutionSpace> source_coords(
"source coordinates",
104 number_source_coords, 3);
107 Kokkos::View<double**, Kokkos::DefaultHostExecutionSpace> target_coords(
"target coordinates", number_target_coords, 3);
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;
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];
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;
133 source_coords(source_index,0) = this_coord[0];
134 source_coords(source_index,1) = 0;
135 source_coords(source_index,2) = 0;
141 for(
int i=0; i<number_target_coords; i++){
144 double rand_dir[3] = {0,0,0};
146 for (
int j=0; j<dimension; ++j) {
148 rand_dir[j] = ((double)rand() / (double) RAND_MAX) - 0.5;
152 for (
int j=0; j<dimension; ++j) {
153 target_coords(i,j) = rand_dir[j];
158 if (do_radius_search) {
161 printf(
"\n1D compressed row neighbor lists:\n");
166 Kokkos::View<int*, Kokkos::DefaultHostExecutionSpace> neighbor_lists(
"neighbor lists",
167 number_target_coords);
168 Kokkos::View<int*, Kokkos::DefaultHostExecutionSpace> number_neighbors_list(
"number of neighbors list",
169 number_target_coords);
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);
175 Kokkos::parallel_for(
"random perturbation of search size", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
176 (0,number_target_coords), [&](
const int i) {
178 random_vec(i) = 0.5 * ((double)rand() / (double) RAND_MAX) - 0.5;
179 epsilon(i) = (h_multiplier + random_vec(i)) * h_spacing;
186 size_t total_num_neighbors = point_cloud_search.generateCRNeighborListsFromRadiusSearch(
true ,
187 target_coords, neighbor_lists, number_neighbors_list, epsilon);
188 printf(
"total num neighbors: %lu\n", total_num_neighbors);
191 Kokkos::resize(neighbor_lists, total_num_neighbors);
194 point_cloud_search.generateCRNeighborListsFromRadiusSearch(
false ,
195 target_coords, neighbor_lists, number_neighbors_list, epsilon);
199 double radius_search_time = timer.seconds();
200 printf(
"nanoflann search time: %f s\n", radius_search_time);
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;
212 double point_cloud_convert_time = timer.seconds();
213 printf(
"point cloud convert time: %f s\n", point_cloud_convert_time);
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) {
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));
227 dist = std::sqrt(dist);
229 if (dist<(h_multiplier + random_vec(i))*h_spacing) {
230 brute_force_neighbor_list_i[j]=dist;
234 double brute_force_search_time = timer.seconds();
235 printf(
"brute force search time: %f s\n", brute_force_search_time);
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) {
249 if (all_found) t_num_passed++;
250 }, Kokkos::Sum<int>(num_passed));
252 if (num_passed != number_target_coords) {
253 printf(
"Brute force neighbor not found in point cloud list\n");
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) {
267 if (all_found) t_num_passed++;
268 }, Kokkos::Sum<int>(num_passed));
270 if (num_passed != number_target_coords) {
271 printf(
"Point cloud neighbor not found in brute force list\n");
275 double check_time = timer.seconds();
276 printf(
"check time: %f s\n", check_time);
280 printf(
"\n2D neighbor lists:\n");
285 Kokkos::View<int**, Kokkos::DefaultHostExecutionSpace> neighbor_lists(
"neighbor lists",
286 number_target_coords, 2);
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);
292 Kokkos::parallel_for(
"random perturbation of search size", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
293 (0,number_target_coords), [&](
const int i) {
295 random_vec(i) = 0.5 * ((double)rand() / (double) RAND_MAX) - 0.5;
296 epsilon(i) = (h_multiplier + random_vec(i)) * h_spacing;
303 auto max_num_neighbors = point_cloud_search.generate2DNeighborListsFromRadiusSearch(
true ,
304 target_coords, neighbor_lists, epsilon);
305 printf(
"max num neighbors: %lu\n", max_num_neighbors);
308 neighbor_lists = Kokkos::View<int**, Kokkos::DefaultHostExecutionSpace>(
"neighbor lists",
309 number_target_coords, 1+max_num_neighbors);
312 max_num_neighbors = point_cloud_search.generate2DNeighborListsFromRadiusSearch(
false ,
313 target_coords, neighbor_lists, epsilon);
315 double radius_search_time = timer.seconds();
316 printf(
"nanoflann search time: %f s\n", radius_search_time);
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;
328 double point_cloud_convert_time = timer.seconds();
329 printf(
"point cloud convert time: %f s\n", point_cloud_convert_time);
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) {
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));
343 dist = std::sqrt(dist);
345 if (dist<(h_multiplier + random_vec(i))*h_spacing) {
346 brute_force_neighbor_list_i[j]=dist;
350 double brute_force_search_time = timer.seconds();
351 printf(
"brute force search time: %f s\n", brute_force_search_time);
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) {
365 if (all_found) t_num_passed++;
366 }, Kokkos::Sum<int>(num_passed));
368 if (num_passed != number_target_coords) {
369 printf(
"Brute force neighbor not found in point cloud list\n");
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) {
383 if (all_found) t_num_passed++;
384 }, Kokkos::Sum<int>(num_passed));
386 if (num_passed != number_target_coords) {
387 printf(
"Point cloud neighbor not found in brute force list\n");
391 double check_time = timer.seconds();
392 printf(
"check time: %f s\n", check_time);
396 printf(
"\n2D neighbor lists converted to 1D compressed row neighbor lists:\n");
401 Kokkos::View<int**, Kokkos::DefaultHostExecutionSpace> neighbor_lists(
"neighbor lists",
402 number_target_coords, 2);
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);
408 Kokkos::parallel_for(
"random perturbation of search size", Kokkos::RangePolicy<Kokkos::DefaultHostExecutionSpace>
409 (0,number_target_coords), [&](
const int i) {
411 random_vec(i) = 0.5 * ((double)rand() / (double) RAND_MAX) - 0.5;
412 epsilon(i) = (h_multiplier + random_vec(i)) * h_spacing;
419 auto max_num_neighbors = point_cloud_search.generate2DNeighborListsFromRadiusSearch(
true ,
420 target_coords, neighbor_lists, epsilon);
421 printf(
"max num neighbors: %lu\n", max_num_neighbors);
424 neighbor_lists = Kokkos::View<int**, Kokkos::DefaultHostExecutionSpace>(
"neighbor lists",
425 number_target_coords, 1+max_num_neighbors);
428 max_num_neighbors = point_cloud_search.generate2DNeighborListsFromRadiusSearch(
false ,
429 target_coords, neighbor_lists, epsilon);
433 double radius_search_time = timer.seconds();
434 printf(
"nanoflann search time: %f s\n", radius_search_time);
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;
446 double point_cloud_convert_time = timer.seconds();
447 printf(
"point cloud convert time: %f s\n", point_cloud_convert_time);
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) {
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));
461 dist = std::sqrt(dist);
463 if (dist<(h_multiplier + random_vec(i))*h_spacing) {
464 brute_force_neighbor_list_i[j]=dist;
468 double brute_force_search_time = timer.seconds();
469 printf(
"brute force search time: %f s\n", brute_force_search_time);
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) {
483 if (all_found) t_num_passed++;
484 }, Kokkos::Sum<int>(num_passed));
486 if (num_passed != number_target_coords) {
487 printf(
"Brute force neighbor not found in point cloud list\n");
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) {
501 if (all_found) t_num_passed++;
502 }, Kokkos::Sum<int>(num_passed));
504 if (num_passed != number_target_coords) {
505 printf(
"Point cloud neighbor not found in brute force list\n");
509 double check_time = timer.seconds();
510 printf(
"check time: %f s\n", check_time);
513 else if (do_knn_search) {
518 printf(
"\n1D compressed row neighbor lists:\n");
524 Kokkos::View<int*, Kokkos::DefaultHostExecutionSpace> neighbor_lists(
"neighbor lists",
525 number_target_coords);
526 Kokkos::View<int*, Kokkos::DefaultHostExecutionSpace> number_neighbors_list(
"number of neighbors list",
527 number_target_coords);
530 Kokkos::View<double*, Kokkos::DefaultHostExecutionSpace> epsilon(
"h supports", number_target_coords);
535 auto total_num_neighbors = point_cloud_search.generateCRNeighborListsFromKNNSearch(
true ,
536 target_coords, neighbor_lists, number_neighbors_list, epsilon, min_neighbors, 1.5 );
537 printf(
"total num neighbors: %lu\n", total_num_neighbors);
540 Kokkos::resize(neighbor_lists, total_num_neighbors);
543 point_cloud_search.generateCRNeighborListsFromKNNSearch(
false ,
544 target_coords, neighbor_lists, number_neighbors_list, epsilon, min_neighbors, 1.5 );
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) {
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));
559 dist = std::sqrt(dist);
560 point_cloud_neighbor_list_i[nla.getNeighborHost(i,j)] = dist;
563 double point_cloud_convert_time = timer.seconds();
564 printf(
"point cloud convert time: %f s\n", point_cloud_convert_time);
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);
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++;
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");
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) {
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));
605 dist = std::sqrt(dist);
608 if (dist<(epsilon(i)/1.5 + eps)) {
609 brute_force_neighbor_list_i[j]=dist;
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");
619 double brute_force_search_time = timer.seconds();
620 printf(
"brute force search time: %f s\n", brute_force_search_time);
625 printf(
"\n2D neighbor lists:\n");
631 Kokkos::View<int**, Kokkos::DefaultHostExecutionSpace> neighbor_lists(
"neighbor lists",
632 number_target_coords, 1+min_neighbors);
635 Kokkos::View<double*, Kokkos::DefaultHostExecutionSpace> epsilon(
"h supports", number_target_coords);
639 auto max_num_neighbors = point_cloud_search.generate2DNeighborListsFromKNNSearch(
true ,
640 target_coords, neighbor_lists, epsilon, min_neighbors, 1.5 );
641 printf(
"max num neighbors: %lu\n", max_num_neighbors);
644 Kokkos::resize(neighbor_lists, neighbor_lists.extent(0), 1+max_num_neighbors);
647 max_num_neighbors = point_cloud_search.generate2DNeighborListsFromKNNSearch(
false ,
648 target_coords, neighbor_lists, epsilon, min_neighbors, 1.5 );
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) {
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));
661 dist = std::sqrt(dist);
662 point_cloud_neighbor_list_i[neighbor_lists(i,j)] = dist;
665 double point_cloud_convert_time = timer.seconds();
666 printf(
"point cloud convert time: %f s\n", point_cloud_convert_time);
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);
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++;
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");
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) {
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));
707 dist = std::sqrt(dist);
710 if (dist<(epsilon(i)/1.5 + eps)) {
711 brute_force_neighbor_list_i[j]=dist;
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");
721 double brute_force_search_time = timer.seconds();
722 printf(
"brute force search time: %f s\n", brute_force_search_time);
727 printf(
"\n2D neighbor lists converted to 1D compressed row neighbor lists:\n");
733 Kokkos::View<int**, Kokkos::DefaultHostExecutionSpace> neighbor_lists(
"neighbor lists",
734 number_target_coords, 1+min_neighbors);
737 Kokkos::View<double*, Kokkos::DefaultHostExecutionSpace> epsilon(
"h supports", number_target_coords);
741 auto max_num_neighbors = point_cloud_search.generate2DNeighborListsFromKNNSearch(
true ,
742 target_coords, neighbor_lists, epsilon, min_neighbors, 1.5 );
743 printf(
"max num neighbors: %lu\n", max_num_neighbors);
746 Kokkos::resize(neighbor_lists, neighbor_lists.extent(0), 1+max_num_neighbors);
749 max_num_neighbors = point_cloud_search.generate2DNeighborListsFromKNNSearch(
false ,
750 target_coords, neighbor_lists, epsilon, min_neighbors, 1.5 );
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) {
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));
765 dist = std::sqrt(dist);
766 point_cloud_neighbor_list_i[nla.getNeighborHost(i,j)] = dist;
769 double point_cloud_convert_time = timer.seconds();
770 printf(
"point cloud convert time: %f s\n", point_cloud_convert_time);
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);
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++;
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");
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) {
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));
811 dist = std::sqrt(dist);
814 if (dist<(epsilon(i)/1.5 + eps)) {
815 brute_force_neighbor_list_i[j]=dist;
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");
825 double brute_force_search_time = timer.seconds();
826 printf(
"brute force search time: %f s\n", brute_force_search_time);
831 printf(
"do_radius_search and do_knn_search both false. Invalid combination.\n");
837 #ifdef COMPADRE_USE_MPI
843 fprintf(stdout,
"Passed test \n");
846 fprintf(stdout,
"Failed test \n");
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...
int main(int argc, char **argv)
NeighborLists< view_type > CreateNeighborLists(view_type number_of_neighbors_list)
CreateNeighborLists allows for the construction of an object of type NeighborLists with template dedu...
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)
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...
NeighborLists< view_type_1d > Convert2DToCompressedRowNeighborLists(view_type_2d neighbor_lists)
Converts 2D neighbor lists to compressed row neighbor lists.