Changeset 33823

Timestamp:

Apr 19, 2012, 5:27:08 PM (14 years ago)

Author:

eugene

Message:

unify the 3 average magnitude calculation steps

Location:

branches/eam_branches/ipp-20120405/Ohana/src/relphot

Files:

• include/relphot.h (modified) (1 diff)
• src/ImageOps.c (modified) (3 diffs)
• src/StarOps.c (modified) (16 diffs)
• src/initialize.c (modified) (1 diff)
• src/reload_catalogs.c (modified) (1 diff)
• src/relphot_objects.c (modified) (2 diffs)
• src/setMrelFinal.c (modified) (3 diffs)

branches/eam_branches/ipp-20120405/Ohana/src/relphot/include/relphot.h

@@ -277 +277 @@
 int           setMmos             PROTO((Catalog *catalog, int Poor, FlatCorrectionTable *flatcorr));
 int           setMrel             PROTO((Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr));
-void          setMrelFinal        PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr));
+void          setMrelFinal        PROTO((Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage));
 int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr));
 int           setMave             PROTO((Catalog *catalog, int Ncatalog));

branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/ImageOps.c

@@ -41 +41 @@
 
 // relationships between the measure,catalog set and the images:
-static IDX_T       **MeasureToImage; // image index from measure,catalog   : MeasureToImage[cat][meas] = ImageIndex
-static IDX_T       **ImageToCatalog; // catalog for given measure on image : ImageCatalog[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex])
-static IDX_T       **ImageToMeasure; // measure for given measure on image : ImageMeasure[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex])
+static IDX_T       **MeasureToImage = NULL; // image index from measure,catalog   : MeasureToImage[cat][meas] = ImageIndex
+static IDX_T       **ImageToCatalog = NULL; // catalog for given measure on image : ImageCatalog[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex])
+static IDX_T       **ImageToMeasure = NULL; // measure for given measure on image : ImageMeasure[ImageIndex][i] = cat (i : 0 < NonImage[ImageIndex])
 
 // MeasureToImage was 'bin'
@@ -322 +322 @@
   off_t i;
 
+  if (!MeasureToImage) return -1;
+
   i = MeasureToImage[cat][meas];
   return (i);
@@ -354 +356 @@
   off_t i;
   short distance;
+
+  if (!MeasureToImage) return -1;
 
   i = MeasureToImage[cat][meas];

branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/StarOps.c

@@ -5 +5 @@
 typedef struct {
   int Nfew;
+  int Ncode;
   int Nsys;
   int Nbad;
@@ -12 +13 @@
   double *list;
   double *dlist;
+  double *aplist;
+  double *daplist;
 } SetMrelInfo;
 
@@ -25 +28 @@
 } ThreadInfo;
 
+// ** internal functions:
 void *setMrel_worker (void *data);
 int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr);
-int setMrel_catalog (Catalog *catalog, int Nc, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
+int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
+int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure);
 
 // we want to allocate the stats list,dlist arrays only once (or once per thread).
