Changeset 38636

Timestamp:

Jul 24, 2015, 8:26:59 AM (11 years ago)

Author:

eugene

Message:

working on loadgalphot

Location:

trunk/Ohana/src/addstar

Files:

• Makefile (modified) (4 diffs)
• include/addstar.h (modified) (1 diff)
• src/find_matches_galphot.c (modified) (2 diffs)
• src/find_matches_starpar.c (modified) (1 diff)
• src/resort_catalog.c (modified) (13 diffs)
• src/resort_threaded.c (modified) (1 diff)
• src/resort_unthreaded.c (modified) (1 diff)

trunk/Ohana/src/addstar/Makefile

@@ -40 +40 @@
 # programs in 'SERVER' use the client-server concept and are out of date
 
-INSTALL = addstar sedstar loadstarpar loadstarpar_client loadICRF loadICRF_client skycells mkcmf dumpskycells findskycell
+INSTALL = addstar sedstar loadgalphot loadstarpar loadstarpar_client loadICRF loadICRF_client skycells mkcmf dumpskycells findskycell
 OTHERS  = load2mass loadwise loadtycho loadsupercos
 SERVER  = addstar_client addstarc addstard addstart
@@ -271 +271 @@
 $(SRC)/loadgalphot_join.$(ARCH).o \
 $(SRC)/loadgalphot_table.$(ARCH).o \
+$(SRC)/resort_catalog.$(ARCH).o \
 $(SRC)/strhash.$(ARCH).o \
 $(SRC)/sortIDs.$(ARCH).o \
@@ -295 +296 @@
 $(SRC)/loadstarpar_save_remote.$(ARCH).o \
 $(SRC)/loadstarpar_table.$(ARCH).o \
+$(SRC)/resort_catalog.$(ARCH).o \
 $(SRC)/psps_ids.$(ARCH).o
 
@@ -312 +314 @@
 $(SRC)/loadstarpar_save_remote.$(ARCH).o \
 $(SRC)/loadstarpar_table.$(ARCH).o \
+$(SRC)/resort_catalog.$(ARCH).o \
 $(SRC)/psps_ids.$(ARCH).o
 

trunk/Ohana/src/addstar/include/addstar.h

@@ -359 +359 @@
 Catalog *addstar_catalog_init (int Nstars);
 
+void resort_catalog_measure (Catalog *catalog);
+void resort_catalog_lensing (Catalog *catalog);
+void resort_catalog_starpar (Catalog *catalog);
+void resort_catalog_galphot (Catalog *catalog);
+
 /** 
     there is an inconsistency to be resolved: fixed structures (like Image)

trunk/Ohana/src/addstar/src/find_matches_galphot.c

@@ -199 +199 @@
   REALLOCATE (catalog[0].galphot, GalPhot, Ngalphot);
  
-  catalog[0].sorted = FALSE;
-
   /* check if the catalog has changed?  if no change, no need to write */
   catalog[0].objID     = objID; // new max value, save on catalog close
@@ -206 +204 @@
   catalog[0].Ngalphot  = Ngalphot;
   catalog[0].Nsecfilt_mem = Nave*Nsecfilt;
+
+  // we need to for
+  catalog[0].sorted = FALSE;
+  resort_catalog_galphot (catalog);
+  catalog[0].sorted = TRUE;
+
   if (VERBOSE) fprintf (stderr, "Nstars, Nave, Ngalphot: %d "OFF_T_FMT" "OFF_T_FMT" ("OFF_T_FMT" matches)\n",  Nstars, Nave, Ngalphot, Nmatch);
 

trunk/Ohana/src/addstar/src/find_matches_starpar.c

@@ -199 +199 @@
   REALLOCATE (catalog[0].starpar, StarPar, Nstarpar);
  
-  catalog[0].sorted = FALSE;
+  resort_catalog_starpar (catalog);
+  catalog[0].sorted = TRUE;
 
   /* check if the catalog has changed?  if no change, no need to write */

trunk/Ohana/src/addstar/src/resort_catalog.c

@@ -2 +2 @@
 
 void SortAveMatch (off_t *MEAS, off_t *AVE, off_t N);
