Index: /branches/eam_branches/ipp-20150625/Ohana/src/opihi/cmd.astro/fitplx.c
===================================================================
--- /branches/eam_branches/ipp-20150625/Ohana/src/opihi/cmd.astro/fitplx.c	(revision 38567)
+++ /branches/eam_branches/ipp-20150625/Ohana/src/opihi/cmd.astro/fitplx.c	(revision 38568)
@@ -7,4 +7,15 @@
 
 typedef struct {
+  double *X;
+  double *Y;
+  double *t;
+  double *pX;
+  double *pY;
+  double *dX;
+  double *dY;
+  int Npts;
+} PlxFitData;
+
+typedef struct {
   double Ro, dRo;
   double Do, dDo;
@@ -17,5 +28,10 @@
   double chisq;
   int Nfit;
+  int getChisq;
 } PlxFit;
+
+int PlxFitDataAlloc (PlxFitData *data, int N);
+void PlxFitDataFree (PlxFitData *data);
+int PlxBootstrapResample (PlxFitData *src, PlxFitData *tgt);
 
 int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE);
@@ -47,7 +63,26 @@
   }
 
+  int Noutlier = 0;
+  float dPsigMax = FLT_MAX;
+  if ((N = get_argument (argc, argv, "-outlier-tests"))) {
+    remove_argument (N, &argc, argv);
+    Noutlier = atoi(argv[N]);
+    remove_argument (N, &argc, argv);
+    dPsigMax = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  Vector *dPvec = NULL;
+  if ((N = get_argument (argc, argv, "-dPsig"))) {
+    if (!Noutlier) { gprint (GP_ERR, "-dPsig requires -outlier-tests to be non-zero\n"); return FALSE; }
+    remove_argument (N, &argc, argv);
+    if (!(dPvec = SelectVector (argv[N], ANYVECTOR, TRUE))) return FALSE;
+    remove_argument (N, &argc, argv);
+  }
+
   if (argc != 6) {
-    gprint (GP_ERR, "USAGE: fitplx (ra) (dR) (dec) (dD) (mjd) [-mask mask]\n");
-    // what about the errors?
+    gprint (GP_ERR, "USAGE: fitplx (ra) (dR) (dec) (dD) (mjd) [-mask mask] [-v] [-vv]\n");
+    gprint (GP_ERR, "  -outlier-tests Nsamples dPsigMax : run Nsample bootstrap-resamples to define the path deviations and reject based on dPsigMax\n");
+    gprint (GP_ERR, "  -dPsig vec : save path deviations in vec\n");
     return (FALSE);
   }
@@ -77,8 +112,8 @@
   }
 
-  N = tvec->Nelements; // XXX check other lengths
+  int Ntotal = tvec->Nelements; // XXX check other lengths
 
   // find mean values to remove
-  double Npts = 0;
+  double Nmean = 0;
   double Tmean = 0;
   double Rmean = 0;
@@ -86,5 +121,5 @@
   double Tmin = +1000000;
   double Tmax = -1000000;
-  for (i = 0; i < N; i++) {
+  for (i = 0; i < Ntotal; i++) {
     if (mask && !mask[i]) continue;
     Rmean += R[i];
@@ -93,9 +128,13 @@
     Tmin = MIN(Tmin, T[i]);
     Tmax = MAX(Tmax, T[i]);
-    Npts += 1.0;
-  }
-  Rmean /= Npts;
-  Dmean /= Npts;
-  Tmean /= Npts;
+    Nmean += 1.0;
+  }
+  Rmean /= Nmean;
+  Dmean /= Nmean;
+  Tmean /= Nmean;
+
+  // Ntotal : all points supplied by user
+  // Nsubset : unmasked points
+  // Nmean : unmasked points (a double only used above)
 
   float Trange = Tmax - Tmin;
@@ -109,13 +148,4 @@
   coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
 
-  double *X, *Y, *t, *pX, *pY, *dX, *dY;
-  ALLOCATE (X, double, N);
-  ALLOCATE (Y, double, N);
-  ALLOCATE (dX, double, N);
-  ALLOCATE (dY, double, N);
-  ALLOCATE (t, double, N);
-  ALLOCATE (pX, double, N);
-  ALLOCATE (pY, double, N);
-
   float pXmin = +2.0;
   float pXmax = -2.0;
