IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38655


Ignore:
Timestamp:
Aug 1, 2015, 10:57:47 AM (11 years ago)
Author:
eugene
Message:

tested bootstrap resample; handle fixed PM; update ref color

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

Legend:

Unmodified
Added
Removed

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

    r38612 r38655  
    294294int    CHIPORDER;
    295295int    CHIPMAP;
     296
     297int    N_BOOTSTRAP_SAMPLES;
    296298
    297299int MaxDensityUse;
    … …  
    710712int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog);
    711713int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange);
    712 int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff, int isStack);
     714int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int isStack);
    713715int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats);
    714716int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int cat, off_t measOff);
    … …  
    722724void FitStatsReset (FitStats *tgt);
    723725FitStats *FitStatsInit (int Nmax, int Nboot);
     726int FitAstromResultSetPM (FitAstromResult *fit, int Nfit, Average *average);

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

    r38604 r38655  
    11# include "relastro.h"
    22
     3// initial values of *fit are ignored
    34int FitPM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
    45

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

    r38612 r38655  
    11# include "relastro.h"
    22
     3// initial values of *fit are ignored
    34int FitPMandPar (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
    45

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

    r38471 r38655  
    185185    // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
    186186
     187    if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES);
     188
    187189    if (DCR_BLUE_COLOR_POS && DCR_BLUE_COLOR_NEG) {
    188190      strextend (&command, "-dcr-blue-color %s %s", DCR_BLUE_COLOR_POS, DCR_BLUE_COLOR_NEG);

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

    r38612 r38655  
    11# include "relastro.h"
    22# define PAR_TOOFEW 5
    3 # define NBOOT 100
    43
    54// This function operates on both Measure and MeasureTiny.  In the big stages, this should
    … …  
    1918  FitStats *sumStatsStack = FitStatsInit (0, 0);
    2019
    21   FitStats *fitStatsChips = FitStatsInit (NmeasureMax, NBOOT);
    22   FitStats *fitStatsStack = FitStatsInit (NmeasureMax, NBOOT);
     20  FitStats *fitStatsChips = FitStatsInit (NmeasureMax, N_BOOTSTRAP_SAMPLES);
     21  FitStats *fitStatsStack = FitStatsInit (NmeasureMax, N_BOOTSTRAP_SAMPLES);
    2322
    2423  int i;
    … …  
    3938      SecFilt *secfilt = &catalog[i].secfilt[j*Nsecfilt];
    4039
     40      int XVERB = FALSE;
     41      XVERB |= (average[0].objID == OBJ_ID_SRC) && (average[0].catID == CAT_ID_SRC);
     42      XVERB |= (average[0].objID == OBJ_ID_DST) && (average[0].catID == CAT_ID_DST);
     43      if ((OBJ_ID_SRC || OBJ_ID_DST) && !XVERB) continue;
     44
    4145      UpdateObjects_Stack(average, secfilt, measure, measureBig, Nsecfilt, fitStatsStack);
    4246      UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, fitStatsChips, i, m);
    … …  
    9094
    9195  // select the measurements to be used in this analysis
    92   UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, cat, measOff, FALSE);
     96  UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, FALSE);
    9397
    9498  // if there are no exposure detections, use the stack position
    … …  
    123127    }
    124128
    125     // XXX handle the case with only 1 or few points
    126 
    127     fitStats->Nfit = 0;
    128     for (k = 0; k < fitStats->NfitAlloc; k++) {
    129       BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
    130       if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue;
    131       fitStats->Nfit ++;
    132     }
    133     BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
    134     BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
    135     BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
    136     BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
    137 
     129    if (fitStats->NfitAlloc == 1) {
     130      // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling:
     131      FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints);
     132    } else {
     133      fitStats->Nfit = 0;
     134      for (k = 0; k < fitStats->NfitAlloc; k++) {
     135        BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
     136        if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue;
     137        fitStats->Nfit ++;
     138      }
     139      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     140      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     141      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
     142      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
     143    }
    138144    FitAstromSetChisq (&fitPM, fitStats->points, fitStats->Npoints, FIT_PM_ONLY);
    139145
    140146    // project Ro, Do back to RA,DEC
    141147    XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &fitStats->coords);
    142     if (XVERB) fprintf (stderr, "project: %f %f : %f %f : %f\n", fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p);
    143148    if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
    144     // XXX : does this make sense at 0,360 boundary?
    145149
    146150    fitPM.p  = fitPM.dp  = 0.0;
    … …  
    168172    }
    169173
    170     fitStats->Nfit = 0;
    171     for (k = 0; k < fitStats->NfitAlloc; k++) {
    172       BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
    173       FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, fitStats->Npoints);
    174       fitStats->Nfit ++;
    175     }
    176     BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
    177     BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
    178     BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
    179     BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
    180     BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX);
    181 
     174    if (fitStats->NfitAlloc == 1) {
     175      // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling:
     176      FitPMandPar (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints);
     177    } else {
     178      fitStats->Nfit = 0;
     179      for (k = 0; k < fitStats->NfitAlloc; k++) {
     180        BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
     181        FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, fitStats->Npoints);
     182        fitStats->Nfit ++;
     183      }
     184      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     185      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     186      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
     187      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
     188      BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX);
     189    }
    182190    FitAstromSetChisq (&fitPar, fitStats->points, fitStats->Npoints, FIT_PM_AND_PAR);
    183191
     192    // project Ro, Do back to RA,DEC
    184193    XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &fitStats->coords);
     194    if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
     195
    185196    average[0].flags |= ID_STAR_FIT_PAR;
    186197    fitStats->Npar ++;
    187 
    188     if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
    189198
    190199    // XXX a hard-wired hack...
    … …  
    199208    // if we only have one point, this is silly...
    200209   
    201     fitStats->Nfit = 0;
    202     for (k = 0; k < fitStats->NfitAlloc; k++) {
    203       BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
    204       FitPosPMfixed (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, fitStats->Npoints);
    205       fitStats->Nfit ++;
    206     }
    207     BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
    208     BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     210    if (fitStats->NfitAlloc == 1) {
     211      FitAstromResultSetPM (&fitPos, 1, average);
     212      FitPosPMfixed (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints);
     213    } else {
     214      fitStats->Nfit = 0;
     215      FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average);
     216      for (k = 0; k < fitStats->NfitAlloc; k++) {
     217        BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
     218        FitPosPMfixed (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, fitStats->Npoints);
     219        fitStats->Nfit ++;
     220      }
     221      BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     222      BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     223    }
    209224    FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE);
    210225
    … …  
    215230  }
    216231
    217   // XXX update the bit flags of which points were used
     232  // update the bit flags of which points were used
    218233  for (k = 0; k < fitStats->Npoints; k++) {
    219234    int Nm = fitStats->points[k].measure;
    … …  
    223238  }
    224239
    225 # if (0)
     240  // we can set the star reference-image color only if we have loaded the image data
    226241  int setRefColor = areImagesMatched();
    227242  if (setRefColor) {
     243    float *C_blue = NULL;
     244    float *C_red = NULL;
     245    ALLOCATE (C_blue, float, fitStats->Npoints);
     246    ALLOCATE (C_red, float, fitStats->Npoints);
     247
     248    int NcBlue = 0;
     249    int NcRed = 0;
     250
     251    for (k = 0; k < fitStats->Npoints; k++) {
     252      int Nm = fitStats->points[k].measure;
     253      float colorBlue = getColorBlue (measOff + Nm, cat);
     254      if (!isnan(colorBlue)) {
     255        C_blue[NcBlue] = colorBlue;
     256        NcBlue++;
     257      }
     258      float colorRed = getColorRed (measOff + Nm, cat);
     259      if (!isnan(colorRed)) {
     260        C_red[NcRed] = colorRed;
     261        NcRed++;
     262      }
     263    }
     264
    228265    // need to reassign here if isfinite()
    229266    float colorMedian;
    230     dsort (C_blue, NcBlue);
     267    fsort (C_blue, NcBlue);
    231268    colorMedian = (NcBlue > 0) ? C_blue[(int)(0.5*NcBlue)] : NAN;
    232269    average[0].refColorBlue = colorMedian;
    233     dsort (C_red, NcRed);
     270    fsort (C_red, NcRed);
    234271    colorMedian = (NcRed > 0) ? C_red[(int)(0.5*NcRed)] : NAN;
    235272    average[0].refColorRed = colorMedian;
    236   }
    237 # endif
     273
     274    free (C_blue);
     275    free (C_red);
     276  }
    238277
    239278  /* choose the result based on the chisq values */
    … …  
    362401
    363402  // select the measurements to be used in this analysis
    364   UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, 0, 0, TRUE);
     403  UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, TRUE);
    365404
    366405  // too few measurements for average position (require 2 values)
    … …  
    421460}
    422461
    423 int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int cat, off_t measOff, int isStack) {
     462int UpdateObjects_SelectMeasures (FitStats *fit, Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int isStack) {
    424463
    425464  // I've already allocated fit->points (and fit->sample) with space for fit->NpointsAlloc entries
    … …  
    512551
    513552    points[Npoints].dT = measure[k].dt;
    514 
    515     // XXX
    516     int setRefColor = 0;
    517     if (setRefColor) {
    518       points[Npoints].C_blue = getColorBlue (measOff+k, cat);
    519       points[Npoints].C_red  = getColorRed (measOff+k, cat);
    520     }
    521553
    522554    points[Npoints].measure = k;
    … …  
    720752      myAbort ("invalid mode");
    721753  }
     754  fit->Nfit = Npoints;
    722755  return (TRUE);
    723756}
     757
     758int FitAstromResultSetPM (FitAstromResult *fit, int Nfit, Average *average) {
     759
     760  int i;
     761
     762  if (USE_GALAXY_MODEL) {
     763    for (i = 0; i < Nfit; i++) {
     764      fit->uR = average->uRgal;
     765      fit->uD = average->uDgal;
     766    }
     767  } else {
     768    for (i = 0; i < Nfit; i++) {
     769      fit->uR = 0.0;
     770      fit->uD = 0.0;
     771    }
     772  }
     773
     774  return TRUE;
     775}

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

    r38441 r38655  
    271271    remove_argument (N, &argc, argv);
    272272    APPLY_OFFSETS = TRUE;
     273  }
     274
     275  N_BOOTSTRAP_SAMPLES = 1;
     276  if ((N = get_argument (argc, argv, "-bootstrap-samples"))) {
     277    remove_argument (N, &argc, argv);
     278    N_BOOTSTRAP_SAMPLES = atoi (argv[N]);
     279    remove_argument (N, &argc, argv);
     280    if ((N_BOOTSTRAP_SAMPLES < 1) || ((N_BOOTSTRAP_SAMPLES > 1) && (N_BOOTSTRAP_SAMPLES < 20))) {
     281      fprintf (stderr, "-bootstrap-samples must be either 1 (no sampling) or >= 20\n");
     282      exit (2);
     283    }
    273284  }
    274285
    … …  
    697708    remove_argument (N, &argc, argv);
    698709    MaxDensityUse = TRUE;
     710  }
     711
     712  N_BOOTSTRAP_SAMPLES = 1;
     713  if ((N = get_argument (argc, argv, "-bootstrap-samples"))) {
     714    remove_argument (N, &argc, argv);
     715    N_BOOTSTRAP_SAMPLES = atoi (argv[N]);
     716    remove_argument (N, &argc, argv);
     717    if ((N_BOOTSTRAP_SAMPLES < 1) || ((N_BOOTSTRAP_SAMPLES > 1) && (N_BOOTSTRAP_SAMPLES < 20))) {
     718      fprintf (stderr, "-bootstrap-samples must be either 1 (no sampling) or >= 20\n");
     719      exit (2);
     720    }
    699721  }
    700722

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

    r38471 r38655  
    176176    // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
    177177
    178     if (MinBadQF > 0.0)        { strextend (&command, "-min-bad-psfqf %f", MinBadQF); }
    179     if (MaxMeanOffset != 10.0) { strextend (&command, "-max-mean-offset  %f", MaxMeanOffset); }
     178    if (MinBadQF > 0.0)          strextend (&command, "-min-bad-psfqf %f", MinBadQF);
     179    if (MaxMeanOffset != 10.0)   strextend (&command, "-max-mean-offset  %f", MaxMeanOffset);
     180    if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES);
    180181
    181182    if (WHERE_A[0]) { strextend (&command, "-D WHERE_A \"%s\"", WHERE_A); }

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

    r38471 r38655  
    174174    // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
    175175
    176     if (MinBadQF > 0.0)        { strextend (&command, "-min-bad-psfqf %f", MinBadQF); }
    177     if (MaxMeanOffset != 10.0) { strextend (&command, "-max-mean-offset  %f", MaxMeanOffset); }
     176    if (MinBadQF > 0.0)          strextend (&command, "-min-bad-psfqf %f", MinBadQF);
     177    if (MaxMeanOffset != 10.0)   strextend (&command, "-max-mean-offset  %f", MaxMeanOffset);
     178    if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES);
    178179
    179180    if (TimeSelect) {

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

    r38471 r38655  
    172172    // XXX note that the above pass in the flag as decimal -- also note that args.c cannot handle 0xHEX values
    173173
    174     if (MinBadQF > 0.0)        { strextend (&command, "-min-bad-psfqf %f", MinBadQF); }
    175     if (MaxMeanOffset != 10.0) { strextend (&command, "-max-mean-offset  %f", MaxMeanOffset); }
     174    if (MinBadQF > 0.0)          strextend (&command, "-min-bad-psfqf %f", MinBadQF);
     175    if (MaxMeanOffset != 10.0)   strextend (&command, "-max-mean-offset  %f", MaxMeanOffset);
     176    if (N_BOOTSTRAP_SAMPLES > 1) strextend (&command, "-bootstrap-samples %d", N_BOOTSTRAP_SAMPLES);
    176177
    177178    if (TimeSelect) {
Note: See TracChangeset for help on using the changeset viewer.