-void resort_catalog_measure (Catalog *catalog);
-void resort_catalog_lensing (Catalog *catalog);
-void resort_catalog_starpar (Catalog *catalog);
-
-void resort_catalog_old (Catalog *catalog) {
-
-  off_t *next_meas, *next_lens;
-  off_t Naves, Nmeas, Nlens;
-
-  if (catalog[0].sorted == TRUE) return;
-
-  INITTIME;
-
-  /* internal counters */
-  Nmeas = catalog[0].Nmeasure;
-  Nlens = catalog[0].Nlensing;
-  Naves = catalog[0].Naverage;
-  
-  /* set up pointers for linked list of measure, missing */
-  next_meas = build_measure_links (catalog[0].average, Naves, catalog[0].measure, Nmeas);
-  next_lens = build_lensing_links (catalog[0].average, Naves, catalog[0].lensing, Nlens);
-
-  catalog[0].sorted = TRUE;
-  catalog[0].measure = sort_measure (catalog[0].average, Naves, catalog[0].measure, Nmeas, next_meas);
-  catalog[0].lensing = sort_lensing (catalog[0].average, Naves, catalog[0].lensing, Nlens, next_lens);
-
-  MARKTIME ("  match time %9.4f sec for %7lld measures, %6lld average\n", dtime, (long long) Nmeas, (long long) Naves);
-
-  return;
-}
 
 void resort_catalog (Catalog *catalog) {
@@ -40 +10 @@
   resort_catalog_lensing (catalog);
   resort_catalog_starpar (catalog);
-
-  // XXX somehow need to set the value of lensing->measure offset here 
-
-
+  resort_catalog_galphot (catalog);
+  catalog[0].sorted = TRUE;
 }
 
@@ -88 +56 @@
       myAssert(average[averageSeq[i]].catID == measure[measureSeq[i]].catID, "object / detection mismatch");
 # if (1)
+      myAssert (average[averageSeq[i]].objID == measure[measureSeq[i]].objID, "object ID mismatch?");
+# endif
+# if (0)
       // myAssert(average[averageSeq[i]].objID == measure[measureSeq[i]].objID, "object / detection mismatch");
       // check if the objID is correct. if not, check if it is byte-swapped (this has happened) and repair if so.  
@@ -110 +81 @@
         }
       }
-# else
+# endif
+# if (0)
       // for reasons I do not understand, the mini dvodbs generated on stsci1X had a handful of detections with an inconsistency between averef and objID.  
       // this happened for 28 detections in the dbs on /data/stsci1?.0/eugene/dvo3pi.20130616, but not at all (as far as I know) in the rest of LAP DVO
@@ -196 +168 @@
   // MARKTIME("  match time %9.4f sec for %7lld measures, %6lld average\n", dtime, (long long) Nmeasure, (long long) Naverage);
 
-  catalog[0].sorted = TRUE;
-
   FREE (measureSeq);
   FREE (averageSeq);
@@ -245 +215 @@
       // earlier formats did not carry the objID or catID, so they are not available (we could assign on load, but we don't)
       myAssert(average[averageSeq[i]].catID == lensing[lensingSeq[i]].catID, "object / detection mismatch");
-# if (1)
       myAssert(average[averageSeq[i]].objID == lensing[lensingSeq[i]].objID, "object / detection mismatch");
-# else
-      // for reasons I do not understand, the mini dvodbs generated on stsci1X had a handful of detections with an inconsistency between averef and objID.  
-      // this happened for 28 detections in the dbs on /data/stsci1?.0/eugene/dvo3pi.20130616, but not at all (as far as I know) in the rest of LAP DVO
-      if (average[averageSeq[i]].objID != lensing[lensingSeq[i]].objID) {
-        fprintf (stderr, "R");
-        lensing[lensingSeq[i]].objID = average[averageSeq[i]].objID; // XXX I don't really like this...
-      }
-# endif
-    }
-  }
-  
-  // check that averageSeq is now in order
-  // for (i = 1; i < Nlensing; i++) {
-  //   if (averageSeq[i] < averageSeq[i-1]) {
-  //     fprintf (stderr, "%d ", (int) i);
-  //   }
-  // }
-  // fprintf (stderr, "\n");
-
+    }
+  }
+  
   SortAveMatch(lensingSeq, averageSeq, Nlensing);
-  // MARKTIME("sort : %f sec\n", dtime);
-
-  // check that averageSeq is now in order
-  // for (i = 1; i < Nlensing; i++) {
-  //   if (averageSeq[i] < averageSeq[i-1]) {
-  //     fprintf (stderr, "%d ", (int) i);
-  //   }
-  // }
-  // fprintf (stderr, "\n");
 
   // copy the lensing entries in the sorted order
@@ -283 +227 @@
     lensingTMP[i] = lensing[j];
   }
-  // MARKTIME("assign lensing : %f sec\n", dtime);
 
   // update the values of average.lensingOffset and average.Nlensing
