Index: trunk/Ohana/src/relastro/Makefile =================================================================== --- trunk/Ohana/src/relastro/Makefile (revision 39245) +++ trunk/Ohana/src/relastro/Makefile (revision 39246) @@ -33,5 +33,7 @@ $(SRC)/FitPM_IRLS.$(ARCH).o \ $(SRC)/FitPMandPar.$(ARCH).o \ +$(SRC)/FitPMandPar_IRLS.$(ARCH).o \ $(SRC)/FitPosPMfixed.$(ARCH).o \ +$(SRC)/BootstrapOps.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o \ $(SRC)/MosaicOps.$(ARCH).o \ @@ -108,5 +110,7 @@ $(SRC)/FitPM_IRLS.$(ARCH).o \ $(SRC)/FitPMandPar.$(ARCH).o \ +$(SRC)/FitPMandPar_IRLS.$(ARCH).o \ $(SRC)/FitPosPMfixed.$(ARCH).o \ +$(SRC)/BootstrapOps.$(ARCH).o \ $(SRC)/ImageOps.$(ARCH).o \ $(SRC)/MosaicOps.$(ARCH).o \ @@ -172,4 +176,5 @@ $(SRC)/FitPMandPar.$(ARCH).o \ $(SRC)/FitPMandPar_IRLS.$(ARCH).o \ +$(SRC)/BootstrapOps.$(ARCH).o \ $(SRC)/mkpolyterm.$(ARCH).o \ $(SRC)/fitpoly.$(ARCH).o Index: trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- trunk/Ohana/src/relastro/include/relastro.h (revision 39245) +++ trunk/Ohana/src/relastro/include/relastro.h (revision 39246) @@ -138,4 +138,6 @@ double *Beta_prev; int Nterms; + int getChisq; + int getError; } FitAstromData; @@ -173,4 +175,5 @@ FitAstromPoint *points; FitAstromPoint *sample; + FitAstromPoint *nomask; int Npoints; int NpointsAlloc; @@ -732,4 +735,6 @@ int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode); int BootstrapResample (FitAstromPoint *sample, FitAstromPoint *points, int Npoints); +int BootstrapSaveUnmasked (FitAstromPoint *nomask, FitAstromPoint *points, int Npoints); + int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog); int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange); Index: trunk/Ohana/src/relastro/src/FitAstromOps.c =================================================================== --- trunk/Ohana/src/relastro/src/FitAstromOps.c (revision 39245) +++ trunk/Ohana/src/relastro/src/FitAstromOps.c (revision 39246) @@ -87,4 +87,5 @@ fitStats->points = NULL; fitStats->sample = NULL; + fitStats->nomask = NULL; fitStats->Npoints = 0; fitStats->NpointsAlloc = Nmax; @@ -92,4 +93,5 @@ ALLOCATE (fitStats->points, FitAstromPoint, Nmax); ALLOCATE (fitStats->sample, FitAstromPoint, Nmax); + ALLOCATE (fitStats->nomask, FitAstromPoint, Nmax); } @@ -139,4 +141,5 @@ FREE (fitStats->points); FREE (fitStats->sample); + FREE (fitStats->nomask); FitAstromDataFree (fitStats->fitdataPos); @@ -160,4 +163,7 @@ ALLOCATE (fit->Beta_prev, double, Nterms); fit->Nterms = Nterms; + + fit->getChisq = TRUE; + fit->getError = TRUE; return fit; @@ -203,5 +209,5 @@ object->u = 0.0; - object->mask = 0; + object->mask = 0; // keep point if mask == 0 return; } Index: trunk/Ohana/src/relastro/src/FitPM_IRLS.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPM_IRLS.c (revision 39245) +++ trunk/Ohana/src/relastro/src/FitPM_IRLS.c (revision 39246) @@ -11,15 +11,8 @@ int i,j; - int p = 4; - int n = 2 * Npoints; - int dof = n - p; - - // data->A,B,Cov,Beta,Beta_prev are allocated outside by FitAstromDataInit() - // points->Wx,Wy,rx,ry,u are elements of FitAstromPoint + int Ndof = 2 * Npoints - data->Nterms; // Convert the measurement errors into initial weights. for (i = 0; i < Npoints; i++) { - // points[i].Wx = 1 / points[i].dX; - // points[i].Wy = 1 / points[i].dY; points[i].Wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1 / SQ(points[i].dX); points[i].Wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1 / SQ(points[i].dY); @@ -39,5 +32,5 @@ sigma_ols += SQ(points[i].rx) + SQ(points[i].ry); } - sigma_ols = sqrt(sigma_ols / dof); + sigma_ols = sqrt(sigma_ols / (float)Ndof); // Save OLS covariance and Beta (solution vector, which is actually also saved in fit) @@ -101,15 +94,34 @@ } - // this section calculates the formal error on the weighted fit using the covariance values + // calculate the weight thresholds to mask the bad points: double Sum_Wx = 0.0; double Sum_Wy = 0.0; - if (1) { + for (i = 0; i < Npoints; i++) { + points[i].Wx = weight_cauchy(points[i].rx / points[i].dX); + points[i].Wy = weight_cauchy(points[i].ry / points[i].dY); + + Sum_Wx += points[i].Wx; + Sum_Wy += points[i].Wy; + } + double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints); + double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints); + + // set a mask (which can be used by the bootstrap resampling analysis) + for (i = 0; i < Npoints; i++) { + // keep if either is above threshold? + // drop if either is below threshold? + // points are marked as keep by default + if ((points[i].Wx < WxThreshold) || (points[i].Wy < WyThreshold)) { + points[i].mask = 1; // keep point if mask == 0 + } + } + + // 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 + if (data->getError) { double ax = 0.0, ay = 0.0; double bx = 0.0, by = 0.0; - double lambda = 0.0; + for (i = 0; i < Npoints; i++) { - points[i].Wx = weight_cauchy(points[i].rx / points[i].dX); - points[i].Wy = weight_cauchy(points[i].ry / points[i].dY); - ax += dpsi_cauchy(points[i].rx / points[i].dX); ay += dpsi_cauchy(points[i].ry / points[i].dY); @@ -117,19 +129,19 @@ bx += SQ(points[i].Wx); by += SQ(points[i].Wy); - - Sum_Wx += points[i].Wx; - Sum_Wy += points[i].Wy; } ax /= 1.0 * Npoints; // mean(psi_dot(r)) ay /= 1.0 * Npoints; - bx /= 1.0 * (Npoints - p); // mean(psi^2(r)) * (N / (N-p)) - by /= 1.0 * (Npoints - p); - - double sigma_robust_x = lambda * sqrt(bx) * sigma_hat * 2.385 / ax; - double sigma_robust_y = lambda * sqrt(by) * sigma_hat * 2.385 / ay; + 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), (n * SQ(sigma_robust_x) + SQ(p * sigma_ols)) / (n + SQ(p))); - double sigma_final_y = MAX(SQ(sigma_robust_y), (n * SQ(sigma_robust_y) + SQ(p * sigma_ols)) / (n + SQ(p))); + 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]); @@ -144,17 +156,4 @@ } - // set a mask (which can be used by the bootstrap resampling analysis) - double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints); - double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints); - - for (i = 0; i < Npoints; i++) { - // keep if either is above threshold? - // drop if either is below threshold? - // points are marked as keep by default - if ((points[i].Wx < WxThreshold) || (points[i].Wy < WyThreshold)) { - points[i].mask = 1; - } - } - // (optionally) add up the chi square for the fit, only counting the unmasked points double chisq = 0.0; @@ -162,16 +161,18 @@ for (i = 0; i < Npoints; i++) { if (points[i].mask) continue; + fit[0].Nfit ++; - double Xf = fit[0].Ro + fit[0].uR*points[i].T; - double Yf = fit[0].Do + fit[0].uD*points[i].T; - double wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dX); - double wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dY); - chisq += SQ(points[i].X - Xf) * wx; - chisq += SQ(points[i].Y - Yf) * wy; - fit[0].Nfit ++; - } - - // the reduced chisq is divided by (Ndof = 2*Nfit - 4) - fit[0].chisq = chisq / (2.0*fit[0].Nfit - 4.0); + if (data->getChisq) { + double Xf = fit[0].Ro + fit[0].uR*points[i].T; + double Yf = fit[0].Do + fit[0].uD*points[i].T; + double wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dX); + double wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dY); + chisq += SQ(points[i].X - Xf) * wx; + chisq += SQ(points[i].Y - Yf) * wy; + } + } + + // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms) + fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms); return (TRUE); Index: trunk/Ohana/src/relastro/src/FitPMandPar.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPMandPar.c (revision 39245) +++ trunk/Ohana/src/relastro/src/FitPMandPar.c (revision 39246) @@ -15,4 +15,7 @@ for (i = 0; i < Npoints; i++) { + if (points[i].mask) continue; // respect the mask if set + fit->Nfit ++; + /* handle case where dX or dY = 0.0 */ wx = 1.0 / SQ(points[i].dX); @@ -89,6 +92,4 @@ fit->dp = sqrt(data->A[4][4]); - fit->Nfit = Npoints; - return (TRUE); } Index: trunk/Ohana/src/relastro/src/FitPMandPar_IRLS.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPMandPar_IRLS.c (revision 39245) +++ trunk/Ohana/src/relastro/src/FitPMandPar_IRLS.c (revision 39246) @@ -11,15 +11,8 @@ int i,j; - int p = 5; - int n = 2 * Npoints; - int dof = n - p; - - // data->A,B,Cov,Beta,Beta_prev are allocated outside by FitAstromDataInit() - // points->Wx,Wy,rx,ry,u are elements of FitAstromPoint + int Ndof = 2 * Npoints - data->Nterms; // Convert the measurement errors into initial weights. for (i = 0; i < Npoints; i++) { - // points[i].Wx = 1 / points[i].dX; - // points[i].Wy = 1 / points[i].dY; points[i].Wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1 / SQ(points[i].dX); points[i].Wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1 / SQ(points[i].dY); @@ -39,6 +32,5 @@ sigma_ols += SQ(points[i].rx) + SQ(points[i].ry); } - sigma_ols = sqrt(sigma_ols / dof); - // sigma_ols = sqrt(sigma_ols / (Npts - 5.0)); + sigma_ols = sqrt(sigma_ols / (float)Ndof); // Save OLS covariance and Beta (solution vector, which is actually also saved in fit) @@ -99,18 +91,37 @@ } if (!converged) { - myAbort ("raise a warning on non-convergence"); - } - - // this section calculates the formal error on the weighted fit using the covariance values + fprintf (stderr, "raise a warning on non-convergence\n"); + } + + // calculate the weight thresholds to mask the bad points: double Sum_Wx = 0.0; double Sum_Wy = 0.0; - if (1) { + for (i = 0; i < Npoints; i++) { + points[i].Wx = weight_cauchy(points[i].rx / points[i].dX); + points[i].Wy = weight_cauchy(points[i].ry / points[i].dY); + + Sum_Wx += points[i].Wx; + Sum_Wy += points[i].Wy; + } + double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints); + double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints); + + // set a mask (which can be used by the bootstrap resampling analysis) + for (i = 0; i < Npoints; i++) { + // keep if either is above threshold? + // drop if either is below threshold? + // points are marked as keep by default + if ((points[i].Wx < WxThreshold) || (points[i].Wy < WyThreshold)) { + points[i].mask = 1; // keep point if mask == 0 + } + } + + // 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 + if (data->getError) { double ax = 0.0, ay = 0.0; double bx = 0.0, by = 0.0; - double lambda = 0.0; + for (i = 0; i < Npoints; i++) { - points[i].Wx = weight_cauchy(points[i].rx / points[i].dX); - points[i].Wy = weight_cauchy(points[i].ry / points[i].dY); - ax += dpsi_cauchy(points[i].rx / points[i].dX); ay += dpsi_cauchy(points[i].ry / points[i].dY); @@ -118,19 +129,19 @@ bx += SQ(points[i].Wx); by += SQ(points[i].Wy); - - Sum_Wx += points[i].Wx; - Sum_Wy += points[i].Wy; } ax /= 1.0 * Npoints; // mean(psi_dot(r)) ay /= 1.0 * Npoints; - bx /= 1.0 * (Npoints - p); // mean(psi^2(r)) * (N / (N-p)) - by /= 1.0 * (Npoints - p); - - double sigma_robust_x = lambda * sqrt(bx) * sigma_hat * 2.385 / ax; - double sigma_robust_y = lambda * sqrt(by) * sigma_hat * 2.385 / ay; + 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), (n * SQ(sigma_robust_x) + SQ(p * sigma_ols)) / (n + SQ(p))); - double sigma_final_y = MAX(SQ(sigma_robust_y), (n * SQ(sigma_robust_y) + SQ(p * sigma_ols)) / (n + SQ(p))); + 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]); @@ -147,17 +158,4 @@ } - // set a mask (which can be used by the bootstrap resampling analysis) - double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints); - double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints); - - for (i = 0; i < Npoints; i++) { - // keep if either is above threshold? - // drop if either is below threshold? - // points are marked as keep by default - if ((points[i].Wx < WxThreshold) || (points[i].Wy < WyThreshold)) { - points[i].mask = 1; - } - } - // (optionally) add up the chi square for the fit, only counting the unmasked points double chisq = 0.0; @@ -165,16 +163,18 @@ for (i = 0; i < Npoints; i++) { if (points[i].mask) continue; + fit[0].Nfit ++; - double Xf = fit[0].Ro + fit[0].uR*points[i].T + fit[0].p*points[i].pR; - double Yf = fit[0].Do + fit[0].uD*points[i].T + fit[0].p*points[i].pD; - double wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dX); - double wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dY); - chisq += SQ(points[i].X - Xf) * wx; - chisq += SQ(points[i].Y - Yf) * wy; - fit[0].Nfit ++; - } - - // the reduced chisq is divided by (Ndof = 2*Nfit - 4) - fit[0].chisq = chisq / (2.0*fit[0].Nfit - 5.0); + if (data->getChisq) { + double Xf = fit[0].Ro + fit[0].uR*points[i].T + fit[0].p*points[i].pR; + double Yf = fit[0].Do + fit[0].uD*points[i].T + fit[0].p*points[i].pD; + double wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dX); + double wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dY); + chisq += SQ(points[i].X - Xf) * wx; + chisq += SQ(points[i].Y - Yf) * wy; + } + } + + // the reduced chisq is divided by Ndof = (2*Nfit - Nterms) + fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms); return (TRUE); Index: trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 39245) +++ trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 39246) @@ -140,4 +140,63 @@ // use a fit with fewer parameters. + // fit the parallax + proper-motion model + // NOTE : we only fit PAR if we have already fitted for proper motion. if we do not fit PM or we fail + // to fit PM, we do not attempt PAR. thus failure to fit PAR falls back to PM-only + if (mode == FIT_PM_AND_PAR) { + if (Trange < PM_DT_MIN) { + mode = FIT_AVERAGE; + goto justPosition; + } + if (parRange < PAR_FACTOR_MIN) { + mode = FIT_PM_ONLY; + goto skipParallax; + } + if (fitStats->Npoints <= PAR_TOOFEW) { + mode = FIT_PM_ONLY; + goto skipParallax; + } + + // we are going to use the IRLS analysis to calculate the mean solution and the masking + // then run N_BOOTSTRAP_SAMPLES to measure the errors + fitStats->fitdataPar->getError = (N_BOOTSTRAP_SAMPLES < 3); + FitPMandPar_IRLS (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints, VERBOSE); + + if (N_BOOTSTRAP_SAMPLES >= 3) { + fitStats->Nfit = 0; + int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); + if (Nnomask < 5) { + myAbort ("handle this case, please"); + } + for (k = 0; k < fitStats->NfitAlloc; k++) { + BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask); + FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask); + fitStats->Nfit ++; + } + // these calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above) + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD); + BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX); + } + + // project Ro, Do back to RA,DEC + XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &fitStats->coords); + if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); + + average[0].flags |= ID_STAR_FIT_PAR; + fitStats->Npar ++; + + // XXX a hard-wired hack... + // unless there is a clear problems (below) with the parallax fit, we will use it + if ((fabs(fitPar.uR) > 4.0) || (fabs(fitPar.uD) > 4.0)) { + mode = FIT_PM_ONLY; + } else { + average[0].flags |= ID_STAR_USE_PAR; + } + } + +skipParallax: + // *** first fit for the proper motion (skip fit if Trange or Npts is too small) *** if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) { @@ -151,25 +210,29 @@ } - if (fitStats->NfitAlloc == 1) { - // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling: - if (1) { - FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints); - } else { - FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints, VERBOSE); + if (average[0].flags & ID_STAR_USE_PAR) { + if (!FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) { + myAbort ("oops"); } } else { - fitStats->Nfit = 0; - for (k = 0; k < fitStats->NfitAlloc; k++) { - BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints); - // if (!FitPM_IRLS (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints, VERBOSE)) continue; - if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue; - fitStats->Nfit ++; - } - BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); - BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); - BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR); - BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD); - } - FitAstromSetChisq (&fitPM, fitStats->points, fitStats->Npoints, FIT_PM_ONLY); + fitStats->fitdataPM->getError = (N_BOOTSTRAP_SAMPLES < 3); + FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints, VERBOSE); + + if (N_BOOTSTRAP_SAMPLES >= 3) { + fitStats->Nfit = 0; + int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); + if (Nnomask < 5) { + myAbort ("handle this case, please"); + } + for (k = 0; k < fitStats->NfitAlloc; k++) { + BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask); + if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue; + fitStats->Nfit ++; + } + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD); + } + } // project Ro, Do back to RA,DEC @@ -177,61 +240,19 @@ if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); - fitPM.p = fitPM.dp = 0.0; average[0].flags |= ID_STAR_FIT_PM; fitStats->Npm ++; // XXX a hard-wired hack... - if ((fabs(fitPM.uR) > 2.0) || (fabs(fitPM.uD) > 2.0)) { + // unless there is a clear problems (below) with the proper-motion fit or we have a parallax fit, we will use pm fit + if ((fabs(fitPM.uR) > 4.0) || (fabs(fitPM.uD) > 4.0)) { mode = FIT_AVERAGE; average[0].flags |= ID_STAR_BAD_PM; + } else { + if (!(average[0].flags & ID_STAR_USE_PAR)) { + average[0].flags |= ID_STAR_USE_PM; + } } } - // fit the parallax + proper-motion model - // NOTE : we only fit PAR if we have already fitted for proper motion. if we do not fit PM or we fail - // to fit PM, we do not attempt PAR. thus failure to fit PAR falls back to PM-only - if (mode == FIT_PM_AND_PAR) { - if (parRange < PAR_FACTOR_MIN) { - mode = FIT_PM_ONLY; - goto justPosition; - } - if (fitStats->Npoints <= PAR_TOOFEW) { - mode = FIT_PM_ONLY; - goto justPosition; - } - - if (fitStats->NfitAlloc == 1) { - // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling: - // FitPMandPar_IRLS (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints); - FitPMandPar (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints); - } else { - fitStats->Nfit = 0; - for (k = 0; k < fitStats->NfitAlloc; k++) { - BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints); - // FitPMandPar_IRLS (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, fitStats->Npoints); - FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, fitStats->Npoints); - fitStats->Nfit ++; - } - BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); - BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); - BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR); - BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD); - BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX); - } - FitAstromSetChisq (&fitPar, fitStats->points, fitStats->Npoints, FIT_PM_AND_PAR); - - // project Ro, Do back to RA,DEC - XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &fitStats->coords); - if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); - - average[0].flags |= ID_STAR_FIT_PAR; - fitStats->Npar ++; - - // XXX a hard-wired hack... - if ((fabs(fitPar.uR) > 2.0) || (fabs(fitPar.uD) > 2.0)) { - mode = FIT_PM_ONLY; - } - } - justPosition: { @@ -239,4 +260,6 @@ // if we only have one point, this is silly... + // XXX I need to rethink here : use a median for the position or weighted average? use IRLS anyway? + if (fitStats->NfitAlloc == 1) { FitAstromResultSetPM (&fitPos, 1, average); @@ -638,103 +661,4 @@ } -int BootstrapResample (FitAstromPoint *sample, FitAstromPoint *points, int Npoints) { - int i; - - // I need to draw Npoints random entries from 'points' with replacement: - for (i = 0; i < Npoints; i++) { - int N = Npoints * drand48(); - sample[i] = points[N]; - } - return TRUE; -} - -// calculate mean and sigma points for the 5 fit parameter -int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode) { - - // generate a histogram for the selected element - double *values = NULL; - ALLOCATE (values, double, Nfit); - - int i; - - for (i = 0; i < Nfit; i++) { - switch (mode) { - case FIT_RESULT_RA: - values[i] = fit[i].Ro; - break; - case FIT_RESULT_DEC: - values[i] = fit[i].Do; - break; - case FIT_RESULT_uR: - values[i] = fit[i].uR; - break; - case FIT_RESULT_uD: - values[i] = fit[i].uD; - break; - case FIT_RESULT_PLX: - values[i] = fit[i].p; - break; - default: - myAbort ("invalid option"); - } - } - - dsort (values, Nfit); - - double median; - if (Nfit % 2) { - int Ncenter = Nfit / 2; - median = values[Ncenter]; - } else { - int Ncenter = Nfit / 2 - 1; - median = 0.5*(values[Ncenter] + values[Ncenter + 1]); - } - - double Slo = VectorFractionInterpolate (values, 0.158655, Nfit); - double Shi = VectorFractionInterpolate (values, 0.841345, Nfit); - double sigma = (Shi - Slo) / 2.0; - - switch (mode) { - case FIT_RESULT_RA: - result->Ro = median; - result->dRo = sigma; - break; - case FIT_RESULT_DEC: - result->Do = median; - result->dDo = sigma; - break; - case FIT_RESULT_uR: - result->uR = median; - result->duR = sigma; - break; - case FIT_RESULT_uD: - result->uD = median; - result->duD = sigma; - break; - case FIT_RESULT_PLX: - result->p = median; - result->dp = sigma; - break; - default: - myAbort ("invalid option"); - } - - return TRUE; -} - -double VectorFractionInterpolate (double *values, float fraction, int Npts) { - - float F = fraction * Npts; - int N = fraction * Npts; - - if (N < 0 ) return NAN; - if (N >= Npts - 2) return NAN; - - // interpolate between N,N+1 - - double S = (F - N) * (values[N+1] - values[N]) + values[N]; - return S; -} - int FitAstromSetChisq (FitAstromResult *fit, FitAstromPoint *points, int Npoints, FitMode mode) { Index: trunk/Ohana/src/relastro/src/fitpm.c =================================================================== --- trunk/Ohana/src/relastro/src/fitpm.c (revision 39245) +++ trunk/Ohana/src/relastro/src/fitpm.c (revision 39246) @@ -12,7 +12,8 @@ # define OUT_ERROR 0.500 /* arcsec */ -# define N_STARS 20000 +# define N_STARS 3000 # define N_POINTS 100 -# define N_OUTLIERS 10 +# define N_OUTLIERS 100 +# define N_BOOTSTRAP 100 # define dcos(THETA) cos(RAD_DEG*THETA) @@ -24,4 +25,5 @@ FIT_PM_NOCLIP, FIT_PLX_IRLS, + FIT_PLX_BOOT, FIT_PLX_NOCLIP, }; @@ -46,4 +48,7 @@ if (!strcasecmp(argv[1], "plx-irls")) { mode = FIT_PLX_IRLS; + } + if (!strcasecmp(argv[1], "plx-boot")) { + mode = FIT_PLX_BOOT; } @@ -62,7 +67,12 @@ gaussdev_init (); - FitStats *fitStats = FitStatsInit (N_POINTS + N_OUTLIERS, 0); - - int i; + FitStats *fitStats = FitStatsInit (N_POINTS + N_OUTLIERS, N_BOOTSTRAP); + fitStats->fitdataPM->getChisq = TRUE; + fitStats->fitdataPM->getError = TRUE; + fitStats->fitdataPar->getChisq = TRUE; + fitStats->fitdataPar->getError = TRUE; + // NOTE: surprisingly, the chisq and error calculations only add ~2-3% to the exec time + + int i, k; int Nstars = N_STARS; @@ -83,4 +93,5 @@ break; case FIT_PLX_IRLS: + case FIT_PLX_BOOT: case FIT_PLX_NOCLIP: plx = 0.5*(drand48() + 0.0); @@ -112,4 +123,22 @@ FitPMandPar_IRLS (&fitPM, fitStats->fitdataPar, fitStats->points, fitStats->Npoints, 0); break; + + case FIT_PLX_BOOT: + FitPMandPar_IRLS (&fitPM, fitStats->fitdataPar, fitStats->points, fitStats->Npoints, 0); + fitStats->Nfit = 0; + int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); + for (k = 0; k < fitStats->NfitAlloc; k++) { + BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask); + FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask); + fitStats->Nfit ++; + } + // these calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above) + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD); + BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX); + break; + default: myAbort("oops"); @@ -118,7 +147,7 @@ XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &fitStats->coords); - fprintf (stdout, "%12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f\n", + fprintf (stdout, "%12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f | %8.6f %8.6f %8.6f %8.6f %8.6f %d | %f\n", Ro, Do, uR, uD, plx, - fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p, Tmean); + fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p, Tmean, fitPM.dRo, fitPM.dDo, fitPM.duR, fitPM.duD, fitPM.dp, fitPM.Nfit, fitPM.chisq); // ok (TRUE, "fitted star");