Changeset 34041 for branches/meh_branches/ppstack_test/Ohana/src/relphot/src/StarOps.c

Timestamp:

Jun 19, 2012, 5:24:19 PM (14 years ago)

Author:

mhuber

Message:

merging latest r34040 trunk changes to branch

Location:

branches/meh_branches/ppstack_test

Files:

• . (modified) (1 prop)
• Ohana/src/relphot/src/StarOps.c (modified) (26 diffs)

branches/meh_branches/ppstack_test/Ohana/src/relphot/src/StarOps.c

@@ -2 +2 @@
 
 static int Nmax;
-static double *list;
-static double *dlist;
-
-// When we rationalize the images/mosaics, we are driving the negative cloud images back
-// to 0.0.  At the same time, we make a guess to the effective impact on all other images,
-// driven by the coupling of common stars.  This array carries the impact of those offsets
-// on each star
-static double *Moffset;
-
+
+typedef struct {
+  int Nfew;
+  int Ncode;
+  int Nsys;
+  int Nbad;
+  int Ncal;
+  int Nmos;
+  int Ngrid;
+  double *list;
+  double *dlist;
+  double *wlist;
+  double *aplist;
+  double *daplist;
+} SetMrelInfo;
+
+enum {THREAD_RUN, THREAD_DONE};
+
+typedef struct {
+  int entry;
+  int state;
+  Catalog *catalog;
+  int Ncatalog;
+  FlatCorrectionTable *flatcorr;
+  SetMrelInfo summary;
+} ThreadInfo;
+
+// ** internal functions:
+void *setMrel_worker (void *data);
+int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr);
+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).
+// this function finds the largest array so we can allocate that max size when needed
 void initMrel (Catalog *catalog, int Ncatalog) {
 
@@ -21 +47 @@
     }
   }
-
-  ALLOCATE (list,    double, MAX (1, Nmax));
-  ALLOCATE (dlist,   double, MAX (1, Nmax));
-  ALLOCATE (Moffset, double, MAX (1, Nmax));
 }  
 
@@ -47 +69 @@
 }
 