@@ -65 +70 @@
 void SetMrelInfoInit (SetMrelInfo *results, int allocLists) {
   results->Nfew  = 0;
+  results->Ncode  = 0;
   results->Nsys  = 0;
   results->Nbad  = 0;
@@ -83 +89 @@
 void SetMrelInfoAccum (SetMrelInfo *summary, SetMrelInfo *results) {
   summary->Nfew  += results->Nfew ;
+  summary->Ncode  += results->Ncode ;
   summary->Nsys  += results->Nsys ;
   summary->Nbad  += results->Nbad ;
@@ -110 +117 @@
 }
 
+// setMrel and setMrelOutput are extremely similar, but have slightly different implications:
+// * setMrel uses the internal Tiny structures only
+// * setMrelOutput skips stars for which there are too few good measurements 
+// * setMrelOutput is meant to be called repeatedly, relaxing the criteria for 'good' on each pass
+// * setMrelOutput updates 2MASS average flags
+// * setMrelOutput updates average EXT flags (PS1 and 2MASS)
+
 int setMrel (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
 
@@ -126 +140 @@
 
   for (i = 0; i < Ncatalog; i++) {
-    setMrel_catalog (catalog, i, flatcorr, &results, Nsecfilt);
+    // pass == -1 for anything other than the final pass
+    setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt);
     SetMrelInfoAccum (&summary, &results);
   }
-  fprintf (stderr, "%d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
+  fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
 
   SetMrelInfoFree (&results);
@@ -135 +150 @@
   return (TRUE);
 }
+
+int setMrelOutput (Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr) {
+
+  int i;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  SetMrelInfo summary, results;
+  SetMrelInfoInit (&summary, FALSE);
+  SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist
+  ALLOCATE (results.aplist, double, Nmax);
+  ALLOCATE (results.daplist, double, Nmax);
+
+  for (i = 0; i < Ncatalog; i++) {
+    setMrel_catalog  (catalog, i, pass, flatcorr, &results, Nsecfilt); // XXX add arguments as needed for options
+    SetMrelInfoAccum (&summary, &results);
+  }
+  fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
+
+  SetMrelInfoFree (&results);
+  free (results.aplist);
+  free (results.daplist);
+  return (TRUE);
+}
+
+int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
+
+  int i;
+
+  SetMrelInfo summary;
+  SetMrelInfoInit (&summary, FALSE);
+
+  pthread_attr_t attr;
+  pthread_attr_init (&attr);
+  pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED);
+  
+  pthread_t *threads;
+  ALLOCATE (threads, pthread_t, NTHREADS);
+
+  ThreadInfo *threadinfo;
+  ALLOCATE (threadinfo, ThreadInfo, NTHREADS);
+
+  // launch N worker threads
+  for (i = 0; i < NTHREADS; i++) {
+    threadinfo[i].entry = i;
+    threadinfo[i].state = THREAD_RUN;
+    threadinfo[i].catalog  =  catalog;
+    threadinfo[i].Ncatalog = Ncatalog;
+    threadinfo[i].flatcorr = flatcorr;
+    SetMrelInfoInit (&threadinfo[i].summary, FALSE);
+    pthread_create (&threads[i], NULL, setMrel_worker, &threadinfo[i]);
+  }
+  pthread_attr_destroy (&attr);
+
+  // wait until all threads have finished
+  while (1) {
+    int allDone = TRUE;
+    for (i = 0; i < NTHREADS; i++) {
+      if (threadinfo[i].state == THREAD_RUN) allDone = FALSE;
+    }
+    if (allDone) {
+      break;
+    }
+    usleep (500000);
+  }
+
+  // all threads are done, free the threads array and grab the info
+  free (threads);
+  
+  // report stats & summary from the threads
+  for (i = 0; i < NTHREADS; i++) {
+    fprintf (stderr, "setMrel thread %d : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", 
+             i, 
+             threadinfo[i].summary.Ncode, 
+             threadinfo[i].summary.Nfew, 
+             threadinfo[i].summary.Nbad, 
+             threadinfo[i].summary.Ncal, 
+             threadinfo[i].summary.Nmos, 
+             threadinfo[i].summary.Ngrid, 
+             threadinfo[i].summary.Nsys);
+    SetMrelInfoAccum (&summary, &threadinfo[i].summary);
+  }
+  fprintf (stderr, "total : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
+  free (threadinfo);
+
+  return TRUE;
+}
+
+void *setMrel_worker (void *data) {
+
+  ThreadInfo *threadinfo = data;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  SetMrelInfo results;
+  SetMrelInfoInit (&results, TRUE); // allocate list, dlist arrays here
+
+  while (1) {
+
+    off_t i = getNextCatalogForThread(threadinfo->Ncatalog);
+    if (i == -1) {
+      threadinfo->state = THREAD_DONE;
+      return NULL;
+    }
+
+    Catalog *catalog = threadinfo->catalog;
+    FlatCorrectionTable *flatcorr = threadinfo->flatcorr;
+
+    // pass == -1 for anything other than the final pass
+    setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt);
+    SetMrelInfoAccum (&threadinfo->summary, &results);
+  }
+
+  SetMrelInfoFree (&results);
+  return NULL;
+}
+
+# define SKIP_THIS_MEAS(REASON) {                               \
+    catalog[Nc].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;     \
+    if (catalog[Nc].measure) {                                  \
+      catalog[Nc].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;   \
+    }                                                           \
+    results->REASON ++;                                         \
+    continue; }
+
+
+// setMrel_catalog is used in 3 different contexts:
+
+// * during the main relphot iterations:
+// ** operations only apply to measureT / averageT
+// ** we are applying the analysis to the bright subset catalog
+// ** we skip any stars found to be bad (STAR_BAD)
+
+// * during the final pass
+// ** operations are applied to all objects
+// ** we have special tests for PS1 extended data, 2MASS extended & good data, & synthetic photometry
+
+// * during the stand-alone '-averages' mode (relphot_objects)
+// ** no image data is loaded
+// ** we only use the provided calibration (measure.Mcal)
+// ** we skip outlier rejection of measurements
 
 // set the Mrel values for the specified catalog
 // NOTE: here 'catalog' is a pointer to a specific catalog, not the root of the array
