Changeset 28223

Timestamp:

Jun 4, 2010, 2:20:19 PM (16 years ago)

Author:

Paul Price

Message:

Add extra diagnostics.

File:

• trunk/psLib/test/types/tap_psTree.c (modified) (5 diffs)

trunk/psLib/test/types/tap_psTree.c

@@ -3 +3 @@
 #include "pstap.h"
 
-#define NUM 1000000                      // Number of points
+#define NUM 1000000                     // Number of points
 #define SPHERICAL_DISTANCE 3.0          // Distance of interest for spherical test
+#define SEED 0                          // Random seed
+
+static double distance(double ra1, double dec1, double ra2, double dec2)
+{
+#if 0
+    // Traditional formula
+    return acos(sin(dec1) * sin(dec2) + cos(dec1) * cos(dec2) * cos(ra1 - ra2));
+#else
+    // Haversine formula: used in psTree
+    double dphi = dec1 - dec2;
+    double sindphi = sin(dphi/2.0);
+    double dlambda = ra1 - ra2;
+    double sindlambda = sin(dlambda/2.0);
+    return 2.0 * asin(sqrt(PS_SQR(sindphi) + cos(dec1) * cos(dec2) * PS_SQR(sindlambda)));
+#endif
+}
 
 int main(int argc, char *argv[])
@@ -54 +70 @@
             ok(closest->data.F64[0] == x->data.F64[bestIndex] &&
                closest->data.F64[1] == y->data.F64[bestIndex],
-               "correst coords: %lf,%lf(%lf) vs %lf,%lf(%lf)",
+               "correct coords: %lf,%lf(%lf) vs %lf,%lf(%lf)",
                closest->data.F64[0], closest->data.F64[1],
                sqrt(PS_SQR(closest->data.F64[0]) + PS_SQR(closest->data.F64[1])),
@@ -77 +93 @@
         psVector *dec = psVectorAlloc(NUM, PS_TYPE_F64);
 
-        psRandom *rng = psRandomAllocSpecific(PS_RANDOM_TAUS, 0);
+        psRandom *rng = psRandomAllocSpecific(PS_RANDOM_TAUS, SEED);
         for (int i = 0; i < NUM; i++) {
             // Using http://mathworld.wolfram.com/SpherePointPicking.html
-            ra->data.F64[i] = psRandomUniform(rng);
+            ra->data.F64[i] = psRandomUniform(rng) * 2.0 * M_PI;
             dec->data.F64[i] = acos(2.0 * psRandomUniform(rng) - 1.0) - M_PI_2;
         }
@@ -92 +108 @@
 
             psVector *coords = psVectorAlloc(2, PS_TYPE_F64);
-            coords->data.F64[0] = psRandomUniform(rng);
+#if 0
+            coords->data.F64[0] = psRandomUniform(rng) * 2.0 * M_PI;
             coords->data.F64[1] = acos(2.0 * psRandomUniform(rng) - 1.0) - M_PI_2;
+#else
+            psVectorInit(coords, 0);
+#endif
 
             psVector *indices = psTreeAllWithin(tree, coords, DEG_TO_RAD(SPHERICAL_DISTANCE));
@@ -102 +122 @@
             ok(psVectorSortInPlace(indices), "sorted indices");
 
+#if 0
+            for (long i = 0; i < indices->n; i++) {
+                long index = indices->data.S64[i];
+                double dist = distance(coords->data.F64[0], coords->data.F64[1],
+                                       ra->data.F64[index], dec->data.F64[index]);
+                diag("%ld (%lf,%lf) is in the list (%lf vs %lf)",
+                     index, ra->data.F64[index], dec->data.F64[index], RAD_TO_DEG(dist), SPHERICAL_DISTANCE);
+            }
+#endif
+
             bool allgood = true;        // All points in the appropriate place?
             double bestDistance = INFINITY; // Distance to best point
             long bestIndex = -1;        // Index of best point
+            long bad = 0;               // Number bad
             for (long i = 0, j = 0; i < NUM; i++) {
-#if 1
-                // Traditional formula
-                double dist = acos(sin(coords->data.F64[1]) * sin(dec->data.F64[i]) +
-                                   cos(coords->data.F64[1]) * cos(dec->data.F64[i]) *
-                                   cos(coords->data.F64[0] - ra->data.F64[i]));
-#else
-                // Haversine formula: used in psTree
-                double dphi = coords->data.F64[1] - dec->data.F64[i];
-                double sindphi = sin(dphi/2.0);
-                double dlambda = coords->data.F64[0] - ra->data.F64[i];
-                double sindlambda = sin(dlambda/2.0);
-                double dist = 2.0 * asin(sqrt(PS_SQR(sindphi) +
-                                              cos(coords->data.F64[1]) * cos(dec->data.F64[i]) *
-                                              PS_SQR(sindlambda)));
-#endif
-                                              if (dist < bestDistance) {
+                double dist = distance(coords->data.F64[0], coords->data.F64[1],
+                                       ra->data.F64[i], dec->data.F64[i]);
+                if (dist < bestDistance) {
                     bestDistance = dist;
                     bestIndex = i;
                 }
-                if (i == indices->data.S64[j]) {
+                if (j < indices->n && i == indices->data.S64[j]) {
                     j++;
                     if (dist > DEG_TO_RAD(SPHERICAL_DISTANCE)) {
-                        diag("%ld is in the list, but shouldn't be (%lf vs %lf)",
-                             i, RAD_TO_DEG(dist), SPHERICAL_DISTANCE);
+                        diag("%ld (%lf,%lf) is in the list, but shouldn't be (%lf vs %lf)",
+                             i, ra->data.F64[i], dec->data.F64[i], RAD_TO_DEG(dist), SPHERICAL_DISTANCE);
                         allgood = false;
+                        bad++;
+                    } else {
+#if 0
+                        diag("%ld (%lf,%lf) correctly identified in the list (%lf vs %lf)",
+                             i, ra->data.F64[i], dec->data.F64[i], RAD_TO_DEG(dist), SPHERICAL_DISTANCE);
+#endif
                     }
                 } else if (dist <= DEG_TO_RAD(SPHERICAL_DISTANCE)) {
-                    diag("%ld is not in the list, but should be (%lf vs %lf)",
-                         i, RAD_TO_DEG(dist), SPHERICAL_DISTANCE);
+                    diag("%ld (%lf,%lf) is not in the list, but should be (%lf vs %lf)",
+                         i, ra->data.F64[i], dec->data.F64[i], RAD_TO_DEG(dist), SPHERICAL_DISTANCE);
                     allgood = false;
+                    bad++;
                 }
             }
-            ok(allgood, "list is acurate");
+            ok(allgood, "list is accurate: %ld bad", bad);
             ok(bestIndex == closeIndex, "correct point: %ld vs %ld", closeIndex, bestIndex);