Changeset 39610 for trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c

Timestamp:

Jun 24, 2016, 5:07:24 PM (10 years ago)

Author:

eugene

Message:

modify irls fitting to apply the multiply the modified weight by the error-based weight

File:

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

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

@@ -119 +119 @@
                          fitdata.t, fitdata.pX, fitdata.pY,
                          fitdata.Wx, fitdata.Wy,
+                         fitdata.Qx, fitdata.Qy,
+                         fitdata.qx, fitdata.qy,
                          fitdata.Npts, max_iterations, outlier_limit, binning_step, VERBOSE)) {
     return FALSE;
@@ -204 +206 @@
     }
 
-    // now calculate median and sigma for each vector
-    VectorRobustStats (pvec,  &fit.p,  &fit.dp);
-    VectorRobustStats (uRvec, &fit.uR, &fit.duR);
-    VectorRobustStats (uDvec, &fit.uD, &fit.duD);
-    VectorRobustStats (Rvec,  &fit.Ro, &fit.dRo);
-    VectorRobustStats (Dvec,  &fit.Do, &fit.dDo);
+    // now calculate robust sigma for each vector
+    VectorRobustStats (pvec,  NULL, &fit.dp);
+    VectorRobustStats (uRvec, NULL, &fit.duR);
+    VectorRobustStats (uDvec, NULL, &fit.duD);
+    VectorRobustStats (Rvec,  NULL, &fit.dRo);
+    VectorRobustStats (Dvec,  NULL, &fit.dDo);
+
+    PlxFitDataFree (&sample);    
   }
 
@@ -250 +254 @@
   set_variable ("Nfit",  fit.Nfit);
 
-
+  PlxFitDataFree (&fitdata);
   
   return (TRUE);
@@ -257 +261 @@
 /* do we want an init function which does the alloc and a clear function to free? */
 int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD,
-                      double *Wx, double *Wy,
+                      double *Wx, double *Wy, double *Qx, double *Qy, double *qx, double *qy,
                       int Npts, int max_iterations, double outlier_limit, double binning_step, int VERBOSE) {
 
@@ -270 +274 @@
 
   double sigma_ols, sigma_hat; // Sigma estimates.
-  //  double *Wx, *Wy;             // Weight vectors.  Not errors.
   double *rx, *ry;             // Deviation from model
   double *u;                   // Deviation magnitude
@@ -305 +308 @@
   // Convert the measurement errors into the initial weights.
   for (i = 0; i < Npts; i++) {
-    Wx[i] = (fabs(dX[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dX[i]);
-    Wy[i] = (fabs(dY[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dY[i]);
+    Qx[i] = (fabs(dX[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dX[i]);
+    Qy[i] = (fabs(dY[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dY[i]);
   }
   
@@ -313 +316 @@
   // Solve OLS equation
   if (binning_step == 0.0) {
-    if (!weighted_LS_PLX(T,pR,pD,X,Wx,Y,Wy,Npts,
-                         A,B,VERBOSE)) {
+    if (!weighted_LS_PLX(T,pR,pD,X,Qx,Y,Qy,Npts,A,B,VERBOSE)) {
       // Handle fail case
       return(FALSE);
@@ -329 +331 @@
     double *WYbin = NULL;
 
-    if (!bin_points_PLX(T, pR, pD, X, Wx, Y, Wy, Npts,
+    if (!bin_points_PLX(T, pR, pD, X, Qx, Y, Qy, Npts,
                         &Tbin, &pRbin, &pDbin, &Xbin, &WXbin, &Ybin, &WYbin, &Nbins,  binning_step) ) {
       return(FALSE);
@@ -396 +398 @@
       Wx[i] = weight_cauchy(rx[i] / dX[i]);
       Wy[i] = weight_cauchy(ry[i] / dY[i]);
+      qx[i] = Qx[i] * Wx[i];
+      qy[i] = Qy[i] * Wy[i];
     }
 
     // Solve
-    if (!weighted_LS_PLX(T,pR,pD,X,Wx,Y,Wy,Npts,A,B,VERBOSE)) {
+    if (!weighted_LS_PLX(T,pR,pD,X,qx,Y,qy,Npts,A,B,VERBOSE)) {
       // Handle fail case
       return(FALSE);