Changeset 39610 for trunk/Ohana/src/opihi

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

Location:

trunk/Ohana/src/opihi

Files:

• cmd.astro/astrom_ops.c (modified) (2 diffs)
• cmd.astro/fitplx.c (modified) (1 diff)
• cmd.astro/fitplx_irls.c (modified) (9 diffs)
• include/astro.h (modified) (2 diffs)

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

@@ -54 +54 @@
   ALLOCATE (data->Wx, double, N);
   ALLOCATE (data->Wy, double, N);
+  ALLOCATE (data->Qx, double, N);
+  ALLOCATE (data->Qy, double, N);
+  ALLOCATE (data->qx, double, N);
+  ALLOCATE (data->qy, double, N);
   ALLOCATE (data->index, int, N);
   return TRUE;
@@ -68 +72 @@
   FREE (data->Wx);
   FREE (data->Wy);
+  FREE (data->Qx);
+  FREE (data->Qy);
+  FREE (data->qx);
+  FREE (data->qy);
   FREE (data->index);
 }

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

@@ -597 +597 @@
   dsort (values, Npts);
 
-  if (Npts % 2) {
-    int Ncenter = Npts / 2;
-    *median = values[Ncenter];
-  } else {
-    int Ncenter = Npts / 2 - 1;
-    *median = 0.5*(values[Ncenter] + values[Ncenter + 1]);
-  }
-
-  double Slo = VectorFractionInterpolate (values, 0.158655, Npts);
-  double Shi = VectorFractionInterpolate (values, 0.841345, Npts);
-
-  *sigma = (Shi - Slo) / 2.0;
-
+  if (median) {
+    if (Npts % 2) {
+      int Ncenter = Npts / 2;
+      *median = values[Ncenter];
+    } else {
+      int Ncenter = Npts / 2 - 1;
+      *median = 0.5*(values[Ncenter] + values[Ncenter + 1]);
+    }
+  }
+    
+  if (sigma) {
+    double Slo = VectorFractionInterpolate (values, 0.158655, Npts);
+    double Shi = VectorFractionInterpolate (values, 0.841345, Npts);
+    *sigma = (Shi - Slo) / 2.0;
+  }
+  
   return TRUE;
 }

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);

trunk/Ohana/src/opihi/include/astro.h

@@ -20 +20 @@
   double *dX;
   double *dY;
-  double *Wx;
-  double *Wy;
+  double *Wx; // modified weight based on distance from fit
+  double *Wy; // modified weight based on distance from fit
+  double *Qx; // raw error-based weight (1/dX^2)
+  double *Qy;
+  double *qx; // modified error-based weight (1/dX^2) * Wx
+  double *qy;
   int *index;
   int Npts;
@@ -59 +63 @@
 
 int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, double binning_step, int VERBOSE);
-int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, double *Wx, double *Wy, int Npts, int max_iterations, double outlier_limit, double binning_step, int VERBOSE);
+int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, 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);
 
 int weighted_LS_PLX (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts, double **A, double **B, int VERBOSE);