Commits

Peter Shinners committed eccbc3b

Code is in working mode

Fixed problems with the bisecting and finding subset

Comments (0)

Files changed (4)

 #nrnr.a: nrnr.c
 
 test: nrnr.c test.c nrnr.h Makefile
-	$(CC) $(CFLAGS) -o test nrnr.c test.c
+	$(CC) $(CFLAGS) -O3 -o test nrnr.c test.c
 
 
 // start is inclusive, end is exclusive
 struct range {
-    const size_t* start;
-    const size_t* end;
+    size_t start;
+    size_t end;
 };
 
 
-#ifndef MIN
-#define MIN(a,b) ((a)<(b) ? (a) : (b))
-#endif
-
-
 
 void SortIndices(size_t* array, int numPoints, const float* points, size_t axis, size_t dimensions);
-struct range BisectIndices(const size_t* array, size_t numPoints, const float* points, size_t axis, size_t dimensions, float minVal, float maxVal);
+struct range BisectIndices(const struct NrNrHood* hood, float minVal, float maxVal, size_t axis);
 
 
 
         SortIndices(array, numPoints, points, d, dimensions);
     }
 
-    // Find inverse mapping to sort index (point to index)
+    // Find inverse mappifng to sort index (point to index)
     // This goes into a temporary set of inverse indices
 
     size_t* inverse = malloc(sizeof(size_t) * numPoints * dimensions);
         }
     }
 
