Changeset 37403 for branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge

Timestamp:

Sep 19, 2014, 4:05:27 PM (12 years ago)

Author:

eugene

Message:

merge changes from trunk

Location:

branches/eam_branches/ps2-tc3-20130727

Files:

• . (modified) (1 prop)
• Ohana (modified) (1 prop)
• Ohana/src/dvomerge/include/dvomerge.h (modified) (1 diff)
• Ohana/src/dvomerge/include/dvoverify.h (modified) (1 diff)
• Ohana/src/dvomerge/src/IDmapIO.c (modified) (1 diff)
• Ohana/src/dvomerge/src/LoadCatalog.c (modified) (1 diff)
• Ohana/src/dvomerge/src/build_links.c (modified) (13 diffs)
• Ohana/src/dvomerge/src/dvo_image_merge_dbs.c (modified) (2 diffs)
• Ohana/src/dvomerge/src/dvomergeContinue.c (modified) (1 diff)
• Ohana/src/dvomerge/src/dvomergeContinue_threaded.c (modified) (1 diff)
• Ohana/src/dvomerge/src/dvomergeCreate.c (modified) (1 diff)
• Ohana/src/dvomerge/src/dvomergeFromList.c (modified) (1 diff)
• Ohana/src/dvomerge/src/dvomergeImageIDs.c (modified) (2 diffs)
• Ohana/src/dvomerge/src/dvomergeUpdate.c (modified) (2 diffs)
• Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c (modified) (2 diffs)
• Ohana/src/dvomerge/src/dvomergeUpdate_threaded.c (modified) (5 diffs)
• Ohana/src/dvomerge/src/dvorepair.c (modified) (1 diff)
• Ohana/src/dvomerge/src/dvorepairCPT.c (modified) (2 diffs)
• Ohana/src/dvomerge/src/dvorepairDeleteImageList.c (modified) (3 diffs)
• Ohana/src/dvomerge/src/dvorepairFixCPT.c (modified) (2 diffs)
• Ohana/src/dvomerge/src/dvorepairFixTables.c (modified) (3 diffs)
• Ohana/src/dvomerge/src/dvorepairImageVsMeasure.c (modified) (1 diff)
• Ohana/src/dvomerge/src/dvorepairImagesVsMeasures.c (modified) (1 diff)
• Ohana/src/dvomerge/src/dvoverify_args.c (modified) (2 diffs)
• Ohana/src/dvomerge/src/dvoverify_catalogs.c (modified) (2 diffs)
• Ohana/src/dvomerge/src/dvoverify_client.c (modified) (1 diff)
• Ohana/src/dvomerge/src/dvoverify_utils.c (modified) (7 diffs)
• Ohana/src/dvomerge/src/merge_catalogs_new.c (modified) (4 diffs)
• Ohana/src/dvomerge/src/merge_catalogs_old.c (modified) (23 diffs)
• Ohana/src/dvomerge/src/replace_match.c (modified) (3 diffs)

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/include/dvomerge.h

@@ -129 +129 @@
 int        merge_catalogs_old     PROTO((SkyRegion *region, Catalog *output, Catalog *input, double RADIUS, int *secflitMap));
 
+off_t     *build_measure_links    PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure));
 off_t     *init_measure_links     PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure));
+int        add_meas_link          PROTO((Average *average, off_t *next, off_t Nmeasure, off_t NMEASURE));
+Measure   *sort_measure           PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure, off_t *next));
+
+off_t     *build_lensing_links    PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing));
+off_t     *init_lensing_links     PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing));
+int        add_lens_link          PROTO((Average *average, off_t *next, off_t Nlensing, off_t NLENSING));
+Lensing   *sort_lensing           PROTO((Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing, off_t *next));
+
 off_t     *init_missing_links     PROTO((Average *average, off_t Naverage, Missing *missing, off_t Nmissing));
-off_t     *build_measure_links    PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure));
-int        add_meas_link          PROTO((Average *average, off_t *next, off_t Nmeasure, off_t NMEASURE));
 int        add_miss_link          PROTO((Average *average, off_t *next, off_t Nmissing));
-Measure   *sort_measure           PROTO((Average *average, off_t Naverage, Measure *measure, off_t Nmeasure, off_t *next));
 Missing   *sort_missing           PROTO((Average *average, off_t Naverage, Missing *missing, off_t Nmissing, off_t *next_miss));
 

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/include/dvoverify.h

@@ -56 +56 @@
 int LoadImageIDs (char *catdir);
 int CheckImageID (Catalog *catalog);
+
+int SaveImageIDsSmall(char *filename);
+int LoadImageIDsSmall (char *filename);

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/IDmapIO.c

@@ -27 +27 @@
   gfits_create_table_header (&theader, "BINTABLE", "IMAGE_ID_MAP");
 
-  gfits_define_bintable_column (&theader, "J", "OLD_IDS", "old image IDs", NULL, 1.0, 1.0*0x8000);
-  gfits_define_bintable_column (&theader, "J", "NEW_IDS", "new image IDs", NULL, 1.0, 1.0*0x8000);
+  gfits_define_bintable_column (&theader, "J", "OLD_IDS", "old image IDs", NULL, 1.0, FT_BZERO_INT32);
+  gfits_define_bintable_column (&theader, "J", "NEW_IDS", "new image IDs", NULL, 1.0, FT_BZERO_INT32);
 
   // generate the output array that carries the data

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/LoadCatalog.c

@@ -8 +8 @@
 
     // always load all of the data (if any exists)
-    catalog[0].catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
+    catalog[0].catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF | LOAD_LENSING | LOAD_LENSOBJ;
 
     catalog[0].catformat = dvo_catalog_catformat (CATFORMAT);  // set the default catformat from config data

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/build_links.c

@@ -1 +1 @@
 # include "dvomerge.h"
+
+/* 
+
+There are two modes for the measure table: sorted and unsorted. 
+
+In sorted mode, all measures associated with a given average are in a single block.  The block
+is pointed to by average->measureOffset and the range is average->Nmeasure
+
+In unsorted mode, it is not possible to go directly from average to measure without scanning.
+In this case, it is necessary to use the value measure->averef to find the corresponding
+average entry.  
+
+Note that average->measureOffset and measure->averef are only valid for a given load of the
+data: they refer to the sequence number in the data blocks.
+
+next_meas is a list of the equivalent sequence of the measure block as if it were sorted.
+
+to find the sequence of measurements for a given average:
+n_0 = average->measureOffset
+n_1 = next_meas[n_0]
+n_i = next_meas[n_i-1]
+
+*/
 
 /* build the initial links assuming the table is sorted, 
@@ -6 +29 @@
 
   off_t i, j, N;
-  off_t *next;
+  off_t *next_meas;
 
   N = 0;
 
-  ALLOCATE (next, off_t, Nmeasure);
+  ALLOCATE (next_meas, off_t, Nmeasure);
   for (i = 0; i < Naverage; i++, N++) {
     for (j = 0; j < average[i].Nmeasure - 1; j++, N++) {
-      next[N] = N + 1;
+      next_meas[N] = N + 1;
       if (N >= Nmeasure) {
         fprintf (stderr, "WARNING: N out of bounds (1)\n");
       }
     }
-    next[N] = -1;
+    next_meas[N] = -1;
     if (N >= Nmeasure) {
       fprintf (stderr, "WARNING: N out of bounds (2)\n");
@@ -28 +51 @@
     }
   }
-  return (next);
-}
-
-/* average[].measureOffset, average[].Nmeasure are valid within an addstar run */
-int add_meas_link (Average *average, off_t *next, off_t Nmeasure, off_t NMEASURE) {
-
-  off_t k, m;
-
-  /* if we have trouble, check validity of next[m] : m < Nmeasure */
-  m = average[0].measureOffset;  
-
-  for (k = 0; k < average[0].Nmeasure - 1; k++)  {
-    m = next[m];
-    if (m >= NMEASURE) {
-      fprintf (stderr, "WARNING: m out of bounds (3)\n");
-    }
-  }
-
-  /* set up references */
-  next[Nmeasure] = -1;
-  if (Nmeasure >= NMEASURE) {
-    fprintf (stderr, "WARNING: Nmeasure out of bounds (1)\n");
-  }
-
-  if (m == -1) {
-    average[0].measureOffset = Nmeasure;
-  } else {
-    next[m] = Nmeasure;
-    if (m >= NMEASURE) {
-      fprintf (stderr, "WARNING: m out of bounds (4)\n");
-    }
-  }
-
-  return (TRUE);
-}
-
-Measure *sort_measure (Average *average, off_t Naverage, Measure *measure, off_t Nmeasure, off_t *next) {
-
-  off_t i, k, n, N;
-  Measure *tmpmeasure;
-
-  /* fix order of Measure (memory intensive, but fast) */
-  N = 0; 
-  ALLOCATE (tmpmeasure, Measure, Nmeasure);
-  for (i = 0; i < Naverage; i++) {
-    n = average[i].measureOffset;
-    average[i].measureOffset = N;
-    for (k = 0; k < average[i].Nmeasure; k++, N++) {
-      tmpmeasure[N] = measure[n]; 
-      tmpmeasure[N].averef = i;
-      n = next[n];
-    }
-  }
-  free (measure);
-  return (tmpmeasure);
+  return (next_meas);
 }
 
