Changeset 41470

Timestamp:

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

Author:

eugene

Message:

update fit1d_irls to include priors on the parameters

Location:

trunk/Ohana/src/opihi

Files:

• cmd.data/fit1d.c (modified) (2 diffs)
• cmd.data/fit1d_irls.c (modified) (8 diffs)
• cmd.data/test/fit1d_irls.sh (modified) (1 diff)
• test/tests.sh (modified) (1 diff)

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

@@ -5 +5 @@
   double **c, **b, *s, X, Y, dY, dY2;
   double ClipNSigma, mean, sigma, maxsigma;
-  int i, j, nterm, mterm, order, Npt, Nmask;
+  int i, j, Npt, Nmask;
   int N, Weight, Quiet, ClipNiter;
   Vector *xvec, *yvec, *dyvec;
@@ -66 +66 @@
 
   /* nterm is number of polynomial terms, starting at x^0 */
-  order = atof (argv[3]);
-  nterm = order + 1;
-  mterm = 2*nterm;
+  int order = atof (argv[3]);
+  int nterm = order + 1;
+  int mterm = 2*order + 1;
 
   ALLOCATE (yfit, opihi_flt, xvec[0].Nelements);

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]);
+        }
+      }
     }
   }

trunk/Ohana/src/opihi/cmd.data/test/fit1d_irls.sh

@@ -346 +346 @@
  end
 end
