Changeset 39640 for trunk/Ohana/src/opihi/cmd.data/fit1d_irls.c

Timestamp:

Jul 28, 2016, 2:14:52 PM (10 years ago)

Author:

eugene

Message:

new api to set vector value en-masse; fix the fit1d_irls mean calculation

File:

• trunk/Ohana/src/opihi/cmd.data/fit1d_irls.c (modified) (8 diffs)

trunk/Ohana/src/opihi/cmd.data/fit1d_irls.c

@@ -2 +2 @@
 
 // These should probably be tunable:
-# define MAX_ITERATIONS 10
 # define FIT_TOLERANCE 1e-4
 # define FLT_TOLERANCE 1e-6
 # define WEIGHT_THRESHOLD 0.3
-# define NBOOTSTRAP 100
 
 int fit_least_squares (double *s, double **b, double **c, int nterm, int mterm, opihi_flt *x, opihi_flt *y, opihi_flt *dy, opihi_flt *wt, char *mask, int Npts);
@@ -32 +30 @@
   }
 
+  Vector *wtVec = NULL;
+  if ((N = get_argument (argc, argv, "-wt"))) {
+    remove_argument (N, &argc, argv);
+    if (!(wtVec = SelectVector (argv[N], ANYVECTOR, TRUE))) return FALSE;
+    remove_argument (N, &argc, argv);
+  }
+
+  Vector *yfitVec = NULL;
+  if ((N = get_argument (argc, argv, "-yfit"))) {
+    remove_argument (N, &argc, argv);
+    if (!(yfitVec = SelectVector (argv[N], ANYVECTOR, TRUE))) return FALSE;
+    remove_argument (N, &argc, argv);
+  }
+  Vector *mkVec = NULL;
+  if ((N = get_argument (argc, argv, "-mask"))) {
+    remove_argument (N, &argc, argv);
+    if (!(mkVec = SelectVector (argv[N], ANYVECTOR, TRUE))) return FALSE;
+    remove_argument (N, &argc, argv);
+  }
+
+  int NBOOTSTRAP = 100;
+  if ((N = get_argument (argc, argv, "-Nboot"))) {
+    remove_argument (N, &argc, argv);
+    NBOOTSTRAP = atoi(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  int MaxIterations = 10;
+  if ((N = get_argument (argc, argv, "-max-iterations"))) {
+    remove_argument (N, &argc, argv);
+    MaxIterations = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  float WeightThreshold = 0.3;
+  if ((N = get_argument (argc, argv, "-weight-threshold"))) {
+    remove_argument (N, &argc, argv);
+    WeightThreshold = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  int UseMean = TRUE;
+  if ((N = get_argument (argc, argv, "-use-median"))) {
+    remove_argument (N, &argc, argv);
+    UseMean = FALSE;
+  }
+
   if ((argc != 4) || !dyvec) {
     // require a y error value
-    gprint (GP_ERR, "USAGE: fit1d x y order -dy wt [-quiet/-q] [-clip Nsigma Niter]\n");
+    gprint (GP_ERR, "USAGE: fit1d x y order -dy wt [-quiet/-q] [-wt wtOut] [-yfit yfit] [-mask maskOut]\n");
     return (FALSE);
   }
@@ -53 +95 @@
   }
 
+  if (mkVec) ResetVector (mkVec, OPIHI_INT, xvec[0].Nelements);
+
   /* nterm is number of polynomial terms, starting at x^0 */
   order = atof (argv[3]);
@@ -87 +131 @@
 
   int converged = FALSE;
-  for (int iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations++) {
+  for (int iterations = 0; !converged && (iterations < MaxIterations); iterations++) {
 
     for (i = 0; i < xvec[0].Nelements; i++) {
@@ -116 +160 @@
   }
       
+  /* XXX TEST 
+  for (i = 0; i < nterm; i++) {
+    snprintf (name, 64, "C%d", i);
+    set_variable (name, b[i][0]);
+    if (!Quiet) gprint (GP_ERR, "%f x^%d ", b[i][0], i);
+    snprintf (name, 64, "dC%d", i);
+    set_variable (name, sqrt(c[i][i]));
+    if (!Quiet) gprint (GP_ERR, "%f x^%d ", sqrt(c[i][i]), i);
+  }
+  if (!Quiet) gprint (GP_ERR, "\n");
+  return TRUE; */
+
+  for (i = 0; i < nterm; i++) {
+    b_old[i][0] = b[i][0];
+  }
+
+
   // calculate the weight thresholds to mask the bad points:
   double Sum_W = 0.0;
+  ALLOCATE_PTR (wtlist, double, xvec[0].Nelements);
   for (i = 0; i < xvec[0].Nelements; i++) {
     wt[i] = weight_cauchy (yoff[i] / dy[i]);
+    wtlist[i] = wt[i];
     Sum_W += wt[i];
   }
-  double WtThreshold = WEIGHT_THRESHOLD * Sum_W / (1.0 * xvec[0].Nelements);
+  dsort (wtlist, xvec[0].Nelements);
+  int midpt = 0.5 * xvec[0].Nelements;
+  double WtMedian = (xvec[0].Nelements % 2) ? wtlist[midpt] : 0.5*(wtlist[midpt] + wtlist[midpt-1]);
+  double WtThreshold = UseMean ? WeightThreshold * Sum_W / (1.0 * xvec[0].Nelements) : WeightThreshold * WtMedian;
+  free (wtlist);
 
   // generate the unmasked subset
@@ -132 +199 @@
   int Npoints = 0;
   for (i = 0; i < xvec[0].Nelements; i++) {
-    if (wt[i] < WtThreshold) continue;
+    if (mkVec) mkVec->elements.Int[i] = 0;
+    if (wt[i] < WtThreshold) {
+      if (mkVec) mkVec->elements.Int[i] = 1;
+      continue;
+    }
      xunmask[Npoints] =  xvec[0].elements.Flt[i];
      yunmask[Npoints] =  yvec[0].elements.Flt[i];
@@ -184 +255 @@
   for (i = 0; i < nterm; i++) {
     snprintf (name, 64, "C%d", i);
-    set_variable (name, b[i][0]);
-    if (!Quiet) gprint (GP_ERR, "%f x^%d ", b[i][0], i);
+    set_variable (name, b_old[i][0]);
+    if (!Quiet) gprint (GP_ERR, "%f x^%d ", b_old[i][0], i);
   }
   if (!Quiet) gprint (GP_ERR, "\n");
@@ -216 +287 @@
   free (s);
   free (yoff);
-  free (yfit);
-  free (wt);
+
+  if (wtVec) {
+    SetVectorValues (wtVec, wt, OPIHI_FLT, xvec[0].Nelements);
+  } else {
+    free (wt);
+  }
+  if (yfitVec) {
+    SetVectorValues (yfitVec, yfit, OPIHI_FLT, xvec[0].Nelements);
+  } else {
+    free (yfit);
+  }
   
   free (xunmask);