Index: /branches/czw_branch/20160809/Ohana/src/relastro/src/UpdateObjects.c
===================================================================
--- /branches/czw_branch/20160809/Ohana/src/relastro/src/UpdateObjects.c	(revision 39782)
+++ /branches/czw_branch/20160809/Ohana/src/relastro/src/UpdateObjects.c	(revision 39783)
@@ -5,8 +5,8 @@
 # define PAR_MIN_NPTS 7
 # define PAR_MIN_NPTS_BOOT 6
-# define PM_MIN_NPTS 4
-# define PM_MIN_NPTS_BOOT 3
-# define POS_MIN_NPTS 3
-# define POS_MIN_NPTS_BOOT 2
+# define PM_MIN_NPTS 5
+# define PM_MIN_NPTS_BOOT 4
+# define POS_MIN_NPTS 4
+# define POS_MIN_NPTS_BOOT 4
 
 typedef enum {
@@ -183,29 +183,25 @@
     if (Trange < PM_DT_MIN) {
       // not enough baseline for proper motion, only set mean position
-      mode = FIT_AVERAGE;
       goto justPosition;
     }
     if (parRange < PAR_FACTOR_MIN) {
       // not enough parallax factor range, skip parallax
-      mode = FIT_PM_ONLY;
       goto skipParallax;
     }
     if (fitStats->Npoints < PAR_MIN_NPTS) {
       // not enough data, skip parallax
-      mode = FIT_PM_ONLY;
       goto skipParallax;
     }
 
     // we are going to use the IRLS analysis to calculate the mean solution and the masking
-    // then run N_BOOTSTRAP_SAMPLES to measure the errors
-    fitStats->fitdataPar->getError = !N_BOOTSTRAP_SAMPLES;
+    // then run N_BOOTSTRAP_SAMPLES to measure the errors.  We first measure the OLS error,
+    // and choose the max of the OLS and bootstrap errors
+    fitStats->fitdataPar->getError = TRUE;
     if (USE_IRLS) {
       if (!FitPMandPar_IRLS (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints)) {
-	mode = FIT_PM_ONLY;
 	goto skipParallax;
       }
     } else {
       if (!FitPMandPar_Basic (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints)) {
-	mode = FIT_PM_ONLY;
 	goto skipParallax;
       }
@@ -217,9 +213,9 @@
       fitStats->Nfit = 0;
       int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
-      if (Nnomask < PAR_MIN_NPTS_BOOT) {
-	// if we do not have enough points to assess parallax error, we cannot do the bootstrap analysis.
-	mode = FIT_PM_ONLY;
-	goto skipParallax;
-      }
+
+      // if we do not have enough points to assess parallax error, skip the bootstrap analysis.
+      // (I think the test above means we never do this skip)
+      if (Nnomask < PAR_MIN_NPTS_BOOT) goto skipParallaxBootstrap;
+
       for (k = 0; k < fitStats->NfitAlloc; k++) {
 	BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
@@ -227,4 +223,5 @@
 	fitStats->Nfit ++;
       }
+
       // These calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above)
       // this call expects the fitted parameters to have the formal error set: it will apply the 
@@ -237,4 +234,5 @@
     }
 
+  skipParallaxBootstrap:
     // project Ro, Do back to RA,DEC
     XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &fitStats->coords);
@@ -255,7 +253,5 @@
     valid = valid && (fabs(fitPar.uD) < 4.0);
     valid = valid && (fabs(fitPar.p) < 2.0);