@@ -95 +64 @@
 
   off_t i, m, k, Nm, averef;
-  off_t *next;
-
-  ALLOCATE (next, off_t, Nmeasure);
+  off_t *next_meas;
+
+  ALLOCATE (next_meas, off_t, Nmeasure);
 
   /* reset the Nm, offset values for average */
@@ -108 +77 @@
     averef = measure[Nm].averef;
     m = average[averef].measureOffset;  
-    next[Nm] = -1;
+    next_meas[Nm] = -1;
 
     if (m == -1) { /* no links yet for source */
@@ -116 +85 @@
     }
 
-    for (k = 0; next[m] != -1; k++) {
-      m = next[m];
+    for (k = 0; next_meas[m] != -1; k++) {
+      m = next_meas[m];
       if (m >= Nmeasure) {
         fprintf (stderr, "WARNING: m out of bounds (1)\n");
@@ -124 +93 @@
 
     average[averef].Nmeasure = k + 2;
-    next[m] = Nm;
+    next_meas[m] = Nm;
     if (m >= Nmeasure) {
       fprintf (stderr, "WARNING: m out of bounds (2)\n");
     }
   }
-  return (next);
-}
+  return (next_meas);
+}
+
+/* average[].measureOffset, average[].Nmeasure are valid within an addstar run */
+int add_meas_link (Average *average, off_t *next_meas, off_t Nmeasure, off_t NMEASURE) {
+
+  off_t k, m;
+
+  /* if we have trouble, check validity of next_meas[m] : m < Nmeasure */
+  m = average[0].measureOffset;  
+
+  for (k = 0; k < average[0].Nmeasure - 1; k++)  {
+    m = next_meas[m];
+    if (m >= NMEASURE) {
+      fprintf (stderr, "WARNING: m out of bounds (3)\n");
+    }
+  }
+
+  /* set up references */
+  next_meas[Nmeasure] = -1;
+  if (Nmeasure >= NMEASURE) {
+    fprintf (stderr, "WARNING: Nmeasure out of bounds (1)\n");
+  }
+
+  if (m == -1) {
+    average[0].measureOffset = Nmeasure;
+  } else {
+    next_meas[m] = Nmeasure;
+    if (m >= NMEASURE) {
+      fprintf (stderr, "WARNING: m out of bounds (4)\n");
+    }
+  }
+
+  return (TRUE);
+}
+
+Measure *sort_measure (Average *average, off_t Naverage, Measure *measure, off_t Nmeasure, off_t *next_meas) {
+
+  off_t i, k, n, N;
+  Measure *tmpmeasure;
+
+  /* fix order of Measure (memory intensive, but fast) */
+  N = 0; 
+  ALLOCATE (tmpmeasure, Measure, Nmeasure);
+  for (i = 0; i < Naverage; i++) {
+    n = average[i].measureOffset;
+    average[i].measureOffset = N;
+    for (k = 0; k < average[i].Nmeasure; k++, N++) {
+      if (n == -1) abort();
+      tmpmeasure[N] = measure[n]; 
+      if (measure[n].averef != i) abort();
+      tmpmeasure[N].averef = i;
+      n = next_meas[n];
+    }
+  }
+  free (measure);
+  return (tmpmeasure);
+}
+
+/*******************************************************************************************/
 
 /* build the initial links assuming the table is sorted */
@@ -136 +163 @@
 
   off_t i, j, N;
-  off_t *next;
+  off_t *next_miss;
 
   N = 0;
 
-  ALLOCATE (next, off_t, Nmissing);
+  ALLOCATE (next_miss, off_t, Nmissing);
   for (i = 0; i < Naverage; i++) {
     for (j = 0; j < average[i].Nmissing - 1; j++, N++) {
-      next[N] = N + 1;
+      next_miss[N] = N + 1;
     }
     if (average[i].Nmissing > 0) {
-      next[N] = -1;
+      next_miss[N] = -1;
       if (N >= Nmissing) {
         fprintf (stderr, "overflow in init_missing_links");
@@ -155 +182 @@
 
   }
-  return (next);
-}
-
-int add_miss_link (Average *average, off_t *next, off_t Nmissing) {
+  return (next_miss);
+}
+
+int add_miss_link (Average *average, off_t *next_miss, off_t Nmissing) {
 
   off_t k, m;
@@ -165 +192 @@
   if (average[0].Nmissing < 1) {
     average[0].missingOffset = Nmissing;
-    next[Nmissing] = -1;
+    next_miss[Nmissing] = -1;
     return (TRUE);
   }
 
   m = average[0].missingOffset;  
-  for (k = 0; k < average[0].Nmissing - 1; k++) m = next[m];
+  for (k = 0; k < average[0].Nmissing - 1; k++) m = next_miss[m];
   /* set up references */
-  next[Nmissing] = -1;
-  next[m] = Nmissing;
+  next_miss[Nmissing] = -1;
+  next_miss[m] = Nmissing;
   return (TRUE);
 }
@@ -180 +207 @@
    we must always save the missing table, if it exists */
 
-Missing *sort_missing (Average *average, off_t Naverage, Missing *missing, off_t Nmissing, off_t *next) {
+Missing *sort_missing (Average *average, off_t Naverage, Missing *missing, off_t Nmissing, off_t *next_miss) {
 
   off_t i, k, n, N;
@@ -193 +220 @@
     for (k = 0; k < average[i].Nmissing; k++, N++) {
       tmpmissing[N] = missing[n]; 
-      n = next[n];
+      n = next_miss[n];
     }
   }
@@ -199 +226 @@
   return (tmpmissing);
 }
+
+/*******************************************************************************************/
+
+/* build the initial links assuming the table is sorted, 
+   not partial, and has a correct set of average[].lensingOffset,Nlensing values */
+off_t *init_lensing_links (Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing) {
+
+  off_t i, j, N;
+  off_t *next_lens;
+
+  N = 0;
+
+  ALLOCATE (next_lens, off_t, Nlensing);
+  for (i = 0; i < Naverage; i++) {
+    if (!average[i].Nlensing) continue;
+    for (j = 0; j < average[i].Nlensing - 1; j++, N++) {
+      next_lens[N] = N + 1;
+      if (N >= Nlensing) {
+        fprintf (stderr, "WARNING: N out of bounds (1)\n");
+      }
+    }
+    next_lens[N] = -1;
+    if (N >= Nlensing) {
+      fprintf (stderr, "WARNING: N out of bounds (2)\n");
+    }
+
+    if (N >= Nlensing) {
+      fprintf (stderr, "overflow in init_lensing_links\n");
+      abort ();
+    }
+    N++;
+  }
+  return (next_lens);
+}
+
+/* construct lensing links which are valid FOR THIS LOAD
+ * - if we have a full load, we will get links which can
+ *   be used by other programs (eg, relphot, etc)
+ * - if we have a partial load, the links are only valid
+ *   for that partial load
+ */ 
+
+off_t *build_lensing_links (Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing) {
+
+  off_t i, m, k, Nm, averef;
+  off_t *next_lens;
+
+  ALLOCATE (next_lens, off_t, Nlensing);
+
+  /* reset the Nm, offset values for average */
+  for (i = 0; i < Naverage; i++) {
+    average[i].lensingOffset = -1;
+    average[i].Nlensing     =  0;
+  }
+
+  for (Nm = 0; Nm < Nlensing; Nm++) {
+    averef = lensing[Nm].averef;
+    m = average[averef].lensingOffset;  
+    next_lens[Nm] = -1;
+
+    if (m == -1) { /* no links yet for source */
+      average[averef].lensingOffset = Nm;
+      average[averef].Nlensing     = 1;
+      continue;
+    }
+
+    for (k = 0; next_lens[m] != -1; k++) {
+      m = next_lens[m];
+      if (m >= Nlensing) {
+        fprintf (stderr, "WARNING: m out of bounds (1)\n");
+      }
+    }
+
+    average[averef].Nlensing = k + 2;
+    next_lens[m] = Nm;
+    if (m >= Nlensing) {
+      fprintf (stderr, "WARNING: m out of bounds (2)\n");
+    }
+  }
+  return (next_lens);
+}
+
+/* average[].lensingOffset, average[].Nlensing are valid within an addstar run */
+int add_lens_link (Average *average, off_t *next_lens, off_t Nlensing, off_t NLENSING) {
+
+  off_t k, m;
+
+  /* if we have trouble, check validity of next_lens[m] : m < Nlensing */
+  m = average[0].lensingOffset;  
+
+  for (k = 0; k < average[0].Nlensing - 1; k++)  {
+    m = next_lens[m];
+    if (m >= NLENSING) {
+      fprintf (stderr, "WARNING: m out of bounds (3)\n");
+    }
+  }
+
+  /* set up references */
+  next_lens[Nlensing] = -1;
+  if (Nlensing >= NLENSING) {
+    fprintf (stderr, "WARNING: Nlensing out of bounds (1)\n");
+  }
+
+  if (m == -1) {
+    average[0].lensingOffset = Nlensing;
+  } else {
+    next_lens[m] = Nlensing;
+    if (m >= NLENSING) {
+      fprintf (stderr, "WARNING: m out of bounds (4)\n");
+    }
+  }
+
+  return (TRUE);
+}
+
+Lensing *sort_lensing (Average *average, off_t Naverage, Lensing *lensing, off_t Nlensing, off_t *next_lens) {
+
+  off_t i, k, n, np, N;
+  Lensing *tmplensing;
+
+  /*
+  for (i = 0; i < Naverage; i++) {
+    if (average[i].Nlensing != 4) {
+      fprintf (stderr, "check %d %d %d\n", (int) i, (int) average[i].Nlensing, (int) average[i].lensingOffset);
+    }
+  }
+  */
+
+  /* fix order of Lensing (memory intensive, but fast) */
+  N = 0; 
+  ALLOCATE (tmplensing, Lensing, Nlensing);
+  for (i = 0; i < Naverage; i++) {
+    if (!average[i].Nlensing) continue;
+    n = average[i].lensingOffset;
+    average[i].lensingOffset = N;
+    for (k = 0; k < average[i].Nlensing; k++, N++) {
+      if (n == -1) {
+        fprintf (stderr, "entry after %d has a problem\n", (int) np);
+        abort();
+      }
+      tmplensing[N] = lensing[n]; 
+      if (lensing[n].averef != i) abort();
+      tmplensing[N].averef = i;
+      np = n;
+      n = next_lens[n];
+    }
+  }
+  free (lensing);
+  return (tmplensing);
+}
+
+

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c

@@ -172 +172 @@
   }
 