-int setMrel_catalog (Catalog *catalog, int Nc, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
+int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
 
   off_t j, k, m;
   int N;
   float Msys, Mcal, Mmos, Mgrid;
-  StatType stats;
-
-  double *list  = results->list;
-  double *dlist = results->dlist;
+  StatType stats, apstats;
+
+  double *list    = results->list;
+  double *dlist   = results->dlist;
+  double *aplist  = results->aplist;
+  double *daplist = results->daplist;
 
   SetMrelInfoInit (results, FALSE); // do not allocate list,dlist arrays
+
+  int isSetMrelFinal = (pass >= 0);
 
   for (j = 0; j < catalog[Nc].Naverage; j++) {
     // XXX accumulate all secfilt values in a single pass?
 
-    int minUbercalDist = 1000;
+    // option for a test print
+    if (FALSE && (catalog[Nc].average[j].objID == 0x46a4) && (catalog[Nc].average[j].catID == 0xf40e)) {
+      fprintf (stderr, "test obj\n");
+      print_measure_set (&catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], catalog[Nc].measure);
+    }
+
+    int GoodPS1 = FALSE;
+    int Good2MASS = FALSE;
+    int Galaxy2MASS = FALSE;
+
+    int NextPS1 = 0;
+    int NpsfPS1 = 0;
 
     int Ns;
@@ -163 +334 @@
       /* calculate the average mag in this SEC photcode for a single star */
 
+      /* star/photcodes already calibrated */
+      if ( isSetMrelFinal && catalog[Nc].found[Nsecfilt*j+Nsec]) continue;  
+      
       // skip bad stars
-      if (catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
-
+      if (!isSetMrelFinal && (catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD)) continue;
+
+      int Ncode = 0;
+      int Next = 0;
+      int haveSynth = FALSE;
+      int haveUbercal = FALSE;
+
+      int minUbercalDist = 1000;
+    
       N = 0;
       m = catalog[Nc].averageT[j].measureOffset;
@@ -174 +355 @@
         if (!code) continue;
         if (code->equiv != thisCode) { continue; }
-
-        if (catalog[Nc].measureT[m].dbFlags & MEAS_BAD) { results->Nbad ++; continue; }
+        Ncode ++;
+
+        if (catalog[Nc].measureT[m].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); 
 
         if (getImageEntry (m, Nc) < 0) {
@@ -185 +367 @@
         } else {
           Mcal  = getMcal  (m, Nc, flatcorr, catalog);
-          if (isnan(Mcal))  { results->Ncal ++; continue; }
+          if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
           Mmos  = getMmos  (m, Nc);
-          if (isnan(Mmos))  { results->Nmos ++; continue; }
+          if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
           Mgrid = getMgrid (m, Nc);
-          if (isnan(Mgrid)) { results->Ngrid++; continue; }
+          if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
         }
 
         Msys = PhotSysTiny (&catalog[Nc].measureT[m], &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt]);
-        if (isnan(Msys)) { results->Nsys++; continue; }
+        if (isnan(Msys)) SKIP_THIS_MEAS(Nsys);
+
         list[N] = Msys - Mcal - Mmos - Mgrid;
 
         int myUbercalDist = getUbercalDist(m,Nc);
         minUbercalDist = MIN(minUbercalDist, myUbercalDist);
+
+        if (isSetMrelFinal) {
+          float Map = PhotAper (&catalog[Nc].measure[m]);
+          aplist[N] = Map - Mcal - Mmos - Mgrid;
+
+          // count the extended detections
+          if ((catalog[Nc].measure[m].photcode >= 10000) && (catalog[Nc].measure[m].photcode <= 10500)) {
+            if (!isnan(catalog[Nc].measure[m].Map)) {
+              if (catalog[Nc].measure[m].M - catalog[Nc].measure[m].Map > 0.5) {
+                Next ++;
+                NextPS1 ++;
+              } else {
+                NpsfPS1 ++;
+              }
+            }
+          }
+          // count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013)
+          if ((catalog[Nc].measure[m].photcode >= 2011) && (catalog[Nc].measure[m].photcode <= 2013)) {
+            if (catalog[Nc].measure[m].photFlags & 0x00c00000) {
+              Next ++;
+              Galaxy2MASS = TRUE;
+            }
+            if (pass == 0) {
+              if (catalog[Nc].measure[m].photFlags & 0x00000007) {
+                Good2MASS = TRUE;
+              } else {
+                // detections without one of these bits should only be used in PASS_1
+                SKIP_THIS_MEAS(Nbad);
+                continue;
+              }
+            }
+          }
+          // ignore SYNTH photocdes until PASS == 4 (where we also accept saturated stars)
+          if ((catalog[Nc].measure[m].photcode >= 3001) && (catalog[Nc].measure[m].photcode <= 3005)) {
+            if (pass < 4) {
+              SKIP_THIS_MEAS(Nbad);
+              continue;
+            }
+            haveSynth = TRUE;
+          }
+        }
 
         // dlist gives the error, which is used as the weight in WT_MEAN.
@@ -225 +449 @@
       }
 
