From 81ab45d2281a3326a02614db210878b81274c880 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 5 Mar 2013 22:32:11 -0500 Subject: [PATCH] Added new test for periodic KDTree --- sample/CMakeLists.txt | 3 + sample/testkd3.cpp | 186 +++++++++++++++++++++++++++++++++++++ src/replicateGenerator.hpp | 2 + 3 files changed, 191 insertions(+) create mode 100644 sample/testkd3.cpp diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt index 765e1b3..d158ae2 100644 --- a/sample/CMakeLists.txt +++ b/sample/CMakeLists.txt @@ -20,6 +20,9 @@ target_link_libraries(testkd ${tolink}) add_executable(testkd2 testkd2.cpp) target_link_libraries(testkd2 ${tolink}) +add_executable(testkd3 testkd3.cpp) +target_link_libraries(testkd3 ${tolink}) + add_executable(testDelaunay testDelaunay.cpp) target_link_libraries(testDelaunay ${tolink}) diff --git a/sample/testkd3.cpp b/sample/testkd3.cpp new file mode 100644 index 0000000..5744985 --- /dev/null +++ b/sample/testkd3.cpp @@ -0,0 +1,186 @@ +/*+ +This is CosmoTool (./sample/testkd2.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ +#include +#include +#include +#include +#include +#define __KD_TREE_NUMNODES +#include "mykdtree.hpp" +#include "kdtree_splitters.hpp" +#include + +#define NTRY 30 +#define NGB 24 +#define ND 3 + +using namespace std; +using namespace CosmoTool; + +typedef KDTree > MyTree; +//typedef KDTree MyTree; +typedef KDCell MyCell; + +static double periodic(double a) +{ + while (a < -0.5) + a+=1; + while (a > 0.5) + a -= 1; + return a; +} + +MyCell *findNearest(MyTree::coords& xc, MyCell *cells, uint32_t Ncells) +{ + MyCell *near2 = 0; + double R2 = INFINITY; + for (int i = 0; i < Ncells; i++) + { + double d2 = 0; + for (int j = 0; j < ND; j++) + { + double delta = periodic(xc[j]-cells[i].coord[j]); + d2 += delta*delta; + } + if (d2 < R2) + { + near2 = &cells[i]; + R2 = d2; + } + } + return near2; +} + +int main() +{ + uint32_t Ncells = 3000000; + MyCell *cells = new MyCell[Ncells]; + + for (int i = 0; i < Ncells; i++) + { + cells[i].active = true; + for (int l = 0; l < ND; l++) + cells[i].coord[l] = drand48(); + } + + // Check timing + clock_t startTimer = clock(); + MyTree tree(cells, Ncells); + clock_t endTimer = clock(); + + tree.setPeriodic(true, 1.0); + + clock_t delta = endTimer-startTimer; + double myTime = delta*1.0/CLOCKS_PER_SEC * 1.0; + + cout << "KDTree build = " << myTime << " s" << endl; + + MyTree::coords *xc = new MyTree::coords[NTRY]; + + cout << "Generating seeds..." << endl; + for (int k = 0; k < NTRY; k++) + { + for (int l = 0; l < ND; l++) + xc[k][l] = drand48(); + } + + // Check consistency + cout << "Check consistency..." << endl; + MyCell **ngb = new MyCell *[NGB]; + double *distances = new double[NGB]; + + ofstream fngb("nearest.txt"); + for (int k = 0; k < NTRY; k++) { + cout << "Seed = " << xc[k][0] << " " << xc[k][1] << " " << xc[k][2] << endl; + tree.getNearestNeighbours(xc[k], NGB, ngb, distances); + int last = -1; + + for (uint32_t i = 0; i < NGB; i++) + { + if (ngb[i] == 0) + continue; + + last = i; + + double d2 = 0; + for (int l = 0; l < 3; l++) + d2 += ({double delta = periodic(xc[k][l] - ngb[i]->coord[l]); delta*delta;}); + fngb << ngb[i]->coord[0] << " " << ngb[i]->coord[1] << " " << ngb[i]->coord[2] << " " << sqrt(d2) << " " << sqrt(distances[i]) << endl; + } + fngb << endl << endl; + double farther_dist = distances[last]; + for (uint32_t i = 0; i < Ncells; i++) + { + bool found = false; + // If the points is not in the list, it means it is farther than the farthest point + for (int j =0; j < NGB; j++) + { + if (&cells[i] == ngb[j]) { + found = true; + break; + } + } + double dist_to_seed = 0; + for (int l = 0; l < 3; l++) + { + double delta = periodic(xc[k][l]-cells[i].coord[l]); + dist_to_seed += delta*delta; + } + double delta_epsilon = fabs(farther_dist/dist_to_seed-1); + if (!found) + { + if (dist_to_seed <= farther_dist && + delta_epsilon > 1000*sqrt(numeric_limits::epsilon())) + { + cout << "SHOULD BE IN (d=" << dist_to_seed << "): Failed iteration " << k << " particle " << i << " (farthest=" << sqrt(farther_dist) << ")" << endl; + cout << "delta_epsilon = " << delta_epsilon << endl; + abort(); + } + } + else + { + if (dist_to_seed > farther_dist && + delta_epsilon > 1000*sqrt(numeric_limits::epsilon())) + { + cout << "SHOULD BE OUT (d=" << sqrt(dist_to_seed) << "): Failed iteration " << k << " particle " << i << " (farthest=" << sqrt(farther_dist) << ")" << endl; + cout << "delta_epsilon = " << delta_epsilon << " epsilon = " << 100*sqrt(numeric_limits::epsilon()) << endl; + abort(); + } + } + } + } + + return 0; +} diff --git a/src/replicateGenerator.hpp b/src/replicateGenerator.hpp index 8f52ace..3b6d3b8 100644 --- a/src/replicateGenerator.hpp +++ b/src/replicateGenerator.hpp @@ -20,6 +20,8 @@ namespace CosmoTool replicate = shift; numFaces = spower(3); std::copy(x, x+N, x_base); + if (!next()) + abort(); } bool next()