+  if (!IDmap->Nmap) {
+    fprintf (stderr, "input database has image IDs, but no Image table\n");
+    return FALSE;
+  }
+
   off_t lastID = IDmap->old[IDmap->Nmap-1];
 
+  // update measure.imageID
   for (i = 0; i < catalog[0].Nmeasure; i++) {
     oldID = catalog[0].measure[i].imageID;
@@ -195 +201 @@
     catalog[0].measure[i].imageID = newID;
   }
+
+  // also update lensing.imageID if lensing exists
+  for (i = 0; i < catalog[0].Nlensing; i++) {
+    oldID = catalog[0].lensing[i].imageID;
+    if (oldID == 0) continue;
+
+    newID = dvo_map_image_ID (IDmap, oldID);
+    if (newID == 0) {
+      if (oldID > lastID) {
+        fprintf (stderr, "problem with image IDs : input out of range\n");
+        fprintf (stderr, "old ID: "OFF_T_FMT", last ID: "OFF_T_FMT"\n", oldID, lastID);
+        exit (2);
+      }
+      if (!IDmap->notFound[oldID]) {
+        fprintf (stderr, "cannot find image ID "OFF_T_FMT"\n",  oldID);
+        IDmap->notFound[oldID] = TRUE;
+        // once we discover an imageID is not found, record that fact so we do not complain for every detection
+      }
+      // optionally exit here? or wait until end to report an error?
+      // exit (2);
+    }
+    catalog[0].lensing[i].imageID = newID;
+  }
   return TRUE;
 }

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvomergeContinue.c

@@ -131 +131 @@
       merge_catalogs_old (&outsky[0].regions[j], &outcatalog, &incatalog, RADIUS, secfiltMap);
 
-      outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
-
       // if we receive a signal which would cause us to exit, wait until the full catalog is written
       SetProtect (TRUE);

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvomergeContinue_threaded.c

@@ -74 +74 @@
       merge_catalogs_old (&threadData->outsky->regions[j], &outcatalog, &incatalog, RADIUS, threadData->secfiltMap);
 
-      outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
-
       // if we receive a signal which would cause us to exit, wait until the full catalog is written
       SetProtect (TRUE);

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvomergeCreate.c

@@ -178 +178 @@
     SkyListFree (inlist);
 
-    outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
     dvo_catalog_save (&outcatalog, VERBOSE);
     dvo_catalog_unlock (&outcatalog);

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvomergeFromList.c

@@ -143 +143 @@
     merge_catalogs_old (&skyregion, &outcatalog, &incatalog, RADIUS, secfiltMap);
     
-    outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
-
     dmhObjectAdd (outstat[0].history, &outcatalog.header, inStats);
 

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvomergeImageIDs.c

@@ -18 +18 @@
   if (inDB.dbstate == LCK_EMPTY) {
     dvo_image_unlock (&inDB); // unlock input
+    IDmap->old = NULL;
+    IDmap->new = NULL;
     IDmap->Nmap = 0;
     return TRUE;
@@ -99 +101 @@
   if (inDB.dbstate == LCK_EMPTY) {
     dvo_image_unlock (&inDB); // unlock input
+    IDmap->old = NULL;
+    IDmap->new = NULL;
     IDmap->Nmap = 0;
     return TRUE;

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvomergeUpdate.c

@@ -13 +13 @@
   int NsecfiltInput, NsecfiltOutput;
 
-  double dtime;
-  struct timeval start, stop;
-  gettimeofday (&start, NULL);
+  INITTIME;
 
   CONTINUE = FALSE;
@@ -121 +119 @@
   }
 
-  gettimeofday (&stop, NULL);
-  dtime = DTIME (stop, start);
-
   if (!status) {
-    fprintf (stderr, "ERROR: elapsed time %9.4f sec\n", dtime);
+    MARKTIME ("ERROR: elapsed time %9.4f sec\n", dtime);
     exit (3);
   }
 
-  fprintf (stderr, "SUCCESS: elapsed time %9.4f sec\n", dtime);
+  MARKTIME ("SUCCESS: elapsed time %9.4f sec\n", dtime);
   exit (0);
 }

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c

@@ -140 +140 @@
     }
 
+    if (!incatalog.sorted) {
+      fprintf (stderr, "ERROR: input catalog %s is not sorted (and must be for dvomerge)\n", filename_input);
+      exit (1);
+    }
+
     dvo_update_image_IDs (IDmap, &incatalog);
 
@@ -158 +163 @@
       LoadCatalog (&outcatalog, outlist[0].regions[j], outcatalog.filename, "w", NsecfiltOutput);
 
+      // if no catalog already exists, use the input catalog to define the format
+      if (outcatalog.Naves_disk == 0) {
+        outcatalog.catformat = incatalog.catformat;
+      }
+
       merge_catalogs_old (outlist[0].regions[j], &outcatalog, &incatalog, RADIUS, secfiltMap);
