IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39364 for trunk


Ignore:
Timestamp:
Feb 21, 2016, 6:00:34 AM (10 years ago)
Author:
eugene
Message:

consolidate irls and non-irls fitting; clean up rules for UpdateObjects

Location:
trunk/Ohana/src/relastro/src
Files:
1 added
11 edited

Legend:

Unmodified
Added
Removed

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

    r39247 r39364  
    2626}
    2727
    28 // calculate mean and sigma points for the 5 fit parameter
     28// calculate sigma values for the 5 fit parameters
    2929int BootstrapRobustStats (FitAstromResult *result, FitAstromResult *fit, int Nfit, int mode) {
    3030

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

    r39246 r39364  
    228228  fit->chisq = NAN;
    229229  fit->Nfit = 0;
     230  fit->converged = FALSE;
    230231 
    231232  return;

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

    r39287 r39364  
    11# include "relastro.h"
     2int FitPM_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
     3int FitPM_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
     4
     5// These should probably be tunable:
     6# define MAX_ITERATIONS 10
     7# define FIT_TOLERANCE 1e-4
     8# define FLT_TOLERANCE 1e-6
     9# define WEIGHT_THRESHOLD 0.3
    210
    311// initial values of *fit are ignored
    4 int FitPM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
     12int FitPM_Basic (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
    513
    614  int i;
    715
     16  // Convert the measurement errors into initial weights.
     17  for (i = 0; i < Npoints; i++) {
     18    points[i].Wx = 1.0 / SQ(points[i].dX);
     19    points[i].Wy = 1.0 / SQ(points[i].dY);
     20  }
     21
     22  int status = FitPM_MinChisq (fit, data, points, Npoints);
     23  if (!status) return FALSE;
     24
     25  FitPM_SetChisq (fit, data, points, Npoints);
     26  return TRUE;
     27}
     28
     29int FitPM_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
     30
     31  int i,j;
     32
     33  int Ndof = 2 * Npoints - data->Nterms;
     34 
     35  // Convert the measurement errors into initial weights.
     36  for (i = 0; i < Npoints; i++) {
     37    points[i].Wx = 1.0 / SQ(points[i].dX);
     38    points[i].Wy = 1.0 / SQ(points[i].dY);
     39  }
     40 
     41  // Solve OLS equation 
     42  if (!FitPM_MinChisq(fit, data, points, Npoints)) {
     43    return(FALSE);
     44  }
     45
     46  // Calculate r vector of residuals and least squares sigma
     47  double sigma_ols = 0.0;
     48  for (i = 0; i < Npoints; i++) {
     49    points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro);
     50    points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do);
     51    sigma_ols += SQ(points[i].rx) + SQ(points[i].ry);
     52  }
     53  sigma_ols = sqrt(sigma_ols / (float)Ndof);
     54
     55  // Save OLS covariance and Beta (solution vector, which is actually also saved in fit)
     56  for (i = 0; i < data->Nterms; i++) {
     57    for (j = 0; j < data->Nterms; j++) {
     58      data->Cov[i][j] = data->A[i][j];
     59    }
     60    data->Beta[i] = data->B[i][0];
     61  }
     62
     63  // Iteratively reweight and solve
     64  double sigma_hat = 0.0; // save for the error model
     65  int converged = FALSE;
     66  int iterations = 0;
     67
     68  // modify the weight based on the distance from the previous fit.  try up to MAX_ITERATIONS.
     69  // at the end "fit", has the last fit parameters
     70  for (iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations ++) {
     71    // Save Beta.
     72    for (i = 0; i < data->Nterms; i ++) {
     73      data->Beta_prev[i] = data->Beta[i];
     74    }
     75
     76    // Assign weights based on the deviation
     77    for (i = 0; i < Npoints; i++) {
     78      points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
     79      points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
     80    }   
     81
     82    // Solve with the new weights
     83    if (!FitPM_MinChisq(fit, data, points, Npoints)) {
     84
     85      // restore the last solution and break
     86      fit->Ro = data->Beta_prev[0];
     87      fit->uR = data->Beta_prev[1];
     88      fit->Do = data->Beta_prev[2];
     89      fit->uD = data->Beta_prev[3];
     90     
     91      // calculate the residuals:
     92      for (i = 0; i < Npoints; i++) {
     93        points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro);
     94        points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do);
     95        points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY));
     96      }
     97      sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745;
     98      break;
     99    }
     100
     101    // store the new Beta.
     102    for (i = 0; i < data->Nterms; i++) {
     103      data->Beta[i] = data->B[i][0];
     104    }
     105
     106    // calculate the residuals:
     107    for (i = 0; i < Npoints; i++) {
     108      points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro);
     109      points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do);
     110      points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY));
     111    }
     112    sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745;
     113   
     114    // Check convergence
     115    converged = TRUE;
     116    for (i = 0; i < data->Nterms; i++) {
     117      // if we are within FIT_TOLERANCE as a fractional error or FLT_TOLERANCE as an absolute error, we are good
     118      if ((fabs(data->Beta[i] - data->Beta_prev[i]) > FIT_TOLERANCE * fabs(data->Beta[i])) &&
     119          (fabs(data->Beta[i] - data->Beta_prev[i]) > FLT_TOLERANCE)) {
     120        converged = FALSE;
     121      }
     122    }
     123  }
     124  fit->converged = converged;
     125
     126  // calculate the weight thresholds to mask the bad points:
     127  double Sum_Wx = 0.0;
     128  double Sum_Wy = 0.0;
     129  for (i = 0; i < Npoints; i++) {
     130    points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
     131    points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
     132   
     133    Sum_Wx += points[i].Wx;
     134    Sum_Wy += points[i].Wy;
     135  }
     136  double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints);
     137  double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints);
     138
     139  // set a mask (which can be used by the bootstrap resampling analysis)
     140  for (i = 0; i < Npoints; i++) {
     141    // keep if either is above threshold?
     142    // drop if either is below threshold?
     143    // points are marked as keep by default
     144    if ((points[i].Wx < WxThreshold) || (points[i].Wy < WyThreshold)) {
     145      points[i].mask = 1; // keep point if mask == 0
     146    }
     147  }
     148
     149  // this section calculates the formal error on the weighted fit using the covariance values
     150  // NOTE EAM: in tests (fitpm.c), they seem to be too large by a factor of ~5.37
     151  if (data->getError) {
     152    double ax = 0.0, ay = 0.0;
     153    double bx = 0.0, by = 0.0;
     154
     155    for (i = 0; i < Npoints; i++) {
     156      ax += dpsi_cauchy(points[i].rx / points[i].dX);
     157      ay += dpsi_cauchy(points[i].ry / points[i].dY);
     158
     159      bx += SQ(points[i].Wx);
     160      by += SQ(points[i].Wy);
     161    }
     162    ax /= 1.0 * Npoints;  // mean(psi_dot(r))
     163    ay /= 1.0 * Npoints;
     164    bx /= 1.0 * (Npoints - data->Nterms); // mean(psi^2(r)) * (N / (N-p))
     165    by /= 1.0 * (Npoints - data->Nterms);
     166 
     167    double lambda_x = 1.0 + (data->Nterms / Npoints) * (1 - ax) / ax;
     168    double lambda_y = 1.0 + (data->Nterms / Npoints) * (1 - ay) / ay;
     169 
     170    double sigma_robust_x = lambda_x * sqrt(bx) * sigma_hat * 2.385 / ax;
     171    double sigma_robust_y = lambda_y * sqrt(by) * sigma_hat * 2.385 / ay;
     172
     173    // This is actually sigma^2, as that's the factor in the covariance (dumouchel 4.1)
     174    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)));
     175    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)));
     176
     177    fit[0].dRo = sqrt(data->Cov[0][0]);
     178    fit[0].duR = sqrt(data->Cov[1][1]);
     179    fit[0].dDo = sqrt(data->Cov[2][2]);
     180    fit[0].duD = sqrt(data->Cov[3][3]);
     181
     182    fit[0].dRo *= sigma_final_x;
     183    fit[0].duR *= sigma_final_x;
     184    fit[0].dDo *= sigma_final_y;
     185    fit[0].duD *= sigma_final_y;
     186  }
     187
     188  FitPM_SetChisq (fit, data, points, Npoints);
     189  return (TRUE);
     190}
     191
     192int FitPM_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
     193
     194  myAssert (data->Nterms == 4, "invalid fit arrays");
     195
     196  int i;
    8197  double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
    9 
    10   myAssert (data->Nterms == 4, "invalid fit arrays");
    11 
    12198  Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
    13199
    14200  for (i = 0; i < Npoints; i++) {
    15201    if (points[i].mask) continue; // respect the mask if set
    16     fit->Nfit ++;
    17 
    18     /* handle case where dX or dY = 0.0 */
    19     wx = 1.0 / SQ(points[i].dX);
    20     wy = 1.0 / SQ(points[i].dY);
     202
     203    wx = points[i].Wx;
     204    wy = points[i].Wy;
    21205
    22206    Wx += wx;
    23207    Wy += wy;
    24208
    25     Tx += points[i].T*wx;
    26     Ty += points[i].T*wy;
    27    
    28     Tx2 += SQ(points[i].T)*wx;
    29     Ty2 += SQ(points[i].T)*wy;
     209    double TWx = points[i].T*wx;
     210    double TWy = points[i].T*wy;
     211
     212    Tx += TWx;
     213    Ty += TWy;
     214   
     215    Tx2 += points[i].T*TWx;
     216    Ty2 += points[i].T*TWy;
    30217   
    31218    Xs += points[i].X*wx;
    32219    Ys += points[i].Y*wy;
    33220
    34     XT += points[i].X*points[i].T*wx;
    35     YT += points[i].Y*points[i].T*wy;
    36   }
    37 
     221    XT += points[i].X*TWx;
     222    YT += points[i].Y*TWy;
     223  }
     224
     225  // X^T W X
    38226  data->A[0][0] = Wx;
    39227  data->A[0][1] = Tx;
     228  data->A[0][2] = 0.0;
     229  data->A[0][3] = 0.0;
    40230
    41231  data->A[1][0] = Tx;
    42232  data->A[1][1] = Tx2;
    43 
     233  data->A[1][2] = 0.0;
     234  data->A[1][3] = 0.0;
     235
     236  data->A[2][0] = 0.0;
     237  data->A[2][1] = 0.0;
    44238  data->A[2][2] = Wy;
    45239  data->A[2][3] = Ty;
    46240
     241  data->A[3][0] = 0.0;
     242  data->A[3][1] = 0.0;
    47243  data->A[3][2] = Ty;
    48244  data->A[3][3] = Ty2;
    49245
     246  // X^T W Y
    50247  data->B[0][0] = Xs;
    51248  data->B[1][0] = XT;
    … …  
    53250  data->B[3][0] = YT;
    54251
    55   dgaussjordan (data->A, data->B, 4, 1);
     252  if (!dgaussjordan (data->A, data->B, 4, 1)) {
     253    return FALSE;
     254  }
    56255
    57256  fit->Ro = data->B[0][0];
    … …  
    60259  fit->uD = data->B[3][0];
    61260  fit->p  = 0.0;
    62  
     261
    63262  fit->dRo = sqrt(data->A[0][0]);
    64263  fit->duR = sqrt(data->A[1][1]);
    … …  
    66265  fit->duD = sqrt(data->A[3][3]);
    67266  fit->dp  = 0.0;
    68  
    69   return (TRUE);
    70 }
     267
     268  return TRUE;
     269}
     270
     271int FitPM_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
     272
     273  int i;
     274
     275  // count unmasked points and (optionally) add up the chi square for the fit
     276  double chisq = 0.0;
     277  fit[0].Nfit = 0;
     278  for (i = 0; i < Npoints; i++) {
     279    if (points[i].mask) continue;
     280    fit[0].Nfit ++;
     281     
     282    if (data->getChisq) {
     283      double Xf = fit[0].Ro + fit[0].uR*points[i].T;
     284      double Yf = fit[0].Do + fit[0].uD*points[i].T;
     285      double wx = 1.0 / SQ(points[i].dX);
     286      double wy = 1.0 / SQ(points[i].dY);
     287      chisq += SQ(points[i].X - Xf) * wx;
     288      chisq += SQ(points[i].Y - Yf) * wy;
     289    }
     290  }
     291   
     292  // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms)
     293  fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms);
     294
     295  return TRUE;
     296}

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

    r39246 r39364  
    2020 
    2121  // Solve OLS equation 
    22   if (!weighted_LS_PM(fit, data, points, Npoints, VERBOSE)) {
    23     myAbort ("handle failures, please!");
     22  if (!FitPM_MinChisq(fit, data, points, Npoints)) {
    2423    return(FALSE);
    2524  }
    … …  
    6261
    6362    // Solve with the new weights
    64     if (!weighted_LS_PM(fit, data, points, Npoints, VERBOSE)) {
    65       myAbort ("handle failures, please!");
    66       return(FALSE);
     63    if (!FitPM_MinChisq(fit, data, points, Npoints)) {
     64
     65      // restore the last solution and break
     66      fit->Ro = data->Beta_prev[0];
     67      fit->uR = data->Beta_prev[1];
     68      fit->Do = data->Beta_prev[2];
     69      fit->uD = data->Beta_prev[3];
     70     
     71      // calculate the residuals:
     72      for (i = 0; i < Npoints; i++) {
     73        points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro);
     74        points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do);
     75        points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY));
     76      }
     77      sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745;
     78      break;
    6779    }
    6880
    … …  
    90102    }
    91103  }
    92   if (!converged) {
    93     myAbort ("raise a warning on non-convergence");
    94   }
     104  fit->converged = converged;
    95105
    96106  // calculate the weight thresholds to mask the bad points:
    … …  
    156166  }
    157167
    158   // (optionally) add up the chi square for the fit, only counting the unmasked points
     168  // count unmasked points and (optionally) add up the chi square for the fit
    159169  double chisq = 0.0;
    160170  fit[0].Nfit = 0;
    … …  
    179189}
    180190
    181 int weighted_LS_PM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) {
     191int FitPM_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
     192
     193  myAssert (data->Nterms == 4, "invalid fit arrays");
    182194
    183195  int i;
    … …  
    186198
    187199  for (i = 0; i < Npoints; i++) {
     200    if (points[i].mask) continue; // respect the mask if set
     201
    188202    wx = points[i].Wx;
    189203    wy = points[i].Wy;
    … …  
    211225  data->A[0][0] = Wx;
    212226  data->A[0][1] = Tx;
     227  data->A[0][2] = 0.0;
     228  data->A[0][3] = 0.0;
    213229
    214230  data->A[1][0] = Tx;
    215231  data->A[1][1] = Tx2;
     232  data->A[1][2] = 0.0;
     233  data->A[1][3] = 0.0;
     234
     235  data->A[2][0] = 0.0;
     236  data->A[2][1] = 0.0;
    216237  data->A[2][2] = Wy;
    217238  data->A[2][3] = Ty;
     239
     240  data->A[3][0] = 0.0;
     241  data->A[3][1] = 0.0;
    218242  data->A[3][2] = Ty;
    219243  data->A[3][3] = Ty2;
    … …  
    226250
    227251  if (!dgaussjordan (data->A, data->B, 4, 1)) {
    228 # if (DEBUG)
    229     if (VERBOSE) fprintf (stderr, "error in fit\n");
    230     int j;
    231     for (i = 0; i < 4; i++) {
    232       for (j = 0; j < 4; j++) {
    233         fprintf (stderr, "%e ", data->A[i][j]);
    234       }
    235       fprintf (stderr, " : %e\n", data->B[i][0]);
    236     }
    237 # endif
    238252    return FALSE;
    239253  }
    … …  
    245259  fit->p  = 0.0;
    246260
     261  fit->dRo = sqrt(data->A[0][0]);
     262  fit->duR = sqrt(data->A[1][1]);
     263  fit->dDo = sqrt(data->A[2][2]);
     264  fit->duD = sqrt(data->A[3][3]);
     265  fit->dp  = 0.0;
     266
    247267  return TRUE;
    248268}
     269

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

    r39246 r39364  
    11# include "relastro.h"
     2int FitPMandPar_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
     3int FitPMandPar_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
     4
     5// These should probably be tunable:
     6# define MAX_ITERATIONS 10
     7# define FIT_TOLERANCE 1e-4
     8# define FLT_TOLERANCE 1e-6
     9# define WEIGHT_THRESHOLD 0.3
    210
    311// initial values of *fit are ignored
    4 int FitPMandPar (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
     12int FitPMandPar_Basic (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
     13
     14  int i;
     15
     16  // Convert the measurement errors into initial weights.
     17  for (i = 0; i < Npoints; i++) {
     18    points[i].Wx = 1.0 / SQ(points[i].dX);
     19    points[i].Wy = 1.0 / SQ(points[i].dY);
     20  }
     21
     22  int status = FitPMandPar_MinChisq (fit, data, points, Npoints);
     23  if (!status) return FALSE;
     24
     25  FitPMandPar_SetChisq (fit, data, points, Npoints);
     26  return TRUE;
     27}
     28
     29int FitPMandPar_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
     30
     31  int i,j;
     32
     33  int Ndof = 2 * Npoints - data->Nterms;
     34 
     35  // Convert the measurement errors into initial weights.
     36  for (i = 0; i < Npoints; i++) {
     37    points[i].Wx = 1.0 / SQ(points[i].dX);
     38    points[i].Wy = 1.0 / SQ(points[i].dY);
     39  }
     40 
     41  // Solve OLS equation: failure here means the chisq matrix is degenerate, give up entirely
     42  if (!FitPMandPar_MinChisq(fit, data, points, Npoints)) {
     43    return(FALSE);
     44  }
     45
     46  // Calculate r vector of residuals and least squares sigma
     47  double sigma_ols = 0.0;
     48  for (i = 0; i < Npoints; i++) {
     49    points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro + points[i].pR * fit->p);
     50    points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do + points[i].pD * fit->p);
     51    sigma_ols += SQ(points[i].rx) + SQ(points[i].ry);
     52  }
     53  sigma_ols = sqrt(sigma_ols / (float)Ndof);
     54
     55  // Save OLS covariance and Beta (solution vector, which is actually also saved in fit)
     56  for (i = 0; i < data->Nterms; i++) {
     57    for (j = 0; j < data->Nterms; j++) {
     58      data->Cov[i][j] = data->A[i][j];
     59    }
     60    data->Beta[i] = data->B[i][0];
     61  }
     62
     63  // Iteratively reweight and solve
     64  double sigma_hat = 0.0; // save for the error model
     65  int converged = FALSE;
     66  int iterations = 0;
     67
     68  // modify the weight based on the distance from the previous fit.  try up to MAX_ITERATIONS.
     69  // at the end "fit", has the last fit parameters
     70  for (iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations ++) {
     71    // Save Beta.
     72    for (i = 0; i < data->Nterms; i ++) {
     73      data->Beta_prev[i] = data->Beta[i];
     74    }
     75
     76    // Assign weights based on the deviation
     77    for (i = 0; i < Npoints; i++) {
     78      points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
     79      points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
     80    }   
     81
     82    // Solve with the new weights
     83    if (!FitPMandPar_MinChisq(fit, data, points, Npoints)) {
     84
     85      // restore the last solution and break
     86      fit->Ro = data->Beta_prev[0];
     87      fit->uR = data->Beta_prev[1];
     88      fit->Do = data->Beta_prev[2];
     89      fit->uD = data->Beta_prev[3];
     90      fit->p  = data->Beta_prev[4];
     91     
     92      // calculate the residuals:
     93      for (i = 0; i < Npoints; i++) {
     94        points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro + points[i].pR * fit->p);
     95        points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do + points[i].pD * fit->p);
     96        points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY));
     97      }
     98      sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745;
     99      break;
     100    }
     101
     102    // store the new Beta.
     103    for (i = 0; i < data->Nterms; i++) {
     104      data->Beta[i] = data->B[i][0];
     105    }
     106
     107    // calculate the residuals:
     108    for (i = 0; i < Npoints; i++) {
     109      points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro + points[i].pR * fit->p);
     110      points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do + points[i].pD * fit->p);
     111      points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY));
     112    }
     113    sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745;
     114
     115    // Check convergence
     116    converged = TRUE;
     117    for (i = 0; i < data->Nterms; i++) {
     118      // if we are within FIT_TOLERANCE as a fractional error or FLT_TOLERANCE as an absolute error, we are good
     119      if ((fabs(data->Beta[i] - data->Beta_prev[i]) > FIT_TOLERANCE * fabs(data->Beta[i])) &&
     120          (fabs(data->Beta[i] - data->Beta_prev[i]) > FLT_TOLERANCE)) {
     121        converged = FALSE;
     122      }
     123    }
     124  }
     125  fit->converged = converged;
     126
     127  // calculate the weight thresholds to mask the bad points:
     128  double Sum_Wx = 0.0;
     129  double Sum_Wy = 0.0;
     130  for (i = 0; i < Npoints; i++) {
     131    points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
     132    points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
     133   
     134    Sum_Wx += points[i].Wx;
     135    Sum_Wy += points[i].Wy;
     136  }
     137  double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints);
     138  double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints);
     139
     140  // set a mask (which can be used by the bootstrap resampling analysis)
     141  for (i = 0; i < Npoints; i++) {
     142    // keep if either is above threshold?
     143    // drop if either is below threshold?
     144    // points are marked as keep by default
     145    if ((points[i].Wx < WxThreshold) || (points[i].Wy < WyThreshold)) {
     146      points[i].mask = 1; // keep point if mask == 0
     147    }
     148  }
     149
     150  // this section calculates the formal error on the weighted fit using the covariance values
     151  // NOTE EAM: in tests (fitpm.c), they seem to be too large by a factor of ~5.37
     152  if (data->getError) {
     153    double ax = 0.0, ay = 0.0;
     154    double bx = 0.0, by = 0.0;
     155
     156    for (i = 0; i < Npoints; i++) {
     157      ax += dpsi_cauchy(points[i].rx / points[i].dX);
     158      ay += dpsi_cauchy(points[i].ry / points[i].dY);
     159
     160      bx += SQ(points[i].Wx);
     161      by += SQ(points[i].Wy);
     162    }
     163    ax /= 1.0 * Npoints;  // mean(psi_dot(r))
     164    ay /= 1.0 * Npoints;
     165    bx /= 1.0 * (Npoints - data->Nterms); // mean(psi^2(r)) * (N / (N-p))
     166    by /= 1.0 * (Npoints - data->Nterms);
     167 
     168    double lambda_x = 1.0 + (data->Nterms / Npoints) * (1 - ax) / ax;
     169    double lambda_y = 1.0 + (data->Nterms / Npoints) * (1 - ay) / ay;
     170 
     171    double sigma_robust_x = lambda_x * sqrt(bx) * sigma_hat * 2.385 / ax;
     172    double sigma_robust_y = lambda_y * sqrt(by) * sigma_hat * 2.385 / ay;
     173
     174    // This is actually sigma^2, as that's the factor in the covariance (dumouchel 4.1)
     175    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)));
     176    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)));
     177
     178    fit[0].dRo = sqrt(data->Cov[0][0]);
     179    fit[0].duR = sqrt(data->Cov[1][1]);
     180    fit[0].dDo = sqrt(data->Cov[2][2]);
     181    fit[0].duD = sqrt(data->Cov[3][3]);
     182    fit[0].dp  = sqrt(data->Cov[4][4]);
     183
     184    fit[0].dRo *= sigma_final_x;
     185    fit[0].duR *= sigma_final_x;
     186    fit[0].dDo *= sigma_final_y;
     187    fit[0].duD *= sigma_final_y;
     188    fit[0].dp  *= sqrt(sigma_final_x * sigma_final_y);
     189  }
     190
     191  FitPMandPar_SetChisq (fit, data, points, Npoints);
     192  return (TRUE);
     193}
     194
     195int FitPMandPar_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
     196
     197  myAssert (data->Nterms == 5, "invalid fit arrays");
    5198
    6199  int i;
    … …  
    9202  double PR, PD, PRT, PDT, PRX, PDY, PR2, PD2;
    10203
    11   myAssert (data->Nterms == 5, "invalid fit arrays");
    12 
    13204  PR = PD = PRT = PDT = PRX = PDY = PR2 = PD2 = 0.0;
    14205  Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
    … …  
    16207  for (i = 0; i < Npoints; i++) {
    17208    if (points[i].mask) continue; // respect the mask if set
    18     fit->Nfit ++;
    19 
    20     /* handle case where dX or dY = 0.0 */
    21     wx = 1.0 / SQ(points[i].dX);
    22     wy = 1.0 / SQ(points[i].dY);
     209
     210    wx = points[i].Wx;
     211    wy = points[i].Wy;
    23212
    24213    Wx += wx;
    25214    Wy += wy;
    26215
    27     Tx += points[i].T*wx;
    28     Ty += points[i].T*wy;
    29    
    30     Tx2 += SQ(points[i].T)*wx;
    31     Ty2 += SQ(points[i].T)*wy;
     216    double TWx = points[i].T*wx;
     217    double TWy = points[i].T*wy;
     218
     219    double XWx = points[i].X*wx;
     220    double YWy = points[i].Y*wy;
     221
     222    Tx += TWx;
     223    Ty += TWy;
     224   
     225    Tx2 += points[i].T*TWx;
     226    Ty2 += points[i].T*TWy;
    32227   
    33228    PR += points[i].pR*wx;
    34229    PD += points[i].pD*wy;
    35230   
    36     PRT += points[i].pR*points[i].T*wx;
    37     PDT += points[i].pD*points[i].T*wy;
    38    
    39     PRX += points[i].pR*points[i].X*wx;
    40     PDY += points[i].pD*points[i].Y*wy;
     231    PRT += points[i].pR*TWx;
     232    PDT += points[i].pD*TWy;
     233   
     234    PRX += points[i].pR*XWx;
     235    PDY += points[i].pD*YWy;
    41236   
    42237    PR2 += SQ(points[i].pR)*wx;
    43238    PD2 += SQ(points[i].pD)*wy;
    44239
    45     Xs += points[i].X*wx;
    46     Ys += points[i].Y*wy;
    47 
    48     XT += points[i].X*points[i].T*wx;
    49     YT += points[i].Y*points[i].T*wy;
     240    Xs += XWx;
     241    Ys += YWy;
     242
     243    XT += points[i].X*TWx;
     244    YT += points[i].Y*TWy;
    50245  }
    51246
    52247  data->A[0][0] = Wx;
    53248  data->A[0][1] = Tx;
     249  data->A[0][2] = 0.0;
     250  data->A[0][3] = 0.0;
    54251  data->A[0][4] = PR;
    55252
    56253  data->A[1][0] = Tx;
    57254  data->A[1][1] = Tx2;
     255  data->A[1][2] = 0.0;
     256  data->A[1][3] = 0.0;
    58257  data->A[1][4] = PRT;
    59258
     259  data->A[2][0] = 0.0;
     260  data->A[2][1] = 0.0;
    60261  data->A[2][2] = Wy;
    61262  data->A[2][3] = Ty;
    62263  data->A[2][4] = PD;
    63264
     265  data->A[3][0] = 0.0;
     266  data->A[3][1] = 0.0;
    64267  data->A[3][2] = Ty;
    65268  data->A[3][3] = Ty2;
    … …  
    78281  data->B[4][0] = PRX + PDY;
    79282
    80   dgaussjordan (data->A, data->B, 5, 1);
     283  if (!dgaussjordan (data->A, data->B, 5, 1)) {
     284    return FALSE;
     285  }
    81286
    82287  fit->Ro = data->B[0][0];
    … …  
    92297  fit->dp  = sqrt(data->A[4][4]);
    93298 
    94   return (TRUE);
     299  return TRUE;
    95300}
     301
     302int FitPMandPar_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
     303
     304  int i;
     305
     306  // count unmasked points and (optionally) add up the chi square for the fit
     307  double chisq = 0.0;
     308  fit[0].Nfit = 0;
     309  for (i = 0; i < Npoints; i++) {
     310    if (points[i].mask) continue;
     311    fit[0].Nfit ++;
     312     
     313    if (data->getChisq) {
     314      double Xf = fit[0].Ro + fit[0].uR*points[i].T + fit[0].p*points[i].pR;
     315      double Yf = fit[0].Do + fit[0].uD*points[i].T + fit[0].p*points[i].pD;
     316      double wx = 1.0 / SQ(points[i].dX);
     317      double wy = 1.0 / SQ(points[i].dY);
     318      chisq += SQ(points[i].X - Xf) * wx;
     319      chisq += SQ(points[i].Y - Yf) * wy;
     320    }
     321  }
     322   
     323  // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms)
     324  fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms);
     325
     326  return TRUE;
     327}

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

    r39246 r39364  
    1919  }
    2020 
    21   // Solve OLS equation
     21  // Solve OLS equation: failure here means the chisq matrix is degenerate, give up entirely
    2222  if (!weighted_LS_PLX(fit, data, points, Npoints, VERBOSE)) {
    23     myAbort ("handle failures, please!");
    2423    return(FALSE);
    2524  }
    … …  
    6362    // Solve with the new weights
    6463    if (!weighted_LS_PLX(fit, data, points, Npoints, VERBOSE)) {
    65       myAbort ("handle failures, please!");
    66       return(FALSE);
     64
     65      // restore the last solution and break
     66      fit->Ro = data->Beta_prev[0];
     67      fit->uR = data->Beta_prev[1];
     68      fit->Do = data->Beta_prev[2];
     69      fit->uD = data->Beta_prev[3];
     70      fit->p  = data->Beta_prev[4];
     71     
     72      // calculate the residuals:
     73      for (i = 0; i < Npoints; i++) {
     74        points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro + points[i].pR * fit->p);
     75        points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do + points[i].pD * fit->p);
     76        points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY));
     77      }
     78      sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745;
     79      break;
    6780    }
    6881
    … …  
    90103    }
    91104  }
    92   if (!converged) {
    93     fprintf (stderr, "raise a warning on non-convergence\n");
    94   }
     105  fit->converged = converged;
    95106
    96107  // calculate the weight thresholds to mask the bad points:
    … …  
    158169  }
    159170
    160   // (optionally) add up the chi square for the fit, only counting the unmasked points
     171  // count unmasked points and (optionally) add up the chi square for the fit
    161172  double chisq = 0.0;
    162173  fit[0].Nfit = 0;
    … …  
    183194int weighted_LS_PLX (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) {
    184195
     196  myAssert (data->Nterms == 5, "invalid fit arrays");
     197
    185198  int i;
    186199  double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
    … …  
    227240  data->A[0][0] = Wx;
    228241  data->A[0][1] = Tx;
     242  data->A[0][2] = 0.0;
     243  data->A[0][3] = 0.0;
    229244  data->A[0][4] = PR;
    230245
    231246  data->A[1][0] = Tx;
    232247  data->A[1][1] = Tx2;
     248  data->A[1][2] = 0.0;
     249  data->A[1][3] = 0.0;
    233250  data->A[1][4] = PRT;
    234251
     252  data->A[2][0] = 0.0;
     253  data->A[2][1] = 0.0;
    235254  data->A[2][2] = Wy;
    236255  data->A[2][3] = Ty;
    237256  data->A[2][4] = PD;
    238257
     258  data->A[3][0] = 0.0;
     259  data->A[3][1] = 0.0;
    239260  data->A[3][2] = Ty;
    240261  data->A[3][3] = Ty2;
    … …  
    254275
    255276  if (!dgaussjordan (data->A, data->B, 5, 1)) {
    256 # if (DEBUG)
    257     if (VERBOSE) fprintf (stderr, "error in fit\n");
    258     int j;
    259     for (i = 0; i < 5; i++) {
    260       for (j = 0; j < 5; j++) {
    261         fprintf (stderr, "%e ", data->A[i][j]);
    262       }
    263       fprintf (stderr, " : %e\n", data->B[i][0]);
    264     }
    265 # endif
    266277    return FALSE;
    267278  }

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

    r38986 r39364  
    11# include "relastro.h"
    2 
    3 int FitPosPMfixed (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
     2int FitPosPMfixed_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
     3int FitPosPMfixed_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
     4
     5// These should probably be tunable:
     6# define MAX_ITERATIONS 10
     7# define FIT_TOLERANCE 1e-4
     8# define FLT_TOLERANCE 1e-6
     9# define WEIGHT_THRESHOLD 0.3
     10
     11// fit->uR,uD are used in the fit
     12int FitPosPMfixed_Basic (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
    413
    514  int i;
    615
     16  // Convert the measurement errors into initial weights.
     17  for (i = 0; i < Npoints; i++) {
     18    points[i].Wx = 1.0 / SQ(points[i].dX);
     19    points[i].Wy = 1.0 / SQ(points[i].dY);
     20  }
     21
     22  int status = FitPosPMfixed_MinChisq (fit, data, points, Npoints);
     23  if (!status) return FALSE;
     24
     25  FitPosPMfixed_SetChisq (fit, data, points, Npoints);
     26  return TRUE;
     27}
     28
     29int FitPosPMfixed_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
     30
     31  int i,j;
     32
     33  int Ndof = 2 * Npoints - data->Nterms;
     34 
     35  // Convert the measurement errors into initial weights.
     36  for (i = 0; i < Npoints; i++) {
     37    points[i].Wx = 1.0 / SQ(points[i].dX);
     38    points[i].Wy = 1.0 / SQ(points[i].dY);
     39  }
     40 
     41  // Solve OLS equation 
     42  if (!FitPosPMfixed_MinChisq(fit, data, points, Npoints)) {
     43    return(FALSE);
     44  }
     45
     46  // Calculate r vector of residuals and least squares sigma
     47  double sigma_ols = 0.0;
     48  for (i = 0; i < Npoints; i++) {
     49    points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro);
     50    points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do);
     51    sigma_ols += SQ(points[i].rx) + SQ(points[i].ry);
     52  }
     53  sigma_ols = sqrt(sigma_ols / (float)Ndof);
     54
     55  // Save OLS covariance and Beta (solution vector, which is actually also saved in fit)
     56  for (i = 0; i < data->Nterms; i++) {
     57    for (j = 0; j < data->Nterms; j++) {
     58      data->Cov[i][j] = data->A[i][j];
     59    }
     60    data->Beta[i] = data->B[i][0];
     61  }
     62
     63  // Iteratively reweight and solve
     64  double sigma_hat = 0.0; // save for the error model
     65  int converged = FALSE;
     66  int iterations = 0;
     67
     68  // modify the weight based on the distance from the previous fit.  try up to MAX_ITERATIONS.
     69  // at the end "fit", has the last fit parameters
     70  for (iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations ++) {
     71    // Save Beta.
     72    for (i = 0; i < data->Nterms; i ++) {
     73      data->Beta_prev[i] = data->Beta[i];
     74    }
     75
     76    // Assign weights based on the deviation
     77    for (i = 0; i < Npoints; i++) {
     78      points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
     79      points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
     80    }   
     81
     82    // Solve with the new weights
     83    if (!FitPosPMfixed_MinChisq(fit, data, points, Npoints)) {
     84
     85      // restore the last solution and break
     86      fit->Ro = data->Beta_prev[0];
     87      fit->Do = data->Beta_prev[1];
     88     
     89      // calculate the residuals:
     90      for (i = 0; i < Npoints; i++) {
     91        points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro);
     92        points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do);
     93        points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY));
     94      }
     95      sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745;
     96      break;
     97    }
     98
     99    // store the new Beta.
     100    for (i = 0; i < data->Nterms; i++) {
     101      data->Beta[i] = data->B[i][0];
     102    }
     103
     104    // calculate the residuals:
     105    for (i = 0; i < Npoints; i++) {
     106      points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro);
     107      points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do);
     108      points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY));
     109    }
     110    sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745;
     111   
     112    // Check convergence
     113    converged = TRUE;
     114    for (i = 0; i < data->Nterms; i++) {
     115      // if we are within FIT_TOLERANCE as a fractional error or FLT_TOLERANCE as an absolute error, we are good
     116      if ((fabs(data->Beta[i] - data->Beta_prev[i]) > FIT_TOLERANCE * fabs(data->Beta[i])) &&
     117          (fabs(data->Beta[i] - data->Beta_prev[i]) > FLT_TOLERANCE)) {
     118        converged = FALSE;
     119      }
     120    }
     121  }
     122  fit->converged = converged;
     123
     124  // calculate the weight thresholds to mask the bad points:
     125  double Sum_Wx = 0.0;
     126  double Sum_Wy = 0.0;
     127  for (i = 0; i < Npoints; i++) {
     128    points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
     129    points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
     130   
     131    Sum_Wx += points[i].Wx;
     132    Sum_Wy += points[i].Wy;
     133  }
     134  double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints);
     135  double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints);
     136
     137  // set a mask (which can be used by the bootstrap resampling analysis)
     138  for (i = 0; i < Npoints; i++) {
     139    // keep if either is above threshold?
     140    // drop if either is below threshold?
     141    // points are marked as keep by default
     142    if ((points[i].Wx < WxThreshold) || (points[i].Wy < WyThreshold)) {
     143      points[i].mask = 1; // keep point if mask == 0
     144    }
     145  }
     146
     147  // this section calculates the formal error on the weighted fit using the covariance values
     148  // NOTE EAM: in tests (fitpm.c), they seem to be too large by a factor of ~5.37
     149  if (data->getError) {
     150    double ax = 0.0, ay = 0.0;
     151    double bx = 0.0, by = 0.0;
     152
     153    for (i = 0; i < Npoints; i++) {
     154      ax += dpsi_cauchy(points[i].rx / points[i].dX);
     155      ay += dpsi_cauchy(points[i].ry / points[i].dY);
     156
     157      bx += SQ(points[i].Wx);
     158      by += SQ(points[i].Wy);
     159    }
     160    ax /= 1.0 * Npoints;  // mean(psi_dot(r))
     161    ay /= 1.0 * Npoints;
     162    bx /= 1.0 * (Npoints - data->Nterms); // mean(psi^2(r)) * (N / (N-p))
     163    by /= 1.0 * (Npoints - data->Nterms);
     164 
     165    double lambda_x = 1.0 + (data->Nterms / Npoints) * (1 - ax) / ax;
     166    double lambda_y = 1.0 + (data->Nterms / Npoints) * (1 - ay) / ay;
     167 
     168    double sigma_robust_x = lambda_x * sqrt(bx) * sigma_hat * 2.385 / ax;
     169    double sigma_robust_y = lambda_y * sqrt(by) * sigma_hat * 2.385 / ay;
     170
     171    // This is actually sigma^2, as that's the factor in the covariance (dumouchel 4.1)
     172    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)));
     173    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)));
     174
     175    fit[0].dRo = sqrt(data->Cov[0][0]);
     176    fit[0].dDo = sqrt(data->Cov[1][1]);
     177
     178    fit[0].dRo *= sigma_final_x;
     179    fit[0].dDo *= sigma_final_y;
     180  }
     181
     182  FitPosPMfixed_SetChisq (fit, data, points, Npoints);
     183  return (TRUE);
     184}
     185
     186int FitPosPMfixed_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
     187
     188  myAssert (data->Nterms == 2, "invalid fit arrays");
     189
     190  int i;
    7191  double wx, wy, Wx, Wy, Tx, Ty, Xs, Ys;
    8 
    9   myAssert (data->Nterms == 2, "invalid fit arrays");
    10 
    11192  Wx = Wy = Tx = Ty = Xs = Ys = 0.0;
    12193
    13194  for (i = 0; i < Npoints; i++) {
    14     /* handle case where dX or dY = 0.0 */
    15     wx = 1.0 / SQ(points[i].dX);
    16     wy = 1.0 / SQ(points[i].dY);
     195    if (points[i].mask) continue; // respect the mask if set
     196
     197    wx = points[i].Wx;
     198    wy = points[i].Wy;
    17199
    18200    Wx += wx;
    … …  
    26208  }
    27209
    28   fit->Ro = (Xs - fit->uR*Tx) / Wx;
    29   fit->Do = (Ys - fit->uD*Ty) / Wy;
    30 
     210  // X^T W X
     211  data->A[0][0] = Wx;
     212  data->A[0][1] = 0.0;
     213
     214  data->A[1][0] = 0.0;
     215  data->A[1][1] = Wy;
     216
     217  // X^T W Y
     218  data->B[0][0] = Xs - fit->uR*Tx;
     219  data->B[1][0] = Ys - fit->uD*Ty;
     220
     221  if (!dgaussjordan (data->A, data->B, 2, 1)) {
     222    return FALSE;
     223  }
     224
     225  fit->Ro = data->B[0][0];
     226  fit->Do = data->B[1][0];
    31227  fit->p  = 0.0;
    32  
    33   // should be ~ 1.0 / Wx
    34   fit->dRo = sqrt(1.0 / Wx);
    35   fit->dDo = sqrt(1.0 / Wy);
    36  
    37   fit->Nfit = Npoints;
    38 
    39   return (TRUE);
    40 }
     228
     229  fit->dRo = sqrt(data->A[0][0]);
     230  fit->dDo = sqrt(data->A[1][1]);
     231  fit->duR = NAN;
     232  fit->duD = NAN;
     233  fit->dp  = 0.0;
     234
     235  return TRUE;
     236}
     237
     238int FitPosPMfixed_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints) {
     239
     240  int i;
     241
     242  // count unmasked points and (optionally) add up the chi square for the fit
     243  double chisq = 0.0;
     244  fit[0].Nfit = 0;
     245  for (i = 0; i < Npoints; i++) {
     246    if (points[i].mask) continue;
     247    fit[0].Nfit ++;
     248     
     249    if (data->getChisq) {
     250      double Xf = fit[0].Ro + fit[0].uR*points[i].T;
     251      double Yf = fit[0].Do + fit[0].uD*points[i].T;
     252      double wx = 1.0 / SQ(points[i].dX);
     253      double wy = 1.0 / SQ(points[i].dY);
     254      chisq += SQ(points[i].X - Xf) * wx;
     255      chisq += SQ(points[i].Y - Yf) * wy;
     256    }
     257  }
     258   
     259  // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms)
     260  fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms);
     261
     262  return TRUE;
     263}

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

    r39287 r39364  
    2020 
    2121  // Solve OLS equation 
    22   if (!weighted_LS_PM(fit, data, points, Npoints, VERBOSE)) {
    23     myAbort ("handle failures, please!");
     22  if (!weighted_LS_POS(fit, data, points, Npoints, VERBOSE)) {
    2423    return(FALSE);
    2524  }
    … …  
    6261
    6362    // Solve with the new weights
    64     if (!weighted_LS_PM(fit, data, points, Npoints, VERBOSE)) {
    65       myAbort ("handle failures, please!");
    66       return(FALSE);
     63    if (!weighted_LS_POS(fit, data, points, Npoints, VERBOSE)) {
     64
     65      // restore the last solution and break
     66      fit->Ro = data->Beta_prev[0];
     67      fit->Do = data->Beta_prev[1];
     68     
     69      // calculate the residuals:
     70      for (i = 0; i < Npoints; i++) {
     71        points[i].rx = points[i].X - (points[i].T * fit->uR + fit->Ro);
     72        points[i].ry = points[i].Y - (points[i].T * fit->uD + fit->Do);
     73        points[i].u = sqrt(SQ(points[i].rx / points[i].dX) + SQ(points[i].ry / points[i].dY));
     74      }
     75      sigma_hat = MedianAbsDeviation(points, Npoints) / 0.6745;
     76      break;
    6777    }
    6878
    … …  
    90100    }
    91101  }
    92   if (!converged) {
    93     myAbort ("raise a warning on non-convergence");
    94   }
     102  fit->converged = converged;
    95103
    96104  // calculate the weight thresholds to mask the bad points:
    … …  
    146154
    147155    fit[0].dRo = sqrt(data->Cov[0][0]);
    148     fit[0].duR = sqrt(data->Cov[1][1]);
    149     fit[0].dDo = sqrt(data->Cov[2][2]);
    150     fit[0].duD = sqrt(data->Cov[3][3]);
     156    fit[0].dDo = sqrt(data->Cov[1][1]);
    151157
    152158    fit[0].dRo *= sigma_final_x;
    153     fit[0].duR *= sigma_final_x;
    154159    fit[0].dDo *= sigma_final_y;
    155     fit[0].duD *= sigma_final_y;
    156   }
    157 
    158   // (optionally) add up the chi square for the fit, only counting the unmasked points
     160  }
     161
     162  // count unmasked points and (optionally) add up the chi square for the fit
    159163  double chisq = 0.0;
    160164  fit[0].Nfit = 0;
    … …  
    179183}
    180184
    181 int weighted_LS_PM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) {
     185int weighted_LS_POS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) {
     186
     187  myAssert (data->Nterms == 2, "invalid fit arrays");
    182188
    183189  int i;
    184   double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
    185   Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
     190  double wx, wy, Wx, Wy, Tx, Ty, Xs, Ys;
     191
     192  Wx = Wy = Tx = Ty = Xs = Ys = 0.0;
    186193
    187194  for (i = 0; i < Npoints; i++) {
    … …  
    192199    Wy += wy;
    193200
    194     double TWx = points[i].T*wx;
    195     double TWy = points[i].T*wy;
    196 
    197     Tx += TWx;
    198     Ty += TWy;
    199    
    200     Tx2 += points[i].T*TWx;
    201     Ty2 += points[i].T*TWy;
     201    Tx += points[i].T*wx;
     202    Ty += points[i].T*wy;
    202203   
    203204    Xs += points[i].X*wx;
    204205    Ys += points[i].Y*wy;
    205 
    206     XT += points[i].X*TWx;
    207     YT += points[i].Y*TWy;
    208206  }
    209207
    210208  // X^T W X
    211209  data->A[0][0] = Wx;
    212   data->A[0][1] = Tx;
    213 
    214   data->A[1][0] = Tx;
    215   data->A[1][1] = Tx2;
    216   data->A[2][2] = Wy;
    217   data->A[2][3] = Ty;
    218   data->A[3][2] = Ty;
    219   data->A[3][3] = Ty2;
     210  data->A[0][1] = 0.0;
     211
     212  data->A[1][0] = 0.0;
     213  data->A[1][1] = Wy;
    220214
    221215  // X^T W Y
    222   data->B[0][0] = Xs;
    223   data->B[1][0] = XT;
    224   data->B[2][0] = Ys;
    225   data->B[3][0] = YT;
    226 
    227   if (!dgaussjordan (data->A, data->B, 4, 1)) {
    228 # if (DEBUG)
    229     if (VERBOSE) fprintf (stderr, "error in fit\n");
    230     int j;
    231     for (i = 0; i < 4; i++) {
    232       for (j = 0; j < 4; j++) {
    233         fprintf (stderr, "%e ", data->A[i][j]);
    234       }
    235       fprintf (stderr, " : %e\n", data->B[i][0]);
    236     }
    237 # endif
     216  data->B[0][0] = Xs - fit->uR*Tx;
     217  data->B[1][0] = Ys - fit->uD*Ty;
     218
     219  if (!dgaussjordan (data->A, data->B, 2, 1)) {
    238220    return FALSE;
    239221  }
    240222
    241223  fit->Ro = data->B[0][0];
    242   fit->uR = data->B[1][0];
    243   fit->Do = data->B[2][0];
    244   fit->uD = data->B[3][0];
     224  fit->Do = data->B[1][0];
    245225  fit->p  = 0.0;
    246226

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

    r39287 r39364  
    11# include "relastro.h"
    22# define PAR_TOOFEW 5
     3
     4typedef enum {
     5  SELECT_MEAS_HAS_DATA  = 1,
     6  SELECT_MEAS_HAS_STACK = 2,
     7  SELECT_MEAS_HAS_2MASS = 4,
     8} SelectMeasureStatus;
    39
    410int DumpObjectsWith2MASS (Catalog *catalog, int Ncatalog);
    … …  
    112118
    113119  int mode = FIT_MODE; // start with the globally-defined fit mode
     120  // mode may be one of FIT_AVERAGE (position only), FIT_PM_ONLY, FIT_PM_AND_PAR
    114121
    115122  int XVERB = FALSE;
    … …  
    135142  FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange);
    136143
    137   // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
     144  // To judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq.
     145  // we start with the parallax fits, and step down to fewer and fewer parameters
    138146
    139147  // if we have too few good detections for the desired fit, or too limited a baseline,
    … …  
    141149
    142150  // 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
    145151  if (mode == FIT_PM_AND_PAR) {
    146152    if (Trange < PM_DT_MIN) {
     153      // not enough baseline for proper motion, only set mean position
    147154      mode = FIT_AVERAGE;
    148155      goto justPosition;
    149156    }
    150157    if (parRange < PAR_FACTOR_MIN) {
     158      // not enough parallax factor range, skip parallax
    151159      mode = FIT_PM_ONLY;
    152160      goto skipParallax;
    153161    }
    154162    if (fitStats->Npoints <= PAR_TOOFEW) {
     163      // not enough data, skip parallax
    155164      mode = FIT_PM_ONLY;
    156165      goto skipParallax;
    … …  
    159168    // we are going to use the IRLS analysis to calculate the mean solution and the masking
    160169    // 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) {
     170    fitStats->fitdataPar->getError = !N_BOOTSTRAP_SAMPLES;
     171    if (!FitPMandPar_IRLS (&fitPar, fitStats->fitdataPar, fitStats->points, fitStats->Npoints)) {
     172      mode = FIT_PM_ONLY;
     173      goto skipParallax;
     174    }
     175
     176    if (N_BOOTSTRAP_SAMPLES) {
    165177      fitStats->Nfit = 0;
    166178      int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
    167179      if (Nnomask < 5) {
    168         myAbort ("handle this case, please");
     180        // if we do not have enough points to fit parallax, we cannot do the bootstrap analysis.
     181        mode = FIT_PM_ONLY;
     182        goto skipParallax;
    169183      }
    170184      for (k = 0; k < fitStats->NfitAlloc; k++) {
    171185        BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
    172         FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask);
     186        if (!FitPMandPar_Basic (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask)) continue;
    173187        fitStats->Nfit ++;
    174188      }
    … …  
    190204    // XXX a hard-wired hack...
    191205    // 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)) {
     206    int valid = TRUE;
     207    valid = valid && isfinite(fitPar.uR);
     208    valid = valid && isfinite(fitPar.uD);
     209    valid = valid && isfinite(fitPar.p);
     210    valid = valid && isfinite(fitPar.duR);
     211    valid = valid && isfinite(fitPar.duD);
     212    valid = valid && isfinite(fitPar.dp);
     213    valid = valid && (fabs(fitPar.uR) < 4.0);
     214    valid = valid && (fabs(fitPar.uD) < 4.0);
     215    valid = valid && (fabs(fitPar.p) < 2.0);
     216    if (!valid) {
    193217      mode = FIT_PM_ONLY;
    194218    } else {
    … …  
    212236    if (average[0].flags & ID_STAR_USE_PAR) {
    213237      // get the best pm fit and chisq given the set of points (mask is respected)
    214       if (!FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
    215         myAbort ("oops");
     238      if (!FitPM_Basic (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
     239        average[0].flags |= ID_STAR_BAD_PM;
     240        goto justPosition;
    216241      }
    217242    } else {
    218       fitStats->fitdataPM->getError = (N_BOOTSTRAP_SAMPLES < 3);
    219       FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints, VERBOSE);
    220 
    221       if (N_BOOTSTRAP_SAMPLES >= 3) {
     243      fitStats->fitdataPM->getError = !N_BOOTSTRAP_SAMPLES;
     244      if (!FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) {
     245        mode = FIT_AVERAGE;
     246        goto justPosition;
     247      }
     248      if (N_BOOTSTRAP_SAMPLES) {
    222249        fitStats->Nfit = 0;
    223250        int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
    224251        if (Nnomask < 5) {
    225           myAbort ("handle this case, please");
     252          mode = FIT_AVERAGE;
     253          goto justPosition;
    226254        }
    227255        for (k = 0; k < fitStats->NfitAlloc; k++) {
    228256          BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
    229           if (!FitPM (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue;
     257          if (!FitPM_Basic (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, fitStats->Npoints)) continue;
    230258          fitStats->Nfit ++;
    231259        }
    … …  
    244272    fitStats->Npm ++;
    245273
    246     // XXX a hard-wired hack...
    247274    // unless there is a clear problems (below) with the proper-motion fit or we have a parallax fit, we will use pm fit
    248     if ((fabs(fitPM.uR) > 4.0) || (fabs(fitPM.uD) > 4.0)) {
     275    int valid = TRUE;
     276    valid = valid && isfinite(fitPM.uR);
     277    valid = valid && isfinite(fitPM.uD);
     278    valid = valid && isfinite(fitPM.duR);
     279    valid = valid && isfinite(fitPM.duD);
     280    valid = valid && (fabs(fitPM.uR) < 4.0);
     281    valid = valid && (fabs(fitPM.uD) < 4.0);
     282    if (!valid) {
    249283      mode = FIT_AVERAGE;
    250284      average[0].flags |= ID_STAR_BAD_PM;
    … …  
    260294    // use bootstrap resampling to check the error distribution
    261295    // if we only have one point, this is silly...
    262    
    263     // XXX I need to rethink here : use a median for the position or weighted average? use IRLS anyway?
    264 
    265     if (fitStats->NfitAlloc == 1) {
     296
     297    FitAstromResultSetPM (&fitPos, 1, average);
     298    if ((average[0].flags & ID_STAR_USE_PAR) || (average[0].flags & ID_STAR_USE_PM)) {
     299      if (!FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
     300        goto doneWithFit;
     301      }
     302    } else {
    266303      FitAstromResultSetPM (&fitPos, 1, average);
    267       FitPosPMfixed (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints);
    268     } else {
    269       fitStats->Nfit = 0;
    270       FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average);
    271       for (k = 0; k < fitStats->NfitAlloc; k++) {
    272         BootstrapResample (fitStats->sample, fitStats->points, fitStats->Npoints);
    273         FitPosPMfixed (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, fitStats->Npoints);
    274         fitStats->Nfit ++;
    275       }
    276       BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
    277       BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
    278     }
    279     FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE);
     304      if (!FitPosPMfixed_IRLS (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) {
     305        goto doneWithFit;
     306      }
     307      if (N_BOOTSTRAP_SAMPLES) {
     308        fitStats->Nfit = 0;
     309        int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
     310        if (Nnomask < 5) {
     311          goto doneWithFit;
     312        }
     313        FitAstromResultSetPM (fitStats->fit, fitStats->NfitAlloc, average);
     314        for (k = 0; k < fitStats->NfitAlloc; k++) {
     315          BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
     316          if (!FitPosPMfixed_Basic (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, fitStats->Npoints)) continue;
     317          fitStats->Nfit ++;
     318        }
     319        BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     320        BootstrapRobustStats (&fitPos, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     321      }
     322      average[0].flags |= ID_STAR_USE_AVE;
     323    }
    280324
    281325    // project Ro, Do back to RA,DEC
    282326    XY_to_RD (&fitPos.Ro, &fitPos.Do, fitPos.Ro, fitPos.Do, &fitStats->coords);
     327    if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
     328
    283329    average[0].flags |= ID_STAR_FIT_AVE;
    284330    fitStats->Nave ++;
    285331  }
    286332
     333 doneWithFit:
    287334  // update the bit flags of which points were used
    288335  for (k = 0; k < fitStats->Npoints; k++) {
    … …  
    333380  /* choose the result based on the chisq values */
    334381  // XXXX for now, just use the mode as the result:
    335   int result = mode;
    336382  FitAstromResult fit;
    337383  FitAstromResultInit (&fit);
    338384
    339   switch (result) {
    340     case FIT_AVERAGE:
    341       average[0].flags |= ID_STAR_USE_AVE;
    342       fit = fitPos;
    343       break;
    344     case FIT_PM_ONLY:
    345       average[0].flags |= ID_STAR_USE_PM;
    346       fit = fitPM;
    347       break;
    348     case FIT_PM_AND_PAR:
    349       average[0].flags |= ID_STAR_USE_PAR;
    350       fit = fitPar;
    351       break;
    352   }
     385  if (average[0].flags & ID_STAR_USE_PAR) {
     386    myAssert ((average[0].flags & (ID_STAR_USE_PM | ID_STAR_USE_AVE)) == 0, "programming error");
     387    fit = fitPar;
     388  }
     389  if (average[0].flags & ID_STAR_USE_PM) {
     390    myAssert ((average[0].flags & (ID_STAR_USE_PAR | ID_STAR_USE_AVE)) == 0, "programming error");
     391    fit = fitPM;
     392  }
     393  if (average[0].flags & ID_STAR_USE_AVE) {
     394    myAssert ((average[0].flags & (ID_STAR_USE_PAR | ID_STAR_USE_PM)) == 0, "programming error");
     395    fit = fitPos;
     396  }
     397
    353398  if (XVERB) {
    354399    fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) plx=(%f +/- %f)\n",
    … …  
    376421  if (!status) {
    377422    fitStats->Nskip ++;
    378     return FALSE; // XXX ??
     423    return FALSE;
    379424  }
    380425
    … …  
    441486int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats) {
    442487
     488  int status;
     489
    443490  // set the default values
    444491  average[0].Rstk  = NAN; // RA in degrees
    … …  
    458505
    459506  // select the measurements to be used in this analysis
    460   UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, TRUE);
    461 
    462   // too few measurements for average position (require 2 values)
    463   if (fitStats->Npoints < 1) return FALSE; // XXX ??
     507  status = UpdateObjects_SelectMeasures (fitStats, average, secfilt, measure, measureBig, TRUE);
     508  if (status & SELECT_MEAS_HAS_STACK) {
     509    // if we have any stack measurements, set the default value to the mean position
     510    // this is set first, so the position may only be the original stack position
     511    // if we exit this function without updating the stack position, it means we
     512    // have a stack measurement, but not a good stack measuremetn
     513    average[0].Rstk  = average[0].R; // RA in degrees
     514    average[0].Dstk  = average[0].D; // DEC in degrees
     515    average[0].dRstk = average[0].dR; // RA scatter in arcsec
     516    average[0].dDstk = average[0].dD; // DEC scatter in arcsec
     517  }
     518
     519  // no stack measurements of acceptable quality (defaults to mean position)
     520  if (fitStats->Npoints < 1) return FALSE;
    464521 
    465522  double Tmean, Trange, parRange;
    466523  FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange);
    467524
    468   FitPosPMfixed (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints);
     525  // the stack positions are not statistically independent...
     526  FitPosPMfixed_Basic (&fitPos, fitStats->fitdataPos, fitStats->points, fitStats->Npoints);
    469527  FitAstromSetChisq (&fitPos, fitStats->points, fitStats->Npoints, FIT_AVERAGE);
    470528
    … …  
    482540                      3600*(average[0].D - fitPos.Do));
    483541
    484   // make sure that the fit succeeded
    485   int status = TRUE;
     542  // check if the fit succeeded
     543  status = TRUE;
    486544  status &= finite(fitPos.Ro);
    487545  status &= finite(fitPos.Do);
    … …  
    522580
    523581  int has2MASS = FALSE;
     582  int hasStack = FALSE;
    524583
    525584  int Npoints = fit->Npoints = 0;
    … …  
    550609    if (isStack) {
    551610      if (!isGPC1stack(measure[k].photcode)) continue;
     611      hasStack = TRUE;
    552612    } else {
    553613      if ( isGPC1stack(measure[k].photcode)) continue;
    … …  
    616676    if (isnan(points[Npoints].dY)) continue;
    617677
     678    // avoid dX,dY == 0.0 so we do not have to constantly test for it
     679    points[Npoints].dX = MAX(points[Npoints].dX, 0.001);
     680    points[Npoints].dY = MAX(points[Npoints].dY, 0.001);
     681
    618682    points[Npoints].dT = measure[k].dt;
    619683
    … …  
    646710
    647711  fit->Npoints = Npoints;
    648   return TRUE;
     712
     713  int status = SELECT_MEAS_HAS_DATA;
     714  if (hasStack) status |= SELECT_MEAS_HAS_STACK;
     715  if (has2MASS) status |= SELECT_MEAS_HAS_2MASS;
     716  return status;
    649717}
    650718

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

    r39225 r39364  
    121121    FIT_MODE = FIT_PM_ONLY;
    122122  }
    123   if ((N = get_argument (argc, argv, "-par"))) {
    124     remove_argument (N, &argc, argv);
    125     FIT_MODE = FIT_PAR_ONLY;
    126   }
    127123  if ((N = get_argument (argc, argv, "-pmpar"))) {
    128124    remove_argument (N, &argc, argv);
    129125    FIT_MODE = FIT_PM_AND_PAR;
    130126  }
     127
     128  // NOTE : this is no longer a valid mode
     129  // if ((N = get_argument (argc, argv, "-par"))) {
     130  //   remove_argument (N, &argc, argv);
     131  //   FIT_MODE = FIT_PAR_ONLY;
     132  // }
    131133
    132134  if ((N = get_argument (argc, argv, "-high-speed"))) {

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

    r39287 r39364  
    1212# define OUT_ERROR 0.500 /* arcsec */
    1313
    14 # define N_STARS 3000
    15 # define N_POINTS 100
    16 # define N_OUTLIERS 5
    17 # define N_BOOTSTRAP 100
     14// # define N_STARS 3000
     15// # define N_POINTS 100
     16// # define N_OUTLIERS 5
     17// # define N_BOOTSTRAP 100
    1818
    1919# define dcos(THETA) cos(RAD_DEG*THETA)
    … …  
    2121
    2222enum {
    23   FIT_PM_NONE,
     23  FIT_POS_NONE,
     24  FIT_POS_IRLS,
     25  FIT_POS_BOOT,
     26  FIT_POS_NOCLIP,
    2427  FIT_PM_IRLS,
     28  FIT_PM_BOOT,
    2529  FIT_PM_NOCLIP,
    2630  FIT_PLX_IRLS,
    … …  
    3135int main (int argc, char **argv) {
    3236 
     37  int N_STARS     = 3000;
     38  int N_POINTS    =  100;
     39  int N_OUTLIERS  =   10;
     40  int N_BOOTSTRAP =  100;
     41
     42  int N;
     43  if ((N = get_argument (argc, argv, "-Nstars"))) {
     44    remove_argument (N, &argc, argv);
     45    N_STARS = atoi (argv[N]);
     46    remove_argument (N, &argc, argv);
     47  }
     48  if ((N = get_argument (argc, argv, "-Npoints"))) {
     49    remove_argument (N, &argc, argv);
     50    N_POINTS = atoi (argv[N]);
     51    remove_argument (N, &argc, argv);
     52  }
     53  if ((N = get_argument (argc, argv, "-Noutliers"))) {
     54    remove_argument (N, &argc, argv);
     55    N_OUTLIERS = atoi (argv[N]);
     56    remove_argument (N, &argc, argv);
     57  }
     58  if ((N = get_argument (argc, argv, "-Nbootstrap"))) {
     59    remove_argument (N, &argc, argv);
     60    N_BOOTSTRAP = atoi (argv[N]);
     61    remove_argument (N, &argc, argv);
     62  }
     63
    3364  if (argc != 2) {
    3465    fprintf (stderr, "USAGE: %s (mode)\n", argv[0]);
     66    fprintf (stderr, "  mode: pos, pos-irls, pos-boot, pm, pm-irls, pm-boot, plx, plx-irls, plx-boot\n");
     67    fprintf (stderr, " options: -Nstars [3000], -Npoints [100], -Noutliers [10], -Nbootstrap [100]\n");
    3568    exit (2);
    3669  }
    3770
    38   int mode = FIT_NONE;
     71  int mode = FIT_POS_NONE;
     72  if (!strcasecmp(argv[1], "pos")) {
     73    mode = FIT_POS_NOCLIP;
     74  }
     75  if (!strcasecmp(argv[1], "pos-irls")) {
     76    mode = FIT_POS_IRLS;
     77  }
     78  if (!strcasecmp(argv[1], "pos-boot")) {
     79    mode = FIT_POS_BOOT;
     80  }
    3981  if (!strcasecmp(argv[1], "pm")) {
    4082    mode = FIT_PM_NOCLIP;
    … …  
    4284  if (!strcasecmp(argv[1], "pm-irls")) {
    4385    mode = FIT_PM_IRLS;
     86  }
     87  if (!strcasecmp(argv[1], "pm-boot")) {
     88    mode = FIT_PM_BOOT;
    4489  }
    4590  if (!strcasecmp(argv[1], "plx")) {
    … …  
    72117  fitStats->fitdataPar->getChisq = TRUE;
    73118  fitStats->fitdataPar->getError = TRUE;
     119  fitStats->fitdataPos->getChisq = TRUE;
     120  fitStats->fitdataPos->getError = TRUE;
    74121  // NOTE: surprisingly, the chisq and error calculations only add ~2-3% to the exec time
    75122
    … …  
    88135    double plx;
    89136    switch (mode) {
     137      case FIT_POS_IRLS:
     138      case FIT_POS_BOOT:
     139      case FIT_POS_NOCLIP:
    90140      case FIT_PM_IRLS:
     141      case FIT_PM_BOOT:
    91142      case FIT_PM_NOCLIP:
    92         plx = 0.0*(drand48() + 0.0);
     143        plx = 0.0;
    93144        break;
    94145      case FIT_PLX_IRLS:
    … …  
    101152    }
    102153
    103 
     154    // generate a single fake star with N_POINTS real points and N_OUTLIERS bad points
    104155    mkstar (fitStats, Ro, Do, uR, uD, plx, N_POINTS, N_OUTLIERS);
    105156
    … …  
    107158    FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange);
    108159
    109     FitAstromResult fitPM;   
    110     FitAstromResultInit (&fitPM);
     160    FitAstromResult fitResult;   
     161    FitAstromResultInit (&fitResult);
     162
     163    int Nnomask;
    111164
    112165    switch (mode) {
    113       case FIT_PM_NOCLIP:
    114         FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints);
    115         break;
    116       case FIT_PM_IRLS:
    117         FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints, 0);
    118         break;
    119       case FIT_PLX_NOCLIP:
    120         FitPMandPar (&fitPM, fitStats->fitdataPar, fitStats->points, fitStats->Npoints);
    121         break;
    122       case FIT_PLX_IRLS:
    123         FitPMandPar_IRLS (&fitPM, fitStats->fitdataPar, fitStats->points, fitStats->Npoints, 0);
    124         break;
    125 
    126       case FIT_PLX_BOOT:
    127         FitPMandPar_IRLS (&fitPM, fitStats->fitdataPar, fitStats->points, fitStats->Npoints, 0);
     166      case FIT_POS_NOCLIP:
     167        fitResult.uR = uR;
     168        fitResult.uD = uD;
     169        if (!FitPosPMfixed_Basic (&fitResult, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) goto escape;
     170        break;
     171      case FIT_POS_IRLS:
     172        fitResult.uR = uR;
     173        fitResult.uD = uD;
     174        if (!FitPosPMfixed_IRLS (&fitResult, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) goto escape;
     175        break;
     176      case FIT_POS_BOOT:
     177        fitResult.uR = uR;
     178        fitResult.uD = uD;
     179        if (!FitPosPMfixed_IRLS (&fitResult, fitStats->fitdataPos, fitStats->points, fitStats->Npoints)) goto escape;
    128180        fitStats->Nfit = 0;
    129         int Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
     181        Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
    130182        for (k = 0; k < fitStats->NfitAlloc; k++) {
    131183          BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
    132           FitPMandPar (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask);
     184          fitStats->fit[k].uR = uR;
     185          fitStats->fit[k].uD = uD;
     186          if (!FitPosPMfixed_Basic (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, Nnomask)) continue;
    133187          fitStats->Nfit ++;
    134188        }
    135189        // 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);
     190        BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     191        BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     192        break;
     193
     194      case FIT_PM_NOCLIP:
     195        if (!FitPM_Basic (&fitResult, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) goto escape;
     196        break;
     197      case FIT_PM_IRLS:
     198        if (!FitPM_IRLS (&fitResult, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) goto escape;
     199        break;
     200      case FIT_PM_BOOT:
     201        if (!FitPM_IRLS (&fitResult, fitStats->fitdataPM, fitStats->points, fitStats->Npoints)) goto escape;
     202        fitStats->Nfit = 0;
     203        Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
     204        for (k = 0; k < fitStats->NfitAlloc; k++) {
     205          BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
     206          if (!FitPM_Basic (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, Nnomask)) continue;
     207          fitStats->Nfit ++;
     208        }
     209        // these calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above)
     210        BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     211        BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     212        BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
     213        BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
     214        break;
     215
     216      case FIT_PLX_NOCLIP:
     217        if (!FitPMandPar_Basic (&fitResult, fitStats->fitdataPar, fitStats->points, fitStats->Npoints)) goto escape;
     218        break;
     219      case FIT_PLX_IRLS:
     220        if (!FitPMandPar_IRLS (&fitResult, fitStats->fitdataPar, fitStats->points, fitStats->Npoints)) goto escape;
     221        break;
     222      case FIT_PLX_BOOT:
     223        if (!FitPMandPar_IRLS (&fitResult, fitStats->fitdataPar, fitStats->points, fitStats->Npoints)) goto escape;
     224        fitStats->Nfit = 0;
     225        Nnomask = BootstrapSaveUnmasked (fitStats->nomask, fitStats->points, fitStats->Npoints);
     226        for (k = 0; k < fitStats->NfitAlloc; k++) {
     227          BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
     228          if (!FitPMandPar_Basic (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask)) continue;
     229          fitStats->Nfit ++;
     230        }
     231        // these calls set the ERRORS on the fit parameters, not the fit values (set in IRLS above)
     232        BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_RA);
     233        BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_DEC);
     234        BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_uR);
     235        BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_uD);
     236        BootstrapRobustStats (&fitResult, fitStats->fit, fitStats->Nfit, FIT_RESULT_PLX);
    141237        break;
    142238
    … …  
    145241    }
    146242
    147     XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &fitStats->coords);
     243    XY_to_RD (&fitResult.Ro, &fitResult.Do, fitResult.Ro, fitResult.Do, &fitStats->coords);
    148244
    149245    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",
    150246             Ro, Do, uR, uD, plx,
    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);
    152 
    153     // ok (TRUE, "fitted star");
     247             fitResult.Ro, fitResult.Do, fitResult.uR, fitResult.uD, fitResult.p, Tmean, fitResult.dRo, fitResult.dDo, fitResult.duR, fitResult.duD, fitResult.dp, fitResult.Nfit, fitResult.chisq);
     248
     249    continue;
     250
     251  escape:
     252    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",
     253             NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, 0.0, NAN, NAN, NAN, NAN, NAN, 0, 0.0);
    154254  }
    155255  // return exit_status();
Note: See TracChangeset for help on using the changeset viewer.