Peter Shinners avatar Peter Shinners committed 24e9586

Cleaned up the finished implementation and test code

Comments (0)

Files changed (4)

 #include <string.h>
 #include <alloca.h>
 #include <stdio.h> // temp
+#include <limits.h>
 #include "nrnr.h"
 
 
+
 // start is inclusive, end is exclusive
 struct range {
-    index_t start;
-    index_t end;
+    unsigned short start;
+    unsigned short end;
 };
 
 
 
-static void SortIndices(index_t* array, index_t numPoints, const float* points, size_t axis, size_t dimensions);
-static struct range BisectIndices(const index_t* array, const float* points, index_t numPoints, size_t dimensions, float minVal, float maxVal, size_t axis);
-
-
-#include <sys/time.h>
-#include <time.h>
-static float GetTime(void)
-{
-	static struct timeval timeval, first_timeval;
-	gettimeofday(&timeval, 0);
-	if(first_timeval.tv_sec == 0) {
-		first_timeval = timeval;
-		return 0;
-	}
-	unsigned int s = (timeval.tv_sec - first_timeval.tv_sec) * 1000 + (timeval.tv_usec - first_timeval.tv_usec) / 1000.0;
-    return s / 1000.0f;
-}
+static void SortIndices(unsigned short* array, unsigned short numPoints, const float* points, size_t axis, size_t dimensions);
+static struct range BisectIndices(const unsigned short* array, const float* points, unsigned short numPoints, size_t dimensions, float minVal, float maxVal, size_t axis);
 
 
 
     // All data is allocated in a single block. Allocate and then
     // aim some pointers into the space.
 
-    size_t allocSize = sizeof(struct NrNrHood) + sizeof(index_t) * numPoints * 6;
+    size_t allocSize = sizeof(struct NrNrHood) + sizeof(unsigned short) * numPoints * 6;
     struct NrNrHood* hood = malloc(allocSize);
     hood->numPoints = numPoints;
     hood->dimensions = dimensions;
     hood->points = points;
-    hood->chunks = (numPoints + NRNRCHUNK - 1) / NRNRCHUNK;
-    hood->indices = (index_t*)(hood + 1);
+    hood->chunks = (numPoints + USHRT_MAX - 1) / USHRT_MAX;
+    hood->indices = (unsigned short*)(hood + 1);
     hood->nextIndices = hood->indices + numPoints * 3;
 
     // Allocate temporary buffer for storing the inverse sorted array
     // Find inverse mappifng to sort index (point to index)
     // This goes into a temporary set of inverse indices
-    size_t inverseSize = NRNRCHUNK;
+    size_t inverseSize = USHRT_MAX;
     if (inverseSize > numPoints) {
         inverseSize = numPoints;
     }
-    index_t* inverse = (index_t*)malloc(sizeof(index_t) * inverseSize);
+    unsigned short* inverse = (unsigned short*)malloc(sizeof(unsigned short) * inverseSize);
 
     // Break the data into small chunks
-    index_t* chunkIndices = hood->indices;
-    index_t* chunkNextIndices = hood->nextIndices;
+    unsigned short* chunkIndices = hood->indices;
+    unsigned short* chunkNextIndices = hood->nextIndices;
     const float* chunkPoints = points;
