IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 17, 2012, 9:51:31 AM (14 years ago)
Author:
eugene
Message:

convert liststats to use a statmode element carried by the structure (instead of the global variable); add a weight term to liststats to allow more complex weighting schemes than constant or inverse-error

Location:
branches/eam_branches/ipp-20120405/Ohana/src/relphot
Files:
8 edited

Legend:

Unmodified
Added
Removed

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

    r33885 r33886  
    3636} Mosaic;
    3737
     38typedef enum {
     39  STATS_NONE,
     40  STATS_MEAN,
     41  STATS_MEDIAN,
     42  STATS_WT_MEAN,
     43  STATS_INNER_MEAN,
     44  STATS_INNER_WTMEAN,
     45  STATS_CHI_INNER_MEAN,
     46  STATS_CHI_INNER_WTMEAN
     47} ListStatsMode;
     48
    3849typedef struct {
    3950  double median;
    … …  
    4859  double total;
    4960  int    Nmeas;
     61  ListStatsMode statmode;
    5062} StatType;
    5163
    … …  
    245257void          initialize          PROTO((int argc, char **argv));
    246258void          initialize_client   PROTO((int argc, char **argv));
    247 void          initstats           PROTO((char *mode));
    248 int           liststats           PROTO((double *value, double *dvalue, int N, StatType *stats));
     259void          liststats_setmode   PROTO((StatType *stats, char *strmode));
     260int           liststats           PROTO((double *value, double *dvalue, double *wvalue, int N, StatType *stats));
    249261Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath));
    250262Catalog      *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog));

  • branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/GridOps.c

    r33651 r33886  
    527527  double *list, *dlist;
    528528  float Msys, Mrel, Mcal, Mmos;
     529
    529530  StatType stats;
     531  liststats_setmode (&stats, STATMODE);
    530532 
    531533  if (!USE_GRID) return;
    … …  
    595597    }
    596598
    597     liststats (list, dlist, N, &stats);
     599    liststats (list, dlist, NULL, N, &stats);
    598600    gridM[i] = stats.mean;
    599601    gridS[i] = stats.sigma;

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

    r33823 r33886  
    399399  float Msys, Mrel, Mmos, Mgrid, Mflat;
    400400  double *list, *dlist, *Mlist, *dMlist;
     401
    401402  StatType stats;
     403  liststats_setmode (&stats, STATMODE);
    402404
    403405  if (FREEZE_IMAGES) return;
    … …  
    518520    }
    519521
    520     liststats (list, dlist, N, &stats);
     522    // no additional weight modification (we treat all stars on an image equally -- note an image is either ubercal-tied or not)
     523    liststats (list, dlist, NULL, N, &stats);
    521524    image[i].Mcal  = stats.mean;
    522525    image[i].dMcal = stats.error;
    … …  
    527530
    528531    // bright end scatter
    529     liststats (Mlist, dMlist, Nbright, &stats);
     532    liststats (Mlist, dMlist, NULL, Nbright, &stats);
    530533    image[i].dMagSys = stats.sigma;
    531534
    … …  
    560563  double *mlist, *slist, *dlist;
    561564  double MaxOffset, MaxScatter, MedOffset;
    562   StatType stats;
    563565
    564566  if (FREEZE_IMAGES) return;
    … …  
    578580    N++;
    579581  }
    580   initstats ("MEAN");
    581   liststats (mlist, dlist, N, &stats);
     582
     583  // use a straight mean to find the global image statistics (no weighting)
     584  StatType stats;
     585  liststats_setmode (&stats, "MEAN");
     586
     587  liststats (mlist, dlist, NULL, N, &stats);
    582588  MaxOffset = MAX (IMAGE_OFFSET, 3*stats.sigma);
    583589  MedOffset = stats.median;
    584   liststats (slist, dlist, N, &stats);
     590
     591  liststats (slist, dlist, NULL, N, &stats);
    585592  MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
    586593  fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
    … …  
    603610
    604611  fprintf (stderr, "%d images marked poor\n", Nmark);
    605   initstats (STATMODE);
    606612  free (mlist);
    607613  free (slist);
    … …  
    724730  double *list, *dlist;
    725731  float Mcal, Mmos, Mgrid;
     732
    726733  StatType stats;
    727 
    728734  bzero (&stats, sizeof (StatType));
    729735  if (FREEZE_IMAGES) return (stats);
    … …  
    755761  }
    756762
    757   liststats (list, dlist, n, &stats);
     763  liststats_setmode (&stats, "MEAN");
     764
     765  liststats (list, dlist, NULL, n, &stats);
    758766  free (list);
    759767  free (dlist);
    … …  
    783791  }
    784792
    785   liststats (list, dlist, n, &stats);
     793  liststats_setmode (&stats, "MEAN");
     794
     795  liststats (list, dlist, NULL, n, &stats);
    786796  free (list);
    787797  free (dlist);
    … …  
    811821  }
    812822
    813   liststats (list, dlist, n, &stats);
     823  liststats_setmode (&stats, "MEAN");
     824
     825  liststats (list, dlist, NULL, n, &stats);
    814826  free (list);
    815827  free (dlist);
    … …  
    839851  }
    840852
    841   liststats (list, dlist, n, &stats);
     853  liststats_setmode (&stats, "MEAN");
     854
     855  liststats (list, dlist, NULL, n, &stats);
    842856  free (list);
    843857  free (dlist);

  • branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/MosaicOps.c

    r33828 r33886  
    764764  Image *imageReal;
    765765  off_t j, NimageReal;
     766
    766767  StatType stats;
     768  liststats_setmode (&stats, STATMODE);
    767769
    768770  double *list   = info->list;
    … …  
    930932    }
    931933  }
    932   liststats (list, dlist, N, &stats);
     934
     935  liststats (list, dlist, NULL, N, &stats);
    933936  if (VERBOSE2 && info->PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
    934937
    … …  
    945948
    946949  // bright end scatter
    947   liststats (Mlist, dMlist, Nbright, &stats);
     950  liststats (Mlist, dMlist, NULL, Nbright, &stats);
    948951  myMosaic[0].dMsys = stats.sigma;
    949952
    … …  
    13101313  }
    13111314
    1312   liststats (list, dlist, n, &stats);
     1315  liststats_setmode (&stats, "MEAN");
     1316
     1317  liststats (list, dlist, NULL, n, &stats);
    13131318  free (list);
    13141319  free (dlist);
    … …  
    13381343  }
    13391344
    1340   liststats (list, dlist, n, &stats);
     1345  liststats_setmode (&stats, "MEAN");
     1346
     1347  liststats (list, dlist, NULL, n, &stats);
    13411348  free (list);
    13421349  free (dlist);
    … …  
    13831390  // fprintf (stderr, "Nmosaic: "OFF_T_FMT", n: "OFF_T_FMT"\n",  Nmosaic,  n);
    13841391
    1385   liststats (list, dlist, n, &stats);
     1392  liststats_setmode (&stats, "MEAN");
     1393
     1394  liststats (list, dlist, NULL, n, &stats);
    13861395  free (list);
    13871396  free (dlist);
    … …  
    14111420  }
    14121421
    1413   liststats (list, dlist, n, &stats);
     1422  liststats_setmode (&stats, "MEAN");
     1423
     1424  liststats (list, dlist, NULL, n, &stats);
    14141425  free (list);
    14151426  free (dlist);
    … …  
    14231434  double *mlist, *slist, *dlist;
    14241435  double MaxOffset, MedOffset, MaxScatter;
    1425   StatType stats;
    14261436
    14271437  if (!MOSAIC_ZEROPT) return;
    … …  
    14421452    N++;
    14431453  }
    1444   initstats ("MEAN");
    1445   liststats (mlist, dlist, N, &stats);
     1454
     1455  StatType stats;
     1456  liststats_setmode (&stats, "MEAN");
     1457
     1458  liststats (mlist, dlist, NULL, N, &stats);
    14461459  MaxOffset = MAX (IMAGE_OFFSET, 2*stats.sigma);
    14471460  MedOffset = stats.median;
    1448   liststats (slist, dlist, N, &stats);
     1461
     1462  liststats (slist, dlist, NULL, N, &stats);
    14491463  MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
    14501464  fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
    … …  
    14761490
    14771491  fprintf (stderr, OFF_T_FMT" mosaics marked poor ("OFF_T_FMT" scatter, "OFF_T_FMT" offset)\n",  Nmark, Nscatter, Noffset);
    1478   initstats (STATMODE);
     1492
    14791493  free (mlist);
    14801494  free (slist);

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

    r33827 r33886  
    1313  double *list;
    1414  double *dlist;
     15  double *wlist;
    1516  double *aplist;
    1617  double *daplist;
    … …  
    7980    ALLOCATE (results->list,  double, Nmax);
    8081    ALLOCATE (results->dlist, double, Nmax);
     82    ALLOCATE (results->wlist, double, Nmax);
    8183  }
    8284}
    … …  
    8587  free (results->list);
    8688  free (results->dlist);
     89  free (results->wlist);
    8790}
    8891
    … …  
    159162  SetMrelInfo summary, results;
    160163  SetMrelInfoInit (&summary, FALSE);
    161   SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist
     164  SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist,wlist
    162165  ALLOCATE (results.aplist, double, Nmax);
    163166  ALLOCATE (results.daplist, double, Nmax);
    … …  
    303306  int N;
    304307  float Msys, Mcal, Mmos, Mgrid;
     308
    305309  StatType stats, apstats;
     310  liststats_setmode (&stats, STATMODE);
     311  liststats_setmode (&apstats, STATMODE);
    306312
    307313  double *list    = results->list;
    308314  double *dlist   = results->dlist;
     315  double *wlist   = results->wlist;
    309316  double *aplist  = results->aplist;
    310317  double *daplist = results->daplist;
    311318
    312   SetMrelInfoInit (results, FALSE); // do not allocate list,dlist arrays
     319  SetMrelInfoInit (results, FALSE); // do not allocate list,dlist,wlist arrays
    313320
    314321  int isSetMrelFinal = (pass >= 0);
    … …  
    390397          aplist[N] = Map - Mcal - Mmos - Mgrid;
    391398
    392           // count the extended detections
     399          // special options for PS1 data
    393400          if ((catalog[Nc].measure[m].photcode >= 10000) && (catalog[Nc].measure[m].photcode <= 10500)) {
     401            // count the extended detections
    394402            if (!isnan(catalog[Nc].measure[m].Map)) {
    395403              if (catalog[Nc].measure[m].M - catalog[Nc].measure[m].Map > 0.5) {
    … …  
    427435        }
    428436
    429         // dlist gives the error, which is used as the weight in WT_MEAN.
    430         // we can modify the resulting weight in a few ways:
     437        // dlist gives the error per measurement, wlist gives the weight
     438        // we can modify the error and weight in a few ways:
    431439        // 1) MIN_ERROR guarantees a floor
    432440        // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
    … …  
    434442        // 4) some reference photcode of some kind can be specified as fixed and have a high weight
    435443        dlist[N] = MAX (hypot(catalog[Nc].measureT[m].dM, code->photomErrSys), MIN_ERROR);
     444        wlist[N] = 1.0;
    436445
    437446        // up-weight the ubercal values (or convergence can take a long time...)
    438447        if (catalog[Nc].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
    439           dlist[N] = MAX (0.1*catalog[Nc].measureT[m].dM, MIN_ERROR);
     448          wlist[N] = 10.0;
    440449        }
    441450
    … …  
    445454        if (refPhotcode) {
    446455          if (code->code == refPhotcode->code) {
    447             // tiny error -> large weight
    448             // dlist[N] = MAX (0.01*catalog[Nc].measureT[m].dM, MIN_ERROR);
    449             dlist[N] = 0.0001;
     456            wlist[N] = 100.0;
    450457          }
    451458        }
    … …  
    469476      }
    470477
    471       liststats (list, dlist, N, &stats);
     478      liststats (list, dlist, wlist, N, &stats);
    472479
    473480      catalog[Nc].secfilt[Nsecfilt*j+Nsec].M  = stats.mean;
    … …  
    481488
    482489      if (isSetMrelFinal) {
    483         liststats (aplist, daplist, N, &apstats);
    484490        catalog[Nc].found[Nsecfilt*j+Nsec] = TRUE;
    485491
    … …  
    490496        catalog[Nc].secfilt[Nsecfilt*j+Nsec].M_80 = 1000 * stats.Upper80;
    491497        catalog[Nc].secfilt[Nsecfilt*j+Nsec].M_20 = 1000 * stats.Lower20;
     498
     499        // NOTE : use the modified weight for apmags as well as psf mags
     500        liststats (aplist, daplist, wlist, N, &apstats);
     501
     502        catalog[Nc].secfilt[Nsecfilt*j+Nsec].Map  = apstats.mean;
    492503
    493504        // NOTE: for 2MASS measurements, Next should be 1, as should N
    … …  
    567578}
    568579
    569 # if (0)
    570 
    571 /* grab Nsec for named photcode */
    572 # define NAMED_PHOTCODE_NSEC(MY_NSEC,NAME)      \
    573   short MY_NSEC = -1;                           \
    574   {                                             \
    575   PhotCode *code = GetPhotcodebyName (NAME);    \
    576   if (code) {                                   \
    577     MY_NSEC = GetPhotcodeNsec (code->equiv);    \
    578   } }
    579 
    580 // For each average object, set the average mags based on existing equiv photometry.
    581 // NOTE: this function operates on the real Measure & Average structures, not the
    582 // MeasureTiny & AverageTiny structures
    583 // NOTE: this function is called on the remote machine -- make sure any (global) options
    584 // are passed to the relphot_client program
    585 int setMave (Catalog *catalog, int Ncatalog) {
    586 
    587   off_t j, k, m, Nmax;
    588   int i, Ns, Nsecfilt, N, Nc;
    589   float Msys;
    590   double *list, *dlist;
    591   StatType stats;
    592   DVOAverageFlags flagBits;
    593 
    594   flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT;
    595 
    596   // pre-allocate a list for stats purposes
    597   Nmax = 0;
    598   for (i = 0; i < Ncatalog; i++) {
    599     for (j = 0; j < catalog[i].Naverage; j++) {
    600       Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
    601     }
    602   }
    603   ALLOCATE (list, double, MAX (1, Nmax));
    604   ALLOCATE (dlist, double, MAX (1, Nmax));
    605 
    606   Nsecfilt = GetPhotcodeNsecfilt ();
    607 
    608   // we want to raise some bits on the 2MASS (JHK) secfilt flags, if we have 2MASS data
    609   NAMED_PHOTCODE_NSEC (Nsec_J, "2MASS_J");
    610   NAMED_PHOTCODE_NSEC (Nsec_H, "2MASS_H");
    611   NAMED_PHOTCODE_NSEC (Nsec_K, "2MASS_K");
    612 
    613   for (i = 0; i < Ncatalog; i++) {
    614     for (j = 0; j < catalog[i].Naverage; j++) {
    615 
    616       // update average photometry for each of the average filters
    617 
    618       // XXX Note that this would be faster if we had an array of results and accumulated
    619       // them in a single pass
    620 
    621       int Next = 0;
    622 
    623       for (Ns = 0; Ns < Nsecfilt; Ns++) {
    624 
    625         PhotCode *code = GetPhotcodebyNsec (Ns);
    626         Nc = code[0].code;
    627        
    628         N = 0;
    629         m = catalog[i].average[j].measureOffset;
    630         for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    631           if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
    632           if (GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode) != Nc) continue;
    633 
    634           Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
    635           if (isnan(Msys)) continue;
    636 
    637           // reject POOR detections (PSF_QF < 0.85) and count extended detections
    638           // XXX only apply this filter for psphot data from GPC1 for now...
    639           if ((catalog[i].measure[m].photcode > 10000) && (catalog[i].measure[m].photcode < 10500)) {
    640             if (catalog[i].measure[m].psfQual < 0.85) continue;
    641             if (!isnan(catalog[i].measure[m].Map)) {
    642               if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) {
    643                 Next ++;
    644               }
    645             }
    646           }
    647 
    648           // XXX make it optional to apply Mcal (do I need an extra field for advisory Mcal?)
    649           list[N] = Msys - catalog[i].measure[m].Mcal;
    650           dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
    651           N++;
    652         }
    653         if (N < 1) continue;
    654 
    655         liststats (list, dlist, N, &stats);
    656 
    657         /* use sigma or error in dM for output? */
    658         catalog[i].secfilt[Nsecfilt*j+Ns].M  = stats.mean;
    659         catalog[i].secfilt[Nsecfilt*j+Ns].dM = stats.error;
    660         catalog[i].secfilt[Nsecfilt*j+Ns].Mstdev = 1000.0*stats.sigma; // Mstdev is in millimags (not enough space for more precision)
    661         catalog[i].secfilt[Nsecfilt*j+Ns].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NAN_S_SHORT;
    662         catalog[i].secfilt[Nsecfilt*j+Ns].Ncode = N;
    663         catalog[i].secfilt[Nsecfilt*j+Ns].Nused = stats.Nmeas;
    664 
    665         catalog[i].secfilt[Nsecfilt*j+Ns].M_80 = 1000 * stats.Upper80;
    666         catalog[i].secfilt[Nsecfilt*j+Ns].M_20 = 1000 * stats.Lower20;
    667 
    668         if ((Next > 0) && (Next > 0.5*N)) {
    669           catalog[i].secfilt[Nsecfilt*j+Ns].flags |= ID_SECF_OBJ_EXT;
    670         }
    671       }
    672 
    673       // update average flags based on the detection stats. 
    674 
    675       int Galaxy2MASS = FALSE; // best guess for galaxy based on 2MASS J measurements (gal_contam == measure.flags[0x00400000 | 0x00800000])
    676       int goodPS1 = FALSE;     // true if any PS1 measurements have psfQual > 0.85
    677       int good2MASS = FALSE;   // true if 2MASS J measurements have significant detections
    678       int nEXT = 0;
    679       int nPSF = 0;     // number of PS1 PSF vs EXT measurements
    680       int have2MASS = FALSE;
    681      
    682       // count, flag good and extended detections
    683       m = catalog[i].average[j].measureOffset;
    684       for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
    685 
    686         // PS1 data :
    687         if ((catalog[i].measure[m].photcode >= 10000) && (catalog[i].measure[m].photcode <= 10500)) {
    688           if (catalog[i].measure[m].psfQual > 0.85) {
    689             goodPS1 = TRUE;
    690             if (!isnan(catalog[i].measure[m].Map)) {
    691               if (catalog[i].measure[m].M - catalog[i].measure[m].Map > 0.5) {
    692                 nEXT ++;
    693               } else {
    694                 nPSF ++;
    695               }
    696             }
    697           }
    698         }
    699          
    700         // 2MASS data J-band flags
    701         if (catalog[i].measure[m].photcode == 2011) {
    702           // only need to do this once (always have JHK triplet; galaxy flag is same for all 3)
    703           have2MASS = TRUE;
    704           if (catalog[i].measure[m].photFlags & 0x00c00000) {
    705             Galaxy2MASS = TRUE;
    706           }
    707           if (catalog[i].measure[m].photFlags & 0x00000007) {
    708             good2MASS = TRUE;
    709             if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_PHOTOM_PASS_0;
    710           } else {
    711             if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_PHOTOM_PASS_1;
    712           }
    713         } 
    714         // 2MASS data H-band flags
    715         if (catalog[i].measure[m].photcode == 2012) {
    716           if (catalog[i].measure[m].photFlags & 0x00000007) {
    717             good2MASS = TRUE;
    718             if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_PHOTOM_PASS_0;
    719           } else {
    720             if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_PHOTOM_PASS_1;
    721           }
    722         } 
    723         // 2MASS data K-band flags
    724         if (catalog[i].measure[m].photcode == 2013) {
    725           if (catalog[i].measure[m].photFlags & 0x00000007) {
    726             good2MASS = TRUE;
    727             if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_PHOTOM_PASS_0;
    728           } else {
    729             if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_PHOTOM_PASS_1;
    730           }
    731         } 
    732       }
    733 
    734       // we attempt to set a few flags here; reset those bits before trying:
    735       catalog[i].average[j].flags &= ~flagBits;
    736 
    737       // XXX set the secfilt bits?
    738       if (nEXT && (nEXT > nPSF)) {
    739         catalog[i].average[j].flags |= ID_OBJ_EXT;
    740       }
    741       if (goodPS1) {
    742         catalog[i].average[j].flags |= ID_OBJ_GOOD;
    743       }
    744       if (Galaxy2MASS) {
    745         catalog[i].average[j].flags |= ID_OBJ_EXT_ALT;
    746         if (Nsec_J > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_J].flags |= ID_SECF_OBJ_EXT;
    747         if (Nsec_H > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_H].flags |= ID_SECF_OBJ_EXT;
    748         if (Nsec_K > -1) catalog[i].secfilt[j*Nsecfilt + Nsec_K].flags |= ID_SECF_OBJ_EXT;
    749       }
    750       if (good2MASS) {
    751         catalog[i].average[j].flags |= ID_OBJ_GOOD_ALT;
    752       }
    753     }
    754   }
    755 
    756   free (list);
    757   free (dlist);
    758   return (TRUE);
    759 }
    760 # endif
    761 
    762580/* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_PHOTOM_NOCAL */
    763581int setMcalOutput (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {
    … …  
    804622  double Chisq, MaxScatter, MaxChisq;
    805623  double *xlist, *slist, *dlist;
     624
    806625  StatType stats;
     626  liststats_setmode (&stats, "MEAN");
    807627
    808628  if (VERBOSE) fprintf (stderr, "marking poor stars\n");
    … …  
    818638  int Nsecfilt = GetPhotcodeNsecfilt ();
    819639
    820   // XX int oldPLOTSTUFF = PLOTSTUFF;
    821   // XX PLOTSTUFF = TRUE;
    822   // XX plot_chisq (catalog, Ncatalog);
    823   // XX PLOTSTUFF = oldPLOTSTUFF;
    824 
    825640  // eliminate bad stars using the stats for a single secfilt at a time
    826   // XXX DEP replace average.flags with secfilt flags
    827641  for (Ns = 0; Ns < Nphotcodes; Ns ++) {
    828642   
    … …  
    843657    }
    844658 
    845     initstats ("MEAN");
    846     liststats (xlist, dlist, Ntot, &stats);
     659    liststats (xlist, dlist, NULL, Ntot, &stats);
    847660    MaxChisq = MAX (STAR_CHISQ, 2*stats.median);
    848     liststats (slist, dlist, Ntot, &stats);
     661
     662    liststats (slist, dlist, NULL, Ntot, &stats);
    849663    MaxScatter = MAX (STAR_SCATTER, 2*stats.median);
    850664    fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq);
    … …  
    870684    }
    871685    fprintf (stderr, "%d stars marked variable (%d scat, %d nan, %d chi), %d total\n", Ndel, Nscat, Nnan, Nchi, Nave);
    872     initstats (STATMODE);
    873686  }
    874687  free (xlist);
    … …  
    876689  free (dlist);
    877690}
     691
     692// clean_measures examines the stats for a single star.  It measures the INNER 50% mean
     693// and sigma, it then re-measures the mean and sigma using all stars with NSIGMA_CLIP (3)
     694// sigma of the INNER 50% mean.  it then flags any measurements which are more than
     695// NSIGMA_REJECT (5) sigma of the mean
    878696
    879697# define NSIGMA_CLIP 3.0
    … …  
    886704  double *tlist, *list, *dlist;
    887705  float Msys, Mcal, Mmos, Mgrid;
    888   StatType stats;
    889706  int Ncal, Nmos, Ngrid, Nfew;
    890707
    … …  
    906723  /* it makes no sense to mark 3-sigma outliers with <5 measurements */
    907724  TOOFEW = MAX (5, STAR_TOOFEW);
     725
     726  // stats structures for inner and full stats
     727  StatType instats, stats;
     728  liststats_setmode (&instats, "INNER_MEAN");
     729  liststats_setmode (&stats, "MEAN");
    908730
    909731  Ndel = Nave = 0;
    … …  
    953775
    954776        // calculated mean of inner 50%
    955         initstats ("INNER_MEAN");
    956         liststats (list, dlist, N, &stats);
    957         stats.sigma = MAX (MIN_ERROR, stats.sigma); /* if measurements agree too well, sigma -> 0.0 */
     777        liststats (list, dlist, NULL, N, &instats);
     778        instats.sigma = MAX (MIN_ERROR, instats.sigma); /* if measurements agree too well, sigma -> 0.0 */
    958779
    959780        // ignore entries > 3sigma from inner mean
    960781        for (k = m = 0; k < N; k++) {
    961           if (fabs (list[k] - stats.median) < NSIGMA_CLIP*stats.sigma) {
     782          if (fabs (list[k] - instats.median) < NSIGMA_CLIP*instats.sigma) {
    962783            list[m] = list[k];
    963784            m++;
    … …  
    965786        }
    966787        // recalculate the mean & sigma of the accepted measurements
    967         initstats ("MEAN");
    968         liststats (list, dlist, m, &stats);
     788        liststats (list, dlist, NULL, m, &stats);
    969789        stats.sigma = MAX (MIN_ERROR, stats.sigma);
    970790
    … …  
    1015835    }
    1016836  }
    1017   initstats (STATMODE);
    1018837  if (VERBOSE) fprintf (stderr, OFF_T_FMT" measures marked poor, "OFF_T_FMT" total\n", Ndel, Nave);
     838
    1019839  free (list);
    1020840  free (dlist);
    1021841  free (ilist);
    1022842  free (tlist);
    1023 
    1024 
    1025843}
    1026844
    … …  
    1072890  }
    1073891
    1074   // fprintf (stderr, "N1: %d, N2: %d, N3: %d, N4: %d, N0: %d\n", N1, N2, N3, N4, N0);
    1075   liststats (list, dlist, n, &stats);
     892  liststats_setmode (&stats, "MEAN");
     893
     894  liststats (list, dlist, NULL, n, &stats);
    1076895  free (list);
    1077896  free (dlist);
    … …  
    1112931  }
    1113932
    1114   liststats (list, dlist, n, &stats);
     933  liststats_setmode (&stats, "MEAN");
     934
     935  liststats (list, dlist, NULL, n, &stats);
    1115936  free (list);
    1116937  free (dlist);
    … …  
    1151972  }
    1152973
    1153   liststats (list, dlist, n, &stats);
     974  liststats_setmode (&stats, "MEAN");
     975
     976  liststats (list, dlist, NULL, n, &stats);
    1154977  free (list);
    1155978  free (dlist);

  • branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/global_stats.c

    r33820 r33886  
    1212
    1313  // struct timeval startTimer, stopTimer;
    14 
    15   initstats ("MEAN");
    1614
    1715  fprintf (stderr, "\n");
    … …  
    6058  fprintf (stderr, "dMcal mosaic:   %7.4f  %7.4f  %7.4f  %7.4f  %7.4f  %6d\n",   msD.median, msD.mean, msD.sigma, msD.min, msD.max, msD.Nmeas);
    6159  fprintf (stderr, "chisq mosaic:   %7.1f  %7.1f  %7.1f  %7.1f  %7.1f  %6d\n",   msX.median, msX.mean, msX.sigma, msX.min, msX.max, msX.Nmeas);
    62  
    63 
    64   initstats (STATMODE);
    65 
    6660}
    6761

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

    r33823 r33886  
    3333    exit (1);
    3434  }
    35 
    36   initstats (STATMODE);
    3735
    3836  IMAGE_BAD = ID_IMAGE_PHOTOM_POOR | ID_IMAGE_PHOTOM_FEW | ID_IMAGE_PHOTOM_SKIP;

  • branches/eam_branches/ipp-20120405/Ohana/src/relphot/src/liststats.c

    r33651 r33886  
    11# include "relphot.h"
    22
    3 enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN,
    4       M_INNER_WTMEAN, M_CHI_INNER_MEAN, M_CHI_INNER_WTMEAN};
     3void liststats_setmode (StatType *stats, char *strmode) {
    54
    6 static int statmode;
     5  stats->statmode = -1;
     6  if (!strcmp (strmode, "MEAN"))             { stats->statmode = STATS_MEAN; return; }
     7  if (!strcmp (strmode, "MEDIAN"))           { stats->statmode = STATS_MEDIAN; return; }
     8  if (!strcmp (strmode, "WT_MEAN"))          { stats->statmode = STATS_WT_MEAN; return; }
     9  if (!strcmp (strmode, "INNER_MEAN"))       { stats->statmode = STATS_INNER_MEAN; return; }
     10  if (!strcmp (strmode, "INNER_WTMEAN"))     { stats->statmode = STATS_INNER_WTMEAN; return; }
     11  if (!strcmp (strmode, "CHI_INNER_MEAN"))   { stats->statmode = STATS_CHI_INNER_MEAN; return; }
     12  if (!strcmp (strmode, "CHI_INNER_WTMEAN")) { stats->statmode = STATS_CHI_INNER_WTMEAN; return; }
    713
    8 void initstats (char *mode) {
    9 
    10   statmode = -1;
    11   if (!strcmp (mode, "MEAN")) statmode = M_MEAN;
    12   if (!strcmp (mode, "MEDIAN")) statmode = M_MEDIAN;
    13   if (!strcmp (mode, "WT_MEAN")) statmode = M_WT_MEAN;
    14   if (!strcmp (mode, "INNER_MEAN")) statmode = M_INNER_MEAN;
    15   if (!strcmp (mode, "INNER_WTMEAN")) statmode = M_INNER_WTMEAN;
    16   if (!strcmp (mode, "CHI_INNER_MEAN")) statmode = M_CHI_INNER_MEAN;
    17   if (!strcmp (mode, "CHI_INNER_WTMEAN")) statmode = M_CHI_INNER_WTMEAN;
    18 
    19   if (statmode == -1) {
    20     fprintf (stderr, "ERROR: invalid stats mode: %s\n", mode);
    21     exit (1);
    22   }
     14  fprintf (stderr, "ERROR: invalid stats mode: %s\n", strmode);
     15  exit (1);
    2316}
    2417
    25 int liststats (double *value, double *dvalue, int N, StatType *stats) {
     18int liststats (double *value, double *dvalue, double *weight, int N, StatType *stats) {
    2619 
    27   int i, ks, ke, Nm;
    28   double Mo, dMo, M, dM, X2, dS, *chi;
     20  int i, ks, ke;
     21  double Mo, dMo, M, dM, Nm, X2, dS, R, W, *chi;
     22
     23  myAssert (stats->statmode != STATS_NONE, "programming error, liststats mode not set");
     24  myAssert (weight, "work out logic for NULL weight");
    2925
    3026  ke = ks = dMo = 0;
    … …  
    3733  if (N < 1) return (FALSE);
    3834
    39   dsortpair (value, dvalue, N);
    40   stats[0].median = value[(int)(0.5*N)];
    41   stats[0].min    = value[0];
    42   stats[0].max    = value[N-1];
    4335  int N80 = MIN (N-1, 0.8*N);
    4436  int N20 = MAX (0, 0.2*N);
     37
     38  if (weight) {
     39    dsortthree (value, dvalue, weight, N);
     40  } else {
     41    dsortpair (value, dvalue, N);
     42  }
     43
     44  // these values do not depend on the errors or weighting scheme
     45  stats[0].median  = value[(int)(0.5*N)];
     46  stats[0].min     = value[0];
     47  stats[0].max     = value[N-1];
    4548  stats[0].Upper80 = value[N80];
    4649  stats[0].Lower20 = value[N20];
    4750
    48   switch (statmode) {
    49   case M_MEDIAN:
    50     ks = 0;
    51     ke = N;
    52     Mo = stats[0].median;
    53     Nm = N;
    54     goto chisq;
    55     break;
    56   case M_MEAN:
    57   case M_WT_MEAN:
    58     ks = 0;
    59     ke = N;
    60     break;
    61   case M_INNER_MEAN:
    62   case M_INNER_WTMEAN:
    63   case M_CHI_INNER_MEAN:
    64   case M_CHI_INNER_WTMEAN:
    65     ks = 0.25*N + 0.50;
    66     ke = 0.75*N + 0.25;
    67     if (N <= 3) {
     51  switch (stats->statmode) {
     52    case STATS_MEDIAN:
    6853      ks = 0;
    6954      ke = N;
    70     }
    71     break;
     55      Mo = stats[0].median;
     56      Nm = N;
     57      goto chisq;
     58      break;
     59    case STATS_MEAN:
     60    case STATS_WT_MEAN:
     61      ks = 0;
     62      ke = N;
     63      break;
     64    case STATS_INNER_MEAN:
     65    case STATS_INNER_WTMEAN:
     66    case STATS_CHI_INNER_MEAN:
     67    case STATS_CHI_INNER_WTMEAN:
     68      ks = 0.25*N + 0.50;
     69      ke = 0.75*N + 0.25;
     70      if (N <= 3) {
     71        ks = 0;
     72        ke = N;
     73      }
     74      break;
     75    case STATS_NONE:
     76      myAbort ("undefined stats");
    7277  }   
    7378
    74   if ((statmode == M_CHI_INNER_MEAN) || (statmode == M_CHI_INNER_WTMEAN)) {
     79  // for these two modes, I need a vector of the chi-square contribution
     80  // I'm actually just using chisq to get the correct sorting order
     81  if ((stats->statmode == STATS_CHI_INNER_MEAN) || (stats->statmode == STATS_CHI_INNER_WTMEAN)) {
    7582    ALLOCATE (chi, double, N);
    7683    for (i = 0; i < N; i++) {
    7784      chi[i] = (value[i] - stats[0].median) / dvalue[i];
    7885    }
    79     dsortthree (chi, value, dvalue, N);
     86    if (weight) {
     87      dsortfour (chi, value, dvalue, weight, N);
     88    } else {
     89      dsortthree (chi, value, dvalue, N);
     90    }
    8091    free (chi);
    8192  }
    8293
    83   /* calculating the per-star offset based on the weighted average */
    84   M = dM = Nm = 0;
    85   if ((statmode == M_WT_MEAN) || (statmode == M_INNER_WTMEAN) || (statmode == M_CHI_INNER_WTMEAN)) {
    86     for (i = ks; i < ke; i++) {
    87       M   += value[i] / SQ (dvalue[i]);
    88       dM  += 1.0 / SQ (dvalue[i]);
    89       Nm  ++; 
    90     }   
    91     Mo = M / dM;
    92     dMo = sqrt (1.0 / dM);
     94  int WeightedMean = FALSE;
     95  WeightedMean |= (stats->statmode == STATS_WT_MEAN);
     96  WeightedMean |= (stats->statmode == STATS_INNER_WTMEAN);
     97  WeightedMean |= (stats->statmode == STATS_CHI_INNER_WTMEAN);
     98
     99  /* calculating the per-star offset based on the desired weighting scheme */
     100  M = dM = Nm = W = R = 0;
     101  if (weight) {
     102    // the weight value is multiplied by whichever nominal weighting scheme is provided
     103    // thus the user should set weight to 1.0 for nominal weight, and 100 for heavy weight (or so)
     104    // and 0.01 for under-weight
     105    if (WeightedMean) {
     106      for (i = ks; i < ke; i++) {
     107        M  += value[i] * weight[i] / SQ(dvalue[i]);
     108        W  +=            weight[i] / SQ(dvalue[i]);
     109        dM += SQ (weight[i] / dvalue[i]);
     110        R  += weight[i] / SQ(dvalue[i]);
     111        Nm += 1.0; 
     112      }
     113      Mo  = M / W;
     114      dMo = sqrt (dM / R);
     115    } else {
     116      for (i = ks; i < ke; i++) {
     117        M  += value[i] * weight[i];
     118        W  +=            weight[i];
     119        dM += SQ (weight[i] * dvalue[i]);
     120        R  += weight[i];
     121        Nm += 1.0; 
     122      }
     123      Mo  = M / W;
     124      dMo = sqrt (dM / R);
     125    }
    93126  } else {
    94     for (i = ks; i < ke; i++) {
    95       M   += value[i];
    96       dM  += SQ (dvalue[i]);
    97       Nm  ++; 
    98     }   
    99     Mo = M / (double) Nm;
    100     dMo = sqrt (dM / (double) Nm);
     127    // NULL weight vector is supplied, revert to standard form (above reverts if weight[i] == 1)
     128    if (WeightedMean) {
     129      // weighted by inverse-variance
     130      for (i = ks; i < ke; i++) {
     131        M   += value[i] / SQ (dvalue[i]);
     132        dM  += 1.0 / SQ (dvalue[i]);
     133        Nm  += 1.0; 
     134      }
     135      Mo = M / dM;
     136      dMo = sqrt (1.0 / dM);
     137    } else {
     138      // pure un-weighted
     139      for (i = ks; i < ke; i++) {
     140        M   += value[i];
     141        dM  += SQ (dvalue[i]);
     142        Nm  += 1.0; 
     143      }
     144      Mo = M / Nm;
     145      dMo = sqrt (dM) / Nm;
     146    }
    101147  }
    102148
    … …  
    111157  }
    112158  X2 = X2 / (Nm - 1);
    113   dS = sqrt (dS / Nm);
     159  dS = sqrt (dS / (Nm - 1));
    114160
    115161  stats[0].mean  = Mo;
    … …  
    122168}
    123169
     170// we could define the weight to be the only scale factor:
     171// \mu      = \sum (value_i * weight_i) / \sum (weight_i)
     172// \sigma^2 = (1/R) \sum (weight_i^2 \sigma_i^2)
     173// R = \sum (weight_i^2)
     174
     175// or, we could define the weight to be a scale factor times the inverse error:
     176// \mu      = \sum (value_i * weight_i / sigma_i) / \sum (weight_i)
     177// \sigma^2 = (1/R) \sum (weight_i^2 \sigma_i^2)
     178// R = \sum (weight_i^2)
Note: See TracChangeset for help on using the changeset viewer.