Changeset 39142 for trunk

Timestamp:

Nov 15, 2015, 10:13:54 AM (11 years ago)

Author:

eugene

Message:

add galphot to the merge options

Location:

trunk/Ohana/src/dvomerge

Files:

• include/dvomerge.h (modified) (1 diff)
• src/build_links.c (modified) (5 diffs)
• src/dvo_image_merge_dbs.c (modified) (1 diff)
• src/merge_catalogs_new.c (modified) (7 diffs)
• src/merge_catalogs_old.c (modified) (9 diffs)

trunk/Ohana/src/dvomerge/include/dvomerge.h

@@ -164 +164 @@
 StarPar   *sort_starpar           PROTO((Average *average, off_t Naverage, StarPar *starpar, off_t Nstarpar, off_t *next));
 
+off_t     *build_galphot_links    PROTO((Average *average, off_t Naverage, GalPhot *galphot, off_t Ngalphot));
+off_t     *init_galphot_links     PROTO((Average *average, off_t Naverage, GalPhot *galphot, off_t Ngalphot));
+int        add_galp_link          PROTO((Average *average, off_t *next, off_t Ngalphot, off_t NGALPHOT));
+GalPhot   *sort_galphot           PROTO((Average *average, off_t Naverage, GalPhot *galphot, off_t Ngalphot, off_t *next));
+
 off_t     *init_missing_links     PROTO((Average *average, off_t Naverage, Missing *missing, off_t Nmissing));
 int        add_miss_link          PROTO((Average *average, off_t *next, off_t Nmissing));

trunk/Ohana/src/dvomerge/src/build_links.c

@@ -23 +23 @@
 
 */
+
+/*** Measure ****************************************************************************************/
 
 /* build the initial links assuming the table is sorted, 
@@ -159 +161 @@
 }
 
-/*******************************************************************************************/
+/*** Missing ****************************************************************************************/
 
 /* build the initial links assuming the table is sorted */
@@ -229 +231 @@
 }
 