-
-      outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
 
       if (outstat[j].missed) {

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvomergeUpdate_threaded.c

@@ -72 +72 @@
       merge_catalogs_old (&threadData->outsky->regions[j], &outcatalog, &incatalog, RADIUS, threadData->secfiltMap);
 
-      outcatalog.catflags = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
-
       // if we receive a signal which would cause us to exit, wait until the full catalog is written
       SetProtect (TRUE);
@@ -115 +113 @@
   int NsecfiltInput, NsecfiltOutput;
 
-  double dtime;
-  struct timeval start, stop;
-  gettimeofday (&start, NULL);
+  INITTIME;
 
   CONTINUE = FALSE;
@@ -134 +130 @@
 
   fprintf (stderr, "WARNING / ERROR : multi-threaded dvomerge does not handle merge tracking yet\n");
+  fprintf (stderr, "WARNING / ERROR : multi-threaded dvomerge does work in a parallel context\n");
   exit (2);
 
@@ -293 +290 @@
   }
 
-  gettimeofday (&stop, NULL);
-  dtime = DTIME (stop, start);
-  fprintf (stderr, "SUCCESS: elapsed time %9.4f sec\n", dtime);
-
+  MARKTIME ("SUCCESS: elapsed time %9.4f sec\n", dtime);
   exit (0);
 
@@ -318 +312 @@
     }
   }
+  MARKTIME ("ERROR: elapsed time %9.4f sec\n", dtime);
   exit (1);
 }

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvorepair.c

@@ -2 +2 @@
 
 int main (int argc, char **argv) {
+
+  fprintf (stderr, "this program needs to be updated to load old format Measure tables (pre PV1_V5) in which dR,dD are saved, not R,D\n");
+  fprintf (stderr, "reminder: relastro can re-construct R,D from X,Y; FtableToMeasure and vice versa could just NAN those values\n");
+  exit (2);
 
   dvorepair_help(argc, argv);

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvorepairCPT.c

@@ -76 +76 @@
   }
 
-  measure = FtableToMeasure (&cpmFtable, &Nmeasure, &catformat);
+  measure = FtableToMeasure (&cpmFtable, NULL, &Nmeasure, &catformat);
   myAssert(measure, "failed to convert ftable to measure data");
 
@@ -126 +126 @@
     average[Nave].measureOffset = -1;
     average[Nave].missingOffset = -1;
-    average[Nave].extendOffset = -1;
+    average[Nave].refColorBlue = NAN;
+    average[Nave].refColorRed = NAN;
 
     average[Nave].objID = measure[i].objID;

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvorepairDeleteImageList.c

@@ -250 +250 @@
         gfits_scan(&cpmHeaderTBL, "NAXIS1", "%d", 1, &NbytesPerRow);
 
-        measure = FtableToMeasure (&cpmFtable, &Nmeasure, &catformat);
+        measure = FtableToMeasure (&cpmFtable, NULL, &Nmeasure, &catformat);
         myAssert(measure, "failed to convert ftable to measure data");
     
@@ -288 +288 @@
 
         // convert internal to external format 
-        if (!MeasureToFtable (&cpmFtable, measureNew, NmeasureNew, catformat)) {
+        if (!MeasureToFtable (&cpmFtable, NULL, measureNew, NmeasureNew, catformat)) {
           myAbort("trouble converting format");
         }
@@ -411 +411 @@
     average[Nave].measureOffset = -1;
     average[Nave].missingOffset = -1;
-    average[Nave].extendOffset = -1;
+    average[Nave].refColorBlue = NAN;
+    average[Nave].refColorRed = NAN;
 
     average[Nave].objID = measure[i].objID;

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvorepairFixCPT.c

@@ -120 +120 @@
   }
 
-  measure = FtableToMeasure (&cpmFtable, &Nmeasure, &catformat);
+  measure = FtableToMeasure (&cpmFtable, NULL, &Nmeasure, &catformat);
   myAssert(measure, "failed to convert ftable to measure data");
 
@@ -127 +127 @@
 
   // convert internal to external format 
-  if (!MeasureToFtable (&cpmFtable, measure, Nmeasure, catformat)) {
+  if (!MeasureToFtable (&cpmFtable, NULL, measure, Nmeasure, catformat)) {
     myAbort("trouble converting format");
   }

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvorepairFixTables.c

@@ -103 +103 @@
 
       // convert internal to external format 
-      if (!MeasureToFtable (&cpmFtable, measureNew, NmeasureNew, catformat)) {
+      if (!MeasureToFtable (&cpmFtable, NULL, measureNew, NmeasureNew, catformat)) {
         myAbort("trouble converting format");
       }
@@ -183 +183 @@
   gfits_scan(&cpmHeaderTBL, "NAXIS2", "%d", 1, &Nrows);
 
-  measure = FtableToMeasure (&cpmFtable, &Nmeasure, &catformat);
+  measure = FtableToMeasure (&cpmFtable, NULL, &Nmeasure, &catformat);
   myAssert(measure, "failed to convert ftable to measure data");
     
@@ -266 +266 @@
     average[Nave].measureOffset = -1;
     average[Nave].missingOffset = -1;
-    average[Nave].extendOffset = -1;
+    average[Nave].refColorBlue = NAN;
+    average[Nave].refColorRef = NAN;
 
     average[Nave].objID = measure[i].objID;

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvorepairImageVsMeasure.c

@@ -128 +128 @@
       gfits_scan(&cpmHeaderTBL, "NAXIS1", "%d", 1, &NbytesPerRow);
 
-      measure = FtableToMeasure (&cpmFtable, &Nmeasure, &catformat);
+      measure = FtableToMeasure (&cpmFtable, NULL, &Nmeasure, &catformat);
       myAssert(measure, "failed to convert ftable to measure data");
     

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvorepairImagesVsMeasures.c

@@ -134 +134 @@
       gfits_scan(&cpmHeaderTBL, "NAXIS1", "%d", 1, &NbytesPerRow);
 
-      measure = FtableToMeasure (&cpmFtable, &Nmeasure, &catformat);
+      measure = FtableToMeasure (&cpmFtable, NULL, &Nmeasure, &catformat);
       myAssert(measure, "failed to convert ftable to measure data");
     

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvoverify_args.c

@@ -56 +56 @@
   }
 
-  /*** dvoverify -parallel (dvoverify_client, actually) ignores this argument
-       dvoverify should save a subset table of just the IDs, and dvoverify_client should load it
-       dvoverify -cpt should have CHECK_IMAGE_ID = F as default ***/
+  /*** dvoverify -cpt should have CHECK_IMAGE_ID = F as default ***/
        
   CHECK_IMAGE_ID = TRUE;
@@ -150 +148 @@
   }
 
+  CHECK_IMAGE_ID = TRUE;
+  if ((N = get_argument (*argc, argv, "-skip-image-ids"))) {
+    CHECK_IMAGE_ID = FALSE;
+    remove_argument (N, argc, argv);
+  }
+
   // restrict to a portion of the sky
   UserPatch.Rmin = 0;

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvoverify_catalogs.c

@@ -47 +47 @@
     }
 
+    sprintf (filename, "%s/%s.cpx", mycatdir, skylist[0].regions[i][0].name);
+    if (!VerifyTableFile (filename)) {
+      Nbad ++;
+      skipIndexCheck = TRUE;
+      AddFailures (filename);
+    }
+
+    sprintf (filename, "%s/%s.cpy", mycatdir, skylist[0].regions[i][0].name);
+    if (!VerifyTableFile (filename)) {
+      Nbad ++;
+      skipIndexCheck = TRUE;
+      AddFailures (filename);
+    }
+
     sprintf (filename, "%s/%s.cpt", mycatdir, skylist[0].regions[i][0].name);
     if (!skipIndexCheck) {
@@ -80 +94 @@
     exit (1);
   }    
+
+  // save the IDlist as a smaller FITS table
+  if (CHECK_IMAGE_ID) {
+    char filename[DVO_MAX_PATH];
+    snprintf (filename, DVO_MAX_PATH, "%s/ImageIDs.fits", CATDIR);
+    if (!SaveImageIDsSmall (filename)) {
+      fprintf (stderr, "ERROR: failure to save image IDs\n");
+      exit (2);
+    }
+  }
 
   int i;

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvoverify_client.c

@@ -22 +22 @@
     // XXX in client mode, we should be reading a reduced table generated by the calling
     // serial program
-    LoadImageIDs (CATDIR);
+    char filename[DVO_MAX_PATH];
+    snprintf (filename, DVO_MAX_PATH, "%s/ImageIDs.fits", CATDIR);
+    if (!LoadImageIDsSmall (filename)) {
+      fprintf (stderr, "ERROR: failure to load image IDs\n");
+      exit (2);
+    }
   }
 

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/dvoverify_utils.c