-    size_t indexStep = hood->dimensions * NRNRCHUNK;
+    size_t indexStep = hood->dimensions * USHRT_MAX;
     for (size_t c = 0; c < hood->chunks; ++c,
                 chunkIndices += indexStep, chunkNextIndices += indexStep, chunkPoints += indexStep) {
-        size_t chunkSize = numPoints - (c * NRNRCHUNK);
-        if (chunkSize > NRNRCHUNK) {
-            chunkSize = NRNRCHUNK;
+        size_t chunkSize = numPoints - (c * USHRT_MAX);
+        if (chunkSize > USHRT_MAX) {
+            chunkSize = USHRT_MAX;
         }
-//printf("CHUNK %zu   chunkSize=%zu\n", c, chunkSize);
 
-        // Initialize indices with list of all points
-        for (index_t i = 0; i < chunkSize; ++i) {
+        for (unsigned short i = 0; i < chunkSize; ++i) {
             chunkIndices[i] = i;
         }
-//printf("presorted 0:"); for (size_t z=0; z < chunkSize; ++z) printf(" %d", (int)chunkIndices[z]); printf("\n");
 
         // Copy ordered point numbers from dimension 0 into all others
         for (size_t d = 1; d < dimensions; ++d) {
-            index_t* initSrc = chunkIndices + (d - 1) * chunkSize;
-            index_t* initDst = chunkIndices + d * chunkSize;
-            memcpy(initDst, initSrc, sizeof(index_t) * chunkSize);
-//printf("presorted %zu:", d); for (size_t z=0; z < chunkSize; ++z) printf(" %zu", initDst[z]); printf("\n");
+            unsigned short* initSrc = chunkIndices + (d - 1) * chunkSize;
+            unsigned short* initDst = chunkIndices + d * chunkSize;
+            memcpy(initDst, initSrc, sizeof(unsigned short) * chunkSize);
         }
 
         // Sort indices for each dimension
         for (size_t d = 0; d < dimensions; ++d) {
-            index_t* array = chunkIndices + d * chunkSize;
+            unsigned short* array = chunkIndices + d * chunkSize;
             SortIndices(array, chunkSize, chunkPoints, d, dimensions);
-//printf("   sorted %zu:", d); for (size_t z=0; z < chunkSize; ++z) printf("  %zu(%.1f)", array[z], chunkPoints[array[z] * dimensions + d]); printf("\n");
-//printf("   sorted %d:", (int)d); for (size_t z=0; z < chunkSize; ++z) printf(" %d", (int)array[z]); printf("\n");
         }
 
         // Build next arrays that map a point from one dimension to the next
         for (size_t d = 0; d < dimensions; ++d) {
 
             // Build an index of the 'inverse' mapping from point number to sorted index
-            index_t* ind = chunkIndices + chunkSize * d;
+            unsigned short* ind = chunkIndices + chunkSize * d;
 
-            for (index_t i = 0; i < chunkSize; ++i) {
+            for (unsigned short i = 0; i < chunkSize; ++i) {
                 inverse[ind[i]] = i;
             }
-//printf("     inverse %zu:", d); for (size_t z=0; z < chunkSize; ++z) printf(" %zu", inverse[z]); printf("\n");
-
 
             // Fill out the next array, move 'ind' to the array coming from
-            index_t* next = chunkNextIndices + chunkSize * d;
+            unsigned short* next = chunkNextIndices + chunkSize * d;
             if (d) {
                 ind = chunkIndices + chunkSize * (d - 1);
             } else {
                 ind = chunkIndices + chunkSize * (dimensions - 1);
             }
-            for (index_t i = 0; i < chunkSize; ++i) {
+            for (unsigned short i = 0; i < chunkSize; ++i) {
                 next[i] = inverse[ind[i]];
             }
-//printf("        next %zu:", d); for (size_t z=0; z < chunkSize; ++z) printf(" %zu", next[z]); printf("\n");
         }
-
     }
 
     free(inverse);
 
 
 
-//size_t NrNrCountBlock(struct NrNrHood* hood, size_t most, float* minimum, float* maximum) {
-//}
-
-
 size_t NrNrSearchBlock(const struct NrNrHood* hood, size_t** results, const float* minimum, const float* maximum) {
     struct range* ranges = alloca(sizeof(ranges) * hood->dimensions);
 
     size_t* workAxis = alloca(sizeof(size_t) * hood->dimensions);;
-    index_t** workNexts = alloca(sizeof(index_t*) * hood->dimensions);
+    unsigned short** workNexts = alloca(sizeof(unsigned short*) * hood->dimensions);
     workAxis[0] = hood->dimensions; // unused (will crash if used)
     workNexts[0] = NULL; // unused (will crash if used)
 
     size_t* found = common;
 
     // Break the data into small chunks
-    index_t* chunkIndices = hood->indices;
-    index_t* chunkNextIndices = hood->nextIndices;
+    unsigned short* chunkIndices = hood->indices;
+    unsigned short* chunkNextIndices = hood->nextIndices;
     const float* chunkPoints = hood->points;
     size_t chunkOffset = 0;
-//printf("Total chunks: %zu\n", hood->chunks);
-    size_t indexStep = hood->dimensions * NRNRCHUNK;
+    size_t totalEmptyChunks = 0;
+    size_t indexStep = hood->dimensions * USHRT_MAX;
     for (size_t c = 0; c < hood->chunks; ++c,
-                chunkIndices += indexStep, chunkNextIndices += indexStep, chunkPoints += indexStep, chunkOffset += NRNRCHUNK) {
-        size_t chunkSize = hood->numPoints - (c * NRNRCHUNK);
+                chunkIndices += indexStep, chunkNextIndices += indexStep, chunkPoints += indexStep, chunkOffset += USHRT_MAX) {
+        size_t chunkSize = hood->numPoints - (c * USHRT_MAX);
         size_t chunkFound = 0;
-        if (chunkSize > NRNRCHUNK) {
-            chunkSize = NRNRCHUNK;
+        if (chunkSize > USHRT_MAX) {
+            chunkSize = USHRT_MAX;
         }
-//printf("Search Chunk %zu  size=%d\n", c, (int)chunkSize);
 
         // Find ranges of indices for each dimension
-
         int emptyBisect = 0;
         for (size_t d = 0; d < hood->dimensions; ++d) {
-            index_t* array = chunkIndices + chunkSize * d;
+            unsigned short* array = chunkIndices + chunkSize * d;
             ranges[d] = BisectIndices(array, chunkPoints, chunkSize, hood->dimensions, minimum[d], maximum[d], d);
-//printf("  bisect %zu: %d-%d, ", d, (int)ranges[d].start, (int)ranges[d].end);
-//for (size_t z=ranges[d].start; z<ranges[d].end; ++z) printf("  %zu(%.1f)", array[z] /*+ chunkOffset*/, chunkPoints[array[z] * hood->dimensions + d]); printf("\n");
-//for (size_t z=ranges[d].start; z<ranges[d].end; ++z) printf(" %d", (int)array[z] /*+ chunkOffset*/); printf("\n");
             if (ranges[d].start == ranges[d].end) {
                 emptyBisect = 1;
-//printf("   chunk found %zu: 0 (empty bisect %zu)\n", c, d);
+                totalEmptyChunks++;
                 break;
             }
         }
         }
 
         // Find dimension with the least number of matches
-
         size_t minAxis = 0;
-        index_t minSize = ranges[0].end - ranges[0].start;
+        unsigned short minSize = ranges[0].end - ranges[0].start;
         for (size_t d = 0; d < hood->dimensions; ++d) {
-            index_t size = ranges[d].end - ranges[d].start;
+            unsigned short size = ranges[d].end - ranges[d].start;
             if (size < minSize) {
                 minSize = size;
                 minAxis = d;
             }
         }
-
         // To eliminate as much code as possible for the search loop,
         // precompute several needed values
 
+        unsigned short* minIndices = chunkIndices + chunkSize * minAxis;
         for (size_t d = 1; d < hood->dimensions; ++d) {
             workAxis[d] = minAxis + d;
             if (workAxis[d] >= hood->dimensions) {
         }
 
         // The main search loop. All searching time is spent here.
-        index_t* minIndices = chunkIndices + chunkSize * minAxis;
 
-        for (index_t minStep = ranges[minAxis].start; minStep < ranges[minAxis].end; ++minStep) {
-            index_t workStep = minStep;
+        for (unsigned short minStep = ranges[minAxis].start; minStep < ranges[minAxis].end; ++minStep) {
+            unsigned short workStep = minStep;
             int match = 1;
-//printf("     axis=%zu  index=%d\n", minAxis, (int)minStep);
             for (size_t d = 1; d < hood->dimensions; ++d) {
-                index_t pt = workNexts[d][workStep];
-//printf("         axis=%zu  index=%d   (%d-%d)\n", workAxis[d], (int)pt, (int)ranges[workAxis[d]].start, (int)ranges[workAxis[d]].end);
+                unsigned short pt = workNexts[d][workStep];
                 if (!(match = (pt >= ranges[workAxis[d]].start && pt < ranges[workAxis[d]].end)))
                     break;
                 workStep = pt;
                 ++chunkFound;
             }
         }
-//printf("   chunk found %zu: %zu (min axis was %d)\n", c, chunkFound, (int)minSize);
     }
 
     // Handle return values
         *results = NULL;
         return 0;
     }
-    // realloc common to the smaller, found size
+
+    // Realloc common to the smaller, found size
     if (totalFound < hood->numPoints) {
         common = realloc(common, sizeof(size_t) * totalFound);
     }
 
 
 
-void SortIndices(index_t* array, index_t numPoints, const float* points, size_t axis, size_t dimensions) {
+void SortIndices(unsigned short* array, unsigned short numPoints, const float* points, size_t axis, size_t dimensions) {
     if (numPoints < 2) {
         return;
     }
 
     // Insertion sort for small ranges
     if (numPoints < 16) {
-        for (index_t i = 1; i < numPoints; ++i) {
-            index_t point = array[i];
+        for (unsigned short i = 1; i < numPoints; ++i) {
+            unsigned short point = array[i];
             float position = points[point * dimensions + axis];
-            index_t j;
+            unsigned short j;
             for (j = i; j > 0 && points[array[j - 1] * dimensions + axis] > position; --j) {
                 array[j] = array[j - 1];
             }
 
     // Quicksort (recursive) for most spans
     float position = points[array[numPoints / 2] * dimensions + axis];
-    index_t *left = array;
-    index_t *right = array + numPoints - 1;
+    unsigned short *left = array;
+    unsigned short *right = array + numPoints - 1;
     while (left <= right) {
         while (points[(*left) * dimensions + axis] < position)
             left++;
         while (points[(*right) * dimensions + axis] > position)
             right--;
         if (left <= right) {
-            index_t t = *left;
+            unsigned short t = *left;
             *left++ = *right;
             *right-- = t;
         }
 
 
 
-struct range BisectIndices(const index_t* array, const float* points, index_t numPoints, size_t dimensions, float minVal, float maxVal, size_t axis) {
+struct range BisectIndices(const unsigned short* array, const float* points, unsigned short numPoints, size_t dimensions, float minVal, float maxVal, size_t axis) {
     struct range range;
-    const index_t *fullStart, *fullEnd;
-    const index_t *start;
-    const index_t *end;
+    const unsigned short *fullStart, *fullEnd;
+    const unsigned short *start;
+    const unsigned short *end;
 
     fullStart = start = array;
     fullEnd = end = fullStart + numPoints;
-#if 0
+
     float startValue = points[array[0] * dimensions + axis];
     if (startValue > maxVal) {
         range.start = range.end = numPoints;
         return range;
     }
-#endif
 
     float finalValue = points[array[numPoints - 1] * dimensions + axis];
-#if 1
     if (finalValue < minVal) {
         range.start = range.end = 0;
         return range;
     }
-#endif
 
     while (start < end) {
-        index_t dist = end - start;
-        const index_t* half = start + (dist / 2);
+        unsigned short dist = end - start;
+        const unsigned short* half = start + (dist / 2);
         float halfValue = points[*half * dimensions + axis];
         if (minVal <= halfValue) {
             end = half;
-#if 0
             if (finalValue < halfValue) {
                 fullEnd = half;
                 finalValue = points[*fullEnd * dimensions + axis];
             }
-#endif
         } else {
             start = half + 1;
         }
     start = end;
     end = fullEnd;
     while (start < end) {
-        index_t dist = end - start;
-        const index_t* half = start + (dist / 2);
+        unsigned short dist = end - start;
+        const unsigned short* half = start + (dist / 2);
         float halfValue = points[*half * dimensions + axis];
         if (maxVal >= halfValue) {
             start = half + 1;
 
 
 #include <stdlib.h>
-#include <inttypes.h>
 
 
-//#define NRNRCHUNK UINT16_MAX
-#define NRNRCHUNK 32000
-
-
-typedef uint16_t index_t;
 
 struct NrNrHood {
     size_t numPoints;
     size_t dimensions;
     size_t chunks;
     const float* points;
-    index_t* indices;
-    index_t* nextIndices;
+    unsigned short* indices;
+    unsigned short* nextIndices;
 };
 
 
 struct NrNrHood* NrNrCreateNeighborhood(size_t numPoints, size_t dimensions, const float* points);
 struct NrNrHood* NrNrFreeNeighborhood(struct NrNrHood* hood);
 
-//size_t NrNrCountBlock(struct NrNrHood* hood, size_t most, float* minimum, float* maximum);
 size_t NrNrSearchBlock(const struct NrNrHood* hood, size_t** results, const float* minimum, const float* maximum);
 size_t* NrNrFreeResults(size_t* results);
 
 
 Moving uint16 from 65k chunks to 32k gave speedup to larger search radius, but had no
 effect on the 10r and 20r searches.
+
+97% of the search time is spent in the main search loop.
+No need to focus on optimizing the bisects.
+
+Searching 2 out of 3 dimensions still takes up most of the time.
+This is good, because it confirms intersecting the first two axis
+takes the majority of time.
+Single axis finds: 3976488  .017 seconds
+Second axis finds: 1594896  .042 seconds
+All three finds:    638857  .054 seconds
+
+
+with uint8 and 255 chunks, .026 time is bisecting and .048 is spent in searching (64%)
+with uint8 and 64 chunks,  .052 time is bisecting and .048 is spent searching (48%)
+Note that there are NONE (0) chunks that end up with 0 points in any axis. LAMEO! (40 radius)
 #include <stdio.h>
-#include "nrnr.h"
-
-
 #include <sys/time.h>
 #include <time.h>
+#include "nrnr.h"
 
 
-#define DOCHECKS 0
-
 
 float GetTime(void);
-void CheckSearch(struct NrNrHood* hood, size_t numResults, const size_t* results, const float* minimum, const float* maximum);
+void RunSearch(int repear, const struct NrNrHood* hood, const float* minimum, const float* maximum);
+void CheckSearch(const struct NrNrHood* hood, size_t numResults, const size_t* results, const float* minimum, const float* maximum);
 
 
 
     float start;
 
     // Create random points
-
     srand(1234);
-    size_t numPoints = 10000000;
+    size_t numPoints = 1000000; // million
     start = GetTime();
     float* points = (float*)malloc(sizeof(float) * 3 * numPoints);
 	for(size_t i=0; i<numPoints; i++) {
 		p[0] = ((float)rand() / RAND_MAX) * 200.0 - 100.0;
 		p[1] = ((float)rand() / RAND_MAX) * 200.0 - 100.0;
 		p[2] = ((float)rand() / RAND_MAX) * 200.0 - 100.0;
-        //printf("  pt %zu: %.1f %.1f %.1f\n",i, p[0], p[1], p[2]);
     }
-	//printf("Create random points (%zu): %.3f sec\n", numPoints, GetTime() - start);
 
     // Create neighborhood
-
     start = GetTime();
     struct NrNrHood* hood = NrNrCreateNeighborhood(numPoints, 3, points);
-	printf("Create neighborhood: %.3f sec\n", GetTime() - start);
+	printf("Create neighborhood: %.2f msec\n", GetTime() - start);
 
-    // Search block
-    start = GetTime();
-#if 0  // 10r
-    float minPos[3] = {-10.0f,10.0f,-30.0f};
-    float maxPos[3] = { 10.0,  30.0f, -10.0f};
-#elif 0 // 20r
-    float minPos[3] = {-20.0f, 0.0f, -40.0f};
-    float maxPos[3] = { 20.0,  40.0f, 00.0f};
-#elif 1 // 40r
-    float minPos[3] = {-40.0f,-20.0f,-60.0f};
-    float maxPos[3] = { 40.0,  60.0f, 20.0f};
-#elif 1 // 80r
-    float minPos[3] = {-80.0f,-60.0f,-100.0f};
-    float maxPos[3] = { 80.0,  100.0f, 60.0f};
-#endif
-
-    //printf("Searching: %.1f - %.1f,  %.1f - %.1f, %.1f - %.1f\n", minPos[0], maxPos[0], minPos[1], maxPos[1], minPos[2],maxPos[2]);
-
-    size_t* results;
-    size_t count = NrNrSearchBlock(hood, &results, minPos, maxPos);
-    printf("Search complete %zu: %.3f sec\n", count, GetTime() - start);
-
-#if 0
-for(size_t z = 0; z < count; ++z) {
-    float* pos = points + results[z] * 3;
-    printf("  found %zu: %.1f %.1f %.1f\n", results[z], pos[0], pos[1], pos[2]);
-}
-#endif
+    // Radius 1 search
+    float minPos1[3] = {-49.0f, -8.0f, 59.0f};
+    float maxPos1[3] = {-47.0f, -6.0f, 61.0f};
+    RunSearch(12000, hood, minPos1, maxPos1);
 
-#if DOCHECKS
-    CheckSearch(hood, count, results, minPos, maxPos);
-#endif
+    // Radius 10 search
+    float minPos10[3] = {-10.0f,10.0f,-30.0f};
+    float maxPos10[3] = { 10.0,  30.0f, -10.0f};
+    RunSearch(3000, hood, minPos10, maxPos10);
 
-    // Cleanup
+    float minPos20[3] = {-20.0f, 0.0f, -40.0f};
+    float maxPos20[3] = { 20.0,  40.0f, 00.0f};
+    RunSearch(1700, hood, minPos20, maxPos20);
 
-    start = GetTime();
-    results = NrNrFreeResults(results);
+    float minPos40[3] = {-40.0f,-20.0f,-60.0f};
+    float maxPos40[3] = { 40.0,  60.0f, 20.0f};
+    RunSearch(800, hood, minPos40, maxPos40);
 
+    float minPos80[3] = {-80.0f,-60.0f,-100.0f};
+    float maxPos80[3] = { 80.0,  100.0f, 60.0f};
+    RunSearch(500, hood, minPos80, maxPos80);
+
+    // Cleanup
+    start = GetTime();
     free(points);
     hood = NrNrFreeNeighborhood(hood);
-    printf("Cleanup: %.3f sec\n", GetTime() - start);
+    printf("Cleanup: %.2f msec\n", GetTime() - start);
 
 	return 0;
 }
 
 
 
-
-
-float GetTime(void)
-{
-	static struct timeval timeval, first_timeval;
-	gettimeofday(&timeval, 0);
-	if(first_timeval.tv_sec == 0) {
-		first_timeval = timeval;
-		return 0;
-	}
-	unsigned int s = (timeval.tv_sec - first_timeval.tv_sec) * 1000 + (timeval.tv_usec - first_timeval.tv_usec) / 1000.0;
-    return s / 1000.0f;
-}
-
-
-
-#if 0
-void CheckSorting(struct NrNrHood* hood) {
-    // Check value ordering
-    for (size_t d = 0; d < hood->dimensions; ++d) {
-        size_t* indices = hood->indices + hood->numPoints * d;
-        float val = hood->points[indices[0] * 3 + d];
-        for (size_t pt = 1; pt < hood->numPoints; ++pt) {
-            float nextVal = hood->points[indices[pt] * 3 + d];
-            if (nextVal < val) {
-                printf("Error sorting, values not in ascending order: axis=%zu, index=%zu\n", d, pt);
-                exit(1);
-            }
-            val = nextVal;
-        }
+void RunSearch(int repeat, const struct NrNrHood* hood, const float* minimum, const float* maximum) {
+    size_t count;
+    size_t* results = NULL;
+    float start = GetTime();
+    for (int r = 0; r < repeat; ++r) {
+        results = NrNrFreeResults(results);
+        count = NrNrSearchBlock(hood, &results, minimum, maximum);
     }
 
-    // Check repeated indices
-    for (size_t d = 0; d < hood->dimensions; ++d) {
-        size_t* indices = hood->indices + hood->numPoints * d;
-        for (size_t i = 0; i < hood->numPoints; ++i) {
-            size_t pt = indices[i];
-            if (pt >= hood->numPoints) {
-                printf("Error sorting, point number out of range: axis=%zu, index=%zu\n", d, i);
-                exit(1);
-            }
-            for (size_t j = i + 1; j < hood->numPoints; ++j) {
-                if (pt == indices[j]) {
-                    printf("Error sorting, point number used more than once: axis=%zu, index=%zu\n", d, i);
-                    exit(1);
-                }
-            }
-        }
-    }
-}
-
-
-void CheckNextIndices(struct NrNrHood* hood) {
-    // Check repeated indices
-    for (size_t d = 0; d < hood->dimensions; ++d) {
-        size_t* next = hood->nextIndices + hood->numPoints * d;
-        for (size_t i = 0; i < hood->numPoints; ++i) {
-            size_t pt = next[i];
-            if (pt >= hood->numPoints) {
-                printf("Error with nexts, point number out of range: axis=%zu, index=%zu\n", d, i);
-                exit(1);
-            }
-            for (size_t j = i + 1; j < hood->numPoints; ++j) {
-                if (pt == next[j]) {
-                    printf("Error with nexts, point number used more than once: axis=%zu, index=%zu\n", d, i);
-                    exit(1);
-                }
-            }
-        }
-    }
+    float radius = (maximum[0] - minimum[0]) / 2.0f;
+    printf("Search found %zu: %.2f msec per search  (repeat=%d radius=%.1f)\n",
+                count, (GetTime() - start) / repeat, repeat, radius);
 
-    // Check roundtrip through dimensions ends with the same
-    for (size_t d = 1; d < hood->dimensions; ++d) {
-        size_t* next = hood->nextIndices + hood->numPoints * d;
-        size_t* indices = hood->indices + hood->numPoints * d;
-        size_t* from;
-        if (d) {
-            from = hood->indices + hood->numPoints * (d - 1);
-        } else {
-            from = hood->indices + hood->numPoints * (hood->dimensions - 1);
-        }
-        for (size_t i = 0; i < hood->numPoints; ++i) {
-            size_t pt = from[i];
-            if (pt != indices[next[i]]) {
-                printf("Error with nexts, incorrect index: axis=%zu, index=%zu\n", d, i);
-                exit(1);
-            }
-        }
+    // Checking is too slow with big pools of matches
+    if (count < 10000) {
+        CheckSearch(hood, count, results, minimum, maximum);
     }
+    results = NrNrFreeResults(results);
 }
-#endif
 
 
-void CheckSearch(struct NrNrHood* hood, size_t numResults, const size_t* results, const float* minimum, const float* maximum) {
+
+void CheckSearch(const struct NrNrHood* hood, size_t numResults, const size_t* results, const float* minimum, const float* maximum) {
     int foundError = 0;
     for (size_t i = 0; i < hood->numPoints; ++i) {
         const float* pos = hood->points + i * hood->dimensions;
         int inside = 1;
+
         for (size_t d = 0; d < hood->dimensions; ++d) {
             if (pos[d] < minimum[d] || pos[d] > maximum[d]) {
                 inside = 0;
         exit(1);
     }
 }
+
+
+
+// Return time in milliseconds
+float GetTime(void)
+{
+	static struct timeval timeval, first_timeval;
+	gettimeofday(&timeval, 0);
+	if(first_timeval.tv_sec == 0) {
+		first_timeval = timeval;
+		return 0;
+	}
+	float s = (timeval.tv_sec - first_timeval.tv_sec) * 1000 + (timeval.tv_usec - first_timeval.tv_usec) / 1000.0;
+    return s;
+}
+
+
+
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.