Changeset 27838 for branches/tap_branches/Ohana/src/addstar/src/find_matches_refstars.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

• . (modified) (1 prop)
• Ohana (modified) (1 prop)
• Ohana/src/addstar/src/find_matches_refstars.c (modified) (17 diffs)

branches/tap_branches/Ohana/src/addstar/src/find_matches_refstars.c

@@ -1 +1 @@
 # include "addstar.h"
 
-int find_matches_refstars (SkyRegion *region, Stars **stars, int Nstars, Catalog *catalog, AddstarClientOptions options) {
-
-  int i, j, k, n, m, N, J;
+int find_matches_refstars (SkyRegion *region, Stars **stars, unsigned int Nstars, Catalog *catalog, AddstarClientOptions options) {
+
+  off_t i, j, n, N, J;
   double RADIUS, RADIUS2;
   double *X1, *Y1, *X2, *Y2;
   double dX, dY, dR;
-  int *N1, *N2,  *next, *next_miss, last, last_miss;
-  int Nave, NAVE, Nmeas, NMEAS, Nmiss, NMISS, Nmatch;
+  off_t *N1, *N2,  *next_meas;
+  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch;
   unsigned int objID, catID;
-  Measure *tmpmeasure;
-  Missing *tmpmissing;
   Coords tcoords;
   int Nsecfilt;
@@ -26 +24 @@
   ALLOCATE (X1, double, Nstars);
   ALLOCATE (Y1, double, Nstars);
-  ALLOCATE (N1, int,   Nstars);
+  ALLOCATE (N1, off_t,  Nstars);
 
   /** allocate local arrays (catalog) **/
@@ -33 +31 @@
   ALLOCATE (X2, double, NAVE);
   ALLOCATE (Y2, double, NAVE);
-  ALLOCATE (N2, int, NAVE);
-  ALLOCATE (catalog[0].found, int, NAVE);
+  ALLOCATE (N2, off_t,  NAVE);
+  ALLOCATE (catalog[0].found, off_t, NAVE);
   REALLOCATE (catalog[0].average, Average, NAVE);
   REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt);
@@ -40 +38 @@
 
   Nmatch = 0;
-  Nmeas = catalog[0].Nmeasure;
-  NMEAS = Nmeas + 1000;
-  ALLOCATE (next, int, NMEAS);
-  REALLOCATE (catalog[0].measure, Measure, NMEAS);
-  
-  Nmiss = catalog[0].Nmissing;
-  NMISS = Nmiss + 1000;
-  if ((NMISS < 1) || (NMISS > 1e10)) {
-      fprintf (stderr, "weird value for NMISS: %d\n", NMISS);
-  }
-  ALLOCATE (next_miss, int, NMISS);
-  REALLOCATE (catalog[0].missing, Missing, NMISS);
+  NMEAS = Nmeas = catalog[0].Nmeasure;
   
   // current max obj ID for this catalog
@@ -82 +69 @@
 
   /* set up pointers for linked list of measurements */
-  for (i = 0; i < Nmeas - 1; i++) {
-    next[i] = i+1;
-  }
-  next[i] = -1;
-  last = i;
-
-  for (i = 0; i < Nmiss - 1; i++) {
-    next_miss[i] = i+1;
-  }
-  next_miss[i] = -1;
-  last_miss = i;
+  if (catalog[0].sorted && (catalog[0].Nmeasure == catalog[0].Nmeas_disk)) {
+    // this version is only valid if we have done a full catalog load, and if the catalog
+    // is sorted while processed
+    next_meas = init_measure_links (catalog[0].average, Nave, catalog[0].measure, Nmeas);
+  } else {
+    next_meas = build_measure_links (catalog[0].average, Nave, catalog[0].measure, Nmeas);
+  }    
 
   /* choose a radius for matches */
@@ -105 +88 @@
   // XXX could use NREFSTAR_GROUP to do this match more quicky
   for (i = j = 0; (i < Nstars) && (j < Nave); ) {
+    if (!finite(X1[i]) || !finite(Y1[i])) { 
+      i++; 
+      continue;
+    }
+    if (!finite(X2[j]) || !finite(Y2[j])) { 
+      j++; 
+      continue;
+    }
     
     dX = X1[i] - X2[j];
-    if (dX <= -2*RADIUS) {
+    if (dX <= -1.02*RADIUS) {
       i++;
       continue;
     }
-    if (dX >= 2*RADIUS) {
+    if (dX >= 1.02*RADIUS) {
       j++;
       continue;
@@ -117 +108 @@
 
     /* negative dX: j is too large, positive dX, i is too large */
-    for (J = j; (dX > -2*RADIUS) && (J < Nave); J++) {
+    for (J = j; (dX > -1.02*RADIUS) && (J < Nave); J++) {
       dX = X1[i] - X2[J];
       dY = Y1[i] - Y2[J];
@@ -123 +114 @@
       if (dR > RADIUS2) continue;
 
+      /*** a match is found, add to average, measure ***/
       Nmatch ++;
       n = N2[J];
       N = N1[i];
-      m = catalog[0].average[n].measureOffset;  
 
       /** in replace mode, search for entry and replace values M, dM, R, D */
-      if (options.replace && replace_match (&catalog[0].average[n], &catalog[0].measure[m], stars[N])) continue;
-
-      /** insert star in measurement list */
-      /* find last measurement of this star */
-      for (k = 0; k < catalog[0].average[n].Nmeasure - 1; k++) m = next[m];
-
-      /* set up references */
-      next[Nmeas] = next[m];
-      next[m] = Nmeas;
-
-      /* last just was moved */
-      if (next[Nmeas] == -1) last = Nmeas;
+      if (options.replace && replace_match (&catalog[0].average[n], catalog[0].measure, stars[N])) continue;
+
+      if (Nmeas == NMEAS) {
+        NMEAS = Nmeas + 1000;
+        REALLOCATE (next_meas, off_t, NMEAS);
+        REALLOCATE (catalog[0].measure, Measure, NMEAS);
+      }
+
+      /* add to end of measurement list */
+      add_meas_link (&catalog[0].average[n], next_meas, Nmeas, NMEAS);
         
       /** add measurements for this star **/
@@ -150 +139 @@
       catalog[0].measure[Nmeas].dR       = 3600.0*(catalog[0].average[n].R - stars[N][0].average.R);
       catalog[0].measure[Nmeas].dD       = 3600.0*(catalog[0].average[n].D - stars[N][0].average.D);
+      catalog[0].measure[Nmeas].dbFlags  = 0;
+      catalog[0].measure[Nmeas].averef   = n;
+      catalog[0].measure[Nmeas].objID    = catalog[0].average[n].objID;
+      catalog[0].measure[Nmeas].catID    = catalog[0].catID;
 
       // rationalize dR:
@@ -162 +155 @@
           catalog[0].measure[Nmeas].dR = 3600.0*(catalog[0].average[n].R - stars[N][0].average.R);
       }
-
-      catalog[0].measure[Nmeas].averef   = n;
-      catalog[0].measure[Nmeas].dbFlags  = 0;
-      catalog[0].measure[Nmeas].objID    = catalog[0].average[n].objID;
-      catalog[0].measure[Nmeas].catID    = catalog[0].catID;
 
       catalog[0].measure[Nmeas].t        = (TIMEREF == 0) ? stars[N][0].measure.t      : TIMEREF; /** careful : time_t vs e_time **/
@@ -204 +192 @@
       catalog[0].average[n].Nmeasure ++;
       Nmeas ++;
-      if (Nmeas == NMEAS) {
-        NMEAS = Nmeas + 1000;
-        REALLOCATE (next, int, NMEAS);
-        REALLOCATE (catalog[0].measure, Measure, NMEAS);
-      }
-
-      update_coords (&catalog[0].average[n], &catalog[0].measure[0], next);
+
+      // update_coords (&catalog[0].average[n], &catalog[0].measure[0], next_meas);
     }
     i++;
@@ -224 +207 @@
 
   for (i = 0; (i < Nstars) && !options.only_match; i+=NREFSTAR_GROUP) {
+    if (Nmeas >= NMEAS - NREFSTAR_GROUP) {
+      NMEAS = Nmeas + 1000;
+      REALLOCATE (next_meas, off_t, NMEAS);
+      REALLOCATE (catalog[0].measure, Measure, NMEAS);
+    }
+    if (Nave == NAVE) {
+      NAVE = Nave + 1000;
+      REALLOCATE (catalog[0].average, Average, NAVE);
+      REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt);
+    }
+
     N = N1[i];
     if (stars[N][0].found > -1) continue;
+    if (!IN_REGION (stars[N][0].average.R, stars[N][0].average.D)) continue;
 
     catalog[0].average[Nave].R             = stars[N][0].average.R;
@@ -260 +255 @@
 
     catalog[0].average[Nave].Xp            = 0;
-    catalog[0].average[Nave].ChiSq         = 0.0;
+    catalog[0].average[Nave].ChiSqAve      = 0.0;
+    catalog[0].average[Nave].ChiSqPM       = 0.0;
+    catalog[0].average[Nave].ChiSqPar      = 0.0;
+    catalog[0].average[Nave].Tmean         = 0;
+    catalog[0].average[Nave].Trange        = 0;
     catalog[0].average[Nave].Npos          = 0;
 
@@ -266 +265 @@
     catalog[0].average[Nave].catID         = catID;
     catalog[0].average[Nave].flags         = 0;
+    if (PSPS_ID) {
+        catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R,
+                                                            catalog[0].average[Nave].D);
+    } else {
+        catalog[0].average[Nave].extID         = 0;
+    }
+
 
     objID ++;
@@ -295 +301 @@
 
       stars[N][0].found = Nmeas;
-      next[last] = Nmeas;
-      next[Nmeas] = -1;
-      last = Nmeas;
-      Nmeas ++;
-      if (Nmeas == NMEAS) {
-        NMEAS = Nmeas + 1000;
-        REALLOCATE (next, int, NMEAS);
-        REALLOCATE (catalog[0].measure, Measure, NMEAS);
-      }
-    }
-
+      next_meas[Nmeas] = -1;  // initial value here, update below
+      Nmeas ++; 
+    }
+
+    for (j = 0; j < NREFSTAR_GROUP - 1; j++) {
+        next_meas[Nmeas - NREFSTAR_GROUP + j] = Nmeas - NREFSTAR_GROUP + j + 1;
+    }
     Nave ++;
-    if (Nave == NAVE) {
-      NAVE = Nave + 1000;
-      REALLOCATE (catalog[0].average, Average, NAVE);
-      REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt);
-    }
   }
       