-int setMrel (Catalog *catalog, int Ncatalog) {
+void SetMrelInfoInit (SetMrelInfo *results, int allocLists) {
+  results->Nfew  = 0;
+  results->Ncode  = 0;
+  results->Nsys  = 0;
+  results->Nbad  = 0;
+  results->Ncal  = 0;
+  results->Nmos  = 0;
+  results->Ngrid = 0;
+  if (allocLists) {
+    ALLOCATE (results->list,  double, Nmax);
+    ALLOCATE (results->dlist, double, Nmax);
+    ALLOCATE (results->wlist, double, Nmax);
+  }
+}
+
+void SetMrelInfoFree (SetMrelInfo *results) {
+  free (results->list);
+  free (results->dlist);
+  free (results->wlist);
+}
+
+void SetMrelInfoAccum (SetMrelInfo *summary, SetMrelInfo *results) {
+  summary->Nfew  += results->Nfew ;
+  summary->Ncode  += results->Ncode ;
+  summary->Nsys  += results->Nsys ;
+  summary->Nbad  += results->Nbad ;
+  summary->Ncal  += results->Ncal ;
+  summary->Nmos  += results->Nmos ;
+  summary->Ngrid += results->Ngrid;
+}
+
+// mutex to lock setMrel_worker operations 
+static pthread_mutex_t setMrel_mutex = PTHREAD_MUTEX_INITIALIZER;
+static int nextCatalog = 0;
+
+// we have an array of catalogs (catalog, Ncatalog).  we need to hand out catalogs one at a time to
+// the worker threads as they need
+off_t getNextCatalogForThread (int Ncatalog) {
+
+  pthread_mutex_lock (&setMrel_mutex);
+  if (nextCatalog >= Ncatalog) {
+    pthread_mutex_unlock (&setMrel_mutex);
+    return (-1);
+  }
+  int thisCatalog = nextCatalog;
+  nextCatalog ++;
+
+  pthread_mutex_unlock (&setMrel_mutex);
+  return (thisCatalog);
+}
+
+// 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) {
+
+  int i;
+
+  if (NTHREADS) {
+    int status = setMrel_threaded (catalog, Ncatalog, flatcorr);
+    return status;
+  }
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  SetMrelInfo summary, results;
+  SetMrelInfoInit (&summary, FALSE);
+  SetMrelInfoInit (&results, TRUE);
+
+  for (i = 0; i < Ncatalog; i++) {
+    // pass == -1 for anything other than the final pass
+    setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt);
+    SetMrelInfoAccum (&summary, &results);
+  }
+  if (VERBOSE2) 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);
+
+  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,wlist
+  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);
+  }
+  if (VERBOSE2) 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);
+
+  // each time this function is called, we cycle through the available catalogs.
+  // make sure we start at 0
+  nextCatalog = 0;;
+
+  // 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++) {
+    if (VERBOSE2) 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);
+  }
+  if (VERBOSE2) 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, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {
 
   off_t j, k, m;
-  int i, N, Nfew, Nsys, Nbad, Ncal, Nmos, Ngrid;
+  int N;
   float Msys, Mcal, Mmos, Mgrid;
-  StatType stats;
-
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-  Nfew = Nsys = Nbad = Ncal = Nmos = Ngrid = 0;
-
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-
-      int Ns;
-      for (Ns = 0; Ns < Nphotcodes; Ns++) {
-
-        int thisCode = photcodes[Ns][0].code;
-        int Nsec = GetPhotcodeNsec(thisCode);
-
-        /* calculate the average value for a single star */
-
-        // skip bad stars
-        if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
-        m = catalog[i].averageT[j].measureOffset;
-
-        N = 0;
-        for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
-          // skip measurements that do not match the current photcode
-          int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
-          if (ecode != thisCode) { continue; }
-
-          if (catalog[i].measureT[m].dbFlags & MEAS_BAD) {
-            Nbad ++;
-            continue;
-          }
-          // XXX allow REF stars (no Image Entry) to be included in the calculation this
-          // should be optionally set, and should allow for REF stars to be downweighted by
-          // more than their reported errors.  how such info is carried is unclear...
-          if (getImageEntry (m, i) < 0) {
-            Mcal = Mmos = Mgrid = 0;
-          } else {
-            Mcal  = getMcal  (m, i);
-            if (isnan(Mcal)) { 
-              Ncal ++;
-              continue;
-            }
-            Mmos  = getMmos  (m, i);
-            if (isnan(Mmos)) {
-              Nmos ++;
-              continue;
-            }
-            Mgrid = getMgrid (m, i);
-            if (isnan(Mgrid)) {
-              Ngrid++;
-              continue;
-            }
-          }
-
-          Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]);
-          if (isnan(Msys)) {
-            Nsys++;
-            continue;
-          }
-          list[N] = Msys - Mcal - Mmos - Mgrid;
-          dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
-
-          // tie down reference photometry if the -refcode (code) option is selected
-          if (refPhotcode) {
-            if (GetPhotcodeEquivCodebyCode(catalog[i].measureT[m].photcode) == refPhotcode[0].equiv) {
-              // increase the weight by a factor of 100:
-              dlist[N] = 0.01*catalog[i].measureT[m].dM;
-            }
-          }
-          N++;
-        }
-
-        // when performing the grid analysis, STAR_TOOFEW will be set to 1;
-        if (N <= STAR_TOOFEW) { /* too few measurements */
-          catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_FEW;
-          Nfew ++;
+
+  StatType stats, apstats;
+  liststats_setmode (&stats, STATMODE);
+  liststats_setmode (&apstats, STATMODE);
+
+  double *list    = results->list;
+  double *dlist   = results->dlist;
+  double *wlist   = results->wlist;
+  double *aplist  = results->aplist;
+  double *daplist = results->daplist;
+
+  SetMrelInfoInit (results, FALSE); // do not allocate list,dlist,wlist arrays
+
+  int isSetMrelFinal = (pass >= 0);
+
+  for (j = 0; j < catalog[Nc].Naverage; j++) {
+    // XXX accumulate all secfilt values in a single pass?
+
+    // 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;
+    for (Ns = 0; Ns < Nphotcodes; Ns++) {
+
+      int thisCode = photcodes[Ns][0].code;
+      int Nsec = GetPhotcodeNsec(thisCode);
+
+      /* 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 (!isSetMrelFinal && (catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD)) continue;
+
+      int Ncode = 0;
+      int Next = 0;
+      int haveSynth = FALSE;
+
+      int forceSynth = FALSE;
+      int forceSynthEntry = -1;
+
+      int haveUbercal = FALSE;
+
+      int minUbercalDist = 1000;
+    
+      N = 0;
+      m = catalog[Nc].averageT[j].measureOffset;
+      for (k = 0; k < catalog[Nc].averageT[j].Nmeasure; k++, m++) {
+
+        // skip measurements that do not match the current photcode
+        PhotCode *code = GetPhotcodebyCode (catalog[Nc].measureT[m].photcode);
+        if (!code) continue;
+        if (code->equiv != thisCode) { continue; }
+        Ncode ++;
+
+        if (catalog[Nc].measureT[m].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); 
+
+        if (getImageEntry (m, Nc) < 0) {
+          // measurements without an image are either external reference photometry or
+          // data for which the associated image has not been loaded (probably because of
+          // overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
+          Mmos = Mgrid = 0;
+          Mcal = catalog[Nc].measureT[m].Mcal; // check that this is zero for loaded REF value
         } else {
-          catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_FEW;
-        }       
-
-        liststats (list, dlist, N, &stats);
-
-        catalog[i].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
-        catalog[i].secfilt[Nsecfilt*j+Nsec].dM = stats.sigma;
-        catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
-      }
-    }
-  }
-  fprintf (stderr, "%d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", Nfew, Nbad, Ncal, Nmos, Ngrid, Nsys);
-
-  return (TRUE);
-}
-
-int setMrelOutput (Catalog *catalog, int Ncatalog, int mark) {
-
-  int i, N;
-  off_t j, k, m, Nmax;
-  float Msys, Mcal, Mmos, Mgrid;
-  double *list, *dlist;
-  StatType stats;
-  int Nsec, Nsecfilt, ecode;
-
-  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));
-
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-      /* skip stars already calibrated */
-      if (catalog[i].found[j]) continue;  
-
-      int Ns;
-      for (Ns = 0; Ns < Nphotcodes; Ns++) {
-        int thisCode = photcodes[Ns][0].code;
-        Nsec = GetPhotcodeNsec(thisCode);
-
-        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
-          ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
-          if (ecode != thisCode) { continue; }
-
-          if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
-
-          // XXX allow REF stars (no Image Entry) to be included in the calculation this
-          // should be optionally set, and should allow for REF stars to be downweighted by
-          // more than their reported errors.  how such info is carried is unclear...
-          if (getImageEntry (m, i) < 0) {
-            Mcal = Mmos = Mgrid = 0;
-          } else {
-            Mcal  = getMcal  (m, i);
-            if (isnan(Mcal)) continue;
-            Mmos  = getMmos  (m, i);
-            if (isnan(Mmos)) continue;
-            Mgrid = getMgrid (m, i);
-            if (isnan(Mgrid)) continue;
-          }
-
-          Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt]);
-          list[N] = Msys - Mcal - Mmos - Mgrid;
-          dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
-          N++;
-        }
-        if (N < 1) continue;
-
-        liststats (list, dlist, N, &stats);
-        if (mark) catalog[i].found[j] = TRUE;
-
-        /* use sigma or error in dM for output? */
-        catalog[i].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
-        catalog[i].secfilt[Nsecfilt*j+Nsec].dM = MAX (stats.error, stats.sigma);
-        catalog[i].secfilt[Nsecfilt*j+Nsec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
-      }
-    }
-  }
-
-  free (list);
-  free (dlist);
-  return (TRUE);
-}
-
-// For each average object, set the average mags based on existing equiv photometry.
-// NOTE: this function operates on the real Measure & Average structures, not the
-// MeasureTiny & AverageTiny structures
-int setMave (Catalog *catalog, int Ncatalog) {
-
-  off_t j, k, m, Nmax;
-  int i, Ns, Nsecfilt, N, Nc;
-  float Msys;
-  double *list, *dlist;
-  StatType stats;
-  PhotCode *code;
-  DVOAverageFlags flagBits;
-
-  flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT;
-
-  // pre-allocate a list for stats purposes
-  Nmax = 0;
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-      Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
-    }
-  }
-  ALLOCATE (list, double, MAX (1, Nmax));
-  ALLOCATE (dlist, double, MAX (1, Nmax));
-
-  Nsecfilt = GetPhotcodeNsecfilt ();
-
-# define PSFQUALSTATS 1
-  for (i = 0; i < Ncatalog; i++) {
-    for (j = 0; j < catalog[i].Naverage; j++) {
-
-      // update average photometry for each of the average filters
-      for (Ns = 0; Ns < Nsecfilt; Ns++) {
-
-        code = GetPhotcodebyNsec (Ns);
-        Nc = code[0].code;
-        
-        N = 0;
-        m = catalog[i].average[j].measureOffset;
-        for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-          if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
-          if (GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode) != Nc) continue;
-
-          Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
-          if (isnan(Msys)) continue;
-
-          // XXX only apply this filter for psphot data from GPC1 for now...
-          if (PSFQUALSTATS && (catalog[i].measure[m].photcode > 10000) && (catalog[i].measure[m].photcode < 10500)) {
-              if (catalog[i].measure[m].psfQual < 0.85) continue;
-          }
-
-          // XXX make it optional to apply Mcal (do I need an extra field for advisory Mcal?)
-          list[N] = Msys - catalog[i].measure[m].Mcal;
-          dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
-          N++;
-        }
-        if (N < 1) continue;
-
-        liststats (list, dlist, N, &stats);
-
-        /* use sigma or error in dM for output? */
-        catalog[i].secfilt[Nsecfilt*j+Ns].M  = stats.mean;
-        catalog[i].secfilt[Nsecfilt*j+Ns].dM = MAX (stats.error, stats.sigma);
-        catalog[i].secfilt[Nsecfilt*j+Ns].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
-        catalog[i].secfilt[Nsecfilt*j+Ns].Ncode = N;
-        catalog[i].secfilt[Nsecfilt*j+Ns].Nused = stats.Nmeas;
-      }
-
-      // update average flags based on the detection stats.  
-      // XXX we need to clean this up -- this is a serious set of hacks...
-      if (PSFQUALSTATS) { 
-        int Galaxy2MASS, GalaxySDSS, goodPS1, nEXT, nPSF, good2MASS;
-
-        Galaxy2MASS = FALSE; // best guess for galaxy based on 2MASS J measurements (gal_contam == measure.flags[0x00400000 | 0x00800000])
-        GalaxySDSS = FALSE;  // best guess for galaxy based on SDSS measurements (XXX need to fix SDSS flags)
-        goodPS1 = FALSE;     // true if any PS1 measurements have psfQual > 0.85
-        good2MASS = FALSE;   // true if 2MASS J measurements have significant detections
-        nEXT = nPSF = 0;     // number of PS1 PSF vs EXT measurements
-
-        m = catalog[i].average[j].measureOffset;
-        for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
-
-          // PS1 data :
-          if ((catalog[i].measure[m].photcode >= 10000) && (catalog[i].measure[m].photcode <= 10500)) {
-            if (catalog[i].measure[m].psfQual > 0.85) {
-              goodPS1 = TRUE;
-              if (!isnan(catalog[i].measure[m].Map)) {
-                if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) {
-                  nEXT ++;
-                } else {
-                  nPSF ++;
-                }
+          Mcal  = getMcal  (m, Nc, flatcorr, catalog);
+          if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
+          Mmos  = getMmos  (m, Nc);
+          if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
+          Mgrid = getMgrid (m, Nc);
+          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)) 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;
+
+          // special options for PS1 data
+          if ((catalog[Nc].measure[m].photcode >= 10000) && (catalog[Nc].measure[m].photcode <= 10500)) {
+            // count the extended detections
+            if (!isnan(catalog[Nc].measure[m].Map)) {
+              if (catalog[Nc].measure[m].M - catalog[Nc].measure[m].Map > 0.5) {
+                Next ++;
+                NextPS1 ++;
+              } else {
+                NpsfPS1 ++;
               }
             }
           }
-          
-          // 2MASS data:
-          if (catalog[i].measure[m].photcode == 2011) {
-            if (catalog[i].measure[m].photFlags & 0x00c00000) {
+          // 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 (catalog[i].measure[m].photFlags & 0x00000007) {
-              good2MASS = 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;
+              }
             }
-          }  
-        }
-
-        // we attempt to set a few flags here; reset those bits before trying:
-        catalog[i].average[j].flags &= ~flagBits;
-
-        if (nEXT && (nEXT > nPSF)) {
-          catalog[i].average[j].flags |= ID_OBJ_EXT;
-        }
-        if (goodPS1) {
-          catalog[i].average[j].flags |= ID_OBJ_GOOD;
-        }
-        if (Galaxy2MASS) {
-          catalog[i].average[j].flags |= ID_OBJ_EXT_ALT;
-        }
-        if (good2MASS) {
-          catalog[i].average[j].flags |= ID_OBJ_GOOD_ALT;
-        }
-      }
-    }
-  }
-
-  free (list);
-  free (dlist);
+          }
+
+          // Blindly accepth the SYNTH mags if we are above saturation, otherwise, 
+          // ignore SYNTH photcodes until PASS == 4 (where we also accept saturated stars)
+          if ((catalog[Nc].measure[m].photcode >= 3001) && (catalog[Nc].measure[m].photcode <= 3005)) {
+            // something of a hack: force object to use synth values if synth mags >>
+            // saturation (3pi instrumental mags < -15)
+            float MaxMagForceSynth = NAN;
+            switch (catalog[Nc].measure[m].photcode) {
+              case 3001:
+                MaxMagForceSynth = 13.64;
+                break;
+              case 3002:
+                MaxMagForceSynth = 13.76;
+                break;
+              case 3003:
+                MaxMagForceSynth = 13.74;
+                break;
+              case 3004:
+                MaxMagForceSynth = 12.94;
+                break;
+              case 3005:
+                MaxMagForceSynth = 12.01;
+                break;
+            }
+            if (catalog[Nc].measureT[m].M < MaxMagForceSynth) {
+              forceSynth = TRUE;
+              forceSynthEntry = N;
+            } else {
+              if (pass < 4) {
+                SKIP_THIS_MEAS(Nbad);
+                continue;
+              }
+              haveSynth = TRUE;
+            }
+          }
+        }
+
+        // dlist gives the error per measurement, wlist gives the weight
+        // we can modify the error and 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[Nc].measureT[m].dM, code->photomErrSys), MIN_ERROR);
+        wlist[N] = 1.0;
+
+        // up-weight the ubercal values (or convergence can take a long time...)
+        if (catalog[Nc].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
+          wlist[N] = 10.0;
+        }
+
+        // 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) {
+            wlist[N] = 100.0;
+          }
+        }
+        N++;
+      }
+
+      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;
+        if (Ncode == 0) {
+          results->Ncode ++;
+        } else {
+          results->Nfew ++;
+        }
+        continue;
+      } else {
+        catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_FEW;
+      } 
+
+      if (forceSynth) {
+        // use the single SYNTH value instead of the other mags here
+        myAssert ((forceSynthEntry < N) && (forceSynthEntry >= 0), "programming error");
+        list[0]  = list[forceSynthEntry];
+        dlist[0] = dlist[forceSynthEntry];
+        wlist[0] = wlist[forceSynthEntry];
+        N = 1;
+      }
+      liststats (list, dlist, wlist, N, &stats);
+
+      catalog[Nc].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
+      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;
+
+      // 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) {
+        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 : use the modified weight for apmags as well as psf mags
+        liststats (aplist, daplist, wlist, N, &apstats);
+
+        catalog[Nc].secfilt[Nsecfilt*j+Nsec].Map  = apstats.mean; 
+
+        // 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 print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
+
+  off_t k;
+
+  int Nsecfilt = GetPhotcodeNsecfilt ();
+
+  off_t m = average[0].measureOffset;
+
+  for (k = 0; k < average[0].Nmeasure; k++, m++) {
+    fprintf (stderr, "meas: %08x\n", measure[m].dbFlags);
+  }
+
+  int Ns;
+  for (Ns = 0; Ns < Nsecfilt; Ns++) {
+    fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags);
+  }
+  return 1;
+}
+
 /* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_PHOTOM_NOCAL */