@@ -178 +178 @@
   catalog.catformat = DVO_FORMAT_UNDEF; // read-only,do not set the catformat
   catalog.catmode   = DVO_MODE_UNDEF; // read-only, do not set the catmode
-  catalog.catflags  = LOAD_AVES | LOAD_MEAS;
+  catalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_LENSING | LOAD_LENSOBJ;
   catalog.Nsecfilt  = 0;
   
@@ -214 +214 @@
   int NmeasureTotal = 0;
   int measureOffsetOK = TRUE;
+  int NlensingTotal = 0;
+  int lensingOffsetOK = TRUE;
+  int NlensobjTotal = 0;
+  int lensobjOffsetOK = TRUE;
   for (i = 0; i < catalog.Naverage; i++) {
     NmeasureTotal += catalog.average[i].Nmeasure;
@@ -227 +231 @@
       fprintf (stderr, "measureOffset + Nmeasure > catalog.Nmeasure : %d %d %d\n", i, catalog.average[i].Nmeasure, (int) catalog.Nmeasure);
     }
+
+    // sum of Nlensing for each object < Nlensing for catalog
+    NlensingTotal += catalog.average[i].Nlensing;
+    if (VERBOSE && !(NlensingTotal <= catalog.Nlensing)) {
+      fprintf (stderr, "NlensingTotal > catalog.Nlensing: %d %d %d\n", i, catalog.average[i].Nlensing, (int) catalog.Nlensing);
+    }
+    // lensingOffset needs to be in range for each object
+    lensingOffsetOK &= (catalog.average[i].lensingOffset < catalog.Nlensing);
+    if (VERBOSE && !(catalog.average[i].lensingOffset < catalog.Nlensing)) {
+      fprintf (stderr, "lensingOffset >= catalog.Nlensing: %d %d %d\n", i, catalog.average[i].lensingOffset, (int) catalog.Nlensing);
+    }
+    lensingOffsetOK &= (catalog.average[i].lensingOffset + catalog.average[i].Nlensing <= catalog.Nlensing);
+    if (VERBOSE && !(catalog.average[i].lensingOffset + catalog.average[i].Nlensing <= catalog.Nlensing)) {
+      fprintf (stderr, "lensingOffset + Nlensing > catalog.Nlensing : %d %d %d\n", i, catalog.average[i].Nlensing, (int) catalog.Nlensing);
+    }
+
+    NlensobjTotal += catalog.average[i].Nlensobj;
+    if (VERBOSE && !(NlensobjTotal <= catalog.Nlensobj)) {
+      fprintf (stderr, "NlensobjTotal > catalog.Nlensobj: %d %d %d\n", i, catalog.average[i].Nlensobj, (int) catalog.Nlensobj);
+    }
+    lensobjOffsetOK &= (catalog.average[i].lensobjOffset < catalog.Nlensobj);
+    if (VERBOSE && !(catalog.average[i].lensobjOffset < catalog.Nlensobj)) {
+      fprintf (stderr, "lensobjOffset >= catalog.Nlensobj: %d %d %d\n", i, catalog.average[i].lensobjOffset, (int) catalog.Nlensobj);
+    }
+    lensobjOffsetOK &= (catalog.average[i].lensobjOffset + catalog.average[i].Nlensobj <= catalog.Nlensobj);
+    if (VERBOSE && !(catalog.average[i].lensobjOffset + catalog.average[i].Nlensobj <= catalog.Nlensobj)) {
+      fprintf (stderr, "lensobjOffset + Nlensobj > catalog.Nlensobj : %d %d %d\n", i, catalog.average[i].Nlensobj, (int) catalog.Nlensobj);
+    }
   }
 
@@ -236 +268 @@
   if (NmeasureTotal != catalog.Nmeasure) {
     fprintf (stderr, "ERROR: catalog %s has an invalid Nmeasure\n", catalog.filename);
+    status = FALSE;
+  }
+
+  if (!lensingOffsetOK) {
+    fprintf (stderr, "ERROR: catalog %s has an invalid lensingOffset\n", catalog.filename);
+    status = FALSE;
+  }
+
+  if (NlensingTotal != catalog.Nlensing) {
+    fprintf (stderr, "ERROR: catalog %s has an invalid Nlensing\n", catalog.filename);
+    status = FALSE;
+  }
+
+  if (!lensobjOffsetOK) {
+    fprintf (stderr, "ERROR: catalog %s has an invalid lensobjOffset\n", catalog.filename);
+    status = FALSE;
+  }
+
+  if (NlensobjTotal != catalog.Nlensobj) {
+    fprintf (stderr, "ERROR: catalog %s has an invalid Nlensobj\n", catalog.filename);
     status = FALSE;
   }
@@ -247 +299 @@
   }
 
-  int objIDsOK = TRUE;
-  int catIDsOK = TRUE;
-  int averefOK = TRUE;
-
-  for (i = 0; i < catalog.Naverage; i++) {
-    m = catalog.average[i].measureOffset;
-    for (j = 0; j < catalog.average[i].Nmeasure; j++) {
-      objIDsOK &= (catalog.average[i].objID == catalog.measure[m+j].objID);
-      catIDsOK &= (catalog.average[i].catID == catalog.measure[m+j].catID);
-      averefOK &= (catalog.measure[m+j].averef == i);
-    }
-  }
+  // check measure <-> average links
+  {
+    int objIDsOK = TRUE;
+    int catIDsOK = TRUE;
+    int averefOK = TRUE;
+
+    for (i = 0; i < catalog.Naverage; i++) {
+      m = catalog.average[i].measureOffset;
+      for (j = 0; j < catalog.average[i].Nmeasure; j++) {
+        objIDsOK &= (catalog.average[i].objID == catalog.measure[m+j].objID);
+        catIDsOK &= (catalog.average[i].catID == catalog.measure[m+j].catID);
+        averefOK &= (catalog.measure[m+j].averef == i);
+      }
+    }
     
-  if (!objIDsOK) {
-    fprintf (stderr, "ERROR: catalog %s has invalid obj IDs\n", catalog.filename);
-    status = FALSE;
-  }
-  if (!catIDsOK) {
-    fprintf (stderr, "ERROR: catalog %s has invalid cat IDs\n", catalog.filename);
-    status = FALSE;
-  }
-  if (!averefOK) {
-    fprintf (stderr, "ERROR: catalog %s has invalid averef values\n", catalog.filename);
-    status = FALSE;
+    if (!objIDsOK) {
+      fprintf (stderr, "ERROR: catalog %s has invalid obj IDs\n", catalog.filename);
+      status = FALSE;
+    }
+    if (!catIDsOK) {
+      fprintf (stderr, "ERROR: catalog %s has invalid cat IDs\n", catalog.filename);
+      status = FALSE;
+    }
+    if (!averefOK) {
+      fprintf (stderr, "ERROR: catalog %s has invalid averef values\n", catalog.filename);
+      status = FALSE;
+    }
+  }
+
+  // check lensing <-> average links
+  {
+    int objIDsOK = TRUE;
+    int catIDsOK = TRUE;
+    int averefOK = TRUE;
+
+    for (i = 0; i < catalog.Naverage; i++) {
+      m = catalog.average[i].lensingOffset;
+      for (j = 0; j < catalog.average[i].Nlensing; j++) {
+        objIDsOK &= (catalog.average[i].objID == catalog.lensing[m+j].objID);
+        catIDsOK &= (catalog.average[i].catID == catalog.lensing[m+j].catID);
+        averefOK &= (catalog.lensing[m+j].averef == i);
+      }
+    }
+    
+    if (!objIDsOK) {
+      fprintf (stderr, "ERROR: catalog %s has invalid lensing obj IDs\n", catalog.filename);
+      status = FALSE;
+    }
+    if (!catIDsOK) {
+      fprintf (stderr, "ERROR: catalog %s has invalid lensing cat IDs\n", catalog.filename);
+      status = FALSE;
+    }
+    if (!averefOK) {
+      fprintf (stderr, "ERROR: catalog %s has invalid lensing averef values\n", catalog.filename);
+      status = FALSE;
+    }
+  }
+
+  // check lensobj <-> average links
+  {
+    int objIDsOK = TRUE;
+    int catIDsOK = TRUE;
+
+    for (i = 0; i < catalog.Naverage; i++) {
+      m = catalog.average[i].lensobjOffset;
+      for (j = 0; j < catalog.average[i].Nlensobj; j++) {
+        objIDsOK &= (catalog.average[i].objID == catalog.lensobj[m+j].objID);
+        catIDsOK &= (catalog.average[i].catID == catalog.lensobj[m+j].catID);
+      }
+    }
+    
+    if (!objIDsOK) {
+      fprintf (stderr, "ERROR: catalog %s has invalid lensobj obj IDs\n", catalog.filename);
+      status = FALSE;
+    }
+    if (!catIDsOK) {
+      fprintf (stderr, "ERROR: catalog %s has invalid lensobj cat IDs\n", catalog.filename);
+      status = FALSE;
+    }
   }
 
@@ -321 +428 @@
     }
   }
