IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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/relastro
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • 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.