Index: trunk/Ohana/src/opihi/cmd.astro/fitplx.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.astro/fitplx.c	(revision 38972)
+++ trunk/Ohana/src/opihi/cmd.astro/fitplx.c	(revision 38986)
@@ -7,4 +7,16 @@
 
 typedef struct {
+  double *X;
+  double *Y;
+  double *t;
+  double *pX;
+  double *pY;
+  double *dX;
+  double *dY;
+  int *index;
+  int Npts;
+} PlxFitData;
+
+typedef struct {
   double Ro, dRo;
   double Do, dDo;
@@ -17,5 +29,17 @@
   double chisq;
   int Nfit;
+  int getChisq;
 } PlxFit;
+
+int VectorRobustStats (Vector *vector, double *median, double *sigma);
+double VectorFractionInterpolate (double *values, float fraction, int Npts);
+
+int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal);
+int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean);
+int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE);
+
+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 +71,33 @@
   }
 
+  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);
+  }
+
+  int Nresample = 0;
+  if ((N = get_argument (argc, argv, "-bootstrap-resample"))) {
+    remove_argument (N, &argc, argv);
+    Nresample = atoi(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,28 +127,11 @@
   }
 
-  N = tvec->Nelements; // XXX check other lengths
-
-  // find mean values to remove
-  double Npts = 0;
-  double Tmean = 0;
-  double Rmean = 0;
-  double Dmean = 0;
-  double Tmin = +1000000;
-  double Tmax = -1000000;
-  for (i = 0; i < N; 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]);
-    Npts += 1.0;
-  }
-  Rmean /= Npts;
-  Dmean /= Npts;
-  Tmean /= Npts;
-
-  float Trange = Tmax - Tmin;
-  // fprintf (stderr, "R,D : %f,%f, T: %f, Trange: %f, Tmin: %f, Tmax: %f\n", Rmean, Dmean, Tmean, Trange, Tmin, Tmax);
+  // Ntotal : all points supplied by user
+  // Nsubset : unmasked points
+  int Ntotal = tvec->Nelements; // XXX check other lengths
+  if (dPvec) ResetVector (dPvec, OPIHI_FLT, Ntotal);
+
+  double Rmean, Dmean, Tmean;
+  PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
 
   /* project coordinates to a plane centered on the object with units of arcsec */
@@ -109,42 +142,101 @@
   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;
-  float pYmin = +2.0;
-  float pYmax = -2.0;
-
-  int n = 0;
-  for (i = 0; i < N; 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++;
-  }
-  float dXRange = pXmax - pXmin;
-  float dYRange = pYmax - pYmin;
-  float parRange = hypot (dXRange, dYRange);
-	
-  // fprintf (stderr, "par factor range: %f\n", parRange);
-
-  PlxFit fit;
-  if (!FitPMandPar (&fit, X, dX, Y, dY, t, pX, pY, n, VERBOSE)) {
+  PlxFitData fitdata;
+  PlxFitDataAlloc (&fitdata, Ntotal);
+  PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
+
+  PlxFit fit; memset (&fit, 0, sizeof(PlxFit));
+
+  // determine dPsig for detections based on Noutlier attempts
+  if (Noutlier) {
+    int clipRetry = TRUE;
+    for (i = 0; clipRetry && (i < 3); i++) {
+      clipRetry = !PlxOutlierClip (&fitdata, mask, Noutlier, dPsigMax, dPvec, VERBOSE);
+
+      // using the new mask values, reset fitdata
+      PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
+      PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
+      if (VERBOSE) fprintf (stderr, "keep %d of %d\n", fitdata.Npts, Ntotal);
+    }
+  }
+
+  for (i = 0; (VERBOSE == 2) && (i < fitdata.Npts); i++) {
+    int n = fitdata.index[i];
+    int maskValue = mask ? mask[n] : 0;
+    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]);
+  }
+
+  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;
   }