@@ -302 +245 @@
     N++;
   }
-  // N++;
   average[currentAve].Nlensing = N;
-  // MARKTIME("update Nlensing : %f sec\n", dtime);
 
   int NlensingTotal = 0;
@@ -330 +271 @@
     fprintf (stderr, "ERROR: catalog %s has an invalid Nlensing\n", catalog[0].filename);
   }
-
-  // MARKTIME("  match time %9.4f sec for %7lld lensings, %6lld average\n", dtime, (long long) Nlensing, (long long) Naverage);
-
-  catalog[0].sorted = TRUE;
 
   FREE (lensingSeq);
@@ -382 +319 @@
       // earlier formats did not carry the objID or catID, so they are not available (we could assign on load, but we don't)
       myAssert(average[averageSeq[i]].catID == starpar[starparSeq[i]].catID, "object / detection mismatch");
-# if (1)
       myAssert(average[averageSeq[i]].objID == starpar[starparSeq[i]].objID, "object / detection mismatch");
-# else
-      // for reasons I do not understand, the mini dvodbs generated on stsci1X had a handful of detections with an inconsistency between averef and objID.  
-      // this happened for 28 detections in the dbs on /data/stsci1?.0/eugene/dvo3pi.20130616, but not at all (as far as I know) in the rest of LAP DVO
-      if (average[averageSeq[i]].objID != starpar[starparSeq[i]].objID) {
-        fprintf (stderr, "R");
-        starpar[starparSeq[i]].objID = average[averageSeq[i]].objID; // XXX I don't really like this...
-      }
-# endif
-    }
-  }
-  
-  // check that averageSeq is now in order
-  // for (i = 1; i < Nstarpar; i++) {
-  //   if (averageSeq[i] < averageSeq[i-1]) {
-  //     fprintf (stderr, "%d ", (int) i);
-  //   }
-  // }
-  // fprintf (stderr, "\n");
-
+    }
+  }
+  
   SortAveMatch(starparSeq, averageSeq, Nstarpar);
-  // MARKTIME("sort : %f sec\n", dtime);
-
-  // check that averageSeq is now in order
-  // for (i = 1; i < Nstarpar; i++) {
-  //   if (averageSeq[i] < averageSeq[i-1]) {
-  //     fprintf (stderr, "%d ", (int) i);
-  //   }
-  // }
-  // fprintf (stderr, "\n");
 
   // copy the starpar entries in the sorted order
@@ -420 +331 @@
     starparTMP[i] = starpar[j];
   }
-  // MARKTIME("assign starpar : %f sec\n", dtime);
 
   // update the values of average.starparOffset and average.Nstarpar
@@ -439 +349 @@
     N++;
   }
-  // N++;
   average[currentAve].Nstarpar = N;
-  // MARKTIME("update Nstarpar : %f sec\n", dtime);
 
   int NstarparTotal = 0;
@@ -468 +376 @@
   }
 
-  // MARKTIME("  match time %9.4f sec for %7lld starpars, %6lld average\n", dtime, (long long) Nstarpar, (long long) Naverage);
-
-  catalog[0].sorted = TRUE;
-
   FREE (starparSeq);
