IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39783


Ignore:
Timestamp:
Oct 20, 2016, 12:11:12 PM (10 years ago)
Author:
eugene
Message:

clean up rules for fitting the objects with par, pm, or pos

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/czw_branch/20160809/Ohana/src/relastro/src/UpdateObjects.c

    r39747 r39783  
    55# define PAR_MIN_NPTS 7
    66# define PAR_MIN_NPTS_BOOT 6
    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
     7# define PM_MIN_NPTS 5
     8# define PM_MIN_NPTS_BOOT 4
     9# define POS_MIN_NPTS 4
     10# define POS_MIN_NPTS_BOOT 4
    1111
    1212typedef enum {
    … …  
    183183    if (Trange < PM_DT_MIN) {
    184184      // not enough baseline for proper motion, only set mean position
    185       mode = FIT_AVERAGE;
    186185      goto justPosition;
    187186    }
    188187    if (parRange < PAR_FACTOR_MIN) {
    189188      // not enough parallax factor range, skip parallax
    190       mode = FIT_PM_ONLY;
    191189      goto skipParallax;
    192190    }
    193191    if (fitStats->Npoints < PAR_MIN_NPTS) {
    194192      // not enough data, skip parallax
    195       mode = FIT_PM_ONLY;
    196193      goto skipParallax;
    197194    }
    198195
    199196    // we are going to use the IRLS analysis to calculate the mean solution and the masking
    200     // then run N_BOOTSTRAP_SAMPLES to measure the errors
    201     fitStats->fitdataPar->getError = !N_BOOTSTRAP_SAMPLES;
     197    // then run N_BOOTSTRAP_SAMPLES to measure the errors.  We first measure the OLS error,
     198    // and choose the max of the OLS and bootstrap errors
     199    fitStats->fitdataPar->getError = TRUE;
    202200    if (USE_IRLS) {
    203201      if (!FitPMandPar_IRLS (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints)) {
    204         mode = FIT_PM_ONLY;
    205202        goto skipParallax;
    206203      }
    207204    } else {
    208205      if (!FitPMandPar_Basic (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints)) {
    209         mode = FIT_PM_ONLY;
    210206        goto skipParallax;
    211207      }
    … …  
    217213      fitStats->Nfit = 0;
    218214      int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
    219       if (Nnomask < PAR_MIN_NPTS_BOOT) {
    220         // if we do not have enough points to assess parallax error, we cannot do the bootstrap analysis.
    221         mode = FIT_PM_ONLY;
    222         goto skipParallax;
    223       }
     215
     216      // if we do not have enough points to assess parallax error, skip the bootstrap analysis.
     217      // (I think the test above means we never do this skip)
     218      if (Nnomask < PAR_MIN_NPTS_BOOT) goto skipParallaxBootstrap;
     219
    224220      for (k = 0; k < fitStats->NfitAlloc; k++) {
    225221        BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
    … …  
    227223        fitStats->Nfit ++;
    228224      }
     225
    229226      // These calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above)
    230227      // this call expects the fitted parameters to have the formal error set: it will apply the
    … …  
    237234    }
    238235
     236  skipParallaxBootstrap:
    239237    // project Ro, Do back to RA,DEC
    240238    XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &fitStats->coords);
    … …  
    255253    valid = valid && (fabs(fitPar.uD) < 4.0);
    256254    valid = valid && (fabs(fitPar.p) < 2.0);
    257     if (!valid) {
    258       mode = FIT_PM_ONLY;
    259     } else {
     255    if (valid) {
    260256      average[0].flags |= ID_OBJ_USE_PAR;
    261257    }
    … …  
    267263  if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) {
    268264    if (Trange < PM_DT_MIN) {
    269       mode = FIT_AVERAGE;
    270265      goto justPosition;
    271266    }
    272267    if (fitStats->Npoints < PM_MIN_NPTS) {
    273       mode = FIT_AVERAGE;
    274268      goto justPosition;
    275269    }
    276270
    277271    // if we have fitted (and accepted) a parallax model, get the best pm fit and chisq
    278     // given the set of points (mask is respected)
    279     if (average[0].flags & ID_OBJ_USE_PAR) {
     272    // given the set of points (mask is respected).  Alternatively, if we do not request
     273    // IRLS fitting, the just use OLS fitting (skip bootstrap)
     274    if ((average[0].flags & ID_OBJ_USE_PAR) || !USE_IRLS) {
    280275      if (!FitPM_Basic (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
    281276        average[0].flags |= ID_OBJ_BAD_PM;
    282277        goto justPosition;
    283278      }
    284     } else {
    285       fitStats->fitdataPM->getError = !N_BOOTSTRAP_SAMPLES;
    286       if (USE_IRLS) {
    287         if (!FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
    288           mode = FIT_AVERAGE;
    289           goto justPosition;
    290         }
    291       } else {
    292         if (!FitPM_Basic (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
    293           mode = FIT_AVERAGE;
    294           goto justPosition;
    295         }
    296       }
     279      goto skipProperMotionBootstrap;
     280    }
     281
     282    // we are going to use the IRLS analysis to calculate the mean solution and the masking
     283    // then run N_BOOTSTRAP_SAMPLES to measure the errors.  We first measure the OLS error,
     284    // and choose the max of the OLS and bootstrap errors
     285    fitStats->fitdataPM->getError = TRUE;
     286    if (!FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
     287      goto justPosition;
     288    }
    297289
    298290    // in the fits above, we have saved the formal error for the unmasked points. 
    299291    // if we do not have enough points for bootstrap, we will keep those errors
    300       if (N_BOOTSTRAP_SAMPLES && (fitPM.Nfit >= PM_MIN_NPTS_BOOT)) {
    301         fitStats->Nfit = 0;
    302         int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
    303         if (Nnomask < PM_MIN_NPTS_BOOT) {
    304           mode = FIT_AVERAGE;
    305           goto justPosition;
    306         }
    307         for (k = 0; k < fitStats->NfitAlloc; k++) {
    308           BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
    309           if (!FitPM_Basic (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, Nnomask)) continue;
    310           fitStats->Nfit ++;
    311         }
     292    if (N_BOOTSTRAP_SAMPLES && (fitPM.Nfit >= PM_MIN_NPTS_BOOT)) {
     293      fitStats->Nfit = 0;
     294      int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
     295       
     296      // if we do not have enough points to assess p.m. error, skip the bootstrap analysis.
     297      // (I think the test above means we never do this skip)
     298      if (Nnomask < PM_MIN_NPTS_BOOT) goto skipProperMotionBootstrap;
     299
     300      for (k = 0; k < fitStats->NfitAlloc; k++) {
     301        BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
     302        if (!FitPM_Basic (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, Nnomask)) continue;
     303        fitStats->Nfit ++;
     304      }
     305
    312306      // These calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above)
    313307      // this call expects the fitted parameters to have the formal error set: it will apply the
    314308      // max of the formal and bootstrap errors
    315         BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
    316         BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
    317         BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
    318         BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
    319       }
    320     }
    321 
     309      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     310      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     311      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
     312      BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
     313    }
     314
     315skipProperMotionBootstrap:
    322316    // project Ro, Do back to RA,DEC
    323317    XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &fitStats->coords);
    … …  
    326320    fitStats->Npm ++;
    327321
     322    // XXX a hard-wired hack...
    328323    // unless there is a clear problems (below) with the proper-motion fit or we have a parallax fit, we will use pm fit
    329324    int valid = TRUE;
    … …  
    335330    valid = valid && (fabs(fitPM.uD) < 4.0);
    336331    if (!valid) {
    337       mode = FIT_AVERAGE;
    338332      average[0].flags |= ID_OBJ_BAD_PM;
    339333    } else {
    … …  
    349343    // if we only have one point, this is silly...
    350344
     345    // set the proper motion (to the galaxy model or average value, if desired; else to 0,0)
    351346    FitAstromResultSetPM (&fitPos, 1, average);
    352     // fprintf (stderr, "fit 1: %f %f : %f %f\n", fitPos.Ro, fitPos.Do, fitPos.uR, fitPos.uD);
    353     if (average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_PM)) {
     347
     348    // if we already have a valid fit (pm or par), use OLS to fit the position
     349    // alternatively, if we do not request IRLS, use OLS
     350    // alternatively, if we do not have enough points, use OLS
     351    if ((average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_PM)) || (fitStats->Npoints < POS_MIN_NPTS) || !USE_IRLS) {
    354352      if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
    355         // if this fails, stick with the PM and/or PAR fit from above
     353        // if this fails, stick with the PM and/or PAR fit from above, or use a single value
    356354        goto doneWithFit;
    357355      }
    358     } else {
    359       if (fitStats->Npoints < POS_MIN_NPTS) {
    360         // I will not try to outlier-reject, just calculate the weighted average
    361         if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
    362           // if we have tried this, we have not masked any points; we will find a single unmasked point below
    363           goto doneWithFit;
    364         }
     356      // if we have not already gotten a good fit, use this fit
     357      if (!(average[0].flags & (ID_OBJ_USE_PAR | ID_OBJ_USE_PM))) {
    365358        average[0].flags |= ID_OBJ_USE_AVE;
    366359        average[0].flags |= ID_OBJ_RAW_AVE;
    367         goto useBasic;
    368       }
    369       if (USE_IRLS) {
    370         if (!FitPosPMfixed_IRLS (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
    371           // if the above fails, we need to clear the masks and try again below
    372           FitPointsClearMasks (fitStats->points, fitStats->Npoints);
    373           // just calculate the weighted average
    374           if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
    375             // if this fails, find a single unmasked point below
    376             goto doneWithFit;
    377           }
    378           average[0].flags |= ID_OBJ_USE_AVE;
    379           average[0].flags |= ID_OBJ_RAW_AVE;
    380           goto useBasic;
    381         }
    382       } else {
    383         if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
    384           // if this fails, find a single unmasked point below
    385           goto doneWithFit;
    386         }
    387         average[0].flags |= ID_OBJ_USE_AVE;
    388         average[0].flags |= ID_OBJ_RAW_AVE;
    389         goto useBasic;
    390       }
    391       // fprintf (stderr, "fit 2: %f %f : %f %f\n", fitPos.Ro, fitPos.Do, fitPos.uR, fitPos.uD);
    392       if (N_BOOTSTRAP_SAMPLES) {
    393         fitStats->Nfit = 0;
    394         int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
    395         if (Nnomask < POS_MIN_NPTS_BOOT) {
    396           // XXX if the above fails, accept the formal error from unmasked points?
    397           FitPointsClearMasks (fitStats->points, fitStats->Npoints);
    398           if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
    399             // if this fails, find a single unmasked point below
    400             goto doneWithFit;
    401           }
    402           average[0].flags |= ID_OBJ_USE_AVE;
    403           average[0].flags |= ID_OBJ_RAW_AVE;
    404           goto useBasic;
    405         }
    406         FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average);
    407         for (k = 0; k < fitStats->NfitAlloc; k++) {
    408           BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
    409           if (!FitPosPMfixed_Basic (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, Nnomask)) continue;
    410           fitStats->Nfit ++;
    411         }
    412         BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
    413         BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     360      }
     361      goto useBasic;
     362    }
     363   
     364    // try the IRLS fitting, otherwise give up and use OLS
     365    if (!FitPosPMfixed_IRLS (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
     366      // if the above fails, we need to clear the masks and try again below
     367      FitPointsClearMasks (fitStats->points, fitStats->Npoints);
     368      // just calculate the weighted average
     369      if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
     370        // if this fails, find a single unmasked point below
     371        goto doneWithFit;
    414372      }
    415373      average[0].flags |= ID_OBJ_USE_AVE;
     374      average[0].flags |= ID_OBJ_RAW_AVE;
     375      goto useBasic;
     376    }
     377    average[0].flags |= ID_OBJ_USE_AVE;
     378
     379    if (N_BOOTSTRAP_SAMPLES && (fitPos.Nfit >= POS_MIN_NPTS_BOOT)) {
     380      fitStats->Nfit = 0;
     381      int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
     382     
     383      if (Nnomask < POS_MIN_NPTS_BOOT) goto useBasic;
     384     
     385      FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average);
     386      for (k = 0; k < fitStats->NfitAlloc; k++) {
     387        BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
     388        if (!FitPosPMfixed_Basic (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, Nnomask)) continue;
     389        fitStats->Nfit ++;
     390      }
     391     
     392      // These calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above)
     393      // this call expects the fitted parameters to have the formal error set: it will apply the
     394      // max of the formal and bootstrap errors
     395      BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     396      BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
    416397    }
    417398
    … …  
    580561  average[0].dP         = fit.dp; // parallax error in arcsec
    581562
     563
    582564  average[0].ChiSqAve   = fitPos.chisq;
    583565  average[0].ChiSqPM    = fitPM.chisq;
    … …  
    873855  int i;
    874856
     857  myAbort ("this should not be called anymore");
     858
    875859  // add up the chi square for the fit
    876860  double chisq = 0.0;
Note: See TracChangeset for help on using the changeset viewer.