Changeset 39181

Timestamp:

Nov 20, 2015, 3:33:29 PM (11 years ago)

Author:

watersc1

Message:

Updated version that can set an output mask based on the weight values, which are passed back up.

File:

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

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

@@ -15 +15 @@
   double *dY;
 
-  double *IRLS_u;
+  double *Wx;
+  double *Wy;
+  
   int *index;
   int Npts;
@@ -32 +34 @@
   int Nfit;
   int getChisq;
-} PlxFit_IRLS;
+} PlxFit;
+
+
+int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE);
 
 int VectorRobustStats_IRLS (Vector *vector, double *median, double *sigma);
@@ -45 +50 @@
 int PlxBootstrapResample_IRLS (PlxFitData *src, PlxFitData *tgt);
 
-int FitPMandPar_IRLS (PlxFit_IRLS *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, double *IRLS_u, int Npts, int max_iterations, int VERBOSE);
+int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, double *Wx, double *Wy, int Npts, int max_iterations, double outlier_limit, int VERBOSE);
 int sun_ecliptic_IRLS (double mjd, double *lambda, double *beta, double *epsilon, double *Radius);
 int ParFactor_IRLS (double *pR, double *pD, double RA, double DEC, double Time);
 
-int IRLS_converged (PlxFit_IRLS *fit);
+int IRLS_converged (PlxFit *fit);
 int Plx_weighted_LS (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts,
                      double **A, double **B, int VERBOSE);
@@ -114 +119 @@
     remove_argument (N, &argc, argv);
   }
+
+  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"))) {
+    remove_argument (N, &argc, argv);
+    outlier_limit = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+        
   
   if (argc != 6) {
@@ -119 +139 @@
     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\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);
   }
@@ -151 +173 @@
   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;
   PlxSetMeanEpoch_IRLS (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
@@ -166 +189 @@
   PlxSetEpochPosition_IRLS (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
 
-  PlxFit_IRLS fit; memset (&fit, 0, sizeof(PlxFit_IRLS));
+  PlxFit fit; memset (&fit, 0, sizeof(PlxFit));
 
   if (0) {
@@ -195 +218 @@
                          fitdata.Y, fitdata.dY, 
                          fitdata.t, fitdata.pX, fitdata.pY,
-                         fitdata.IRLS_u,
-                         fitdata.Npts, max_iterations, VERBOSE)) {
+                         fitdata.Wx, fitdata.Wy,
+                         fitdata.Npts, max_iterations, outlier_limit, VERBOSE)) {
     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);
+        }
+      }
+    }
+  }
+
+  
   if (Nresample){
     PlxFitData sample;
     PlxFitDataAlloc_IRLS (&sample, fitdata.Npts);
 
-    PlxFit_IRLS *testfit = NULL;
-    ALLOCATE (testfit, PlxFit_IRLS, Nresample);
+    PlxFit *testfit = NULL;
+    ALLOCATE (testfit, PlxFit, Nresample);
 
     int Ngood = 0;
@@ -215 +264 @@
       // fit the sample
       testfit[Ngood].getChisq = FALSE;
-      if (!FitPMandPar_IRLS (&testfit[Ngood], 
-                             sample.X, sample.dX, 
-                             sample.Y, sample.dY, sample.t, 
-                             sample.pX, sample.pY,
-                             sample.IRLS_u, sample.Npts, 1, VERBOSE)) continue;
+      if (!FitPMandPar (&testfit[Ngood], 
+                        sample.X, sample.dX, 
+                        sample.Y, sample.dY, sample.t, 
+                        sample.pX, sample.pY,
+                        sample.Npts, VERBOSE)) continue;
       Ngood ++;
     }
@@ -292 +341 @@
   set_variable ("Nfit",  fit.Nfit);
 
+
+  
   return (TRUE);
 }
 
 /* do we want an init function which does the alloc and a clear function to free? */