+
+  for (i = 0; i < catalog[0].Naverage; i++) {
+    m = catalog[0].average[i].lensingOffset;
+    for (j = 0; j < catalog[0].average[i].Nlensing; j++) {
+      id = catalog[0].lensing[m+j].imageID;
+      if (id == 0) continue; // detections from ref photcodes can (should) have unset image IDs
+      if (id > maxID) {
+        Nfail ++;
+        continue;
+        // is this sufficient to catch IDs set without an image table?
+      }
+      if (IDlist) {
+        if (IDlist[id] < 0) {
+          Nfail ++;
+          continue;
+        }
+      } else {
+        if (id > 0) {
+          Nfail ++;
+          continue;
+        }
+      }
+    }
+  }
+
   return Nfail;
 }
@@ -382 +514 @@
   return TRUE;
 }
+
+# define GET_COLUMN(OUT,NAME,TYPE) \
+  OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \
+  myAssert (!strcmp(type, #TYPE), "wrong column type");
+
+// STATUS is value expected for success
+# define CHECK_STATUS(STATUS,MSG,...)                                   \
+  if (!(STATUS)) {                                                      \
+    fprintf (stderr, MSG, __VA_ARGS__);                                 \
+    return FALSE;                                                       \
+  }
+
+// write out the IDmap data for clients to read
+int SaveImageIDsSmall(char *filename) {
+
+  Header header;
+  Header theader;
+  Matrix matrix;
+  FTable ftable;
+
+  gfits_init_header (&header);
+  header.extend = TRUE;
+  gfits_create_header (&header);
+  gfits_create_matrix (&header, &matrix);
+
+  gfits_create_table_header (&theader, "BINTABLE", "IMAGE_IDS");
+
+  gfits_define_bintable_column (&theader, "J", "IMAGE_IDS", "image IDs", NULL, 1.0, FT_BZERO_INT32);
+
+  // generate the output array that carries the data
+  gfits_create_table (&theader, &ftable);
+
+  // add the columns to the output array
+  gfits_set_bintable_column (&theader, &ftable, "IMAGE_IDS", IDlist, maxID + 1);
+
+  FILE *f = fopen (filename, "w");
+  if (!f) {
+    fprintf (stderr, "ERROR: cannot open image ID file for output %s\n", filename);
+    return FALSE;
+  }
+
+  int status;
+  status = gfits_fwrite_header  (f, &header);
+  CHECK_STATUS (status, "ERROR: cannot write header for image ID file %s\n", filename);
+
+  status = gfits_fwrite_matrix  (f, &matrix);
+  CHECK_STATUS (status, "ERROR: cannot write matrix for image ID file %s\n", filename);
+
+  status = gfits_fwrite_Theader (f, &theader);
+  CHECK_STATUS (status, "ERROR: cannot write table header for image ID file %s\n", filename);
+
+  status = gfits_fwrite_table  (f, &ftable);
+  CHECK_STATUS (status, "ERROR: cannot write table data for image ID file %s\n", filename);
+
+  gfits_free_header (&header);
+  gfits_free_matrix (&matrix);
+  gfits_free_header (&theader);
+  gfits_free_table (&ftable);
+
+  int fd = fileno (f);
+
+  status = fflush (f);
+  CHECK_STATUS (!status, "ERROR: cannot flush file image ID file %s\n", filename);
+
+  status = fsync (fd);
+  CHECK_STATUS (!status, "ERROR: cannot flush file image ID file %s\n", filename);
+
+  status = fclose (f);
+  CHECK_STATUS (!status, "ERROR: problem closing image ID file file %s\n", filename);
+
+  return TRUE;
+}
+
+int LoadImageIDsSmall (char *filename) {
+
+  int Ncol;
+  off_t Nrow;
+  Header header;
+  Header theader;
+  Matrix matrix;
+  FTable ftable;
+
+  FILE *f = fopen (filename, "r");
+  if (!f) {
+    fprintf (stderr, "ERROR: cannot open image subset file %s\n", filename);
+    return FALSE;
+  }
+
+  /* load in PHU segment (ignore) */
+  if (!gfits_fread_header (f, &header)) {
+    if (VERBOSE) fprintf (stderr, "can't read image subset header\n");
+    fclose (f);
+    return FALSE;
+  }
+  if (!gfits_fread_matrix (f, &matrix, &header)) {
+    if (VERBOSE) fprintf (stderr, "can't read image subset matrix\n");
+    gfits_free_header (&header);
+    fclose (f);
+    return FALSE;
+  }
+
+  ftable.header = &theader;
+
+  // load data for this header 
+  if (!gfits_load_header (f, &theader)) {
+    fclose (f);
+    return FALSE;
+  }
+
+  // read the fits table bytes
+  if (!gfits_fread_ftable_data (f, &ftable, FALSE)) {
+    fclose (f);
+    return FALSE;
+  }
+  fclose (f);
+
+  char type[16];
+
+  GET_COLUMN (IDlist, "IMAGE_IDS", int);
+  maxID = Nrow - 1;
+  fprintf (stderr, "loaded data for %lld images\n", (long long) Nrow);
+
+  return TRUE;
+}
+

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/merge_catalogs_new.c

@@ -8 +8 @@
 // input entries always define new objects
 
-#define notyet 1
-
 int merge_catalogs_new (SkyRegion *region, Catalog *output, Catalog *input, int *secfiltMap) {
   
   off_t i, j, offset;
-  off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, NsecfiltIn, NsecfiltOut, Nm;
+  off_t NAVERAGE, NMEASURE, NLENSING, Naverage, Nmeasure, Nlensing, NsecfiltIn, NsecfiltOut, Nm;
 
   Naverage = output[0].Naverage;
@@ -63 +61 @@
       }
       output[0].average[Naverage].Nmeasure = Nm;
+
+      Nm = 0;
+      for (j = 0; j < input[0].average[i].Nlensing; j++) {
+          offset = input[0].average[i].lensingOffset + j;
+
+          output[0].lensing[Nlensing] = input[0].lensing[offset];
+          output[0].lensing[Nlensing].averef = Naverage;
+
+          Nlensing ++;
+          Nm ++;
+          if (Nlensing == NLENSING) {
+              NLENSING += 1000;
+              REALLOCATE (output[0].lensing, Lensing, NLENSING);
+          }
+      }
+      output[0].average[Naverage].Nlensing = Nm;
+
       Naverage ++;
       if (Naverage == NAVERAGE) {
@@ -72 +87 @@
   REALLOCATE (output[0].average, Average, MAX (Naverage, 1));
   REALLOCATE (output[0].measure, Measure, MAX (Nmeasure, 1));
+  REALLOCATE (output[0].lensing, Lensing, MAX (Nlensing, 1));
   REALLOCATE (output[0].secfilt, SecFilt, NsecfiltOut*MAX (Naverage, 1));
   output[0].Naverage = Naverage;
-  output[0].Nmeasure = Nmeasure;
+  output[0].Nlensing = Nlensing;
   output[0].Nsecf_mem = Naverage * NsecfiltOut;
 
@@ -84 +100 @@
   
   if (VERBOSE) {
-      fprintf (stderr, OFF_T_FMT": using "OFF_T_FMT" stars ("OFF_T_FMT" measures) for catalog\n", 
+      fprintf (stderr, OFF_T_FMT": using "OFF_T_FMT" stars ("OFF_T_FMT" measures, "OFF_T_FMT" lensing) for catalog\n", 
                 i, 
                 output[0].Naverage, 
-                output[0].Nmeasure);
+                output[0].Nmeasure, 
+                output[0].Nlensing);
   }
   return (TRUE);

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/merge_catalogs_old.c

@@ -1 +1 @@
 # include "dvomerge.h"
 # define PSPS_ID TRUE
-
-# define MARKTIME(MSG,...) {                    \
-    float dtime;                                \
-    gettimeofday (&stop, (void *) NULL);        \
-    dtime = DTIME (stop, start);                \
-    fprintf (stderr, MSG, __VA_ARGS__);         \
-    gettimeofday (&start, (void *) NULL);       \
-  }
 
 # define IN_REGION(R,D) (                                       \
@@ -19 +11 @@
 
   off_t i, j, k, Nin, offset, J, Jmin, status, Nstars;
-  double RADIUS2, Rmin, Rin, Din;
+  double RADIUS2, Rmin;
   double *X1, *Y1, *X2, *Y2;
   double dX, dY, dR;
-  off_t *N1, *N2,  *next_meas;
-  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch;
+  off_t *N1, *N2, *next_meas, *next_lens;
+  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch, Nlens, NLENS;
   int NsecfiltIn;
   int NsecfiltOut;
@@ -29 +21 @@
   Coords tcoords;
   
-  // struct timeval start, stop;
-  // gettimeofday (&start, (void *) NULL);
+  // INITTIME;
 
   NsecfiltOut = output[0].Nsecfilt;
@@ -40 +31 @@
   ALLOCATE (Y1, double, input[0].Naverage);
   ALLOCATE (N1, off_t,  input[0].Naverage);
-  ALLOCATE (input[0].found, off_t, input[0].Naverage);
+  if (!input[0].found_t) {
+    ALLOCATE (input[0].found_t, off_t, input[0].Naverage);
+  } else {
+    REALLOCATE (input[0].found_t, off_t, input[0].Naverage);
+  }
 
   /** allocate local arrays (catalog) **/
@@ -47 +42 @@
   ALLOCATE (Y2, double, NAVE);
   ALLOCATE (N2, off_t,    NAVE);
-  ALLOCATE (output[0].found, off_t, NAVE);
+  if (!output[0].found_t) {
+    ALLOCATE (output[0].found_t, off_t, NAVE);
+  } else {
+    REALLOCATE (output[0].found_t, off_t, NAVE);
+  }
   /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */
 
@@ -53 +52 @@
   Nmatch = 0;
   NMEAS = Nmeas = output[0].Nmeasure;
+  NLENS = Nlens = output[0].Nlensing;
 
   // current max obj ID for this catalog
@@ -87 +87 @@
     N1[Nstars] = i;
     Nstars ++;
-    input[0].found[i] = -1;
+    input[0].found_t[i] = -1;
   }
   if (Nstars < 1) {
     if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", output[0].filename);
-    free (output[0].found);
-    free (input[0].found);
     free (X1);
     free (Y1);
@@ -106 +104 @@
   for (i = 0; i < Nave; i++) {
     RD_to_XY (&X2[i], &Y2[i], output[0].average[i].R, output[0].average[i].D, &tcoords);
-    output[0].found[i] = -1;
+    output[0].found_t[i] = -1;
     N2[i] = i;
   }
@@ -115 +113 @@
     exit (3);
   }
-
 
   /* set up pointers for linked list of measure */
@@ -122 +119 @@
     // is sorted while processed
     next_meas = init_measure_links (output[0].average, Nave, output[0].measure, Nmeas);
+    next_lens = init_lensing_links (output[0].average, Nave, output[0].lensing, Nlens);
   } else {
     next_meas = build_measure_links (output[0].average, Nave, output[0].measure, Nmeas);
+    next_lens = build_lensing_links (output[0].average, Nave, output[0].lensing, Nlens);
   }    
 