+  FREE (averageSeq);
+
+  return;
+}
+
+void resort_catalog_galphot (Catalog *catalog) {
+
+  off_t Naverage, Ngalphot;
+  GalPhot *galphot;
+  Average *average;
+  off_t i, j, N, currentAve;
+
+  off_t *galphotSeq = NULL;
+  off_t *averageSeq = NULL;
+  GalPhot *galphotTMP = NULL;
+
+  // struct timeval start, stop;
+  // gettimeofday (&start, NULL);
+
+  /* internal counters */
+  Ngalphot = catalog[0].Ngalphot;
+  Naverage = catalog[0].Naverage;
+
+  if (!Ngalphot) return;
+
+  galphot = catalog[0].galphot;
+  average = catalog[0].average;
+  
+  myAbort("relphot_catalog_galphot needs to handle average entries with Ngalphot = 0");
+
+  // we have a table of average objects and an unsorted table of measurements.  each measurement
+  // has a reference to the average object sequence (as well as an ID)
+  // galphot[i].averef -> average[averef]
+  // galphot[i].objID = average[averef].objID
+  // galphot[i].catID = average[averef].catID
+
+  // we want a sorted galphot array with all averef entries in sequence
+
+  ALLOCATE (galphotSeq, off_t,   Ngalphot);
+  ALLOCATE (averageSeq, off_t,   Ngalphot);
+
+  for (i = 0; i < Ngalphot; i++) {
+    galphotSeq[i] = i;
+    averageSeq[i] = galphot[i].averef;
+    
+    if (catalog[0].catformat >= DVO_FORMAT_PS1_V1) {
+      // earlier formats did not carry the objID or catID, so they are not available (we could assign on load, but we don't)
+      myAssert(average[averageSeq[i]].catID == galphot[galphotSeq[i]].catID, "object / detection mismatch");
+      myAssert(average[averageSeq[i]].objID == galphot[galphotSeq[i]].objID, "object / detection mismatch");
+    }
+  }
+  
+  SortAveMatch(galphotSeq, averageSeq, Ngalphot);
+
+  // copy the galphot entries in the sorted order
+  ALLOCATE (galphotTMP, GalPhot, Ngalphot);
+  for (i = 0; i < Ngalphot; i++) {
+    j = galphotSeq[i];
+    galphotTMP[i] = galphot[j];
+  }
+
+  // update the values of average.galphotOffset and average.Ngalphot
+  FREE(galphot);
+  catalog[0].galphot = galphotTMP;
+
+  N = 0;
+  currentAve = averageSeq[0];
+  average[currentAve].galphotOffset = 0;
+  for (i = 0; i < Ngalphot; i++) {
+    if (averageSeq[i] != currentAve) {
+      // we have hit the next entry in the list
+      average[currentAve].Ngalphot = N;
+      N = 0;
+      currentAve = averageSeq[i];
+      average[currentAve].galphotOffset = i;
+    }
+    N++;
+  }
+  average[currentAve].Ngalphot = N;
+
+  int NgalphotTotal = 0;
+  int galphotOffsetOK = TRUE;
+  for (i = 0; i < Naverage; i++) {
+    NgalphotTotal += catalog[0].average[i].Ngalphot;
+    if (VERBOSE && !(NgalphotTotal <= catalog[0].Ngalphot)) {
+      fprintf (stderr, "too few galphot: %d %d %d\n", (int) i, NgalphotTotal, (int) catalog[0].Ngalphot);
+    }
+    galphotOffsetOK &= (catalog[0].average[i].galphotOffset < catalog[0].Ngalphot);
+    if (VERBOSE && !(catalog[0].average[i].galphotOffset < catalog[0].Ngalphot)) {
+      fprintf (stderr, "offset too large: %d %d %d\n", (int) i, catalog[0].average[i].Ngalphot, (int) catalog[0].Ngalphot);
+    }
+    galphotOffsetOK &= (catalog[0].average[i].galphotOffset + catalog[0].average[i].Ngalphot <= catalog[0].Ngalphot);
+    if (VERBOSE && !(catalog[0].average[i].galphotOffset + catalog[0].average[i].Ngalphot <= catalog[0].Ngalphot)) {
+      fprintf (stderr, "orrset + Ngalphot too large: %d + %d > %d %d\n", (int) i, catalog[0].average[i].galphotOffset, catalog[0].average[i].Ngalphot, (int) catalog[0].Ngalphot);
+    }
+  }
+
+  if (!galphotOffsetOK) {
+    fprintf (stderr, "ERROR: catalog %s has an invalid galphotOffset\n", catalog[0].filename);
+  }
+
+  if (NgalphotTotal != catalog[0].Ngalphot) {
+    fprintf (stderr, "ERROR: catalog %s has an invalid Ngalphot\n", catalog[0].filename);
+  }
+
+  FREE (galphotSeq);
   FREE (averageSeq);
 

trunk/Ohana/src/addstar/src/resort_threaded.c

@@ -72 +72 @@
     // this is an overloaded value to mean 'force sort'
     if (threadData->forcesort) catalog.sorted = FALSE;
-
-    if (OLD_RESORT) { 
-      resort_catalog_old (&catalog);
-    } else {
-      resort_catalog (&catalog);
-    }
+    resort_catalog (&catalog);
 
     /* report total updated values */

trunk/Ohana/src/addstar/src/resort_unthreaded.c

@@ -42 +42 @@
     }
 
-    // this is an overloaded value to mean 'force sort'
+    // to force the sort, we need to set sorted 'false'
     if (ForceSort) catalog.sorted = FALSE;
-
-    if (OLD_RESORT) { 
-      resort_catalog_old (&catalog);
-    } else {
-      resort_catalog (&catalog);
-    }
+    resort_catalog (&catalog);
 
     // report total updated values