-      // when performing the grid analysis, STAR_TOOFEW will be set to 1;
-      if (N <= STAR_TOOFEW) { /* too few measurements */
+      int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1;
+
+      // when performing the grid analysis, STAR_TOOFEW should be set to 1;
+      if (N < Nminmeas) { /* too few measurements */
         // fprintf (f, "%10.6f %10.6f %d %d %d\n", catalog[Nc].averageT[j].R, catalog[Nc].averageT[j].D, catalog[Nc].measureT[catalog[Nc].averageT[j].measureOffset].imageID, N, STAR_TOOFEW); 
         catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW;
-        results->Nfew ++;
+        if (Ncode == 0) {
+          results->Ncode ++;
+        } else {
+          results->Nfew ++;
+        }
+        continue;
       } else {
         catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_FEW;
@@ -235 +466 @@
 
       liststats (list, dlist, N, &stats);
-        
+
       catalog[Nc].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
-      catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM = stats.sigma;
+      catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM = stats.error;
       catalog[Nc].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq + 1e-4) : NAN_S_SHORT;
-        
-      catalog[Nc].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist;
-    }
-  }
-
-  // values which need to be passed out: Nfew, Nbad, Ncal, Nmos, Ngrid, Nsys
+
+      // when running -averages, we have no information about the images, so we cannot set this
+      if (minUbercalDist > -1) {
+        catalog[Nc].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist;
+      }
+
+      if (isSetMrelFinal) {
+        liststats (aplist, daplist, N, &apstats);
+        catalog[Nc].found[Nsecfilt*j+Nsec] = TRUE;
+
+        catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mstdev = 1000.0*stats.sigma; // Mstdev is in millimags (not enough space for more precision)
+        catalog[Nc].secfilt[Nsecfilt*j+Nsec].Ncode = Ncode;
+        catalog[Nc].secfilt[Nsecfilt*j+Nsec].Nused = stats.Nmeas;
+
+        catalog[Nc].secfilt[Nsecfilt*j+Nsec].M_80 = 1000 * stats.Upper80;
+        catalog[Nc].secfilt[Nsecfilt*j+Nsec].M_20 = 1000 * stats.Lower20;
+
+        // NOTE: for 2MASS measurements, Next should be 1, as should N
+        if ((Next > 0) && (Next > 0.5*N)) {
+          catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_OBJ_EXT;
+        }
+
+        switch (pass) {
+          case 0:
+            catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_0;
+            GoodPS1 = TRUE;
+            break;
+          case 1:
+            catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1;
+            GoodPS1 = TRUE;
+            break;
+          case 2:
+            catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2;
+            GoodPS1 = TRUE;
+            break;
+          case 3:
+            catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_3;
+            break;
+          case 4:
+            catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_4;
+            break;
+        }
+        if (haveSynth) {
+          catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_SYNTH;
+        }       
+        if (haveUbercal) {
+          catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_UBERCAL;
+        }       
+      }
+    }
+
+    if (isSetMrelFinal) {
+      DVOAverageFlags flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT;
+
+      // we attempt to set a few flags here; reset those bits before trying:
+      catalog[Nc].average[j].flags &= ~flagBits;
+
+      if (NextPS1 && (NextPS1 > NpsfPS1)) {
+        catalog[Nc].average[j].flags |= ID_OBJ_EXT;
+      }
+      if (GoodPS1) {
+        catalog[Nc].average[j].flags |= ID_OBJ_GOOD;
+      }
+      if (Galaxy2MASS) {
+        catalog[Nc].average[j].flags |= ID_OBJ_EXT_ALT;
+      }
+      if (Good2MASS) {
+        catalog[Nc].average[j].flags |= ID_OBJ_GOOD_ALT;
+      }
+    }
+  }
   return (TRUE);
