IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39610 for trunk


Ignore:
Timestamp:
Jun 24, 2016, 5:07:24 PM (10 years ago)
Author:
eugene
Message:

modify irls fitting to apply the multiply the modified weight by the error-based weight

Location:
trunk/Ohana/src
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/opihi/cmd.astro/astrom_ops.c

    r39228 r39610  
    5454  ALLOCATE (data->Wx, double, N);
    5555  ALLOCATE (data->Wy, double, N);
     56  ALLOCATE (data->Qx, double, N);
     57  ALLOCATE (data->Qy, double, N);
     58  ALLOCATE (data->qx, double, N);
     59  ALLOCATE (data->qy, double, N);
    5660  ALLOCATE (data->index, int, N);
    5761  return TRUE;
    … …  
    6872  FREE (data->Wx);
    6973  FREE (data->Wy);
     74  FREE (data->Qx);
     75  FREE (data->Qy);
     76  FREE (data->qx);
     77  FREE (data->qy);
    7078  FREE (data->index);
    7179}

  • trunk/Ohana/src/opihi/cmd.astro/fitplx.c

    r39234 r39610  
    597597  dsort (values, Npts);
    598598
    599   if (Npts % 2) {
    600     int Ncenter = Npts / 2;
    601     *median = values[Ncenter];
    602   } else {
    603     int Ncenter = Npts / 2 - 1;
    604     *median = 0.5*(values[Ncenter] + values[Ncenter + 1]);
    605   }
    606 
    607   double Slo = VectorFractionInterpolate (values, 0.158655, Npts);
    608   double Shi = VectorFractionInterpolate (values, 0.841345, Npts);
    609 
    610   *sigma = (Shi - Slo) / 2.0;
    611 
     599  if (median) {
     600    if (Npts % 2) {
     601      int Ncenter = Npts / 2;
     602      *median = values[Ncenter];
     603    } else {
     604      int Ncenter = Npts / 2 - 1;
     605      *median = 0.5*(values[Ncenter] + values[Ncenter + 1]);
     606    }
     607  }
     608   
     609  if (sigma) {
     610    double Slo = VectorFractionInterpolate (values, 0.158655, Npts);
     611    double Shi = VectorFractionInterpolate (values, 0.841345, Npts);
     612    *sigma = (Shi - Slo) / 2.0;
     613  }
     614 
    612615  return TRUE;
    613616}

  • trunk/Ohana/src/opihi/cmd.astro/fitplx_irls.c

    r39596 r39610  
    119119                         fitdata.t, fitdata.pX, fitdata.pY,
    120120                         fitdata.Wx, fitdata.Wy,
     121                         fitdata.Qx, fitdata.Qy,
     122                         fitdata.qx, fitdata.qy,
    121123                         fitdata.Npts, max_iterations, outlier_limit, binning_step, VERBOSE)) {
    122124    return FALSE;
    … …  
    204206    }
    205207
    206     // now calculate median and sigma for each vector
    207     VectorRobustStats (pvec,  &fit.p,  &fit.dp);
    208     VectorRobustStats (uRvec, &fit.uR, &fit.duR);
    209     VectorRobustStats (uDvec, &fit.uD, &fit.duD);
    210     VectorRobustStats (Rvec,  &fit.Ro, &fit.dRo);
    211     VectorRobustStats (Dvec,  &fit.Do, &fit.dDo);
     208    // now calculate robust sigma for each vector
     209    VectorRobustStats (pvec,  NULL, &fit.dp);
     210    VectorRobustStats (uRvec, NULL, &fit.duR);
     211    VectorRobustStats (uDvec, NULL, &fit.duD);
     212    VectorRobustStats (Rvec,  NULL, &fit.dRo);
     213    VectorRobustStats (Dvec,  NULL, &fit.dDo);
     214
     215    PlxFitDataFree (&sample);   
    212216  }
    213217
    … …  
    250254  set_variable ("Nfit",  fit.Nfit);
    251255
    252 
     256  PlxFitDataFree (&fitdata);
    253257 
    254258  return (TRUE);
    … …  
    257261/* do we want an init function which does the alloc and a clear function to free? */
    258262int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD,
    259                       double *Wx, double *Wy,
     263                      double *Wx, double *Wy, double *Qx, double *Qy, double *qx, double *qy,
    260264                      int Npts, int max_iterations, double outlier_limit, double binning_step, int VERBOSE) {
    261265
    … …  
    270274
    271275  double sigma_ols, sigma_hat; // Sigma estimates.
    272   //  double *Wx, *Wy;             // Weight vectors.  Not errors.
    273276  double *rx, *ry;             // Deviation from model
    274277  double *u;                   // Deviation magnitude
    … …  
    305308  // Convert the measurement errors into the initial weights.
    306309  for (i = 0; i < Npts; i++) {
    307     Wx[i] = (fabs(dX[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dX[i]);
    308     Wy[i] = (fabs(dY[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dY[i]);
     310    Qx[i] = (fabs(dX[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dX[i]);
     311    Qy[i] = (fabs(dY[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dY[i]);
    309312  }
    310313 
    … …  
    313316  // Solve OLS equation
    314317  if (binning_step == 0.0) {
    315     if (!weighted_LS_PLX(T,pR,pD,X,Wx,Y,Wy,Npts,
    316                          A,B,VERBOSE)) {
     318    if (!weighted_LS_PLX(T,pR,pD,X,Qx,Y,Qy,Npts,A,B,VERBOSE)) {
    317319      // Handle fail case
    318320      return(FALSE);
    … …  
    329331    double *WYbin = NULL;
    330332
    331     if (!bin_points_PLX(T, pR, pD, X, Wx, Y, Wy, Npts,
     333    if (!bin_points_PLX(T, pR, pD, X, Qx, Y, Qy, Npts,
    332334                        &Tbin, &pRbin, &pDbin, &Xbin, &WXbin, &Ybin, &WYbin, &Nbins,  binning_step) ) {
    333335      return(FALSE);
    … …  
    396398      Wx[i] = weight_cauchy(rx[i] / dX[i]);
    397399      Wy[i] = weight_cauchy(ry[i] / dY[i]);
     400      qx[i] = Qx[i] * Wx[i];
     401      qy[i] = Qy[i] * Wy[i];
    398402    }
    399403
    400404    // Solve
    401     if (!weighted_LS_PLX(T,pR,pD,X,Wx,Y,Wy,Npts,A,B,VERBOSE)) {
     405    if (!weighted_LS_PLX(T,pR,pD,X,qx,Y,qy,Npts,A,B,VERBOSE)) {
    402406      // Handle fail case
    403407      return(FALSE);

  • trunk/Ohana/src/opihi/include/astro.h

    r39585 r39610  
    2020  double *dX;
    2121  double *dY;
    22   double *Wx;
    23   double *Wy;
     22  double *Wx; // modified weight based on distance from fit
     23  double *Wy; // modified weight based on distance from fit
     24  double *Qx; // raw error-based weight (1/dX^2)
     25  double *Qy;
     26  double *qx; // modified error-based weight (1/dX^2) * Wx
     27  double *qy;
    2428  int *index;
    2529  int Npts;
    … …  
    5963
    6064int FitPMonly_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts, double binning_step, int VERBOSE);
    61 int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, double *Wx, double *Wy, int Npts, int max_iterations, double outlier_limit, double binning_step, int VERBOSE);
     65int FitPMandPar_IRLS (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, double *Wx, double *Wy, double *Qx, double *Qy, double *qx, double *qy, int Npts, int max_iterations, double outlier_limit, double binning_step, int VERBOSE);
    6266
    6367int weighted_LS_PLX (double *T, double *pR, double *pD, double *X, double *WX, double *Y, double *WY, int Npts, double **A, double **B, int VERBOSE);

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

    r39600 r39610  
    163163  double D, dD;
    164164  double T, dT;
    165   double Wx, Wy;
     165  double Qx, Qy; // unmodified error-based weight (1/dX^2)
     166  double qx, qy; // IRLS-modified error-based weight (1/dX^2)
     167  double Wx, Wy; // IRLS weight factor
    166168  double rx, ry;
    167169  double u;

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

    r39457 r39610  
    227227  object->Wy     = 1.0;
    228228
     229  object->Qx     = 0.0;
     230  object->Qy     = 0.0;
     231  object->qx     = 0.0;
     232  object->qy     = 0.0;
     233
    229234  object->rx     = 0.0;
    230235  object->ry     = 0.0;

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

    r39371 r39610  
    11# include "relastro.h"
     2# define OLD_METHOD 0
     3
    24int FitPMandPar_MinChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
    35int FitPMandPar_SetChisq (FitAstromResult *fit, FitAstromData *data, FitAstromPoint *points, int Npoints);
    … …  
    1618  // Convert the measurement errors into initial weights.
    1719  for (i = 0; i < Npoints; i++) {
     20# if (OLD_METHOD)
    1821    points[i].Wx = 1.0 / SQ(points[i].dX);
    1922    points[i].Wy = 1.0 / SQ(points[i].dY);
     23# else
     24    points[i].Wx = 1.0;
     25    points[i].Wy = 1.0;
     26    points[i].Qx = 1.0 / SQ(points[i].dX);
     27    points[i].Qy = 1.0 / SQ(points[i].dY);
     28    points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0
     29    points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0
     30# endif
    2031  }
    2132
    … …  
    3445  // Convert the measurement errors into initial weights.
    3546  for (i = 0; i < Npoints; i++) {
     47# if (OLD_METHOD)
    3648    points[i].Wx = 1.0 / SQ(points[i].dX);
    3749    points[i].Wy = 1.0 / SQ(points[i].dY);
     50    points[i].qx = points[i].Wx; // Wx, Wy start out at 1.0
     51    points[i].qy = points[i].Wy; // Wx, Wy start out at 1.0
     52# else
     53    points[i].Qx = 1.0 / SQ(points[i].dX);
     54    points[i].Qy = 1.0 / SQ(points[i].dY);
     55    points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0
     56    points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0
     57# endif
    3858  }
    3959 
    … …  
    7797      points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
    7898      points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
     99      points[i].qx = points[i].Wx * points[i].Qx;
     100      points[i].qy = points[i].Wy * points[i].Qy;
    79101    }   
    80102
    … …  
    209231    Nfit ++;
    210232
     233# if (OLD_METHOD)   
    211234    wx = points[i].Wx;
    212235    wy = points[i].Wy;
     236# else
     237    wx = points[i].qx;
     238    wy = points[i].qy;
     239# endif
    213240
    214241    Wx += wx;

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

    r39364 r39610  
    11# include "relastro.h"
     2# define OLD_METHOD 1
    23
    34// These should probably be tunable:
    … …  
    1516  // Convert the measurement errors into initial weights.
    1617  for (i = 0; i < Npoints; i++) {
     18# if (OLD_METHOD)
    1719    points[i].Wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1 / SQ(points[i].dX);
    1820    points[i].Wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1 / SQ(points[i].dY);
     21# else
     22    points[i].Qx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1 / SQ(points[i].dX);
     23    points[i].Qy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1 / SQ(points[i].dY);
     24# endif
     25    points[i].qx = points[i].Wx * points[i].Qx; // Wx, Wy start out at 1.0
     26    points[i].qy = points[i].Wy * points[i].Qy; // Wx, Wy start out at 1.0
    1927  }
    2028 
    … …  
    5866      points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
    5967      points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
     68      points[i].qx = points[i].Wx * points[i].Qx;
     69      points[i].qy = points[i].Wy * points[i].Qy;
    6070    }   
    6171
    … …  
    138148      ay += dpsi_cauchy(points[i].ry / points[i].dY);
    139149
    140       bx += SQ(points[i].Wx);
    141       by += SQ(points[i].Wy);
     150      bx += SQ(points[i].qx);
     151      by += SQ(points[i].qy);
    142152    }
    143153    ax /= 1.0 * Npoints;  // mean(psi_dot(r))
    … …  
    207217    // if (VERBOSE == 2) fprintf (stderr, "%f %f : %f %f : %f : %f %f\n", X[i], WX[i], Y[i], WY[i], T[i], pR[i], pD[i]);
    208218
     219# if (OLD_METHOD)
    209220    wx = points[i].Wx;
    210221    wy = points[i].Wy;
     222# else
     223    wx = points[i].Wx;
     224    wy = points[i].Wy;
     225# endif
    211226
    212227    Wx += wx;

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

    r39457 r39610  
    121121  for (i = 0; i < Nstars; i++) {
    122122
    123     double Ro  = 360.0*(drand48() - 0.5);
     123    double Ro  = 360.0*drand48();
    124124    double Phi = 2.0*(drand48() - 0.5);
    125125    double Do  = DEG_RAD * asin(Phi);

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

    r39371 r39610  
    155155  for (i = 0; i < Nstars; i++) {
    156156
    157     double Ro  = 360.0*(drand48_cnt() - 0.5);
     157    double Ro  = 360.0*drand48_cnt();
    158158    double Phi = 2.0*(drand48_cnt() - 0.5);
    159159    double Do  = DEG_RAD * asin(Phi);
    … …  
    181181    }
    182182
    183     int Npoints = N_POINTS*drand48_cnt(); // minimum of 0 points
    184     int Noutliers_max = Npoints * F_OUTLIERS;
    185     int Noutliers = Noutliers_max*drand48_cnt(); // minimum of 0 points
     183    // int Npoints = N_POINTS*drand48_cnt(); // minimum of 0 points
     184    // int Noutliers_max = Npoints * F_OUTLIERS;
     185    // int Noutliers = Noutliers_max*drand48_cnt(); // minimum of 0 points
     186
     187    int Npoints = N_POINTS; // minimum of 0 points
     188    int Noutliers = Npoints * F_OUTLIERS;
    186189
    187190    // generate a single fake star with N_POINTS real points and N_OUTLIERS bad points
    … …  
    287290  escape:
    288291    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",
    289              NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, NAN, 0.0, NAN, NAN, NAN, NAN, NAN, 0, 0.0);
     292             Ro, Do, uR, uD, plx,
     293             NAN, NAN, NAN, NAN, NAN, 0.0, NAN, NAN, NAN, NAN, NAN, 0, 0.0);
    290294  }
    291295  // return exit_status();

  • trunk/Ohana/src/relastro/test/mana.sh

    r39375 r39610  
    9191  # foreach outliers 0 3 10 30 100
    9292
    93   foreach outliers 10
    94     exec ../bin/fitpm.lin64 -Nstars 3000 -Noutliers $outliers -Npoints 100 -Nbootstrap 100 $1 > test.$1.dat
     93  foreach outliers 0
     94    echo ../bin/fitpm.lin64 -Nstars 3000 -Npoints 100 -Nbootstrap 100 $1
     95    exec ../bin/fitpm.lin64 -Nstars 3000 -Npoints 100 -Nbootstrap 100 $1 > test.$1.dat
    9596    reload test.$1.dat 0
    9697
    … …  
    107108  data $1
    108109  if ($2)
    109     read Npts 1 Nout 2 Ro 3 Do 4 uRo 5 uDo 6 Po 7 Rx 9 Dx 10 uRx 11 uDx 12 Px 13 Tx 15 dRo 17 dDo 18 duRo 19 duDo 20 dPo 21 Nfit 22 chisq 26
    110     # read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 24
     110    # read Npts 1 Nout 2 Ro 3 Do 4 uRo 5 uDo 6 Po 7 Rx 9 Dx 10 uRx 11 uDx 12 Px 13 Tx 15 dRo 17 dDo 18 duRo 19 duDo 20 dPo 21 Nfit 22 chisq 26
     111    read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 24
    111112  else
    112113    read Ro 1 Do 2 uRo 3 uDo 4 Po 5 Rx 7 Dx 8 uRx 9 uDx 10 Px 11 Tx 13 dRo 15 dDo 16 duRo 17 duDo 18 dPo 19 Nfit 20 chisq 22
Note: See TracChangeset for help on using the changeset viewer.