@@ -185 +184 @@
       REALLOCATE (output[0].measure, Measure, NMEAS);
     }
+    if (Nlens + input[0].average[N].Nmeasure >= NLENS) {
+      NLENS = Nlens + input[0].average[N].Nmeasure + 1000;
+      REALLOCATE (next_lens, off_t, NLENS);
+      REALLOCATE (output[0].lensing, Lensing, NLENS);
+    }
 
     // 4) average properties from the input and the output db need to be properly merged.
@@ -194 +198 @@
       if (REPLACE_BY_PHOTCODE) {
         // index to first measure for this object
+        // XXX this does not support lensing measurements
         int Mout = output[0].average[n].measureOffset;  
         if (replace_match(&output[0].average[n], &output[0].measure[Mout], &input[0].average[N], &input[0].measure[offset])) {
-          input[0].found[N] = Mout;
+          input[0].found_t[N] = Mout;
           continue;
         }
@@ -206 +211 @@
       output[0].measure[Nmeas] = input[0].measure[offset];
 
-      Rin = input[0].average[N].R - input[0].measure[offset].dR / 3600.0;
-      Din = input[0].average[N].D - input[0].measure[offset].dD / 3600.0;
-
-      /** dR,dD now represent arcsec **/
-      output[0].measure[Nmeas].dR       = 3600.0*(output[0].average[n].R - Rin);  // XXX update these based on choice of astromety
-      output[0].measure[Nmeas].dD       = 3600.0*(output[0].average[n].D - Din);  // XXX update these based on choice of astromety
+      // old code: find R,D using average_in[0], the get offset relative to average_out[0].  no longer
+      // needed since we carry around R,D
+      // Rin = input[0].average[N].R - input[0].measure[offset].dR / 3600.0;
+      // Din = input[0].average[N].D - input[0].measure[offset].dD / 3600.0;
+      // output[0].measure[Nmeas].dR = 3600.0*(output[0].average[n].R - Rin);
+      // output[0].measure[Nmeas].dD = 3600.0*(output[0].average[n].D - Din);
+
       output[0].measure[Nmeas].dbFlags  = 0;  // XXX why reset these?
       output[0].measure[Nmeas].averef   = n;
@@ -221 +227 @@
       // fprintf (stderr, "Nave : "OFF_T_FMT", Nmeas : "OFF_T_FMT", dR: %f, dD: %f, catID: %d\n",  n,  Nmeas, output[0].measure[Nmeas].dR, output[0].measure[Nmeas].dD, output[0].measure[i].catID);
 
-      // rationalize dR
-      if (output[0].measure[Nmeas].dR > +180.0*3600.0) {
+      float dRoff = dvoOffsetR(&output[0].measure[Nmeas], &output[0].average[n]);
+
+      // rationalize R
+      if (dRoff > +180.0*3600.0) {
         // average on high end of boundary, move star up
-        Rin += 360.0;
-        output[0].measure[Nmeas].dR = 3600.0*(output[0].average[n].R - Rin);
-      }
-      if (output[0].measure[Nmeas].dR < -180.0*3600.0) {
+        output[0].measure[Nmeas].R += 360.0;
+        dRoff -= 360.0*3600.0;
+      }
+      if (dRoff < -180.0*3600.0) {
         // average on low end of boundary, move star down
-        Rin -= 360.0;
-        output[0].measure[Nmeas].dR = 3600.0*(output[0].average[n].R - Rin);
-      }
-      if (fabs(output[0].measure[Nmeas].dR) > 10*RADIUS) {
-        // ok take declination into account and check again.
-        double cosD = cos(RAD_DEG*Din);
-        if (fabs(output[0].measure[Nmeas].dR*cosD) > 10*RADIUS) {
-
-            fprintf (stderr, "error: %10.6f,%10.6f vs %10.6f,%10.6f (%f,%f vs %f,%f)\n", 
-             output[0].average[n].R, output[0].average[n].D, Rin, Din,
-             X1[i], X2[Jmin], Y1[i], Y2[Jmin]);
-
-            // XXX abort on this? -- this is a bad failure...
+        output[0].measure[Nmeas].R -= 360.0;
+        dRoff += 360.0*3600.0;
+      }
+      if (fabs(dRoff) > 10*RADIUS) {
+        // take declination into account and check again.
+        double cosD = cos(RAD_DEG*output[0].average[n].D);
+        if (fabs(dRoff*cosD) > 10*RADIUS) {
+          fprintf (stderr, "error: %10.6f,%10.6f vs %10.6f,%10.6f (%f,%f vs %f,%f)\n", 
+                   output[0].average[n].R, output[0].average[n].D, 
+                   output[0].measure[Nmeas].R, output[0].measure[Nmeas].D,
+                   X1[i], X2[Jmin], Y1[i], Y2[Jmin]);
+          // XXX abort on this? -- this is a bad failure...
         }
       }
-      input[0].found[N] = Nmeas;
+      input[0].found_t[N] = Nmeas;
       output[0].average[n].Nmeasure ++;
       Nmeas ++;
+    }
+
+    // if lensing measurements exist, add them too
+    if (output[0].lensing) {
+      for (Nin = 0; Nin < input[0].average[N].Nlensing; Nin++) {
+        /* add to end of lensing list */
+        add_lens_link (&output[0].average[n], next_lens, Nlens, NLENS);
+        
+        // set the new lensing
+        off_t lensoff = input[0].average[N].lensingOffset + Nin;
+        output[0].lensing[Nlens] = input[0].lensing[lensoff];
+
+        output[0].lensing[Nlens].averef   = n;
+        output[0].lensing[Nlens].objID    = output[0].average[n].objID;
+        output[0].lensing[Nlens].catID    = output[0].catID;
+        output[0].average[n].Nlensing ++;
+        Nlens ++;
+      }
     }
 
@@ -273 +298 @@
     /* Nm is updated, but not written out in -update mode (for existing entries)
        Nm is recalculated in build_meas_links if loaded table is not sorted */
-    output[0].found[n] = Nmeas;
+    output[0].found_t[n] = Nmeas;
     i++;
   }
@@ -289 +314 @@
       REALLOCATE (output[0].measure, Measure, NMEAS);
     }
