IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39241 for trunk


Ignore:
Timestamp:
Dec 9, 2015, 3:41:43 PM (11 years ago)
Author:
eugene
Message:

added IRLS versions of FitPM and FitPMandPar, added test to prove they are working

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

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relastro/Makefile

    r39239 r39241  
    170170$(SRC)/FitPM.$(ARCH).o               \
    171171$(SRC)/FitPM_IRLS.$(ARCH).o               \
     172$(SRC)/FitPMandPar.$(ARCH).o               \
     173$(SRC)/FitPMandPar_IRLS.$(ARCH).o               \
    172174$(SRC)/mkpolyterm.$(ARCH).o            \
    173175$(SRC)/fitpoly.$(ARCH).o

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

    r39238 r39241  
    750750
    751751int FitPM_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE);
     752int FitPMandPar_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE);
    752753
    753754double MedianAbsDeviation(FitAstromPoint *points, int Npoints);
    754755
    755756int weighted_LS_PM (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE);
     757int weighted_LS_PLX (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE);
    756758
    757759double weight_cauchy (double x);

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

    r39239 r39241  
    225225  return;
    226226}
     227
     228double weight_cauchy (double x) {
     229  double r = x / 2.385;
     230  return (1.0 / (1.0 + SQ(r)));
     231}
     232
     233// dpsi = (d/dx) (x * weight(x))
     234double dpsi_cauchy (double x) {
     235  double r2 = SQ(x / 2.385);
     236  return ((1.0 - r2) / (SQ(1 + r2)));
     237}
     238
     239
     240// median absolute deviation
     241// MAD = median(abs(x - median(x)))
     242double MedianAbsDeviation(FitAstromPoint *points, int Npoints) {
     243
     244  double *x;
     245  double median = 0.0;
     246  int i;
     247 
     248  ALLOCATE(x, double, Npoints);
     249  for (i = 0; i < Npoints; i++) {
     250    x[i] = points[i].u;
     251  }
     252  dsort(x, Npoints);
     253
     254  if (Npoints % 2) {
     255    median = 0.5*(x[(int)(0.5*Npoints)] + x[(int)(0.5*Npoints) - 1]);
     256  } else {
     257    median = x[(int)(0.5*Npoints)];
     258  }
     259
     260  for (i = 0; i < Npoints; i++ ) {
     261    x[i] = fabs(x[i] - median);
     262  }
     263  dsort(x, Npoints);
     264
     265  if (Npoints % 2) {
     266    median = 0.5*(x[(int)(0.5*Npoints)] + x[(int)(0.5*Npoints) - 1]);
     267  } else {
     268    median = x[(int)(0.5*Npoints)];
     269  }
     270
     271  return median;
     272}

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

    r39238 r39241  
    44# define MAX_ITERATIONS 10
    55# define FIT_TOLERANCE 1e-4
     6# define FLT_TOLERANCE 1e-6
    67# define WEIGHT_THRESHOLD 0.3
    78
    8 /* do we want an init function which does the alloc and a clear function to free? */
    99int FitPM_IRLS (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints, int VERBOSE) {
    1010
    1111  int i,j;
    1212
    13   int dof = 2 * Npoints - 4;
    1413  int p   = 4;
    1514  int n   = 2 * Npoints;
     15  int dof = n - p;
    1616 
    1717  // data->A,B,Cov,Beta,Beta_prev are allocated outside by FitAstromDataInit()
    … …  
    2020  // Convert the measurement errors into initial weights.
    2121  for (i = 0; i < Npoints; i++) {
    22     points[i].Wx = 1 / points[i].dX;
    23     points[i].Wy = 1 / points[i].dY;
     22    // points[i].Wx = 1 / points[i].dX;
     23    // points[i].Wy = 1 / points[i].dY;
     24    points[i].Wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1 / SQ(points[i].dX);
     25    points[i].Wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1 / SQ(points[i].dY);
    2426  }
    2527 
    … …  
    4042
    4143  // Save OLS covariance and Beta (solution vector, which is actually also saved in fit)
    42   for (i = 0; i < 4; i++) {
    43     for (j = 0; j < 4; j++) {
     44  for (i = 0; i < data->Nterms; i++) {
     45    for (j = 0; j < data->Nterms; j++) {
    4446      data->Cov[i][j] = data->A[i][j];
    4547    }
    … …  
    4749  }
    4850
    49   // Iterately reweight and solve
     51  // Iteratively reweight and solve
    5052  double sigma_hat = 0.0; // save for the error model
    5153  int converged = FALSE;
    … …  
    5658  for (iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations ++) {
    5759    // Save Beta.
    58     for (i = 0; i < 4; i ++) {
     60    for (i = 0; i < data->Nterms; i ++) {
    5961      data->Beta_prev[i] = data->Beta[i];
    6062    }
    … …  
    7375
    7476    // store the new Beta.
    75     for (i = 0; i < 4; i++) {
     77    for (i = 0; i < data->Nterms; i++) {
    7678      data->Beta[i] = data->B[i][0];
    7779    }
    … …  
    8789    // Check convergence
    8890    converged = TRUE;
    89     for (i = 0; i < 4; i++) {
    90       if (fabs(data->Beta[i] - data->Beta_prev[i]) > FIT_TOLERANCE * fabs(data->Beta[i])) {
     91    for (i = 0; i < data->Nterms; i++) {
     92      // if we are within FIT_TOLERANCE as a fractional error or FLT_TOLERANCE as an absolute error, we are good
     93      if ((fabs(data->Beta[i] - data->Beta_prev[i]) > FIT_TOLERANCE * fabs(data->Beta[i])) &&
     94          (fabs(data->Beta[i] - data->Beta_prev[i]) > FLT_TOLERANCE)) {
    9195        converged = FALSE;
    9296      }
    … …  
    134138    fit[0].duD = sqrt(data->Cov[3][3]);
    135139
    136     fit[9].dRo *= sigma_final_x;
    137     fit[9].duR *= sigma_final_x;
    138 
    139     fit[9].dDo *= sigma_final_y;
    140     fit[9].duD *= sigma_final_y;
     140    fit[0].dRo *= sigma_final_x;
     141    fit[0].duR *= sigma_final_x;
     142    fit[0].dDo *= sigma_final_y;
     143    fit[0].duD *= sigma_final_y;
    141144  }
    142145
    … …  
    154157  }
    155158
    156   // add up the chi square for the fit, only counting the unmasked points
     159  // (optionally) add up the chi square for the fit, only counting the unmasked points
    157160  double chisq = 0.0;
    158161  fit[0].Nfit = 0;
    159162  for (i = 0; i < Npoints; i++) {
    160163    if (points[i].mask) continue;
    161 
     164     
    162165    double Xf = fit[0].Ro + fit[0].uR*points[i].T;
    163166    double Yf = fit[0].Do + fit[0].uD*points[i].T;
    164     chisq += SQ(points[i].X - Xf) / SQ(points[i].dX);
    165     chisq += SQ(points[i].Y - Yf) / SQ(points[i].dY);
     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;
    166171    fit[0].Nfit ++;
    167172  }
    168 
     173   
    169174  // the reduced chisq is divided by (Ndof = 2*Nfit - 4)
    170175  fit[0].chisq = chisq / (2.0*fit[0].Nfit - 4.0);
     176
    171177  return (TRUE);
    172178}
    … …  
    175181
    176182  int i;
    177   double Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
     183  double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
    178184  Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
    179185
    180186  for (i = 0; i < Npoints; i++) {
    181     Wx += points[i].Wx;
    182     Wy += points[i].Wy;
    183 
    184     double TWx = points[i].T*points[i].Wx;
    185     double TWy = points[i].T*points[i].Wy;
     187    wx = points[i].Wx;
     188    wy = points[i].Wy;
     189
     190    Wx += wx;
     191    Wy += wy;
     192
     193    double TWx = points[i].T*wx;
     194    double TWy = points[i].T*wy;
    186195
    187196    Tx += TWx;
    … …  
    191200    Ty2 += points[i].T*TWy;
    192201   
    193     Xs += points[i].X*points[i].Wx;
    194     Ys += points[i].Y*points[i].Wy;
     202    Xs += points[i].X*wx;
     203    Ys += points[i].Y*wy;
    195204
    196205    XT += points[i].X*TWx;
    … …  
    237246  return TRUE;
    238247}
    239 
    240 double weight_cauchy (double x) {
    241   double r = x / 2.385;
    242   return (1.0 / (1.0 + SQ(r)));
    243 }
    244 
    245 // dpsi = (d/dx) (x * weight(x))
    246 double dpsi_cauchy (double x) {
    247   double r2 = SQ(x / 2.385);
    248   return ((1.0 - r2) / (SQ(1 + r2)));
    249 }
    250 
    251 
    252 // median absolute deviation
    253 // MAD = median(abs(x - median(x)))
    254 double MedianAbsDeviation(FitAstromPoint *points, int Npoints) {
    255 
    256   double *x;
    257   double median = 0.0;
    258   int i;
    259  
    260   ALLOCATE(x, double, Npoints);
    261   for (i = 0; i < Npoints; i++) {
    262     x[i] = points[i].u;
    263   }
    264   dsort(x, Npoints);
    265 
    266   if (Npoints % 2) {
    267     median = 0.5*(x[(int)(0.5*Npoints)] + x[(int)(0.5*Npoints) - 1]);
    268   } else {
    269     median = x[(int)(0.5*Npoints)];
    270   }
    271 
    272   for (i = 0; i < Npoints; i++ ) {
    273     x[i] = fabs(x[i] - median);
    274   }
    275   dsort(x, Npoints);
    276 
    277   if (Npoints % 2) {
    278     median = 0.5*(x[(int)(0.5*Npoints)] + x[(int)(0.5*Npoints) - 1]);
    279   } else {
    280     median = x[(int)(0.5*Npoints)];
    281   }
    282 
    283   return median;
    284 }

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

    r39239 r39241  
    11# include "relastro.h"
    22
    3 double rnd_gauss (double mean, double sigma);
    4 void gauss_init (int Nbin);
    5 double gaussian (double x, double mean, double sigma);
    63int mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad);
    74
    … …  
    1512# define OUT_ERROR 0.500 /* arcsec */
    1613
     14# define N_STARS 20000
    1715# define N_POINTS 100
    1816# define N_OUTLIERS 10
    … …  
    2119# define dsin(THETA) sin(RAD_DEG*THETA)
    2220
     21enum {
     22  FIT_PM_NONE,
     23  FIT_PM_IRLS,
     24  FIT_PM_NOCLIP,
     25  FIT_PLX_IRLS,
     26  FIT_PLX_NOCLIP,
     27};
     28
    2329int main (int argc, char **argv) {
    2430 
     31  if (argc != 2) {
     32    fprintf (stderr, "USAGE: %s (mode)\n", argv[0]);
     33    exit (2);
     34  }
     35
     36  int mode = FIT_NONE;
     37  if (!strcasecmp(argv[1], "pm")) {
     38    mode = FIT_PM_NOCLIP;
     39  }
     40  if (!strcasecmp(argv[1], "pm-irls")) {
     41    mode = FIT_PM_IRLS;
     42  }
     43  if (!strcasecmp(argv[1], "plx")) {
     44    mode = FIT_PLX_NOCLIP;
     45  }
     46  if (!strcasecmp(argv[1], "plx-irls")) {
     47    mode = FIT_PLX_IRLS;
     48  }
     49
    2550  // plan_tests (14);
    2651
    … …  
    3459    srand48(A);
    3560  }
    36   gauss_init (2048);
     61
     62  gaussdev_init ();
    3763
    3864  FitStats *fitStats = FitStatsInit (N_POINTS + N_OUTLIERS, 0);
    … …  
    4066  int i;
    4167
    42   // XXX try a single star to see if it looks good:
    43   if (0) {
    44     mkstar (fitStats, 3.0, 2.0, 0.2, 0.4, 0.1, N_POINTS, N_OUTLIERS);
    45     for (i = 0; i < fitStats->Npoints; i++) {
    46       fprintf (stdout, "%f %f\n", 3600*(fitStats->points[i].R - 3.0), 3600*(fitStats->points[i].D - 2.0));
    47     }
    48     exit (0);
    49   }
    50 
    51   int Nstars = 1000;
     68  int Nstars = N_STARS;
    5269  for (i = 0; i < Nstars; i++) {
    5370
    … …  
    5875    double uR  = 0.5*(drand48() - 0.5);
    5976    double uD  = 0.5*(drand48() - 0.5);
    60     double plx = 0.0*(drand48() + 0.0);
     77
     78    double plx;
     79    switch (mode) {
     80      case FIT_PM_IRLS:
     81      case FIT_PM_NOCLIP:
     82        plx = 0.0*(drand48() + 0.0);
     83        break;
     84      case FIT_PLX_IRLS:
     85      case FIT_PLX_NOCLIP:
     86        plx = 0.5*(drand48() + 0.0);
     87        break;
     88      default:
     89        myAbort("oops");
     90    }
     91
    6192
    6293    mkstar (fitStats, Ro, Do, uR, uD, plx, N_POINTS, N_OUTLIERS);
    … …  
    6899    FitAstromResultInit (&fitPM);
    69100
    70     if (0) {
    71       FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints);
    72     } else {
    73       FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints, 0);
     101    switch (mode) {
     102      case FIT_PM_NOCLIP:
     103        FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints);
     104        break;
     105      case FIT_PM_IRLS:
     106        FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints, 0);
     107        break;
     108      case FIT_PLX_NOCLIP:
     109        FitPMandPar (&fitPM, fitStats->fitdataPar, fitStats->points, fitStats->Npoints);
     110        break;
     111      case FIT_PLX_IRLS:
     112        FitPMandPar_IRLS (&fitPM, fitStats->fitdataPar, fitStats->points, fitStats->Npoints, 0);
     113        break;
     114      default:
     115        myAbort("oops");
    74116    }
    75117
    … …  
    83125  }
    84126  // return exit_status();
     127  gaussdev_free();
     128
    85129  exit (0);
    86130}
    … …  
    103147    ParFactor (&pR, &pD, Ro, Do, points[i].T);
    104148
    105     double dR = rnd_gauss(0.0, POS_ERROR);
    106     double dD = rnd_gauss(0.0, POS_ERROR);
     149    double dR = gaussdev_rnd(0.0, POS_ERROR);
     150    double dD = gaussdev_rnd(0.0, POS_ERROR);
    107151
    108152    points[i].R = Ro + (dR + plx*pR + uR*(points[i].T - MJD_REF_YRS))/3600/dcos(Do);
    … …  
    135179  return TRUE;
    136180}
    137 
    138 static int Ngaussint = 0;
    139 static double *gaussint;
    140 
    141 extern double drand48();
    142 
    143 double gaussian (double x, double mean, double sigma) {
    144 
    145   double f;
    146 
    147   f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
    148 
    149   return (f);
    150 
    151 }
    152 
    153 /* integrate a gaussian from -5 sigma to +5 sigma */
    154 void gauss_init (int Nbin) {
    155  
    156   int i;
    157   double val, x, dx, dx1, dx2, dx3, df;
    158   double mean, sigma;
    159  
    160   /* no need to generate this if it already exists */
    161   if (Ngaussint == Nbin) return;
    162 
    163   // A = time(NULL);
    164   // // XXX this is expensive if called a lot (1 sec min)
    165   // // for (B = 0; A == time(NULL); B++);
    166   // B = A + 10000;
    167   // srand48(B);
    168  
    169   Ngaussint = Nbin;
    170   ALLOCATE (gaussint, double, Ngaussint + 1);
    171 
    172   val = 0;
    173   dx = 1.0 / Ngaussint;
    174   dx1 = dx / 3.0;
    175   dx2 = 2.0*dx/3.0;
    176   dx3 = dx;
    177   mean = 0.0;
    178   sigma = 1.0;
    179  
    180   for (i = 0, x = -7.0; (i < Ngaussint) && (x < 7.0); x += dx)  {
    181     df = (3.0*gaussian(x    , mean, sigma) +
    182           9.0*gaussian(x+dx1, mean, sigma) +
    183           9.0*gaussian(x+dx2, mean, sigma) +
    184           3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0);
    185     val += df;
    186     if (val > (i + 0.5) / (double) Ngaussint) {
    187       gaussint[i] = x + dx / 2.0;
    188       i++;
    189     }
    190   }
    191 }
    192 
    193 double rnd_gauss (double mean, double sigma) {
    194  
    195   int i;
    196   double y;
    197  
    198   y = drand48();
    199   i = Ngaussint*y;
    200   y = gaussint[i]*sigma + mean;
    201  
    202   return (y);
    203  
    204 }
    205  
Note: See TracChangeset for help on using the changeset viewer.