+
+macro try.priors
+ create x 0 100
+ set y = 5.2 + 2.3*x
+ set ytru = 5.2 + 2.3*x
+ gaussdev dyoff ytru[] 0.0 5.0
+ set dy = 5.0 + zero(ytru)
+ set yobs = ytru + dyoff
+ yobs[1] = 210
+ $P0v = 5.0
+ $dP0 = 0.5
+
+ fit1d_irls x yobs 1 -dy dy -yfit yfit -max-iterations 0 -use-priors
+# fit1d_irls x yobs 1 -dy dy -yfit yfit -max-iterations 5
+# fit1d_irls  x yobs 1 -dy dy -yfit yfit -max-iterations 0
+
+ lim x yobs; clear; box; 
+ plot -c blue -x line -lw 2 x ytru
+ plot x yobs -pt cir -op 0.5
+ plot -c red -x line -lw 2 x yfit -op 0.5
+end
+ 
+macro run.priors
+ create x 0 100
+ set y = 5.2 + 2.3*x
+ set ytru = 5.2 + 2.3*x
+ gaussdev dyoff ytru[] 0.0 5.0
+ set dy = 5.0 + zero(ytru)
+ set yobs = ytru + dyoff
+
+ $P0v = NAN
+ $P1v = NAN
+
+ $n = 0
+
+ if ($n == 0)
+   $P0v = 4.5
+   $dP0 = 1.0
+ end
+ if ($n == 1)
+   $P1v = 2.3
+   $dP1 = 0.5
+ end
+
+ set yobsS = yobs
+
+ delete -q C0norm C0irls C0prio C0irpi
+ delete -q N0norm N0irls N0prio N0irpi
+ for i 0 5000
+
+   set yobs = yobsS
+   for j 0 10
+     $Npt = int(yobs[]*rnd(0))
+     yobs[$Npt] = 1000*rnd(0) - 500
+   end
+   fit1d_irls -q x yobs 1 -dy dy -yfit yfit -mask ymsk -max-iterations 0
+   concat $C$n C0norm; vstat -q ymsk; concat $TOTAL N0norm
+   fit1d_irls -q x yobs 1 -dy dy -yfit yfit -mask ymsk -max-iterations 0 -use-priors
+   concat $C$n C0prio; vstat -q ymsk; concat $TOTAL N0prio
+   fit1d_irls -q x yobs 1 -dy dy -yfit yfit -mask ymsk -max-iterations 5
+   concat $C$n C0irls; vstat -q ymsk; concat $TOTAL N0irls
+   fit1d_irls -q x yobs 1 -dy dy -yfit yfit -mask ymsk -max-iterations 5 -use-priors
+   concat $C$n C0irpi; vstat -q ymsk; concat $TOTAL N0irpi
+ end
+
+  histogram C0norm nCnorm -30 30 0.25 -range dC
+  histogram C0prio nCprio -30 30 0.25 -range dC
+  histogram C0irls nCirls -30 30 0.25 -range dC
+  histogram C0irpi nCirpi -30 30 0.25 -range dC
+
+  lim dC nCirls; clear; box; 
+
+  plot dC nCnorm -x hist -op 0.6 -lw 3 -c black  ; vstat C0norm 
+  plot dC nCprio -x hist -op 0.6 -lw 3 -c blue   ; vstat C0prio
+  plot dC nCirls -x hist -op 0.6 -lw 3 -c red    ; vstat C0irls
+  plot dC nCirpi -x hist -op 0.6 -lw 3 -c blue60 ; vstat C0irpi
+
+  if ($n == 0)
+    line -c red -lt dot $P0v 0 to $P0v 1000
+    line -c blue -lt dot 5.2 0 to 5.2 1000
+  else
+    line -c red -lt dot $P1v 0 to $P1v 1000
+    line -c blue -lt dot 2.3 0 to 2.3 1000
+  end
+end
+ 
+macro run.2dpriors.fixmeas
+ $C0tru = +3
+ $C1tru = -4
+ $Sigma =  5
+
+ create x -25 25
+ set ytru = $C0tru + $C1tru*x
+ gaussdev yoff ytru[] 0.0 $Sigma
+
+ set dy = $Sigma + zero(ytru)
+ set yobs = ytru + yoff
+
+ $P0v = NAN
+ $P1v = NAN
+
+ $P0v = $C0tru + 0.5
+ $dP0 = 1.0
+
+ $P1v = $C1tru + 0.2
+ $dP1 = 0.5
+
+ set yobsSave = yobs
+
+ delete -q N0norm N0irls N0prio N0irpi
+ delete -q C0norm C0irls C0prio C0irpi
+ delete -q C1norm C1irls C1prio C1irpi
+
+ for i 0 5000
+
+   set yobs = yobsSave
+
+   # inject outliers
+   for j 0 10
+     $Npt = int(yobs[]*rnd(0))
+     yobs[$Npt] = 1000*rnd(0) - 500
+   end
+
+   fit1d_irls -q x yobs 1 -dy dy -yfit yfit -mask ymsk -max-iterations 0
+   concat $C0 C0norm; concat $C1 C1norm; vstat -q ymsk; concat $TOTAL N0norm
+
+   fit1d_irls -q x yobs 1 -dy dy -yfit yfit -mask ymsk -max-iterations 0 -use-priors
+   concat $C0 C0prio; concat $C1 C1prio; vstat -q ymsk; concat $TOTAL N0prio
+
+   fit1d_irls -q x yobs 1 -dy dy -yfit yfit -mask ymsk -max-iterations 5
+   concat $C0 C0irls; concat $C1 C1irls; vstat -q ymsk; concat $TOTAL N0irls
+
+   fit1d_irls -q x yobs 1 -dy dy -yfit yfit -mask ymsk -max-iterations 5 -use-priors
+   concat $C0 C0irpi; concat $C1 C1irpi; vstat -q ymsk; concat $TOTAL N0irpi
+ end
+
+ histogram C0norm nC0norm -30 30 0.25 -range dC
+ histogram C0prio nC0prio -30 30 0.25 -range dC
+ histogram C0irls nC0irls -30 30 0.25 -range dC
+ histogram C0irpi nC0irpi -30 30 0.25 -range dC
+
+ histogram C1norm nC1norm -30 30 0.25 -range dC
+ histogram C1prio nC1prio -30 30 0.25 -range dC
+ histogram C1irls nC1irls -30 30 0.25 -range dC
+ histogram C1irpi nC1irpi -30 30 0.25 -range dC
+
+ dev -n 0
+ lim dC nC0irpi; clear; box; 
+ line -c red  -lw 2 -lt dot $P0v   0 to $P0v   1000
+ line -c blue -lw 2 -lt dot $C0tru 0 to $C0tru 1000
+ plot dC nC0norm -x hist -op 0.6 -lw 3 -c black  ; vstat C0norm 
+ plot dC nC0prio -x hist -op 0.6 -lw 3 -c blue   ; vstat C0prio
+ plot dC nC0irls -x hist -op 0.6 -lw 3 -c red    ; vstat C0irls
+ plot dC nC0irpi -x hist -op 0.6 -lw 3 -c blue60 ; vstat C0irpi
+
+ dev -n 1
+ lim dC nC1irpi; clear; box; 
+ line -c red  -lw 2 -lt dot $P1v   0 to $P1v   1000
+ line -c blue -lw 2 -lt dot $C1tru 0 to $C1tru 1000
+ plot dC nC1norm -x hist -op 0.6 -lw 3 -c black  ; vstat C1norm 
+ plot dC nC1prio -x hist -op 0.6 -lw 3 -c blue   ; vstat C1prio
+ plot dC nC1irls -x hist -op 0.6 -lw 3 -c red    ; vstat C1irls
+ plot dC nC1irpi -x hist -op 0.6 -lw 3 -c blue60 ; vstat C1irpi
+
+end
+ 
+macro run.2dpriors
+ $C0tru = +3
+ $C1tru = -4
+ $Sigma =  5
+
+ create x -25 25
+ set ytru = $C0tru + $C1tru*x
+
+ set dy = $Sigma + zero(ytru)
+
+ $P0v = NAN
+ $P1v = NAN
+
+ $P0v = $C0tru + 0.5
+ $dP0 = 1.0
+
+ $P1v = $C1tru + 0.2
+ $dP1 = 0.5
+
+
+ delete -q N0norm N0irls N0prio N0irpi
+ delete -q C0norm C0irls C0prio C0irpi
+ delete -q C1norm C1irls C1prio C1irpi
+
+ for i 0 5000
+
+   gaussdev yoff ytru[] 0.0 $Sigma
+   set yobs = ytru + yoff
+
+   # inject outliers
+   for j 0 10
+     $Npt = int(yobs[]*rnd(0))
+     yobs[$Npt] = 1000*rnd(0) - 500
+   end
+
+   fit1d_irls -q x yobs 1 -dy dy -yfit yfit -mask ymsk -max-iterations 0
+   concat $C0 C0norm; concat $C1 C1norm; vstat -q ymsk; concat $TOTAL N0norm
+
+   fit1d_irls -q x yobs 1 -dy dy -yfit yfit -mask ymsk -max-iterations 0 -use-priors
+   concat $C0 C0prio; concat $C1 C1prio; vstat -q ymsk; concat $TOTAL N0prio
+
+   fit1d_irls -q x yobs 1 -dy dy -yfit yfit -mask ymsk -max-iterations 5
+   concat $C0 C0irls; concat $C1 C1irls; vstat -q ymsk; concat $TOTAL N0irls
+
+   fit1d_irls -q x yobs 1 -dy dy -yfit yfit -mask ymsk -max-iterations 5 -use-priors
+   concat $C0 C0irpi; concat $C1 C1irpi; vstat -q ymsk; concat $TOTAL N0irpi
+ end
+end
+
+macro plot.2dpriors 
+
+ if (1)
+   histogram C0norm nC0norm -30 30 0.25 -range dC
+   histogram C0prio nC0prio -30 30 0.25 -range dC
+   histogram C0irls nC0irls -30 30 0.25 -range dC
+   histogram C0irpi nC0irpi -30 30 0.25 -range dC
+   
+   histogram C1norm nC1norm -30 30 0.25 -range dC
+   histogram C1prio nC1prio -30 30 0.25 -range dC
+   histogram C1irls nC1irls -30 30 0.25 -range dC
+   histogram C1irpi nC1irpi -30 30 0.25 -range dC
+   
+   dev -n 0
+   lim dC nC0irpi; clear; box; 
+   line -c red  -lw 2 -lt dot $P0v   0 to $P0v   1000
+   line -c blue -lw 2 -lt dot $C0tru 0 to $C0tru 1000
+   plot dC nC0norm -x hist -op 0.6 -lw 3 -c black  ; vstat C0norm 
+   plot dC nC0prio -x hist -op 0.6 -lw 3 -c blue         ; vstat C0prio
+   plot dC nC0irls -x hist -op 0.6 -lw 3 -c red    ; vstat C0irls
+   plot dC nC0irpi -x hist -op 0.6 -lw 3 -c blue60 ; vstat C0irpi
+   
+   dev -n 1
+   lim dC nC1irpi; clear; box; 
+   line -c red  -lw 2 -lt dot $P1v   0 to $P1v   1000
+   line -c blue -lw 2 -lt dot $C1tru 0 to $C1tru 1000
+   plot dC nC1norm -x hist -op 0.6 -lw 3 -c black  ; vstat C1norm 
+   plot dC nC1prio -x hist -op 0.6 -lw 3 -c blue         ; vstat C1prio
+   plot dC nC1irls -x hist -op 0.6 -lw 3 -c red    ; vstat C1irls
+   plot dC nC1irpi -x hist -op 0.6 -lw 3 -c blue60 ; vstat C1irpi
+ end
+
+ dev -n 2
+ $C0range = 50; $C1range = 5
+ resize 1800 1800
+ clear -s
+ section a0 0.0 0.0 0.5 0.5
+ section a1 0.5 0.0 0.5 0.5
+ section a2 0.0 0.5 0.5 0.5
+ section a3 0.5 0.5 0.5 0.5
+
+ section a0; subplot.2dprior C0norm C1norm
+ section a1; subplot.2dprior C0prio C1prio
+ section a2; subplot.2dprior C0irls C1irls
+ section a3; subplot.2dprior C0irpi C1irpi
+
+ dev -n 3
+ $C0range = 8; $C1range = 2
+ resize 1800 1800
+ clear -s
+ section a0 0.0 0.0 0.5 0.5
+ section a1 0.5 0.0 0.5 0.5
+ section a2 0.0 0.5 0.5 0.5
+ section a3 0.5 0.5 0.5 0.5
+
+ section a0; subplot.2dprior C0norm C1norm
+ section a1; subplot.2dprior C0prio C1prio
+ section a2; subplot.2dprior C0irls C1irls
+ section a3; subplot.2dprior C0irpi C1irpi
+end
+
+macro subplot.2dprior 
+ if ($0 != 3)
+   echo "USAGE: subplot.2dprior (vec) (vec)"
+   break
+ end
+    
+ lim {$C0tru - $C0range} {$C0tru + $C0range} {$C1tru - $C1range} {$C1tru + $C1range}
+ box; label -x C0 -y C1 
+
+ plot -op 0.1 -pt cir -sz 1 -c blue $1 $2
+
+ line -c blue -lw 2 -lt dot $C0tru -1000 to $C0tru 1000
+ line -c red  -lw 2 -lt dot $P0v   -1000 to $P0v   1000
+ line -c blue -lw 2 -lt dot -1000 $C1tru to 1000 $C1tru
+ line -c red  -lw 2 -lt dot -1000   $P1v to 1000 $P1v
+end

trunk/Ohana/src/opihi/test/tests.sh

@@ -89 +89 @@
     else
       $Npass ++
+      if ($VERBOSE > 2)
+        echo " $tests:$Tk : success"
+      end
     end
     $Ntest ++