-}
-
-int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
-
-  int i;
-
-  SetMrelInfo summary;
-  SetMrelInfoInit (&summary, FALSE);
-
-  pthread_attr_t attr;
-  pthread_attr_init (&attr);
-  pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED);
-  
-  pthread_t *threads;
-  ALLOCATE (threads, pthread_t, NTHREADS);
-
-  ThreadInfo *threadinfo;
-  ALLOCATE (threadinfo, ThreadInfo, NTHREADS);
-
-  // launch N worker threads
-  for (i = 0; i < NTHREADS; i++) {
-    threadinfo[i].entry = i;
-    threadinfo[i].state = THREAD_RUN;
-    threadinfo[i].catalog  =  catalog;
-    threadinfo[i].Ncatalog = Ncatalog;
-    threadinfo[i].flatcorr = flatcorr;
-    SetMrelInfoInit (&threadinfo[i].summary, FALSE);
-    pthread_create (&threads[i], NULL, setMrel_worker, &threadinfo[i]);
-  }
-  pthread_attr_destroy (&attr);
-
-  // wait until all threads have finished
-  while (1) {
-    int allDone = TRUE;
-    for (i = 0; i < NTHREADS; i++) {
-      if (threadinfo[i].state == THREAD_RUN) allDone = FALSE;
-    }
-    if (allDone) {
-      break;
-    }
-    usleep (500000);
-  }
-
-  // all threads are done, free the threads array and grab the info
-  free (threads);
-  
-  // report stats & summary from the threads
-  for (i = 0; i < NTHREADS; i++) {
-    fprintf (stderr, "setMrel thread %d : %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", 
-             i, 
-             threadinfo[i].summary.Nfew, 
-             threadinfo[i].summary.Nbad, 
-             threadinfo[i].summary.Ncal, 
-             threadinfo[i].summary.Nmos, 
-             threadinfo[i].summary.Ngrid, 
-             threadinfo[i].summary.Nsys);
-    SetMrelInfoAccum (&summary, &threadinfo[i].summary);
-  }
-  fprintf (stderr, "total : %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
-  free (threadinfo);
-
-  return TRUE;
-}
-
-void *setMrel_worker (void *data) {
-
-  ThreadInfo *threadinfo = data;
-
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-
-  SetMrelInfo results;
-  SetMrelInfoInit (&results, TRUE); // allocate list, dlist arrays here
-
-  while (1) {
-
-    off_t i = getNextCatalogForThread(threadinfo->Ncatalog);
-    if (i == -1) {
-      threadinfo->state = THREAD_DONE;
-      return NULL;
-    }
-
-    Catalog *catalog = threadinfo->catalog;
-    FlatCorrectionTable *flatcorr = threadinfo->flatcorr;
-
-    setMrel_catalog (catalog, i, flatcorr, &results, Nsecfilt);
-    SetMrelInfoAccum (&threadinfo->summary, &results);
-  }
-
-  SetMrelInfoFree (&results);
-  return NULL;
 }
 