-/*******************************************************************************************/
+/*** Lensing ****************************************************************************************/
 
 /* build the initial links assuming the table is sorted, 
@@ -384 +386 @@
 }
 
-/********  StarPar ********************************************************************************/
+/*** StarPar ********************************************************************************/
 
 /* build the initial links assuming the table is sorted, 
@@ -535 +537 @@
 }
 
-
+/*** GalPhot ****************************************************************************************/
+
+/* build the initial links assuming the table is sorted, 
+   not partial, and has a correct set of average[].galphotOffset,Ngalphot values */
+off_t *init_galphot_links (Average *average, off_t Naverage, GalPhot *galphot, off_t Ngalphot) {
+
+  off_t i, j, N;
+  off_t *next_galp;
+
+  N = 0;
+
+  ALLOCATE (next_galp, off_t, Ngalphot);
+  for (i = 0; i < Naverage; i++) {
+    if (!average[i].Ngalphot) continue;
+    off_t m = average[i].galphotOffset;
+    myAssert (galphot[m].averef == i, "not sorted");
+    for (j = 0; j < average[i].Ngalphot - 1; j++, N++) {
+      myAssert (galphot[m+j+1].averef == i, "not sorted");
+      next_galp[N] = N + 1;
+      if (N >= Ngalphot) {
+        fprintf (stderr, "WARNING: N out of bounds (1)\n");
+      }
+    }
+    next_galp[N] = -1;
+    if (N >= Ngalphot) {
+      fprintf (stderr, "WARNING: N out of bounds (2)\n");
+    }
+
+    if (N >= Ngalphot) {
+      fprintf (stderr, "overflow in init_galphot_links\n");
+      abort ();
+    }
+    N++;
+  }
+  return (next_galp);
+}
+
+/* construct galphot 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_galphot_links (Average *average, off_t Naverage, GalPhot *galphot, off_t Ngalphot) {
+
+  off_t i, m, k, Nm, averef;
+  off_t *next_galp;
+
+  ALLOCATE (next_galp, off_t, Ngalphot);
+
+  /* reset the Nm, offset values for average */
+  for (i = 0; i < Naverage; i++) {
+    average[i].galphotOffset = -1;
+    average[i].Ngalphot     =  0;
+  }
+
+  for (Nm = 0; Nm < Ngalphot; Nm++) {
+    averef = galphot[Nm].averef;
+    m = average[averef].galphotOffset;  
+    next_galp[Nm] = -1;
+
+    if (m == -1) { /* no links yet for source */
+      average[averef].galphotOffset = Nm;
+      average[averef].Ngalphot     = 1;
+      continue;
+    }
+
+    for (k = 0; next_galp[m] != -1; k++) {
+      m = next_galp[m];
+      if (m >= Ngalphot) {
+        fprintf (stderr, "WARNING: m out of bounds (1)\n");
+      }
+    }
+
+    average[averef].Ngalphot = k + 2;
+    next_galp[m] = Nm;
+    if (m >= Ngalphot) {
+      fprintf (stderr, "WARNING: m out of bounds (2)\n");
+    }
+  }
+  return (next_galp);
+}
+
+/* average[].galphotOffset, average[].Ngalphot are valid within an addstar run */
+int add_galp_link (Average *average, off_t *next_galp, off_t Ngalphot, off_t NGALPHOT) {
+
+  off_t k, m;
+
+  /* if we have trouble, check validity of next_galp[m] : m < Ngalphot */
+  m = average[0].galphotOffset;  
+
+  for (k = 0; k < average[0].Ngalphot - 1; k++)  {
+    m = next_galp[m];
+    if (m >= NGALPHOT) {
+      fprintf (stderr, "WARNING: m out of bounds (3)\n");
+    }
+  }
+
+  /* set up references */
+  next_galp[Ngalphot] = -1;
+  if (Ngalphot >= NGALPHOT) {
+    fprintf (stderr, "WARNING: Ngalphot out of bounds (1)\n");
+  }
+
+  // if Ngalphot is 0, m may have been mis-set; add to the end
+  if ((average[0].Ngalphot == 0) || (m == -1)) {
+    average[0].galphotOffset = Ngalphot;
+  } else {
+    next_galp[m] = Ngalphot;
+    if (m >= NGALPHOT) {
+      fprintf (stderr, "WARNING: m out of bounds (4)\n");
+    }
+  }
+
+  return (TRUE);
+}
+
+GalPhot *sort_galphot (Average *average, off_t Naverage, GalPhot *galphot, off_t Ngalphot, off_t *next_galp) {
+
+  off_t i, k, n, np, N;
+  GalPhot *tmpgalphot;
+
+  /*
+  for (i = 0; i < Naverage; i++) {
+    if (average[i].Ngalphot != 4) {
+      fprintf (stderr, "check %d %d %d\n", (int) i, (int) average[i].Ngalphot, (int) average[i].galphotOffset);
+    }
+  }
+  */
+
+  /* fix order of GalPhot (memory intensive, but fast) */
+  np = -1;
+  N = 0; 
+  ALLOCATE (tmpgalphot, GalPhot, Ngalphot);
+  for (i = 0; i < Naverage; i++) {
+    if (!average[i].Ngalphot) continue;
+    n = average[i].galphotOffset;
+    average[i].galphotOffset = N;
+    for (k = 0; k < average[i].Ngalphot; k++, N++) {
+      if (n == -1) {
+        fprintf (stderr, "entry after %d has a problem\n", (int) np);
+        abort();
+      }
+      tmpgalphot[N] = galphot[n]; 
+      if (galphot[n].averef != i) abort();
+      tmpgalphot[N].averef = i;
+      np = n;
+      n = next_galp[n];
+    }
+  }
+  free (galphot);
+  return (tmpgalphot);
+}
+

trunk/Ohana/src/dvomerge/src/dvo_image_merge_dbs.c

@@ -234 +234 @@
     catalog[0].lensing[i].imageID = newID;
   }
+
+  // NOTE galphot also carries an imageID, but it is actually the EXTERNAL image ID reference.  Do NOT update it.
+
   return TRUE;
 }

trunk/Ohana/src/dvomerge/src/merge_catalogs_new.c

@@ -11 +11 @@
   
   off_t i, j, offset;
-  off_t NAVERAGE, NMEASURE, NLENSING, NSTARPAR, Naverage, Nmeasure, Nlensing, Nstarpar, NsecfiltIn, NsecfiltOut, Nm;
+  off_t NAVERAGE, NMEASURE, NLENSING, NSTARPAR, NGALPHOT, Naverage, Nmeasure, Nlensing, Nstarpar, Ngalphot, NsecfiltIn, NsecfiltOut, Nm;
 
   Naverage = output[0].Naverage;
