IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39636 for trunk


Ignore:
Timestamp:
Jul 23, 2016, 12:46:44 PM (10 years ago)
Author:
eugene
Message:

add liststats_irls, magStatsByRankingIRLS; add code to set mask bit for objects masked by the IRLS fitting; add irls test code; use same method for chip and warp

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

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relphot/Makefile

    r39513 r39636  
    2020relphot_client: $(BIN)/relphot_client.$(ARCH)
    2121
     22test_liststats: $(BIN)/test_liststats.$(ARCH)
     23
    2224install: $(DESTBIN)/relphot $(DESTBIN)/relphot_client
    2325
    2426RELPHOT = \
     27$(SRC)/StatDataSetOps.$(ARCH).o  \
    2528$(SRC)/ConfigInit.$(ARCH).o      \
    2629$(SRC)/GridOps.$(ARCH).o         \
    … …  
    7982
    8083RELPHOT_CLIENT = \
     84$(SRC)/StatDataSetOps.$(ARCH).o  \
    8185$(SRC)/ConfigInit.$(ARCH).o      \
    8286$(SRC)/GridOps.$(ARCH).o         \
    … …  
    120124$(RELPHOT_CLIENT): $(INC)/relphot.h
    121125$(BIN)/relphot_client.$(ARCH): $(RELPHOT_CLIENT)
     126
     127TEST_LISTSTATS = \
     128$(SRC)/StatDataSetOps.$(ARCH).o  \
     129$(SRC)/liststats.$(ARCH).o       \
     130$(SRC)/test_liststats.$(ARCH).o
     131
     132$(TEST_LISTSTATS): $(INC)/relphot.h
     133$(BIN)/test_liststats.$(ARCH): $(TEST_LISTSTATS)
     134

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

    r39517 r39636  
    119119  int    *ranking;            // weights to use for mean mags
    120120  int    *measSeq;            // weights to use for mean mags
     121  int    *msklist;            // mask modifications
    121122  int     Nlist;
    122123} StatDataSet;
    … …  
    428429int           liststats           PROTO((double *value, double *dvalue, double *wvalue, int N, StatType *stats));
    429430int           liststats_init      PROTO((StatType *stats));
     431int           liststats_irls      PROTO((StatDataSet *dataset, int Npoints, StatType *stats));
     432
    430433Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int hostID, char *hostpath, char *syncfile));
    431434Catalog      *load_catalogs_parallel PROTO((SkyList *sky, int *Ncatalog, char *syncfile));

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

    r39522 r39636  
    5252  value = catalog[cat].secfilt[Nsecfilt*ave+Nsec].M;
    5353  return (value);
    54 }
    55 
    56 StatDataSet *StatDataSetAlloc (int Nsecfilt, int Nmax) {
    57 
    58   int i;
    59 
    60   StatDataSet *dataset = NULL;
    61   ALLOCATE (dataset, StatDataSet, Nsecfilt);
    62   for (i = 0; i < Nsecfilt; i++) {
    63     ALLOCATE (dataset[i].flxlist, double, Nmax);
    64     ALLOCATE (dataset[i].wgtlist, double, Nmax);
    65     ALLOCATE (dataset[i].errlist, double, Nmax);
    66     ALLOCATE (dataset[i].ranking,    int, Nmax);
    67     ALLOCATE (dataset[i].measSeq,    int, Nmax);
    68   } 
    69   return dataset;
    70 }
    71 
    72 void StatDataSetFree (StatDataSet *dataset, int Nsecfilt) {
    73 
    74   int i;
    75 
    76   for (i = 0; i < Nsecfilt; i++) {
    77     FREE (dataset[i].flxlist);
    78     FREE (dataset[i].wgtlist);
    79     FREE (dataset[i].errlist);
    80     FREE (dataset[i].ranking);
    81     FREE (dataset[i].measSeq);
    82   } 
    83   FREE (dataset);
    8454}
    8555

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

    r38466 r39636  
    188188}
    189189
     190// These should probably be tunable:
     191# define MAX_ITERATIONS 10
     192# define FIT_TOLERANCE 1e-4
     193# define FLT_TOLERANCE 1e-6
     194# define WEIGHT_THRESHOLD 0.3
     195# define NBOOTSTRAP 100
     196
     197int fit_least_squares (double *fit, double *y, double *dy, double *wgt, double *wt, int Npts);
     198double VectorFractionInterpolate (double *values, float fraction, int Npts);
     199double weight_cauchy (double x);
     200
     201// this is a zero-order fit (constant value only)
     202int liststats_irls (StatDataSet *dataset, int Npoints, StatType *stats) {
     203
     204  double value;
     205
     206  liststats_init (stats);
     207
     208  // OLS
     209  if (!fit_least_squares (&value, dataset->flxlist, dataset->errlist, dataset->wgtlist, NULL, Npoints)) return FALSE;
     210 
     211  // XXX add to dataset elements?
     212  ALLOCATE_PTR (wt, double, Npoints);
     213
     214  int converged = FALSE;
     215  for (int iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations++) {
     216   
     217    for (int i = 0; i < Npoints; i++) {
     218      // we are only including the formal error, not the weight in the definition of wt[]
     219      wt[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
     220    }
     221   
     222    double oldValue = value;
     223    if (!fit_least_squares (&value, dataset->flxlist, dataset->errlist, dataset->wgtlist, wt, Npoints)) {
     224      value = oldValue;
     225      break;
     226    }
     227   
     228    converged = TRUE;
     229    if ((fabs(value - oldValue) > FIT_TOLERANCE * fabs(value)) &&
     230        (fabs(value - oldValue) > FLT_TOLERANCE)) {
     231      converged = FALSE;
     232    }
     233  }
     234  stats->mean = value;
     235     
     236  // calculate the weight thresholds to mask the bad points:
     237  double Sum_W = 0.0;
     238  for (int i = 0; i < Npoints; i++) {
     239    wt[i] = weight_cauchy ((dataset->flxlist[i] - value) / dataset->errlist[i]);
     240    Sum_W += wt[i];
     241  }
     242  double WtThreshold = WEIGHT_THRESHOLD * Sum_W / (1.0 * Npoints);
     243
     244  // generate the unmasked subset
     245  // XXX add these to the dataset elements?
     246  ALLOCATE_PTR ( ykeep, double, Npoints);
     247  ALLOCATE_PTR (dykeep, double, Npoints);
     248  ALLOCATE_PTR (wtkeep, double, Npoints);
     249
     250  // save unmasked points
     251  int Nkeep = 0;
     252  stats->min = NAN;
     253  stats->max = NAN;
     254  double dChi = 0.0;
     255  double dSig = 0.0;
     256  for (int i = 0; i < Npoints; i++) {
     257    if ((wt[i] < WtThreshold) || !isfinite(dataset->flxlist[i])) {
     258      dataset->msklist[i] = TRUE; // mark the masked points
     259      continue;
     260    }
     261     ykeep[Nkeep] = dataset->flxlist[i];
     262    dykeep[Nkeep] = dataset->errlist[i];
     263    wtkeep[Nkeep] = dataset->wgtlist[i]; // externally-supplied weight
     264    Nkeep ++;
     265   
     266    // record mean, error, chisq, min, max, sigma, Nmeas (unmasked points)
     267    stats->min = isfinite(stats->min) ? MIN(dataset->flxlist[i], stats->min) : dataset->flxlist[i];
     268    stats->max = isfinite(stats->max) ? MAX(dataset->flxlist[i], stats->max) : dataset->flxlist[i];
     269    double dValue2 = SQ(dataset->flxlist[i] - value);
     270    dChi += dValue2 / SQ (dataset->errlist[i]);
     271    dSig += dValue2;
     272  }
     273  stats->Nmeas = Nkeep;
     274  stats->chisq = dChi / (Nkeep - 1);
     275  stats->sigma = sqrt (dSig / (Nkeep - 1));
     276
     277  // bootstrap resampling to generate the errorbars
     278  // XXX add these to the dataset elements?
     279  ALLOCATE_PTR (ysample,  double, Nkeep);
     280  ALLOCATE_PTR (dysample, double, Nkeep);
     281  ALLOCATE_PTR (wtsample, double, Nkeep);
     282  ALLOCATE_PTR (bvalue,   double, NBOOTSTRAP); // vector to save the bootstrap values
     283
     284  int Nboot = 0;
     285  for (int iboot = 0; iboot < NBOOTSTRAP; iboot++) {
     286   
     287    // resample
     288    for (int i = 0; i < Nkeep; i++) {
     289      // I need to draw Npoints random entries from 'points' with replacement:
     290      int N = Nkeep * drand48();
     291       ysample[i] =  ykeep[N];
     292      dysample[i] = dykeep[N];
     293      wtsample[i] = wtkeep[N];
     294    }
     295
     296    if (!fit_least_squares (&value, ysample, dysample, wtsample, NULL, Nkeep)) continue;
     297
     298    bvalue[Nboot] = value;
     299    Nboot ++;
     300  }
     301
     302  dsort (bvalue, Nboot);
     303 
     304  double Slo = VectorFractionInterpolate (bvalue, 0.158655, Nboot);
     305  double Shi = VectorFractionInterpolate (bvalue, 0.841345, Nboot);
     306  stats->error = (Shi - Slo) / 2.0;
     307
     308  free (bvalue);
     309  free ( ysample);
     310  free (dysample);
     311  free (wtsample);
     312  free ( ykeep);
     313  free (dykeep);
     314  free (wtkeep);
     315  free (wt);
     316
     317  return TRUE;
     318}
     319
     320// wgt is externally-supplied weight, wt is optional
     321int fit_least_squares (double *fit, double *y, double *dy, double *wgt, double *wt, int Npts) {
     322
     323  int i;
     324
     325  double S0 = 0;
     326  double S1 = 0;
     327
     328  /* perform linear fit */
     329  for (i = 0; i < Npts; i++, y++, dy++, wgt++) {
     330    if (!finite(*y)) continue;
     331
     332    // wt is optional
     333    double dY = wt ? wt[i] * (*wgt) / SQ(*dy) : (*wgt) / SQ(*dy);
     334
     335    S0 +=    dY;
     336    S1 += *y*dY;
     337  }
     338  if (S0 == 0.0) return FALSE;
     339  *fit = S1 / S0;
     340  return TRUE;
     341}
     342
     343double weight_cauchy (double x) {
     344  double r = x / 2.385;
     345  return (1.0 / (1.0 + SQ(r)));
     346}
     347
     348double VectorFractionInterpolate (double *values, float fraction, int Npts) {
     349
     350  float F = fraction * Npts;
     351  int   N = fraction * Npts;
     352
     353  if (N < 0        ) return NAN;
     354  if (N >= Npts - 2) return NAN;
     355
     356  // interpolate between N,N+1
     357   
     358  double S = (F - N) * (values[N+1] - values[N]) + values[N];
     359  return S;
     360}
     361
    190362// we could define the weight to be the only scale factor:
    191363// \mu      = \sum (value_i * weight_i) / \sum (weight_i)

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

    r39635 r39636  
    1010
    1111int magStatsByRankingClipped (StatDataSet *dataset, StatType *stats);
    12 int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags);
     12int magStatsByRankingIRLS (StatDataSet *dataset, StatType *stats);
     13int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags keepflag, DVOMeasureFlags maskflag);
    1314void GetPhotFlagStats (uint32_t *photFlagUpper, uint32_t *photFlagLower, uint32_t *photflag_list, int Nphotflag);
    1415void sort_entry_by_offset (double *offset, int *entry, int N);
    1516void sort_StatDataSet (StatDataSet *dataset);
     17int liststats_irls (StatDataSet *dataset, int Npoints, StatType *stats);
     18
     19# define MAG_STATS_BY_RANKING magStatsByRankingIRLS
    1620
    1721# define UBERCAL_WEIGHT 100.0
    … …  
    255259        results->aperData[Nsec].ranking[Nap] = measureRank[k];
    256260        results->aperData[Nsec].measSeq[Nap] = k;
     261        results->aperData[Nsec].msklist[Nap] = 0;
    257262        results->aperData[Nsec].Nlist ++;
    258263      }
    … …  
    270275        results->kronData[Nsec].ranking[Nkron] = measureRank[k];
    271276        results->kronData[Nsec].measSeq[Nkron] = k;
     277        results->kronData[Nsec].msklist[Nkron] = 0;
    272278        results->kronData[Nsec].Nlist ++;
    273279      }
    … …  
    283289      results->psfData[Nsec].ranking[Npsf] = measureRank ? measureRank[k] : 0;
    284290      results->psfData[Nsec].measSeq[Npsf] = k;
     291      results->psfData[Nsec].msklist[Npsf] = 0;
    285292      results->psfData[Nsec].Nlist ++;
    286293    }
    … …  
    362369    // if too few valid measurements meet the minimum criteria, go to the next entry
    363370    StatType *psfstats = &results->psfstats;
    364     int Nranking = magStatsByRanking (&results->psfData[Nsec], psfstats);
     371    int Nranking = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats);
    365372    if (Nranking < Nminmeas) {
    366373      secfilt[Nsec].flags |= ID_OBJ_FEW;
    … …  
    370377      secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
    371378    }
    372     int minRank = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;
     379    int minRankPSF = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;
    373380
    374381    // when running -averages, we have no information about the images, so we cannot set this
    … …  
    383390
    384391      // mark the measurements matching this ranking
    385       markMeasureByRanking (&results->psfData[Nsec], measure, minRank, ID_MEAS_PHOTOM_PSF);
     392      markMeasureByRanking (&results->psfData[Nsec], measure, minRankPSF, ID_MEAS_PHOTOM_PSF, ID_MEAS_MASKED_PSF);
    386393
    387394      if (Nranking) {
    … …  
    391398        secfilt[Nsec].Mmin   = psfstats->min;
    392399      }
    393 
    394400      secfilt[Nsec].psfQfMax     = results->psfQfMax[Nsec];
    395401      secfilt[Nsec].psfQfPerfMax = results->psfQfPerfMax[Nsec];
    … …  
    399405
    400406      StatType *apstats = &results->apstats;
    401       Nranking = magStatsByRanking (&results->aperData[Nsec], apstats);
     407      Nranking = MAG_STATS_BY_RANKING (&results->aperData[Nsec], apstats);
    402408      if (Nranking) {
    403409        secfilt[Nsec].Map     = apstats->mean;
    … …  
    406412        secfilt[Nsec].NusedAp = Nranking;
    407413      }
    408       markMeasureByRanking (&results->aperData[Nsec], measure, minRank, ID_MEAS_PHOTOM_APER);
     414      int minRankAper = (Nranking > 0) ? results->aperData[Nsec].ranking[0] : 10;
     415      markMeasureByRanking (&results->aperData[Nsec], measure, minRankAper, ID_MEAS_PHOTOM_APER, ID_MEAS_MASKED_APER);
    409416
    410417      StatType *kronstats = &results->kronstats;
    411       Nranking = magStatsByRanking (&results->kronData[Nsec], kronstats);
     418      Nranking = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats);
    412419      if (Nranking) {
    413420        secfilt[Nsec].Mkron     = kronstats->mean;
    … …  
    416423        secfilt[Nsec].NusedKron = Nranking;
    417424      }
    418       markMeasureByRanking (&results->kronData[Nsec], measure, minRank, ID_MEAS_PHOTOM_KRON);
     425      int minRankKron = (Nranking > 0) ? results->kronData[Nsec].ranking[0] : 10;
     426      markMeasureByRanking (&results->kronData[Nsec], measure, minRankKron, ID_MEAS_PHOTOM_KRON, ID_MEAS_MASKED_KRON);
    419427
    420428      // does this object appear extended in > 50% of measurements?
    … …  
    423431      }
    424432
    425       switch (minRank) {
     433      switch (minRankPSF) {
    426434        case 0:
    427435          secfilt[Nsec].flags |= ID_SECF_RANK_0;
    … …  
    853861      results->psfData[Nsec].ranking[Npsf] = measureRank ? measureRank[k] : 0;
    854862      results->psfData[Nsec].measSeq[Npsf] = k;
     863      results->psfData[Nsec].msklist[Npsf] = 0;
    855864      results->psfData[Nsec].Nlist ++;
    856865    }     
    … …  
    871880      results->aperData[Nsec].ranking[Naper] = measureRank ? measureRank[k] : 0;
    872881      results->aperData[Nsec].measSeq[Naper] = k;
     882      results->aperData[Nsec].msklist[Naper] = 0;
    873883      results->aperData[Nsec].Nlist ++;
    874884    }
    … …  
    884894      results->kronData[Nsec].ranking[Nkron] = measureRank ? measureRank[k] : 0;
    885895      results->kronData[Nsec].measSeq[Nkron] = k;
     896      results->kronData[Nsec].msklist[Nkron] = 0;
    886897      results->kronData[Nsec].Nlist ++;
    887898    }
    … …  
    905916    }
    906917
     918    int Nranking, minRank;
     919
    907920    // if too few valid measurements meet the minimum criteria, go to the next entry
    908921    StatType *psfstats = &results->psfstats;
    909     int Nranking = magStatsByRankingClipped (&results->psfData[Nsec], psfstats);
     922    Nranking = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats);
    910923    if (Nranking) {
    911924      secfilt[Nsec].FpsfWrp  = psfstats->mean;
    … …  
    915928      secfilt[Nsec].MpsfWrp  = isnan(secfilt[Nsec].FpsfWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FpsfWrp); // 8.9 since flux is in Jy
    916929    }
    917     int minRank = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;
    918     markMeasureByRanking (&results->psfData[Nsec], measure, minRank, ID_MEAS_WARP_USED);
     930    minRank = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;
     931    markMeasureByRanking (&results->psfData[Nsec], measure, minRank, ID_MEAS_WARP_USED | ID_MEAS_PHOTOM_PSF, ID_MEAS_MASKED_PSF);
    919932
    920933    // if too few valid measurements meet the minimum criteria, go to the next entry
    921934    StatType *apstats = &results->apstats;
    922     Nranking = magStatsByRankingClipped (&results->aperData[Nsec], apstats);
     935    Nranking = MAG_STATS_BY_RANKING (&results->aperData[Nsec], apstats);
    923936    if (Nranking) {
    924937      secfilt[Nsec].FapWrp     = apstats->mean;
    … …  
    928941      secfilt[Nsec].MapWrp     = isnan(secfilt[Nsec].FapWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FapWrp); // 8.9 since flux is in Jy
    929942    }
     943    minRank = (Nranking > 0) ? results->aperData[Nsec].ranking[0] : 10;
     944    markMeasureByRanking (&results->aperData[Nsec], measure, minRank, ID_MEAS_WARP_USED | ID_MEAS_PHOTOM_APER, ID_MEAS_MASKED_APER);
    930945
    931946    // if too few valid measurements meet the minimum criteria, go to the next entry
    932947    StatType *kronstats = &results->kronstats;
    933     Nranking = magStatsByRankingClipped (&results->kronData[Nsec], kronstats);
     948    Nranking = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats);
    934949    if (Nranking) {
    935950      secfilt[Nsec].FkronWrp     = kronstats->mean;
    … …  
    939954      secfilt[Nsec].MkronWrp     = isnan(secfilt[Nsec].FkronWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FkronWrp); // 8.9 since flux is in Jy
    940955    }
     956    minRank = (Nranking > 0) ? results->kronData[Nsec].ranking[0] : 10;
     957    markMeasureByRanking (&results->kronData[Nsec], measure, minRank, ID_MEAS_WARP_USED | ID_MEAS_PHOTOM_KRON, ID_MEAS_MASKED_KRON);
    941958
    942959    secfilt[Nsec].Nwarp     = results->Nmeas[Nsec];
    … …  
    967984}
    968985
    969 # if (0)
    970986int magStatsByRankingIRLS (StatDataSet *dataset, StatType *stats) {
    971987
    … …  
    9851001  for (i = 0; (i < dataset->Nlist) && (dataset->ranking[i] == minRank); i++, Nranking++);
    9861002
    987   liststats (dataset->flxlist, dataset->errlist, dataset->wgtlist, Nranking, stats);
     1003  liststats_irls (dataset, Nranking, stats);
    9881004  return (Nranking);
    9891005}
    990 # endif
    9911006
    9921007// outlier warp measurements are driving bad mean values. 
    … …  
    10611076}
    10621077
    1063 int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags) {
     1078int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags keepflag, DVOMeasureFlags maskflag) {
    10641079  int i;
    10651080
    … …  
    10671082    int k = dataset->measSeq[i];
    10681083    if (dataset->ranking[i] > minrank) {
    1069       measure[k].dbFlags &= ~flags;
     1084      measure[k].dbFlags &= ~keepflag;
    10701085    } else {
    1071       measure[k].dbFlags |= flags;
    1072     }
     1086      measure[k].dbFlags |= keepflag;
     1087    }
     1088    if (dataset->msklist[i]) measure[k].dbFlags |= maskflag;
    10731089  }
    10741090  return TRUE;
Note: See TracChangeset for help on using the changeset viewer.