IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39246 for trunk


Ignore:
Timestamp:
Dec 11, 2015, 9:12:57 AM (11 years ago)
Author:
eugene
Message:

adding IRLS versions of fitpm and fitplx

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

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relastro/Makefile

    r39241 r39246  
    3333$(SRC)/FitPM_IRLS.$(ARCH).o          \
    3434$(SRC)/FitPMandPar.$(ARCH).o         \
     35$(SRC)/FitPMandPar_IRLS.$(ARCH).o               \
    3536$(SRC)/FitPosPMfixed.$(ARCH).o       \
     37$(SRC)/BootstrapOps.$(ARCH).o        \
    3638$(SRC)/ImageOps.$(ARCH).o            \
    3739$(SRC)/MosaicOps.$(ARCH).o           \
    … …  
    108110$(SRC)/FitPM_IRLS.$(ARCH).o          \
    109111$(SRC)/FitPMandPar.$(ARCH).o         \
     112$(SRC)/FitPMandPar_IRLS.$(ARCH).o               \
    110113$(SRC)/FitPosPMfixed.$(ARCH).o       \
     114$(SRC)/BootstrapOps.$(ARCH).o        \
    111115$(SRC)/ImageOps.$(ARCH).o            \
    112116$(SRC)/MosaicOps.$(ARCH).o           \
    … …  
    172176$(SRC)/FitPMandPar.$(ARCH).o               \
    173177$(SRC)/FitPMandPar_IRLS.$(ARCH).o               \
     178$(SRC)/BootstrapOps.$(ARCH).o        \
    174179$(SRC)/mkpolyterm.$(ARCH).o            \
    175180$(SRC)/fitpoly.$(ARCH).o

  • trunk/Ohana/src/relastro/include/relastro.h

    r39241 r39246  
    138138  double *Beta_prev;
    139139  int Nterms;
     140  int getChisq;
     141  int getError;
    140142} FitAstromData;
    141143
    … …  
    173175  FitAstromPoint *points;
    174176  FitAstromPoint *sample;
     177  FitAstromPoint *nomask;
    175178  int Npoints;
    176179  int NpointsAlloc;
    … …  
    732735int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode);
    733736int BootstrapResample (FitAstromPoint *sample, FitAstromPoint *points, int Npoints);
     737int BootstrapSaveUnmasked (FitAstromPoint *nomask, FitAstromPoint *points, int Npoints);
     738
    734739int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog);
    735740int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange);

  • trunk/Ohana/src/relastro/src/FitAstromOps.c

    r39241 r39246  
    8787  fitStats->points = NULL;
    8888  fitStats->sample = NULL;
     89  fitStats->nomask = NULL;
    8990  fitStats->Npoints = 0;
    9091  fitStats->NpointsAlloc = Nmax;
    … …  
    9293    ALLOCATE (fitStats->points, FitAstromPoint, Nmax);
    9394    ALLOCATE (fitStats->sample, FitAstromPoint, Nmax);
     95    ALLOCATE (fitStats->nomask, FitAstromPoint, Nmax);
    9496  }
    9597
    … …  
    139141  FREE (fitStats->points);
    140142  FREE (fitStats->sample);
     143  FREE (fitStats->nomask);
    141144
    142145  FitAstromDataFree (fitStats->fitdataPos);
    … …  
    160163  ALLOCATE (fit->Beta_prev, double, Nterms);
    161164  fit->Nterms = Nterms;
     165
     166  fit->getChisq = TRUE;
     167  fit->getError = TRUE;
    162168
    163169  return fit;
    … …  
    203209  object->u      = 0.0;
    204210
    205   object->mask   = 0;
     211  object->mask   = 0; // keep point if mask == 0
    206212  return;
    207213}

  • trunk/Ohana/src/relastro/src/FitPM_IRLS.c

    r39241 r39246  
    1111  int i,j;
    1212
    13   int p   = 4;
    14   int n   = 2 * Npoints;
    15   int dof = n - p;
    16  
    17   // data->A,B,Cov,Beta,Beta_prev are allocated outside by FitAstromDataInit()
    18   // points->Wx,Wy,rx,ry,u are elements of FitAstromPoint
     13  int Ndof = 2 * Npoints - data->Nterms;
    1914 
    2015  // Convert the measurement errors into initial weights.
    2116  for (i = 0; i < Npoints; i++) {
    22     // points[i].Wx = 1 / points[i].dX;
    23     // points[i].Wy = 1 / points[i].dY;
    2417    points[i].Wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1 / SQ(points[i].dX);
    2518    points[i].Wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1 / SQ(points[i].dY);
    … …  
    3932    sigma_ols += SQ(points[i].rx) + SQ(points[i].ry);
    4033  }
    41   sigma_ols = sqrt(sigma_ols / dof);
     34  sigma_ols = sqrt(sigma_ols / (float)Ndof);
    4235
    4336  // Save OLS covariance and Beta (solution vector, which is actually also saved in fit)
    … …  
    10194  }
    10295
    103   // this section calculates the formal error on the weighted fit using the covariance values
     96  // calculate the weight thresholds to mask the bad points:
    10497  double Sum_Wx = 0.0;
    10598  double Sum_Wy = 0.0;
    106   if (1) {
     99  for (i = 0; i < Npoints; i++) {
     100    points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
     101    points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
     102   
     103    Sum_Wx += points[i].Wx;
     104    Sum_Wy += points[i].Wy;
     105  }
     106  double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints);
     107  double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints);
     108
     109  // set a mask (which can be used by the bootstrap resampling analysis)
     110  for (i = 0; i < Npoints; i++) {
     111    // keep if either is above threshold?
     112    // drop if either is below threshold?
     113    // points are marked as keep by default
     114    if ((points[i].Wx < WxThreshold) || (points[i].Wy < WyThreshold)) {
     115      points[i].mask = 1; // keep point if mask == 0
     116    }
     117  }
     118
     119  // this section calculates the formal error on the weighted fit using the covariance values
     120  // NOTE EAM: in tests (fitpm.c), they seem to be too large by a factor of ~5.37
     121  if (data->getError) {
    107122    double ax = 0.0, ay = 0.0;
    108123    double bx = 0.0, by = 0.0;
    109     double lambda = 0.0;
     124
    110125    for (i = 0; i < Npoints; i++) {
    111       points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
    112       points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
    113    
    114126      ax += dpsi_cauchy(points[i].rx / points[i].dX);
    115127      ay += dpsi_cauchy(points[i].ry / points[i].dY);
    … …  
    117129      bx += SQ(points[i].Wx);
    118130      by += SQ(points[i].Wy);
    119 
    120       Sum_Wx += points[i].Wx;
    121       Sum_Wy += points[i].Wy;
    122131    }
    123132    ax /= 1.0 * Npoints;  // mean(psi_dot(r))
    124133    ay /= 1.0 * Npoints;
    125     bx /= 1.0 * (Npoints - p); // mean(psi^2(r)) * (N / (N-p))
    126     by /= 1.0 * (Npoints - p);
    127  
    128     double sigma_robust_x = lambda * sqrt(bx) * sigma_hat * 2.385 / ax;
    129     double sigma_robust_y = lambda * sqrt(by) * sigma_hat * 2.385 / ay;
     134    bx /= 1.0 * (Npoints - data->Nterms); // mean(psi^2(r)) * (N / (N-p))
     135    by /= 1.0 * (Npoints - data->Nterms);
     136 
     137    double lambda_x = 1.0 + (data->Nterms / Npoints) * (1 - ax) / ax;
     138    double lambda_y = 1.0 + (data->Nterms / Npoints) * (1 - ay) / ay;
     139 
     140    double sigma_robust_x = lambda_x * sqrt(bx) * sigma_hat * 2.385 / ax;
     141    double sigma_robust_y = lambda_y * sqrt(by) * sigma_hat * 2.385 / ay;
    130142
    131143    // This is actually sigma^2, as that's the factor in the covariance (dumouchel 4.1)
    132     double sigma_final_x  = MAX(SQ(sigma_robust_x), (n * SQ(sigma_robust_x) + SQ(p * sigma_ols)) / (n + SQ(p)));
    133     double sigma_final_y  = MAX(SQ(sigma_robust_y), (n * SQ(sigma_robust_y) + SQ(p * sigma_ols)) / (n + SQ(p)));
     144    double sigma_final_x  = MAX(SQ(sigma_robust_x), (2 * Npoints * SQ(sigma_robust_x) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms)));
     145    double sigma_final_y  = MAX(SQ(sigma_robust_y), (2 * Npoints * SQ(sigma_robust_y) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms)));
    134146
    135147    fit[0].dRo = sqrt(data->Cov[0][0]);
    … …  
    144156  }
    145157
    146   // set a mask (which can be used by the bootstrap resampling analysis)
    147   double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints);
    148   double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints);
    149 
    150   for (i = 0; i < Npoints; i++) {
    151     // keep if either is above threshold?
    152     // drop if either is below threshold?
    153     // points are marked as keep by default
    154     if ((points[i].Wx < WxThreshold) || (points[i].Wy < WyThreshold)) {
    155       points[i].mask = 1;
    156     }
    157   }
    158 
    159158  // (optionally) add up the chi square for the fit, only counting the unmasked points
    160159  double chisq = 0.0;
    … …  
    162161  for (i = 0; i < Npoints; i++) {
    163162    if (points[i].mask) continue;
     163    fit[0].Nfit ++;
    164164     
    165     double Xf = fit[0].Ro + fit[0].uR*points[i].T;
    166     double Yf = fit[0].Do + fit[0].uD*points[i].T;
    167     double wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dX);
    168     double wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dY);
    169     chisq += SQ(points[i].X - Xf) * wx;
    170     chisq += SQ(points[i].Y - Yf) * wy;
    171     fit[0].Nfit ++;
    172   }
    173    
    174   // the reduced chisq is divided by (Ndof = 2*Nfit - 4)
    175   fit[0].chisq = chisq / (2.0*fit[0].Nfit - 4.0);
     165    if (data->getChisq) {
     166      double Xf = fit[0].Ro + fit[0].uR*points[i].T;
     167      double Yf = fit[0].Do + fit[0].uD*points[i].T;
     168      double wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dX);
     169      double wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dY);
     170      chisq += SQ(points[i].X - Xf) * wx;
     171      chisq += SQ(points[i].Y - Yf) * wy;
     172    }
     173  }
     174   
     175  // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms)
     176  fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms);
    176177
    177178  return (TRUE);

  • trunk/Ohana/src/relastro/src/FitPMandPar.c

    r38986 r39246  
    1515
    1616  for (i = 0; i < Npoints; i++) {
     17    if (points[i].mask) continue; // respect the mask if set
     18    fit->Nfit ++;
     19
    1720    /* handle case where dX or dY = 0.0 */
    1821    wx = 1.0 / SQ(points[i].dX);
    … …  
    8992  fit->dp  = sqrt(data->A[4][4]);
    9093 
    91   fit->Nfit = Npoints;
    92 
    9394  return (TRUE);
    9495}

  • trunk/Ohana/src/relastro/src/FitPMandPar_IRLS.c

    r39241 r39246  
    1111  int i,j;
    1212
    13   int p   = 5;
    14   int n   = 2 * Npoints;
    15   int dof = n - p;
    16  
    17   // data->A,B,Cov,Beta,Beta_prev are allocated outside by FitAstromDataInit()
    18   // points->Wx,Wy,rx,ry,u are elements of FitAstromPoint
     13  int Ndof = 2 * Npoints - data->Nterms;
    1914 
    2015  // Convert the measurement errors into initial weights.
    2116  for (i = 0; i < Npoints; i++) {
    22     // points[i].Wx = 1 / points[i].dX;
    23     // points[i].Wy = 1 / points[i].dY;
    2417    points[i].Wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1 / SQ(points[i].dX);
    2518    points[i].Wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1 / SQ(points[i].dY);
    … …  
    3932    sigma_ols += SQ(points[i].rx) + SQ(points[i].ry);
    4033  }
    41   sigma_ols = sqrt(sigma_ols / dof);
    42   // sigma_ols = sqrt(sigma_ols / (Npts - 5.0));
     34  sigma_ols = sqrt(sigma_ols / (float)Ndof);
    4335 
    4436  // Save OLS covariance and Beta (solution vector, which is actually also saved in fit)
    … …  
    9991  }
    10092  if (!converged) {
    101     myAbort ("raise a warning on non-convergence");
    102   }
    103 
    104   // this section calculates the formal error on the weighted fit using the covariance values
     93    fprintf (stderr, "raise a warning on non-convergence\n");
     94  }
     95
     96  // calculate the weight thresholds to mask the bad points:
    10597  double Sum_Wx = 0.0;
    10698  double Sum_Wy = 0.0;
    107   if (1) {
     99  for (i = 0; i < Npoints; i++) {
     100    points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
     101    points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
     102   
     103    Sum_Wx += points[i].Wx;
     104    Sum_Wy += points[i].Wy;
     105  }
     106  double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints);
     107  double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints);
     108
     109  // set a mask (which can be used by the bootstrap resampling analysis)
     110  for (i = 0; i < Npoints; i++) {
     111    // keep if either is above threshold?
     112    // drop if either is below threshold?
     113    // points are marked as keep by default
     114    if ((points[i].Wx < WxThreshold) || (points[i].Wy < WyThreshold)) {
     115      points[i].mask = 1; // keep point if mask == 0
     116    }
     117  }
     118
     119  // this section calculates the formal error on the weighted fit using the covariance values
     120  // NOTE EAM: in tests (fitpm.c), they seem to be too large by a factor of ~5.37
     121  if (data->getError) {
    108122    double ax = 0.0, ay = 0.0;
    109123    double bx = 0.0, by = 0.0;
    110     double lambda = 0.0;
     124
    111125    for (i = 0; i < Npoints; i++) {
    112       points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
    113       points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
    114    
    115126      ax += dpsi_cauchy(points[i].rx / points[i].dX);
    116127      ay += dpsi_cauchy(points[i].ry / points[i].dY);
    … …  
    118129      bx += SQ(points[i].Wx);
    119130      by += SQ(points[i].Wy);
    120 
    121       Sum_Wx += points[i].Wx;
    122       Sum_Wy += points[i].Wy;
    123131    }
    124132    ax /= 1.0 * Npoints;  // mean(psi_dot(r))
    125133    ay /= 1.0 * Npoints;
    126     bx /= 1.0 * (Npoints - p); // mean(psi^2(r)) * (N / (N-p))
    127     by /= 1.0 * (Npoints - p);
    128  
    129     double sigma_robust_x = lambda * sqrt(bx) * sigma_hat * 2.385 / ax;
    130     double sigma_robust_y = lambda * sqrt(by) * sigma_hat * 2.385 / ay;
     134    bx /= 1.0 * (Npoints - data->Nterms); // mean(psi^2(r)) * (N / (N-p))
     135    by /= 1.0 * (Npoints - data->Nterms);
     136 
     137    double lambda_x = 1.0 + (data->Nterms / Npoints) * (1 - ax) / ax;
     138    double lambda_y = 1.0 + (data->Nterms / Npoints) * (1 - ay) / ay;
     139 
     140    double sigma_robust_x = lambda_x * sqrt(bx) * sigma_hat * 2.385 / ax;
     141    double sigma_robust_y = lambda_y * sqrt(by) * sigma_hat * 2.385 / ay;
    131142
    132143    // This is actually sigma^2, as that's the factor in the covariance (dumouchel 4.1)
    133     double sigma_final_x  = MAX(SQ(sigma_robust_x), (n * SQ(sigma_robust_x) + SQ(p * sigma_ols)) / (n + SQ(p)));
    134     double sigma_final_y  = MAX(SQ(sigma_robust_y), (n * SQ(sigma_robust_y) + SQ(p * sigma_ols)) / (n + SQ(p)));
     144    double sigma_final_x  = MAX(SQ(sigma_robust_x), (2 * Npoints * SQ(sigma_robust_x) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms)));
     145    double sigma_final_y  = MAX(SQ(sigma_robust_y), (2 * Npoints * SQ(sigma_robust_y) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms)));
    135146
    136147    fit[0].dRo = sqrt(data->Cov[0][0]);
    … …  
    147158  }
    148159
    149   // set a mask (which can be used by the bootstrap resampling analysis)
    150   double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints);
    151   double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints);
    152 
    153   for (i = 0; i < Npoints; i++) {
    154     // keep if either is above threshold?
    155     // drop if either is below threshold?
    156     // points are marked as keep by default
    157     if ((points[i].Wx < WxThreshold) || (points[i].Wy < WyThreshold)) {
    158       points[i].mask = 1;
    159     }
    160   }
    161 
    162160  // (optionally) add up the chi square for the fit, only counting the unmasked points
    163161  double chisq = 0.0;
    … …  
    165163  for (i = 0; i < Npoints; i++) {
    166164    if (points[i].mask) continue;
     165    fit[0].Nfit ++;
    167166     
    168     double Xf = fit[0].Ro + fit[0].uR*points[i].T + fit[0].p*points[i].pR;
    169     double Yf = fit[0].Do + fit[0].uD*points[i].T + fit[0].p*points[i].pD;
    170     double wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dX);
    171     double wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dY);
    172     chisq += SQ(points[i].X - Xf) * wx;
    173     chisq += SQ(points[i].Y - Yf) * wy;
    174     fit[0].Nfit ++;
    175   }
    176    
    177   // the reduced chisq is divided by (Ndof = 2*Nfit - 4)
    178   fit[0].chisq = chisq / (2.0*fit[0].Nfit - 5.0);
     167    if (data->getChisq) {
     168      double Xf = fit[0].Ro + fit[0].uR*points[i].T + fit[0].p*points[i].pR;
     169      double Yf = fit[0].Do + fit[0].uD*points[i].T + fit[0].p*points[i].pD;
     170      double wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dX);
     171      double wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dY);
     172      chisq += SQ(points[i].X - Xf) * wx;
     173      chisq += SQ(points[i].Y - Yf) * wy;
     174    } 
     175  }
     176   
     177  // the reduced chisq is divided by Ndof = (2*Nfit - Nterms)
     178  fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms);
    179179 
    180180  return (TRUE);

  • trunk/Ohana/src/relastro/src/UpdateObjects.c

    r39239 r39246  
    140140  // use a fit with fewer parameters.
    141141
     142  // fit the parallax + proper-motion model
     143  // NOTE : we only fit PAR if we have already fitted for proper motion. if we do not fit PM or we fail
     144  // to fit PM, we do not attempt PAR.  thus failure to fit PAR falls back to PM-only
     145  if (mode == FIT_PM_AND_PAR) {
     146    if (Trange < PM_DT_MIN) {
     147      mode = FIT_AVERAGE;
     148      goto justPosition;
     149    }
     150    if (parRange < PAR_FACTOR_MIN) {
     151      mode = FIT_PM_ONLY;
     152      goto skipParallax;
     153    }
     154    if (fitStats->Npoints <= PAR_TOOFEW) {
     155      mode = FIT_PM_ONLY;
     156      goto skipParallax;
     157    }
     158
     159    // we are going to use the IRLS analysis to calculate the mean solution and the masking
     160    // then run N_BOOTSTRAP_SAMPLES to measure the errors
     161    fitStats->fitdataPar->getError = (N_BOOTSTRAP_SAMPLES < 3);
     162    FitPMandPar_IRLS (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints, VERBOSE);
     163
     164    if (N_BOOTSTRAP_SAMPLES >= 3) {
     165      fitStats->Nfit = 0;
     166      int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
     167      if (Nnomask < 5) {
     168        myAbort ("handle this case, please");
     169      }
     170      for (k = 0; k < fitStats->NfitAlloc; k++) {
     171        BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
     172        FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask);
     173        fitStats->Nfit ++;
     174      }
     175      // these calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above)
     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    }
     182
     183    // project Ro, Do back to RA,DEC
     184    XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &fitStats->coords);
     185    if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
     186
     187    average[0].flags |= ID_STAR_FIT_PAR;
     188    fitStats->Npar ++;
     189
     190    // XXX a hard-wired hack...
     191    // unless there is a clear problems (below) with the parallax fit, we will use it
     192    if ((fabs(fitPar.uR) > 4.0) || (fabs(fitPar.uD) > 4.0)) {
     193      mode = FIT_PM_ONLY;
     194    } else {
     195      average[0].flags |= ID_STAR_USE_PAR;
     196    }
     197  }       
     198
     199skipParallax:
     200
    142201  // *** first fit for the proper motion (skip fit if Trange or Npts is too small) ***
    143202  if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) {
    … …  
    151210    }
    152211
    153     if (fitStats->NfitAlloc == 1) {
    154       // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling:
    155       if (1) {
    156         FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints);
    157       } else {
    158         FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints, VERBOSE);
     212    if (average[0].flags & ID_STAR_USE_PAR) {
     213      if (!FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
     214        myAbort ("oops");
    159215      }
    160216    } else {
    161       fitStats->Nfit = 0;
    162       for (k = 0; k < fitStats->NfitAlloc; k++) {
    163         BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
    164         // if (!FitPM_IRLS (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints, VERBOSE)) continue;
    165         if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue;
    166         fitStats->Nfit ++;
    167       }
    168       BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
    169       BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
    170       BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
    171       BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
    172     }
    173     FitAstromSetChisq (&fitPM, fitStats->points, fitStats->Npoints, FIT_PM_ONLY);
     217      fitStats->fitdataPM->getError = (N_BOOTSTRAP_SAMPLES < 3);
     218      FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints, VERBOSE);
     219
     220      if (N_BOOTSTRAP_SAMPLES >= 3) {
     221        fitStats->Nfit = 0;
     222        int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
     223        if (Nnomask < 5) {
     224          myAbort ("handle this case, please");
     225        }
     226        for (k = 0; k < fitStats->NfitAlloc; k++) {
     227          BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
     228          if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue;
     229          fitStats->Nfit ++;
     230        }
     231        BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     232        BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     233        BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
     234        BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
     235      }
     236    }
    174237
    175238    // project Ro, Do back to RA,DEC
    … …  
    177240    if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
    178241
    179     fitPM.p  = fitPM.dp  = 0.0;
    180242    average[0].flags |= ID_STAR_FIT_PM;
    181243    fitStats->Npm ++;
    182244
    183245    // XXX a hard-wired hack...
    184     if ((fabs(fitPM.uR) > 2.0) || (fabs(fitPM.uD) > 2.0)) {
     246    // unless there is a clear problems (below) with the proper-motion fit or we have a parallax fit, we will use pm fit
     247    if ((fabs(fitPM.uR) > 4.0) || (fabs(fitPM.uD) > 4.0)) {
    185248      mode = FIT_AVERAGE;
    186249      average[0].flags |= ID_STAR_BAD_PM;
     250    } else {
     251      if (!(average[0].flags & ID_STAR_USE_PAR)) {
     252        average[0].flags |= ID_STAR_USE_PM;
     253      }
    187254    }
    188255  }
    189256 
    190   // fit the parallax + proper-motion model
    191   // NOTE : we only fit PAR if we have already fitted for proper motion. if we do not fit PM or we fail
    192   // to fit PM, we do not attempt PAR.  thus failure to fit PAR falls back to PM-only
    193   if (mode == FIT_PM_AND_PAR) {
    194     if (parRange < PAR_FACTOR_MIN) {
    195       mode = FIT_PM_ONLY;
    196       goto justPosition;
    197     }
    198     if (fitStats->Npoints <= PAR_TOOFEW) {
    199       mode = FIT_PM_ONLY;
    200       goto justPosition;
    201     }
    202 
    203     if (fitStats->NfitAlloc == 1) {
    204       // if N_BOOTSTRAP_SAMPLES = 1, no bootstrap resampling:
    205       // FitPMandPar_IRLS (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints);
    206       FitPMandPar (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints);
    207     } else {
    208       fitStats->Nfit = 0;
    209       for (k = 0; k < fitStats->NfitAlloc; k++) {
    210         BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
    211         // FitPMandPar_IRLS (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, fitStats->Npoints);
    212         FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, fitStats->Npoints);
    213         fitStats->Nfit ++;
    214       }
    215       BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
    216       BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
    217       BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
    218       BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
    219       BootstrapRobustStats (&fitPar, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX);
    220     }
    221     FitAstromSetChisq (&fitPar, fitStats->points, fitStats->Npoints, FIT_PM_AND_PAR);
    222 
    223     // project Ro, Do back to RA,DEC
    224     XY_to_RD (&fitPar.Ro, &fitPar.Do, fitPar.Ro, fitPar.Do, &fitStats->coords);
    225     if (fabs(fitPar.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
    226 
    227     average[0].flags |= ID_STAR_FIT_PAR;
    228     fitStats->Npar ++;
    229 
    230     // XXX a hard-wired hack...
    231     if ((fabs(fitPar.uR) > 2.0) || (fabs(fitPar.uD) > 2.0)) {
    232       mode = FIT_PM_ONLY;
    233     }
    234   }       
    235 
    236257justPosition:
    237258  {
    … …  
    239260    // if we only have one point, this is silly...
    240261   
     262    // XXX I need to rethink here : use a median for the position or weighted average? use IRLS anyway?
     263
    241264    if (fitStats->NfitAlloc == 1) {
    242265      FitAstromResultSetPM (&fitPos, 1, average);
    … …  
    638661}
    639662
    640 int BootstrapResample (FitAstromPoint *sample, FitAstromPoint *points, int Npoints) {
    641   int i;
    642 
    643   // I need to draw Npoints random entries from 'points' with replacement:
    644   for (i = 0; i < Npoints; i++) {
    645     int N = Npoints * drand48();
    646     sample[i] = points[N];
    647   }
    648   return TRUE;
    649 }
    650 
    651 // calculate mean and sigma points for the 5 fit parameter
    652 int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode) {
    653 
    654   // generate a histogram for the selected element
    655   double *values = NULL;
    656   ALLOCATE (values, double, Nfit);
    657  
    658   int i;
    659 
    660   for (i = 0; i < Nfit; i++) {
    661     switch (mode) {
    662       case FIT_RESULT_RA:
    663         values[i] = fit[i].Ro;
    664         break;
    665       case FIT_RESULT_DEC:
    666         values[i] = fit[i].Do;
    667         break;
    668       case FIT_RESULT_uR:
    669         values[i] = fit[i].uR;
    670         break;
    671       case FIT_RESULT_uD:
    672         values[i] = fit[i].uD;
    673         break;
    674       case FIT_RESULT_PLX:
    675         values[i] = fit[i].p;
    676         break;
    677       default:
    678         myAbort ("invalid option");
    679     }
    680   }
    681 
    682   dsort (values, Nfit);
    683 
    684   double median;
    685   if (Nfit % 2) {
    686     int Ncenter = Nfit / 2;
    687     median = values[Ncenter];
    688   } else {
    689     int Ncenter = Nfit / 2 - 1;
    690     median = 0.5*(values[Ncenter] + values[Ncenter + 1]);
    691   }
    692 
    693   double Slo = VectorFractionInterpolate (values, 0.158655, Nfit);
    694   double Shi = VectorFractionInterpolate (values, 0.841345, Nfit);
    695   double sigma = (Shi - Slo) / 2.0;
    696 
    697   switch (mode) {
    698     case FIT_RESULT_RA:
    699       result->Ro = median;
    700       result->dRo = sigma;
    701       break;
    702     case FIT_RESULT_DEC:
    703       result->Do = median;
    704       result->dDo = sigma;
    705       break;
    706     case FIT_RESULT_uR:
    707       result->uR = median;
    708       result->duR = sigma;
    709       break;
    710     case FIT_RESULT_uD:
    711       result->uD = median;
    712       result->duD = sigma;
    713       break;
    714     case FIT_RESULT_PLX:
    715       result->p = median;
    716       result->dp = sigma;
    717       break;
    718     default:
    719       myAbort ("invalid option");
    720   }
    721 
    722   return TRUE;
    723 }
    724 
    725 double VectorFractionInterpolate (double *values, float fraction, int Npts) {
    726 
    727   float F = fraction * Npts;
    728   int   N = fraction * Npts;
    729 
    730   if (N < 0        ) return NAN;
    731   if (N >= Npts - 2) return NAN;
    732 
    733   // interpolate between N,N+1
    734    
    735   double S = (F - N) * (values[N+1] - values[N]) + values[N];
    736   return S;
    737 }
    738 
    739663int FitAstromSetChisq (FitAstromResult *fit, FitAstromPoint *points, int Npoints, FitMode mode) {
    740664

  • trunk/Ohana/src/relastro/src/fitpm.c

    r39241 r39246  
    1212# define OUT_ERROR 0.500 /* arcsec */
    1313
    14 # define N_STARS 20000
     14# define N_STARS 3000
    1515# define N_POINTS 100
    16 # define N_OUTLIERS 10
     16# define N_OUTLIERS 100
     17# define N_BOOTSTRAP 100
    1718
    1819# define dcos(THETA) cos(RAD_DEG*THETA)
    … …  
    2425  FIT_PM_NOCLIP,
    2526  FIT_PLX_IRLS,
     27  FIT_PLX_BOOT,
    2628  FIT_PLX_NOCLIP,
    2729};
    … …  
    4648  if (!strcasecmp(argv[1], "plx-irls")) {
    4749    mode = FIT_PLX_IRLS;
     50  }
     51  if (!strcasecmp(argv[1], "plx-boot")) {
     52    mode = FIT_PLX_BOOT;
    4853  }
    4954
    … …  
    6267  gaussdev_init ();
    6368
    64   FitStats *fitStats = FitStatsInit (N_POINTS + N_OUTLIERS, 0);
    65 
    66   int i;
     69  FitStats *fitStats = FitStatsInit (N_POINTS + N_OUTLIERS, N_BOOTSTRAP);
     70  fitStats->fitdataPM->getChisq  = TRUE;
     71  fitStats->fitdataPM->getError  = TRUE;
     72  fitStats->fitdataPar->getChisq = TRUE;
     73  fitStats->fitdataPar->getError = TRUE;
     74  // NOTE: surprisingly, the chisq and error calculations only add ~2-3% to the exec time
     75
     76  int i, k;
    6777
    6878  int Nstars = N_STARS;
    … …  
    8393        break;
    8494      case FIT_PLX_IRLS:
     95      case FIT_PLX_BOOT:
    8596      case FIT_PLX_NOCLIP:
    8697        plx = 0.5*(drand48() + 0.0);
    … …  
    112123        FitPMandPar_IRLS (&fitPM, fitStats->fitdataPar, fitStats->points, fitStats->Npoints, 0);
    113124        break;
     125
     126      case FIT_PLX_BOOT:
     127        FitPMandPar_IRLS (&fitPM, fitStats->fitdataPar, fitStats->points, fitStats->Npoints, 0);
     128        fitStats->Nfit = 0;
     129        int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
     130        for (k = 0; k < fitStats->NfitAlloc; k++) {
     131          BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
     132          FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask);
     133          fitStats->Nfit ++;
     134        }
     135        // these calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above)
     136        BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     137        BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     138        BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
     139        BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
     140        BootstrapRobustStats (&fitPM, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX);
     141        break;
     142
    114143      default:
    115144        myAbort("oops");
    … …  
    118147    XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &fitStats->coords);
    119148
    120     fprintf (stdout, "%12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f\n",
     149    fprintf (stdout, "%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\n",
    121150             Ro, Do, uR, uD, plx,
    122              fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p, Tmean);
     151             fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p, Tmean, fitPM.dRo, fitPM.dDo, fitPM.duR, fitPM.duD, fitPM.dp, fitPM.Nfit, fitPM.chisq);
    123152
    124153    // ok (TRUE, "fitted star");
Note: See TracChangeset for help on using the changeset viewer.