@@ -17 +17 @@
   Nlensing = output[0].Nlensing;
   Nstarpar = output[0].Nstarpar;
+  Ngalphot = output[0].Ngalphot;
 
   NsecfiltOut = output[0].Nsecfilt;
@@ -25 +26 @@
   NLENSING = Nlensing + 1000;
   NSTARPAR = Nstarpar + 1000;
+  NGALPHOT = Ngalphot + 1000;
 
   REALLOCATE (output[0].average, Average, NAVERAGE);
@@ -31 +33 @@
   REALLOCATE (output[0].lensing, Lensing, NLENSING);
   REALLOCATE (output[0].starpar, StarPar, NSTARPAR);
+  REALLOCATE (output[0].galphot, GalPhot, NGALPHOT);
 
   // add all of these entries to the existing catalog
@@ -100 +103 @@
       output[0].average[Naverage].Nstarpar = Nm;
 
+      Nm = 0;
+      for (j = 0; j < input[0].average[i].Ngalphot; j++) {
+          offset = input[0].average[i].galphotOffset + j;
+
+          output[0].galphot[Ngalphot] = input[0].galphot[offset];
+          output[0].galphot[Ngalphot].averef = Naverage;
+
+          Ngalphot ++;
+          Nm ++;
+          if (Ngalphot == NGALPHOT) {
+              NGALPHOT += 1000;
+              REALLOCATE (output[0].galphot, GalPhot, NGALPHOT);
+          }
+      }
+      output[0].average[Naverage].Ngalphot = Nm;
+
       Naverage ++;
       if (Naverage == NAVERAGE) {
@@ -111 +130 @@
   REALLOCATE (output[0].lensing, Lensing, MAX (Nlensing, 1));
   REALLOCATE (output[0].starpar, StarPar, MAX (Nstarpar, 1));
+  REALLOCATE (output[0].galphot, GalPhot, MAX (Ngalphot, 1));
   REALLOCATE (output[0].secfilt, SecFilt, NsecfiltOut*MAX (Naverage, 1));
   output[0].Naverage = Naverage;
   output[0].Nlensing = Nlensing;
   output[0].Nstarpar = Nstarpar;
+  output[0].Ngalphot = Ngalphot;
   output[0].Nsecfilt_mem = Naverage * NsecfiltOut;
 
@@ -124 +145 @@
   
   if (VERBOSE) {
-      fprintf (stderr, OFF_T_FMT": using "OFF_T_FMT" stars ("OFF_T_FMT" measures, "OFF_T_FMT" lensing) for catalog\n", 
+      fprintf (stderr, OFF_T_FMT": using "OFF_T_FMT" stars ("OFF_T_FMT" measures, "OFF_T_FMT" lensing, "OFF_T_FMT" starpar, "OFF_T_FMT" galphot) for catalog\n", 
                 i, 
                 output[0].Naverage, 
                 output[0].Nmeasure, 
-                output[0].Nlensing);
+                output[0].Nlensing, 
+                output[0].Nstarpar, 
+                output[0].Ngalphot);
   }
   return (TRUE);

trunk/Ohana/src/dvomerge/src/merge_catalogs_old.c

@@ -14 +14 @@
   double *X1, *Y1, *X2, *Y2;
   double dX, dY, dR;
-  off_t *N1, *N2, *next_meas, *next_lens, *next_star;
-  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch, Nlens, NLENS, Nstar, NSTAR;
+  off_t *N1, *N2, *next_meas, *next_lens, *next_star, *next_galp;
+  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch, Nlens, NLENS, Nstar, NSTAR, Ngalp, NGALP;
   int NsecfiltIn;
   int NsecfiltOut;
@@ -199 +199 @@
       NGALP = Ngalp + input[0].average[N].Ngalphot + 1000;
       REALLOCATE (next_galp, off_t, NGALP);
-      REALLOCATE (output[0].galphot, Galphot, NGALP);
+      REALLOCATE (output[0].galphot, GalPhot, NGALP);
     }
 
