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

Timestamp:

Jan 12, 2021, 4:16:47 PM (6 years ago)

Author:

eugene

Message:

update fit1d_irls to include priors on the parameters

File:

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

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

@@ -10 +10 @@
 double weight_cauchy (double x);
 double VectorFractionInterpolate (double *values, float fraction, int Npts);
+
+// using file-globals to carry values into 'fit_least_squares'
+double *priorValue = NULL;
+double *priorStdev = NULL;
 
 int fit1d_irls (int argc, char **argv) {
@@ -73 +77 @@
     UseMean = FALSE;
   }
+  int UsePriors = FALSE;
+  if ((N = get_argument (argc, argv, "-use-priors"))) {
+    remove_argument (N, &argc, argv);
+    UsePriors = TRUE;
+  }
 
   if ((argc != 4) || !dyvec) {
     // require a y error value
-    gprint (GP_ERR, "USAGE: fit1d x y order -dy wt [-quiet/-q] [-wt wtOut] [-yfit yfit] [-mask maskOut]\n");
+    gprint (GP_ERR, "USAGE: fit1d_irls x y order [-dy wt] [-quiet/-q] [-wt wtOut] [-yfit yfit] [-mask maskOut] [-Nboot N] [-max-iterations N] [-weight-threshold f] [-use-median] [-use-priors]\n");
+    gprint (GP_ERR, "  [-dy wt]        : optional error vector (use value of 1 by default)\n");
+    gprint (GP_ERR, "  [-quiet/-q]     : reduce verbosity\n");
+    gprint (GP_ERR, "  [-wt wtOut]     : save the weights in supplied vector\n");
+    gprint (GP_ERR, "  [-yfit yfit]    : save the fitted values in supplied vector\n");
+    gprint (GP_ERR, "  [-mask maskOut] : save the masked / unmasked state (weight below threshold\n");
+    gprint (GP_ERR, "  [-Nboot N]      : number of bootstrap samples\n");
+    gprint (GP_ERR, "  [-max-iterations N] : max number of iterations\n");
+    gprint (GP_ERR, "  [-weight-threshold f] : limit for masked values\n");
+    gprint (GP_ERR, "  [-use-median]   : use median (as opposed to mean) to calculate weight threshold\n");
+    gprint (GP_ERR, "  [-use-priors]   : include priors for parameters (specified as prior value $P0v and stdev $dP0, $P1v, $dP1, etc).  Unspecified or NAN values are ignored.\n");
     return (FALSE);
   }
@@ -100 +119 @@
   order = atof (argv[3]);
   nterm = order + 1;
-  mterm = 2*nterm;
+  mterm = 2*order + 1;
 
   ALLOCATE_PTR (yoff, opihi_flt, xvec[0].Nelements);
@@ -119 +138 @@
   }
 
+  if (UsePriors) {
+    ALLOCATE (priorValue, double, nterm);
+    ALLOCATE (priorStdev, double, nterm);
+    for (i = 0; i < nterm; i++) {
+      int isFound;
+      double tmpValue;
+      snprintf (name, 64, "P%dv", i);
+      tmpValue = get_double_variable (name, &isFound);
+      priorValue[i] = isFound ? tmpValue : NAN;
+      snprintf (name, 64, "dP%d", i);
+      tmpValue = get_double_variable (name, &isFound);
+      priorStdev[i] = isFound ? tmpValue : NAN;
+      if (!isFound) priorValue[i] = NAN; // need both Pnv and dPn 
+    }
+  }
+
   // define initial weights as 1.0
   for (i = 0; i < xvec[0].Nelements; i++) { 
@@ -137 +172 @@
     }
 
+    // save this solution
+    for (i = 0; i < nterm; i++) {
+      b_old[i][0] = b[i][0];
+    }
+
     if (!fit_least_squares (s, b, c, nterm, mterm, xvec[0].elements.Flt, yvec[0].elements.Flt, dyvec[0].elements.Flt, wt, NULL, xvec[0].Nelements)) {
-      // restore
+      // restore the saved solution
       for (i = 0; i < nterm; i++) {
         b[i][0] = b_old[i][0];
       }
       break;
-    }
-
-    // save this solution
-    for (i = 0; i < nterm; i++) {
-      b_old[i][0] = b[i][0];
     }
 
@@ -172 +207 @@
   return TRUE; */
 
+  // save the solution. we re-use b[][] below in the bootstrap resampling
   for (i = 0; i < nterm; i++) {
     b_old[i][0] = b[i][0];
   }
-
 
   // calculate the weight thresholds to mask the bad points:
@@ -307 +342 @@
   free (dysample);
 
+  if (priorValue) {
+    free (priorValue); priorValue = NULL;
+    free (priorStdev); priorStdev = NULL;
+  }
+
   return (TRUE);
 
@@ -356 +396 @@
     for (j = 0; j < nterm; j++) {
       c[i][j] = s[i + j];
+      if (priorValue && (i == j)) {
+        if (!isnan(priorValue[i])) {
+          c[i][i] += 1.0 / SQ(priorStdev[i]);
+          b[i][0] += priorValue[i] / SQ(priorStdev[i]);
+        }
+      }
     }
   }