+
+  if (Nresample){
+    PlxFitData sample;
+    PlxFitDataAlloc (&sample, fitdata.Npts);
+
+    PlxFit *testfit = NULL;
+    ALLOCATE (testfit, PlxFit, Nresample);
+
+    int Ngood = 0;
+    for (i = 0; i < Nresample; i++) {
+      PlxBootstrapResample (&fitdata, &sample);
+      
+      if (i % 100000 == 99999) fprintf (stderr, ".");
+
+      // fit the sample
+      testfit[Ngood].getChisq = FALSE;
+      if (!FitPMandPar (&testfit[Ngood], 
+			sample.X, sample.dX, 
+			sample.Y, sample.dY, sample.t, 
+			sample.pX, sample.pY, sample.Npts, VERBOSE)) continue;
+      Ngood ++;
+    }
+
+    Vector *pvec, *uRvec, *uDvec, *Rvec, *Dvec;
+
+    // save the Nresample histograms
+    if ((pvec  = SelectVector ("plxVector", ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "plxVector");
+    if ((uRvec = SelectVector ("uRVector",  ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "uDVector");
+    if ((uDvec = SelectVector ("uDVector",  ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "uRVector");
+    if ((Rvec  = SelectVector ("RoVector",  ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "RoVector");
+    if ((Dvec  = SelectVector ("DoVector",  ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "DoVector");
+    
+    ResetVector ( pvec, OPIHI_FLT, Ngood);
+    ResetVector (uRvec, OPIHI_FLT, Ngood);
+    ResetVector (uDvec, OPIHI_FLT, Ngood);
+    ResetVector ( Rvec, OPIHI_FLT, Ngood);
+    ResetVector ( Dvec, OPIHI_FLT, Ngood);
+    
+    for (i = 0; i < Ngood; i++) {
+      pvec->elements.Flt[i]  = testfit[i].p;
+      uRvec->elements.Flt[i] = testfit[i].uR;
+      uDvec->elements.Flt[i] = testfit[i].uD;
+      Rvec->elements.Flt[i]  = testfit[i].Ro;
+      Dvec->elements.Flt[i]  = testfit[i].Do;
+    }
+
+    // now calculate median and sigma for each vector
+    VectorRobustStats (pvec,  &fit.p,  &fit.dp);
+    VectorRobustStats (uRvec, &fit.uR, &fit.duR);
+    VectorRobustStats (uDvec, &fit.uD, &fit.duD);
+    VectorRobustStats (Rvec,  &fit.Ro, &fit.dRo);
+    VectorRobustStats (Dvec,  &fit.Do, &fit.dDo);
+  }
+
+  // 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 +244,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);
@@ -164,12 +256,10 @@
   set_variable ("uR",   fit.uR);
   set_variable ("uD",   fit.uD);
-  set_variable ("duR",   fit.duR);
-  set_variable ("duD",   fit.duD);
+  set_variable ("duR",  fit.duR);
+  set_variable ("duD",  fit.duD);
   set_variable ("plx",  fit.p);
   set_variable ("dplx", fit.dp);
   
   set_variable ("Tmean",  Tmean);
-  set_variable ("Trange", Trange);
-  set_variable ("Prange", parRange);
 
   set_variable ("chisq", fit.chisq);
@@ -293,19 +383,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 +441,289 @@
   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);
+  ALLOCATE (data->index, int, 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);
+  FREE (data->index);
+}
+
+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;
+}
+
+int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal) {
+
+  int i;
+
+  // find mean values to remove
+  double Nmean = 0;
+  *Tmean = 0;
+  *Rmean = 0;
+  *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;
+  
+  double Trange = Tmax - Tmin;
+
+  // fprintf (stderr, "R,D : %f,%f, T: %f, Trange: %f, Tmin: %f, Tmax: %f\n", *Rmean, *Dmean, *Tmean, Trange, Tmin, Tmax);
+
+  set_variable ("Trange", Trange);
+  return TRUE;
+}
+
+// generate the fit values (projected X,Y; parallax factors; 
+int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean) {
+
+  int i;
+
+  float pXmin = +2.0;
+  float pXmax = -2.0;
+  float pYmin = +2.0;
+  float pYmax = -2.0;
+
+  int Nsubset = 0;
+  for (i = 0; i < Ntotal; i++) {
+    if (mask && !mask[i]) continue;
+    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]);
+    fitdata->index[Nsubset] = i;
+    Nsubset++;
+  }
+  fitdata->Npts = Nsubset;
+  float dXRange = pXmax - pXmin;
+  float dYRange = pYmax - pYmin;
+  float parRange = hypot (dXRange, dYRange);
+
+  set_variable ("Prange", parRange);
+  // fprintf (stderr, "par factor range: %f\n", parRange);
+
+  return TRUE;
+}
+
+/* Outlier clipping based on bootstrap-resampling tests of the plx path
+ * generate Noutlier resampled datasets
+ * fit the Noutlier plx paths
+ * determine and save the distribution of dXsig and dYsig for each point
+ * sort the resulting distributions and find dPsig (median point) for each measurement
+ * find the 90% point of dPsig : if > dPsigMax, only clip the 10% most deviant points
+ * set the dPvec values if desired
+ * -- mask is modified, dPvec values are set
+ * -- fitdata is unchanged
+ */
+
+# define MAX_REJECT 0.1
+
+int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE) {
+
+  int i, n;
+
+  PlxFit testfit;
+  testfit.getChisq = FALSE;
+
+  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);
+  }
+
+  int Nsamples = 0;
+  for (n = 0; n < Noutlier; n++) {
+    // bootstrap resample (fitdata -> sample)
+    PlxBootstrapResample (fitdata, &sample);
+      
+    if (n % 100000 == 99999) 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]);
+  }
+
+  // make a copy of dPsig[] and check if > 10% are > dPsigMax
+  double *dPsigSort;
+  ALLOCATE (dPsigSort, double, fitdata->Npts);
+  for (i = 0; i < fitdata->Npts; i++) {
+    dPsigSort[i] = dPsig[i];
+  }
+  dsort (dPsigSort, fitdata->Npts);
+  int Nmax = (1.0 - MAX_REJECT)*fitdata->Npts;
+
+  int completeClip = TRUE;
+  if (dPsigSort[Nmax] > dPsigMax) {
+    if (VERBOSE) fprintf (stderr, "too many outliers: %f at 90\n", dPsigSort[Nmax]);
+    dPsigMax = dPsigSort[Nmax];
+    completeClip = FALSE;
+  }
+
+  for (i = 0; i < fitdata->Npts; i++) {
+    if (dPsig[i] < dPsigMax) continue;
+    int n = fitdata->index[i];
+    // fprintf (stderr, "clip %d: %f : %f\n", i, fitdata->t[i], dPsig[i]);
+    mask[n] = 0; // mask these points
+  }
+
+  // only set dPvec if we have completed the clipping?
+  if (dPvec) {
+    for (i = 0; i < dPvec->Nelements; i++) {
+      dPvec->elements.Flt[i] = NAN;
+    }
+    for (i = 0; i < fitdata->Npts; i++) {
+      int n = fitdata->index[i];
+      dPvec->elements.Flt[n] = dPsig[i];
+    }
+  }
+
+  free (dPsig);
+  free (dPsigSort);
+  
+  for (i = 0; i < fitdata->Npts; i++) {
+    free (dXsig[i]);
+    free (dYsig[i]);
+  }
+  free (dXsig);
+  free (dYsig);
+
+  return completeClip;
+}
+
+int VectorRobustStats (Vector *vector, double *median, double *sigma) {
+
+  // warn if vector->Nelements > 1000? 10000?)
+  // warn if vector is not float
+
+  // we need to copy the vector to avoid changing the sort order
+  double *values = NULL;
+  ALLOCATE (values, double, vector->Nelements);
+
+  int i;
+  int Npts = 0;
+  for (i = 0; i < vector->Nelements; i++) {
+    if (!isfinite(vector->elements.Flt[i])) continue;
+    values[Npts] = vector->elements.Flt[i];
+    Npts++;
+  }
+
+  dsort (values, Npts);
+
+  if (Npts % 2) {
+    int Ncenter = Npts / 2;
+    *median = values[Ncenter];
+  } else {
+    int Ncenter = Npts / 2 - 1;
+    *median = 0.5*(values[Ncenter] + values[Ncenter + 1]);
+  }
+
+  double Slo = VectorFractionInterpolate (values, 0.158655, Npts);
+  double Shi = VectorFractionInterpolate (values, 0.841345, Npts);
+
+  *sigma = (Shi - Slo) / 2.0;
+
+  return TRUE;
+}
+
+double VectorFractionInterpolate (double *values, float fraction, int Npts) {
+
+  float F = fraction * Npts;
+  int   N = fraction * Npts;
+
+  if (N < 0        ) return NAN;
+  if (N >= Npts - 2) return NAN;
+
+  // interpolate between N,N+1
+    
+  double S = (F - N) * (values[N+1] - values[N]) + values[N];
+  return S;
+}
Index: trunk/Ohana/src/opihi/cmd.basic/date.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.basic/date.c	(revision 38972)
+++ trunk/Ohana/src/opihi/cmd.basic/date.c	(revision 38986)
@@ -3,5 +3,6 @@
 int date (int argc, char **argv) {
   
-  int N, SECONDS, REFTIME;
+  int N, SECONDS;
+  double REFTIME;
   struct timeval now;
   char *tstring = NULL;
@@ -16,8 +17,8 @@
   }
 