-int setMcalOutput (Catalog *catalog, int Ncatalog) {
+int setMcalOutput (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
 
   int i;
@@ -368 +633 @@
       for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
         if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
-        Mcal  = getMcal  (m, i);
+        Mcal  = getMcal  (m, i, flatcorr, catalog);
         if (isnan(Mcal)) continue;
         Mmos  = getMmos  (m, i);
@@ -375 +640 @@
         if (isnan(Mgrid)) continue;
 
+        // note that measurements for which the image is not selected will not be modified
+        // (that is a good thing! -- we keep a prior calibration)
+
         // set the output calibration
         catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid;
+
+        if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
+          myAssert (isfinite(catalog[i].measure[m].Mcal), "oops, broke an ubercal mag");
+        }
       }
     }
@@ -385 +657 @@
 void clean_stars (Catalog *catalog, int Ncatalog) {
 
-  int i, j, Ndel, Nave, Ntot, mark, Ns;
+  int i, j, Ndel, Nave, Ntot, mark, Ns, Nscat, Nchi, Nnan;
   float dM, Xm;
   double Chisq, MaxScatter, MaxChisq;
   double *xlist, *slist, *dlist;
+
   StatType stats;
+  liststats_setmode (&stats, "MEAN");
 
   if (VERBOSE) fprintf (stderr, "marking poor stars\n");
@@ -404 +678 @@
 
   // eliminate bad stars using the stats for a single secfilt at a time
-  // XXX DEP replace average.flags with secfilt flags
   for (Ns = 0; Ns < Nphotcodes; Ns ++) {
     
@@ -423 +696 @@
     }
   
-    initstats ("MEAN");
-    liststats (xlist, dlist, Ntot, &stats);
+    liststats (xlist, dlist, NULL, Ntot, &stats);
     MaxChisq = MAX (STAR_CHISQ, 2*stats.median);
-    liststats (slist, dlist, Ntot, &stats);
+
+    liststats (slist, dlist, NULL, Ntot, &stats);
     MaxScatter = MAX (STAR_SCATTER, 2*stats.median);
     fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq);
 