@@ -357 +563 @@
 }
 
-# define MARK_SKIP_MEAS \
-  catalog[i].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM; \
-  catalog[i].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
-
-// setMrel and setMrelOutput are extremely similar, but have slightly different implications:
-// * setMrel uses the internal Tiny structures only
-// * setMrelOutput skips stars for which there are too few good measurements 
-// * setMrelOutput is meant to be called repeatedly, relaxing the criteria for 'good' on each pass
-
-// setMave is also similar to the above.  but:
-// * setMave is called by relphot -update-objects
-// * setMave excludes all detections with (PSF_QF < 0.85), setMrelOutput allows these for PASS > 2
-// * setMave updates 2MASS average flags
-// * setMave updates average EXT flags (PS1 and 2MASS)
-int setMrelOutput (Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr) {
-
-  int i, N;
-  off_t j, k, m, Nmax;
-  float Msys, Mcal, Mmos, Mgrid, Map;
-  double *list, *dlist, *aplist, *daplist;
-  StatType stats, apstats;
-  int Nsec;
-
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-
-  /* Nmeasure is now different, need to reallocate */
-  Nmax = 0;
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-      Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure);
-    }
-  }
-  ALLOCATE (list, double, MAX (1, Nmax));
-  ALLOCATE (dlist, double, MAX (1, Nmax));
-
-  ALLOCATE (aplist, double, MAX (1, Nmax));
-  ALLOCATE (daplist, double, MAX (1, Nmax));
-
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-
-      if (FALSE && (catalog[i].average[j].objID == 0x46a4) && (catalog[i].average[j].catID == 0xf40e)) {
-        fprintf (stderr, "test obj\n");
-        print_measure_set (&catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt], catalog[i].measure);
-      }
-
-      int minUbercalDist = 1000;
-
-      int Ns;
-      for (Ns = 0; Ns < Nphotcodes; Ns++) {
-
-        int thisCode = photcodes[Ns][0].code;
-        Nsec = GetPhotcodeNsec(thisCode);
-
-        /* star/photcodes already calibrated */
-        if (catalog[i].found[Nsecfilt*j+Nsec]) continue;  
-
-        int Ncode = 0;
-        int Next = 0;
-        int haveSynth = FALSE;
-        int haveUbercal = FALSE;
-        int isRefPhot = FALSE;
-
-        N = 0;
-        m = catalog[i].averageT[j].measureOffset;
-        for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
-
-          // skip measurements that do not match the current photcode
-          PhotCode *code = GetPhotcodebyCode (catalog[i].measureT[m].photcode);
-          if (!code) continue;
-          if (code->equiv != thisCode) { continue; }
-          Ncode ++;
-
-          if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
-
-          if (getImageEntry (m, i) < 0) {
-          XXXX : likely problem -- compare with setMrel_catalog;
-            // these detetions have no image (eg, ref values such as 2MASS)
-            Mcal = Mmos = Mgrid = 0;
-            isRefPhot = TRUE;
-          } else {
-            Mcal  = getMcal  (m, i, flatcorr, catalog);
-            if (isnan(Mcal)) {
-              MARK_SKIP_MEAS;
-              continue;
-            }
-            Mmos  = getMmos  (m, i);
-            if (isnan(Mmos)) {
-              MARK_SKIP_MEAS;
-              continue;
-            }
-            Mgrid = getMgrid (m, i);
-            if (isnan(Mgrid)) {
-              MARK_SKIP_MEAS;
-              continue;
-            }
-          }
-
-          Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]);
-          if (isnan(Msys)) { 
-            if (!isRefPhot) {
-              MARK_SKIP_MEAS;
-            }
-            continue; 
-          }
-          list[N] = Msys - Mcal - Mmos - Mgrid;
-
-          Map = PhotAper (&catalog[i].measure[m]);
-          aplist[N] = Map - Mcal - Mmos - Mgrid;
-
-          int myUbercalDist = getUbercalDist(m,i);
-          minUbercalDist = MIN(minUbercalDist, myUbercalDist);
-
-          // count the extended detections
-          if ((catalog[i].measure[m].photcode >= 10000) && (catalog[i].measure[m].photcode <= 10500)) {
-            if (!isnan(catalog[i].measure[m].Map)) {
-              if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) {
-                Next ++;
-              }
-            }
-          }
-          // count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013)
-          if ((catalog[i].measure[m].photcode >= 2011) && (catalog[i].measure[m].photcode <= 2013)) {
-            if (catalog[i].measure[m].photFlags & 0x00c00000) {
-              Next ++;
-            }
-            if ((pass == 0) && !(catalog[i].measure[m].photFlags & 0x00000007)) {
-              // detections without one of these bits should only be used in PASS_1
-              MARK_SKIP_MEAS;
-              continue;
-            }
-          }
-          // ignore SYNTH photocdes until PASS == 4 (where we also accept saturated stars)
-          if ((catalog[i].measure[m].photcode >= 3001) && (catalog[i].measure[m].photcode <= 3005)) {
-            if (pass < 4) {
-              MARK_SKIP_MEAS;
-              continue;
-            }
-            haveSynth = TRUE;
-          }
-
-          // dlist gives the error, which is used as the weight in WT_MEAN.
-          // we can modify the resulting weight in a few ways:
-          // 1) MIN_ERROR guarantees a floor
-          // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
-          // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
-          // 4) some reference photcode of some kind can be specified as fixed and have a high weight
-          dlist[N] = MAX (hypot(catalog[i].measureT[m].dM, code->photomErrSys), MIN_ERROR);
-
-          // up-weight the ubercal values (or convergence can take a long time...) (XXX make this optional?)
-          if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
-            dlist[N] = MAX (0.1*catalog[i].measureT[m].dM, MIN_ERROR);
-            haveUbercal = TRUE;
-          }
-
-          // tie down reference photometry if the -refcode (code) option is selected
-          // eg, -refcode g_SDSS
-          // this probably makes no sense in the context of multifilter analysis
-          if (refPhotcode) {
-            if (code->code == refPhotcode->code) {
-              // tiny error -> large weight
-              // dlist[N] = MAX (0.01*catalog[i].measureT[m].dM, MIN_ERROR);
-              dlist[N] = 0.0001;
-            }
-          }
-          N++;
-        }
-        if (N < 1) continue;
-
-        for (m = 0; m < N; m++) {
-          daplist[m] = dlist[m];
-        }
-
-        // XXX force WT_MEAN or MEAN here?
-        liststats (list, dlist, N, &stats);
-        liststats (aplist, daplist, N, &apstats);
-        catalog[i].found[Nsecfilt*j+Nsec] = TRUE;
-        switch (pass) {
-          case 0:
-            catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_0;
-            break;
-          case 1:
-            catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1;
-            break;
-          case 2:
-            catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2;
-            break;
-          case 3:
-            catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_3;
-            break;
-          case 4:
-            catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_4;
-            break;
-        }
-
-        if (haveSynth) {
-          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_SYNTH;
-        }       
-        if (haveUbercal) {
-          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_UBERCAL;
-        }       
-
-        /* use sigma or error in dM for output? */
-        catalog[i].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
-        catalog[i].secfilt[Nsecfilt*j+Nsec].Map = apstats.mean;
-        catalog[i].secfilt[Nsecfilt*j+Nsec].dM = stats.error;
-        catalog[i].secfilt[Nsecfilt*j+Nsec].Mstdev = 1000.0*stats.sigma; // Mstdev is in millimags (not enough space for more precision)
-        catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
-        catalog[i].secfilt[Nsecfilt*j+Nsec].Ncode = Ncode;
-        catalog[i].secfilt[Nsecfilt*j+Nsec].Nused = stats.Nmeas;
-
-        catalog[i].secfilt[Nsecfilt*j+Nsec].M_80 = 1000 * stats.Upper80;
-        catalog[i].secfilt[Nsecfilt*j+Nsec].M_20 = 1000 * stats.Lower20;
-        catalog[i].secfilt[Nsecfilt*j+Nsec].ubercalDist = minUbercalDist;
-
-        if ((Next > 0) && (Next > 0.5*N)) {
-          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_OBJ_EXT;
-        }
-      }
-    }
-  }
-
-  free (list);
-  free (dlist);
-
-  free (aplist);
-  free (daplist);
-  return (TRUE);
-}
+# if (0)
 
 /* grab Nsec for named photcode */