@@ -123,18 +153,30 @@
   float pYmax = -2.0;
 
-  int n = 0;
-  for (i = 0; i < N; i++) {
+  PlxFit fit; memset (&fit, 0, sizeof(PlxFit));
+  PlxFitData fitdata;
+  PlxFitDataAlloc (&fitdata, Ntotal);
+
+  // generate the fit values (projected X,Y; parallax factors; 
+
+  // save an index so we can supplied dPsig for the unmasked points
+  int *index;  
+  ALLOCATE (index, int, Ntotal);
+  
+  int Nsubset = 0;
+  for (i = 0; i < Ntotal; i++) {
     if (mask && !mask[i]) continue;
-    RD_to_XY (&X[n], &Y[n], R[i], D[i], &coords);
-    dX[n] = dR[i];
-    dY[n] = dD[i];
-    t[n] = (T[i] - Tmean) / 365.25;
-    ParFactor (&pX[n], &pY[n], R[i], D[i], T[i]);
-    pXmin = MIN (pXmin, pX[n]);
-    pXmax = MAX (pXmax, pX[n]);
-    pYmin = MIN (pYmin, pY[n]);
-    pYmax = MAX (pYmax, pY[n]);
-    n++;
-  }
+    RD_to_XY (&fitdata.X[Nsubset], &fitdata.Y[Nsubset], R[i], D[i], &coords);
+    fitdata.dX[Nsubset] = dR[i];
+    fitdata.dY[Nsubset] = dD[i];
+    fitdata.t[Nsubset] = (T[i] - Tmean) / 365.25;
+    ParFactor (&fitdata.pX[Nsubset], &fitdata.pY[Nsubset], R[i], D[i], T[i]);
+    pXmin = MIN (pXmin, fitdata.pX[Nsubset]);
+    pXmax = MAX (pXmax, fitdata.pX[Nsubset]);
+    pYmin = MIN (pYmin, fitdata.pY[Nsubset]);
+    pYmax = MAX (pYmax, fitdata.pY[Nsubset]);
+    index[Nsubset] = i;
+    Nsubset++;
+  }
+  fitdata.Npts = Nsubset;
   float dXRange = pXmax - pXmin;
   float dYRange = pYmax - pYmin;
@@ -143,8 +185,128 @@
   // fprintf (stderr, "par factor range: %f\n", parRange);
 