-    Ndel = Nave = 0;
+    Ndel = Nave = Nscat = Nnan = Nchi = 0;
     for (i = 0; i < Ncatalog; i++) {
       for (j = 0; j < catalog[i].Naverage; j++) {
@@ -440 +713 @@
           catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_POOR;
           Ndel ++;
+          if (dM > MaxScatter) { Nscat ++; }
+          if (Xm == NAN_S_SHORT) { Nnan ++; }
+          if (Chisq > MaxChisq) { Nchi ++; }
         } else {
           catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_POOR;
@@ -446 +722 @@
       }
     }
-    fprintf (stderr, "%d stars marked variable, %d total\n", Ndel, Nave);
-    initstats (STATMODE);
+    fprintf (stderr, "%d stars marked variable (%d scat, %d nan, %d chi), %d total\n", Ndel, Nscat, Nnan, Nchi, Nave);
   }
   free (xlist);
@@ -454 +729 @@
 }
 
+// clean_measures examines the stats for a single star.  It measures the INNER 50% mean
+// and sigma, it then re-measures the mean and sigma using all stars with NSIGMA_CLIP (3)
+// sigma of the INNER 50% mean.  it then flags any measurements which are more than
+// NSIGMA_REJECT (5) sigma of the mean
+
 # define NSIGMA_CLIP 3.0
 # define NSIGMA_REJECT 5.0