-  REFTIME = 0;
+  REFTIME = 0.0;
   if ((N = get_argument (argc, argv, "-reftime"))) {
     remove_argument (N, &argc, argv);
-    REFTIME = atoi (argv[N]);
+    REFTIME = atof (argv[N]);
     remove_argument (N, &argc, argv);
   } 
@@ -37,8 +38,9 @@
   gettimeofday (&now, NULL);
   if (SECONDS) {
+    double nowSec = now.tv_sec + 1e-6*now.tv_usec;
     if (varName) {
-      set_int_variable (varName, now.tv_sec - REFTIME);
+      set_variable (varName, nowSec - REFTIME);
     } else {
-      gprint (GP_ERR, "%d\n", (int) now.tv_sec - REFTIME);
+      gprint (GP_ERR, "%.12g\n", nowSec - REFTIME);
     }
   } else {
Index: trunk/Ohana/src/opihi/cmd.data/box.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/box.c	(revision 38972)
+++ trunk/Ohana/src/opihi/cmd.data/box.c	(revision 38986)
@@ -48,5 +48,5 @@
     if (strlen (graphmode.ticks) != 4) { goto usage; }
     for (i = 0; i < strlen (graphmode.ticks); i++) {
-      if ((graphmode.ticks[i] != '0') && (graphmode.ticks[i] != '1') && (graphmode.ticks[i] != '2')) { goto usage; }
+      if ((graphmode.ticks[i] != '0') && (graphmode.ticks[i] != '1') && (graphmode.ticks[i] != '2') && (graphmode.ticks[i] != '3')) { goto usage; }
     }
   }
