Changeset 39371
- Timestamp:
- Feb 23, 2016, 10:55:44 AM (10 years ago)
- Location:
- trunk/Ohana/src/relastro
- Files:
-
- 9 edited
-
include/relastro.h (modified) (2 diffs)
-
src/ConfigInit.c (modified) (1 diff)
-
src/FitPM.c (modified) (1 diff)
-
src/FitPMandPar.c (modified) (2 diffs)
-
src/FitPosPMfixed.c (modified) (2 diffs)
-
src/UpdateObjects.c (modified) (9 diffs)
-
src/fitobj2.c (modified) (7 diffs)
-
src/fitpm.c (modified) (10 diffs)
-
test/mana.sh (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/Ohana/src/relastro/include/relastro.h
r39365 r39371 310 310 int SHOW_PARAMS; 311 311 char STATMODE[32]; 312 int POS_TOOFEW;313 int PM_TOOFEW;312 // int POS_TOOFEW; 313 // int PM_TOOFEW; 314 314 double PM_DT_MIN; 315 315 double PAR_FACTOR_MIN; … … 759 759 int FitPMandPar_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); 760 760 int FitPosPMfixed_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints); 761 int FitPosPMfixed_Single (FitAstromResult *fit, FitAstromPoint *points, int Npoints); 761 762 762 763 double MedianAbsDeviation(FitAstromPoint *points, int Npoints); -
trunk/Ohana/src/relastro/src/ConfigInit.c
r39225 r39371 28 28 29 29 if (!ScanConfig (config, "PM_DT_MIN", "%lf", 0, &PM_DT_MIN)) PM_DT_MIN = 0.25; 30 if (!ScanConfig (config, "PM_TOOFEW", "%d", 0, &PM_TOOFEW)) PM_TOOFEW = 4;31 if (!ScanConfig (config, "POS_TOOFEW", "%d", 0, &POS_TOOFEW)) POS_TOOFEW = 1;30 // if (!ScanConfig (config, "PM_TOOFEW", "%d", 0, &PM_TOOFEW)) PM_TOOFEW = 4; 31 // if (!ScanConfig (config, "POS_TOOFEW", "%d", 0, &POS_TOOFEW)) POS_TOOFEW = 1; 32 32 if (!ScanConfig (config, "PAR_FACTOR_MIN", "%lf", 0, &PAR_FACTOR_MIN)) PAR_FACTOR_MIN = 0.2; 33 33 if (!ScanConfig (config, "RELASTRO_MAP_NX", "%d", 0, &NX_MAP)) NX_MAP = 5; -
trunk/Ohana/src/relastro/src/FitPM.c
r39370 r39371 20 20 } 21 21 22 int status = FitPM_MinChisq (fit, data, points, Npoints); 23 if (!status) return FALSE; 24 25 FitPM_SetChisq (fit, data, points, Npoints); 22 if (!FitPM_MinChisq (fit, data, points, Npoints)) return FALSE; 23 if (!FitPM_SetChisq (fit, data, points, Npoints)) return FALSE; 26 24 return TRUE; 27 25 } -
trunk/Ohana/src/relastro/src/FitPMandPar.c
r39370 r39371 20 20 } 21 21 22 int status = FitPMandPar_MinChisq (fit, data, points, Npoints); 23 if (!status) return FALSE; 24 25 FitPMandPar_SetChisq (fit, data, points, Npoints); 22 if (!FitPMandPar_MinChisq (fit, data, points, Npoints)) return FALSE; 23 24 if (!FitPMandPar_SetChisq (fit, data, points, Npoints)) return FALSE; 26 25 return TRUE; 27 26 } … … 331 330 // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms) 332 331 fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms); 333 334 335 332 return TRUE; 336 333 } -
trunk/Ohana/src/relastro/src/FitPosPMfixed.c
r39370 r39371 9 9 # define WEIGHT_THRESHOLD 0.3 10 10 11 int FitPosPMfixed_Single (FitAstromResult *fit, FitAstromPoint *points, int Npoints) { 12 13 int i; 14 15 // find a single unmasked point 16 int found = -1; 17 for (i = 0; (found == -1) && (i < Npoints); i++) { 18 if (points[i].mask) continue; // respect the mask if set 19 fit->Ro = points[i].X - fit->uR*points[i].T; 20 fit->Do = points[i].Y - fit->uD*points[i].T; 21 fit->dRo = points[i].dX; 22 fit->dDo = points[i].dY; 23 24 fit->p = 0.0; 25 fit->duR = 0.0; 26 fit->duD = 0.0; 27 fit->dp = 0.0; 28 29 found = i; 30 } 31 if (found == -1) return FALSE; 32 for (i = 0; (i < Npoints); i++) { 33 if (i == found) continue; 34 points[i].mask = TRUE; 35 } 36 return TRUE; 37 } 38 11 39 // fit->uR,uD are used in the fit 12 40 int FitPosPMfixed_Basic (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) { … … 20 48 } 21 49 22 int status = FitPosPMfixed_MinChisq (fit, data, points, Npoints); 23 if (!status) return FALSE; 24 25 FitPosPMfixed_SetChisq (fit, data, points, Npoints); 50 if (!FitPosPMfixed_MinChisq (fit, data, points, Npoints)) { 51 if (!FitPosPMfixed_Single (fit, points, Npoints)) return FALSE; 52 return TRUE; 53 } 54 55 if (!FitPosPMfixed_SetChisq (fit, data, points, Npoints)) { 56 if (!FitPosPMfixed_Single (fit, points, Npoints)) return FALSE; 57 return TRUE; 58 } 26 59 return TRUE; 27 60 } -
trunk/Ohana/src/relastro/src/UpdateObjects.c
r39370 r39371 1 1 # include "relastro.h" 2 # define PAR_TOOFEW 5 2 3 // for good time sampling, these values are probably conservative 4 // but with tti pairs, these may not be sufficient... 5 # define PAR_MIN_NPTS 5 6 # define PAR_MIN_NPTS_BOOT 4 7 # define PM_MIN_NPTS 4 8 # define PM_MIN_NPTS_BOOT 3 9 # define POS_MIN_NPTS 3 10 # define POS_MIN_NPTS_BOOT 2 3 11 4 12 typedef enum { … … 160 168 goto skipParallax; 161 169 } 162 if (fitStats->Npoints < = PAR_TOOFEW) {170 if (fitStats->Npoints < PAR_MIN_NPTS) { 163 171 // not enough data, skip parallax 164 172 mode = FIT_PM_ONLY; … … 177 185 fitStats->Nfit = 0; 178 186 int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); 179 if (Nnomask < 5) {180 // if we do not have enough points to fit parallax, we cannot do the bootstrap analysis.187 if (Nnomask < PAR_MIN_NPTS_BOOT) { 188 // if we do not have enough points to assess parallax error, we cannot do the bootstrap analysis. 181 189 mode = FIT_PM_ONLY; 182 190 goto skipParallax; … … 229 237 goto justPosition; 230 238 } 231 if (fitStats->Npoints < = PM_TOOFEW) {239 if (fitStats->Npoints < PM_MIN_NPTS) { 232 240 mode = FIT_AVERAGE; 233 241 goto justPosition; … … 249 257 fitStats->Nfit = 0; 250 258 int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); 251 if (Nnomask < 5) {259 if (Nnomask < PM_MIN_NPTS_BOOT) { 252 260 mode = FIT_AVERAGE; 253 261 goto justPosition; … … 297 305 FitAstromResultSetPM (&fitPos, 1, average); 298 306 // fprintf (stderr, "fit 1: %f %f : %f %f\n", fitPos.Ro, fitPos.Do, fitPos.uR, fitPos.uD); 299 if ( (average[0].flags & ID_STAR_USE_PAR) || (average[0].flags &ID_STAR_USE_PM)) {307 if (average[0].flags & (ID_STAR_USE_PAR | ID_STAR_USE_PM)) { 300 308 if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { 301 309 goto doneWithFit; 302 310 } 303 311 } else { 312 if (fitStats->Npoints < POS_MIN_NPTS) { 313 // I will not try to outlier-reject, just calculate the weighted average 314 if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { 315 goto doneWithFit; 316 } 317 // XXX Maybe add a flag here? 318 average[0].flags |= ID_STAR_USE_AVE; 319 goto useBasic; 320 } 304 321 if (!FitPosPMfixed_IRLS (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { 305 322 goto doneWithFit; … … 309 326 fitStats->Nfit = 0; 310 327 int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints); 311 if (Nnomask < 5) { 312 goto doneWithFit; 328 if (Nnomask < POS_MIN_NPTS_BOOT) { 329 if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) { 330 goto doneWithFit; 331 } 332 // XXX Maybe add a flag here? 333 average[0].flags |= ID_STAR_USE_AVE; 334 goto useBasic; 313 335 } 314 336 FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average); … … 324 346 } 325 347 348 useBasic: 326 349 // project Ro, Do back to RA,DEC 327 350 XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords); … … 333 356 334 357 doneWithFit: 358 // if no valid fit has been found, try to use a single unmasked point: 359 if (!(average[0].flags & (ID_STAR_USE_PAR | ID_STAR_USE_PM | ID_STAR_USE_AVE))) { 360 if (!FitPosPMfixed_Single (&fitPos, fitStats->points, fitStats->Npoints)) { 361 fitStats->Nskip ++; 362 return FALSE; 363 } 364 // project Ro, Do back to RA,DEC 365 XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords); 366 if (fabs(fitPos.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n"); 367 368 average[0].flags |= ID_STAR_USE_AVE; 369 average[0].flags |= ID_STAR_FIT_AVE; 370 fitStats->Nave ++; 371 } 372 335 373 // update the bit flags of which points were used 336 374 for (k = 0; k < fitStats->Npoints; k++) { -
trunk/Ohana/src/relastro/src/fitobj2.c
r39370 r39371 1 1 # include "relastro.h" 2 3 static int Nrand = 0; 4 double drand48_cnt () { 5 double value = drand48(); 6 Nrand ++; 7 return value; 8 } 2 9 3 10 Catalog *mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad); … … 109 116 long A = now.tv_sec + now.tv_usec * 1000000; 110 117 srand48(A); 111 //srand48(1);118 srand48(1); 112 119 } 113 120 … … 121 128 for (i = 0; i < Nstars; i++) { 122 129 123 double Ro = 360.0*(drand48 () - 0.5);124 double Phi = 2.0*(drand48 () - 0.5);130 double Ro = 360.0*(drand48_cnt() - 0.5); 131 double Phi = 2.0*(drand48_cnt() - 0.5); 125 132 double Do = DEG_RAD * asin(Phi); 126 133 127 double uR = 0.5*(drand48() - 0.5); 128 double uD = 0.5*(drand48() - 0.5); 129 130 double plx = (FIT_MODE == FIT_PM_AND_PAR) ? 0.5*(drand48() + 0.0) : 0.0; 131 132 int Npoints = (N_POINTS - 2)*drand48() + 2; // minimum of 2 points 133 int Noutliers = Npoints * F_OUTLIERS; 134 double uR = 0.5*(drand48_cnt() - 0.5); 135 double uD = 0.5*(drand48_cnt() - 0.5); 136 137 double plx = 0.5*(drand48_cnt() + 0.0); 138 if (FIT_MODE == FIT_PM_AND_PAR) plx = 0.0; 139 140 int Npoints = N_POINTS*drand48_cnt(); // minimum of 0 points 141 int Noutliers_max = Npoints * F_OUTLIERS; 142 int Noutliers = Noutliers_max*drand48_cnt(); // minimum of 0 points 134 143 135 144 // generate a single fake star with N_POINTS real points and N_OUTLIERS bad points … … 140 149 double Tmean = ohana_sec_to_mjd(catalog->average->Tmean); 141 150 double Tyear = (Tmean - MJD_J2000) / 365.25; 142 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",151 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", 143 152 Npoints, Noutliers, Ro, Do, uR, uD, plx, 144 153 catalog->average->R, catalog->average->D, catalog->average->uR, catalog->average->uD, catalog->average->P, Tyear, 145 154 catalog->average->dR, catalog->average->dD, catalog->average->duR, catalog->average->duD, catalog->average->dP, catalog->average->Npos, 146 catalog->average->ChiSqAve, catalog->average->ChiSqPM, catalog->average->ChiSqPar); 155 catalog->average->ChiSqAve, catalog->average->ChiSqPM, catalog->average->ChiSqPar, catalog->average->flags); 156 if (i > 0) { 157 // fprintf (stderr, ".\n"); 158 } 147 159 148 160 continue; … … 187 199 for (i = 0; i < Npoints; i++) { 188 200 // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined) 189 double Tyears = MJD_MIN_YRS + drand48 ()*(MJD_MAX_YRS - MJD_MIN_YRS);201 double Tyears = MJD_MIN_YRS + drand48_cnt()*(MJD_MAX_YRS - MJD_MIN_YRS); 190 202 double Tmjd = 365.25*Tyears + MJD_J2000; 191 203 … … 214 226 for (i = Npoints; i < Ntotal; i++) { 215 227 // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined) 216 double Tyears = MJD_MIN_YRS + drand48 ()*(MJD_MAX_YRS - MJD_MIN_YRS);228 double Tyears = MJD_MIN_YRS + drand48_cnt()*(MJD_MAX_YRS - MJD_MIN_YRS); 217 229 double Tmjd = 365.25*Tyears + MJD_J2000; 218 230 … … 222 234 ParFactor (&pR, &pD, Ro, Do, Tyears); 223 235 224 double dR = OUT_ERROR*(drand48 () - 0.5);225 double dD = OUT_ERROR*(drand48 () - 0.5);236 double dR = OUT_ERROR*(drand48_cnt() - 0.5); 237 double dD = OUT_ERROR*(drand48_cnt() - 0.5); 226 238 227 239 catalog->measureT[i].R = Ro + (dR + plx*pR + uR*(Tyears - MJD_REF_YRS))/3600/dcos(Do); -
trunk/Ohana/src/relastro/src/fitpm.c
r39370 r39371 1 1 # include "relastro.h" 2 3 static int Nrand = 0; 4 double drand48_cnt () { 5 double value = drand48(); 6 Nrand ++; 7 return value; 8 } 2 9 3 10 int mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad); … … 41 48 int N_STARS = 3000; 42 49 int N_POINTS = 100; 43 int N_OUTLIERS = 10;50 // int N_OUTLIERS = 10; 44 51 int N_BOOTSTRAP = 100; 52 float F_OUTLIERS = 0.1; 45 53 46 54 int N; … … 55 63 remove_argument (N, &argc, argv); 56 64 } 57 if ((N = get_argument (argc, argv, "-Noutliers"))) {58 remove_argument (N, &argc, argv);59 N_OUTLIERS = atoi (argv[N]);60 remove_argument (N, &argc, argv);61 }65 // if ((N = get_argument (argc, argv, "-Noutliers"))) { 66 // remove_argument (N, &argc, argv); 67 // N_OUTLIERS = atoi (argv[N]); 68 // remove_argument (N, &argc, argv); 69 // } 62 70 if ((N = get_argument (argc, argv, "-Nbootstrap"))) { 63 71 remove_argument (N, &argc, argv); … … 132 140 ohana_gaussdev_init (); 133 141 134 FitStats *fitStats = FitStatsInit (N_POINTS + N_OUTLIERS, N_BOOTSTRAP); 142 FitStats *fitStats = FitStatsInit (N_POINTS * (1 + 2*F_OUTLIERS), N_BOOTSTRAP); 143 // FitStats *fitStats = FitStatsInit (N_POINTS + N_OUTLIERS, N_BOOTSTRAP); 135 144 fitStats->fitdataPM->getChisq = TRUE; 136 145 fitStats->fitdataPM->getError = TRUE; … … 146 155 for (i = 0; i < Nstars; i++) { 147 156 148 double Ro = 360.0*(drand48 () - 0.5);149 double Phi = 2.0*(drand48 () - 0.5);157 double Ro = 360.0*(drand48_cnt() - 0.5); 158 double Phi = 2.0*(drand48_cnt() - 0.5); 150 159 double Do = DEG_RAD * asin(Phi); 151 160 152 double uR = 0.5*(drand48() - 0.5); 153 double uD = 0.5*(drand48() - 0.5); 154 155 double plx; 161 double uR = 0.5*(drand48_cnt() - 0.5); 162 double uD = 0.5*(drand48_cnt() - 0.5); 163 164 double plx = 0.5*(drand48_cnt() + 0.0); 165 156 166 switch (mode) { 157 167 case FIT_POS_IRLS: … … 166 176 case FIT_PLX_BOOT: 167 177 case FIT_PLX_NOCLIP: 168 plx = 0.5*(drand48() + 0.0);169 178 break; 170 179 default: … … 172 181 } 173 182 183 int Npoints = N_POINTS*drand48_cnt(); // minimum of 0 points 184 int Noutliers_max = Npoints * F_OUTLIERS; 185 int Noutliers = Noutliers_max*drand48_cnt(); // minimum of 0 points 186 174 187 // generate a single fake star with N_POINTS real points and N_OUTLIERS bad points 175 mkstar (fitStats, Ro, Do, uR, uD, plx, N _POINTS, N_OUTLIERS);188 mkstar (fitStats, Ro, Do, uR, uD, plx, Npoints, Noutliers); 176 189 177 190 double Tmean, Trange, parRange; … … 266 279 Ro, Do, uR, uD, plx, 267 280 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); 281 if (i > 0) { 282 // fprintf (stderr, ".\n"); 283 } 268 284 269 285 continue; … … 291 307 292 308 // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined) 293 points[i].T = MJD_MIN_YRS + drand48 ()*(MJD_MAX_YRS - MJD_MIN_YRS);309 points[i].T = MJD_MIN_YRS + drand48_cnt()*(MJD_MAX_YRS - MJD_MIN_YRS); 294 310 295 311 double pR, pD; … … 310 326 311 327 // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined) 312 points[i].T = MJD_MIN_YRS + drand48 ()*(MJD_MAX_YRS - MJD_MIN_YRS);328 points[i].T = MJD_MIN_YRS + drand48_cnt()*(MJD_MAX_YRS - MJD_MIN_YRS); 313 329 314 330 double pR, pD; 315 331 ParFactor (&pR, &pD, Ro, Do, points[i].T); 316 332 317 double dR = OUT_ERROR*(drand48 () - 0.5) + RA_BIAS;318 double dD = OUT_ERROR*(drand48 () - 0.5) + DEC_BIAS;333 double dR = OUT_ERROR*(drand48_cnt() - 0.5) + RA_BIAS; 334 double dD = OUT_ERROR*(drand48_cnt() - 0.5) + DEC_BIAS; 319 335 320 336 points[i].R = Ro + (dR + plx*pR + uR*(points[i].T - MJD_REF_YRS))/3600/dcos(Do); -
trunk/Ohana/src/relastro/test/mana.sh
r39370 r39371 131 131 132 132 data $1 133 #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 26133 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 134 134 # 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 135 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 22135 # 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 136 136 137 137 foreach f R D uD uR P … … 167 167 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 168 168 169 end 170 171 macro checkstats 172 if ($0 != 2) 173 echo "USAGE: checkstats (file)" 174 break 175 end 176 177 teststats $1 178 179 delete dRbin oRbin 180 for i 0 100 181 subset tmp = dRf if (Npts == $i) 182 vstat tmp -q 183 concat $SIGMA dRbin 184 concat $MEAN oRbin 185 end 186 187 create n 0 dRbin[] 188 set dRmod = 0.010 / sqrt (n) 189 lim n -0.001 0.015; clear; box; plot n dRbin 190 plot -c red -pt 7 n dRmod 169 191 end 170 192
Note:
See TracChangeset
for help on using the changeset viewer.