-void clean_measures (Catalog *catalog, int Ncatalog, int final) {
+void clean_measures (Catalog *catalog, int Ncatalog, int final, FlatCorrectionTable *flatcorr) {
 
   off_t j, k, m, Nmax, Ndel, Nave;
@@ -463 +743 @@
   double *tlist, *list, *dlist;
   float Msys, Mcal, Mmos, Mgrid;
-  StatType stats;
   int Ncal, Nmos, Ngrid, Nfew;
 
@@ -484 +763 @@
   TOOFEW = MAX (5, STAR_TOOFEW);
 
+  // stats structures for inner and full stats
+  StatType instats, stats;
+  liststats_setmode (&instats, "INNER_MEAN");
+  liststats_setmode (&stats, "MEAN");
+
   Ndel = Nave = 0;
   Ncal = Nmos = Ngrid = Nfew = 0;
@@ -492 +776 @@
       for (Ns = 0; Ns < Nphotcodes; Ns++) {
         
-        /* on final processing, skip stars already measured */
-        if (final && catalog[i].found[j]) continue;  
-
         int thisCode = photcodes[Ns][0].code;
         int Nsec = GetPhotcodeNsec(thisCode);
         
+        /* on final processing, skip stars already measured */
+        if (final && catalog[i].found[Nsecfilt*j + Nsec]) continue;  
+
         /* skip bad stars to prevent them from becoming good (on inner sample) */
         if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
@@ -509 +793 @@
           if (ecode != thisCode) { continue; }
 
-          /* if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; */
-          Mcal  = getMcal  (m, i);
+          // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
+          Mcal  = getMcal  (m, i, flatcorr, catalog);
           if (isnan(Mcal)) { Ncal ++; continue; }
           Mmos  = getMmos  (m, i);
@@ -522 +806 @@
           N++;
         }
-        if (N <= TOOFEW) { Nfew ++; continue; }
+        if (N <= TOOFEW) { 
+          Nfew ++; 
+          continue; 
+        }
 
         /* 3-sigma clip based on stats of inner 50% */
 
         // calculated mean of inner 50%
-        initstats ("INNER_MEAN");
-        liststats (list, dlist, N, &stats);
-        stats.sigma = MAX (MIN_ERROR, stats.sigma); /* if measurements agree too well, sigma -> 0.0 */
+        liststats (list, dlist, NULL, N, &instats);
+        instats.sigma = MAX (MIN_ERROR, instats.sigma); /* if measurements agree too well, sigma -> 0.0 */
 
         // ignore entries > 3sigma from inner mean
         for (k = m = 0; k < N; k++) {
-          if (fabs (list[k] - stats.median) < NSIGMA_CLIP*stats.sigma) {
+          if (fabs (list[k] - instats.median) < NSIGMA_CLIP*instats.sigma) {
             list[m] = list[k];
             m++;
@@ -539 +825 @@
         }
         // recalculate the mean & sigma of the accepted measurements
-        initstats ("MEAN");
-        liststats (list, dlist, m, &stats);
+        liststats (list, dlist, NULL, m, &stats);
         stats.sigma = MAX (MIN_ERROR, stats.sigma);
 
@@ -553 +838 @@
           if (ecode != thisCode) { continue; }
 
-          /* if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue; */
-          Mcal  = getMcal  (m, i);
+          // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
+          Mcal  = getMcal  (m, i, flatcorr, catalog);
           if (isnan(Mcal)) continue;
           Mmos  = getMmos  (m, i);
@@ -572 +857 @@
         /* mark bad measures (> 3 sigma deviant) */
         for (k = 0; k < N; k++) {
-          if (fabs (list[k] - stats.median) > NSIGMA_REJECT*stats.sigma) {
+          // treat the scatter of the star as a systematic term.  this is a bit of an
+          // over-estimage (a perfect Gauss distribution with perfect errors would have
+          // mySigma = sqrt(2) too large)
+          float mySigma = hypot (stats.sigma, dlist[k]);
+          if (fabs (list[k] - stats.median) > NSIGMA_REJECT*mySigma) {
             catalog[i].measureT[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
             if (final) {
@@ -585 +874 @@
     }
   }
-  initstats (STATMODE);
   if (VERBOSE) fprintf (stderr, OFF_T_FMT" measures marked poor, "OFF_T_FMT" total\n", Ndel, Nave);
+
   free (list);
   free (dlist);
@@ -593 +882 @@
 }
 
-StatType statsStarN (Catalog *catalog, int Ncatalog, int Nsec, int seccode) {
+StatType statsStarN (Catalog *catalog, int Ncatalog, int Nsec, int seccode, FlatCorrectionTable *flatcorr) {
 
   off_t j, k, m, Ntot;
@@ -625 +914 @@
         int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
         if (ecode != seccode) { continue;}
-        Mcal = getMcal  (m, i);
+        Mcal = getMcal  (m, i, flatcorr, catalog);
         if (isnan(Mcal)) { continue;}
         Mmos = getMmos  (m, i);
@@ -640 +929 @@
   }
 
-  // fprintf (stderr, "N1: %d, N2: %d, N3: %d, N4: %d, N0: %d\n", N1, N2, N3, N4, N0);
-  liststats (list, dlist, n, &stats);
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (list, dlist, NULL, n, &stats);
   free (list);
   free (dlist);
@@ -680 +970 @@
   }
 
-  liststats (list, dlist, n, &stats);
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (list, dlist, NULL, n, &stats);
   free (list);
   free (dlist);
@@ -712 +1004 @@
 
       dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
+      if (isnan(dM)) continue;
       list[n] = dM;
       dlist[n] = 1;
@@ -718 +1011 @@
   }
 
-  liststats (list, dlist, n, &stats);
+  liststats_setmode (&stats, "MEAN");
+
+  liststats (list, dlist, NULL, n, &stats);
   free (list);
   free (dlist);