@@ -139,4 +139,62 @@
     remove_argument (N, &argc, argv);
     graphmode.padYp = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if ((N = get_argument (argc, argv, "-fminor"))) {
+    remove_argument (N, &argc, argv);
+    graphmode.fMinorXm = atof(argv[N]);
+    graphmode.fMinorXp = atof(argv[N]);
+    graphmode.fMinorYm = atof(argv[N]);
+    graphmode.fMinorYp = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-xfminor"))) {
+    remove_argument (N, &argc, argv);
+    graphmode.fMinorXm = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "+xfminor"))) {
+    remove_argument (N, &argc, argv);
+    graphmode.fMinorXp = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-yfminor"))) {
+    remove_argument (N, &argc, argv);
+    graphmode.fMinorYm = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "+yfminor"))) {
+    remove_argument (N, &argc, argv);
+    graphmode.fMinorYp = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  if ((N = get_argument (argc, argv, "-flabel"))) {
+    remove_argument (N, &argc, argv);
+    graphmode.fLabelRangeXm = atof(argv[N]);
+    graphmode.fLabelRangeXp = atof(argv[N]);
+    graphmode.fLabelRangeYm = atof(argv[N]);
+    graphmode.fLabelRangeYp = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-xflabel"))) {
+    remove_argument (N, &argc, argv);
+    graphmode.fLabelRangeXm = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "+xflabel"))) {
+    remove_argument (N, &argc, argv);
+    graphmode.fLabelRangeXp = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-yflabel"))) {
+    remove_argument (N, &argc, argv);
+    graphmode.fLabelRangeYm = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "+yflabel"))) {
+    remove_argument (N, &argc, argv);
+    graphmode.fLabelRangeYp = atof(argv[N]);
     remove_argument (N, &argc, argv);
   }
@@ -174,4 +232,12 @@
   gprint (GP_ERR, "         alternatively, set each axis independently with:\n");
   gprint (GP_ERR, "        -xpad, -ypad, +xpad, +ypad\n");