-int FitPMandPar_IRLS (PlxFit_IRLS *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD,
-                      double *IRLS_u,
-                      int Npts, int max_iterations, int VERBOSE) {
+int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD,
+                      double *Wx, double *Wy,
+                      int Npts, int max_iterations, double outlier_limit, int VERBOSE) {
 
   int i,j;
@@ -310 +361 @@
 
   double sigma_ols, sigma_hat; // Sigma estimates.
-  double *Wx, *Wy;             // Weight vectors.  Not errors.
+  //  double *Wx, *Wy;             // Weight vectors.  Not errors.
   double *rx, *ry;             // Deviation from model
-  //  double *u;                   // Deviation magnitude
+  double *u;                   // Deviation magnitude
   int converged;
   int iterations;
@@ -338 +389 @@
   ALLOCATE(Beta, double, 5);
   ALLOCATE(Beta_prev, double, 5);
-  ALLOCATE(Wx, double, Npts);
-  ALLOCATE(Wy, double, Npts);
+  //  ALLOCATE(Wx, double, Npts);
+  //  ALLOCATE(Wy, double, Npts);
   ALLOCATE(rx,  double, Npts);
   ALLOCATE(ry,  double, Npts);
+  ALLOCATE(u,  double, Npts);
   
   // Convert the measurement errors into the initial weights.
@@ -416 +468 @@
       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]);
-      IRLS_u[i] = sqrt(SQ(rx[i] / dX[i]) + SQ(ry[i] / dY[i]));
+      u[i] = sqrt(SQ(rx[i] / dX[i]) + SQ(ry[i] / dY[i]));
     }
 
     // Calculate sigma_hat from distribution of residual magnitudes
-    sigma_hat = Plx_MAD(IRLS_u,Npts) / 0.6745;
+    sigma_hat = Plx_MAD(u,Npts) / 0.6745;
 
     // Check convergence
@@ -511 +563 @@
     double wx,wy;
     for (i = 0; i < Npts; i++) {
-      if ((Wx[i] > 0.1 * Sum_Wx / (1.0 * Npts))||
-          (Wy[i] > 0.1 * Sum_Wy / (1.0 * Npts))) {
+      if ((Wx[i] > outlier_limit * Sum_Wx / (1.0 * Npts))||
+          (Wy[i] > outlier_limit * Sum_Wy / (1.0 * Npts))) {
         
         Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].p*pR[i];
@@ -584 +636 @@
   ALLOCATE (data->pX, double, N);
   ALLOCATE (data->pY, double, N);
-  ALLOCATE (data->IRLS_u, double, N);
+  ALLOCATE (data->Wx, double, N);
+  ALLOCATE (data->Wy, double, N);
   ALLOCATE (data->index, int, N);
   return TRUE;
@@ -597 +650 @@
   FREE (data->pX);
   FREE (data->pY);
-  FREE (data->IRLS_u);
+  FREE (data->Wx);
+  FREE (data->Wy);
   FREE (data->index);
 }
@@ -614 +668 @@
     tgt->pX[i] = src->pX[N];
     tgt->pY[i] = src->pY[N];
-    tgt->IRLS_u[i] = src->IRLS_u[N];
+    tgt->Wx[i] = src->Wx[N];
+    tgt->Wy[i] = src->Wy[N];
   }
   return TRUE;
@@ -674 +729 @@
     pYmax = MAX (pYmax, fitdata->pY[Nsubset]);
 
-    fitdata->IRLS_u[Nsubset] = 0.0;
+    fitdata->Wx[Nsubset] = 1.0;
+    fitdata->Wy[Nsubset] = 1.0;    
     fitdata->index[Nsubset] = i;
     Nsubset++;
@@ -706 +762 @@
   int i, n;
 
-  PlxFit_IRLS testfit;
+  PlxFit testfit;
   testfit.getChisq = FALSE;
 
@@ -728 +784 @@
 
     // fit the sample
-    if (!FitPMandPar_IRLS (&testfit, 
-                           sample.X, sample.dX, 
-                           sample.Y, sample.dY, sample.t, 
-                           sample.pX, sample.pY, sample.IRLS_u, sample.Npts, 1, VERBOSE)) continue;
+    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);