-  free (catalog[0].found);
   REALLOCATE (catalog[0].average, Average, Nave);
 
-  /* fix order of Measure (memory intensive, but fast) */
-  N = 0; 
-  ALLOCATE (tmpmeasure, Measure, Nmeas);
-  for (i = 0; i < Nave; i++) {
-    n = catalog[0].average[i].measureOffset;
-    catalog[0].average[i].measureOffset = N;
-    for (k = 0; k < catalog[0].average[i].Nmeasure; k++, N++) {
-      tmpmeasure[N] = catalog[0].measure[n]; 
-      tmpmeasure[N].averef = i;
-      n = next[n];
-    }
-  }
-  free (catalog[0].measure);
-  catalog[0].measure = tmpmeasure;
-    
-  /* fix order of Missing (memory intensive, but fast) */
-  N = 0; 
-  ALLOCATE (tmpmissing, Missing, Nmiss);
-  for (i = 0; i < Nave; i++) {
-    if (catalog[0].average[i].Nmissing > 0) {
-      n = catalog[0].average[i].missingOffset;
-      catalog[0].average[i].missingOffset = N;
-      for (k = 0; k < catalog[0].average[i].Nmissing; k++, N++) {
-        tmpmissing[N] = catalog[0].missing[n]; 
-        n = next_miss[n];
-      }
-    }
-  }
-  free (catalog[0].missing);
-  catalog[0].missing = tmpmissing;
+  // XXX allow for unsorted output?
+  catalog[0].measure = sort_measure (catalog[0].average, Nave, catalog[0].measure, Nmeas, next_meas);
 
   /* note stars which have been found in this catalog */
@@ -357 +325 @@
   }
 
+  catalog[0].objID    = objID; // new max value, save on catalog close
+  catalog[0].Naverage = Nave;
+  catalog[0].Nmeasure = Nmeas;
+  catalog[0].Nsecf_mem = Nave*Nsecfilt;
+  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: %lld %lld %lld, (%lld matches)\n", (long long) Nstars, (long long) Nave, (long long) Nmeas, (long long) Nmatch);
+
+  free (catalog[0].found);
   free (X1);
   free (Y1);
@@ -363 +338 @@
   free (Y2);
   free (N2);
-  free (next);
-  free (next_miss);
-
-  catalog[0].objID    = objID; // new max value, save on catalog close
-  catalog[0].Naverage = Nave;
-  catalog[0].Nmeasure = Nmeas;
-  catalog[0].Nmissing = Nmiss;
-  catalog[0].Nsecf_mem = Nave*Nsecfilt;
-  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nmiss: %d %d %d %d, (%d matches)\n", Nstars, Nave, Nmeas, Nmiss, Nmatch);
+  free (next_meas);
+
   return (Nmatch);
 }