@@ -222 +222 @@
       if (REPLACE_BY_PHOTCODE) {
         // index to first measure for this object
-        // XXX this does not support lensing measurements
+        // XXX this does not support lensing, starpar, or galphot measurements
         if (replace_match (&output[0].average[n], output[0].measure, next_meas, &input[0].average[N], &input[0].measure[offset])) {
           input[0].found_t[N] = output[0].average[n].measureOffset;  
@@ -315 +315 @@
     }
 
+    // if galphot measurements exist, add them too
+    if (output[0].galphot) {
+      for (Nin = 0; Nin < input[0].average[N].Ngalphot; Nin++) {
+        /* add to end of galphot list */
+        add_galp_link (&output[0].average[n], next_galp, Ngalp, NGALP);
+        
+        // set the new galphot
+        off_t galpoff = input[0].average[N].galphotOffset + Nin;
+        output[0].galphot[Ngalp] = input[0].galphot[galpoff];
+
+        output[0].galphot[Ngalp].averef   = n;
+        output[0].galphot[Ngalp].objID    = output[0].average[n].objID;
+        output[0].galphot[Ngalp].catID    = output[0].catID;
+        output[0].average[n].Ngalphot ++;
+        Ngalp ++;
+      }
+    }
+
     // update the average properties to reflect the incoming entries:
     // if the original value is NAN but the input value is not, accept the input:
@@ -365 +383 @@
       REALLOCATE (output[0].starpar, StarPar, NSTAR);
     }
+    if (Ngalp + input[0].average[N].Ngalphot >= NGALP) {
+      NGALP = Ngalp + input[0].average[N].Ngalphot + 1000;
+      REALLOCATE (next_galp, off_t, NGALP);
+      REALLOCATE (output[0].galphot, GalPhot, NGALP);
+    }
     if (Nave >= NAVE) {
       NAVE = Nave + 1000;
@@ -488 +511 @@
     }
 
+    /** add galphot for this input average object **/
+    if (output[0].galphot) {
+      output[0].average[Nave].galphotOffset  = Ngalp;
+      for (Nin = 0; Nin < input[0].average[N].Ngalphot; Nin ++) {
+        // supply the galphot values from this detection
+        off_t galpoff = input[0].average[N].galphotOffset + Nin;
+        output[0].galphot[Ngalp]           = input[0].galphot[galpoff];
+
+        // the following galphot elements cannot be set until here:
+        output[0].galphot[Ngalp].averef   = Nave;
+        output[0].galphot[Ngalp].objID    = output[0].average[Nave].objID;
+        output[0].galphot[Ngalp].catID    = output[0].catID;
+
+        // as we add galphot, update Ngalphot to match
+        output[0].average[Nave].Ngalphot ++;
+
+        /* we set next[Ngalp] to -1 here, and update correctly below */
+        next_galp[Ngalp] = -1;
+        Ngalp ++;
+      }
+      int Ngroup = input[0].average[N].Ngalphot;
+      for (j = 0; j < Ngroup - 1; j++) {
+        next_galp[Ngalp - Ngroup + j] = Ngalp - Ngroup + j + 1;
+      }
+    }
+
     Nave ++;
   }
@@ -497 +546 @@
   REALLOCATE (output[0].lensing, Lensing, Nlens);
   REALLOCATE (output[0].starpar, StarPar, Nstar);
+  REALLOCATE (output[0].galphot, GalPhot, Ngalp);
  
 # define NOSORT 0
@@ -506 +556 @@
     output[0].lensing = sort_lensing (output[0].average, Nave, output[0].lensing, Nlens, next_lens);
     output[0].starpar = sort_starpar (output[0].average, Nave, output[0].starpar, Nstar, next_star);
+    output[0].galphot = sort_galphot (output[0].average, Nave, output[0].galphot, Ngalp, next_galp);
   }
 
@@ -523 +574 @@
   output[0].Nlensing = Nlens;
   output[0].Nstarpar = Nstar;
+  output[0].Ngalphot = Ngalp;
   output[0].Nsecfilt_mem = Nave*NsecfiltOut;
-  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);
+  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nlens, Ngalp: "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT" "OFF_T_FMT", ("OFF_T_FMT" matches)\n",  Nstars,  Nave,  Nmeas,  Nlens, Ngalp, Nmatch);
 
   free (next_meas);
   free (next_lens);
   free (next_star);
+  free (next_galp);
 
   free (X2);