Changeset 39236 for trunk

Timestamp:

Dec 7, 2015, 7:52:53 PM (11 years ago)

Author:

eugene

Message:

fix fitplx_irls to generate accurate plx fits

File:

• trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c (modified) (12 diffs)

trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c

@@ -25 +25 @@
   }
 
-  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"))) {
@@ -42 +32 @@
   }
 
-  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);
-  }
-
   int max_iterations = 10;
   if ((N = get_argument (argc, argv, "-max-iterations"))) {
@@ -57 +39 @@
   }
 
-  Vector *outMask = NULL;
-  if ((N = get_argument (argc, argv, "-out-mask"))) {
-    remove_argument (N, &argc, argv);
-    if (!(outMask = SelectVector(argv[N], ANYVECTOR, TRUE))) return FALSE;
-    remove_argument (N, &argc, argv);
-  }
-
   double outlier_limit = 0.1;
   if ((N = get_argument (argc, argv, "-outlier-limit"))) {
@@ -74 +49 @@
   if (argc != 6) {
     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");
     gprint (GP_ERR, "  -max-iterations : maximum number of IRLS iterations to run (default 10)\n");
-    gprint (GP_ERR, "  -outMask vect : save outlier mask in vec\n");
     gprint (GP_ERR, "  -outlier-limit : fraction of average weight to reject on (default 0.1)\n");
     return (FALSE);
@@ -109 +81 @@
   // Nsubset : unmasked points
   int Ntotal = tvec->Nelements; // XXX check other lengths
-  if (dPvec) ResetVector (dPvec, OPIHI_FLT, Ntotal);
-  if (outMask) ResetVector (outMask, OPIHI_INT, Ntotal);
   
   double Rmean, Dmean, Tmean;
@@ -125 +95 @@
   PlxFitDataAlloc (&fitdata, Ntotal);
   PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
+  // fitdata only contains the points which are not masked
 
   PlxFit fit; memset (&fit, 0, sizeof(PlxFit));
 
-  if (0) {
-    // This should be accomplished via the IRLS call
-    // 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;
+    int maskValue = mask ? mask[n] : 1;
     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]);
   }
@@ -159 +114 @@
     return FALSE;
   }
-
-  // Set the output mask based on the input mask and the outlier limits.
-  if (outMask) {
-    double Sum_Wx = 0;
-    double Sum_Wy = 0;
-    int *outMask_V = outMask->elements.Int;
-    for (i = 0; i < fitdata.Npts; i++) {
-      Sum_Wx += fitdata.Wx[i];
-      Sum_Wy += fitdata.Wy[i];
-    }
-    for (i = 0; i < Ntotal; i++) {
-      outMask_V[i] = mask ? mask[i] : 0;
-    }
-    
-    for (i = 0; i < fitdata.Npts; i++) {
-      if ((fitdata.Wx[i] > outlier_limit * Sum_Wx / (1.0 * fitdata.Npts))||
-          (fitdata.Wy[i] > outlier_limit * Sum_Wy / (1.0 * fitdata.Npts))) {
-        int n = fitdata.index[i];
-        outMask_V[n] = 1;
-
-        if (VERBOSE == 2) {
-          fprintf (stderr, "%f %f : %f %d : %f %f %f : %f %f %f %f\n", R[n], D[n], T[n], outMask_V[n], fitdata.t[i], fitdata.X[i], fitdata.Y[i], fitdata.Wx[i], fitdata.Wy[i], Sum_Wx, Sum_Wy);
-        }
-      }
-    }
-  }
-
-  
+  // FitPMandPar_IRLS sets the values of Wx,Wy based on the fit distance
+
+  // update the mask based on the input mask and the outlier limits.
+  double Sum_Wx = 0;
+  double Sum_Wy = 0;
+  mask = mvec->elements.Int;
+  for (i = 0; i < fitdata.Npts; i++) {
+    Sum_Wx += fitdata.Wx[i];
+    Sum_Wy += fitdata.Wy[i];
+  }
+  for (i = 0; i < fitdata.Npts; i++) {
+    // fitdata only includes the previously unmasked points
+    if ((fitdata.Wx[i] < outlier_limit * Sum_Wx / (1.0 * fitdata.Npts))||
+        (fitdata.Wy[i] < outlier_limit * Sum_Wy / (1.0 * fitdata.Npts))) {
+      int n = fitdata.index[i];
+      mask[n] = 0;
+      
+      if (VERBOSE == 2) {
+          fprintf (stderr, "%f %f : %f %d : %f %f %f : %f %f %f %f\n", R[n], D[n], T[n], mask[n], fitdata.t[i], fitdata.X[i], fitdata.Y[i], fitdata.Wx[i], fitdata.Wy[i], Sum_Wx, Sum_Wy);
+      }
+    }
+  }
+
   if (Nresample){
+    // now that the mask has been updated, we need to recalculate mean epoch and positions
+    // XXX make this conditional on actually masking unmasked points above
+    PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
+    PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
+
     PlxFitData sample;
     PlxFitDataAlloc (&sample, fitdata.Npts);
@@ -319 +274 @@
   }
 
+  // XXX EAM : I suggest making these static as well like the A and B matrices
   // Allocate the IRLS specific objects
   ALLOCATE (Cov, double *, 5);
@@ -327 +283 @@
   ALLOCATE(Beta, double, 5);
   ALLOCATE(Beta_prev, double, 5);
-  //  ALLOCATE(Wx, double, Npts);
-  //  ALLOCATE(Wy, double, Npts);
   ALLOCATE(rx,  double, Npts);
   ALLOCATE(ry,  double, Npts);
@@ -348 +302 @@
   sigma_ols = 0.0;
   for (i = 0; i < Npts; i++) {
+    // double Xf = B[0][0] + B[1][0]*T[i] + B[4][0]*pR[i];
+    // double Yf = B[2][0] + B[3][0]*T[i] + B[4][0]*pD[i];
+
     rx[i] = X[i] - (T[i] * B[1][0] + B[0][0] + B[4][0] * pR[i]);
     ry[i] = Y[i] - (T[i] * B[3][0] + B[2][0] + B[4][0] * pD[i]);
+
+
     //    u[i] = r[i] /
     sigma_ols += SQ(rx[i]) + SQ(ry[i]);
-
   }
   sigma_ols = sqrt(sigma_ols / (Npts - 5.0));
@@ -390 +348 @@
 
     // Solve
-    if (!weighted_LS_PLX(T,pD,pR,X,Wx,Y,Wy,Npts,A,B,VERBOSE)) {
+    if (!weighted_LS_PLX(T,pR,pD,X,Wx,Y,Wy,Npts,A,B,VERBOSE)) {
       // Handle fail case
       return(FALSE);
@@ -478 +436 @@
     }
   }
-  
-  
 
   // Finish