IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 41671


Ignore:
Timestamp:
Jun 29, 2021, 11:01:38 AM (5 years ago)
Author:
eugene
Message:

add option to calculate extended variability stats

Location:
trunk/Ohana/src/relphot
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relphot/include/relphot.h

    r41662 r41671  
    163163
    164164typedef enum {
    165   STATS_NONE,
    166   STATS_MEAN,
    167   STATS_MEDIAN,
    168   STATS_WT_MEAN,
    169   STATS_INNER_MEAN,
    170   STATS_INNER_WTMEAN,
    171   STATS_CHI_INNER_MEAN,
    172   STATS_CHI_INNER_WTMEAN
     165  // these modes are primary and mutually exclusive:
     166  STATS_NONE               = 0x0000,
     167  STATS_MEAN               = 0x0001,
     168  STATS_MEDIAN             = 0x0002,
     169  STATS_WT_MEAN            = 0x0003,
     170  STATS_INNER_MEAN         = 0x0004,
     171  STATS_INNER_WTMEAN       = 0x0005,
     172  STATS_CHI_INNER_MEAN     = 0x0006,
     173  STATS_CHI_INNER_WTMEAN   = 0x0007,
     174  // these modes are additional options
     175  STATS_PRIMARY            = 0x0007, // use this to mask out the optional bits
     176  STATS_VARSTATS           = 0x0010,
    173177} ListStatsMode;
    174178
    … …  
    185189  double Lower20;
    186190  double Lower10;
     191  double Upper90Nsig;
     192  double Upper80Nsig;
     193  double Lower20Nsig;
     194  double Lower10Nsig;
    187195  double total;
    188196  int    Nmeas;
    … …  
    419427double IMFIT_SYS_SIGMA_LIM;
    420428double CLOUD_TOLERANCE;
     429int    VARIABILITY_STATS;
    421430
    422431int    PARALLEL;

  • trunk/Ohana/src/relphot/src/args.c

    r41648 r41671  
    434434    remove_argument (N, &argc, argv);
    435435    KEEP_UBERCAL = FALSE;
     436  }
     437
     438  VARIABILITY_STATS = FALSE;
     439  if ((N = get_argument (argc, argv, "-varstats"))) {
     440    remove_argument (N, &argc, argv);
     441    VARIABILITY_STATS = TRUE;
    436442  }
    437443
    … …  
    865871
    866872
     873  VARIABILITY_STATS = FALSE;
     874  if ((N = get_argument (argc, argv, "-varstats"))) {
     875    remove_argument (N, &argc, argv);
     876    VARIABILITY_STATS = TRUE;
     877  }
     878
    867879  MIN_ERROR = 0.001;
    868880  if ((N = get_argument (argc, argv, "-minerror"))) {

  • trunk/Ohana/src/relphot/src/liststats.c

    r41670 r41671  
    33void liststats_setmode (StatType *stats, char *strmode) {
    44
    5   stats->statmode = -1;
     5  stats->statmode = STATS_NONE;
    66  if (!strcmp (strmode, "MEAN"))             { stats->statmode = STATS_MEAN; return; }
    77  if (!strcmp (strmode, "MEDIAN"))           { stats->statmode = STATS_MEDIAN; return; }
    … …  
    1717
    1818int liststats_init (StatType *stats) {
    19   stats->median  = NAN;
    20   stats->mean    = NAN;
    21   stats->sigma   = NAN;
    22   stats->error   = NAN;
    23   stats->chisq   = NAN;
    24   stats->min     = NAN;
    25   stats->max     = NAN;
    26   stats->Upper80 = NAN;
    27   stats->Lower20 = NAN;
    28   stats->Upper90 = NAN;
    29   stats->Lower10 = NAN;
    30   stats->total   = NAN;
    31   stats->Nmeas = 0;
     19  stats->median      = NAN;
     20  stats->mean        = NAN;
     21  stats->sigma       = NAN;
     22  stats->error       = NAN;
     23  stats->chisq       = NAN;
     24  stats->min         = NAN;
     25  stats->max         = NAN;
     26  stats->Upper80     = NAN;
     27  stats->Lower20     = NAN;
     28  stats->Upper90     = NAN;
     29  stats->Lower10     = NAN;
     30  stats->Upper80Nsig = NAN;
     31  stats->Lower20Nsig = NAN;
     32  stats->Upper90Nsig = NAN;
     33  stats->Lower10Nsig = NAN;
     34  stats->total       = NAN;
     35  stats->Nmeas       = 0;
    3236  return TRUE;
    3337}
    … …  
    3943
    4044  myAssert (stats->statmode != STATS_NONE, "programming error, liststats mode not set");
     45  ListStatsMode myMode = stats->statmode & STATS_PRIMARY; // exclude the option bits
    4146
    4247  liststats_init (stats);
    … …  
    7580  stats[0].Lower10 = value[N10];
    7681
    77   switch (stats->statmode) {
     82  switch (myMode) {
    7883    case STATS_MEDIAN:
    7984      ks = 0;
    … …  
    100105      break;
    101106    case STATS_NONE:
     107    default:
    102108      myAbort ("undefined stats");
    103109  }   
    … …  
    105111  // for these two modes, I need a vector of the chi-square contribution
    106112  // I'm actually just using chisq to get the correct sorting order
    107   if ((stats->statmode == STATS_CHI_INNER_MEAN) || (stats->statmode == STATS_CHI_INNER_WTMEAN)) {
     113  if ((myMode == STATS_CHI_INNER_MEAN) || (myMode == STATS_CHI_INNER_WTMEAN)) {
    108114    ALLOCATE (chi, double, N);
    109115    for (i = 0; i < N; i++) {
    … …  
    119125
    120126  int WeightedMean = FALSE;
    121   WeightedMean |= (stats->statmode == STATS_WT_MEAN);
    122   WeightedMean |= (stats->statmode == STATS_INNER_WTMEAN);
    123   WeightedMean |= (stats->statmode == STATS_CHI_INNER_WTMEAN);
     127  WeightedMean |= (myMode == STATS_WT_MEAN);
     128  WeightedMean |= (myMode == STATS_INNER_WTMEAN);
     129  WeightedMean |= (myMode == STATS_CHI_INNER_WTMEAN);
    124130  if (!dvalue) WeightedMean = FALSE;  // warn the user?
    125131
    … …  
    310316     
    311317  // calculate the weight thresholds to mask the bad points:
    312   // double Sum_W = 0.0;
    313318  for (int i = 0; i < Npoints; i++) {
    314319    dataset->wtvals[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
    315320    dataset->wtlist[i] = dataset->wtvals[i];
    316     // Sum_W += dataset->wtvals[i];
    317321  }
    318322  dsort (dataset->wtlist, Npoints);
    319323  double WtMedian = (Npoints % 2) ? dataset->wtlist[midpt] : 0.5*(dataset->wtlist[midpt] + dataset->wtlist[midpt-1]);
    320   // double WtThreshold = WEIGHT_THRESHOLD * Sum_W / (1.0 * Npoints);
    321324  double WtThreshold = WEIGHT_THRESHOLD * WtMedian;
    322325
    … …  
    348351  stats->sigma = sqrt (dSig / (Nkeep - 1));
    349352
     353  // if percentile ranges are desired, calculate them
     354  if (stats->statmode & STATS_VARSTATS) {
     355    // save ykeep values in a vector to be sorted
     356    for (int i = 0; i < Nkeep; i++) {
     357      dataset->wtvals[i] = dataset->ykeep[i];
     358    }
     359    dsort (dataset->wtvals, Nkeep);
     360
     361    int N90 = MIN (Nkeep-1, 0.9*Nkeep);
     362    int N80 = MIN (Nkeep-1, 0.8*Nkeep);
     363    int N20 = MAX (0, 0.2*Nkeep);
     364    int N10 = MAX (0, 0.1*Nkeep);
     365
     366    stats->Upper90 = dataset->wtvals[N90];
     367    stats->Upper80 = dataset->wtvals[N80];
     368    stats->Lower20 = dataset->wtvals[N20];
     369    stats->Lower10 = dataset->wtvals[N10];
     370   
     371    int NkeepMid = 0.5 * Nkeep;
     372    double ClipMedian = (Nkeep % 2) ? dataset->wtlist[NkeepMid] : 0.5*(dataset->wtlist[NkeepMid] + dataset->wtlist[NkeepMid-1]);
     373
     374    // save ykeep values in a vector to be sorted
     375    for (int i = 0; i < Nkeep; i++) {
     376      dataset->wtvals[i] = (dataset->ykeep[i] - ClipMedian) / dataset->dykeep[i];
     377    }
     378    dsort (dataset->wtvals, Nkeep);
     379
     380    stats->Upper90Nsig = dataset->wtvals[N90];
     381    stats->Upper80Nsig = dataset->wtvals[N80];
     382    stats->Lower20Nsig = dataset->wtvals[N20];
     383    stats->Lower10Nsig = dataset->wtvals[N10];
     384  }
     385
    350386  int Nboot = 0;
    351387  for (int iboot = 0; iboot < NBOOTSTRAP; iboot++) {

  • trunk/Ohana/src/relphot/src/reload_catalogs.c

    r41647 r41671  
    237237    if (USE_BASIC_CHECK)   {  strextend (&command, "-basic-image-search"); }
    238238    // if (USE_ALL_IMAGES)    { strextend (&command, "-use-all-images"); }
     239    if (VARIABILITY_STATS) { strextend (&command, "-varstats"); }
    239240    if (USE_MCAL_PSF_FOR_STACK_APER) { strextend (&command, "-use-mcal-psf-for-stack-aper"); }
    240241

  • trunk/Ohana/src/relphot/src/setMrelCatalog.c

    r41647 r41671  
    11# include "relphot.h"
     2
     3# define UPPER90_VAL MpsfStk
     4# define UPPER80_VAL FpsfStk
     5# define LOWER20_VAL MkronStk
     6# define LOWER10_VAL FkronStk
     7# define UPPER90_SIG MpsfWrp
     8# define UPPER80_SIG FpsfWrp
     9# define LOWER20_SIG MkronWrp
     10# define LOWER10_SIG FkronWrp
     11
    212
    313# if (0)
    … …  
    435445    StatType *psfstats = &results->psfstats;
    436446
     447    // if isSetMrelFinal and varstats are requested, request percentile ranges
     448    if (isSetMrelFinal && VARIABILITY_STATS) psfstats->statmode |= STATS_VARSTATS;
     449
    437450    int NrankingPSF = useOLS ? magStatsByRanking (&results->psfData[Nsec], psfstats) : magStatsByRankingIRLS (&results->psfData[Nsec], psfstats);
    438451    if (NrankingPSF < Nminmeas) {
    … …  
    469482        secfilt[Nsec].Mmax     = psfstats->max;
    470483        secfilt[Nsec].Mmin     = psfstats->min;
     484      }
     485      if (VARIABILITY_STATS) {
     486        // XXX NOTE: these are overloaded on warp and stack mags and fluxes
     487        secfilt[Nsec].UPPER90_VAL = psfstats->Upper90;
     488        secfilt[Nsec].UPPER80_VAL = psfstats->Upper80;
     489        secfilt[Nsec].LOWER20_VAL = psfstats->Lower20;
     490        secfilt[Nsec].LOWER10_VAL = psfstats->Lower10;
     491        secfilt[Nsec].UPPER90_SIG = psfstats->Upper90Nsig;
     492        secfilt[Nsec].UPPER80_SIG = psfstats->Upper80Nsig;
     493        secfilt[Nsec].LOWER20_SIG = psfstats->Lower20Nsig;
     494        secfilt[Nsec].LOWER10_SIG = psfstats->Lower10Nsig;
    471495      }
    472496      secfilt[Nsec].psfQfMax     = results->psfQfMax[Nsec];
Note: See TracChangeset for help on using the changeset viewer.