@@ -776 +754 @@
   return (TRUE);
 }
+# endif
 
 /* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_PHOTOM_NOCAL */
@@ -800 +779 @@
         if (isnan(Mgrid)) continue;
 
-        // note that measurements for which the image is not selected will not up modified
+        // note that measurements for which the image is not selected will not be modified
         // (that is a good thing! -- we keep a prior calibration)
 

branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/initialize.c

@@ -12 +12 @@
     PhotcodeList = strcreate (argv[1]);
     photcodes = ParsePhotcodeList (PhotcodeList, &Nphotcodes, TRUE); // require SEC photcodes
-  }
+  } else {
+    char tmpline1[256];
+    int Ns;
+    Nphotcodes = GetPhotcodeNsecfilt ();
+    ALLOCATE (photcodes, PhotCode *, Nphotcodes);
+    ALLOCATE (PhotcodeList, char, 256);
+    for (Ns = 0; Ns < Nphotcodes; Ns++) {
+      photcodes[Ns] = GetPhotcodebyNsec (Ns);
+      if (Ns > 0) {
+        snprintf (tmpline1, 256, "%s,%s", PhotcodeList, photcodes[Ns][0].name);
+      } else {
+        snprintf (tmpline1, 256, "%s", photcodes[Ns][0].name);
+      }
+      strcpy (PhotcodeList, tmpline1);
+    }
+  }    
+    
   if (USE_GRID && (Nphotcodes > 1)) {
     fprintf (stderr, "grid correction analysis currently can only operate on a single photcode\n");

branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/reload_catalogs.c

@@ -86 +86 @@
     TIMESTAMP(time4);
 
-    setMrelFinal (&catalog, flatcorr);
+    initMrel (&catalog, 1);
+    setMrelFinal (&catalog, flatcorr, FALSE);
     TIMESTAMP(time5);
 

branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/relphot_objects.c

@@ -90 +90 @@
     }
 