+    if (Nlens + input[0].average[N].Nlensing >= NLENS) {
+      NLENS = Nlens + input[0].average[N].Nlensing + 1000;
+      REALLOCATE (next_lens, off_t, NLENS);
+      REALLOCATE (output[0].lensing, Lensing, NLENS);
+    }
     if (Nave >= NAVE) {
       NAVE = Nave + 1000;
@@ -295 +325 @@
     }
 
-    if (input[0].found[N] >= 0) continue;
+    if (input[0].found_t[N] >= 0) continue;
     if (!IN_REGION (input[0].average[N].R, input[0].average[N].D)) continue;
 
@@ -348 +378 @@
 
       /* we set next[Nmeas] to -1 here, and update correctly below */
-      input[0].found[N] = Nmeas;
+      input[0].found_t[N] = Nmeas;
       next_meas[Nmeas] = -1;
       Nmeas ++;
@@ -356 +386 @@
       next_meas[Nmeas - Ngroup + j] = Nmeas - Ngroup + j + 1;
     }
+
+    /** add lensing for this input average object **/
+    if (output[0].lensing) {
+      output[0].average[Nave].lensingOffset  = Nlens;
+      for (Nin = 0; Nin < input[0].average[N].Nlensing; Nin ++) {
+        // supply the lensing values from this detection
+        off_t lensoff = input[0].average[N].lensingOffset + Nin;
+        output[0].lensing[Nlens]           = input[0].lensing[lensoff];
+
+        // the following lensing elements cannot be set until here:
+        output[0].lensing[Nlens].averef   = Nave;
+        output[0].lensing[Nlens].objID    = output[0].average[Nave].objID;
+        output[0].lensing[Nlens].catID    = output[0].catID;
+
+        // as we add lensing, update Nlensing to match
+        output[0].average[Nave].Nlensing ++;
+
+        /* we set next[Nlens] to -1 here, and update correctly below */
+        next_lens[Nlens] = -1;
+        Nlens ++;
+      }
+      int Ngroup = input[0].average[N].Nlensing;
+      for (j = 0; j < Ngroup - 1; j++) {
+        next_lens[Nlens - Ngroup + j] = Nlens - Ngroup + j + 1;
+      }
+    }
+
     Nave ++;
   }
@@ -363 +420 @@
   REALLOCATE (output[0].average, Average, Nave);
   REALLOCATE (output[0].measure, Measure, Nmeas);
+  REALLOCATE (output[0].lensing, Lensing, Nlens);
  
 # define NOSORT 0
@@ -370 +428 @@
     output[0].sorted = TRUE;
     output[0].measure = sort_measure (output[0].average, Nave, output[0].measure, Nmeas, next_meas);
+    output[0].lensing = sort_lensing (output[0].average, Nave, output[0].lensing, Nlens, next_lens);
   }
 
   /* note stars which have been found in this catalog */
   for (i = 0; i < input[0].Naverage; i++) {
-    if (input[0].found[i] > -1) {
-      input[0].found[i] = -2;
+    if (input[0].found_t[i] > -1) {
+      input[0].found_t[i] = -2;
     } else {
-      input[0].found[i] = -3;
+      input[0].found_t[i] = -3;
     }
   }
@@ -385 +444 @@
   output[0].Naverage = Nave;
   output[0].Nmeasure = Nmeas;
+  output[0].Nlensing = Nlens;
   output[0].Nsecf_mem = Nave*NsecfiltOut;
-  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nmatch);
+  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nlens: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nlens, Nmatch);
 
   free (next_meas);
-
-  free (output[0].found);
+  free (next_lens);
+
   free (X2);
   free (Y2);
   free (N2);
-  free (input[0].found);
   free (X1);
   free (Y1);
@@ -411 +470 @@
    images have boundaries which are lines in pixels coords, but curve in RA and DEC
    
-   output[0].found[Ncat] but stars[Nstar].found
+   output[0].found_t[Ncat] but stars[Nstar].found
    
 */

branches/eam_branches/ps2-tc3-20130727/Ohana/src/dvomerge/src/replace_match.c

@@ -7 +7 @@
   unsigned int averef;
   unsigned int catID;
-  double Rin, Din;
 
   // find the matching photcode in the object's list of measurements
@@ -18 +17 @@
     measure_out[Nout] = measure_in[0];
 
-    Rin = average_in[0].R - measure_in[0].dR / 3600.0;
-    Din = average_in[0].D - measure_in[0].dD / 3600.0;
+    // old code: find R,D using average_in[0], the get offset relative to average_out[0].  no longer
+    // needed since we carry around R,D
+    // double Rin = average_in[0].R - measure_in[0].dR / 3600.0;
+    // double Din = average_in[0].D - measure_in[0].dD / 3600.0;
+    // measure_out[Nout].dR = 3600.0*(average_out[0].R - Rin);
+    // measure_out[Nout].dD = 3600.0*(average_out[0].D - Din);
 
-    /** dR,dD now represent arcsec **/
-    measure_out[Nout].dR       = 3600.0*(average_out[0].R - Rin);
-    measure_out[Nout].dD       = 3600.0*(average_out[0].D - Din);
     measure_out[Nout].dbFlags  = 0;  // XXX why reset these?
     measure_out[Nout].averef   = averef;
@@ -29 +29 @@
     measure_out[Nout].catID    = catID;
 
+    float dRoff = dvoOffsetR(&measure_out[Nout], average_out);
+
     // rationalize dR
-    if (measure_out[Nout].dR > +180.0*3600.0) {
+    if (dRoff > +180.0*3600.0) {
       // average on high end of boundary, move star up
-      Rin += 360.0;
-      measure_out[Nout].dR = 3600.0*(average_out[0].R - Rin);
+      measure_out[Nout].R += 360.0;
+      dRoff -= 360.0*3600.0;
     }
-    if (measure_out[Nout].dR < -180.0*3600.0) {
+    if (dRoff < -180.0*3600.0) {
       // average on low end of boundary, move star down
-      Rin -= 360.0;
-      measure_out[Nout].dR = 3600.0*(average_out[0].R - Rin);
+      measure_out[Nout].R -= 360.0;
+      dRoff += 360.0*3600.0;
     }
 
     // warn on surprisingly distant detections
-    if (fabs(measure_out[Nout].dR) > 10*RADIUS) {
+    if (fabs(dRoff) > 10*RADIUS) {
       // ok take declination into account and check again.
-      double cosD = cos(RAD_DEG*Din);
-      if (fabs(measure_out[Nout].dR*cosD) > 10*RADIUS) {
-
+      double cosD = cos(RAD_DEG*average_out[0].D);
+      if (fabs(dRoff*cosD) > 10*RADIUS) {
         fprintf (stderr, "surprisingly distant detection: %10.6f,%10.6f vs %10.6f,%10.6f\n", 
-                 average_out[0].R, average_out[0].D, Rin, Din);
+                 average_out[0].R, average_out[0].D, measure_out[Nout].R, measure_out[Nout].D);
       }
     }