diff --git a/sample/testkd2.cpp b/sample/testkd2.cpp index 2b79408..bc65899 100644 --- a/sample/testkd2.cpp +++ b/sample/testkd2.cpp @@ -7,13 +7,14 @@ #include "mykdtree.hpp" #include "kdtree_splitters.hpp" -#define NTRY 100 +#define NTRY 10 #define ND 3 using namespace std; using namespace CosmoTool; typedef KDTree > MyTree; +//typedef KDTree MyTree; typedef KDCell MyCell; MyCell *findNearest(MyTree::coords& xc, MyCell *cells, uint32_t Ncells) @@ -39,7 +40,7 @@ MyCell *findNearest(MyTree::coords& xc, MyCell *cells, uint32_t Ncells) int main() { - uint32_t Ncells = 100000; + uint32_t Ncells = 10000000; MyCell *cells = new MyCell[Ncells]; for (int i = 0; i < Ncells; i++) @@ -49,8 +50,16 @@ int main() cells[i].coord[l] = drand48(); } + // Check timing + clock_t startTimer = clock(); MyTree tree(cells, Ncells); - + clock_t endTimer = clock(); + + 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; @@ -69,31 +78,39 @@ int main() for (int k = 0; k < NTRY; k++) { cout << "Seed = " << xc[k][0] << " " << xc[k][1] << " " << xc[k][2] << endl; tree.getNearestNeighbours(xc[k], 12, ngb, distances); + int last = -1; for (uint32_t i = 0; i < 12; i++) { + if (ngb[i] == 0) + continue; + + last = i; + double d2 = 0; for (int l = 0; l < 3; l++) d2 += ({double delta = xc[k][l] - ngb[i]->coord[l]; delta*delta;}); fngb << ngb[i]->coord[0] << " " << ngb[i]->coord[1] << " " << ngb[i]->coord[2] << " " << sqrt(d2) << endl; } - fngb << endl << endl; - double farther_dist = distances[11]; + 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 farther point - for (int j =0; j < 12; j++) + // If the points is not in the list, it means it is farther than the farthest point + for (int j =0; j < 12; j++) { if (&cells[i] == ngb[j]) { found = true; break; } } - double dist_to_seed = 0; - for (int l = 0; l < 3; l++) - { double delta = xc[k][l]-cells[i].coord[l]; - dist_to_seed += delta*delta; } + double dist_to_seed = 0; + for (int l = 0; l < 3; l++) + { + double delta = xc[k][l]-cells[i].coord[l]; + dist_to_seed += delta*delta; + } if (!found) { if (dist_to_seed <= farther_dist) diff --git a/src/kdtree_splitters.hpp b/src/kdtree_splitters.hpp index 4594298..4092cde 100644 --- a/src/kdtree_splitters.hpp +++ b/src/kdtree_splitters.hpp @@ -12,56 +12,119 @@ namespace CosmoTool typedef typename KDDef::KDCoordinates coords; typedef typename KDDef::CoordType ctype; + + void check_splitting(KDCell **cells, uint32_t Ncells, int axis, uint32_t split_index, ctype midCoord) + { + ctype delta = std::numeric_limits::max(); + assert(split_index < Ncells); + assert(axis < N); + for (uint32_t i = 0; i < split_index; i++) + { + assert(cells[i]->coord[axis] <= midCoord); + delta = min(midCoord-cells[i]->coord[axis], delta); + } + for (uint32_t i = split_index+1; i < Ncells; i++) + { + assert(cells[i]->coord[axis] > midCoord); + delta = min(cells[i]->coord[axis]-midCoord, delta); + } + assert(delta >= 0); + assert (std::abs(cells[split_index]->coord[axis]-midCoord) <= delta); + } + void operator()(KDCell **cells, uint32_t Ncells, uint32_t& split_index, int axis, coords minBound, coords maxBound) { + if (Ncells == 1) + { + split_index = 0; + return; + } + ctype midCoord = 0.5*(maxBound[axis]+minBound[axis]); uint32_t below = 0, above = Ncells-1; - ctype delta_max = std::abs(cells[0]->coord[axis]-midCoord); - uint32_t idx_max = 0; + ctype delta_min = std::numeric_limits::max(); + uint32_t idx_min = std::numeric_limits::max(); while (below < above) { ctype delta = cells[below]->coord[axis]-midCoord; if (delta > 0) { - if (delta < delta_max) + if (delta < delta_min) { - delta_max = delta; - idx_max = above; + delta_min = delta; + idx_min = above; } std::swap(cells[below], cells[above--]); } else { - if (-delta < delta_max) + if (-delta < delta_min) { - delta_max = -delta; - idx_max = below; + delta_min = -delta; + idx_min = below; } below++; } } - if (idx_max != above) + // Last iteration + { + ctype delta = cells[below]->coord[axis]-midCoord; + if (delta > 0) + { + if (delta < delta_min) + { + delta_min = delta; + idx_min = above; + } + } + else + { + if (-delta < delta_min) + { + delta_min = -delta; + idx_min = above; + } + } + } + + if (idx_min != above) { - bool cond1 = cells[idx_max]->coord[axis] > midCoord; + bool cond1 = cells[idx_min]->coord[axis] > midCoord; bool cond2 = cells[above]->coord[axis] > midCoord; if ((cond1 && cond2) || (!cond1 && !cond2)) { split_index = above; - std::swap(cells[above], cells[idx_max]); + std::swap(cells[above], cells[idx_min]); } else if (cond2) { - split_index = above-1; - std::swap(cells[above-1], cells[idx_max]); + if (above >= 1) + { + split_index = above-1; + std::swap(cells[above-1], cells[idx_min]); + } + else + split_index = 0; + assert(split_index >= 0); } else { - split_index = above+1; - std::swap(cells[above+1], cells[idx_max]); + if (above+1 < Ncells) + { + split_index = above+1; + std::swap(cells[above+1], cells[idx_min]); + } + else + split_index = Ncells-1; + + assert(split_index < Ncells); } } else split_index = above; + + + // check_splitting(cells, Ncells, axis, split_index, midCoord); } }; diff --git a/src/mykdtree.tcc b/src/mykdtree.tcc index 6f5eebf..a5e3abb 100644 --- a/src/mykdtree.tcc +++ b/src/mykdtree.tcc @@ -390,11 +390,11 @@ namespace CosmoTool { // If not it is in 1. go = node->children[1]; other = node->children[0]; - if (go == 0) - { - go = other; - other = 0; - } + // if (go == 0) + // { + // go = other; + //other = 0; + //} } if (go != 0) @@ -407,8 +407,8 @@ namespace CosmoTool { computeDistance(node->value, info.x); info.queue.push(node->value, thisR2); info.traversed++; - if (go == 0) - return; + // if (go == 0) + // return; // Now we found the best. We check whether the hypersphere // intersect the hyperplane of the other branch