Index: /trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- /trunk/Ohana/src/relastro/include/relastro.h (revision 39370) +++ /trunk/Ohana/src/relastro/include/relastro.h (revision 39371) @@ -310,6 +310,6 @@ int SHOW_PARAMS; char STATMODE[32]; -int POS_TOOFEW; -int PM_TOOFEW; +// int POS_TOOFEW; +// int PM_TOOFEW; double PM_DT_MIN; double PAR_FACTOR_MIN; @@ -759,4 +759,5 @@ int FitPMandPar_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); int FitPosPMfixed_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); +int FitPosPMfixed_Single (FitAstromResult *fit, FitAstromPoint *points, int Npoints); double MedianAbsDeviation(FitAstromPoint *points, int Npoints); Index: /trunk/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- /trunk/Ohana/src/relastro/src/ConfigInit.c (revision 39370) +++ /trunk/Ohana/src/relastro/src/ConfigInit.c (revision 39371) @@ -28,6 +28,6 @@ if (!ScanConfig (config, "PM_DT_MIN", "%lf", 0, &PM_DT_MIN)) PM_DT_MIN = 0.25; - if (!ScanConfig (config, "PM_TOOFEW", "%d", 0, &PM_TOOFEW)) PM_TOOFEW = 4; - if (!ScanConfig (config, "POS_TOOFEW", "%d", 0, &POS_TOOFEW)) POS_TOOFEW = 1; + // if (!ScanConfig (config, "PM_TOOFEW", "%d", 0, &PM_TOOFEW)) PM_TOOFEW = 4; + // if (!ScanConfig (config, "POS_TOOFEW", "%d", 0, &POS_TOOFEW)) POS_TOOFEW = 1; if (!ScanConfig (config, "PAR_FACTOR_MIN", "%lf", 0, &PAR_FACTOR_MIN)) PAR_FACTOR_MIN = 0.2; if (!ScanConfig (config, "RELASTRO_MAP_NX", "%d", 0, &NX_MAP)) NX_MAP = 5; Index: /trunk/Ohana/src/relastro/src/FitPM.c =================================================================== --- /trunk/Ohana/src/relastro/src/FitPM.c (revision 39370) +++ /trunk/Ohana/src/relastro/src/FitPM.c (revision 39371) @@ -20,8 +20,6 @@ } - int status = FitPM_MinChisq (fit, data, points, Npoints); - if (!status) return FALSE; - - FitPM_SetChisq (fit, data, points, Npoints); + if (!FitPM_MinChisq (fit, data, points, Npoints)) return FALSE; + if (!FitPM_SetChisq (fit, data, points, Npoints)) return FALSE; return TRUE; } Index: /trunk/Ohana/src/relastro/src/FitPMandPar.c =================================================================== --- /trunk/Ohana/src/relastro/src/FitPMandPar.c (revision 39370) +++ /trunk/Ohana/src/relastro/src/FitPMandPar.c (revision 39371) @@ -20,8 +20,7 @@ } - int status = FitPMandPar_MinChisq (fit, data, points, Npoints); - if (!status) return FALSE; - - FitPMandPar_SetChisq (fit, data, points, Npoints); + if (!FitPMandPar_MinChisq (fit, data, points, Npoints)) return FALSE; + + if (!FitPMandPar_SetChisq (fit, data, points, Npoints)) return FALSE; return TRUE; } @@ -331,6 +330,4 @@ // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms) fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms); - - return TRUE; } Index: /trunk/Ohana/src/relastro/src/FitPosPMfixed.c =================================================================== --- /trunk/Ohana/src/relastro/src/FitPosPMfixed.c (revision 39370) +++ /trunk/Ohana/src/relastro/src/FitPosPMfixed.c (revision 39371) @@ -9,4 +9,32 @@ # define WEIGHT_THRESHOLD 0.3 +int FitPosPMfixed_Single (FitAstromResult *fit, FitAstromPoint *points, int Npoints) { + + int i; + + // find a single unmasked point + int found = -1; + for (i = 0; (found == -1) && (i < Npoints); i++) { + if (points[i].mask) continue; // respect the mask if set + fit->Ro = points[i].X - fit->uR*points[i].T; + fit->Do = points[i].Y - fit->uD*points[i].T; + fit->dRo = points[i].dX; + fit->dDo = points[i].dY; + + fit->p = 0.0; + fit->duR = 0.0; + fit->duD = 0.0; + fit->dp = 0.0; + + found = i; + } + if (found == -1) return FALSE; + for (i = 0; (i < Npoints); i++) { + if (i == found) continue; + points[i].mask = TRUE; + } + return TRUE; +} + // fit->uR,uD are used in the fit int FitPosPMfixed_Basic (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) { @@ -20,8 +48,13 @@ } - int status = FitPosPMfixed_MinChisq (fit, data, points, Npoints); - if (!status) return FALSE; - - FitPosPMfixed_SetChisq (fit, data, points, Npoints); + if (!FitPosPMfixed_MinChisq (fit, data, points, Npoints)) { + if (!FitPosPMfixed_Single (fit, points, Npoints)) return FALSE; + return TRUE; + } + + if (!FitPosPMfixed_SetChisq (fit, data, points, Npoints)) { + if (!FitPosPMfixed_Single (fit, points, Npoints)) return FALSE; + return TRUE; + } return TRUE; } Index: /trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- /trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 39370) +++ /trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 39371) @@ -1,4 +1,12 @@ # include "relastro.h" -# define PAR_TOOFEW 5 + +// for good time sampling, these values are probably conservative +// but with tti pairs, these may not be sufficient... +# define PAR_MIN_NPTS 5 +# define PAR_MIN_NPTS_BOOT 4 +# define PM_MIN_NPTS 4 +# define PM_MIN_NPTS_BOOT 3 +# define POS_MIN_NPTS 3 +# define POS_MIN_NPTS_BOOT 2 typedef enum { @@ -160,5 +168,5 @@ goto skipParallax; } - if (fitStats->Npoints <= PAR_TOOFEW) { + if (fitStats->Npoints < PAR_MIN_NPTS) { // not enough data, skip parallax mode = FIT_PM_ONLY; @@ -177,6 +185,6 @@ fitStats->Nfit = 0; int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); - if (Nnomask < 5) { - // if we do not have enough points to fit parallax, we cannot do the bootstrap analysis. + 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; @@ -229,5 +237,5 @@ goto justPosition; } - if (fitStats->Npoints <= PM_TOOFEW) { + if (fitStats->Npoints < PM_MIN_NPTS) { mode = FIT_AVERAGE; goto justPosition; @@ -249,5 +257,5 @@ fitStats->Nfit = 0; int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); - if (Nnomask < 5) { + if (Nnomask < PM_MIN_NPTS_BOOT) { mode = FIT_AVERAGE; goto justPosition; @@ -297,9 +305,18 @@ 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_STAR_USE_PAR) || (average[0].flags & ID_STAR_USE_PM)) { + if (average[0].flags & (ID_STAR_USE_PAR | ID_STAR_USE_PM)) { if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { 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)) { + goto doneWithFit; + } + // XXX Maybe add a flag here? + average[0].flags |= ID_STAR_USE_AVE; + goto useBasic; + } if (!FitPosPMfixed_IRLS (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { goto doneWithFit; @@ -309,6 +326,11 @@ fitStats->Nfit = 0; int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); - if (Nnomask < 5) { - goto doneWithFit; + if (Nnomask < POS_MIN_NPTS_BOOT) { + if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { + goto doneWithFit; + } + // XXX Maybe add a flag here? + average[0].flags |= ID_STAR_USE_AVE; + goto useBasic; } FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average); @@ -324,4 +346,5 @@ } + useBasic: // project Ro, Do back to RA,DEC XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords); @@ -333,4 +356,19 @@ doneWithFit: + // if no valid fit has been found, try to use a single unmasked point: + if (!(average[0].flags & (ID_STAR_USE_PAR | ID_STAR_USE_PM | ID_STAR_USE_AVE))) { + if (!FitPosPMfixed_Single (&fitPos, fitStats->points, fitStats->Npoints)) { + fitStats->Nskip ++; + return FALSE; + } + // project Ro, Do back to RA,DEC + XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords); + if (fabs(fitPos.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); + + average[0].flags |= ID_STAR_USE_AVE; + average[0].flags |= ID_STAR_FIT_AVE; + fitStats->Nave ++; + } + // update the bit flags of which points were used for (k = 0; k < fitStats->Npoints; k++) { Index: /trunk/Ohana/src/relastro/src/fitobj2.c =================================================================== --- /trunk/Ohana/src/relastro/src/fitobj2.c (revision 39370) +++ /trunk/Ohana/src/relastro/src/fitobj2.c (revision 39371) @@ -1,3 +1,10 @@ # include "relastro.h" + +static int Nrand = 0; +double drand48_cnt () { + double value = drand48(); + Nrand ++; + return value; +} Catalog *mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad); @@ -109,5 +116,5 @@ long A = now.tv_sec + now.tv_usec * 1000000; srand48(A); - // srand48(1); + srand48(1); } @@ -121,15 +128,17 @@ for (i = 0; i < Nstars; i++) { - double Ro = 360.0*(drand48() - 0.5); - double Phi = 2.0*(drand48() - 0.5); + double Ro = 360.0*(drand48_cnt() - 0.5); + double Phi = 2.0*(drand48_cnt() - 0.5); double Do = DEG_RAD * asin(Phi); - double uR = 0.5*(drand48() - 0.5); - double uD = 0.5*(drand48() - 0.5); - - double plx = (FIT_MODE == FIT_PM_AND_PAR) ? 0.5*(drand48() + 0.0) : 0.0; - - int Npoints = (N_POINTS - 2)*drand48() + 2; // minimum of 2 points - int Noutliers = Npoints * F_OUTLIERS; + double uR = 0.5*(drand48_cnt() - 0.5); + double uD = 0.5*(drand48_cnt() - 0.5); + + double plx = 0.5*(drand48_cnt() + 0.0); + if (FIT_MODE == FIT_PM_AND_PAR) plx = 0.0; + + int Npoints = N_POINTS*drand48_cnt(); // minimum of 0 points + int Noutliers_max = Npoints * F_OUTLIERS; + int Noutliers = Noutliers_max*drand48_cnt(); // minimum of 0 points // generate a single fake star with N_POINTS real points and N_OUTLIERS bad points @@ -140,9 +149,12 @@ double Tmean = ohana_sec_to_mjd(catalog->average->Tmean); double Tyear = (Tmean - MJD_J2000) / 365.25; - fprintf (stdout, "%3d %3d %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f | %8.6f %8.6f %8.6f %8.6f %8.6f %d | %f %f %f\n", + fprintf (stdout, "%3d %3d %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f | %8.6f %8.6f %8.6f %8.6f %8.6f %d | %6.2f %6.2f %6.2f 0x%08x\n", Npoints, Noutliers, Ro, Do, uR, uD, plx, catalog->average->R, catalog->average->D, catalog->average->uR, catalog->average->uD, catalog->average->P, Tyear, catalog->average->dR, catalog->average->dD, catalog->average->duR, catalog->average->duD, catalog->average->dP, catalog->average->Npos, - catalog->average->ChiSqAve, catalog->average->ChiSqPM, catalog->average->ChiSqPar); + catalog->average->ChiSqAve, catalog->average->ChiSqPM, catalog->average->ChiSqPar, catalog->average->flags); + if (i > 0) { + // fprintf (stderr, ".\n"); + } continue; @@ -187,5 +199,5 @@ for (i = 0; i < Npoints; i++) { // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined) - double Tyears = MJD_MIN_YRS + drand48()*(MJD_MAX_YRS - MJD_MIN_YRS); + double Tyears = MJD_MIN_YRS + drand48_cnt()*(MJD_MAX_YRS - MJD_MIN_YRS); double Tmjd = 365.25*Tyears + MJD_J2000; @@ -214,5 +226,5 @@ for (i = Npoints; i < Ntotal; i++) { // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined) - double Tyears = MJD_MIN_YRS + drand48()*(MJD_MAX_YRS - MJD_MIN_YRS); + double Tyears = MJD_MIN_YRS + drand48_cnt()*(MJD_MAX_YRS - MJD_MIN_YRS); double Tmjd = 365.25*Tyears + MJD_J2000; @@ -222,6 +234,6 @@ ParFactor (&pR, &pD, Ro, Do, Tyears); - double dR = OUT_ERROR*(drand48() - 0.5); - double dD = OUT_ERROR*(drand48() - 0.5); + double dR = OUT_ERROR*(drand48_cnt() - 0.5); + double dD = OUT_ERROR*(drand48_cnt() - 0.5); catalog->measureT[i].R = Ro + (dR + plx*pR + uR*(Tyears - MJD_REF_YRS))/3600/dcos(Do); Index: /trunk/Ohana/src/relastro/src/fitpm.c =================================================================== --- /trunk/Ohana/src/relastro/src/fitpm.c (revision 39370) +++ /trunk/Ohana/src/relastro/src/fitpm.c (revision 39371) @@ -1,3 +1,10 @@ # include "relastro.h" + +static int Nrand = 0; +double drand48_cnt () { + double value = drand48(); + Nrand ++; + return value; +} int mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad); @@ -41,6 +48,7 @@ int N_STARS = 3000; int N_POINTS = 100; - int N_OUTLIERS = 10; + // int N_OUTLIERS = 10; int N_BOOTSTRAP = 100; + float F_OUTLIERS = 0.1; int N; @@ -55,9 +63,9 @@ remove_argument (N, &argc, argv); } - if ((N = get_argument (argc, argv, "-Noutliers"))) { - remove_argument (N, &argc, argv); - N_OUTLIERS = atoi (argv[N]); - remove_argument (N, &argc, argv); - } + // if ((N = get_argument (argc, argv, "-Noutliers"))) { + // remove_argument (N, &argc, argv); + // N_OUTLIERS = atoi (argv[N]); + // remove_argument (N, &argc, argv); + // } if ((N = get_argument (argc, argv, "-Nbootstrap"))) { remove_argument (N, &argc, argv); @@ -132,5 +140,6 @@ ohana_gaussdev_init (); - FitStats *fitStats = FitStatsInit (N_POINTS + N_OUTLIERS, N_BOOTSTRAP); + FitStats *fitStats = FitStatsInit (N_POINTS * (1 + 2*F_OUTLIERS), N_BOOTSTRAP); + // FitStats *fitStats = FitStatsInit (N_POINTS + N_OUTLIERS, N_BOOTSTRAP); fitStats->fitdataPM->getChisq = TRUE; fitStats->fitdataPM->getError = TRUE; @@ -146,12 +155,13 @@ for (i = 0; i < Nstars; i++) { - double Ro = 360.0*(drand48() - 0.5); - double Phi = 2.0*(drand48() - 0.5); + double Ro = 360.0*(drand48_cnt() - 0.5); + double Phi = 2.0*(drand48_cnt() - 0.5); double Do = DEG_RAD * asin(Phi); - double uR = 0.5*(drand48() - 0.5); - double uD = 0.5*(drand48() - 0.5); - - double plx; + double uR = 0.5*(drand48_cnt() - 0.5); + double uD = 0.5*(drand48_cnt() - 0.5); + + double plx = 0.5*(drand48_cnt() + 0.0); + switch (mode) { case FIT_POS_IRLS: @@ -166,5 +176,4 @@ case FIT_PLX_BOOT: case FIT_PLX_NOCLIP: - plx = 0.5*(drand48() + 0.0); break; default: @@ -172,6 +181,10 @@ } + int Npoints = N_POINTS*drand48_cnt(); // minimum of 0 points + int Noutliers_max = Npoints * F_OUTLIERS; + int Noutliers = Noutliers_max*drand48_cnt(); // minimum of 0 points + // generate a single fake star with N_POINTS real points and N_OUTLIERS bad points - mkstar (fitStats, Ro, Do, uR, uD, plx, N_POINTS, N_OUTLIERS); + mkstar (fitStats, Ro, Do, uR, uD, plx, Npoints, Noutliers); double Tmean, Trange, parRange; @@ -266,4 +279,7 @@ Ro, Do, uR, uD, plx, fitResult.Ro, fitResult.Do, fitResult.uR, fitResult.uD, fitResult.p, Tmean, fitResult.dRo, fitResult.dDo, fitResult.duR, fitResult.duD, fitResult.dp, fitResult.Nfit, fitResult.chisq); + if (i > 0) { + // fprintf (stderr, ".\n"); + } continue; @@ -291,5 +307,5 @@ // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined) - points[i].T = MJD_MIN_YRS + drand48()*(MJD_MAX_YRS - MJD_MIN_YRS); + points[i].T = MJD_MIN_YRS + drand48_cnt()*(MJD_MAX_YRS - MJD_MIN_YRS); double pR, pD; @@ -310,11 +326,11 @@ // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined) - points[i].T = MJD_MIN_YRS + drand48()*(MJD_MAX_YRS - MJD_MIN_YRS); + points[i].T = MJD_MIN_YRS + drand48_cnt()*(MJD_MAX_YRS - MJD_MIN_YRS); double pR, pD; ParFactor (&pR, &pD, Ro, Do, points[i].T); - double dR = OUT_ERROR*(drand48() - 0.5) + RA_BIAS; - double dD = OUT_ERROR*(drand48() - 0.5) + DEC_BIAS; + double dR = OUT_ERROR*(drand48_cnt() - 0.5) + RA_BIAS; + double dD = OUT_ERROR*(drand48_cnt() - 0.5) + DEC_BIAS; points[i].R = Ro + (dR + plx*pR + uR*(points[i].T - MJD_REF_YRS))/3600/dcos(Do); Index: /trunk/Ohana/src/relastro/test/mana.sh =================================================================== --- /trunk/Ohana/src/relastro/test/mana.sh (revision 39370) +++ /trunk/Ohana/src/relastro/test/mana.sh (revision 39371) @@ -131,7 +131,7 @@ data $1 - # read Npts 1 Nout 2 Ro 3 Do 4 uRo 5 uDo 6 Po 7 Rx 9 Dx 10 uRx 11 uDx 12 Px 13 Tx 15 dRo 17 dDo 18 duRo 19 duDo 20 dPo 21 Nfit 22 chisq 26 + read Npts 1 Nout 2 Ro 3 Do 4 uRo 5 uDo 6 Po 7 Rx 9 Dx 10 uRx 11 uDx 12 Px 13 Tx 15 dRo 17 dDo 18 duRo 19 duDo 20 dPo 21 Nfit 22 chisq 26 # read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 24 - read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 22 + # read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 22 foreach f R D uD uR P @@ -167,4 +167,26 @@ fprintf "%5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %6.2f %5.2f" $value_dRf $value_dDf $value_duR $value_duD $value_dP $value_dRo $value_dDo $value_duRo $value_duDo $value_dPo $sigma_dRo $sigma_dDo $sigma_duRo $sigma_duDo $sigma_dPo $mean_Nfit $mean_chisq +end + +macro checkstats + if ($0 != 2) + echo "USAGE: checkstats (file)" + break + end + + teststats $1 + + delete dRbin oRbin + for i 0 100 + subset tmp = dRf if (Npts == $i) + vstat tmp -q + concat $SIGMA dRbin + concat $MEAN oRbin + end + + create n 0 dRbin[] + set dRmod = 0.010 / sqrt (n) + lim n -0.001 0.015; clear; box; plot n dRbin + plot -c red -pt 7 n dRmod end