-  PlxFit fit;
-  if (!FitPMandPar (&fit, X, dX, Y, dY, t, pX, pY, n, VERBOSE)) {
+  // determine dPsig for detections based on Noutlier attempts
+  if (Noutlier) {
+    PlxFit testfit;
+    PlxFitData sample;
+    PlxFitDataAlloc (&sample, fitdata.Npts);
+
+    double **dXsig, **dYsig;
+    ALLOCATE (dXsig, double *, fitdata.Npts);
+    ALLOCATE (dYsig, double *, fitdata.Npts);
+    for (i = 0; i < fitdata.Npts; i++) {
+      ALLOCATE (dXsig[i], double, Noutlier);
+      ALLOCATE (dYsig[i], double, Noutlier);
+    }
+
+    testfit.getChisq = FALSE;
+
+    int n;
+    int Nsamples = 0;
+    for (n = 0; n < Noutlier; n++) {
+      // bootstrap resample (fitdata -> sample)
+      PlxBootstrapResample (&fitdata, &sample);
+      
+      if (n % 1000 == 999) fprintf (stderr, ".");
+
+      // fit the sample
+      if (!FitPMandPar (&testfit, 
+			sample.X, sample.dX, 
+			sample.Y, sample.dY, sample.t, 
+			sample.pX, sample.pY, sample.Npts, VERBOSE)) continue;
+
+      // fprintf (stderr, "%f +/- %f | %f %f\n", testfit.p, testfit.dp, testfit.uR, testfit.uD);
+
+      // find the distances to the path
+      for (i = 0; i < fitdata.Npts; i++) {
+	double Xf = testfit.Ro + testfit.uR*fitdata.t[i] + testfit.p*fitdata.pX[i];
+	double Yf = testfit.Do + testfit.uD*fitdata.t[i] + testfit.p*fitdata.pY[i];
+	dXsig[i][Nsamples] = fabs(fitdata.X[i] - Xf) / fitdata.dX[i];
+	dYsig[i][Nsamples] = fabs(fitdata.Y[i] - Yf) / fitdata.dY[i];
+        // 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]);
+      }
+      Nsamples ++;
+    }
+
+    double *dPsig;
+    ALLOCATE (dPsig, double, fitdata.Npts);
+    
+    for (i = 0; i < fitdata.Npts; i++) {
+      dsort (dXsig[i], Nsamples);
+      dsort (dYsig[i], Nsamples);
+
+      // choose the median values
+      double dXsigMedian, dYsigMedian;
+      if (Nsamples % 2) {
+	int Ncenter = Nsamples / 2;
+	dXsigMedian = dXsig[i][Ncenter];
+	dYsigMedian = dYsig[i][Ncenter];
+      } else {
+	int Ncenter = Nsamples / 2 - 1;
+	dXsigMedian = 0.5*(dXsig[i][Ncenter] + dXsig[i][Ncenter + 1]);
+	dYsigMedian = 0.5*(dYsig[i][Ncenter] + dYsig[i][Ncenter + 1]);
+      }
+      // XXX replace with hypotenuse?
+      dPsig[i] = 0.5*(dXsigMedian + dYsigMedian);
+      // 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]);
+    }
+
+    int Nout = 0;
+    for (i = 0; i < fitdata.Npts; i++) {
+      if (dPsig[i] > dPsigMax) {
+	fprintf (stderr, "clip %d: %f : %f\n", i, fitdata.t[i], dPsig[i]);
+	continue;
+      }
+      sample.X [Nout] = fitdata.X [i];
+      sample.Y [Nout] = fitdata.Y [i];
+      sample.dX[Nout] = fitdata.dX[i];
+      sample.dY[Nout] = fitdata.dY[i];
+      sample.t [Nout] = fitdata.t [i];
+      sample.pX[Nout] = fitdata.pX[i];
+      sample.pY[Nout] = fitdata.pY[i];
+      Nout ++;
+    }
+    sample.Npts = Nout;
+    fprintf (stderr, "keep %d of %d\n", sample.Npts, fitdata.Npts);
+
+    if (dPvec) {
+        ResetVector (dPvec, OPIHI_FLT, Ntotal);
+	for (i = 0; i < Ntotal; i++) {
+ 	  dPvec->elements.Flt[i] = NAN;
+	}
+	for (i = 0; i < fitdata.Npts; i++) {
+	  int n = index[i];
+ 	  dPvec->elements.Flt[n] = dPsig[i];
+	}
+	dPvec->Nelements = Ntotal;
+    }
+    free (dPsig);
+
+    fitdata = sample;
+  }
+
+  for (i = 0; VERBOSE && (i < fitdata.Npts); i++) {
+    int n = index[i];
+    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]);
+  }
+
+  fit.getChisq = TRUE;
+  if (!FitPMandPar (&fit, 
+		    fitdata.X, fitdata.dX, 
+		    fitdata.Y, fitdata.dY, 
+		    fitdata.t, fitdata.pX, fitdata.pY, fitdata.Npts, VERBOSE)) {
     return FALSE;
   }
+
+  // fprintf (stderr, "%f +/- %f | %f %f\n", fit.p, fit.dp, fit.uR, fit.uD);
+
+/*
+  FILE *f = fopen ("test.pf.dat", "w");
+  for (i = 0; i < Ntotal; i++) {
+    double Xf = fit.Ro + fit.uR*fitdata.t[i] + fit.p*fitdata.pX[i];
+    double Yf = fit.Do + fit.uD*fitdata.t[i] + fit.p*fitdata.pY[i];
+    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]);
+  }
+  fclose (f);
+*/
 
   // fprintf (stderr, "Roff, Doff: %f, %f; dRo, dDo: %f, %f\n", fit.Ro, fit.Do, fit.dRo, fit.dDo);