-    setMave (&catalog, 1);
+    populate_tiny_values(&catalog, DVO_TV_MEASURE | DVO_TV_AVERAGE);
+
+    initMrel (&catalog, 1);
+    setMrelFinal (&catalog, NULL, TRUE);
 
     if (!UPDATE) {
       dvo_catalog_unlock (&catalog);
+      free_tiny_values(&catalog);
       dvo_catalog_free (&catalog);
       continue;
@@ -105 +109 @@
     dvo_catalog_save (&catalog, VERBOSE); 
     dvo_catalog_unlock (&catalog);
+    free_tiny_values(&catalog);
     dvo_catalog_free (&catalog);
   }

branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/setMrelFinal.c

@@ -5 +5 @@
 // output dbFlags values
 
-void setMrelFinal (Catalog *catalog, FlatCorrectionTable *flatcorr) {
+void setMrelFinal (Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage) {
 
   off_t i;
@@ -94 +94 @@
   }
 
-  clean_measures (catalog, 1, TRUE, flatcorr);    /* mark outliers ID_MEAS_POOR_PHOTOM */
+  // XXX make this optional? (do not clean for -averages?)
+  if (!simpleAverage) clean_measures (catalog, 1, TRUE, flatcorr);    /* mark outliers ID_MEAS_POOR_PHOTOM */
   for (i = 0; i < 5; i++) {
     skip_measurements (catalog, i, flatcorr);       /* set ID_MEAS_SKIP for measures to be skipped */
     setMrelOutput  (catalog, 1, i, flatcorr);    /* set Mrel using remaining measures */
   }
-  setMcalOutput (catalog, 1, flatcorr);
+  if (!simpleAverage) setMcalOutput (catalog, 1, flatcorr);
 
   /* clear ID_STAR_POOR, ID_STAR_FEW values before writing ??? */
@@ -258 +259 @@
     }
   }
-  if (VERBOSE2) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n",
+  if (VERBOSE) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n",
                         pass, Ntot, Ntry, Nskip, Nkeep);
 }