IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38604


Ignore:
Timestamp:
Jul 19, 2015, 6:26:28 AM (11 years ago)
Author:
eugene
Message:

big re-work of UpdateObjects to allow bootstrap resampling, outlier rejection, and fixed-pm fitting

Location:
branches/eam_branches/ipp-20150625/Ohana/src/relastro
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20150625/Ohana/src/relastro/include/relastro.h

    r38599 r38604  
    2828
    2929typedef enum {TARGET_NONE, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS} FitTarget;
     30
     31typedef enum {
     32  FIT_RESULT_RA,
     33  FIT_RESULT_DEC,
     34  FIT_RESULT_uR,
     35  FIT_RESULT_uD,
     36  FIT_RESULT_PLX,
     37} FitAstromResultMode;
    3038
    3139typedef enum {
    … …  
    119127  double  p, dp;
    120128
    121   int getChisq;
    122129  double chisq;
    123130  int Nfit;
    … …  
    140147  double C_blue;
    141148  double C_red;
    142 } FitAstromPoints;
     149} FitAstromPoint;
    143150
    144151typedef struct {
    … …  
    149156  off_t Noffset;
    150157
     158  double *values;
    151159  FitAstromResult *fit; // use bootstrap resampling to generate Nfit fits to measure the stats
    152160  int Nfit;
    153161  int NfitAlloc;
    154162
    155   FitAstromPoints *points;
    156   FitAstromPoints *sample;
     163  FitAstromPoint *points;
     164  FitAstromPoint *sample;
    157165  int Npoints;
    158166  int NpointsAlloc;
    159167
    160   FitAstromData *fitPos;
    161   FitAstromData *fitPM;
    162   FitAstromData *fitPar;
     168  FitAstromData *fitdataPos;
     169  FitAstromData *fitdataPM;
     170  FitAstromData *fitdataPar;
    163171
    164172  Coords coords;
    … …  
    491499int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon, double *Radius);
    492500int ParFactor (double *pR, double *pD, double RA, double Dec, double Time, double Tmean);
    493 int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, int XVERB);
    494 int FitPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *pR, double *pD, int Npts);
    495 int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int XVERB);
     501int FitPM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
     502int FitPMandPar (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
     503int FitPosPMfixed (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
    496504
    497505Mosaic *getMosaicForImage (off_t N);
    … …  
    694702int client_logger_init (char *dirname);
    695703int client_logger_message (char *format,...);
     704
     705int FitAstromSetChisq (FitAstromResult *fit, FitAstromPoint *points, int Npoints, FitMode mode);
     706double VectorFractionInterpolate (double *values, float fraction, int Npts);
     707int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode);
     708int BootstrapResample (FitAstromPoints *sample, FitAstromPoints *points, int Npoints);
     709int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog);
     710int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange);
     711int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff);
     712int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats);
     713int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int cat, off_t measOff);
     714
     715void FitAstromResultInit (FitAstromResult *fit);
     716void FitAstromPointInit (FitAstromObject *object);
     717void FitAstromDataFree (FitAstromData *fit);
     718FitAstromData *FitAstromDataInit (int Nterms);
     719void FitStatsFree (FitStats *fitStats);
     720void FitStatsSum (FitStats *src, FitStats *tgt);
     721void FitStatsReset (FitStats *tgt);
     722FitStats *FitStatsInit (int Nmax, int Nboot);

  • branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitAstromOps.c

    r38599 r38604  
    1212  fitStats->Noffset = 0;
    1313
    14   // container to hold the fit results (Nboot > 1 for bootstrap resampling)
    15   fitStats->fit = NULL;
     14  fitStats->values = NULL; // pre-allocated array for median & robust sigma
     15  fitStats->fit = NULL; // container to hold the fit results (Nboot > 1 for bootstrap resampling)
    1616  fitStats->Nfit = 0;
    1717  fitStats->NfitAlloc = Nboot;
    1818  if (Nboot > 0) {
    1919    ALLOCATE (fitStats->fit, FitAstromResult, Nboot);
     20    ALLOCATE (fitStats->values, double, Nboot);
    2021  }
    2122
    … …  
    3132
    3233  // pre-allocated fit matrices for the 3 fit options
    33   fitStats->fitPos = NULL;
    34   fitStats->fitPM  = NULL;
    35   fitStats->fitPar = NULL;
     34  fitStats->fitdataPos = NULL;
     35  fitStats->fitdataPM  = NULL;
     36  fitStats->fitdataPar = NULL;
    3637
    3738  if (Nmax > 0) {
    38     fitStats->fitPos = FitAstromDataInit (2);
    39     fitStats->fitPM  = FitAstromDataInit (4);
    40     fitStats->fitPar = FitAstromDataInit (5);
     39    fitStats->fitdataPos = FitAstromDataInit (2);
     40    fitStats->fitdataPM  = FitAstromDataInit (4);
     41    fitStats->fitdataPar = FitAstromDataInit (5);
    4142  }
    4243
    … …  
    7273
    7374  FREE (fitStats->fit);
     75  FREE (fitStats->values);
    7476  FREE (fitStats->points);
    7577  FREE (fitStats->sample);
    7678
    77   FitAstromDataFree (fitStats->fitPos);
    78   FitAstromDataFree (fitStats->fitPM);
    79   FitAstromDataFree (fitStats->fitPar);
     79  FitAstromDataFree (fitStats->fitdataPos);
     80  FitAstromDataFree (fitStats->fitdataPM);
     81  FitAstromDataFree (fitStats->fitdataPar);
    8082
    8183  free (fitStats);
    … …  
    105107}
    106108
    107 } 
    108 
    109 void FitAstromObjectInit (FitAstromObject *object) {
     109void FitAstromPointInit (FitAstromObject *object) {
    110110  object->X      = 0.0;
    111111  object->Y      = 0.0;
    … …  
    125125}
    126126
    127 void freeObjectData () {
     127void FitAstromResultInit (FitAstromResult *fit) {
    128128
    129   free (R);
    130   free (D);
    131   free (T);
    132   free (X);
    133   free (Y);
     129  fit->Ro  = 0.0;
     130  fit->dRo = 0.0;
     131  fit->Do  = 0.0;
     132  fit->dDo = 0.0;
     133  fit->uR  = 0.0;
     134  fit->duR = 0.0;
     135  fit->uD  = 0.0;
     136  fit->duD = 0.0;
     137  fit->p   = 0.0;
     138  fit->dp  = 0.0;
    134139
    135   free (dR);
    136   free (dD);
    137   free (dT);
    138   free (dX);
    139   free (dY);
    140 
    141   free (pX);
    142   free (pY);
    143 
    144   free (C_blue);
    145   free (C_red);
    146 } 
    147 
     140  fit->chisq = NAN;
     141  fit->Nfit = 0;
     142 
     143  return;
     144}

  • branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPM.c

    r38599 r38604  
    6565  fit->Nfit = Npoints;
    6666
    67   // add up the chi square for the fit
    68   if (fit->getChisq) {
    69     double chisq = 0.0;
    70     for (i = 0; i < Npoints; i++) {
    71       double Xf = fit->Ro + fit->uR*points[i].T;
    72       double Yf = fit->Do + fit->uD*points[i].T;
    73       chisq += SQ(points[i].X - Xf) / SQ(points[i].dX);
    74       chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY);
    75     }
    76     // the reduced chisq is divided by (Ndof = 2*Npts - 4)
    77     fit->chisq = chisq / (2.0*Npts - 4.0);
    78   }
    79 
    8067  return (TRUE);
    8168}

  • branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPMandPar.c

    r38599 r38604  
    9090  fit->Nfit = Npoints;
    9191
    92   // add up the chi square for the fit
    93   if (fit->getChisq) {
    94     double chisq = 0.0;
    95     for (i = 0; i < Npoints; i++) {
    96       double Xf = fit->Ro + fit->uR*points[i].T + fit->p*points[i].pR;
    97       double Yf = fit->Do + fit->uD*points[i].T + fit->p*points[i].pD;
    98       chisq += SQ(points[i].X - Xf) / SQ(points[i].dX);
    99       chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY);
    100     }
    101     // the reduced chisq is divided by (Ndof = 2*Npts - 5)
    102     fit->chisq = chisq / (2.0*Npts - 5.0);
    103   }
    104 
    10592  return (TRUE);
    10693}

  • branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/FitPosPMfixed.c

    r38599 r38604  
    3737  fit->Nfit = Npoints;
    3838
    39   // add up the chi square for the fit
    40   if (fit->getChisq) {
    41     double chisq = 0.0;
    42     for (i = 0; i < Npoints; i++) {
    43       double Xf = fit->Ro + fit->uR*points[i].T;
    44       double Yf = fit->Do + fit->uD*points[i].T;
    45       chisq += SQ(points[i].X - Xf) / SQ(points[i].dX);
    46       chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY);
    47     }
    48     // the reduced chisq is divided by (Ndof = 2*Npts - 4)
    49     fit->chisq = chisq / (2.0*Npts - 4.0);
    50   }
    51 
    5239  return (TRUE);
    5340}

  • branches/eam_branches/ipp-20150625/Ohana/src/relastro/src/UpdateObjects.c

    r38599 r38604  
    22# define PAR_TOOFEW 5
    33# define NBOOT 100
    4 
    5 int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset);
    6 int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats);
    74
    85// This function operates on both Measure and MeasureTiny.  In the big stages, this should
    … …  
    7168  /* calculate the average value of R,D for a single star */
    7269
    73   PMFit fit;    memset (&fit,    0, sizeof(fit));
    74   PMFit fitAve; memset (&fitAve, 0, sizeof(fitAve)); fitAve.chisq = NAN;
    75   PMFit fitPM;  memset (&fitPM,  0, sizeof(fitPM));  fitPM.chisq  = NAN;
    76   PMFit fitPAR; memset (&fitPAR, 0, sizeof(fitPAR)); fitPAR.chisq = NAN;
     70  FitAstromResult fitPos, fitPM, fitPar;
     71  FitAstromResultInit (&fitPos);
     72  FitAstromResultInit (&fitPM);
     73  FitAstromResultInit (&fitPar);
    7774
    7875  // if we fail to fit the astrometry for some reason, we need to set/reset these
    … …  
    9693
    9794  // select the measurements to be used in this analysis
    98   FitAstromObject *points = UpdateObjects_SelectMeasures (average, secfilt, measure, measureBig, cat, measOff, &Npoints);
     95  UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, cat, measOff);
    9996
    10097  // if there are no exposure detections, use the stack position
    101   if (Npoints < 1) {
     98  if (fitStats->Npoints < 1) {
    10299    if (isfinite(average[0].Rstk) && isfinite(average[0].Dstk)) {
    103100      average[0].R  = average[0].Rstk;
    … …  
    111108
    112109  double Tmean, Trange, parRange;
    113   UpdateObjects_Project (points, Npoints, &Tmean, &Trange, &parRange);
     110  UpdateObjects_Project (fitStats->points, fitStats->Npoints, &Tmean, &Trange, &parRange);
    114111
    115112  // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
    … …  
    129126    }
    130127
     128    // XXX handle the case with only 1 or few points
     129
    131130    fitStats->Nfit = 0;
    132131    for (k = 0; k < fitStats->NfitAlloc; k++) {
    133       BootstrapResample (sample, points, Npoints);
    134       if (!FitPM (&fitStats->fit[k], fitdata, sample, Npoints)) continue;
     132      BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
     133      // XXX init fitStats->fit[k]? set fitStats->fit[k].getChisq?
     134      if (!FitPM (&fitStats->fit[k], fitStats->fitPM, fitStats->sample, fitStats->Npoints)) continue;
    135135      fitStats->Nfit ++;
    136136    }
    137     // XXX convert the fit distributions of uR,uD to (uRo, uDo) +/- (duRo, duDo)
    138 
    139     if (XVERB) fprintf (stderr, "fitted PM:  %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq);
     137    BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     138    BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     139    BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
     140    BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
     141
     142    FitAstromSetChisq (&fitPM, fitStats->points, fitStats->Npoints, FIT_PM_ONLY);
    140143
    141144    // project Ro, Do back to RA,DEC
    … …  
    171174    fitStats->Nfit = 0;
    172175    for (k = 0; k < fitStats->NfitAlloc; k++) {
    173       BootstrapResample (sample, points, Npoints);
    174       FitPMandPar (&fitStats->fit[k], fitsdata, sample, Npoints);
     176      BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
     177      FitPMandPar (&fitStats->fit[k], fitStats->fitPar, fitStats->sample, fitStats->Npoints);
    175178      fitStats->Nfit ++;
    176179    }
    177     // XXX convert the fit distributions of uR,uD to (uRo, uDo) +/- (duRo, duDo)
    178 
    179     if (XVERB) fprintf (stderr, "fitted PM+PAR:  %f - %f : %f %f : %f %f : %f %f : %f vs %f vs %f\n", Tmin, Tmax, fitPAR.Ro, fitPAR.Do, fitPAR.uR, fitPAR.uD, fitPAR.p, fitPAR.dp, fitPAR.chisq, fitPM.chisq, fitAve.chisq);
    180 
    181     XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords);
     180    BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     181    BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     182    BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
     183    BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
     184    BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX);
     185
     186    FitAstromSetChisq (&fitPar, fitStats->points, fitStats->Npoints, FIT_PM_AND_PAR);
     187
     188    XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &coords);
    182189    average[0].flags |= ID_STAR_FIT_PAR;
    183190    fitStats->Npar ++;
    184191
    185     if (fabs(fitPAR.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
     192    if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
    186193
    187194    // XXX a hard-wired hack...
    188     if ((fabs(fitPAR.uR) > 2.0) || (fabs(fitPAR.uD) > 2.0)) {
     195    if ((fabs(fitPar.uR) > 2.0) || (fabs(fitPar.uD) > 2.0)) {
    189196      mode = FIT_PM_ONLY;
    190197    }
    … …  
    198205    fitStats->Nfit = 0;
    199206    for (k = 0; k < fitStats->NfitAlloc; k++) {
    200       BootstrapResample (sample, points, Npoints);
    201       FitPosPMfixed (&fitStats->fit[k], fitsdata, sample, Npoints);
     207      BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
     208      FitPosPMfixed (&fitStats->fit[k], fitStats->fitPos, fitStats->sample, fitStats->Npoints);
    202209      fitStats->Nfit ++;
    203210    }
    204 
    205     if (XVERB) fprintf (stderr, "average: %f %f\n", fitAve.Ro, fitAve.Do);
     211    BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     212    BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     213    FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE);
    206214
    207215    // project Ro, Do back to RA,DEC
    208     XY_to_RD (&fitAve.Ro, &fitAve.Do, fitAve.Ro, fitAve.Do, &coords);
     216    XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &coords);
    209217    average[0].flags |= ID_STAR_FIT_AVE;
    210218    fitStats->Nave ++;
    211 
    212     // XXX fitAve.dRo = statsR.sigma;
    213     // XXX fitAve.dDo = statsD.sigma;
    214     // XXX
    215     // XXX fitAve.chisq = (N > 1) ? 0.5 * (statsR.chisq + statsD.chisq) : NAN;
    216     // XXX fitAve.Nfit = N;
    217     // XXX
    218     // XXX fitAve.uR = fitAve.duR = 0.0;
    219     // XXX fitAve.uD = fitAve.duD = 0.0;
    220     // XXX fitAve.p  = fitAve.dp  = 0.0;
    221219  }
    222220
    223221  // XXX update the bit flags of which points were used
    224   for (k = 0; k < Npoints; k++) {
    225     int Nm = points[k].measure;
     222  for (k = 0; k < fitStats->Npoints; k++) {
     223    int Nm = fitStats->points[k].measure;
    226224    measure[Nm].dbFlags |= ID_MEAS_USED_OBJ;
    227225    if (measureBig) { measureBig[Nm].dbFlags |= ID_MEAS_USED_OBJ; }
    228226  }
    229227
     228# if (0)
    230229  if (setRefColor) {
    231230    // need to reassign here if isfinite()
    … …  
    238237    average[0].refColorRed = colorMedian;
    239238  }
     239# endif
    240240
    241241  /* choose the result based on the chisq values */
    242242  // XXXX for now, just use the mode as the result:
    243243  int result = mode;
    244 
     244  FitAstromResult fit;
    245245  switch (result) {
    246246    case FIT_AVERAGE:
    247247      average[0].flags |= ID_STAR_USE_AVE;
    248       fit = fitAve;
     248      fit = fitPos;
    249249      break;
    250250    case FIT_PM_ONLY:
    … …  
    254254    case FIT_PM_AND_PAR:
    255255      average[0].flags |= ID_STAR_USE_PAR;
    256       fit = fitPAR;
     256      fit = fitPar;
    257257      break;
    258258  }
    … …  
    314314                      average[0].uR,
    315315                      average[0].uD,
    316                       fitAve.chisq, fitPM.chisq, fitPAR.chisq);
     316                      fitPos.chisq, fitPM.chisq, fitPar.chisq);
    317317
    318318  average[0].R          = fit.Ro; // RA in degrees
    … …  
    329329  average[0].dP         = fit.dp; // parallax error in arcsec
    330330
    331   average[0].ChiSqAve   = fitAve.chisq;
     331  average[0].ChiSqAve   = fitPos.chisq;
    332332  average[0].ChiSqPM    = fitPM.chisq;
    333   average[0].ChiSqPar   = fitPAR.chisq;
     333  average[0].ChiSqPar   = fitPar.chisq;
    334334
    335335  average[0].Tmean      = (Tmean * 86400 * 365.25) + T2000;
    … …  
    479479}
    480480
    481 int BootstrapResample (FitAstromPoints *sample, FitAstromPoints *points, int Npoints) {
    482   int i;
    483 
    484   // I need to draw Npoints random entries from 'points' with replacement:
    485   for (i = 0; i < Npoints; i++) {
    486     int N = Npoints * drand48();
    487     sample[i] = points[N];
    488   }
    489   return TRUE;
    490 }
    491 
    492 FitAstromObject *UpdateObjects_SelectMeasures (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff, int *npoints) {
    493 
    494   // XXX move this above, but just for re-working now:
    495   int Npoints = 0;
    496   int NPOINTS = average->Nmeasure;
    497   FitAstromObject *points = NULL;
    498   ALLOCATE (points, FitAstromObject, NPOINTS);
     481int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff) {
     482
     483  // I've already allocated fit->points (and fit->sample) with space for fit->NpointsAlloc entries
     484
     485  int Npoints = fit->Npoints = 0;
     486  FitAstromPoint *points = fit->points;
    499487
    500488  // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
    … …  
    551539    }
    552540
    553     FitAstromObjectInit (&points[Npoints]);
     541    FitAstromPointInit (&points[Npoints]);
    554542
    555543    points[Npoints].R = Ri;
    … …  
    586574    points[Npoints].measure = k;
    587575    Npoints++;
     576
     577    myAssert (Npoints <= fit->NpointsAlloc, "oops");
    588578  } // loop over measurements : average[0].Nmeasure
    589579
    590   *npoints = Npoints;
    591   return points;
    592 }
    593 
    594 int FitAstromObject_Project (FitAstromObject *points, int Npoints, double *Tmean, double *Trange, double *parRange) {
     580  fit->Npoints = Npoints;
     581  return TRUE;
     582}
     583
     584int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange) {
    595585
    596586  int k;
     587
     588  FitAstromPoint *points = fitStats->points;
    597589
    598590  // find Tmin & Tmax from the list of accepted measurements
    … …  
    632624
    633625  /* we need to do the fit in a locally linear space; choose a ref coordinate */
    634   coords.crval1 = points[0].R;
    635   coords.crval2 = points[0].D;
     626  fitStats->coords.crval1 = points[0].R;
     627  fitStats->coords.crval2 = points[0].D;
    636628
    637629  // project all of the R,D coordinates to a plane centered on this coordinate. set
    638630  // the times to be relative to Tmean
    639631  for (k = 0; k < Npoints; k++) {
    640     RD_to_XY (&points[k].X, &points[k].Y, points[k].R, points[k].D, &coords);
     632    RD_to_XY (&points[k].X, &points[k].Y, points[k].R, points[k].D, &fitStats->coords);
    641633    points[k].T -= Tmean;
    642634  }       
    … …  
    658650}
    659651
    660 /* fitting proper-motion and parallax:
    661 
    662    given a source at position r,d, at a time t, we need to calculate a vector (pr,pd)
    663 
    664    let x,y be the coordinate in the linearized frame with y parallel to DEC lines
    665 
    666    L,B are the ecliptic longitude and latitude of the object,
    667    dL and dB are the offsets in the L and B directions
    668 
    669    dL = sin(t - topp)
    670    dB = cos(t - topp)*sin(B)
    671 
    672    these need to be rotated to the R,D frame to yield pR,pD.  Then, the equation of motion
    673    for the source in the x,y frame is:
    674 
    675    x = Ro + uR * (t - to) + p * pR
    676    y = Do + uD * (t - to) + p * pD
    677 
    678    the unknowns in these equations are Ro, uR, Do, uD, and p
    679 
    680    XXX think through the concepts for the pole a bit better.  all objects near the pole
    681    move the same way with the same phase.  choose a projection center and define dL,dB relative
    682    to that center point coordinate system?
    683 
    684 */
     652int BootstrapResample (FitAstromPoints *sample, FitAstromPoints *points, int Npoints) {
     653  int i;
     654
     655  // I need to draw Npoints random entries from 'points' with replacement:
     656  for (i = 0; i < Npoints; i++) {
     657    int N = Npoints * drand48();
     658    sample[i] = points[N];
     659  }
     660  return TRUE;
     661}
     662
     663// calculate mean and sigma points for the 5 fit parameter
     664int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode) {
     665
     666  // generate a histogram for the selected element
     667  double *values = NULL;
     668  ALLOCATE (values, double, Nfit);
     669 
     670  for (i = 0; i < Nfit; i++) {
     671    switch (mode) {
     672      case FIT_RESULT_RA:
     673        values[i] = fit[i].Ro;
     674        break;
     675      case FIT_RESULT_DEC:
     676        values[i] = fit[i].Do;
     677        break;
     678      case FIT_RESULT_uR:
     679        values[i] = fit[i].uR;
     680        break;
     681      case FIT_RESULT_uD:
     682        values[i] = fit[i].uD;
     683        break;
     684      case FIT_RESULT_PLX:
     685        values[i] = fit[i].p;
     686        break;
     687      default:
     688        myAbort ("invalid option");
     689    }
     690  }
     691
     692  dsort (values, Nfit);
     693
     694  double median;
     695  if (Nfit % 2) {
     696    int Ncenter = Nfit / 2;
     697    median = values[Ncenter];
     698  } else {
     699    int Ncenter = Nfit / 2 - 1;
     700    median = 0.5*(values[Ncenter] + values[Ncenter + 1]);
     701  }
     702
     703  double Slo = VectorFractionInterpolate (values, 0.158655, Nfit);
     704  double Shi = VectorFractionInterpolate (values, 0.841345, Nfit);
     705  double sigma = (Shi - Slo) / 2.0;
     706
     707  switch (mode) {
     708    case FIT_RESULT_RA:
     709      result->Ro = median;
     710      result->dRo = sigma;
     711      break;
     712    case FIT_RESULT_DEC:
     713      result->Do = median;
     714      result->dDo = sigma;
     715      break;
     716    case FIT_RESULT_uR:
     717      result->uR = median;
     718      result->duR = sigma;
     719      break;
     720    case FIT_RESULT_uD:
     721      result->uD = median;
     722      result->duD = sigma;
     723      break;
     724    case FIT_RESULT_PLX:
     725      result->p = median;
     726      result->dp = sigma;
     727      break;
     728    default:
     729      myAbort ("invalid option");
     730  }
     731
     732  return TRUE;
     733}
     734
     735double VectorFractionInterpolate (double *values, float fraction, int Npts) {
     736
     737  float F = fraction * Npts;
     738  int   N = fraction * Npts;
     739
     740  if (N < 0        ) return NAN;
     741  if (N >= Npts - 2) return NAN;
     742
     743  // interpolate between N,N+1
     744   
     745  double S = (F - N) * (values[N+1] - values[N]) + values[N];
     746  return S;
     747}
     748
     749int FitAstromSetChisq (FitAstromResult *fit, FitAstromPoint *points, int Npoints, FitMode mode) {
     750
     751  // add up the chi square for the fit
     752  double chisq = 0.0;
     753  for (i = 0; i < Npoints; i++) {
     754    double Xf = fit->Ro + fit->uR*points[i].T + fit->p*points[i].pR;
     755    double Yf = fit->Do + fit->uD*points[i].T + fit->p*points[i].pD;
     756    chisq += SQ(points[i].X - Xf) / SQ(points[i].dX);
     757    chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY);
     758  }
     759  switch (mode) {
     760    case FIT_AVERAGE:
     761      fit->chisq = chisq / (2.0*Npts - 2.0);
     762      break;
     763    case FIT_PM_ONLY:
     764      fit->chisq = chisq / (2.0*Npts - 4.0);
     765      break;
     766    case FIT_PM_AND_PAR:
     767      fit->chisq = chisq / (2.0*Npts - 5.0);
     768      break;
     769    default:
     770      myAbort ("invalid mode");
     771  }
     772  return (TRUE);
     773}
Note: See TracChangeset for help on using the changeset viewer.