@@ -152,5 +314,5 @@
   XY_to_RD (&Rmean, &Dmean, fit.Ro, fit.Do, &coords);
   if (VERBOSE) {
-    fprintf (stderr, "Ro, Do: %f, %f +/- %f, %f\n", Rmean, Dmean, fit.dRo, fit.dDo);
+    fprintf (stderr, "Ro, Do: %f, %f +/- %f, %f (%f, %f)\n", Rmean, Dmean, fit.dRo, fit.dDo, fit.Ro, fit.Do);
     fprintf (stderr, "uR, uD: %f, %f; duR, duD: %f, %f\n", fit.uR, fit.uD, fit.duR, fit.duD);
     fprintf (stderr, "par: %f +/- %f\n", fit.p, fit.dp);
@@ -293,19 +455,21 @@
   fit[0].dp  = sqrt(A[4][4]);
   
-  // add up the chi square for the fit
-  chisq = 0.0;
-  for (i = 0; i < Npts; i++) {
-    Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].p*pR[i];
-    Yf = fit[0].Do + fit[0].uD*T[i] + fit[0].p*pD[i];
-    wx = (fabs(dX[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dX[i]);
-    wy = (fabs(dY[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dY[i]);
-    chisq += SQ(X[i] - Xf) * wx;
-    chisq += SQ(Y[i] - Yf) * wy;
-    // 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);
-  }
+  // (optionally) add up the chi square for the fit
+  if (fit->getChisq) {
+    chisq = 0.0;
+    for (i = 0; i < Npts; i++) {
+      Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].p*pR[i];
+      Yf = fit[0].Do + fit[0].uD*T[i] + fit[0].p*pD[i];
+      wx = (fabs(dX[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dX[i]);
+      wy = (fabs(dY[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dY[i]);
+      chisq += SQ(X[i] - Xf) * wx;
+      chisq += SQ(Y[i] - Yf) * wy;
+      // 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);
+    }
+    // the reduced chisq is divided by (Ndof = 2*Npts - 5)
+    fit[0].chisq = chisq / (2.0*Npts - 5.0);
+  }
+  
   fit[0].Nfit = Npts;
-
-  // the reduced chisq is divided by (Ndof = 2*Npts - 5)
-  fit[0].chisq = chisq / (2.0*Npts - 5.0);
   return (TRUE);
 }
@@ -349,2 +513,71 @@
   return TRUE;
 }
+
+// allocate arrays but not the container
+int PlxFitDataAlloc (PlxFitData *data, int N) {
+
+  data->Npts = N;
+  ALLOCATE (data->X, double, N);
+  ALLOCATE (data->Y, double, N);
+  ALLOCATE (data->dX, double, N);
+  ALLOCATE (data->dY, double, N);
+  ALLOCATE (data->t, double, N);
+  ALLOCATE (data->pX, double, N);
+  ALLOCATE (data->pY, double, N);
+  return TRUE;
+}
+
+void PlxFitDataFree (PlxFitData *data) {
+  FREE (data->X);
+  FREE (data->Y);
+  FREE (data->dX);
+  FREE (data->dY);
+  FREE (data->t);
+  FREE (data->pX);
+  FREE (data->pY);
+}
+
+int PlxBootstrapResample (PlxFitData *src, PlxFitData *tgt) {
+  int i;
+  tgt->Npts = src->Npts;
+  for (i = 0; i < src->Npts; i++) {
+    int N = tgt->Npts * drand48();
+    // int N = i;
+    tgt->X [i] = src->X [N];
+    tgt->Y [i] = src->Y [N];
+    tgt->dX[i] = src->dX[N];
+    tgt->dY[i] = src->dY[N];
+    tgt->t [i] = src->t [N];
+    tgt->pX[i] = src->pX[N];
+    tgt->pY[i] = src->pY[N];
+  }
+  return TRUE;
+}
+
+# if (0)
+int PlxSetMeanEpoch () {
+
+  // find mean values to remove
+  double Nmean = 0;
+  double Tmean = 0;
+  double Rmean = 0;
+  double Dmean = 0;
+  double Tmin = +1000000;
+  double Tmax = -1000000;
+  for (i = 0; i < Ntotal; i++) {
+    if (mask && !mask[i]) continue;
+    Rmean += R[i];
+    Dmean += D[i];
+    Tmean += T[i];
+    Tmin = MIN(Tmin, T[i]);
+    Tmax = MAX(Tmax, T[i]);
+    Nmean += 1.0;
+  }
+  Rmean /= Nmean;
+  Dmean /= Nmean;
+  Tmean /= Nmean;
+
+  float Trange = Tmax - Tmin;
+
+}
+# endif