IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39371


Ignore:
Timestamp:
Feb 23, 2016, 10:55:44 AM (10 years ago)
Author:
eugene
Message:

updates to UpdateObjects: fix rules for number of detections (+ associated tests)

Location:
trunk/Ohana/src/relastro
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • trunk/Ohana/src/relastro/include/relastro.h

    r39365 r39371  
    310310int    SHOW_PARAMS;
    311311char   STATMODE[32];
    312 int    POS_TOOFEW;
    313 int    PM_TOOFEW;
     312// int    POS_TOOFEW;
     313// int    PM_TOOFEW;
    314314double PM_DT_MIN;
    315315double PAR_FACTOR_MIN;
     
    759759int FitPMandPar_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
    760760int FitPosPMfixed_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
     761int FitPosPMfixed_Single (FitAstromResult *fit, FitAstromPoint *points, int Npoints);
    761762
    762763double MedianAbsDeviation(FitAstromPoint *points, int Npoints);
  • trunk/Ohana/src/relastro/src/ConfigInit.c

    r39225 r39371  
    2828
    2929  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; 
    3232  if (!ScanConfig (config, "PAR_FACTOR_MIN",         "%lf", 0, &PAR_FACTOR_MIN))   PAR_FACTOR_MIN = 0.2;
    3333  if (!ScanConfig (config, "RELASTRO_MAP_NX",        "%d",  0, &NX_MAP))           NX_MAP = 5;
  • trunk/Ohana/src/relastro/src/FitPM.c

    r39370 r39371  
    2020  }
    2121
    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;
    2624  return TRUE;
    2725}
  • trunk/Ohana/src/relastro/src/FitPMandPar.c

    r39370 r39371  
    2020  }
    2121
    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;
    2625  return TRUE;
    2726}
     
    331330  // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms)
    332331  fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms);
    333 
    334 
    335332  return TRUE;
    336333}
  • trunk/Ohana/src/relastro/src/FitPosPMfixed.c

    r39370 r39371  
    99# define WEIGHT_THRESHOLD 0.3
    1010
     11int 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
    1139// fit->uR,uD are used in the fit
    1240int FitPosPMfixed_Basic (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
     
    2048  }
    2149
    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  }
    2659  return TRUE;
    2760}
  • trunk/Ohana/src/relastro/src/UpdateObjects.c

    r39370 r39371  
    11# 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
    311
    412typedef enum {
     
    160168      goto skipParallax;
    161169    }
    162     if (fitStats->Npoints <= PAR_TOOFEW) {
     170    if (fitStats->Npoints < PAR_MIN_NPTS) {
    163171      // not enough data, skip parallax
    164172      mode = FIT_PM_ONLY;
     
    177185      fitStats->Nfit = 0;
    178186      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.
    181189        mode = FIT_PM_ONLY;
    182190        goto skipParallax;
     
    229237      goto justPosition;
    230238    }
    231     if (fitStats->Npoints <= PM_TOOFEW) {
     239    if (fitStats->Npoints < PM_MIN_NPTS) {
    232240      mode = FIT_AVERAGE;
    233241      goto justPosition;
     
    249257        fitStats->Nfit = 0;
    250258        int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
    251         if (Nnomask < 5) {
     259        if (Nnomask < PM_MIN_NPTS_BOOT) {
    252260          mode = FIT_AVERAGE;
    253261          goto justPosition;
     
    297305    FitAstromResultSetPM (&fitPos, 1, average);
    298306    // 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)) {
    300308      if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
    301309        goto doneWithFit;
    302310      }
    303311    } 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      }
    304321      if (!FitPosPMfixed_IRLS (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
    305322        goto doneWithFit;
     
    309326        fitStats->Nfit = 0;
    310327        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;
    313335        }
    314336        FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average);
     
    324346    }
    325347
     348  useBasic:
    326349    // project Ro, Do back to RA,DEC
    327350    XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords);
     
    333356
    334357 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
    335373  // update the bit flags of which points were used
    336374  for (k = 0; k < fitStats->Npoints; k++) {
  • trunk/Ohana/src/relastro/src/fitobj2.c

    r39370 r39371  
    11# include "relastro.h"
     2
     3static int Nrand = 0;
     4double drand48_cnt () {
     5  double value = drand48();
     6  Nrand ++;
     7  return value;
     8}
    29
    310Catalog *mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad);
     
    109116    long A = now.tv_sec + now.tv_usec * 1000000;
    110117    srand48(A);
    111     // srand48(1);
     118    srand48(1);
    112119  }
    113120
     
    121128  for (i = 0; i < Nstars; i++) {
    122129
    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);
    125132    double Do  = DEG_RAD * asin(Phi);
    126133
    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
    134143
    135144    // generate a single fake star with N_POINTS real points and N_OUTLIERS bad points
     
    140149    double Tmean = ohana_sec_to_mjd(catalog->average->Tmean);
    141150    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",
    143152             Npoints, Noutliers, Ro, Do, uR, uD, plx,
    144153             catalog->average->R,  catalog->average->D,  catalog->average->uR,  catalog->average->uD,  catalog->average->P, Tyear,
    145154             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    }
    147159
    148160    continue;
     
    187199  for (i = 0; i < Npoints; i++) {
    188200    // 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);
    190202    double Tmjd   = 365.25*Tyears + MJD_J2000;
    191203
     
    214226  for (i = Npoints; i < Ntotal; i++) {
    215227    // 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);
    217229    double Tmjd   = 365.25*Tyears + MJD_J2000;
    218230
     
    222234    ParFactor (&pR, &pD, Ro, Do, Tyears);
    223235
    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);
    226238
    227239    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  
    11# include "relastro.h"
     2
     3static int Nrand = 0;
     4double drand48_cnt () {
     5  double value = drand48();
     6  Nrand ++;
     7  return value;
     8}
    29
    310int mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad);
     
    4148  int N_STARS     = 3000;
    4249  int N_POINTS    =  100;
    43   int N_OUTLIERS  =   10;
     50  // int N_OUTLIERS  =   10;
    4451  int N_BOOTSTRAP =  100;
     52  float F_OUTLIERS  =  0.1;
    4553
    4654  int N;
     
    5563    remove_argument (N, &argc, argv);
    5664  }
    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  // }
    6270  if ((N = get_argument (argc, argv, "-Nbootstrap"))) {
    6371    remove_argument (N, &argc, argv);
     
    132140  ohana_gaussdev_init ();
    133141
    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);
    135144  fitStats->fitdataPM->getChisq  = TRUE;
    136145  fitStats->fitdataPM->getError  = TRUE;
     
    146155  for (i = 0; i < Nstars; i++) {
    147156
    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);
    150159    double Do  = DEG_RAD * asin(Phi);
    151160
    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
    156166    switch (mode) {
    157167      case FIT_POS_IRLS:
     
    166176      case FIT_PLX_BOOT:
    167177      case FIT_PLX_NOCLIP:
    168         plx = 0.5*(drand48() + 0.0);
    169178        break;
    170179      default:
     
    172181    }
    173182
     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
    174187    // 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);
    176189
    177190    double Tmean, Trange, parRange;
     
    266279             Ro, Do, uR, uD, plx,
    267280             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    }
    268284
    269285    continue;
     
    291307
    292308    // 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);
    294310   
    295311    double pR, pD;
     
    310326
    311327    // 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);
    313329   
    314330    double pR, pD;
    315331    ParFactor (&pR, &pD, Ro, Do, points[i].T);
    316332
    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;
    319335
    320336    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  
    131131
    132132  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 26
     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 26
    134134  # 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 22
     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 22
    136136
    137137  foreach f R D uD uR P
     
    167167  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
    168168
     169end
     170
     171macro 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
    169191end
    170192
Note: See TracChangeset for help on using the changeset viewer.