+  gprint (GP_ERR, "  \n");
+  gprint (GP_ERR, "  -fminor : set the number of minor ticks per major (all axes)\n");
+  gprint (GP_ERR, "         alternatively, set each axis independently with:\n");
+  gprint (GP_ERR, "        -xfminor, -yfminor, +xfminor, +yfminor\n");
+  gprint (GP_ERR, "  \n");
+  gprint (GP_ERR, "  -flabel : set the fraction of axis over which major ticks have labels (all axes)\n");
+  gprint (GP_ERR, "         alternatively, set each axis independently with:\n");
+  gprint (GP_ERR, "        -xflabel, -yflabel, +xflabel, +yflabel\n");
 
   return (FALSE);
Index: trunk/Ohana/src/opihi/cmd.data/device.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/device.c	(revision 38972)
+++ trunk/Ohana/src/opihi/cmd.data/device.c	(revision 38986)
@@ -6,4 +6,10 @@
   char *name;;
   /* set / get current graphics device */
+
+  int QUIET = FALSE;
+  if ((N = get_argument (argc, argv, "-q"))) {
+    remove_argument (N, &argc, argv);
+    QUIET = TRUE;
+  }
 
   name = NULL;
@@ -23,5 +29,5 @@
     if (!GetGraph (NULL, &kapa, name)) return (FALSE);
   }
-  gprint (GP_ERR, "kapa %s\n", name); 
+  if (!QUIET) gprint (GP_ERR, "kapa %s\n", name); 
 
   return (TRUE);
Index: trunk/Ohana/src/opihi/cmd.data/read_vectors.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/read_vectors.c	(revision 38972)
+++ trunk/Ohana/src/opihi/cmd.data/read_vectors.c	(revision 38986)
@@ -511,5 +511,6 @@
     if (!gfits_fread_ftable_data (f, &table, padIfShort)) ESCAPE ("error reading table for extension %d\n", Nextend);
   } else {
-    if (!gfits_fread_ftable_range (f, padIfShort, &table, start, Nrows)) ESCAPE ("error reading table for extension %d\n", Nextend);
+    // arg3 (FALSE) : seek to this segment start
+    if (!gfits_fread_ftable_range (f, padIfShort, FALSE, &table, start, Nrows)) ESCAPE ("error reading table for extension %d\n", Nextend);
   }
 
Index: trunk/Ohana/src/opihi/cmd.data/uniq.c
===================================================================
--- trunk/Ohana/src/opihi/cmd.data/uniq.c	(revision 38972)
+++ trunk/Ohana/src/opihi/cmd.data/uniq.c	(revision 38986)
@@ -1,11 +1,24 @@
 # include "data.h"