-    if (!valid) {
-      mode = FIT_PM_ONLY;
-    } else {
+    if (valid) {
       average[0].flags |= ID_OBJ_USE_PAR;
     }
@@ -267,57 +263,55 @@
   if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) {
     if (Trange < PM_DT_MIN) {
-      mode = FIT_AVERAGE;
       goto justPosition;
     }
     if (fitStats->Npoints < PM_MIN_NPTS) {
-      mode = FIT_AVERAGE;
       goto justPosition;
     }
 
     // if we have fitted (and accepted) a parallax model, get the best pm fit and chisq
-    // given the set of points (mask is respected)
-    if (average[0].flags & ID_OBJ_USE_PAR) {
+    // given the set of points (mask is respected).  Alternatively, if we do not request
+    // IRLS fitting, the just use OLS fitting (skip bootstrap)
+    if ((average[0].flags & ID_OBJ_USE_PAR) || !USE_IRLS) {
       if (!FitPM_Basic (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
 	average[0].flags |= ID_OBJ_BAD_PM;
 	goto justPosition;
       }
-    } else {
-      fitStats->fitdataPM->getError = !N_BOOTSTRAP_SAMPLES;
-      if (USE_IRLS) {
-	if (!FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
-	  mode = FIT_AVERAGE;
-	  goto justPosition;
-	}
-      } else {
-	if (!FitPM_Basic (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
-	  mode = FIT_AVERAGE;
-	  goto justPosition;
-	}
-      }
+      goto skipProperMotionBootstrap;
+    } 
+
+    // we are going to use the IRLS analysis to calculate the mean solution and the masking
+    // then run N_BOOTSTRAP_SAMPLES to measure the errors.  We first measure the OLS error,
+    // and choose the max of the OLS and bootstrap errors
+    fitStats->fitdataPM->getError = TRUE;
+    if (!FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
+      goto justPosition;
+    }
 
     // in the fits above, we have saved the formal error for the unmasked points.  
     // if we do not have enough points for bootstrap, we will keep those errors
-      if (N_BOOTSTRAP_SAMPLES && (fitPM.Nfit >= PM_MIN_NPTS_BOOT)) {
-	fitStats->Nfit = 0;
-	int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
-	if (Nnomask < PM_MIN_NPTS_BOOT) {
-	  mode = FIT_AVERAGE;
-	  goto justPosition;
-	}
-	for (k = 0; k < fitStats->NfitAlloc; k++) {
-	  BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
-	  if (!FitPM_Basic (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, Nnomask)) continue;
-	  fitStats->Nfit ++;
-	}
+    if (N_BOOTSTRAP_SAMPLES && (fitPM.Nfit >= PM_MIN_NPTS_BOOT)) {
+      fitStats->Nfit = 0;
+      int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
+	
+      // if we do not have enough points to assess p.m. error, skip the bootstrap analysis.
+      // (I think the test above means we never do this skip)
+      if (Nnomask < PM_MIN_NPTS_BOOT) goto skipProperMotionBootstrap;
+
+      for (k = 0; k < fitStats->NfitAlloc; k++) {
+	BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
+	if (!FitPM_Basic (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, Nnomask)) continue;
+	fitStats->Nfit ++;
+      }
+
       // These calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above)
       // this call expects the fitted parameters to have the formal error set: it will apply the 
       // max of the formal and bootstrap errors
-	BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
-	BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
-	BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
-	BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
-      }
-    }
-
+      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
+      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
+      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
+      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
+    }
+
+skipProperMotionBootstrap:
     // project Ro, Do back to RA,DEC
     XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &fitStats->coords);
@@ -326,4 +320,5 @@
     fitStats->Npm ++;
 
+    // XXX a hard-wired hack...
     // unless there is a clear problems (below) with the proper-motion fit or we have a parallax fit, we will use pm fit
     int valid = TRUE;
@@ -335,5 +330,4 @@
     valid = valid && (fabs(fitPM.uD) < 4.0);
     if (!valid) {
-      mode = FIT_AVERAGE;
       average[0].flags |= ID_OBJ_BAD_PM;
     } else {
@@ -349,69 +343,56 @@
     // if we only have one point, this is silly...
 
+    // set the proper motion (to the galaxy model or average value, if desired; else to 0,0)
     FitAstromResultSetPM (&fitPos, 1, average);
-    // fprintf (stderr, "fit 1: %f %f : %f %f\n", fitPos.Ro, fitPos.Do, fitPos.uR, fitPos.uD);
-    if (average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_PM)) {
+
+    // if we already have a valid fit (pm or par), use OLS to fit the position
+    // alternatively, if we do not request IRLS, use OLS
+    // alternatively, if we do not have enough points, use OLS
+    if ((average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_PM)) || (fitStats->Npoints < POS_MIN_NPTS) || !USE_IRLS) {
       if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { 
-	// if this fails, stick with the PM and/or PAR fit from above
+	// if this fails, stick with the PM and/or PAR fit from above, or use a single value
 	goto doneWithFit;
       }
-    } else {
-      if (fitStats->Npoints < POS_MIN_NPTS) {
-	// I will not try to outlier-reject, just calculate the weighted average
-	if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { 
-	  // if we have tried this, we have not masked any points; we will find a single unmasked point below
-	  goto doneWithFit;
-	}
+      // if we have not already gotten a good fit, use this fit
+      if (!(average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_PM))) {
 	average[0].flags |= ID_OBJ_USE_AVE;
 	average[0].flags |= ID_OBJ_RAW_AVE;
-	goto useBasic;
-      }
-      if (USE_IRLS) {
-	if (!FitPosPMfixed_IRLS (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
-	  // if the above fails, we need to clear the masks and try again below
-	  FitPointsClearMasks (fitStats->points, fitStats->Npoints); 
-	  // just calculate the weighted average
-	  if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { 
-	    // if this fails, find a single unmasked point below
-	    goto doneWithFit;
-	  }
-	  average[0].flags |= ID_OBJ_USE_AVE;
-	  average[0].flags |= ID_OBJ_RAW_AVE;
-	  goto useBasic;
-	}
-      } else {
-	if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { 
-	  // if this fails, find a single unmasked point below
-	  goto doneWithFit;
-	}
-	average[0].flags |= ID_OBJ_USE_AVE;
-	average[0].flags |= ID_OBJ_RAW_AVE;
-	goto useBasic;
-      }
-      // fprintf (stderr, "fit 2: %f %f : %f %f\n", fitPos.Ro, fitPos.Do, fitPos.uR, fitPos.uD);
-      if (N_BOOTSTRAP_SAMPLES) {
-	fitStats->Nfit = 0;
-	int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
-	if (Nnomask < POS_MIN_NPTS_BOOT) {
-	  // XXX if the above fails, accept the formal error from unmasked points?
-	  FitPointsClearMasks (fitStats->points, fitStats->Npoints); 
-	  if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { 
-	    // if this fails, find a single unmasked point below
-	    goto doneWithFit;
-	  }
-	  average[0].flags |= ID_OBJ_USE_AVE;
-	  average[0].flags |= ID_OBJ_RAW_AVE;
-	  goto useBasic;
-	}
-	FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average);
-	for (k = 0; k < fitStats->NfitAlloc; k++) {
-	  BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
-	  if (!FitPosPMfixed_Basic (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, Nnomask)) continue;
-	  fitStats->Nfit ++;
-	}
-	BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
-	BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
+      }
+      goto useBasic;
+    } 
+    
+    // try the IRLS fitting, otherwise give up and use OLS
+    if (!FitPosPMfixed_IRLS (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
+      // if the above fails, we need to clear the masks and try again below
+      FitPointsClearMasks (fitStats->points, fitStats->Npoints); 
+      // just calculate the weighted average
+      if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { 
+	// if this fails, find a single unmasked point below
+	goto doneWithFit;
       }
       average[0].flags |= ID_OBJ_USE_AVE;
+      average[0].flags |= ID_OBJ_RAW_AVE;
+      goto useBasic;
+    }
+    average[0].flags |= ID_OBJ_USE_AVE;
+
+    if (N_BOOTSTRAP_SAMPLES && (fitPos.Nfit >= POS_MIN_NPTS_BOOT)) {
+      fitStats->Nfit = 0;
+      int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
+      
+      if (Nnomask < POS_MIN_NPTS_BOOT) goto useBasic;
+      
+      FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average);
+      for (k = 0; k < fitStats->NfitAlloc; k++) {
+	BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
+	if (!FitPosPMfixed_Basic (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, Nnomask)) continue;
+	fitStats->Nfit ++;
+      }
+      
+      // These calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above)
+      // this call expects the fitted parameters to have the formal error set: it will apply the 
+      // max of the formal and bootstrap errors
+      BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
+      BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     }
 
@@ -580,4 +561,5 @@
   average[0].dP         = fit.dp; // parallax error in arcsec
 
+
   average[0].ChiSqAve   = fitPos.chisq;
   average[0].ChiSqPM    = fitPM.chisq;
@@ -873,4 +855,6 @@
   int i;
 
+  myAbort ("this should not be called anymore");
+
   // add up the chi square for the fit
   double chisq = 0.0;