IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38574


Ignore:
Timestamp:
Jul 6, 2015, 8:11:18 PM (11 years ago)
Author:
eugene
Message:

add clipping and bootstrap-resampling measurement of the errors

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20150625/Ohana/src/opihi/cmd.astro/fitplx.c

    r38568 r38574  
    1414  double *dX;
    1515  double *dY;
     16  int *index;
    1617  int Npts;
    1718} PlxFitData;
    … …  
    3132} PlxFit;
    3233
     34int VectorRobustStats (Vector *vector, double *median, double *sigma);
     35double VectorFractionInterpolate (double *values, float fraction, int Npts);
     36
     37int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal);
     38int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean);
     39int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE);
     40
    3341int PlxFitDataAlloc (PlxFitData *data, int N);
    3442void PlxFitDataFree (PlxFitData *data);
    … …  
    7078    remove_argument (N, &argc, argv);
    7179    dPsigMax = atof(argv[N]);
     80    remove_argument (N, &argc, argv);
     81  }
     82
     83  int Nresample = 0;
     84  if ((N = get_argument (argc, argv, "-bootstrap-resample"))) {
     85    remove_argument (N, &argc, argv);
     86    Nresample = atoi(argv[N]);
    7287    remove_argument (N, &argc, argv);
    7388  }
    … …  
    112127  }
    113128
    114   int Ntotal = tvec->Nelements; // XXX check other lengths
    115 
    116   // find mean values to remove
    117   double Nmean = 0;
    118   double Tmean = 0;
    119   double Rmean = 0;
    120   double Dmean = 0;
    121   double Tmin = +1000000;
    122   double Tmax = -1000000;
    123   for (i = 0; i < Ntotal; i++) {
    124     if (mask && !mask[i]) continue;
    125     Rmean += R[i];
    126     Dmean += D[i];
    127     Tmean += T[i];
    128     Tmin = MIN(Tmin, T[i]);
    129     Tmax = MAX(Tmax, T[i]);
    130     Nmean += 1.0;
    131   }
    132   Rmean /= Nmean;
    133   Dmean /= Nmean;
    134   Tmean /= Nmean;
    135 
    136129  // Ntotal : all points supplied by user
    137130  // Nsubset : unmasked points
    138   // Nmean : unmasked points (a double only used above)
    139 
    140   float Trange = Tmax - Tmin;
    141   // fprintf (stderr, "R,D : %f,%f, T: %f, Trange: %f, Tmin: %f, Tmax: %f\n", Rmean, Dmean, Tmean, Trange, Tmin, Tmax);
     131  int Ntotal = tvec->Nelements; // XXX check other lengths
     132  if (dPvec) ResetVector (dPvec, OPIHI_FLT, Ntotal);
     133
     134  double Rmean, Dmean, Tmean;
     135  PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
    142136
    143137  /* project coordinates to a plane centered on the object with units of arcsec */
    … …  
    148142  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
    149143
    150   float pXmin = +2.0;
    151   float pXmax = -2.0;
    152   float pYmin = +2.0;
    153   float pYmax = -2.0;
    154 
    155   PlxFit fit; memset (&fit, 0, sizeof(PlxFit));
    156144  PlxFitData fitdata;
    157145  PlxFitDataAlloc (&fitdata, Ntotal);
    158 
    159   // generate the fit values (projected X,Y; parallax factors;
    160 
    161   // save an index so we can supplied dPsig for the unmasked points
    162   int *index; 
    163   ALLOCATE (index, int, Ntotal);
    164  
    165   int Nsubset = 0;
    166   for (i = 0; i < Ntotal; i++) {
    167     if (mask && !mask[i]) continue;
    168     RD_to_XY (&fitdata.X[Nsubset], &fitdata.Y[Nsubset], R[i], D[i], &coords);
    169     fitdata.dX[Nsubset] = dR[i];
    170     fitdata.dY[Nsubset] = dD[i];
    171     fitdata.t[Nsubset] = (T[i] - Tmean) / 365.25;
    172     ParFactor (&fitdata.pX[Nsubset], &fitdata.pY[Nsubset], R[i], D[i], T[i]);
    173     pXmin = MIN (pXmin, fitdata.pX[Nsubset]);
    174     pXmax = MAX (pXmax, fitdata.pX[Nsubset]);
    175     pYmin = MIN (pYmin, fitdata.pY[Nsubset]);
    176     pYmax = MAX (pYmax, fitdata.pY[Nsubset]);
    177     index[Nsubset] = i;
    178     Nsubset++;
    179   }
    180   fitdata.Npts = Nsubset;
    181   float dXRange = pXmax - pXmin;
    182   float dYRange = pYmax - pYmin;
    183   float parRange = hypot (dXRange, dYRange);
    184        
    185   // fprintf (stderr, "par factor range: %f\n", parRange);
     146  PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
     147
     148  PlxFit fit; memset (&fit, 0, sizeof(PlxFit));
    186149
    187150  // determine dPsig for detections based on Noutlier attempts
    188151  if (Noutlier) {
    189     PlxFit testfit;
    190     PlxFitData sample;
    191     PlxFitDataAlloc (&sample, fitdata.Npts);
    192 
    193     double **dXsig, **dYsig;
    194     ALLOCATE (dXsig, double *, fitdata.Npts);
    195     ALLOCATE (dYsig, double *, fitdata.Npts);
    196     for (i = 0; i < fitdata.Npts; i++) {
    197       ALLOCATE (dXsig[i], double, Noutlier);
    198       ALLOCATE (dYsig[i], double, Noutlier);
    199     }
    200 
    201     testfit.getChisq = FALSE;
    202 
    203     int n;
    204     int Nsamples = 0;
    205     for (n = 0; n < Noutlier; n++) {
    206       // bootstrap resample (fitdata -> sample)
    207       PlxBootstrapResample (&fitdata, &sample);
    208      
    209       if (n % 1000 == 999) fprintf (stderr, ".");
    210 
    211       // fit the sample
    212       if (!FitPMandPar (&testfit,
    213                         sample.X, sample.dX,
    214                         sample.Y, sample.dY, sample.t,
    215                         sample.pX, sample.pY, sample.Npts, VERBOSE)) continue;
    216 
    217       // fprintf (stderr, "%f +/- %f | %f %f\n", testfit.p, testfit.dp, testfit.uR, testfit.uD);
    218 
    219       // find the distances to the path
    220       for (i = 0; i < fitdata.Npts; i++) {
    221         double Xf = testfit.Ro + testfit.uR*fitdata.t[i] + testfit.p*fitdata.pX[i];
    222         double Yf = testfit.Do + testfit.uD*fitdata.t[i] + testfit.p*fitdata.pY[i];
    223         dXsig[i][Nsamples] = fabs(fitdata.X[i] - Xf) / fitdata.dX[i];
    224         dYsig[i][Nsamples] = fabs(fitdata.Y[i] - Yf) / fitdata.dY[i];
    225         // fprintf (stderr, "%f : %f %f : %f %f : %f %f : %f %f %f\n", T[i], Xf, Yf, fitdata.X[i], fitdata.Y[i], fitdata.dX[i], fitdata.dY[i], fitdata.t[i], fitdata.pX[i], fitdata.pY[i]);
    226       }
    227       Nsamples ++;
    228     }
    229 
    230     double *dPsig;
    231     ALLOCATE (dPsig, double, fitdata.Npts);
    232    
    233     for (i = 0; i < fitdata.Npts; i++) {
    234       dsort (dXsig[i], Nsamples);
    235       dsort (dYsig[i], Nsamples);
    236 
    237       // choose the median values
    238       double dXsigMedian, dYsigMedian;
    239       if (Nsamples % 2) {
    240         int Ncenter = Nsamples / 2;
    241         dXsigMedian = dXsig[i][Ncenter];
    242         dYsigMedian = dYsig[i][Ncenter];
    243       } else {
    244         int Ncenter = Nsamples / 2 - 1;
    245         dXsigMedian = 0.5*(dXsig[i][Ncenter] + dXsig[i][Ncenter + 1]);
    246         dYsigMedian = 0.5*(dYsig[i][Ncenter] + dYsig[i][Ncenter + 1]);
    247       }
    248       // XXX replace with hypotenuse?
    249       dPsig[i] = 0.5*(dXsigMedian + dYsigMedian);
    250       // fprintf (stderr, "%d %10.6f %10.6f %10.6f  %f %f : %f\n", i, R[i], D[i], T[i], dXsig[i][Ncenter], dYsig[i][Ncenter], dPsig[i]);
    251     }
    252 
    253     int Nout = 0;
    254     for (i = 0; i < fitdata.Npts; i++) {
    255       if (dPsig[i] > dPsigMax) {
    256         fprintf (stderr, "clip %d: %f : %f\n", i, fitdata.t[i], dPsig[i]);
    257         continue;
    258       }
    259       sample.X [Nout] = fitdata.X [i];
    260       sample.Y [Nout] = fitdata.Y [i];
    261       sample.dX[Nout] = fitdata.dX[i];
    262       sample.dY[Nout] = fitdata.dY[i];
    263       sample.t [Nout] = fitdata.t [i];
    264       sample.pX[Nout] = fitdata.pX[i];
    265       sample.pY[Nout] = fitdata.pY[i];
    266       Nout ++;
    267     }
    268     sample.Npts = Nout;
    269     fprintf (stderr, "keep %d of %d\n", sample.Npts, fitdata.Npts);
    270 
    271     if (dPvec) {
    272         ResetVector (dPvec, OPIHI_FLT, Ntotal);
    273         for (i = 0; i < Ntotal; i++) {
    274           dPvec->elements.Flt[i] = NAN;
    275         }
    276         for (i = 0; i < fitdata.Npts; i++) {
    277           int n = index[i];
    278           dPvec->elements.Flt[n] = dPsig[i];
    279         }
    280         dPvec->Nelements = Ntotal;
    281     }
    282     free (dPsig);
    283 
    284     fitdata = sample;
    285   }
    286 
    287   for (i = 0; VERBOSE && (i < fitdata.Npts); i++) {
    288     int n = index[i];
    289     fprintf (stderr, "%f %f : %f %d : %f %f %f\n", R[n], D[n], T[n], mask[n], fitdata.t[i], fitdata.X[i], fitdata.Y[i]);
     152    int clipRetry = TRUE;
     153    for (i = 0; clipRetry && (i < 3); i++) {
     154      clipRetry = !PlxOutlierClip (&fitdata, mask, Noutlier, dPsigMax, dPvec, VERBOSE);
     155
     156      // using the new mask values, reset fitdata
     157      PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
     158      PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
     159      if (VERBOSE) fprintf (stderr, "keep %d of %d\n", fitdata.Npts, Ntotal);
     160    }
     161  }
     162
     163  for (i = 0; (VERBOSE == 2) && (i < fitdata.Npts); i++) {
     164    int n = fitdata.index[i];
     165    int maskValue = mask ? mask[n] : 0;
     166    fprintf (stderr, "%f %f : %f %d : %f %f %f\n", R[n], D[n], T[n], maskValue, fitdata.t[i], fitdata.X[i], fitdata.Y[i]);
    290167  }
    291168
    … …  
    298175  }
    299176
     177  if (Nresample){
     178    PlxFitData sample;
     179    PlxFitDataAlloc (&sample, fitdata.Npts);
     180
     181    PlxFit *testfit = NULL;
     182    ALLOCATE (testfit, PlxFit, Nresample);
     183
     184    int Ngood = 0;
     185    for (i = 0; i < Nresample; i++) {
     186      PlxBootstrapResample (&fitdata, &sample);
     187     
     188      if (i % 100000 == 99999) fprintf (stderr, ".");
     189
     190      // fit the sample
     191      testfit[Ngood].getChisq = FALSE;
     192      if (!FitPMandPar (&testfit[Ngood],
     193                        sample.X, sample.dX,
     194                        sample.Y, sample.dY, sample.t,
     195                        sample.pX, sample.pY, sample.Npts, VERBOSE)) continue;
     196      Ngood ++;
     197    }
     198
     199    Vector *pvec, *uRvec, *uDvec, *Rvec, *Dvec;
     200
     201    // save the Nresample histograms
     202    if ((pvec  = SelectVector ("plxVector", ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "plxVector");
     203    if ((uRvec = SelectVector ("uRVector",  ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "uDVector");
     204    if ((uDvec = SelectVector ("uDVector",  ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "uRVector");
     205    if ((Rvec  = SelectVector ("RoVector",  ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "RoVector");
     206    if ((Dvec  = SelectVector ("DoVector",  ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "DoVector");
     207   
     208    ResetVector ( pvec, OPIHI_FLT, Ngood);
     209    ResetVector (uRvec, OPIHI_FLT, Ngood);
     210    ResetVector (uDvec, OPIHI_FLT, Ngood);
     211    ResetVector ( Rvec, OPIHI_FLT, Ngood);
     212    ResetVector ( Dvec, OPIHI_FLT, Ngood);
     213   
     214    for (i = 0; i < Ngood; i++) {
     215      pvec->elements.Flt[i]  = testfit[i].p;
     216      uRvec->elements.Flt[i] = testfit[i].uR;
     217      uDvec->elements.Flt[i] = testfit[i].uD;
     218      Rvec->elements.Flt[i]  = testfit[i].Ro;
     219      Dvec->elements.Flt[i]  = testfit[i].Do;
     220    }
     221
     222    // now calculate median and sigma for each vector
     223    VectorRobustStats (pvec,  &fit.p,  &fit.dp);
     224    VectorRobustStats (uRvec, &fit.uR, &fit.duR);
     225    VectorRobustStats (uDvec, &fit.uD, &fit.duD);
     226    VectorRobustStats (Rvec,  &fit.Ro, &fit.dRo);
     227    VectorRobustStats (Dvec,  &fit.Do, &fit.dDo);
     228  }
     229
    300230  // fprintf (stderr, "%f +/- %f | %f %f\n", fit.p, fit.dp, fit.uR, fit.uD);
    301231
    … …  
    326256  set_variable ("uR",   fit.uR);
    327257  set_variable ("uD",   fit.uD);
    328   set_variable ("duR",   fit.duR);
    329   set_variable ("duD",   fit.duD);
     258  set_variable ("duR",  fit.duR);
     259  set_variable ("duD",  fit.duD);
    330260  set_variable ("plx",  fit.p);
    331261  set_variable ("dplx", fit.dp);
    332262 
    333263  set_variable ("Tmean",  Tmean);
    334   set_variable ("Trange", Trange);
    335   set_variable ("Prange", parRange);
    336264
    337265  set_variable ("chisq", fit.chisq);
    … …  
    525453  ALLOCATE (data->pX, double, N);
    526454  ALLOCATE (data->pY, double, N);
     455  ALLOCATE (data->index, int, N);
    527456  return TRUE;
    528457}
    … …  
    536465  FREE (data->pX);
    537466  FREE (data->pY);
     467  FREE (data->index);
    538468}
    539469
    … …  
    555485}
    556486
    557 # if (0)
    558 int PlxSetMeanEpoch () {
     487int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal) {
     488
     489  int i;
    559490
    560491  // find mean values to remove
    561492  double Nmean = 0;
    562   double Tmean = 0;
    563   double Rmean = 0;
    564   double Dmean = 0;
     493  *Tmean = 0;
     494  *Rmean = 0;
     495  *Dmean = 0;
    565496  double Tmin = +1000000;
    566497  double Tmax = -1000000;
    567498  for (i = 0; i < Ntotal; i++) {
    568499    if (mask && !mask[i]) continue;
    569     Rmean += R[i];
    570     Dmean += D[i];
    571     Tmean += T[i];
     500    *Rmean += R[i];
     501    *Dmean += D[i];
     502    *Tmean += T[i];
    572503    Tmin = MIN(Tmin, T[i]);
    573504    Tmax = MAX(Tmax, T[i]);
    574505    Nmean += 1.0;
    575506  }
    576   Rmean /= Nmean;
    577   Dmean /= Nmean;
    578   Tmean /= Nmean;
    579 
    580   float Trange = Tmax - Tmin;
    581 
    582 }
    583 # endif
     507  *Rmean /= Nmean;
     508  *Dmean /= Nmean;
     509  *Tmean /= Nmean;
     510 
     511  double Trange = Tmax - Tmin;
     512
     513  // fprintf (stderr, "R,D : %f,%f, T: %f, Trange: %f, Tmin: %f, Tmax: %f\n", *Rmean, *Dmean, *Tmean, Trange, Tmin, Tmax);
     514
     515  set_variable ("Trange", Trange);
     516  return TRUE;
     517}
     518
     519// generate the fit values (projected X,Y; parallax factors;
     520int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean) {
     521
     522  int i;
     523
     524  float pXmin = +2.0;
     525  float pXmax = -2.0;
     526  float pYmin = +2.0;
     527  float pYmax = -2.0;
     528
     529  int Nsubset = 0;
     530  for (i = 0; i < Ntotal; i++) {
     531    if (mask && !mask[i]) continue;
     532    RD_to_XY (&fitdata->X[Nsubset], &fitdata->Y[Nsubset], R[i], D[i], coords);
     533    fitdata->dX[Nsubset] = dR[i];
     534    fitdata->dY[Nsubset] = dD[i];
     535    fitdata->t[Nsubset] = (T[i] - Tmean) / 365.25;
     536    ParFactor (&fitdata->pX[Nsubset], &fitdata->pY[Nsubset], R[i], D[i], T[i]);
     537    pXmin = MIN (pXmin, fitdata->pX[Nsubset]);
     538    pXmax = MAX (pXmax, fitdata->pX[Nsubset]);
     539    pYmin = MIN (pYmin, fitdata->pY[Nsubset]);
     540    pYmax = MAX (pYmax, fitdata->pY[Nsubset]);
     541    fitdata->index[Nsubset] = i;
     542    Nsubset++;
     543  }
     544  fitdata->Npts = Nsubset;
     545  float dXRange = pXmax - pXmin;
     546  float dYRange = pYmax - pYmin;
     547  float parRange = hypot (dXRange, dYRange);
     548
     549  set_variable ("Prange", parRange);
     550  // fprintf (stderr, "par factor range: %f\n", parRange);
     551
     552  return TRUE;
     553}
     554
     555/* Outlier clipping based on bootstrap-resampling tests of the plx path
     556 * generate Noutlier resampled datasets
     557 * fit the Noutlier plx paths
     558 * determine and save the distribution of dXsig and dYsig for each point
     559 * sort the resulting distributions and find dPsig (median point) for each measurement
     560 * find the 90% point of dPsig : if > dPsigMax, only clip the 10% most deviant points
     561 * set the dPvec values if desired
     562 * -- mask is modified, dPvec values are set
     563 * -- fitdata is unchanged
     564 */
     565
     566# define MAX_REJECT 0.1
     567
     568int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE) {
     569
     570  int i, n;
     571
     572  PlxFit testfit;
     573  testfit.getChisq = FALSE;
     574
     575  PlxFitData sample;
     576  PlxFitDataAlloc (&sample, fitdata->Npts);
     577
     578  double **dXsig, **dYsig;
     579  ALLOCATE (dXsig, double *, fitdata->Npts);
     580  ALLOCATE (dYsig, double *, fitdata->Npts);
     581  for (i = 0; i < fitdata->Npts; i++) {
     582    ALLOCATE (dXsig[i], double, Noutlier);
     583    ALLOCATE (dYsig[i], double, Noutlier);
     584  }
     585
     586  int Nsamples = 0;
     587  for (n = 0; n < Noutlier; n++) {
     588    // bootstrap resample (fitdata -> sample)
     589    PlxBootstrapResample (fitdata, &sample);
     590     
     591    if (n % 100000 == 99999) fprintf (stderr, ".");
     592
     593    // fit the sample
     594    if (!FitPMandPar (&testfit,
     595                      sample.X, sample.dX,
     596                      sample.Y, sample.dY, sample.t,
     597                      sample.pX, sample.pY, sample.Npts, VERBOSE)) continue;
     598
     599    // fprintf (stderr, "%f +/- %f | %f %f\n", testfit.p, testfit.dp, testfit.uR, testfit.uD);
     600
     601    // find the distances to the path
     602    for (i = 0; i < fitdata->Npts; i++) {
     603      double Xf = testfit.Ro + testfit.uR*fitdata->t[i] + testfit.p*fitdata->pX[i];
     604      double Yf = testfit.Do + testfit.uD*fitdata->t[i] + testfit.p*fitdata->pY[i];
     605      dXsig[i][Nsamples] = fabs(fitdata->X[i] - Xf) / fitdata->dX[i];
     606      dYsig[i][Nsamples] = fabs(fitdata->Y[i] - Yf) / fitdata->dY[i];
     607      // fprintf (stderr, "%f : %f %f : %f %f : %f %f : %f %f %f\n", T[i], Xf, Yf, fitdata->X[i], fitdata->Y[i], fitdata->dX[i], fitdata->dY[i], fitdata->t[i], fitdata->pX[i], fitdata->pY[i]);
     608    }
     609    Nsamples ++;
     610  }
     611
     612  double *dPsig;
     613  ALLOCATE (dPsig, double, fitdata->Npts);
     614   
     615  for (i = 0; i < fitdata->Npts; i++) {
     616    dsort (dXsig[i], Nsamples);
     617    dsort (dYsig[i], Nsamples);
     618
     619    // choose the median values
     620    double dXsigMedian, dYsigMedian;
     621    if (Nsamples % 2) {
     622      int Ncenter = Nsamples / 2;
     623      dXsigMedian = dXsig[i][Ncenter];
     624      dYsigMedian = dYsig[i][Ncenter];
     625    } else {
     626      int Ncenter = Nsamples / 2 - 1;
     627      dXsigMedian = 0.5*(dXsig[i][Ncenter] + dXsig[i][Ncenter + 1]);
     628      dYsigMedian = 0.5*(dYsig[i][Ncenter] + dYsig[i][Ncenter + 1]);
     629    }
     630    // XXX replace with hypotenuse?
     631    dPsig[i] = 0.5*(dXsigMedian + dYsigMedian);
     632    // fprintf (stderr, "%d %10.6f %10.6f %10.6f  %f %f : %f\n", i, R[i], D[i], T[i], dXsig[i][Ncenter], dYsig[i][Ncenter], dPsig[i]);
     633  }
     634
     635  // make a copy of dPsig[] and check if > 10% are > dPsigMax
     636  double *dPsigSort;
     637  ALLOCATE (dPsigSort, double, fitdata->Npts);
     638  for (i = 0; i < fitdata->Npts; i++) {
     639    dPsigSort[i] = dPsig[i];
     640  }
     641  dsort (dPsigSort, fitdata->Npts);
     642  int Nmax = (1.0 - MAX_REJECT)*fitdata->Npts;
     643
     644  int completeClip = TRUE;
     645  if (dPsigSort[Nmax] > dPsigMax) {
     646    if (VERBOSE) fprintf (stderr, "too many outliers: %f at 90\n", dPsigSort[Nmax]);
     647    dPsigMax = dPsigSort[Nmax];
     648    completeClip = FALSE;
     649  }
     650
     651  for (i = 0; i < fitdata->Npts; i++) {
     652    if (dPsig[i] < dPsigMax) continue;
     653    int n = fitdata->index[i];
     654    // fprintf (stderr, "clip %d: %f : %f\n", i, fitdata->t[i], dPsig[i]);
     655    mask[n] = 0; // mask these points
     656  }
     657
     658  // only set dPvec if we have completed the clipping?
     659  if (dPvec) {
     660    for (i = 0; i < dPvec->Nelements; i++) {
     661      dPvec->elements.Flt[i] = NAN;
     662    }
     663    for (i = 0; i < fitdata->Npts; i++) {
     664      int n = fitdata->index[i];
     665      dPvec->elements.Flt[n] = dPsig[i];
     666    }
     667  }
     668
     669  free (dPsig);
     670  free (dPsigSort);
     671 
     672  for (i = 0; i < fitdata->Npts; i++) {
     673    free (dXsig[i]);
     674    free (dYsig[i]);
     675  }
     676  free (dXsig);
     677  free (dYsig);
     678
     679  return completeClip;
     680}
     681
     682int VectorRobustStats (Vector *vector, double *median, double *sigma) {
     683
     684  // warn if vector->Nelements > 1000? 10000?)
     685  // warn if vector is not float
     686
     687  // we need to copy the vector to avoid changing the sort order
     688  double *values = NULL;
     689  ALLOCATE (values, double, vector->Nelements);
     690
     691  int i;
     692  int Npts = 0;
     693  for (i = 0; i < vector->Nelements; i++) {
     694    if (!isfinite(vector->elements.Flt[i])) continue;
     695    values[Npts] = vector->elements.Flt[i];
     696    Npts++;
     697  }
     698
     699  dsort (values, Npts);
     700
     701  if (Npts % 2) {
     702    int Ncenter = Npts / 2;
     703    *median = values[Ncenter];
     704  } else {
     705    int Ncenter = Npts / 2 - 1;
     706    *median = 0.5*(values[Ncenter] + values[Ncenter + 1]);
     707  }
     708
     709  double Slo = VectorFractionInterpolate (values, 0.158655, Npts);
     710  double Shi = VectorFractionInterpolate (values, 0.841345, Npts);
     711
     712  *sigma = (Shi - Slo) / 2.0;
     713
     714  return TRUE;
     715}
     716
     717double VectorFractionInterpolate (double *values, float fraction, int Npts) {
     718
     719  float F = fraction * Npts;
     720  int   N = fraction * Npts;
     721
     722  if (N < 0        ) return NAN;
     723  if (N >= Npts - 2) return NAN;
     724
     725  // interpolate between N,N+1
     726   
     727  double S = (F - N) * (values[N+1] - values[N]) + values[N];
     728  return S;
     729}
Note: See TracChangeset for help on using the changeset viewer.