+// NOTE: if there are only a few uniq values, the old algorithm is not bad.  
+// for 10000 uniq values, 30M points take ~20sec in the new algorithm, 
+// 3M points takess 45 sec in the old method.
 
 int uniq (int argc, char **argv) {
   
-  int Nnew, i, j, found;
+  int Nnew, i, N;
   Vector *ivec, *ovec;
 
+  Vector *cvec = NULL;
+  if ((N = get_argument (argc, argv, "-c"))) {
+    remove_argument (N, &argc, argv);
+    if ((cvec = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) {
+      gprint (GP_ERR, "invalid vector %s\n", argv[N]);
+      return FALSE;
+    }
+    remove_argument (N, &argc, argv);
+  }
+
   if (argc != 3) {
-    gprint (GP_ERR, "USAGE: uniq (in) (out)\n");
+    gprint (GP_ERR, "USAGE: uniq (in) (out) -c count\n");
     return (FALSE);
   }
@@ -16,37 +29,69 @@
   /* allocate the maximum possible needed */
   ResetVector (ovec, ivec->type, ivec->Nelements);
+  if (cvec) {
+    ResetVector (cvec, OPIHI_INT, ivec->Nelements);
+  }
 
   Nnew = 0;
 
   if (ivec->type == OPIHI_FLT) {
-    opihi_flt *v1 = ivec[0].elements.Flt;
-    for (i = 0; i < ivec[0].Nelements; i++, v1++) {
-      opihi_flt *v2 = ovec[0].elements.Flt;
-      found = FALSE;
-      for (j = 0; !found && (j < Nnew); j++, v2++) {
-	if (*v1 == *v2) found = TRUE;
+    // copy the input data to a temporary array to avoid damaging it with sort
+    opihi_flt *indata = NULL;
+    ALLOCATE (indata, opihi_flt, ivec[0].Nelements);
+    memcpy (indata, ivec->elements.Flt, ivec[0].Nelements*sizeof(opihi_flt));
+
+    dsort (indata, ivec->Nelements);
+
+    Nnew = 0;
+    opihi_flt *vtgt = ovec[0].elements.Flt;
+
+    opihi_flt *vsrc = indata;
+
+    for (i = 0; i < ivec->Nelements; Nnew++) {
+      vtgt[Nnew] = *vsrc;
+      int Ndup = 0;
+      opihi_flt lastValue = *vsrc;
+      while ((i < ivec->Nelements) && (*vsrc == lastValue)) {
+	i++;
+	vsrc ++;
+	Ndup ++;
       }
-      if (!found) {
-	ovec[0].elements.Flt[Nnew] = *v1;
-	Nnew ++;
+      if (cvec) {
+	cvec->elements.Int[Nnew] = Ndup;
       }
     }
+    free (indata);
   } else {
-    opihi_int *v1 = ivec[0].elements.Int;
-    for (i = 0; i < ivec[0].Nelements; i++, v1++) {
-      opihi_int *v2 = ovec[0].elements.Int;
-      found = FALSE;
-      for (j = 0; !found && (j < Nnew); j++, v2++) {
-	if (*v1 == *v2) found = TRUE;
+    // copy the input data to a temporary array to avoid damaging it with sort
+    opihi_int *indata = NULL;
+    ALLOCATE (indata, opihi_int, ivec[0].Nelements);
+    memcpy (indata, ivec->elements.Int, ivec[0].Nelements*sizeof(opihi_int));
+
+    isort (indata, ivec->Nelements);
+
+    Nnew = 0;
+    opihi_int *vtgt = ovec[0].elements.Int;
+
+    opihi_int *vsrc = indata;
+
+    for (i = 0; i < ivec->Nelements; Nnew++) {
+      vtgt[Nnew] = *vsrc;
+      int Ndup = 0;
+      opihi_int lastValue = *vsrc;
+      while ((i < ivec->Nelements) && (*vsrc == lastValue)) {
+	i++;
+	vsrc ++;
+	Ndup ++;
       }
-      if (!found) {
-	ovec[0].elements.Int[Nnew] = *v1;
-	Nnew ++;
+      if (cvec) {
+	cvec->elements.Int[Nnew] = Ndup;
       }
     }
+    free (indata);
   }
 
   // free up extra memory
   ResetVector (ovec, ivec->type, Nnew);
+  if (cvec) ResetVector (cvec, OPIHI_INT, Nnew);
 
   return (TRUE);
Index: trunk/Ohana/src/opihi/dvo/mextract.c
===================================================================
--- trunk/Ohana/src/opihi/dvo/mextract.c	(revision 38972)
+++ trunk/Ohana/src/opihi/dvo/mextract.c	(revision 38986)
@@ -213,5 +213,5 @@
   }
 
-  // int needLensing = dbFieldNeedLensing (fields, Nfields);
+  int needLensing = dbFieldNeedLensing (fields, Nfields);
   int needStarpar = dbFieldNeedStarpar (fields, Nfields, FALSE);
 
@@ -234,5 +234,5 @@
     catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i];
     catalog.catflags = DVO_LOAD_AVERAGE | DVO_LOAD_MEASURE | DVO_LOAD_SECFILT;
-    // catalog.catflags |= needLensing ? DVO_LOAD_LENSING : DVO_SKIP_LENSING;
+    catalog.catflags |= needLensing ? DVO_LOAD_LENSING : DVO_SKIP_LENSING;
     catalog.catflags |= needStarpar ? DVO_LOAD_STARPAR : DVO_SKIP_STARPAR;
     catalog.Nsecfilt = Nsecfilt;
@@ -274,6 +274,6 @@
 	StarPar *starpar = needStarpar ? &catalog.starpar[Nstarpar] : NULL;
 
-	// int Nlensing = average->lensobjOffset;
-	// Lensobj *lensobj = needLensobj ? &catalog.lensobj[m] : NULL;
+	int Nlensing = average->lensingOffset;
+	Lensing *lensing = needLensing ? &catalog.lensing[Nlensing] : NULL;
 
 	int Nsec = j*Nsecfilt;
@@ -281,5 +281,5 @@
 
 	for (n = 0; n < Nfields; n++) {
-	  values[n] = dbExtractMeasures (average, secfilt, &catalog.measure[m], NULL, starpar, &fields[n]);
+	  values[n] = dbExtractMeasures (average, secfilt, &catalog.measure[m], lensing, starpar, &fields[n]);
 	}
 	// fprintf (stderr, "object: ave: %f, cat: %f, averef %d\n", fields[n].name, values[2], values[3], catalog.measure[m].averef);
Index: trunk/Ohana/src/opihi/lib.data/gaussian.c
===================================================================
--- trunk/Ohana/src/opihi/lib.data/gaussian.c	(revision 38972)
+++ trunk/Ohana/src/opihi/lib.data/gaussian.c	(revision 38986)
@@ -20,5 +20,4 @@
  
   int i;
-  long A, B;
   double val, x, dx, dx1, dx2, dx3, df;
   double mean, sigma;
@@ -27,7 +26,6 @@
   if (Ngaussint == Nbin) return;
 
-  A = time(NULL);
-  for (B = 0; A == time(NULL); B++);
-  srand48(B);
+  long A = time(NULL);
+  srand48(A);
  
   Ngaussint = Nbin;
Index: trunk/Ohana/src/opihi/lib.shell/ConfigInit.c
===================================================================
--- trunk/Ohana/src/opihi/lib.shell/ConfigInit.c	(revision 38972)
+++ trunk/Ohana/src/opihi/lib.shell/ConfigInit.c	(revision 38986)
@@ -38,4 +38,9 @@
   }
 
+  // XXX this is a bit dangerous : answer may have an arbitrary
+  // length, but we do not know the available space in 'ptr'
+  // we really should be passing in the buffer size and limiting the copy
+
+
   if (!strcmp (mode, "%s"))  strcpy ((char *) ptr, answer);
   if (!strcmp (mode, "%d"))  *(int *) ptr       = atoi (answer);
Index: trunk/Ohana/src/opihi/lib.shell/SocketOps.c
===================================================================
--- trunk/Ohana/src/opihi/lib.shell/SocketOps.c	(revision 38972)
+++ trunk/Ohana/src/opihi/lib.shell/SocketOps.c	(revision 38986)
@@ -23,5 +23,9 @@
   status = gethostname (myHostname, HOST_NAME_MAX);
 
-  fprintf (stderr, "target host: %s, real host: %s\n", hostname, myHostname);
+  if (strcmp (hostname, myHostname)) {
+    fprintf (stderr, "target host: %s, real host: %s\n", hostname, myHostname);
+    fprintf (stderr, "please run on the correct host\n");
+    exit (2);
+  }
 
   GetPortRange (&start, &stop, portinfo);
@@ -197,4 +201,9 @@
 
   host = gethostbyname (hostname);
+  if (!host) {
+    gprint (GP_ERR, "cannot connect to pantasks server %s\n", hostname);
+    exit (3);
+  }
+
   bzero (hostip, 80);
   for (i = 0; i < host[0].h_length; i++) {
Index: trunk/Ohana/src/opihi/lib.shell/gprint.c
===================================================================
--- trunk/Ohana/src/opihi/lib.shell/gprint.c	(revision 38972)
+++ trunk/Ohana/src/opihi/lib.shell/gprint.c	(revision 38986)
@@ -255,5 +255,5 @@
   if (file == NULL) {
     // XXX this is a problem: we are leaving open the old file
-    fprintf (stderr, "cannot open file %s\n", stream[0].name);
+    fprintf (stderr, "gprint cannot open file %s\n", stream[0].name);
     free (stream[0].name);
     file = (dest == GP_LOG) ? stdout : stderr;
Index: trunk/Ohana/src/opihi/pantasks/CheckJobs.c
===================================================================
--- trunk/Ohana/src/opihi/pantasks/CheckJobs.c	(revision 38972)
+++ trunk/Ohana/src/opihi/pantasks/CheckJobs.c	(revision 38986)
@@ -1,3 +1,4 @@
 # include "pantasks.h"
+static int Ncheck = 0;
 
 float CheckJobs () {
@@ -12,6 +13,5 @@
   float time_running, next_timeout;
 
-  // int Ncheck;
-  // Ncheck = 0;
+  Ncheck ++;
 
   // actual maximum delay is controlled in job_threads.c
@@ -21,5 +21,4 @@
   /** test all jobs: ready to test?  finished? **/
   while ((job = NextJob ()) != NULL) {
-    // Ncheck ++;
 
     task = job[0].task;
@@ -220,5 +219,5 @@
     SetTaskTimer (&job[0].last);
   }
-  // fprintf (stderr, "check %d jobs\n", Ncheck);
+
   JobTaskUnlock();
   return (next_timeout);
Index: trunk/Ohana/src/opihi/pantasks/CheckTasks.c
===================================================================
--- trunk/Ohana/src/opihi/pantasks/CheckTasks.c	(revision 38972)
+++ trunk/Ohana/src/opihi/pantasks/CheckTasks.c	(revision 38986)
@@ -1,3 +1,4 @@
 # include "pantasks.h"
+static int Ncheck = 0; 
 
 float CheckTasks () {
@@ -7,4 +8,6 @@
   int status;
   float time_running, next_timeout, fuzz;
+
+  Ncheck ++;
 
   // actual maximum delay is controlled in job_threads.c
Index: trunk/Ohana/src/opihi/pantasks/TaskOps.c
===================================================================
--- trunk/Ohana/src/opihi/pantasks/TaskOps.c	(revision 38972)
+++ trunk/Ohana/src/opihi/pantasks/TaskOps.c	(revision 38986)
@@ -22,5 +22,7 @@
 void FreeTasks () {
   int i;
-  for (i = 0; i < Ntasks; i++) {
+  int ntasks = Ntasks;
+  Ntasks = 0;
+  for (i = 0; i < ntasks; i++) {
     FreeTask (tasks[i]);
   }
Index: trunk/Ohana/src/opihi/pantasks/pantasks_server.c.in
===================================================================
--- trunk/Ohana/src/opihi/pantasks/pantasks_server.c.in	(revision 38972)
+++ trunk/Ohana/src/opihi/pantasks/pantasks_server.c.in	(revision 38986)
@@ -19,5 +19,5 @@
   
   char hostname[256], portinfo[256];
-  char log_stdout[1024], log_stderr[1024];
+  char log_stdout[1024], log_stderr[1024], tmpname[1024];
   pthread_t JobsAndTasksThread;
   pthread_t clientsThread;
@@ -34,6 +34,7 @@
   stdin = freopen ("/dev/zero", "r", stdin);
 
-  // this block redirects the actual stderr, stdout steams to the output files
-  if (VarConfig ("PANTASKS_SERVER_STDOUT", "%s", log_stdout) != NULL) {
+  // this block redirects the actual stderr, stdout streams to the output files
+  if (VarConfig ("PANTASKS_SERVER_STDOUT", "%s", tmpname) != NULL) {
+    strcpy (log_stdout, tmpname);
     if (strcmp (log_stdout, "stdout")) {
       stdout = freopen (log_stdout, "a", stdout);
@@ -44,5 +45,6 @@
     }
   }
-  if (VarConfig ("PANTASKS_SERVER_STDERR", "%s", log_stderr) != NULL) {
+  if (VarConfig ("PANTASKS_SERVER_STDERR", "%s", tmpname) != NULL) {
+    strcpy (log_stderr, tmpname);
     if (strcmp (log_stderr, "stderr")) {
       stderr = freopen (log_stderr, "a", stderr);
Index: trunk/Ohana/src/opihi/pantasks/test/local.sh
===================================================================
--- trunk/Ohana/src/opihi/pantasks/test/local.sh	(revision 38972)
+++ trunk/Ohana/src/opihi/pantasks/test/local.sh	(revision 38986)
@@ -1,4 +1,5 @@
 
 ## a basic test of memory allocation
+if (not($?VERBOSE)) set VERBOSE = 0
 
 exec rm -f tmp.txt
@@ -38,4 +39,5 @@
  $startmem = $word:1
 end
+
 macro memcheck
  list word -x "ps -p $PID -o rss"