-	printf("  invx:"); for (size_t i=0; i<MIN(numPoints,10); ++i) {printf(" %zu", inverse[i]);} printf("\n");
-	printf("  invy:"); for (size_t i=0; i<MIN(numPoints,10); ++i) {printf(" %zu", inverse[numPoints + i]);} printf("\n");
-	printf("  invz:"); for (size_t i=0; i<MIN(numPoints,10); ++i) {printf(" %zu", inverse[numPoints * 2 + i]);} printf("\n");
-
     // Build next arrays that map a point from one dimension to the next
     // Allows jumping from dim0 index N to dim1, and then to dim2, ...
 
     for (size_t d = 0; d < dimensions; ++d) {
         size_t* next = hood->nextIndices + numPoints * d;
-        size_t* ind = hood->indices + numPoints * d;
-        size_t* inv;
+        size_t* inv = inverse + numPoints * d;
+        size_t* ind;
         if (d) {
-            inv = inverse + numPoints * (d - 1);
+            ind = hood->indices + numPoints * (d - 1);
         } else {
-            inv = inverse + numPoints * (dimensions - 1);
+            ind = hood->indices + numPoints * (dimensions - 1);
         }
         for (size_t i = 0; i < numPoints; ++i) {
             next[i] = inv[ind[i]];
 
 size_t NrNrSearchBlock(const struct NrNrHood* hood, size_t** results, const float* minimum, const float* maximum) {
     // Find ranges of indices for each dimension
-
     struct range* ranges = alloca(sizeof(ranges) * hood->dimensions);
+
     for (size_t d = 0; d < hood->dimensions; ++d) {
         size_t* array = hood->indices + hood->numPoints * d;
-        ranges[d] = BisectIndices(array, hood->numPoints, hood->points, d, hood->dimensions, minimum[d], maximum[d]);
+        ranges[d] = BisectIndices(hood, minimum[d], maximum[d], d);
     }
 
     // Find dimension with the least number of matches
     size_t* common = (size_t*)malloc(sizeof(size_t) * minSize);
     size_t* found = common;
 
-#if 0
-    for (size_t* minStep = ranges[minAxis].start; minStep < ranges[minAxis].end; ++minStep) {
-        unsigned d;
-        size_t* workPt = hood->indices + hood->numPoints * minAxis - minStep;
+    for (size_t minStep = ranges[minAxis].start; minStep < ranges[minAxis].end; ++minStep) {
+        size_t d;
+        size_t workStep = minStep;
         for (d = 1; d < hood->dimensions; ++d) {
             size_t workAxis = minAxis + d;
             if (workAxis >= hood->dimensions) {
                 workAxis -= hood->dimensions;
             }
-
-            workPt = hood->nextIndices[hood->numPoints * workAxis + workPt];
-            if (workPt < ranges[workAxis].start || workPt > ranges[workAxis].end) {
+            size_t pt = hood->nextIndices[hood->numPoints * workAxis + workStep];
+            if (pt < ranges[workAxis].start || pt > ranges[workAxis].end) {
                 break;
             }
+            workStep = pt;
         }
         if (d == hood->dimensions) {
-            *found = *minStep;
+            *found = hood->indices[hood->numPoints * minAxis + minStep];
             ++found;
         }
     }
-#endif
 
     // Handle return values
 
 
 
 
-struct range BisectIndices(const size_t* array, size_t numPoints, const float* points, size_t axis, size_t dimensions, float minVal, float maxVal) {
-    const size_t* fullStart = array;
-    const size_t* fullEnd = fullStart + numPoints;
-    const size_t* start = fullStart;
-    const size_t* end = fullEnd;
+struct range BisectIndices(const struct NrNrHood* hood, float minVal, float maxVal, size_t axis) {
     struct range range;
+    size_t *fullStart, *fullEnd;
+    size_t *start;
+    size_t *end;
 
-    // Find min index
+    fullStart = start = hood->indices + hood->numPoints * axis;
+    fullEnd = end = fullStart + hood->numPoints;
+
+    // TODO, early out when maxVal < start point or minVal > end point
     while (start < end) {
         size_t dist = end - start;
-        const size_t* half = start + (dist / 2);
-        float halfValue = array[*half * 3 + axis];
-        if (minVal < halfValue) {
+        size_t* half = start + (dist / 2);
+        float halfValue = hood->points[*half * hood->dimensions + axis];
+        if (minVal <= halfValue) {
             end = half;
         } else {
             start = half + 1;
         }
     }
-    range.start = start;
+    range.start = start - fullStart;
 
-    // Find max index
+    start = end;
     end = fullEnd;
     while (start < end) {
         size_t dist = end - start;
-        const size_t* half = start + (dist / 2);
-        float halfValue = points[*half * 3 + axis];
-        if (maxVal > halfValue) {
+        size_t* half = start + (dist / 2);
+        float halfValue = hood->points[*half * hood->dimensions + axis];
+        if (maxVal >= halfValue) {
             start = half + 1;
         } else {
             end = half;
         }
     }
-    range.end = end;
+    range.end = start - fullStart;
 
     return range;
 }
-
+nrnr
+Create neighborhood: 15.615 sec
+50r = Search block 5120872: 0.296 sec
+40r = Search block 639449: 0.114 sec
+20r = Search block 79971: 0.047 sec
+10r = Search block 20031: 0.021 sec
+Cleanup: 0.012 sec
+3282352maxresident
+
+40r 4d: sig11
+40r 3d: .119 sec
+40r 2d: .057 sec
+40r 1d: .023 sec
+
+
+
+kdtree
+inserting 10000000 random vectors... 29.000 sec
+80r = range query returned 2680049 items in 1.02700 sec
+40r = range query returned 334901 items in 0.16200 sec
+20r = range query returned 41902 items in 0.02000 sec
+10r = range query returned 5175 items in 0.00300 sec
+cleaned up in 2.32100 sec
+3169120maxresident
+
+
+
+nanoflann
+Create Kdtree: 12.735sec
+80r = Search: 2680049, 0.511999sec
+40r = Search: 334901, 0.0760002sec
+20r = Search: 41902, 0.00699997sec
+10r = Search: 5175, 0.000999451sec
+Cleanup: 0.0279999sec
+1180688maxresident
+
+
 #include <time.h>
 
 
+#define DOCHECKS 0
+
 
 float GetTime(void);
 void CheckSorting(struct NrNrHood* hood);
 void CheckNextIndices(struct NrNrHood* hood);
-
-
-
-#ifndef MIN
-#define MIN(a,b) ((a)<(b) ? (a) : (b))
-#endif
+void CheckSearch(struct NrNrHood* hood, size_t numResults, const size_t* results, const float* minimum, const float* maximum);
 
 
 
 
     // Create random points
 
-    start = GetTime();
     srand(12345);
-    size_t numPoints = 5; //0000000;
+    size_t numPoints = 10000000;
+    start = GetTime();
     float* points = (float*)malloc(sizeof(float) * 3 * numPoints);
 	for(size_t i=0; i<numPoints; i++) {
 		float* p = points + i * 3;
 		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;
-
+		p[3] = ((float)rand() / RAND_MAX) * 200.0 - 100.0;
+        //printf("  pt %zu: %.2f %.2f %.2f\n",i, p[0], p[1], p[2]);
     }
 	printf("Create random points (%zu): %.3f sec\n", numPoints, GetTime() - start);
 
     struct NrNrHood* hood = NrNrCreateNeighborhood(numPoints, 3, points);
 	printf("Create neighborhood: %.3f sec\n", GetTime() - start);
 
-
-	printf("  x:"); for (size_t i=0; i<MIN(numPoints,10); ++i) {printf(" %zu=%.3f", hood->indices[i], points[hood->indices[i] * 3]);} printf("\n");
-	printf("  y:"); for (size_t i=0; i<MIN(numPoints,10); ++i) {printf(" %zu=%.3f", hood->indices[numPoints + i], points[hood->indices[numPoints + i] * 3 + 1]);} printf("\n");
-	printf("  z:"); for (size_t i=0; i<MIN(numPoints,10); ++i) {printf(" %zu=%.3f", hood->indices[numPoints * 2 + i], points[hood->indices[numPoints * 2 + i] * 3 + 2]);} printf("\n");
+#if DOCHECKS
+    printf("Checking sorting\n");
     CheckSorting(hood);
-
-	printf(" nx:"); for (size_t i=0; i<MIN(numPoints,10); ++i) {printf(" %zu", hood->nextIndices[i]);} printf("\n");
-	printf(" ny:"); for (size_t i=0; i<MIN(numPoints,10); ++i) {printf(" %zu", hood->nextIndices[numPoints + i]);} printf("\n");
-	printf(" nz:"); for (size_t i=0; i<MIN(numPoints,10); ++i) {printf(" %zu", hood->nextIndices[numPoints * 2 + i]);} printf("\n");
+    printf("Checking nexts\n");
     CheckNextIndices(hood);
+#endif
 
     // Search block
 
     start = GetTime();
-    float minPos[3] = {0.0f, 20.0f, -20.0f};
-    float maxPos[3] = {minPos[0] + 40.0, minPos[1] + 40.0, minPos[2] + 40.0};
+    float minPos[3] = {-40.0f,-20.0f,-60.0f};
+    float maxPos[3] = { 40.0,  60.0f, 20.0f};
+//    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 block %zu: %.3f sec\n", count, GetTime() - start);
+#if DOCHECKS
+    CheckSearch(hood, count, results, minPos, maxPos);
+#endif
 
     // Cleanup
 
             }
         }
     }
+}
 
+
+void CheckSearch(struct NrNrHood* hood, size_t numResults, const size_t* results, const float* minimum, const float* maximum) {
+    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;
+                break;
+            }
+        }
+
+        int inresults = 0;
+        for (size_t j = 0; j < numResults; ++j) {
+            if (results[j] == i) {
+                inresults = 1;
+                break;
+            }
+        }
+        if (inside && !inresults) {
+            printf("Error with results, points should be inside: index=%zu, %.2f %.2f %.2f\n", i,
+                        hood->points[i * hood->dimensions], hood->points[i * hood->dimensions + 1], hood->points[i * hood->dimensions + 2]);
+            exit(1);
+        } else if (!inside && inresults) {
+            printf("Error with results, points should not be inside: index=%zu, %.2f %.2f %.2f\n", i,
+                        hood->points[i * hood->dimensions], hood->points[i * hood->dimensions + 1], hood->points[i * hood->dimensions + 2]);
+            exit(1);
+        }
+    }
 }