IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38568


Ignore:
Timestamp:
Jul 4, 2015, 7:42:08 AM (11 years ago)
Author:
eugene
Message:

fitplx does outlier clipping ok

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r37807 r38568  
    77
    88typedef struct {
     9  double *X;
     10  double *Y;
     11  double *t;
     12  double *pX;
     13  double *pY;
     14  double *dX;
     15  double *dY;
     16  int Npts;
     17} PlxFitData;
     18
     19typedef struct {
    920  double Ro, dRo;
    1021  double Do, dDo;
    … …  
    1728  double chisq;
    1829  int Nfit;
     30  int getChisq;
    1931} PlxFit;
     32
     33int PlxFitDataAlloc (PlxFitData *data, int N);
     34void PlxFitDataFree (PlxFitData *data);
     35int PlxBootstrapResample (PlxFitData *src, PlxFitData *tgt);
    2036
    2137int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE);
    … …  
    4763  }
    4864
     65  int Noutlier = 0;
     66  float dPsigMax = FLT_MAX;
     67  if ((N = get_argument (argc, argv, "-outlier-tests"))) {
     68    remove_argument (N, &argc, argv);
     69    Noutlier = atoi(argv[N]);
     70    remove_argument (N, &argc, argv);
     71    dPsigMax = atof(argv[N]);
     72    remove_argument (N, &argc, argv);
     73  }
     74
     75  Vector *dPvec = NULL;
     76  if ((N = get_argument (argc, argv, "-dPsig"))) {
     77    if (!Noutlier) { gprint (GP_ERR, "-dPsig requires -outlier-tests to be non-zero\n"); return FALSE; }
     78    remove_argument (N, &argc, argv);
     79    if (!(dPvec = SelectVector (argv[N], ANYVECTOR, TRUE))) return FALSE;
     80    remove_argument (N, &argc, argv);
     81  }
     82
    4983  if (argc != 6) {
    50     gprint (GP_ERR, "USAGE: fitplx (ra) (dR) (dec) (dD) (mjd) [-mask mask]\n");
    51     // what about the errors?
     84    gprint (GP_ERR, "USAGE: fitplx (ra) (dR) (dec) (dD) (mjd) [-mask mask] [-v] [-vv]\n");
     85    gprint (GP_ERR, "  -outlier-tests Nsamples dPsigMax : run Nsample bootstrap-resamples to define the path deviations and reject based on dPsigMax\n");
     86    gprint (GP_ERR, "  -dPsig vec : save path deviations in vec\n");
    5287    return (FALSE);
    5388  }
    … …  
    77112  }
    78113
    79   N = tvec->Nelements; // XXX check other lengths
     114  int Ntotal = tvec->Nelements; // XXX check other lengths
    80115
    81116  // find mean values to remove
    82   double Npts = 0;
     117  double Nmean = 0;
    83118  double Tmean = 0;
    84119  double Rmean = 0;
    … …  
    86121  double Tmin = +1000000;
    87122  double Tmax = -1000000;
    88   for (i = 0; i < N; i++) {
     123  for (i = 0; i < Ntotal; i++) {
    89124    if (mask && !mask[i]) continue;
    90125    Rmean += R[i];
    … …  
    93128    Tmin = MIN(Tmin, T[i]);
    94129    Tmax = MAX(Tmax, T[i]);
    95     Npts += 1.0;
    96   }
    97   Rmean /= Npts;
    98   Dmean /= Npts;
    99   Tmean /= Npts;
     130    Nmean += 1.0;
     131  }
     132  Rmean /= Nmean;
     133  Dmean /= Nmean;
     134  Tmean /= Nmean;
     135
     136  // Ntotal : all points supplied by user
     137  // Nsubset : unmasked points
     138  // Nmean : unmasked points (a double only used above)
    100139
    101140  float Trange = Tmax - Tmin;
    … …  
    109148  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
    110149
    111   double *X, *Y, *t, *pX, *pY, *dX, *dY;
    112   ALLOCATE (X, double, N);
    113   ALLOCATE (Y, double, N);
    114   ALLOCATE (dX, double, N);
    115   ALLOCATE (dY, double, N);
    116   ALLOCATE (t, double, N);
    117   ALLOCATE (pX, double, N);
    118   ALLOCATE (pY, double, N);
    119 
    120150  float pXmin = +2.0;
    121151  float pXmax = -2.0;
    … …  
    123153  float pYmax = -2.0;
    124154
    125   int n = 0;
    126   for (i = 0; i < N; i++) {
     155  PlxFit fit; memset (&fit, 0, sizeof(PlxFit));
     156  PlxFitData fitdata;
     157  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++) {
    127167    if (mask && !mask[i]) continue;
    128     RD_to_XY (&X[n], &Y[n], R[i], D[i], &coords);
    129     dX[n] = dR[i];
    130     dY[n] = dD[i];
    131     t[n] = (T[i] - Tmean) / 365.25;
    132     ParFactor (&pX[n], &pY[n], R[i], D[i], T[i]);
    133     pXmin = MIN (pXmin, pX[n]);
    134     pXmax = MAX (pXmax, pX[n]);
    135     pYmin = MIN (pYmin, pY[n]);
    136     pYmax = MAX (pYmax, pY[n]);
    137     n++;
    138   }
     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;
    139181  float dXRange = pXmax - pXmin;
    140182  float dYRange = pYmax - pYmin;
    … …  
    143185  // fprintf (stderr, "par factor range: %f\n", parRange);
    144186
    145   PlxFit fit;
    146   if (!FitPMandPar (&fit, X, dX, Y, dY, t, pX, pY, n, VERBOSE)) {
     187  // determine dPsig for detections based on Noutlier attempts
     188  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]);
     290  }
     291
     292  fit.getChisq = TRUE;
     293  if (!FitPMandPar (&fit,
     294                    fitdata.X, fitdata.dX,
     295                    fitdata.Y, fitdata.dY,
     296                    fitdata.t, fitdata.pX, fitdata.pY, fitdata.Npts, VERBOSE)) {
    147297    return FALSE;
    148298  }
     299
     300  // fprintf (stderr, "%f +/- %f | %f %f\n", fit.p, fit.dp, fit.uR, fit.uD);
     301
     302/*
     303  FILE *f = fopen ("test.pf.dat", "w");
     304  for (i = 0; i < Ntotal; i++) {
     305    double Xf = fit.Ro + fit.uR*fitdata.t[i] + fit.p*fitdata.pX[i];
     306    double Yf = fit.Do + fit.uD*fitdata.t[i] + fit.p*fitdata.pY[i];
     307    fprintf (f, "%f : %f %f : %f %f : %f : %f %f : %f %f\n", T[i], R[i], D[i], Xf, Yf, fitdata.t[i], fitdata.X[i], fitdata.Y[i], fitdata.pX[i], fitdata.pY[i]);
     308  }
     309  fclose (f);
     310*/
    149311
    150312  // fprintf (stderr, "Roff, Doff: %f, %f; dRo, dDo: %f, %f\n", fit.Ro, fit.Do, fit.dRo, fit.dDo);
    … …  
    152314  XY_to_RD (&Rmean, &Dmean, fit.Ro, fit.Do, &coords);
    153315  if (VERBOSE) {
    154     fprintf (stderr, "Ro, Do: %f, %f +/- %f, %f\n", Rmean, Dmean, fit.dRo, fit.dDo);
     316    fprintf (stderr, "Ro, Do: %f, %f +/- %f, %f (%f, %f)\n", Rmean, Dmean, fit.dRo, fit.dDo, fit.Ro, fit.Do);
    155317    fprintf (stderr, "uR, uD: %f, %f; duR, duD: %f, %f\n", fit.uR, fit.uD, fit.duR, fit.duD);
    156318    fprintf (stderr, "par: %f +/- %f\n", fit.p, fit.dp);
    … …  
    293455  fit[0].dp  = sqrt(A[4][4]);
    294456 
    295   // add up the chi square for the fit
    296   chisq = 0.0;
    297   for (i = 0; i < Npts; i++) {
    298     Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].p*pR[i];
    299     Yf = fit[0].Do + fit[0].uD*T[i] + fit[0].p*pD[i];
    300     wx = (fabs(dX[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dX[i]);
    301     wy = (fabs(dY[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dY[i]);
    302     chisq += SQ(X[i] - Xf) * wx;
    303     chisq += SQ(Y[i] - Yf) * wy;
    304     // if (VERBOSE) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq);
    305   }
     457  // (optionally) add up the chi square for the fit
     458  if (fit->getChisq) {
     459    chisq = 0.0;
     460    for (i = 0; i < Npts; i++) {
     461      Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].p*pR[i];
     462      Yf = fit[0].Do + fit[0].uD*T[i] + fit[0].p*pD[i];
     463      wx = (fabs(dX[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dX[i]);
     464      wy = (fabs(dY[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dY[i]);
     465      chisq += SQ(X[i] - Xf) * wx;
     466      chisq += SQ(Y[i] - Yf) * wy;
     467      // if (VERBOSE) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq);
     468    }
     469    // the reduced chisq is divided by (Ndof = 2*Npts - 5)
     470    fit[0].chisq = chisq / (2.0*Npts - 5.0);
     471  }
     472 
    306473  fit[0].Nfit = Npts;
    307 
    308   // the reduced chisq is divided by (Ndof = 2*Npts - 5)
    309   fit[0].chisq = chisq / (2.0*Npts - 5.0);
    310474  return (TRUE);
    311475}
    … …  
    349513  return TRUE;
    350514}
     515
     516// allocate arrays but not the container
     517int PlxFitDataAlloc (PlxFitData *data, int N) {
     518
     519  data->Npts = N;
     520  ALLOCATE (data->X, double, N);
     521  ALLOCATE (data->Y, double, N);
     522  ALLOCATE (data->dX, double, N);
     523  ALLOCATE (data->dY, double, N);
     524  ALLOCATE (data->t, double, N);
     525  ALLOCATE (data->pX, double, N);
     526  ALLOCATE (data->pY, double, N);
     527  return TRUE;
     528}
     529
     530void PlxFitDataFree (PlxFitData *data) {
     531  FREE (data->X);
     532  FREE (data->Y);
     533  FREE (data->dX);
     534  FREE (data->dY);
     535  FREE (data->t);
     536  FREE (data->pX);
     537  FREE (data->pY);
     538}
     539
     540int PlxBootstrapResample (PlxFitData *src, PlxFitData *tgt) {
     541  int i;
     542  tgt->Npts = src->Npts;
     543  for (i = 0; i < src->Npts; i++) {
     544    int N = tgt->Npts * drand48();
     545    // int N = i;
     546    tgt->X [i] = src->X [N];
     547    tgt->Y [i] = src->Y [N];
     548    tgt->dX[i] = src->dX[N];
     549    tgt->dY[i] = src->dY[N];
     550    tgt->t [i] = src->t [N];
     551    tgt->pX[i] = src->pX[N];
     552    tgt->pY[i] = src->pY[N];
     553  }
     554  return TRUE;
     555}
     556
     557# if (0)
     558int PlxSetMeanEpoch () {
     559
     560  // find mean values to remove
     561  double Nmean = 0;
     562  double Tmean = 0;
     563  double Rmean = 0;
     564  double Dmean = 0;
     565  double Tmin = +1000000;
     566  double Tmax = -1000000;
     567  for (i = 0; i < Ntotal; i++) {
     568    if (mask && !mask[i]) continue;
     569    Rmean += R[i];
     570    Dmean += D[i];
     571    Tmean += T[i];
     572    Tmin = MIN(Tmin, T[i]);
     573    Tmax = MAX(Tmax, T[i]);
     574    Nmean += 1.0;
     575  }
     576  Rmean /= Nmean;
     577  Dmean /= Nmean;
     578  Tmean /= Nmean;
     579
     580  float Trange = Tmax - Tmin;
     581
     582}
     583# endif
Note: See TracChangeset for help on using the changeset viewer.