Changeset 39731 for branches/czw_branch/20160809/Ohana/src/relastro/src/FitPMandPar.c

Timestamp:

Oct 6, 2016, 9:44:08 AM (10 years ago)

Author:

eugene

Message:

split galaxy motion model option into two options: galaxy motion model + apply proper motion; remove deprecated chunk of code; add options to skip HSC, CFH, PS1_CHIP, or PS1_STACK images (leaving them static); calculate the error based on the ordinary least squares and regular poisson weights (when doing bootstrap error analysis, choose the MAX); calculate dXpixSys,dYpixSys in arcsec based on plate scale; add option to avoid changing measurements for specific cameras; StarMap tests are now done in real arcsec using plate scales; option to fit the chips without modifying the mean objects; remove excessively verbose memory checks; add tests for HSC and CFH photcodes

File:

• branches/czw_branch/20160809/Ohana/src/relastro/src/FitPMandPar.c (modified) (5 diffs)

branches/czw_branch/20160809/Ohana/src/relastro/src/FitPMandPar.c

@@ -25 +25 @@
   }
 
+  fit->useWeight = FALSE; // Ordinary Least Squares
   if (!FitPMandPar_MinChisq (fit, data, points, Npoints)) return FALSE;
   if (!FitPMandPar_SetChisq (fit, data, points, Npoints)) return FALSE;
@@ -47 +48 @@
   
   // Solve OLS equation: failure here means the chisq matrix is degenerate, give up entirely
+  fit->useWeight = FALSE; // Ordinary Least Squares
   if (!FitPMandPar_MinChisq(fit, data, points, Npoints)) {
     return(FALSE);
@@ -90 +92 @@
 
     // Solve with the new weights
+    fit->useWeight = TRUE; // Reweighted Least Squares
     if (!FitPMandPar_MinChisq(fit, data, points, Npoints)) {
 
@@ -157 +160 @@
   }
 
-  // this section calculates the formal error on the weighted fit using the covariance values
-  // NOTE EAM: in tests (fitpm.c), they seem to be too large by a factor of ~5.37
+  // this section calculates the formal error on the regular (unweighted) fit using the covariance values
+  // NOTE 20160929 : use only the unmasked points to calculate the error
   if (data->getError) {
-    double ax = 0.0, ay = 0.0;
-    double bx = 0.0, by = 0.0;
-
-    for (i = 0; i < Npoints; i++) {
-      ax += dpsi_cauchy(points[i].rx / points[i].dX);
-      ay += dpsi_cauchy(points[i].ry / points[i].dY);
-
-      bx += SQ(points[i].Wx);
-      by += SQ(points[i].Wy);
-    }
-    ax /= 1.0 * Npoints;  // mean(psi_dot(r))
-    ay /= 1.0 * Npoints; 
-    bx /= 1.0 * (Npoints - data->Nterms); // mean(psi^2(r)) * (N / (N-p))
-    by /= 1.0 * (Npoints - data->Nterms);
-  
-    double lambda_x = 1.0 + (data->Nterms / Npoints) * (1 - ax) / ax;
-    double lambda_y = 1.0 + (data->Nterms / Npoints) * (1 - ay) / ay;
-  
-    double sigma_robust_x = lambda_x * sqrt(bx) * sigma_hat * 2.385 / ax;
-    double sigma_robust_y = lambda_y * sqrt(by) * sigma_hat * 2.385 / ay;
-
-    // This is actually sigma^2, as that's the factor in the covariance (dumouchel 4.1)
-    double sigma_final_x  = MAX(SQ(sigma_robust_x), (2 * Npoints * SQ(sigma_robust_x) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms)));
-    double sigma_final_y  = MAX(SQ(sigma_robust_y), (2 * Npoints * SQ(sigma_robust_y) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms)));
-
-    fit[0].dRo = sqrt(data->Cov[0][0]);
-    fit[0].duR = sqrt(data->Cov[1][1]);
-    fit[0].dDo = sqrt(data->Cov[2][2]);
-    fit[0].duD = sqrt(data->Cov[3][3]);
-    fit[0].dp  = sqrt(data->Cov[4][4]);
-
-    fit[0].dRo *= sigma_final_x;
-    fit[0].duR *= sigma_final_x;
-    fit[0].dDo *= sigma_final_y;
-    fit[0].duD *= sigma_final_y;
-    fit[0].dp  *= sqrt(sigma_final_x * sigma_final_y);
+    FitAstromResult fitErrors;
+    FitAstromResultInit (&fitErrors);
+    fitErrors.useWeight = FALSE; 
+
+    FitPMandPar_MinChisq(&fitErrors, data, points, Npoints);
+
+    // we use the errors from a simple OLS, ignoring masked points
+    fit[0].dRo = fitErrors.dRo;
+    fit[0].duR = fitErrors.duR;
+    fit[0].dDo = fitErrors.dDo;
+    fit[0].duD = fitErrors.duD;
+    fit[0].dp  = fitErrors.dp;
   }
 
@@ -219 +198 @@
     Nfit ++;
 
-    wx = points[i].qx;
-    wy = points[i].qy;
+    if (fit->useWeight) {
+      wx = points[i].qx;
+      wy = points[i].qy;
+    } else {
+      wx = points[i].Qx;
+      wy = points[i].Qy;
+    }
 
     Wx += wx;