Index: trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- trunk/Ohana/src/relastro/include/relastro.h (revision 39693) +++ trunk/Ohana/src/relastro/include/relastro.h (revision 39926) @@ -33,5 +33,5 @@ typedef enum {OP_NONE, OP_IMAGES, OP_HIGH_SPEED, OP_MERGE_SOURCE, OP_UPDATE_OBJECTS, OP_UPDATE_OFFSETS, OP_LOAD_OBJECTS, OP_HPM, OP_PARALLEL_REGIONS, OP_PARALLEL_IMAGES, OP_REPAIR_STACKS, OP_REPAIR_WARPS, OP_REPAIR_OBJECT_ID} RelastroOp; -typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS} FitTarget; +typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, SET_CHIPS, TARGET_MOSAICS} FitTarget; typedef enum { @@ -144,4 +144,6 @@ int Nfit; int converged; + int useWeight; + } FitAstromResult; @@ -296,5 +298,5 @@ } CheckMeasureResult; -# define ID_MEAS_OBJECT_HAS_2MASS ID_MEAS_POOR_PHOTOM +// # define ID_MEAS_OBJECT_HAS_2MASS ID_MEAS_POOR_PHOTOM /* global variables set in parameter file */ @@ -315,4 +317,5 @@ int PARALLEL_REGIONS_MANUAL; char *MANUAL_UNIQUER; +int CATCH_UP; int HOST_ID; @@ -352,4 +355,7 @@ int VERBOSE2; +float TEST_SCALE; +char *GALAXY_MODEL; + int USE_FIXED_PIXCOORDS; int USE_GALAXY_MODEL; @@ -363,5 +369,18 @@ int USE_IMAGE_COORDS_FOR_REPAIR; int USE_ALL_IMAGES; +int KEEP_ALL_IMAGES_RA; int CHECK_MEASURE_TO_IMAGE; + +int SKIP_PS1_CHIP; +int SKIP_PS1_STACK; +int SKIP_HSC; +int SKIP_CFH; + +int UPDATE_PS1_STACK_MEASURE; +int UPDATE_PS1_CHIP_MEASURE; +int UPDATE_HSC_MEASURE; +int UPDATE_CFH_MEASURE; + +int APPLY_PROPER_MOTION; int RESET; @@ -383,4 +402,9 @@ int CHIPMAP; +int *ChipMapLoop; +int *ChipOrderLoop; +char *ChipMapLoopStr; +char *ChipOrderLoopStr; + int N_BOOTSTRAP_SAMPLES; @@ -407,7 +431,11 @@ float *LoopWeight2MASS; +char *LoopWeight2MASSstr; + float *LoopWeightTycho; -char *LoopWeight2MASSstr; char *LoopWeightTychostr; + +float *LoopWeightGAIA; +char *LoopWeightGAIAstr; int ImagSelect; @@ -433,4 +461,6 @@ FitMode FIT_MODE; +int USE_IRLS; +int ALLOW_IRLS; RelastroOp RELASTRO_OP; @@ -498,5 +528,5 @@ void initGridBins PROTO((Catalog *catalog, int Ncatalog)); void initImageBins PROTO((Catalog *catalog, int Ncatalog, int FULLINIT)); -void initImages PROTO((Image *input, off_t *line_number, off_t N)); +void initImages PROTO((Image *input, off_t *line_number, off_t N, int isSubset)); void freeImages PROTO((char *dbImagePtr)); void initMosaicBins PROTO((Catalog *catalog, int Ncatalog)); @@ -744,4 +774,6 @@ int isGPC1stack (int photcode); int isGPC1warp (int photcode); +int isHSCchip (int photcode); +int isCFHchip (int photcode); int save_astrom_table (); Index: trunk/Ohana/src/relastro/src/BootstrapOps.c =================================================================== --- trunk/Ohana/src/relastro/src/BootstrapOps.c (revision 39693) +++ trunk/Ohana/src/relastro/src/BootstrapOps.c (revision 39926) @@ -77,21 +77,21 @@ case FIT_RESULT_RA: // result->Ro = median; - result->dRo = sigma; + result->dRo = MAX(sigma, result->dRo); break; case FIT_RESULT_DEC: // result->Do = median; - result->dDo = sigma; + result->dDo = MAX(sigma, result->dDo); break; case FIT_RESULT_uR: // result->uR = median; - result->duR = sigma; + result->duR = MAX(sigma, result->duR); break; case FIT_RESULT_uD: // result->uD = median; - result->duD = sigma; + result->duD = MAX(sigma, result->duD); break; case FIT_RESULT_PLX: // result->p = median; - result->dp = sigma; + result->dp = MAX(sigma, result->dp); break; default: Index: trunk/Ohana/src/relastro/src/BrightCatalog.c =================================================================== --- trunk/Ohana/src/relastro/src/BrightCatalog.c (revision 39693) +++ trunk/Ohana/src/relastro/src/BrightCatalog.c (revision 39926) @@ -353,8 +353,8 @@ /*** MeasureTiny ***/ { - ohana_memcheck (1); + // ohana_memcheck (1); gfits_create_table_header (&theader, "BINTABLE", "MEASURE_TINY"); - ohana_memcheck (1); + // ohana_memcheck (1); gfits_define_bintable_column (&theader, "D", "RA", "ra", "degrees", 1.0, 0.0); @@ -429,6 +429,6 @@ gfits_set_bintable_column (&theader, &ftable, "RA", R, catalog->Nmeasure); - fprintf (stderr, "--------------- after set_bintable RA --------------"); - ohana_memdump_file (stderr, TRUE); + // fprintf (stderr, "--------------- after set_bintable RA --------------"); + // ohana_memdump_file (stderr, TRUE); gfits_set_bintable_column (&theader, &ftable, "DEC", D, catalog->Nmeasure); Index: trunk/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- trunk/Ohana/src/relastro/src/ConfigInit.c (revision 39693) +++ trunk/Ohana/src/relastro/src/ConfigInit.c (revision 39926) @@ -84,11 +84,4 @@ SetZeroPoint (25.0); - if (USE_GALAXY_MODEL) { - if (!InitGalaxyModel ("FEAST-HIPPARCOS")) { - fprintf (stderr, "failed to init galaxy model\n"); - exit (2); - } - } - FreeConfigFile(); free (config); Index: trunk/Ohana/src/relastro/src/FitAstromOps.c =================================================================== --- trunk/Ohana/src/relastro/src/FitAstromOps.c (revision 39693) +++ trunk/Ohana/src/relastro/src/FitAstromOps.c (revision 39926) @@ -257,4 +257,8 @@ fit->converged = FALSE; + // this is an input value + // if true, use the IRLS modified weight + fit->useWeight = FALSE; + return; } Index: trunk/Ohana/src/relastro/src/FitChip.c =================================================================== --- trunk/Ohana/src/relastro/src/FitChip.c (revision 39693) +++ trunk/Ohana/src/relastro/src/FitChip.c (revision 39926) @@ -78,5 +78,5 @@ order_use = MIN(MIN(order_use, CHIPMAP), 6); // can only go up to 6th order map (can be user limited) } else { - order_use = MIN(order_use, 3); // can only go up to 3rd order for polynomials + order_use = MIN(MIN(order_use, CHIPORDER), 3); // can only go up to 3rd order for polynomials } @@ -206,9 +206,26 @@ if (VERBOSE2) fprintf (stderr, "fit sigma: %f (%f, %f) : full: %f (%f, %f), scatter limit: %f (%d full, %d bright, %d fit, %d all) (%d %d %d %d %d)\n", dRsig, dLsig, dMsig, dRsigFull, dLsigFull, dMsigFull, dRmax, NstatFull, Nstat, fit[0].Npts, Nmatch, nMask1, nMask2, nMask3, nMask4, nMask5); - image[0].dXpixSys = dLsig; - image[0].dYpixSys = dMsig; + // need to convert dLsig, dMsig back to pixel scale (or up to arcsec) + + float plateScale; + if (image[0].coords.mosaic) { + // NOTE: for the full pixel to sky plate scale, use this: + // float plateScaleX = 3600.0*image[0].coords.mosaic->cdelt1*image[0].coords.cdelt1; + // float plateScaleY = 3600.0*image[0].coords.mosaic->cdelt2*image[0].coords.cdelt2; + + // since we are compare L,M values, just need to compensate for focal plate to sky: + float plateScaleX = 3600.0*fabs(image[0].coords.mosaic->cdelt1); + float plateScaleY = 3600.0*fabs(image[0].coords.mosaic->cdelt2); + plateScale = 0.5*(plateScaleX + plateScaleY); + } else { + // since we are compare L,M values, just need to compensate for arcsec vs degrees: + plateScale = 3600.0; + } + + image[0].dXpixSys = plateScale*dLsig; + image[0].dYpixSys = plateScale*dMsig; image[0].nFitAstrom = fit[0].Npts; - // fprintf (stderr, "%s %6.3f %4d %4d\n", image[0].name, image[0].refColor, Ncolor, image[0].nFitAstrom); + if (VERBOSE2) fprintf (stderr, "%s | %6.3f %6.3f | %4d %4d | %6.3f %6.3f\n", image[0].name, image[0].refColorRed, image[0].refColorBlue, Ncolor, image[0].nFitAstrom, image[0].dXpixSys, image[0].dYpixSys); if (fit) fit_free (fit); Index: trunk/Ohana/src/relastro/src/FitPM.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPM.c (revision 39693) +++ trunk/Ohana/src/relastro/src/FitPM.c (revision 39926) @@ -25,4 +25,5 @@ } + fit->useWeight = FALSE; // Ordinary Least Squares if (!FitPM_MinChisq (fit, data, points, Npoints)) return FALSE; if (!FitPM_SetChisq (fit, data, points, Npoints)) return FALSE; @@ -47,4 +48,5 @@ // Solve OLS equation + fit->useWeight = FALSE; // Ordinary Least Squares if (!FitPM_MinChisq(fit, data, points, Npoints)) { return(FALSE); @@ -69,5 +71,5 @@ // Iteratively reweight and solve - double sigma_hat = 0.0; // save for the error model + // double sigma_hat = 0.0; // save for the error model int converged = FALSE; int iterations = 0; @@ -90,4 +92,5 @@ // Solve with the new weights + fit->useWeight = TRUE; // Reweighted Least Squares if (!FitPM_MinChisq(fit, data, points, Npoints)) { @@ -104,5 +107,5 @@ points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY)); } - sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; + // sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; break; } @@ -119,5 +122,5 @@ points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY)); } - sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; + // sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; // Check convergence @@ -159,38 +162,15 @@ // 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; - - 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].dRo *= sigma_final_x; - fit[0].duR *= sigma_final_x; - fit[0].dDo *= sigma_final_y; - fit[0].duD *= sigma_final_y; + FitAstromResult fitErrors; + FitAstromResultInit (&fitErrors); + fitErrors.useWeight = FALSE; + + FitPM_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; } @@ -212,6 +192,11 @@ 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; Index: trunk/Ohana/src/relastro/src/FitPM_IRLS.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPM_IRLS.c (revision 39693) +++ (revision ) @@ -1,269 +1,0 @@ -# include "relastro.h" - -// These should probably be tunable: -# define MAX_ITERATIONS 10 -# define FIT_TOLERANCE 1e-4 -# define FLT_TOLERANCE 1e-6 -# define WEIGHT_THRESHOLD 0.3 - -int FitPM_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) { - - int i,j; - - int Ndof = 2 * Npoints - data->Nterms; - - // Convert the measurement errors into initial weights. - for (i = 0; i < Npoints; i++) { - 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); - } - - // Solve OLS equation - if (!FitPM_MinChisq(fit, data, points, Npoints)) { - return(FALSE); - } - - // Calculate r vector of residuals and least squares sigma - double sigma_ols = 0.0; - for (i = 0; i < Npoints; i++) { - points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro); - points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do); - sigma_ols += SQ(points[i].rx) + SQ(points[i].ry); - } - sigma_ols = sqrt(sigma_ols / (float)Ndof); - - // Save OLS covariance and Beta (solution vector, which is actually also saved in fit) - for (i = 0; i < data->Nterms; i++) { - for (j = 0; j < data->Nterms; j++) { - data->Cov[i][j] = data->A[i][j]; - } - data->Beta[i] = data->B[i][0]; - } - - // Iteratively reweight and solve - double sigma_hat = 0.0; // save for the error model - int converged = FALSE; - int iterations = 0; - - // modify the weight based on the distance from the previous fit. try up to MAX_ITERATIONS. - // at the end "fit", has the last fit parameters - for (iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations ++) { - // Save Beta. - for (i = 0; i < data->Nterms; i ++) { - data->Beta_prev[i] = data->Beta[i]; - } - - // Assign weights based on the deviation - 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); - } - - // Solve with the new weights - if (!FitPM_MinChisq(fit, data, points, Npoints)) { - - // restore the last solution and break - fit->Ro = data->Beta_prev[0]; - fit->uR = data->Beta_prev[1]; - fit->Do = data->Beta_prev[2]; - fit->uD = data->Beta_prev[3]; - - // calculate the residuals: - for (i = 0; i < Npoints; i++) { - points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro); - points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do); - points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY)); - } - sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; - break; - } - - // store the new Beta. - for (i = 0; i < data->Nterms; i++) { - data->Beta[i] = data->B[i][0]; - } - - // calculate the residuals: - for (i = 0; i < Npoints; i++) { - points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro); - points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do); - points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY)); - } - sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; - - // Check convergence - converged = TRUE; - for (i = 0; i < data->Nterms; i++) { - // if we are within FIT_TOLERANCE as a fractional error or FLT_TOLERANCE as an absolute error, we are good - if ((fabs(data->Beta[i] - data->Beta_prev[i]) > FIT_TOLERANCE * fabs(data->Beta[i])) && - (fabs(data->Beta[i] - data->Beta_prev[i]) > FLT_TOLERANCE)) { - converged = FALSE; - } - } - } - fit->converged = converged; - - // calculate the weight thresholds to mask the bad points: - double Sum_Wx = 0.0; - double Sum_Wy = 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); - - 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; - - 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].dRo *= sigma_final_x; - fit[0].duR *= sigma_final_x; - fit[0].dDo *= sigma_final_y; - fit[0].duD *= sigma_final_y; - } - - // count unmasked points and (optionally) add up the chi square for the fit - double chisq = 0.0; - fit[0].Nfit = 0; - for (i = 0; i < Npoints; i++) { - if (points[i].mask) continue; - fit[0].Nfit ++; - - 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); -} - -int FitPM_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) { - - myAssert (data->Nterms == 4, "invalid fit arrays"); - - int i; - double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT; - Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0; - - for (i = 0; i < Npoints; i++) { - if (points[i].mask) continue; // respect the mask if set - - wx = points[i].Wx; - wy = points[i].Wy; - - Wx += wx; - Wy += wy; - - double TWx = points[i].T*wx; - double TWy = points[i].T*wy; - - Tx += TWx; - Ty += TWy; - - Tx2 += points[i].T*TWx; - Ty2 += points[i].T*TWy; - - Xs += points[i].X*wx; - Ys += points[i].Y*wy; - - XT += points[i].X*TWx; - YT += points[i].Y*TWy; - } - - // X^T W X - data->A[0][0] = Wx; - data->A[0][1] = Tx; - data->A[0][2] = 0.0; - data->A[0][3] = 0.0; - - data->A[1][0] = Tx; - data->A[1][1] = Tx2; - data->A[1][2] = 0.0; - data->A[1][3] = 0.0; - - data->A[2][0] = 0.0; - data->A[2][1] = 0.0; - data->A[2][2] = Wy; - data->A[2][3] = Ty; - - data->A[3][0] = 0.0; - data->A[3][1] = 0.0; - data->A[3][2] = Ty; - data->A[3][3] = Ty2; - - // X^T W Y - data->B[0][0] = Xs; - data->B[1][0] = XT; - data->B[2][0] = Ys; - data->B[3][0] = YT; - - if (!dgaussjordan (data->A, data->B, 4, 1)) { - return FALSE; - } - - fit->Ro = data->B[0][0]; - fit->uR = data->B[1][0]; - fit->Do = data->B[2][0]; - fit->uD = data->B[3][0]; - fit->p = 0.0; - - fit->dRo = sqrt(data->A[0][0]); - fit->duR = sqrt(data->A[1][1]); - fit->dDo = sqrt(data->A[2][2]); - fit->duD = sqrt(data->A[3][3]); - fit->dp = 0.0; - - return TRUE; -} - Index: trunk/Ohana/src/relastro/src/FitPMandPar.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPMandPar.c (revision 39693) +++ trunk/Ohana/src/relastro/src/FitPMandPar.c (revision 39926) @@ -25,4 +25,5 @@ } + 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,4 +48,5 @@ // 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); @@ -69,5 +71,5 @@ // Iteratively reweight and solve - double sigma_hat = 0.0; // save for the error model + // double sigma_hat = 0.0; // save for the error model int converged = FALSE; int iterations = 0; @@ -90,4 +92,5 @@ // Solve with the new weights + fit->useWeight = TRUE; // Reweighted Least Squares if (!FitPMandPar_MinChisq(fit, data, points, Npoints)) { @@ -105,5 +108,5 @@ points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY)); } - sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; + // sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; break; } @@ -120,5 +123,5 @@ points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY)); } - sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; + // sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; // Check convergence @@ -157,43 +160,19 @@ } - // 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,6 +198,11 @@ 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; Index: trunk/Ohana/src/relastro/src/FitPMandPar_IRLS.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPMandPar_IRLS.c (revision 39693) +++ (revision ) @@ -1,302 +1,0 @@ -# include "relastro.h" -# define OLD_METHOD 1 - -// These should probably be tunable: -# define MAX_ITERATIONS 10 -# define FIT_TOLERANCE 1e-4 -# define FLT_TOLERANCE 1e-6 -# define WEIGHT_THRESHOLD 0.3 - -int FitPMandPar_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) { - - int i,j; - - int Ndof = 2 * Npoints - data->Nterms; - - // Convert the measurement errors into initial weights. - for (i = 0; i < Npoints; i++) { -# if (OLD_METHOD) - 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); -# else - points[i].Qx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1 / SQ(points[i].dX); - points[i].Qy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1 / SQ(points[i].dY); -# endif - points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0 - points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0 - } - - // Solve OLS equation: failure here means the chisq matrix is degenerate, give up entirely - if (!weighted_LS_PLX(fit, data, points, Npoints, VERBOSE)) { - return(FALSE); - } - - // Calculate r vector of residuals and least squares sigma - double sigma_ols = 0.0; - for (i = 0; i < Npoints; i++) { - points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro + points[i].pR * fit->p); - points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do + points[i].pD * fit->p); - sigma_ols += SQ(points[i].rx) + SQ(points[i].ry); - } - sigma_ols = sqrt(sigma_ols / (float)Ndof); - - // Save OLS covariance and Beta (solution vector, which is actually also saved in fit) - for (i = 0; i < data->Nterms; i++) { - for (j = 0; j < data->Nterms; j++) { - data->Cov[i][j] = data->A[i][j]; - } - data->Beta[i] = data->B[i][0]; - } - - // Iteratively reweight and solve - double sigma_hat = 0.0; // save for the error model - int converged = FALSE; - int iterations = 0; - - // modify the weight based on the distance from the previous fit. try up to MAX_ITERATIONS. - // at the end "fit", has the last fit parameters - for (iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations ++) { - // Save Beta. - for (i = 0; i < data->Nterms; i ++) { - data->Beta_prev[i] = data->Beta[i]; - } - - // Assign weights based on the deviation - 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); - points[i].qx = points[i].Wx * points[i].Qx; - points[i].qy = points[i].Wy * points[i].Qy; - } - - // Solve with the new weights - if (!weighted_LS_PLX(fit, data, points, Npoints, VERBOSE)) { - - // restore the last solution and break - fit->Ro = data->Beta_prev[0]; - fit->uR = data->Beta_prev[1]; - fit->Do = data->Beta_prev[2]; - fit->uD = data->Beta_prev[3]; - fit->p = data->Beta_prev[4]; - - // calculate the residuals: - for (i = 0; i < Npoints; i++) { - points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro + points[i].pR * fit->p); - points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do + points[i].pD * fit->p); - points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY)); - } - sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; - break; - } - - // store the new Beta. - for (i = 0; i < data->Nterms; i++) { - data->Beta[i] = data->B[i][0]; - } - - // calculate the residuals: - for (i = 0; i < Npoints; i++) { - points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro + points[i].pR * fit->p); - points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do + points[i].pD * fit->p); - points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY)); - } - sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; - - // Check convergence - converged = TRUE; - for (i = 0; i < data->Nterms; i++) { - // if we are within FIT_TOLERANCE as a fractional error or FLT_TOLERANCE as an absolute error, we are good - if ((fabs(data->Beta[i] - data->Beta_prev[i]) > FIT_TOLERANCE * fabs(data->Beta[i])) && - (fabs(data->Beta[i] - data->Beta_prev[i]) > FLT_TOLERANCE)) { - converged = FALSE; - } - } - } - fit->converged = converged; - - // calculate the weight thresholds to mask the bad points: - double Sum_Wx = 0.0; - double Sum_Wy = 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); - - 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; - - 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].qx); - by += SQ(points[i].qy); - } - 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); - } - - // count unmasked points and (optionally) add up the chi square for the fit - double chisq = 0.0; - fit[0].Nfit = 0; - for (i = 0; i < Npoints; i++) { - if (points[i].mask) continue; - fit[0].Nfit ++; - - 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); -} - -int weighted_LS_PLX (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) { - - myAssert (data->Nterms == 5, "invalid fit arrays"); - - int i; - double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT; - double PR, PD, PRT, PDT, PRX, PDY, PR2, PD2; - - PR = PD = PRT = PDT = PRX = PDY = PR2 = PD2 = 0.0; - Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0; - - for (i = 0; i < Npoints; i++) { - - // if (VERBOSE == 2) fprintf (stderr, "%f %f : %f %f : %f : %f %f\n", X[i], WX[i], Y[i], WY[i], T[i], pR[i], pD[i]); - -# if (OLD_METHOD) - wx = points[i].Wx; - wy = points[i].Wy; -# else - wx = points[i].Wx; - wy = points[i].Wy; -# endif - - Wx += wx; - Wy += wy; - - Tx += points[i].T*wx; - Ty += points[i].T*wy; - - Tx2 += SQ(points[i].T)*wx; - Ty2 += SQ(points[i].T)*wy; - - PR += points[i].pR*wx; - PD += points[i].pD*wy; - - PRT += points[i].pR*points[i].T*wx; - PDT += points[i].pD*points[i].T*wy; - - PRX += points[i].pR*points[i].X*wx; - PDY += points[i].pD*points[i].Y*wy; - - PR2 += SQ(points[i].pR)*wx; - PD2 += SQ(points[i].pD)*wy; - - Xs += points[i].X*wx; - Ys += points[i].Y*wy; - - XT += points[i].X*points[i].T*wx; - YT += points[i].Y*points[i].T*wy; - } - - data->A[0][0] = Wx; - data->A[0][1] = Tx; - data->A[0][2] = 0.0; - data->A[0][3] = 0.0; - data->A[0][4] = PR; - - data->A[1][0] = Tx; - data->A[1][1] = Tx2; - data->A[1][2] = 0.0; - data->A[1][3] = 0.0; - data->A[1][4] = PRT; - - data->A[2][0] = 0.0; - data->A[2][1] = 0.0; - data->A[2][2] = Wy; - data->A[2][3] = Ty; - data->A[2][4] = PD; - - data->A[3][0] = 0.0; - data->A[3][1] = 0.0; - data->A[3][2] = Ty; - data->A[3][3] = Ty2; - data->A[3][4] = PDT; - - data->A[4][0] = PR; - data->A[4][1] = PRT; - data->A[4][2] = PD; - data->A[4][3] = PDT; - data->A[4][4] = PR2 + PD2; - - data->B[0][0] = Xs; - data->B[1][0] = XT; - data->B[2][0] = Ys; - data->B[3][0] = YT; - data->B[4][0] = PRX + PDY; - - if (!dgaussjordan (data->A, data->B, 5, 1)) { - return FALSE; - } - - fit->Ro = data->B[0][0]; - fit->uR = data->B[1][0]; - fit->Do = data->B[2][0]; - fit->uD = data->B[3][0]; - fit->p = data->B[4][0]; - - return TRUE; -} Index: trunk/Ohana/src/relastro/src/FitPosPMfixed.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPosPMfixed.c (revision 39693) +++ trunk/Ohana/src/relastro/src/FitPosPMfixed.c (revision 39926) @@ -54,4 +54,5 @@ } + fit->useWeight = FALSE; // Ordinary Least Squares if (!FitPosPMfixed_MinChisq (fit, data, points, Npoints)) { if (!FitPosPMfixed_Single (fit, points, Npoints)) return FALSE; @@ -83,4 +84,5 @@ // Solve OLS equation + fit->useWeight = FALSE; // Ordinary Least Squares if (!FitPosPMfixed_MinChisq(fit, data, points, Npoints)) { return(FALSE); @@ -105,5 +107,5 @@ // Iteratively reweight and solve - double sigma_hat = 0.0; // save for the error model + // double sigma_hat = 0.0; // save for the error model int converged = FALSE; int iterations = 0; @@ -126,4 +128,5 @@ // Solve with the new weights + fit->useWeight = TRUE; // Reweighted Least Squares if (!FitPosPMfixed_MinChisq(fit, data, points, Npoints)) { @@ -138,5 +141,5 @@ points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY)); } - sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; + // sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; break; } @@ -153,5 +156,5 @@ points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY)); } - sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; + // sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; // Check convergence @@ -193,34 +196,16 @@ // 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; - - 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].dDo = sqrt(data->Cov[1][1]); - - fit[0].dRo *= sigma_final_x; - fit[0].dDo *= sigma_final_y; + FitAstromResult fitErrors; + FitAstromResultInit (&fitErrors); + fitErrors.useWeight = FALSE; + + FitPosPMfixed_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; } @@ -242,6 +227,11 @@ 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; Index: trunk/Ohana/src/relastro/src/FitPosPMfixed_IRLS.c =================================================================== --- trunk/Ohana/src/relastro/src/FitPosPMfixed_IRLS.c (revision 39693) +++ (revision ) @@ -1,228 +1,0 @@ -# include "relastro.h" - -// These should probably be tunable: -# define MAX_ITERATIONS 10 -# define FIT_TOLERANCE 1e-4 -# define FLT_TOLERANCE 1e-6 -# define WEIGHT_THRESHOLD 0.3 - -int FitPosPMfixed_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) { - - int i,j; - - int Ndof = 2 * Npoints - data->Nterms; - - // Convert the measurement errors into initial weights. - for (i = 0; i < Npoints; i++) { - 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); - } - - // Solve OLS equation - if (!weighted_LS_POS(fit, data, points, Npoints, VERBOSE)) { - return(FALSE); - } - - // Calculate r vector of residuals and least squares sigma - double sigma_ols = 0.0; - for (i = 0; i < Npoints; i++) { - points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro); - points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do); - sigma_ols += SQ(points[i].rx) + SQ(points[i].ry); - } - sigma_ols = sqrt(sigma_ols / (float)Ndof); - - // Save OLS covariance and Beta (solution vector, which is actually also saved in fit) - for (i = 0; i < data->Nterms; i++) { - for (j = 0; j < data->Nterms; j++) { - data->Cov[i][j] = data->A[i][j]; - } - data->Beta[i] = data->B[i][0]; - } - - // Iteratively reweight and solve - double sigma_hat = 0.0; // save for the error model - int converged = FALSE; - int iterations = 0; - - // modify the weight based on the distance from the previous fit. try up to MAX_ITERATIONS. - // at the end "fit", has the last fit parameters - for (iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations ++) { - // Save Beta. - for (i = 0; i < data->Nterms; i ++) { - data->Beta_prev[i] = data->Beta[i]; - } - - // Assign weights based on the deviation - 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); - } - - // Solve with the new weights - if (!weighted_LS_POS(fit, data, points, Npoints, VERBOSE)) { - - // restore the last solution and break - fit->Ro = data->Beta_prev[0]; - fit->Do = data->Beta_prev[1]; - - // calculate the residuals: - for (i = 0; i < Npoints; i++) { - points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro); - points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do); - points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY)); - } - sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; - break; - } - - // store the new Beta. - for (i = 0; i < data->Nterms; i++) { - data->Beta[i] = data->B[i][0]; - } - - // calculate the residuals: - for (i = 0; i < Npoints; i++) { - points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro); - points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do); - points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY)); - } - sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745; - - // Check convergence - converged = TRUE; - for (i = 0; i < data->Nterms; i++) { - // if we are within FIT_TOLERANCE as a fractional error or FLT_TOLERANCE as an absolute error, we are good - if ((fabs(data->Beta[i] - data->Beta_prev[i]) > FIT_TOLERANCE * fabs(data->Beta[i])) && - (fabs(data->Beta[i] - data->Beta_prev[i]) > FLT_TOLERANCE)) { - converged = FALSE; - } - } - } - fit->converged = converged; - - // calculate the weight thresholds to mask the bad points: - double Sum_Wx = 0.0; - double Sum_Wy = 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); - - 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; - - 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].dDo = sqrt(data->Cov[1][1]); - - fit[0].dRo *= sigma_final_x; - fit[0].dDo *= sigma_final_y; - } - - // count unmasked points and (optionally) add up the chi square for the fit - double chisq = 0.0; - fit[0].Nfit = 0; - for (i = 0; i < Npoints; i++) { - if (points[i].mask) continue; - fit[0].Nfit ++; - - 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); -} - -int weighted_LS_POS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) { - - myAssert (data->Nterms == 2, "invalid fit arrays"); - - int i; - double wx, wy, Wx, Wy, Tx, Ty, Xs, Ys; - - Wx = Wy = Tx = Ty = Xs = Ys = 0.0; - - for (i = 0; i < Npoints; i++) { - wx = points[i].Wx; - wy = points[i].Wy; - - Wx += wx; - Wy += wy; - - Tx += points[i].T*wx; - Ty += points[i].T*wy; - - Xs += points[i].X*wx; - Ys += points[i].Y*wy; - } - - // X^T W X - data->A[0][0] = Wx; - data->A[0][1] = 0.0; - - data->A[1][0] = 0.0; - data->A[1][1] = Wy; - - // X^T W Y - data->B[0][0] = Xs - fit->uR*Tx; - data->B[1][0] = Ys - fit->uD*Ty; - - if (!dgaussjordan (data->A, data->B, 2, 1)) { - return FALSE; - } - - fit->Ro = data->B[0][0]; - fit->Do = data->B[1][0]; - fit->p = 0.0; - - return TRUE; -} Index: trunk/Ohana/src/relastro/src/GetAstromError.c =================================================================== --- trunk/Ohana/src/relastro/src/GetAstromError.c (revision 39693) +++ trunk/Ohana/src/relastro/src/GetAstromError.c (revision 39926) @@ -48,4 +48,13 @@ if (isnan(code[0].astromErrSys)) return NAN; + if (measure[0].photcode == 1030) { + if (mode == ERROR_MODE_RA) { + dPobs = pow(10.0, (6.0 * measure[0].dXccd - 3.0)); // dXccd is a value in pixels + } + if (mode == ERROR_MODE_DEC) { + dPobs = pow(10.0, (6.0 * measure[0].dYccd - 3.0)); // dXccd is a value in pixels + } + } + AS = code[0].astromErrScale; MS = code[0].astromErrMagScale; @@ -69,7 +78,11 @@ // to match the 2MASS / Tycho / ICRS reference frame. As Nloop gets higher, the weight // needs to drop to allow the ps1 measurements to drive the solution + int isGAIA = USE_GALAXY_MODEL && !isImage && (measure[0].photcode == 1030); int is2MASS = USE_GALAXY_MODEL && !isImage && (measure[0].photcode >= 2011) && (measure[0].photcode <= 2013); int isTycho = USE_GALAXY_MODEL && !isImage && (measure[0].photcode >= 2020) && (measure[0].photcode <= 2021); + + int hasGAIA = USE_GALAXY_MODEL && isImage && (measure[0].dbFlags & ID_MEAS_OBJECT_HAS_GAIA); int has2MASS = USE_GALAXY_MODEL && isImage && (measure[0].dbFlags & ID_MEAS_OBJECT_HAS_2MASS); + int hasTycho = USE_GALAXY_MODEL && isImage && (measure[0].dbFlags & ID_MEAS_OBJECT_HAS_TYCHO); // modest hack: if the object has 2MASS or Tycho, we set this internal bit and adjust the @@ -79,6 +92,16 @@ dPtotal = dPtotal / LoopWeight2MASS[Nloop]; } + if (hasGAIA && LoopWeightGAIA && (Nloop >= 0)) { + dPtotal = dPtotal / LoopWeightGAIA[Nloop]; + } + if (hasTycho && LoopWeightTycho && (Nloop >= 0)) { + dPtotal = dPtotal / LoopWeightTycho[Nloop]; + } + if (is2MASS && LoopWeight2MASS && (Nloop >= 0)) { dPtotal = dPtotal / LoopWeight2MASS[Nloop]; + } + if (isGAIA && LoopWeightGAIA && (Nloop >= 0)) { + dPtotal = dPtotal / LoopWeightGAIA[Nloop]; } if (isTycho && LoopWeightTycho && (Nloop >= 0)) { @@ -120,4 +143,13 @@ // do not raise an exception, just send back the result if (isnan(code[0].astromErrSys)) return NAN; + + if (measure[0].photcode == 1030) { + if (mode == ERROR_MODE_RA) { + dPobs = pow(10.0, (6.0 * measure[0].dXccd - 3.0)); // dXccd is a value in pixels + } + if (mode == ERROR_MODE_DEC) { + dPobs = pow(10.0, (6.0 * measure[0].dYccd - 3.0)); // dXccd is a value in pixels + } + } AS = code[0].astromErrScale; @@ -143,7 +175,11 @@ // to match the 2MASS / Tycho / ICRS reference frame. As Nloop gets higher, the weight // needs to drop to allow the ps1 measurements to drive the solution + int isGAIA = USE_GALAXY_MODEL && !isImage && (measure[0].photcode == 1030); int is2MASS = USE_GALAXY_MODEL && !isImage && (measure[0].photcode >= 2011) && (measure[0].photcode <= 2013); int isTycho = USE_GALAXY_MODEL && !isImage && (measure[0].photcode >= 2020) && (measure[0].photcode <= 2021); + + int hasGAIA = USE_GALAXY_MODEL && isImage && (measure[0].dbFlags & ID_MEAS_OBJECT_HAS_GAIA); int has2MASS = USE_GALAXY_MODEL && isImage && (measure[0].dbFlags & ID_MEAS_OBJECT_HAS_2MASS); + int hasTycho = USE_GALAXY_MODEL && isImage && (measure[0].dbFlags & ID_MEAS_OBJECT_HAS_TYCHO); // modest hack: if the object has 2MASS or Tycho, we set this internal bit and adjust the @@ -153,6 +189,16 @@ dPtotal = dPtotal / LoopWeight2MASS[Nloop]; } + if (hasGAIA && LoopWeightGAIA && (Nloop >= 0)) { + dPtotal = dPtotal / LoopWeightGAIA[Nloop]; + } + if (hasTycho && LoopWeightTycho && (Nloop >= 0)) { + dPtotal = dPtotal / LoopWeightTycho[Nloop]; + } + if (is2MASS && LoopWeight2MASS) { dPtotal = dPtotal / LoopWeight2MASS[Nloop]; + } + if (isGAIA && LoopWeightGAIA && (Nloop >= 0)) { + dPtotal = dPtotal / LoopWeightGAIA[Nloop]; } if (isTycho && LoopWeightTycho) { Index: trunk/Ohana/src/relastro/src/ImageOps.c =================================================================== --- trunk/Ohana/src/relastro/src/ImageOps.c (revision 39693) +++ trunk/Ohana/src/relastro/src/ImageOps.c (revision 39926) @@ -6,4 +6,5 @@ static Image *image; // list of available images static off_t Nimage; // number of available images +static int isImageSubset; // if we read only a subset of the rows from the Image FITS, LineNumber tells us to which row @@ -78,8 +79,9 @@ } -void initImages (Image *input, off_t *line_number, off_t N) { +void initImages (Image *input, off_t *line_number, off_t N, int isSubset) { off_t i; + isImageSubset = isSubset; image = input; LineNumber = line_number; @@ -136,5 +138,5 @@ // we call gfits_db_free as well as this function. sometimes those point at the same // memory location, in which case we should only do the free once. - if ((void *) dbImagePtr != (void *) image) free (image); + if (((void *) dbImagePtr != (void *) image) && isImageSubset) free (image); free_astrom_table(); } @@ -793,6 +795,13 @@ ref[i].D = catalog[c].average[n].D; + int XVERB = FALSE; + XVERB |= (catalog[c].average[n].objID == OBJ_ID_SRC) && (catalog[c].average[n].catID == CAT_ID_SRC); + XVERB |= (catalog[c].average[n].objID == OBJ_ID_DST) && (catalog[c].average[n].catID == CAT_ID_DST); + if (XVERB) { + fprintf (stderr, "found test object\n"); + } + // if we are applying the galaxy model, move the reference position... - if (USE_GALAXY_MODEL) { + if (APPLY_PROPER_MOTION) { // apply proper-motion from average position to measure epoch: float dTime = (measure[0].t - catalog[c].average[n].Tmean) / (86400*365.25) ; // time relative to Tmean in years @@ -934,371 +943,2 @@ } -# if (0) -/** lifted from relphot/StarOps.clean_measures */ -void FlagOutliers2D(Catalog *catalog); - -// operates on Full values (not tiny) -void FlagOutliers (Catalog *catalog) { - - // XXX FlagOutliers is now just using FlagOutliers2D - FlagOutliers2D(catalog); - return; - - int Ndel, Nave; - off_t i, j, k, m, N, Nmax, TOOFEW, Nsecfilt; - double Ns, theta, x, y; - double *R, *D, *dR, *dD; - StatType statsR, statsD; - - Nsecfilt = GetPhotcodeNsecfilt(); - assert(catalog[0].Nsecfilt == Nsecfilt); - - if (VERBOSE2) fprintf (stderr, "marking poor measures\n"); - Nmax = 0; - for (i = 0; i < catalog[0].Naverage; i++) { - Nmax = MAX (Nmax, catalog[0].average[i].Nmeasure); - } - - ALLOCATE (R, double, Nmax); - ALLOCATE (D, double, Nmax); - ALLOCATE (dR, double, Nmax); - ALLOCATE (dD, double, Nmax); - - /* it makes no sense to mark 3-sigma outliers with <5 measurements */ - TOOFEW = MAX (5, SRC_MEAS_TOOFEW); - - Ns = CLIP_THRESH; - Ndel = Nave = 0; - - /* loop over each object in the catalog */ - for (j = 0; j < catalog[0].Naverage; j++) { - - // pointer to this set of measurements - m = catalog[0].average[j].measureOffset; - Measure *measure = &catalog[0].measure[m]; - - /* accumulate list of valid measurements */ - N = 0; - for (k = 0; k < catalog[0].average[j].Nmeasure; k++) { - // skip measurements based on user selected criteria - if (!MeasFilterTest(&measure[k], FALSE)) continue; - R[N] = measure[k].R; - D[N] = measure[k].D; - dR[N] = GetAstromError (&measure[k], ERROR_MODE_RA); - dD[N] = GetAstromError (&measure[k], ERROR_MODE_DEC); - if (isnan(R[N]) || isnan(D[N])) continue; - N++; - } - if (N <= TOOFEW) continue; - - /* 3-sigma clip based on stats of inner 50% */ - initstats ("MEAN"); - liststats (R, dR, N, &statsR); - liststats (D, dD, N, &statsD); - - statsR.sigma = MAX (MIN_ERROR, statsR.sigma); - statsD.sigma = MAX (MIN_ERROR, statsD.sigma); - - /* compare per-object distance to this standard deviation, and flag outliers*/ - N = 0; - for (k = 0; k < catalog[0].average[j].Nmeasure; k++) { - // reset flag on each invocation - measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; - - // skip measurements based on user selected criteria - if (!MeasFilterTest(&measure[k], FALSE)) continue; - - x = measure[k].R - statsR.median; - y = measure[k].D - statsD.median; - theta = atan2(y,x); - if ((x*x + y*y) > (SQR(statsR.sigma * Ns * cos(theta)) + - SQR(statsD.sigma * Ns * sin(theta)))) { - measure[k].dbFlags |= ID_MEAS_POOR_ASTROM; - Ndel++; - } - N++; - Nave ++; - } - - // examine results - // relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, catalog[0].average[j].Nmeasure, statsR, statsD, Ns); - } - - if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave); - free (R); - free (dR); - free (D); - free (dD); -} - - -/** an alternative outlier rejection scheme */ -void FlagOutliers2D (Catalog *catalog) { - - int Ndel, Nave; - off_t i, j, k, m, N, Nmax, TOOFEW, Nsecfilt; - double *index; - double Ns, theta, x, y; - double *R, *D, *dR, *dD, *d2; - StatType statsR, statsD; - - // XXX we are not going to use this for now - return; - - Nsecfilt = GetPhotcodeNsecfilt(); - assert(catalog[0].Nsecfilt == Nsecfilt); - - if (VERBOSE2) fprintf (stderr, "marking poor measures\n"); - Nmax = 0; - for (i = 0; i < catalog[0].Naverage; i++) { - Nmax = MAX (Nmax, catalog[0].average[i].Nmeasure); - } - - ALLOCATE (R, double, Nmax); - ALLOCATE (D, double, Nmax); - ALLOCATE (dR, double, Nmax); - ALLOCATE (dD, double, Nmax); - ALLOCATE (d2, double, Nmax); - ALLOCATE (index, double, Nmax); - - /* it makes no sense to mark 3-sigma outliers with <5 measurements */ - TOOFEW = MAX (5, SRC_MEAS_TOOFEW); - - Ns = CLIP_THRESH; - Ndel = Nave = 0; - - /* loop over each object in the catalog */ - for (j = 0; j < catalog[0].Naverage; j++) { - - // pointer to this set of measurements - m = catalog[0].average[j].measureOffset; - Measure *measure = &catalog[0].measure[m]; - - /* accumulate list of valid measurements */ - N = 0; - for (k = 0; k < catalog[0].average[j].Nmeasure; k++) { - - // reset flag on each invocation - measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; - - // skip measurements based on user selected criteria - if (!MeasFilterTest(&measure[k], FALSE)) continue; - R[N] = measure[k].R; - D[N] = measure[k].D; - dR[N] = GetAstromError(&measure[k], ERROR_MODE_RA); - dD[N] = GetAstromError(&measure[k], ERROR_MODE_DEC); - if (isnan(R[N]) || isnan(D[N])) continue; - N++; - } - if (N <= TOOFEW) continue; - - /* calculate mean of all points*/ - initstats ("MEAN"); - liststats (R, dR, N, &statsR); - liststats (D, dD, N, &statsD); - statsR.sigma = MAX (MIN_ERROR, statsR.sigma); - statsD.sigma = MAX (MIN_ERROR, statsD.sigma); - - /* calculate deviations of all points*/ - N = 0; - for (k = 0; k < catalog[0].average[j].Nmeasure; k++) { - // skip bad measurements - if (!MeasFilterTest(&measure[k], FALSE)) continue; - x = measure[k].R - statsR.median; - y = measure[k].D - statsD.median; - theta = atan2(y,x); - d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) + - SQR(statsD.sigma * Ns * sin(theta))); - index[N] = k; - N++; - } - - // sort d2 - dsortpair(d2, index, N); - N = (N/2 > (N-1)) ? N/2 : N-1; - - // recalculate image center, sigma based on closest 50% of points - for (k = 0; k < N; k++) { - off_t ind = (off_t) index[k]; - R[k] = measure[ind].R; - D[k] = measure[ind].D; - dR[k] = GetAstromError(&measure[ind], ERROR_MODE_RA); - dD[k] = GetAstromError(&measure[ind], ERROR_MODE_DEC); - } - liststats (R, dR, N, &statsR); - liststats (D, dD, N, &statsD); - statsR.sigma = MAX (MIN_ERROR, statsR.sigma); - statsD.sigma = MAX (MIN_ERROR, statsD.sigma); - - // use these new statistics to flag outliers - N = 0; - for (k = 0; k < catalog[0].average[j].Nmeasure; k++) { - //skip bad measurements - if (!MeasFilterTest(&measure[k], FALSE)) continue; - x = measure[k].R - statsR.median; - y = measure[k].D - statsD.median; - theta = atan2(y,x); - d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) + - SQR(statsD.sigma * Ns * sin(theta))); - if ((d2[N]) > 1) { - measure[k].dbFlags |= ID_MEAS_POOR_ASTROM; - Ndel ++; - } - N++; - Nave++; - } // done rejecting outliers - - // examine results - // relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, catalog[0].average[j].Nmeasure, statsR, statsD, Ns); - - } // done looping over objects - - if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave); - free (R); - free (dR); - free (D); - free (dD); - free (d2); - free (index); -} - -/** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ -// we only optionally apply the sigma limit: for object averages, this should not be used -int MeasFilterTestTiny(MeasureTiny *measure, int applySigmaLim) { - int found, k; - long mask; - PhotCode *code; - float mag; - - if (!finite(measure[0].R) || !finite(measure[0].D)) return FALSE; - if (!finite(measure[0].M)) return FALSE; //XXX is this necessary for all relastro tasks? - if (!finite(measure[0].dM)) return FALSE; //XXX is this necessary for all relastro tasks? - - if ((MinBadQF > 0.0) && (isGPC1chip(measure[0].photcode) || isGPC1stack(measure[0].photcode))) { - if (!isfinite(measure[0].psfQF)) { return FALSE; }; - if (measure[0].psfQF < MinBadQF) { return FALSE; }; - } - - /* select measurements by photcode, or equiv photcode, if specified */ - if (NphotcodesKeep > 0) { - found = FALSE; - for (k = 0; (k < NphotcodesKeep) && !found; k++) { - if (photcodesKeep[k][0].code == measure[0].photcode) found = TRUE; - if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE; - } - if (!found) return FALSE; - } - - if (NphotcodesSkip > 0) { - found = FALSE; - for (k = 0; (k < NphotcodesSkip) && !found; k++) { - if (photcodesSkip[k][0].code == measure[0].photcode) found = TRUE; - if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE; - } - if (found) return FALSE; - } - - /* select measurements by time */ - if (TimeSelect) { - if (measure[0].t < TSTART) return FALSE; - if (measure[0].t > TSTOP) return FALSE; - } - - /* select measurements by quality */ - if (PhotFlagSelect) { - if (PhotFlagBad) { - mask = PhotFlagBad; - } else { - code = GetPhotcodebyCode (measure[0].photcode); - if (!code) return FALSE; - mask = code[0].astromBadMask; - } - if (mask & measure[0].photFlags) return FALSE; - } - - /* select measurements by measurement error */ - // this is a bit convoluted: applySigmaLim is only TRUE when this function is - // called by bcatalog. for UpdateObjects, it is FALSE - if (applySigmaLim && (SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) { - return FALSE; - } - - /* select measurements by mag limit */ - if (ImagSelect) { - mag = PhotInstTiny (measure, MAG_CLASS_PSF); - if (mag < ImagMin || mag > ImagMax) return FALSE; - } - - return TRUE; -} - -# define SUPER_VERBOSE 0 - -/** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ -// we only optionally apply the sigma limit: for object averages, this should not be used (should it?) -int MeasFilterTest(Measure *measure, int applySigmaLim) { - int found, k; - long mask; - PhotCode *code; - float mag; - - if (!finite(measure[0].R) || !finite(measure[0].D)) { if (SUPER_VERBOSE) fprintf (stderr, "filter 1\n"); return FALSE; }; - if (!finite(measure[0].M)) { if (SUPER_VERBOSE) fprintf (stderr, "filter 2\n"); return FALSE; }; //XXX is this necessary for all relastro tasks? - if (!finite(measure[0].dM)) { if (SUPER_VERBOSE) fprintf (stderr, "filter 3\n"); return FALSE; }; //XXX is this necessary for all relastro tasks? - - /* select measurements by photcode, or equiv photcode, if specified */ - if (NphotcodesKeep > 0) { - found = FALSE; - for (k = 0; (k < NphotcodesKeep) && !found; k++) { - if (photcodesKeep[k][0].code == measure[0].photcode) found = TRUE; - if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE; - } - if (!found) { if (SUPER_VERBOSE) fprintf (stderr, "filter 4\n"); return FALSE; }; - } - - if (NphotcodesSkip > 0) { - found = FALSE; - for (k = 0; (k < NphotcodesSkip) && !found; k++) { - if (photcodesSkip[k][0].code == measure[0].photcode) found = TRUE; - if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE; - } - if (found) { if (SUPER_VERBOSE) fprintf (stderr, "filter 5\n"); return FALSE; }; - } - - if ((MinBadQF > 0.0) && isGPC1chip(measure[0].photcode)) { - if (!isfinite(measure[0].psfQF)) { if (SUPER_VERBOSE) fprintf (stderr, "filter 6\n"); return FALSE; }; - if (measure[0].psfQF < MinBadQF) { if (SUPER_VERBOSE) fprintf (stderr, "filter 7\n"); return FALSE; }; - } - - /* select measurements by time */ - if (TimeSelect) { - if (measure[0].t < TSTART) { if (SUPER_VERBOSE) fprintf (stderr, "filter 8\n"); return FALSE; }; - if (measure[0].t > TSTOP) { if (SUPER_VERBOSE) fprintf (stderr, "filter 9\n"); return FALSE; }; - } - - /* select measurements by quality */ - if (PhotFlagSelect) { - if (PhotFlagBad) { - mask = PhotFlagBad; - } else { - code = GetPhotcodebyCode (measure[0].photcode); - if (!code) { if (SUPER_VERBOSE) fprintf (stderr, "filter 10\n"); return FALSE; }; - mask = code[0].astromBadMask; - } - if (mask & measure[0].photFlags) { if (SUPER_VERBOSE) fprintf (stderr, "filter 11\n"); return FALSE; }; - } - - /* select measurements by measurement error */ - if (applySigmaLim && (SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) { - { if (SUPER_VERBOSE) fprintf (stderr, "filter 12\n"); return FALSE; }; - } - - /* select measurements by mag limit */ - if (ImagSelect) { - mag = PhotInst (measure, MAG_CLASS_PSF); - if (mag < ImagMin || mag > ImagMax) { if (SUPER_VERBOSE) fprintf (stderr, "filter 13\n"); return FALSE; }; - } - - return TRUE; -} -# endif Index: trunk/Ohana/src/relastro/src/StarMaps.c =================================================================== --- trunk/Ohana/src/relastro/src/StarMaps.c (revision 39693) +++ trunk/Ohana/src/relastro/src/StarMaps.c (revision 39926) @@ -166,4 +166,19 @@ dLmax = dMmax = 0.0; + float plateScale; + if (images[N].coords.mosaic) { + // NOTE: for the full pixel to sky plate scale, use this: + // float plateScaleX = 3600.0*images[N].coords.mosaic->cdelt1*images[N].coords.cdelt1; + // float plateScaleY = 3600.0*images[N].coords.mosaic->cdelt2*images[N].coords.cdelt2; + + // since we are compare L,M values, just need to compensate for focal plate to sky: + float plateScaleX = 3600.0*fabs(images[N].coords.mosaic->cdelt1); + float plateScaleY = 3600.0*fabs(images[N].coords.mosaic->cdelt2); + plateScale = 0.5*(plateScaleX + plateScaleY); + } else { + // since we are compare L,M values, just need to compensate for arcsec vs degrees: + plateScale = 3600.0; + } + for (i = 0; i < starmap[N].Npoints; i++) { @@ -171,6 +186,6 @@ XY_to_LM (&L, &M, starmap[N].points[i].X, starmap[N].points[i].Y, &images[N].coords); - starmap[N].points[i].dL = starmap[N].points[i].L - L; - starmap[N].points[i].dM = starmap[N].points[i].M - M; + starmap[N].points[i].dL = plateScale*(starmap[N].points[i].L - L); + starmap[N].points[i].dM = plateScale*(starmap[N].points[i].M - M); dLmax = MAX(fabs(starmap[N].points[i].dL), dLmax); Index: trunk/Ohana/src/relastro/src/UpdateChips.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateChips.c (revision 39693) +++ trunk/Ohana/src/relastro/src/UpdateChips.c (revision 39926) @@ -42,4 +42,8 @@ AstromErrorSetLoop (Nloop, TRUE); + // if ChipMapLoop or ChipOrderLoop is set use that to define the value of CHIPMAP and/or CHIPORDER this loop + if (ChipMapLoop) { CHIPMAP = ChipMapLoop[Nloop]; } + if (ChipOrderLoop) { CHIPORDER = ChipOrderLoop[Nloop]; } + if (NTHREADS) { UpdateChips_threaded (catalog, Ncatalog); @@ -60,9 +64,18 @@ for (i = 0; i < Nimage; i++) { - VERBOSE_IMAGE = !strcmp(image[i].name, "o5745g0516o.356887.cm.982631.smf[XY54]"); + VERBOSE = FALSE; + VERBOSE_IMAGE |= !strcmp(image[i].name, "o5745g0516o.356887.cm.982631.smf[XY45]"); + VERBOSE_IMAGE |= !strcmp(image[i].name, "o5745g0526o.356899.cm.982643.smf[XY45]"); + VERBOSE_IMAGE |= !strcmp(image[i].name, "o5748g0436o.358811.cm.982690.smf[XY34]"); // XXX looks like everything below is thread safe : we can unroll this into a set of // helper functions that grab the next available chip.... + // allow certain cameras to stay static + if (SKIP_PS1_CHIP && isGPC1chip (image[i].photcode)) { Nskip ++; mode[i] = 0; continue; } + if (SKIP_PS1_STACK && isGPC1stack(image[i].photcode)) { Nskip ++; mode[i] = 0; continue; } + if (SKIP_HSC && isHSCchip (image[i].photcode)) { Nskip ++; mode[i] = 0; continue; } + if (SKIP_CFH && isCFHchip (image[i].photcode)) { Nskip ++; mode[i] = 0; continue; } + /* skip all except WRP images */ if (strcmp(&image[i].coords.ctype[4], "-WRP")) { @@ -160,5 +173,5 @@ setImageRaw (catalog, Ncatalog, i, raw, Nraw, MODE_MOSAIC); if (USE_GALAXY_MODEL) { - // XXX DEPRECATE? + // the image calibration was calculated using a galaxy motion model image[i].flags |= ID_IMAGE_ASTROM_GMM; } @@ -281,4 +294,10 @@ } + // allow certain cameras to stay static + if (SKIP_PS1_CHIP && isGPC1chip (image[i].photcode)) { threadinfo->Nskip ++; threadinfo->mode[i] = 0; continue; } + if (SKIP_PS1_STACK && isGPC1stack(image[i].photcode)) { threadinfo->Nskip ++; threadinfo->mode[i] = 0; continue; } + if (SKIP_HSC && isHSCchip (image[i].photcode)) { threadinfo->Nskip ++; threadinfo->mode[i] = 0; continue; } + if (SKIP_CFH && isCFHchip (image[i].photcode)) { threadinfo->Nskip ++; threadinfo->mode[i] = 0; continue; } + /* skip all except WRP images */ if (strcmp(&image[i].coords.ctype[4], "-WRP")) { Index: trunk/Ohana/src/relastro/src/UpdateMeasures.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateMeasures.c (revision 39693) +++ trunk/Ohana/src/relastro/src/UpdateMeasures.c (revision 39926) @@ -115,5 +115,32 @@ // only modify the chip coordinates - if (isGPC1chip (measureT->photcode)) { + if (UPDATE_PS1_STACK_MEASURE && isGPC1stack (measureT->photcode)) { + measureT->R = R; + measureT->D = D; + if (measureB) { + measureB->R = R; + measureB->D = D; + } + } + // only modify the chip coordinates + if (UPDATE_PS1_CHIP_MEASURE && isGPC1chip (measureT->photcode)) { + measureT->R = R; + measureT->D = D; + if (measureB) { + measureB->R = R; + measureB->D = D; + } + } + // only modify the chip coordinates + if (UPDATE_HSC_MEASURE && isGPC1chip (measureT->photcode)) { + measureT->R = R; + measureT->D = D; + if (measureB) { + measureB->R = R; + measureB->D = D; + } + } + // only modify the chip coordinates + if (UPDATE_CFH_MEASURE && isGPC1chip (measureT->photcode)) { measureT->R = R; measureT->D = D; Index: trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 39693) +++ trunk/Ohana/src/relastro/src/UpdateObjectOffsets.c (revision 39926) @@ -167,10 +167,14 @@ strextend (&command, "relastro_client -update-offsets"); strextend (&command, "-hostID %d", group->hosts[i][0].hostID); + strextend (&command, "-hostdir %s", group->hosts[i][0].pathname); + strextend (&command, "-D CATDIR %s", CATDIR); - strextend (&command, "-hostdir %s", group->hosts[i][0].pathname); strextend (&command, "-region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); strextend (&command, "-statmode %s", STATMODE); strextend (&command, "-minerror %f", MIN_ERROR); + strextend (&command, "-D RELASTRO_SIGMA_LIM %f", SIGMA_LIM); + strextend (&command, "-D RELASTRO_SRC_MEAS_TOOFEW %d", SRC_MEAS_TOOFEW); + if (FIT_MODE == FIT_PM_ONLY) strextend (&command, "-pm"); if (FIT_MODE == FIT_PAR_ONLY) strextend (&command, "-par"); @@ -180,21 +184,11 @@ if (VERBOSE2) strextend (&command, "-vv"); if (RESET) strextend (&command, "-reset"); - if (UPDATE) strextend (&command, "-update"); - - if (RESET_BAD_IMAGES) strextend (&command, "-reset-bad-images"); if (ImagSelect) strextend (&command, "-instmag %f %f", ImagMin, ImagMax); if (MaxDensityUse) strextend (&command, "-max-density %f", MaxDensityValue); - - if (USE_BASIC_CHECK) strextend (&command, "-basic-image-search"); if (FlagOutlier) strextend (&command, "-clip %d", CLIP_THRESH); if (ExcludeBogus) strextend (&command, "-exclude-bogus %f", ExcludeBogusRadius); - if (USE_ALL_IMAGES) strextend (&command, "-use-all-images"); if (USE_FIXED_PIXCOORDS) strextend (&command, "-D USE_FIXED_PIXCOORDS 1"); - - if (REPAIR_STACKS) strextend (&command, "-repair-stacks-on-update"); - if (CHECK_MEASURE_TO_IMAGE) strextend (&command, "-check-measures"); - if (PHOTCODE_KEEP_LIST) strextend (&command, "+photcode %s", PHOTCODE_KEEP_LIST); if (PHOTCODE_SKIP_LIST) strextend (&command, "-photcode %s", PHOTCODE_SKIP_LIST); @@ -204,6 +198,4 @@ // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values - if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES); - if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) { strextend (&command, "-dcr-blue-color %s %s", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG); @@ -213,6 +205,20 @@ } + if (REPAIR_STACKS) strextend (&command, "-repair-stacks-on-update"); + if (CHECK_MEASURE_TO_IMAGE) strextend (&command, "-check-measures"); + + if (UPDATE_PS1_STACK_MEASURE) strextend (&command, "-update-ps1-stack"); + if (UPDATE_PS1_CHIP_MEASURE) strextend (&command, "-update-ps1-chip"); + if (UPDATE_HSC_MEASURE) strextend (&command, "-update-hsc"); + if (UPDATE_CFH_MEASURE) strextend (&command, "-update-cfh"); + + if (UPDATE) strextend (&command, "-update"); + if (RESET_BAD_IMAGES) strextend (&command, "-reset-bad-images"); + if (USE_BASIC_CHECK) strextend (&command, "-basic-image-search"); + if (USE_ALL_IMAGES) strextend (&command, "-use-all-images"); + if (MinBadQF > 0.0) strextend (&command, "-min-bad-psfqf %f", MinBadQF); if (MaxMeanOffset != 10.0) strextend (&command, "-max-mean-offset %f", MaxMeanOffset); + if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES); if (TimeSelect) { @@ -300,10 +306,14 @@ strextend (&command, "relastro_client -update-offsets"); strextend (&command, "-hostID %d", table->hosts[i].hostID); + strextend (&command, "-hostdir %s", table->hosts[i].pathname); + strextend (&command, "-D CATDIR %s", CATDIR); - strextend (&command, "-hostdir %s", table->hosts[i].pathname); strextend (&command, "-region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); strextend (&command, "-statmode %s", STATMODE); strextend (&command, "-minerror %f", MIN_ERROR); + strextend (&command, "-D RELASTRO_SIGMA_LIM %f", SIGMA_LIM); + strextend (&command, "-D RELASTRO_SRC_MEAS_TOOFEW %d", SRC_MEAS_TOOFEW); + if (FIT_MODE == FIT_PM_ONLY) strextend (&command, "-pm"); if (FIT_MODE == FIT_PAR_ONLY) strextend (&command, "-par"); @@ -313,21 +323,11 @@ if (VERBOSE2) strextend (&command, "-vv"); if (RESET) strextend (&command, "-reset"); - if (UPDATE) strextend (&command, "-update"); - - if (RESET_BAD_IMAGES) strextend (&command, "-reset-bad-images"); if (ImagSelect) strextend (&command, "-instmag %f %f", ImagMin, ImagMax); if (MaxDensityUse) strextend (&command, "-max-density %f", MaxDensityValue); - - if (USE_BASIC_CHECK) strextend (&command, "-basic-image-search"); if (FlagOutlier) strextend (&command, "-clip %d", CLIP_THRESH); if (ExcludeBogus) strextend (&command, "-exclude-bogus %f", ExcludeBogusRadius); - if (USE_ALL_IMAGES) strextend (&command, "-use-all-images"); if (USE_FIXED_PIXCOORDS) strextend (&command, "-D USE_FIXED_PIXCOORDS 1"); - - if (REPAIR_STACKS) strextend (&command, "-repair-stacks-on-update"); - if (CHECK_MEASURE_TO_IMAGE) strextend (&command, "-check-measures"); - if (PHOTCODE_KEEP_LIST) strextend (&command, "+photcode %s", PHOTCODE_KEEP_LIST); if (PHOTCODE_SKIP_LIST) strextend (&command, "-photcode %s", PHOTCODE_SKIP_LIST); @@ -337,6 +337,4 @@ // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values - if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES); - if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) { strextend (&command, "-dcr-blue-color %s %s", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG); @@ -345,4 +343,19 @@ strextend (&command, "-dcr-red-color %s %s", DCR_RED_COLOR_POS, DCR_RED_COLOR_NEG); } + + if (REPAIR_STACKS) strextend (&command, "-repair-stacks-on-update"); + if (CHECK_MEASURE_TO_IMAGE) strextend (&command, "-check-measures"); + + if (UPDATE_PS1_STACK_MEASURE) strextend (&command, "-update-ps1-stack"); + if (UPDATE_PS1_CHIP_MEASURE) strextend (&command, "-update-ps1-chip"); + if (UPDATE_HSC_MEASURE) strextend (&command, "-update-hsc"); + if (UPDATE_CFH_MEASURE) strextend (&command, "-update-cfh"); + + if (UPDATE) strextend (&command, "-update"); + if (RESET_BAD_IMAGES) strextend (&command, "-reset-bad-images"); + if (USE_BASIC_CHECK) strextend (&command, "-basic-image-search"); + if (USE_ALL_IMAGES) strextend (&command, "-use-all-images"); + + if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES); if (MinBadQF > 0.0) strextend (&command, "-min-bad-psfqf %f", MinBadQF); Index: trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 39693) +++ trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 39926) @@ -5,13 +5,15 @@ # define PAR_MIN_NPTS 7 # define PAR_MIN_NPTS_BOOT 6 -# define PM_MIN_NPTS 4 -# define PM_MIN_NPTS_BOOT 3 -# define POS_MIN_NPTS 3 -# define POS_MIN_NPTS_BOOT 2 +# define PM_MIN_NPTS 5 +# define PM_MIN_NPTS_BOOT 4 +# define POS_MIN_NPTS 4 +# define POS_MIN_NPTS_BOOT 4 typedef enum { - SELECT_MEAS_HAS_DATA = 1, - SELECT_MEAS_HAS_STACK = 2, - SELECT_MEAS_HAS_2MASS = 4, + SELECT_MEAS_HAS_DATA = 0x01, + SELECT_MEAS_HAS_STACK = 0x02, + SELECT_MEAS_HAS_2MASS = 0x04, + SELECT_MEAS_HAS_GAIA = 0x08, + SELECT_MEAS_HAS_TYCHO = 0x10, } SelectMeasureStatus; @@ -181,34 +183,39 @@ if (Trange < PM_DT_MIN) { // not enough baseline for proper motion, only set mean position - mode = FIT_AVERAGE; goto justPosition; } if (parRange < PAR_FACTOR_MIN) { // not enough parallax factor range, skip parallax - mode = FIT_PM_ONLY; goto skipParallax; } if (fitStats->Npoints < PAR_MIN_NPTS) { // not enough data, skip parallax - 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; - if (!FitPMandPar_IRLS (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints)) { - mode = FIT_PM_ONLY; - goto skipParallax; - } - - if (N_BOOTSTRAP_SAMPLES) { + // then run N_BOOTSTRAP_SAMPLES to measure the errors. We first measure the OLS error, + // and choose the max of the OLS and bootstrap errors + fitStats->fitdataPar->getError = TRUE; + if (USE_IRLS) { + if (!FitPMandPar_IRLS (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints)) { + goto skipParallax; + } + } else { + if (!FitPMandPar_Basic (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints)) { + goto skipParallax; + } + } + + // in the fits above, we have saved the formal error for the unmasked points. + // if we do not have enough points for bootstrap, we will keep those errors + if (N_BOOTSTRAP_SAMPLES && (fitPar.Nfit >= PAR_MIN_NPTS_BOOT)) { fitStats->Nfit = 0; int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); - if (Nnomask < PAR_MIN_NPTS_BOOT) { - // if we do not have enough points to assess parallax error, we cannot do the bootstrap analysis. - mode = FIT_PM_ONLY; - goto skipParallax; - } + + // if we do not have enough points to assess parallax error, skip the bootstrap analysis. + // (I think the test above means we never do this skip) + if (Nnomask < PAR_MIN_NPTS_BOOT) goto skipParallaxBootstrap; + for (k = 0; k < fitStats->NfitAlloc; k++) { BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask); @@ -216,5 +223,8 @@ fitStats->Nfit ++; } - // these calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above) + + // These calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above) + // this call expects the fitted parameters to have the formal error set: it will apply the + // max of the formal and bootstrap errors BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); @@ -224,4 +234,5 @@ } + skipParallaxBootstrap: // project Ro, Do back to RA,DEC XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &fitStats->coords); @@ -242,7 +253,5 @@ valid = valid && (fabs(fitPar.uD) < 4.0); valid = valid && (fabs(fitPar.p) < 2.0); - if (!valid) { - mode = FIT_PM_ONLY; - } else { + if (valid) { average[0].flags |= ID_OBJ_USE_PAR; } @@ -254,44 +263,55 @@ if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) { if (Trange < PM_DT_MIN) { - mode = FIT_AVERAGE; goto justPosition; } if (fitStats->Npoints < PM_MIN_NPTS) { - mode = FIT_AVERAGE; goto justPosition; } // if we have fitted (and accepted) a parallax model, get the best pm fit and chisq - // given the set of points (mask is respected) - if (average[0].flags & ID_OBJ_USE_PAR) { + // given the set of points (mask is respected). Alternatively, if we do not request + // IRLS fitting, the just use OLS fitting (skip bootstrap) + if ((average[0].flags & ID_OBJ_USE_PAR) || !USE_IRLS) { if (!FitPM_Basic (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) { average[0].flags |= ID_OBJ_BAD_PM; goto justPosition; } - } else { - fitStats->fitdataPM->getError = !N_BOOTSTRAP_SAMPLES; - if (!FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) { - mode = FIT_AVERAGE; - goto justPosition; - } - if (N_BOOTSTRAP_SAMPLES) { - fitStats->Nfit = 0; - int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); - if (Nnomask < PM_MIN_NPTS_BOOT) { - mode = FIT_AVERAGE; - goto justPosition; - } - for (k = 0; k < fitStats->NfitAlloc; k++) { - BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask); - if (!FitPM_Basic (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, Nnomask)) 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); - } - } - + goto skipProperMotionBootstrap; + } + + // 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. We first measure the OLS error, + // and choose the max of the OLS and bootstrap errors + fitStats->fitdataPM->getError = TRUE; + if (!FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) { + goto justPosition; + } + + // in the fits above, we have saved the formal error for the unmasked points. + // if we do not have enough points for bootstrap, we will keep those errors + if (N_BOOTSTRAP_SAMPLES && (fitPM.Nfit >= PM_MIN_NPTS_BOOT)) { + fitStats->Nfit = 0; + int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); + + // if we do not have enough points to assess p.m. error, skip the bootstrap analysis. + // (I think the test above means we never do this skip) + if (Nnomask < PM_MIN_NPTS_BOOT) goto skipProperMotionBootstrap; + + for (k = 0; k < fitStats->NfitAlloc; k++) { + BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask); + if (!FitPM_Basic (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, Nnomask)) continue; + fitStats->Nfit ++; + } + + // These calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above) + // this call expects the fitted parameters to have the formal error set: it will apply the + // max of the formal and bootstrap errors + 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); + } + +skipProperMotionBootstrap: // project Ro, Do back to RA,DEC XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &fitStats->coords); @@ -300,4 +320,5 @@ fitStats->Npm ++; + // XXX a hard-wired hack... // unless there is a clear problems (below) with the proper-motion fit or we have a parallax fit, we will use pm fit int valid = TRUE; @@ -309,5 +330,4 @@ valid = valid && (fabs(fitPM.uD) < 4.0); if (!valid) { - mode = FIT_AVERAGE; average[0].flags |= ID_OBJ_BAD_PM; } else { @@ -323,59 +343,56 @@ // if we only have one point, this is silly... + // set the proper motion (to the galaxy model or average value, if desired; else to 0,0) FitAstromResultSetPM (&fitPos, 1, average); - // fprintf (stderr, "fit 1: %f %f : %f %f\n", fitPos.Ro, fitPos.Do, fitPos.uR, fitPos.uD); - if (average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_PM)) { + + // if we already have a valid fit (pm or par), use OLS to fit the position + // alternatively, if we do not request IRLS, use OLS + // alternatively, if we do not have enough points, use OLS + if ((average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_PM)) || (fitStats->Npoints < POS_MIN_NPTS) || !USE_IRLS) { if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { - // if this fails, stick with the PM and/or PAR fit from above + // if this fails, stick with the PM and/or PAR fit from above, or use a single value goto doneWithFit; } - } else { - if (fitStats->Npoints < POS_MIN_NPTS) { - // I will not try to outlier-reject, just calculate the weighted average - if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { - // if we have tried this, we have not masked any points; we will find a single unmasked point below - goto doneWithFit; - } + // if we have not already gotten a good fit, use this fit + if (!(average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_PM))) { average[0].flags |= ID_OBJ_USE_AVE; average[0].flags |= ID_OBJ_RAW_AVE; - goto useBasic; - } - if (!FitPosPMfixed_IRLS (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { - // if the above fails, we need to clear the masks and try again below - FitPointsClearMasks (fitStats->points, fitStats->Npoints); - // just calculate the weighted average - if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { - // if this fails, find a single unmasked point below - goto doneWithFit; - } - average[0].flags |= ID_OBJ_USE_AVE; - average[0].flags |= ID_OBJ_RAW_AVE; - goto useBasic; - } - // fprintf (stderr, "fit 2: %f %f : %f %f\n", fitPos.Ro, fitPos.Do, fitPos.uR, fitPos.uD); - if (N_BOOTSTRAP_SAMPLES) { - fitStats->Nfit = 0; - int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); - if (Nnomask < POS_MIN_NPTS_BOOT) { - // if the above fails, we need to clear the masks and try again below - FitPointsClearMasks (fitStats->points, fitStats->Npoints); - if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { - // if this fails, find a single unmasked point below - goto doneWithFit; - } - average[0].flags |= ID_OBJ_USE_AVE; - average[0].flags |= ID_OBJ_RAW_AVE; - goto useBasic; - } - FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average); - for (k = 0; k < fitStats->NfitAlloc; k++) { - BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask); - if (!FitPosPMfixed_Basic (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, Nnomask)) continue; - fitStats->Nfit ++; - } - BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); - BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); + } + goto useBasic; + } + + // try the IRLS fitting, otherwise give up and use OLS + if (!FitPosPMfixed_IRLS (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { + // if the above fails, we need to clear the masks and try again below + FitPointsClearMasks (fitStats->points, fitStats->Npoints); + // just calculate the weighted average + if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { + // if this fails, find a single unmasked point below + goto doneWithFit; } average[0].flags |= ID_OBJ_USE_AVE; + average[0].flags |= ID_OBJ_RAW_AVE; + goto useBasic; + } + average[0].flags |= ID_OBJ_USE_AVE; + + if (N_BOOTSTRAP_SAMPLES && (fitPos.Nfit >= POS_MIN_NPTS_BOOT)) { + fitStats->Nfit = 0; + int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); + + if (Nnomask < POS_MIN_NPTS_BOOT) goto useBasic; + + FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average); + for (k = 0; k < fitStats->NfitAlloc; k++) { + BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask); + if (!FitPosPMfixed_Basic (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, Nnomask)) continue; + fitStats->Nfit ++; + } + + // These calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above) + // this call expects the fitted parameters to have the formal error set: it will apply the + // max of the formal and bootstrap errors + BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA); + BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC); } @@ -544,4 +561,5 @@ average[0].dP = fit.dp; // parallax error in arcsec + average[0].ChiSqAve = fitPos.chisq; average[0].ChiSqPM = fitPM.chisq; @@ -666,4 +684,6 @@ int has2MASS = FALSE; + int hasGAIA = FALSE; + int hasTycho = FALSE; int hasStack = FALSE; if (stackEntry) *stackEntry = -1; @@ -774,7 +794,7 @@ Npoints++; - if ((measure[k].photcode >= 2011) && (measure[k].photcode <= 2013)) { - has2MASS = TRUE; - } + hasGAIA = (measure[k].photcode == 1030); + has2MASS = ((measure[k].photcode >= 2011) && (measure[k].photcode <= 2013)); + hasTycho = ((measure[k].photcode >= 2020) && (measure[k].photcode <= 2021)); myAssert (Npoints <= fit->NpointsAlloc, "oops"); @@ -796,4 +816,14 @@ measure[k].dbFlags &= ~ID_MEAS_OBJECT_HAS_2MASS; } + if (hasGAIA) { + measure[k].dbFlags |= ID_MEAS_OBJECT_HAS_GAIA; + } else { + measure[k].dbFlags &= ~ID_MEAS_OBJECT_HAS_GAIA; + } + if (hasTycho) { + measure[k].dbFlags |= ID_MEAS_OBJECT_HAS_TYCHO; + } else { + measure[k].dbFlags &= ~ID_MEAS_OBJECT_HAS_TYCHO; + } } @@ -803,4 +833,6 @@ if (hasStack) status |= SELECT_MEAS_HAS_STACK; if (has2MASS) status |= SELECT_MEAS_HAS_2MASS; + if (hasGAIA) status |= SELECT_MEAS_HAS_GAIA; + if (hasTycho) status |= SELECT_MEAS_HAS_TYCHO; return status; } @@ -822,4 +854,6 @@ int i; + + myAbort ("this should not be called anymore"); // add up the chi square for the fit @@ -852,8 +886,13 @@ int i; - if (USE_GALAXY_MODEL) { + if (APPLY_PROPER_MOTION) { for (i = 0; i < Nfit; i++) { - fit[i].uR = average->uRgal; - fit[i].uD = average->uDgal; + if (USE_GALAXY_MODEL) { + fit[i].uR = average->uRgal; + fit[i].uD = average->uDgal; + } else { + fit[i].uR = average->uR; + fit[i].uD = average->uD; + } } } else { Index: trunk/Ohana/src/relastro/src/UpdateStacks.c =================================================================== --- trunk/Ohana/src/relastro/src/UpdateStacks.c (revision 39693) +++ trunk/Ohana/src/relastro/src/UpdateStacks.c (revision 39926) @@ -3,4 +3,6 @@ // NOTE: we only measure the systematic floor of the astrometric scatter per stack, no change to the calibration int UpdateStacks (Catalog *catalog, int Ncatalog) { + + if (SKIP_PS1_STACK) return TRUE; off_t Nimage; @@ -46,6 +48,6 @@ // XXX: I need to convert dLsig, dMsig from degrees to pixels - dLsig *= 3600.0 / 0.25; - dMsig *= 3600.0 / 0.25; + dLsig *= 3600.0; + dMsig *= 3600.0; image[i].dXpixSys = dLsig; Index: trunk/Ohana/src/relastro/src/args.c =================================================================== --- trunk/Ohana/src/relastro/src/args.c (revision 39693) +++ trunk/Ohana/src/relastro/src/args.c (revision 39926) @@ -1,8 +1,9 @@ # include "relastro.h" -void usage (void); -void usage_client (void); +void usage (int argc, char **argv); +void usage_client (int argc, char **argv); void usage_merge_source (void); void usage_merge_source_id (char *name); float *ParseLoopWeights (char *rawlist); +int *ParseLoopOrder (char *rawlist, int minValue); int args (int argc, char **argv) { @@ -39,5 +40,5 @@ remove_argument (N, &argc, argv); - if (argc != 1) usage (); + if (argc != 1) usage (argc, argv); return TRUE; } @@ -76,4 +77,11 @@ remove_argument (N, &argc, argv); CHECK_MEASURE_TO_IMAGE = TRUE; + } + + // catch-up mode : for a re-run, allow the sync file to be ahead of the desired location: + CATCH_UP = FALSE; + if ((N = get_argument (argc, argv, "-catch-up"))) { + remove_argument (N, &argc, argv); + CATCH_UP = TRUE; } @@ -105,8 +113,8 @@ remove_argument (N, &argc, argv); RELASTRO_OP = OP_PARALLEL_IMAGES; - if (N >= argc) usage(); + if (N >= argc) usage (argc, argv); IMAGE_TABLE = strcreate (argv[N]); remove_argument (N, &argc, argv); - if (!REGION_FILE) usage(); + if (!REGION_FILE) usage (argc, argv); } @@ -114,4 +122,81 @@ remove_argument (N, &argc, argv); RELASTRO_OP = OP_IMAGES; + } + + // used to decide if changes to the image parameters get applied to the measures + APPLY_OFFSETS = FALSE; + if ((N = get_argument (argc, argv, "-apply-offsets"))) { + remove_argument (N, &argc, argv); + APPLY_OFFSETS = TRUE; + } + + APPLY_PROPER_MOTION = FALSE; + if ((N = get_argument (argc, argv, "-apply-proper-motion"))) { + remove_argument (N, &argc, argv); + APPLY_PROPER_MOTION = TRUE; + } + + SKIP_PS1_CHIP = FALSE; + if ((N = get_argument (argc, argv, "-skip-ps1-chip"))) { + remove_argument (N, &argc, argv); + SKIP_PS1_CHIP = TRUE; + } + SKIP_PS1_STACK = FALSE; + if ((N = get_argument (argc, argv, "-skip-ps1-stack"))) { + remove_argument (N, &argc, argv); + SKIP_PS1_STACK = TRUE; + } + SKIP_HSC = FALSE; + if ((N = get_argument (argc, argv, "-skip-hsc"))) { + remove_argument (N, &argc, argv); + SKIP_HSC = TRUE; + } + SKIP_CFH = FALSE; + if ((N = get_argument (argc, argv, "-skip-cfh"))) { + remove_argument (N, &argc, argv); + SKIP_CFH = TRUE; + } + + UPDATE_PS1_STACK_MEASURE = FALSE; + if ((N = get_argument (argc, argv, "-update-ps1-stack"))) { + remove_argument (N, &argc, argv); + UPDATE_PS1_STACK_MEASURE = TRUE; + } + UPDATE_PS1_CHIP_MEASURE = FALSE; + if ((N = get_argument (argc, argv, "-update-ps1-chip"))) { + remove_argument (N, &argc, argv); + UPDATE_PS1_CHIP_MEASURE = TRUE; + } + UPDATE_HSC_MEASURE = FALSE; + if ((N = get_argument (argc, argv, "-update-hsc"))) { + remove_argument (N, &argc, argv); + UPDATE_HSC_MEASURE = TRUE; + } + UPDATE_CFH_MEASURE = FALSE; + if ((N = get_argument (argc, argv, "-update-cfh"))) { + remove_argument (N, &argc, argv); + UPDATE_CFH_MEASURE = TRUE; + } + if (RELASTRO_OP == OP_UPDATE_OFFSETS) { + if (!UPDATE_PS1_STACK_MEASURE && !UPDATE_PS1_CHIP_MEASURE && !UPDATE_HSC_MEASURE && !UPDATE_CFH_MEASURE) { + fprintf (stderr, "for -update-offsets, need to select at least one of -update-ps1-stack, -update-ps1-chip, -update-hsc, -update-cfh\n"); + exit (2); + } + } + if ((RELASTRO_OP == OP_PARALLEL_IMAGES) || (RELASTRO_OP == OP_IMAGES)) { + if (APPLY_OFFSETS && !UPDATE_PS1_STACK_MEASURE && !UPDATE_PS1_CHIP_MEASURE && !UPDATE_HSC_MEASURE && !UPDATE_CFH_MEASURE) { + fprintf (stderr, "for [-images or -parallel-images] with -apply-offsets, need to select at least one of -update-ps1-stack, -update-ps1-chip, -update-hsc, -update-cfh\n"); + exit (2); + } + } + + // for fitting objects, this is always the same as 'ALLOW_IRLS' below, but for fitting + // images, this is set to FALSE while doing the fit for the image parameters + USE_IRLS = TRUE; + ALLOW_IRLS = TRUE; + if ((N = get_argument (argc, argv, "-no-irls"))) { + remove_argument (N, &argc, argv); + USE_IRLS = FALSE; + ALLOW_IRLS = FALSE; } @@ -120,5 +205,5 @@ remove_argument (N, &argc, argv); RELASTRO_OP = OP_PARALLEL_REGIONS; - if (!REGION_FILE) usage(); + if (!REGION_FILE) usage (argc, argv); if ((N = get_argument (argc, argv, "-parallel-regions-manual"))) { remove_argument (N, &argc, argv); @@ -128,22 +213,22 @@ if ((N = get_argument (argc, argv, "-testobj1"))) { - if (N > argc - 3) usage (); + if (N > argc - 3) usage (argc, argv); remove_argument (N, &argc, argv); OBJ_ID_SRC = strtol(argv[N], &endptr, 0); - if (*endptr) usage (); + if (*endptr) usage (argc, argv); remove_argument (N, &argc, argv); CAT_ID_SRC = strtol(argv[N], &endptr, 0); - if (*endptr) usage (); + if (*endptr) usage (argc, argv); remove_argument (N, &argc, argv); } if ((N = get_argument (argc, argv, "-testobj2"))) { - if (N > argc - 3) usage (); + if (N > argc - 3) usage (argc, argv); remove_argument (N, &argc, argv); OBJ_ID_DST = strtol(argv[N], &endptr, 0); - if (*endptr) usage (); + if (*endptr) usage (argc, argv); remove_argument (N, &argc, argv); CAT_ID_DST = strtol(argv[N], &endptr, 0); - if (*endptr) usage (); + if (*endptr) usage (argc, argv); remove_argument (N, &argc, argv); } @@ -167,5 +252,5 @@ if ((N = get_argument (argc, argv, "-high-speed"))) { // XXX include a parallax / no-parallax option - if (N >= argc - 4) usage(); + if (N >= argc - 4) usage (argc, argv); RELASTRO_OP = OP_HIGH_SPEED; remove_argument (N, &argc, argv); @@ -181,5 +266,5 @@ if ((N = get_argument (argc, argv, "-hpm"))) { - if (N >= argc - 2) usage(); + if (N >= argc - 2) usage (argc, argv); RELASTRO_OP = OP_HPM; remove_argument (N, &argc, argv); @@ -212,4 +297,8 @@ FIT_TARGET = TARGET_MOSAICS; } + if ((N = get_argument (argc, argv, "-set-chips"))) { + remove_argument (N, &argc, argv); + FIT_TARGET = SET_CHIPS; + } FlagOutlier = FALSE; @@ -231,7 +320,7 @@ } - if (RELASTRO_OP == OP_NONE) usage(); - - if (((RELASTRO_OP == OP_IMAGES) || (RELASTRO_OP == OP_PARALLEL_REGIONS) || (RELASTRO_OP == OP_PARALLEL_IMAGES)) && (FIT_TARGET == TARGET_NONE)) usage(); + if (RELASTRO_OP == OP_NONE) usage (argc, argv); + + if (((RELASTRO_OP == OP_IMAGES) || (RELASTRO_OP == OP_PARALLEL_REGIONS) || (RELASTRO_OP == OP_PARALLEL_IMAGES)) && (FIT_TARGET == TARGET_NONE)) usage (argc, argv); /* specify portion of the sky : allow default of all sky? */ @@ -303,4 +392,10 @@ } + KEEP_ALL_IMAGES_RA = FALSE; + if ((N = get_argument (argc, argv, "-keep-all-images-ra"))) { + remove_argument (N, &argc, argv); + KEEP_ALL_IMAGES_RA = TRUE; + } + USE_BASIC_CHECK = FALSE; if ((N = get_argument (argc, argv, "-basic-image-search"))) { @@ -315,10 +410,4 @@ remove_argument (N, &argc, argv); MaxDensityUse = TRUE; - } - - APPLY_OFFSETS = FALSE; - if ((N = get_argument (argc, argv, "-apply-offsets"))) { - remove_argument (N, &argc, argv); - APPLY_OFFSETS = TRUE; } @@ -439,19 +528,4 @@ } - CHIPORDER = 0; - if ((N = get_argument (argc, argv, "-chiporder"))) { - remove_argument (N, &argc, argv); - CHIPORDER = atoi(argv[N]); - remove_argument (N, &argc, argv); - } - - CHIPMAP = 0; - if ((N = get_argument (argc, argv, "-chipmap"))) { - remove_argument (N, &argc, argv); - CHIPMAP = atoi(argv[N]); - remove_argument (N, &argc, argv); - - } - SAVEPLOT = FALSE; PLOTSTUFF = FALSE; @@ -580,4 +654,38 @@ } + // e.g., -chiporderloop 3,4,5 + // NOTE: this must come after -nloop above + ChipOrderLoop = NULL; + ChipOrderLoopStr = NULL; + CHIPORDER = 1; + if ((N = get_argument (argc, argv, "-chiporder"))) { + remove_argument (N, &argc, argv); + CHIPORDER = atoi(argv[N]); + remove_argument (N, &argc, argv); + } + if ((N = get_argument (argc, argv, "-chiporderloop"))) { + remove_argument (N, &argc, argv); + ChipOrderLoopStr = strcreate(argv[N]); + ChipOrderLoop = ParseLoopOrder (argv[N], 1); + remove_argument (N, &argc, argv); + } + + ChipMapLoop = NULL; + ChipMapLoopStr = NULL; + CHIPMAP = 0; + if ((N = get_argument (argc, argv, "-chipmap"))) { + remove_argument (N, &argc, argv); + CHIPMAP = atoi(argv[N]); + remove_argument (N, &argc, argv); + + } + if ((N = get_argument (argc, argv, "-chipmaploop"))) { + remove_argument (N, &argc, argv); + ChipMapLoopStr = strcreate(argv[N]); + ChipMapLoop = ParseLoopOrder (argv[N], 0); + remove_argument (N, &argc, argv); + } + + // e.g., -loop-weights-2mass 1000,300,300,200,200,100 // NOTE: this must come after -nloop above @@ -591,4 +699,5 @@ } LoopWeightTycho = NULL; + LoopWeightTychostr = NULL; if ((N = get_argument (argc, argv, "-loop-weights-tycho"))) { remove_argument (N, &argc, argv); @@ -597,4 +706,28 @@ remove_argument (N, &argc, argv); } + LoopWeightGAIA = NULL; + LoopWeightGAIAstr = NULL; + if ((N = get_argument (argc, argv, "-loop-weights-gaia"))) { + remove_argument (N, &argc, argv); + LoopWeightGAIAstr = strcreate(argv[N]); + LoopWeightGAIA = ParseLoopWeights (argv[N]); + remove_argument (N, &argc, argv); + } + + GALAXY_MODEL = NULL; + if ((N = get_argument (argc, argv, "-galaxy-model"))) { + remove_argument (N, &argc, argv); + GALAXY_MODEL = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + if (!GALAXY_MODEL) GALAXY_MODEL = strcreate ("FEAST-HIPPARCOS"); + + // for testing, allow the galaxy model scale to be non-unity + TEST_SCALE = 1.0; + if ((N = get_argument (argc, argv, "-testing"))) { + remove_argument (N, &argc, argv); + TEST_SCALE = atof(argv[N]); + remove_argument (N, &argc, argv); + } NTHREADS = 0; @@ -605,5 +738,5 @@ } - if (argc != 1) usage (); + if (argc != 1) usage (argc, argv); return TRUE; } @@ -614,4 +747,13 @@ FREE (LoopWeight2MASSstr); FREE (LoopWeightTychostr); + FREE (LoopWeightGAIAstr); + FREE (LoopWeight2MASS); + FREE (LoopWeightTycho); + FREE (LoopWeightGAIA); + + FREE (ChipMapLoop); + FREE (ChipMapLoopStr); + FREE (ChipOrderLoop); + FREE (ChipOrderLoopStr); FREE (PHOTCODE_SKIP_LIST); @@ -628,4 +770,5 @@ FREE (BCATALOG); FREE (HOSTDIR); + FREE (GALAXY_MODEL); // these are set in initialize @@ -673,5 +816,5 @@ remove_argument (N, &argc, argv); } - if (!HOST_ID) usage_client(); + if (!HOST_ID) usage_client (argc, argv); HOSTDIR = NULL; @@ -681,5 +824,5 @@ remove_argument (N, &argc, argv); } - if (!HOSTDIR) usage_client(); + if (!HOSTDIR) usage_client (argc, argv); if ((N = get_argument (argc, argv, "-load-objects"))) { @@ -726,4 +869,52 @@ } + SKIP_PS1_CHIP = FALSE; + if ((N = get_argument (argc, argv, "-skip-ps1-chip"))) { + remove_argument (N, &argc, argv); + SKIP_PS1_CHIP = TRUE; + } + SKIP_PS1_STACK = FALSE; + if ((N = get_argument (argc, argv, "-skip-ps1-stack"))) { + remove_argument (N, &argc, argv); + SKIP_PS1_STACK = TRUE; + } + SKIP_HSC = FALSE; + if ((N = get_argument (argc, argv, "-skip-hsc"))) { + remove_argument (N, &argc, argv); + SKIP_HSC = TRUE; + } + SKIP_CFH = FALSE; + if ((N = get_argument (argc, argv, "-skip-cfh"))) { + remove_argument (N, &argc, argv); + SKIP_CFH = TRUE; + } + + UPDATE_PS1_STACK_MEASURE = FALSE; + if ((N = get_argument (argc, argv, "-update-ps1-stack"))) { + remove_argument (N, &argc, argv); + UPDATE_PS1_STACK_MEASURE = TRUE; + } + UPDATE_PS1_CHIP_MEASURE = FALSE; + if ((N = get_argument (argc, argv, "-update-ps1-chip"))) { + remove_argument (N, &argc, argv); + UPDATE_PS1_CHIP_MEASURE = TRUE; + } + UPDATE_HSC_MEASURE = FALSE; + if ((N = get_argument (argc, argv, "-update-hsc"))) { + remove_argument (N, &argc, argv); + UPDATE_HSC_MEASURE = TRUE; + } + UPDATE_CFH_MEASURE = FALSE; + if ((N = get_argument (argc, argv, "-update-cfh"))) { + remove_argument (N, &argc, argv); + UPDATE_CFH_MEASURE = TRUE; + } + if (RELASTRO_OP == OP_UPDATE_OFFSETS) { + if (!UPDATE_PS1_STACK_MEASURE && !UPDATE_PS1_CHIP_MEASURE && !UPDATE_HSC_MEASURE && !UPDATE_CFH_MEASURE) { + fprintf (stderr, "for -update-offsets, need to select at least one of -update-ps1-stack, -update-ps1-chip, -update-hsc, -update-cfh\n"); + exit (2); + } + } + // check for object fitting modes if ((N = get_argument (argc, argv, "-pm"))) { @@ -742,5 +933,5 @@ if ((N = get_argument (argc, argv, "-high-speed"))) { // XXX include a parallax / no-parallax option - if (N >= argc - 5) usage_client(); + if (N >= argc - 5) usage_client (argc, argv); RELASTRO_OP = OP_HIGH_SPEED; remove_argument (N, &argc, argv); @@ -756,5 +947,5 @@ if ((N = get_argument (argc, argv, "-hpm"))) { - if (N >= argc - 3) usage(); + if (N >= argc - 3) usage_client (argc, argv); RELASTRO_OP = OP_HPM; remove_argument (N, &argc, argv); @@ -766,22 +957,22 @@ if ((N = get_argument (argc, argv, "-testobj1"))) { - if (N > argc - 3) usage (); + if (N > argc - 3) usage_client (argc, argv); remove_argument (N, &argc, argv); OBJ_ID_SRC = strtol(argv[N], &endptr, 0); - if (*endptr) usage (); + if (*endptr) usage_client (argc, argv); remove_argument (N, &argc, argv); CAT_ID_SRC = strtol(argv[N], &endptr, 0); - if (*endptr) usage (); + if (*endptr) usage_client (argc, argv); remove_argument (N, &argc, argv); } if ((N = get_argument (argc, argv, "-testobj2"))) { - if (N > argc - 3) usage (); + if (N > argc - 3) usage_client (argc, argv); remove_argument (N, &argc, argv); OBJ_ID_DST = strtol(argv[N], &endptr, 0); - if (*endptr) usage (); + if (*endptr) usage_client (argc, argv); remove_argument (N, &argc, argv); CAT_ID_DST = strtol(argv[N], &endptr, 0); - if (*endptr) usage (); + if (*endptr) usage_client (argc, argv); remove_argument (N, &argc, argv); } @@ -805,5 +996,21 @@ } - if (RELASTRO_OP == OP_NONE) usage_client(); + GALAXY_MODEL = NULL; + if ((N = get_argument (argc, argv, "-galaxy-model"))) { + remove_argument (N, &argc, argv); + GALAXY_MODEL = strcreate(argv[N]); + remove_argument (N, &argc, argv); + } + if (!GALAXY_MODEL) GALAXY_MODEL = strcreate ("FEAST-HIPPARCOS"); + + // for testing, allow the galaxy model scale to be non-unity + TEST_SCALE = 1.0; + if ((N = get_argument (argc, argv, "-testing"))) { + remove_argument (N, &argc, argv); + TEST_SCALE = atof(argv[N]); + remove_argument (N, &argc, argv); + } + + if (RELASTRO_OP == OP_NONE) usage_client (argc, argv); /* specify portion of the sky : allow default of all sky? */ @@ -1032,5 +1239,5 @@ } - if (argc != 1) usage_client (); + if (argc != 1) usage_client (argc, argv); return TRUE; } @@ -1049,5 +1256,6 @@ FREE(HIGH_SPEED_DIR); FREE(BCATALOG); - FREE (HOSTDIR); + FREE(HOSTDIR); + FREE(GALAXY_MODEL); // these are set in initialize @@ -1066,5 +1274,5 @@ } -void usage () { +void usage (int argc, char **argv) { fprintf (stderr, "ERROR: USAGE: relastro -images -update-simple [options]\n"); fprintf (stderr, " OR: relastro -images -update-chips [options]\n"); @@ -1116,8 +1324,15 @@ fprintf (stderr, " -v\n"); fprintf (stderr, " \n"); + + fprintf (stderr, "remaining args: "); + for (int i = 0; i < argc; i++) { + fprintf (stderr, "%s ", argv[i]); + } + fprintf (stderr, "\n"); + exit (2); } -void usage_client () { +void usage_client (int argc, char **argv) { fprintf (stderr, "ERROR: USAGE: relastro_client -load\n"); fprintf (stderr, " OR: relastro_client -update-offsets\n"); @@ -1150,4 +1365,11 @@ fprintf (stderr, " -v\n"); fprintf (stderr, " \n"); + + fprintf (stderr, "remaining args: "); + for (int i = 0; i < argc; i++) { + fprintf (stderr, "%s ", argv[i]); + } + fprintf (stderr, "\n"); + exit (2); } @@ -1189,8 +1411,49 @@ } + // this sets the last loops to match the last value... while (Nloop < NLOOP) { weights[Nloop] = weights[Nloop - 1]; Nloop ++; } + return weights; } + +int *ParseLoopOrder (char *rawlist, int minValue) { + + int *orders = NULL; + ALLOCATE (orders, int, NLOOP); + + int Nloop = 0; + + /* parse the comma-separated list of photcodes */ + char *myList = strcreate(rawlist); + char *list = myList; + char *entry = NULL; + char *ptr = NULL; + while ((Nloop < NLOOP) && ((entry = strtok_r (list, ",", &ptr)) != NULL)) { + list = NULL; // pass NULL on successive strtok_r calls + + orders[Nloop] = atoi(entry); + if (orders[Nloop] < minValue) { + fprintf (stderr, "order cannot be < %d: %s\n", minValue, rawlist); + exit (3); + } + + Nloop ++; + } + free (myList); + + if (Nloop == 0) { + fprintf (stderr, "syntax error parsing orders: %s\n", rawlist); + exit (3); + } + + // this sets the last loops to match the last value... + while (Nloop < NLOOP) { + orders[Nloop] = orders[Nloop - 1]; + Nloop ++; + } + + return orders; +} Index: trunk/Ohana/src/relastro/src/assign_images.c =================================================================== --- trunk/Ohana/src/relastro/src/assign_images.c (revision 39693) +++ trunk/Ohana/src/relastro/src/assign_images.c (revision 39926) @@ -40,5 +40,5 @@ // register the image array with ImageOps.c for later getimageByID calls - initImages (image, NULL, Nimage); + initImages (image, NULL, Nimage, FALSE); if (VERBOSE) fprintf (stderr, "finding images\n"); @@ -72,4 +72,10 @@ for (j = 0; j < Nimage; j++) { + + // allow certain cameras to stay static + if (SKIP_PS1_CHIP && isGPC1chip (image[j].photcode)) continue; + if (SKIP_PS1_STACK && isGPC1stack(image[j].photcode)) continue; + if (SKIP_HSC && isHSCchip (image[j].photcode)) continue; + if (SKIP_CFH && isCFHchip (image[j].photcode)) continue; /* select images by photcode, or equiv photcode, if specified */ Index: trunk/Ohana/src/relastro/src/bcatalog.c =================================================================== --- trunk/Ohana/src/relastro/src/bcatalog.c (revision 39693) +++ trunk/Ohana/src/relastro/src/bcatalog.c (revision 39926) @@ -72,4 +72,6 @@ int myNskip1 = 0, myNskip2 = 0, myNskip3 = 0, myNskip4 = 0, myNskip5 = 0, myNskip6 = 0; + + int NgaiaObject = 0; /* exclude stars not in range or with too few measurements */ @@ -187,4 +189,10 @@ if (isGPC1warp(catalog[0].measure[offset].photcode)) continue; + // allow certain cameras to stay static + if (SKIP_PS1_CHIP && isGPC1chip(catalog[0].measure[offset].photcode)) continue; + if (SKIP_PS1_STACK && isGPC1stack(catalog[0].measure[offset].photcode)) continue; + if (SKIP_HSC && isHSCchip(catalog[0].measure[offset].photcode)) continue; + if (SKIP_CFH && isCFHchip(catalog[0].measure[offset].photcode)) continue; + // filter objects based on user supplied criteria, including SIGMA_LIM if (!MeasFilterTest(&catalog[0].measure[offset], TRUE)) { @@ -250,4 +258,6 @@ } + if (catalog[0].measure[offset].photcode == 1030) { NgaiaObject ++; } + CopyMeasureToTiny (&subcatalog[0].measureT[Nmeasure], &catalog[0].measure[offset]); // subcatalog[0].measure[Nmeasure] = catalog[0].measure[offset]; @@ -281,5 +291,5 @@ } } - fprintf (stderr, "skips: %d %d %d %d %d %d\n", myNskip1, myNskip2, myNskip3, myNskip4, myNskip5, myNskip6); + fprintf (stderr, "skips: %d %d %d %d %d %d, Ngaia: %d\n", myNskip1, myNskip2, myNskip3, myNskip4, myNskip5, myNskip6, NgaiaObject); REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1)); REALLOCATE (subcatalog[0].measureT, MeasureTiny, MAX (Nmeasure, 1)); Index: trunk/Ohana/src/relastro/src/extra.c =================================================================== --- trunk/Ohana/src/relastro/src/extra.c (revision 39693) +++ trunk/Ohana/src/relastro/src/extra.c (revision 39926) @@ -39,2 +39,26 @@ return FALSE; } + +// for now (20160925) I need to identify HSC chips explicitly. generalize in the future +int isHSCchip (int photcode) { + + if ((photcode >= 20000) && (photcode <= 20111)) return TRUE; // g-band + if ((photcode >= 21000) && (photcode <= 21111)) return TRUE; // r-band + if ((photcode >= 22000) && (photcode <= 22111)) return TRUE; // i-band + if ((photcode >= 23000) && (photcode <= 23111)) return TRUE; // z-band + if ((photcode >= 24000) && (photcode <= 24111)) return TRUE; // y-band + + return FALSE; +} + +// for now (20160925) I need to identify CFH chips explicitly. generalize in the future +int isCFHchip (int photcode) { + + if ((photcode >= 100) && (photcode <= 152)) return TRUE; // g-band + if ((photcode >= 200) && (photcode <= 252)) return TRUE; // r-band + if ((photcode >= 300) && (photcode <= 352)) return TRUE; // i-band + if ((photcode >= 400) && (photcode <= 452)) return TRUE; // z-band + if ((photcode >= 500) && (photcode <= 552)) return TRUE; // y-band + + return FALSE; +} Index: trunk/Ohana/src/relastro/src/initialize.c =================================================================== --- trunk/Ohana/src/relastro/src/initialize.c (revision 39693) +++ trunk/Ohana/src/relastro/src/initialize.c (revision 39926) @@ -11,4 +11,11 @@ ConfigInit (&argc, argv); args (argc, argv); + + if (USE_GALAXY_MODEL) { + if (!InitGalaxyModel (GALAXY_MODEL)) { + fprintf (stderr, "failed to init galaxy model %s\n", GALAXY_MODEL); + exit (2); + } + } if (RELASTRO_OP == OP_MERGE_SOURCE) return; Index: trunk/Ohana/src/relastro/src/launch_region_hosts.c =================================================================== --- trunk/Ohana/src/relastro/src/launch_region_hosts.c (revision 39693) +++ trunk/Ohana/src/relastro/src/launch_region_hosts.c (revision 39926) @@ -81,10 +81,26 @@ strextend (&command, "-region-hosts %s", REGION_FILE); strextend (&command, "-region-hostID %d", host->hostID); + strextend (&command, "-D CATDIR %s", CATDIR); strextend (&command, "-region %f %f %f %f", host->RminCat, host->RmaxCat, host->DminCat, host->DmaxCat); strextend (&command, "-statmode %s", STATMODE); strextend (&command, "-minerror %f", MIN_ERROR); - strextend (&command, "-nloop %d", NLOOP); - strextend (&command, "-threads %d", NTHREADS); + + strextend (&command, "-D RELASTRO_SIGMA_LIM %f", SIGMA_LIM); + strextend (&command, "-D RELASTRO_SRC_MEAS_TOOFEW %d", SRC_MEAS_TOOFEW); + + strextend (&command, " -D RELASTRO_MIN_DISTANCE_MOD %f", MIN_DISTANCE_MOD); + strextend (&command, " -D RELASTRO_MAX_DISTANCE_MOD %f", MAX_DISTANCE_MOD); + strextend (&command, " -D RELASTRO_MAX_DISTANCE_MOD_ERR %f", MAX_DISTANCE_MOD_ERR); + + strextend (&command, "-D USE_GALAXY_MODEL %d", USE_GALAXY_MODEL); + strextend (&command, "-D USE_ICRF_CORRECT %d", USE_ICRF_CORRECT); + + strextend (&command, "-D RELASTRO_DPOS_MAX %f", DPOS_MAX); + strextend (&command, "-D ADDSTAR_RADIUS %f", ADDSTAR_RADIUS); + + strextend (&command, "-D USE_ICRF_LOCAL %d", USE_ICRF_LOCAL); + strextend (&command, "-D USE_ICRF_SHFIT %d", USE_ICRF_SHFIT); + strextend (&command, "-D USE_ICRF_POLE %d", USE_ICRF_POLE); switch (FIT_TARGET) { @@ -98,4 +114,7 @@ strextend (&command, "-update-mosaics"); break; + case SET_CHIPS: + strextend (&command, "-set-chips"); + break; case TARGET_NONE: abort(); @@ -105,28 +124,46 @@ if (VERBOSE2) strextend (&command, "-vv"); if (RESET) strextend (&command, "-reset"); + + if (ImagSelect) strextend (&command, "-instmag %f %f", ImagMin, ImagMax); + if (MaxDensityUse) strextend (&command, "-max-density %f", MaxDensityValue); + if (FlagOutlier) strextend (&command, "-clip %d", CLIP_THRESH); + if (ExcludeBogus) strextend (&command, "-exclude-bogus %f", ExcludeBogusRadius); + + if (USE_FIXED_PIXCOORDS) strextend (&command, "-D USE_FIXED_PIXCOORDS 1"); + if (PHOTCODE_KEEP_LIST) strextend (&command, "+photcode %s", PHOTCODE_KEEP_LIST); + if (PHOTCODE_SKIP_LIST) strextend (&command, "-photcode %s", PHOTCODE_SKIP_LIST); + if (PhotFlagSelect) strextend (&command, "+photflags"); + if (PhotFlagBad) strextend (&command, "+photflagbad %d", PhotFlagBad); + if (PhotFlagPoor) strextend (&command, "+photflagpoor %d", PhotFlagPoor); + + if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) { + strextend (&command, "-dcr-blue-color %s %s", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG); + } + if (DCR_RED_COLOR_POS && DCR_RED_COLOR_NEG) { + strextend (&command, "-dcr-red-color %s %s", DCR_RED_COLOR_POS, DCR_RED_COLOR_NEG); + } + + if (TEST_SCALE != 1.0) strextend (&command, "-testing %f", TEST_SCALE); + + if (SKIP_PS1_CHIP) strextend (&command, "-skip-ps1-chip"); + if (SKIP_PS1_STACK) strextend (&command, "-skip-ps1-stack"); + if (SKIP_HSC) strextend (&command, "-skip-hsc"); + if (SKIP_CFH) strextend (&command, "-skip-cfh"); + + strextend (&command, "-nloop %d", NLOOP); + strextend (&command, "-threads %d", NTHREADS); + + if (PHOTCODE_RESET_LIST) strextend (&command, "-reset-to-photcode %s", PHOTCODE_RESET_LIST); + if (UPDATE) strextend (&command, "-update"); if (PARALLEL) strextend (&command, "-parallel"); if (PARALLEL_MANUAL) strextend (&command, "-parallel-manual"); if (PARALLEL_SERIAL) strextend (&command, "-parallel-serial"); - if (PHOTCODE_KEEP_LIST) strextend (&command, "+photcode %s", PHOTCODE_KEEP_LIST); - if (PHOTCODE_SKIP_LIST) strextend (&command, "-photcode %s", PHOTCODE_SKIP_LIST); - if (PHOTCODE_RESET_LIST) strextend (&command, "-reset-to-photcode %s", PHOTCODE_RESET_LIST); - - if (MaxDensityUse) strextend (&command, "-max-density %f", MaxDensityValue); - if (ImagSelect) strextend (&command, "-instmag %f %f", ImagMin, ImagMax); - if (ExcludeBogus) strextend (&command, "-exclude-bogus %f", ExcludeBogusRadius); - - if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) { - strextend (&command, "-dcr-blue-color %s %s", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG); - } - if (DCR_RED_COLOR_POS && DCR_RED_COLOR_NEG) { - strextend (&command, "-dcr-red-color %s %s", DCR_RED_COLOR_POS, DCR_RED_COLOR_NEG); - } - - if (PhotFlagSelect) strextend (&command, "+photflags"); - if (PhotFlagBad) strextend (&command, "+photflagbad %d", PhotFlagBad); - if (PhotFlagPoor) strextend (&command, "+photflagpoor %d", PhotFlagPoor); - + + strextend (&command, "-chiporder %d", CHIPORDER); if (CHIPMAP) strextend (&command, "-chipmap %d", CHIPMAP); + if (ChipMapLoop) strextend (&command, "-chipmaploop %s", ChipMapLoopStr); + if (ChipOrderLoop) strextend (&command, "-chiporderloop %s", ChipOrderLoopStr); + if (RESET_IMAGES) strextend (&command, "-reset-images"); @@ -136,18 +173,6 @@ if (LoopWeight2MASS) { strextend (&command, "-loop-weights-2mass %s", LoopWeight2MASSstr); } if (LoopWeightTycho) { strextend (&command, "-loop-weights-tycho %s", LoopWeightTychostr); } - - strextend (&command, "-D RELASTRO_SRC_MEAS_TOOFEW %d", SRC_MEAS_TOOFEW); - strextend (&command, "-D RELASTRO_SIGMA_LIM %f", SIGMA_LIM); - strextend (&command, "-D RELASTRO_DPOS_MAX %f", DPOS_MAX); - strextend (&command, "-D ADDSTAR_RADIUS %f", ADDSTAR_RADIUS); - - strextend (&command, "-D USE_GALAXY_MODEL %d", USE_GALAXY_MODEL); - - strextend (&command, "-D USE_ICRF_CORRECT %d", USE_ICRF_CORRECT); - strextend (&command, "-D USE_ICRF_LOCAL %d", USE_ICRF_LOCAL); - strextend (&command, "-D USE_ICRF_SHFIT %d", USE_ICRF_SHFIT); - strextend (&command, "-D USE_ICRF_POLE %d", USE_ICRF_POLE); - - if (USE_FIXED_PIXCOORDS) strextend (&command, "-D USE_FIXED_PIXCOORDS 1"); + if (LoopWeightGAIA) { strextend (&command, "-loop-weights-gaia %s", LoopWeightGAIAstr); } + if (APPLY_PROPER_MOTION) strextend (&command, "-apply-proper-motion"); if (TimeSelect) { Index: trunk/Ohana/src/relastro/src/load_catalogs.c =================================================================== --- trunk/Ohana/src/relastro/src/load_catalogs.c (revision 39693) +++ trunk/Ohana/src/relastro/src/load_catalogs.c (revision 39926) @@ -178,37 +178,34 @@ strextend (&command, "relastro_client -load-objects %s", table->hosts[i].results); strextend (&command, " -hostID %d", table->hosts[i].hostID); + strextend (&command, " -hostdir %s", table->hosts[i].pathname); + strextend (&command, " -D CATDIR %s", CATDIR); - strextend (&command, " -hostdir %s", table->hosts[i].pathname); strextend (&command, " -region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); strextend (&command, " -statmode %s", STATMODE); strextend (&command, " -minerror %f", MIN_ERROR); + strextend (&command, " -D RELASTRO_SIGMA_LIM %f", SIGMA_LIM); + strextend (&command, " -D RELASTRO_SRC_MEAS_TOOFEW %d", SRC_MEAS_TOOFEW); + strextend (&command, " -D RELASTRO_MIN_DISTANCE_MOD %f", MIN_DISTANCE_MOD); strextend (&command, " -D RELASTRO_MAX_DISTANCE_MOD %f", MAX_DISTANCE_MOD); strextend (&command, " -D RELASTRO_MAX_DISTANCE_MOD_ERR %f", MAX_DISTANCE_MOD_ERR); + strextend (&command, "-D USE_GALAXY_MODEL %d", USE_GALAXY_MODEL); + strextend (&command, "-D USE_ICRF_CORRECT %d", USE_ICRF_CORRECT); + if (FIT_MODE == FIT_PM_ONLY) strextend (&command, "-pm"); if (FIT_MODE == FIT_PAR_ONLY) strextend (&command, "-par"); if (FIT_MODE == FIT_PM_AND_PAR) strextend (&command, "-pmpar"); - if (VERBOSE) strextend (&command, "-v"); - if (VERBOSE2) strextend (&command, "-vv"); - if (RESET) strextend (&command, "-reset"); - if (ImagSelect) strextend (&command, "-instmag %f %f", ImagMin, ImagMax); - if (MaxDensityUse) strextend (&command, "-max-density %f", MaxDensityValue); - if (FlagOutlier) strextend (&command, "-clip %d", CLIP_THRESH); - if (ExcludeBogus) strextend (&command, "-exclude-bogus %f", ExcludeBogusRadius); - - if (USE_ICRF_CORRECT) strextend (&command, "-D USE_ICRF_CORRECT %d", USE_ICRF_CORRECT); - if (USE_GALAXY_MODEL) strextend (&command, "-D USE_GALAXY_MODEL %d", USE_GALAXY_MODEL); - - if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) { - strextend (&command, "-dcr-blue-color %s %s", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG); - } - if (DCR_RED_COLOR_POS && DCR_RED_COLOR_NEG) { - strextend (&command, "-dcr-red-color %s %s", DCR_RED_COLOR_POS, DCR_RED_COLOR_NEG); - } - - if (USE_ALL_IMAGES) strextend (&command, "-use-all-images"); + if (VERBOSE) strextend (&command, "-v"); + if (VERBOSE2) strextend (&command, "-vv"); + if (RESET) strextend (&command, "-reset"); + + if (ImagSelect) strextend (&command, "-instmag %f %f", ImagMin, ImagMax); + if (MaxDensityUse) strextend (&command, "-max-density %f", MaxDensityValue); + if (FlagOutlier) strextend (&command, "-clip %d", CLIP_THRESH); + if (ExcludeBogus) strextend (&command, "-exclude-bogus %f", ExcludeBogusRadius); + if (USE_FIXED_PIXCOORDS) strextend (&command, "-D USE_FIXED_PIXCOORDS 1"); if (PHOTCODE_KEEP_LIST) strextend (&command, "+photcode %s", PHOTCODE_KEEP_LIST); @@ -217,4 +214,19 @@ if (PhotFlagBad) strextend (&command, "+photflagbad %d", PhotFlagBad); if (PhotFlagPoor) strextend (&command, "+photflagpoor %d", PhotFlagPoor); + + if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) { + strextend (&command, "-dcr-blue-color %s %s", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG); + } + if (DCR_RED_COLOR_POS && DCR_RED_COLOR_NEG) { + strextend (&command, "-dcr-red-color %s %s", DCR_RED_COLOR_POS, DCR_RED_COLOR_NEG); + } + + if (SKIP_PS1_CHIP) strextend (&command, "-skip-ps1-chip"); + if (SKIP_PS1_STACK) strextend (&command, "-skip-ps1-stack"); + if (SKIP_HSC) strextend (&command, "-skip-hsc"); + if (SKIP_CFH) strextend (&command, "-skip-cfh"); + + if (USE_ALL_IMAGES) strextend (&command, "-use-all-images"); + // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values @@ -282,4 +294,6 @@ CatalogSplitter *catalogs = BrightCatalogSplitInit (Nsecfilt); + ohana_memstats (TRUE); + for (i = 0; i < table->Nhosts; i++) { @@ -296,4 +310,6 @@ free (bcatalog->secfilt); free (bcatalog); + + ohana_memstats (TRUE); } @@ -320,4 +336,5 @@ BrightCatalogSplitFree (catalogs); + ohana_memstats (TRUE); return (catalog); Index: trunk/Ohana/src/relastro/src/load_images.c =================================================================== --- trunk/Ohana/src/relastro/src/load_images.c (revision 39693) +++ trunk/Ohana/src/relastro/src/load_images.c (revision 39926) @@ -66,5 +66,5 @@ } - initImages (subset, LineNumber, Nsubset); + initImages (subset, LineNumber, Nsubset, !USE_ALL_IMAGES); MARKTIME(" init images: %f sec\n", dtime); Index: trunk/Ohana/src/relastro/src/relastro_images.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_images.c (revision 39693) +++ trunk/Ohana/src/relastro/src/relastro_images.c (revision 39926) @@ -16,4 +16,7 @@ int finalPassMode = FIT_MODE; // start with the globally-defined fit mode FIT_MODE = FIT_AVERAGE; + + int RESET_ON_UPDATE = RESET; + RESET = TRUE; // we need to reset when we load the bright catalog subset /* lock and load the image db table */ @@ -65,4 +68,6 @@ // XXX NOTE : for 2mass reset, photcodesKeep should now limit to 2MASS measurements + USE_IRLS = FALSE; // do not use IRLS yet -- leads to excessive outlier rejections in the loops + /* major modes */ switch (FIT_TARGET) { @@ -90,4 +95,11 @@ break; + case SET_CHIPS: + // we just want to fit the selected chips to the mean positions + UpdateChips (catalog, Ncatalog, 0); // measure.X,Y -> R,D, fit image.coords + MARKTIME("update chips: %f sec\n", dtime); + + break; + case TARGET_MOSAICS: for (i = 0; i < NLOOP; i++) { @@ -102,12 +114,4 @@ } - if (!UPDATE) { - freeStarMaps(); - gfits_db_free (&db); - ohana_memcheck (VERBOSE); - ohana_memdump (VERBOSE); - exit (0); - } - // free the image / measurement pointers freeImageBins (Ncatalog); @@ -117,4 +121,12 @@ free (catalog); freeMosaics (); + + if (!UPDATE) { + freeStarMaps(); + dvo_image_unlock (&db); + freeImages (db.ftable.buffer); + gfits_db_free (&db); + return TRUE; + } // If we did NOT use all images, then we applied the measured corrections to a subset of @@ -154,4 +166,6 @@ // iterate over catalogs to make detection coordinates consistant if (APPLY_OFFSETS) { + USE_IRLS = ALLOW_IRLS; // now that we have fitted the images, choose the user's option + RESET = RESET_ON_UPDATE; UpdateObjectOffsets (skylist, 0, NULL); } Index: trunk/Ohana/src/relastro/src/relastro_objects.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_objects.c (revision 39693) +++ trunk/Ohana/src/relastro/src/relastro_objects.c (revision 39926) @@ -159,8 +159,13 @@ strextend (&command, "relastro_client -update-objects"); strextend (&command, "-hostID %d", table->hosts[i].hostID); + strextend (&command, "-hostdir %s", table->hosts[i].pathname); + strextend (&command, "-D CATDIR %s", CATDIR); - strextend (&command, "-hostdir %s", table->hosts[i].pathname); strextend (&command, "-region %f %f %f %f", UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax); strextend (&command, "-statmode %s", STATMODE); + strextend (&command, "-minerror %f", MIN_ERROR); + + strextend (&command, "-D RELASTRO_SIGMA_LIM %f", SIGMA_LIM); + strextend (&command, "-D RELASTRO_SRC_MEAS_TOOFEW %d", SRC_MEAS_TOOFEW); if (FIT_MODE == FIT_PM_ONLY) { strextend (&command, "-pm"); } @@ -171,10 +176,10 @@ if (VERBOSE2) { strextend (&command, "-vv"); } if (RESET) { strextend (&command, "-reset"); } - if (UPDATE) { strextend (&command, "-update"); } + if (ImagSelect) { strextend (&command, "-instmag %f %f", ImagMin, ImagMax); } if (MaxDensityUse) { strextend (&command, "-max-density %f", MaxDensityValue); } if (FlagOutlier) { strextend (&command, "-clip %d", CLIP_THRESH); } - - if (USE_ALL_IMAGES) { strextend (&command, "-use-all-images"); } + if (ExcludeBogus) strextend (&command, "-exclude-bogus %f", ExcludeBogusRadius); + if (USE_FIXED_PIXCOORDS) { strextend (&command, "-D USE_FIXED_PIXCOORDS 1"); } if (PHOTCODE_KEEP_LIST) { strextend (&command, "+photcode %s", PHOTCODE_KEEP_LIST); } @@ -185,4 +190,14 @@ // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values + if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) { + strextend (&command, "-dcr-blue-color %s %s", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG); + } + if (DCR_RED_COLOR_POS && DCR_RED_COLOR_NEG) { + strextend (&command, "-dcr-red-color %s %s", DCR_RED_COLOR_POS, DCR_RED_COLOR_NEG); + } + + if (UPDATE) { strextend (&command, "-update"); } + if (USE_ALL_IMAGES) { strextend (&command, "-use-all-images"); } + if (MinBadQF > 0.0) strextend (&command, "-min-bad-psfqf %f", MinBadQF); if (MaxMeanOffset != 10.0) strextend (&command, "-max-mean-offset %f", MaxMeanOffset); Index: trunk/Ohana/src/relastro/src/relastro_parallel_images.c =================================================================== --- trunk/Ohana/src/relastro/src/relastro_parallel_images.c (revision 39693) +++ trunk/Ohana/src/relastro/src/relastro_parallel_images.c (revision 39926) @@ -61,5 +61,5 @@ initMosaics (image, Nimage); - initImages (image, NULL, Nimage); + initImages (image, NULL, Nimage, FALSE); SkyTable *sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE); @@ -103,4 +103,8 @@ client_logger_message ("starting the loops: %s\n", myHostName); + + RESET = TRUE; // we need to reset when we load the bright catalog subset + FIT_MODE = FIT_AVERAGE; // we need to only fit the average + USE_IRLS = FALSE; // do not use IRLS yet -- leads to excessive outlier rejections in the loops /* major modes */ Index: trunk/Ohana/src/relastro/src/select_images.c =================================================================== --- trunk/Ohana/src/relastro/src/select_images.c (revision 39693) +++ trunk/Ohana/src/relastro/src/select_images.c (revision 39926) @@ -121,4 +121,10 @@ } + // allow certain cameras to stay static + if (SKIP_PS1_CHIP && isGPC1chip (timage[i].photcode)) continue; + if (SKIP_PS1_STACK && isGPC1stack(timage[i].photcode)) continue; + if (SKIP_HSC && isHSCchip (timage[i].photcode)) continue; + if (SKIP_CFH && isCFHchip (timage[i].photcode)) continue; + /* select images by photcode, or equiv photcode, if specified */ if (NphotcodesKeep > 0) { @@ -187,4 +193,6 @@ if (DmaxImage < DminSkyRegion) continue; + if (KEEP_ALL_IMAGES_RA) goto found_it; + // the sky region RA is defined to be 0 - 360.0 if (RminImage > RmaxSkyRegion) continue; Index: trunk/Ohana/src/relastro/src/syncfile.c =================================================================== --- trunk/Ohana/src/relastro/src/syncfile.c (revision 39693) +++ trunk/Ohana/src/relastro/src/syncfile.c (revision 39926) @@ -39,4 +39,5 @@ int loop; sscanf (message, "%*s %d", &loop); + if (CATCH_UP && (nloop < loop)) return TRUE; // if I am behind this machine, I need to catch up! if (loop != nloop) { usleep (2000000);