IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39370 for trunk


Ignore:
Timestamp:
Feb 23, 2016, 5:36:56 AM (10 years ago)
Author:
eugene
Message:

more relastro work on outlier clipping

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

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relastro/Makefile

    r39365 r39370  
    200200$(BIN)/fitobj.$(ARCH): $(FITOBJ)
    201201
     202fitobj2: $(BIN)/fitobj2.$(ARCH)
     203
     204FITOBJ2 = \
     205$(SRC)/fitobj2.$(ARCH).o              \
     206$(SRC)/ParFactor.$(ARCH).o           \
     207$(SRC)/FitAstromOps.$(ARCH).o        \
     208$(SRC)/FitPosPMfixed.$(ARCH).o               \
     209$(SRC)/FitPM.$(ARCH).o               \
     210$(SRC)/FitPMandPar.$(ARCH).o               \
     211$(SRC)/UpdateObjects.$(ARCH).o               \
     212$(SRC)/BootstrapOps.$(ARCH).o        \
     213$(SRC)/mkpolyterm.$(ARCH).o            \
     214$(SRC)/GetAstromError.$(ARCH).o            \
     215$(SRC)/ImageOps.$(ARCH).o            \
     216$(SRC)/MosaicOps.$(ARCH).o            \
     217$(SRC)/extra.$(ARCH).o            \
     218$(SRC)/plotstuff.$(ARCH).o            \
     219$(SRC)/fitpoly.$(ARCH).o
     220
     221$(FITOBJ2): $(INC)/relastro.h
     222$(BIN)/fitobj2.$(ARCH): $(FITOBJ2)
     223
    202224# $(SRC)/FitPM_IRLS.$(ARCH).o          \
    203225# $(SRC)/FitPMandPar_IRLS.$(ARCH).o               \

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

    r39364 r39370  
    186186  }
    187187
    188   FitPM_SetChisq (fit, data, points, Npoints);
     188  if (!FitPM_SetChisq (fit, data, points, Npoints)) return FALSE;
    189189  return (TRUE);
    190190}
    … …  
    198198  Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
    199199
     200  int Nfit = 0;
    200201  for (i = 0; i < Npoints; i++) {
    201202    if (points[i].mask) continue; // respect the mask if set
     203    Nfit ++;
    202204
    203205    wx = points[i].Wx;
    … …  
    222224    YT += points[i].Y*TWy;
    223225  }
     226  if (Nfit < 3) return FALSE;
    224227
    225228  // X^T W X
    … …  
    290293  }
    291294   
     295  if (fit[0].Nfit < 3) {
     296    fit[0].chisq = NAN;
     297    return FALSE;
     298  }
     299
    292300  // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms)
    293301  fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms);

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

    r39364 r39370  
    189189  }
    190190
    191   FitPMandPar_SetChisq (fit, data, points, Npoints);
     191  if (!FitPMandPar_SetChisq (fit, data, points, Npoints)) return FALSE;
    192192  return (TRUE);
    193193}
    … …  
    205205  Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
    206206
     207  int Nfit = 0;
    207208  for (i = 0; i < Npoints; i++) {
    208209    if (points[i].mask) continue; // respect the mask if set
     210    Nfit ++;
    209211
    210212    wx = points[i].Wx;
    … …  
    244246    YT += points[i].Y*TWy;
    245247  }
     248  if (Nfit < 3) return FALSE;
    246249
    247250  data->A[0][0] = Wx;
    … …  
    321324  }
    322325   
     326  if (fit[0].Nfit < 3) {
     327    fit[0].chisq = NAN;
     328    return FALSE;
     329  }
     330
    323331  // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms)
    324332  fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms);
    325333
     334
    326335  return TRUE;
    327336}

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

    r39367 r39370  
    180180  }
    181181
    182   FitPosPMfixed_SetChisq (fit, data, points, Npoints);
     182  if (!FitPosPMfixed_SetChisq (fit, data, points, Npoints)) return FALSE;
    183183  return (TRUE);
    184184}
    … …  
    192192  Wx = Wy = Tx = Ty = Xs = Ys = 0.0;
    193193
     194  int Nfit = 0;
    194195  for (i = 0; i < Npoints; i++) {
    195196    if (points[i].mask) continue; // respect the mask if set
     197    Nfit ++;
    196198
    197199    wx = points[i].Wx;
    … …  
    207209    Ys += points[i].Y*wy;
    208210  }
     211  if (Nfit < 2) return FALSE;
    209212
    210213  // X^T W X
    … …  
    257260  }
    258261   
     262  if (fit[0].Nfit < 2) {
     263    fit[0].chisq = NAN;
     264    return FALSE;
     265  }
     266
    259267  // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms)
    260268  fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms);

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

    r39367 r39370  
    766766  if (USE_GALAXY_MODEL) {
    767767    for (i = 0; i < Nfit; i++) {
    768       fit->uR = average->uRgal;
    769       fit->uD = average->uDgal;
     768      fit[i].uR = average->uRgal;
     769      fit[i].uD = average->uDgal;
    770770    }
    771771  } else {
    772772    for (i = 0; i < Nfit; i++) {
    773       fit->uR = 0.0;
    774       fit->uD = 0.0;
     773      fit[i].uR = 0.0;
     774      fit[i].uD = 0.0;
    775775    }
    776776  }

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

    r39367 r39370  
    99# define MJD_REF_YRS 12.4148 /* 2012/06/01,00:00:00 */
    1010
    11 # define MJD_J2000   51544.0 /* 2010/01/01,00:00:00 */
     11# define MJD_J2000   51544.5 /* 2000/01/01,12:00:00 */
    1212
    1313# define POS_ERROR 0.010 /* arcsec */
    … …  
    109109    long A = now.tv_sec + now.tv_usec * 1000000;
    110110    srand48(A);
     111    // srand48(1);
    111112  }
    112113
    … …  
    179180  catalog->average->uRgal = uR;
    180181  catalog->average->uDgal = uD;
     182
     183  // int T2000 = ohana_date_to_sec ("2000/01/01,12:00:00");
    181184  for (i = 0; i < Npoints; i++) {
    182185    // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined)
    183     double Tmjd = MJD_MIN_MJD + drand48()*(MJD_MAX_MJD - MJD_MIN_MJD);
     186    double Tyears = MJD_MIN_YRS + drand48()*(MJD_MAX_YRS - MJD_MIN_YRS);
     187    double Tmjd   = 365.25*Tyears + MJD_J2000;
    184188
    185189    catalog->measureT[i].t = ohana_mjd_to_sec (Tmjd);
    186     double Tyears = (Tmjd - MJD_J2000) / 365.25;
    187190   
    188191    double pR, pD;
    … …  
    208211  for (i = Npoints; i < Ntotal; i++) {
    209212    // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined)
    210     double Tmjd = MJD_MIN_MJD + drand48()*(MJD_MAX_MJD - MJD_MIN_MJD);
     213    double Tyears = MJD_MIN_YRS + drand48()*(MJD_MAX_YRS - MJD_MIN_YRS);
     214    double Tmjd   = 365.25*Tyears + MJD_J2000;
    211215
    212216    catalog->measureT[i].t = ohana_mjd_to_sec (Tmjd);
    213     double Tyears = (Tmjd - MJD_J2000) / 365.25;
    214217   
    215218    double pR, pD;

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

    r39364 r39370  
    3333};
    3434
     35float ERROR_BIAS = 1.0;
     36float RA_BIAS = 0.0;
     37float DEC_BIAS = 0.0;
     38
    3539int main (int argc, char **argv) {
    3640 
    … …  
    5963    remove_argument (N, &argc, argv);
    6064    N_BOOTSTRAP = atoi (argv[N]);
     65    remove_argument (N, &argc, argv);
     66  }
     67  if ((N = get_argument (argc, argv, "-error-bias"))) {
     68    remove_argument (N, &argc, argv);
     69    ERROR_BIAS = atof (argv[N]);
     70    remove_argument (N, &argc, argv);
     71  }
     72  if ((N = get_argument (argc, argv, "-ra-bias"))) {
     73    remove_argument (N, &argc, argv);
     74    RA_BIAS = atof (argv[N]);
     75    remove_argument (N, &argc, argv);
     76  }
     77  if ((N = get_argument (argc, argv, "-dec-bias"))) {
     78    remove_argument (N, &argc, argv);
     79    DEC_BIAS = atof (argv[N]);
    6180    remove_argument (N, &argc, argv);
    6281  }
    … …  
    108127    long A = now.tv_sec + now.tv_usec * 1000000;
    109128    srand48(A);
     129    srand48(1);
    110130  }
    111131
    … …  
    184204          fitStats->fit[k].uR = uR;
    185205          fitStats->fit[k].uD = uD;
    186           if (!FitPosPMfixed_Basic (&fitStats->fit[k], fitStats->fitdataPos, fitStats->sample, Nnomask)) continue;
     206          if (!FitPosPMfixed_Basic (&fitStats->fit[fitStats->Nfit], fitStats->fitdataPos, fitStats->sample, Nnomask)) continue;
    187207          fitStats->Nfit ++;
    188208        }
    … …  
    204224        for (k = 0; k < fitStats->NfitAlloc; k++) {
    205225          BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
    206           if (!FitPM_Basic (&fitStats->fit[k], fitStats->fitdataPM, fitStats->sample, Nnomask)) continue;
     226          if (!FitPM_Basic (&fitStats->fit[fitStats->Nfit], fitStats->fitdataPM, fitStats->sample, Nnomask)) continue;
    207227          fitStats->Nfit ++;
    208228        }
    … …  
    226246        for (k = 0; k < fitStats->NfitAlloc; k++) {
    227247          BootstrapResample (fitStats->sample, fitStats->nomask, Nnomask);
    228           if (!FitPMandPar_Basic (&fitStats->fit[k], fitStats->fitdataPar, fitStats->sample, Nnomask)) continue;
     248          if (!FitPMandPar_Basic (&fitStats->fit[fitStats->Nfit], fitStats->fitdataPar, fitStats->sample, Nnomask)) continue;
    229249          fitStats->Nfit ++;
    230250        }
    … …  
    282302    points[i].D = Do + (dD + plx*pD + uD*(points[i].T - MJD_REF_YRS))/3600;
    283303
    284     points[i].dX = POS_ERROR;
    285     points[i].dY = POS_ERROR;
     304    points[i].dX = POS_ERROR*ERROR_BIAS;
     305    points[i].dY = POS_ERROR*ERROR_BIAS;
    286306  }
    287307
    … …  
    295315    ParFactor (&pR, &pD, Ro, Do, points[i].T);
    296316
    297     double dR = OUT_ERROR*(drand48() - 0.5);
    298     double dD = OUT_ERROR*(drand48() - 0.5);
     317    double dR = OUT_ERROR*(drand48() - 0.5) +  RA_BIAS;
     318    double dD = OUT_ERROR*(drand48() - 0.5) + DEC_BIAS;
    299319
    300320    points[i].R = Ro + (dR + plx*pR + uR*(points[i].T - MJD_REF_YRS))/3600/dcos(Do);
    301321    points[i].D = Do + (dD + plx*pD + uD*(points[i].T - MJD_REF_YRS))/3600;
    302322
    303     points[i].dX = POS_ERROR;
    304     points[i].dY = POS_ERROR;
     323    points[i].dX = POS_ERROR*ERROR_BIAS;
     324    points[i].dY = POS_ERROR*ERROR_BIAS;
    305325  }
    306326  fitStats->Npoints = Ntotal;

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

    r39367 r39370  
    3030end
    3131
     32macro testobj
     33  if ($0 != 1)
     34    echo "USAGE: testobj"
     35    break
     36  end
     37
     38  foreach mode pos pm plx
     39    exec ../bin/fitobj.lin64 $mode > test.$mode.dat
     40    if ("$mode" == "plx")
     41      echo "is plx"
     42      reload test.$mode.dat 1
     43    else
     44      reload test.$mode.dat 0
     45    end
     46
     47    fprintf "%5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %6.2f %5.2f" $value_dRf $value_dDf $value_duR $value_duD $value_dP $value_dRo $value_dDo $value_duRo $value_duDo $value_dPo $sigma_dRo $sigma_dDo $sigma_duRo $sigma_duDo $sigma_dPo $mean_Nfit $mean_chisq
     48  end
     49end
     50
    3251macro testnboot
    3352  if ($0 != 2)
    … …  
    3857  foreach nboot 20 30 100 300 1000
    3958    exec ../bin/fitpm.lin64 -Nstars 3000 -Noutliers 20 -Npoints 100 -Nbootstrap $nboot $1 > test.$1.dat
    40     reload test.$1.dat
     59    reload test.$1.dat 0
    4160
    4261    fprintf "%4d : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %6.2f %5.2f" $nboot $value_dRf $value_dDf $value_duR $value_duD $value_dP $value_dRo $value_dDo $value_duRo $value_duDo $value_dPo $sigma_dRo $sigma_dDo $sigma_duRo $sigma_duDo $sigma_dPo $mean_Nfit $mean_chisq
    … …  
    5069  end
    5170
    52   foreach outliers 0 3 10 30 100
     71  # foreach outliers 0 3 10 30 100
     72
     73  foreach outliers 10
    5374    exec ../bin/fitpm.lin64 -Nstars 3000 -Noutliers $outliers -Npoints 100 -Nbootstrap 100 $1 > test.$1.dat
    54     reload test.$1.dat
     75    reload test.$1.dat 0
    5576
    5677    fprintf "%4d : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %6.2f %5.2f" $outliers $value_dRf $value_dDf $value_duR $value_duD $value_dP $value_dRo $value_dDo $value_duRo $value_duDo $value_dPo $sigma_dRo $sigma_dDo $sigma_duRo $sigma_duDo $sigma_dPo $mean_Nfit $mean_chisq
    … …  
    5980
    6081macro reload
    61   if ($0 != 2)
    62     echo "USAGE: reload (file)"
     82  if ($0 != 3)
     83    echo "USAGE: reload (file) (isObjPLX)"
    6384    break
    6485  end
    6586
    6687  data $1
    67   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
     88  if ($2)
     89    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
     90    # 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
     91  else
     92    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
     93  end
    6894
    6995  foreach f R D uD uR P
    … …  
    98124end
    99125
     126macro teststats
     127  if ($0 != 2)
     128    echo "USAGE: teststats (file)"
     129    break
     130  end
     131
     132  data $1
     133  # 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
     134  # 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
     135  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
     136
     137  foreach f R D uD uR P
     138    set d$f = $f\o - $f\x
     139  end
     140
     141  set Rf = Rx - uRx*(Tx - $Tref)/3600.0/dcos(Dx)
     142  set Df = Dx - uDx*(Tx - $Tref)/3600.0
     143
     144  set dRf = 3600*(Ro - Rf)*dcos(Dx)
     145  set dDf = 3600*(Do - Df)
     146
     147  # echo "fit scatter"
     148  vstat -q dRf;  $value_dRf =  $SIGMA*1000
     149  vstat -q dDf;  $value_dDf =  $SIGMA*1000
     150  vstat -q duR;  $value_duR =  $SIGMA*1000
     151  vstat -q duD;  $value_duD =  $SIGMA*1000
     152  vstat -q dP ;  $value_dP  =  $SIGMA*1000
     153
     154  # echo "fit errors"
     155  vstat -q dRo;  $value_dRo  = $MEAN*1000; $sigma_dRo  = $SIGMA*1000
     156  vstat -q dDo;  $value_dDo  = $MEAN*1000; $sigma_dDo  = $SIGMA*1000
     157  vstat -q duRo; $value_duRo = $MEAN*1000; $sigma_duRo = $SIGMA*1000
     158  vstat -q duDo; $value_duDo = $MEAN*1000; $sigma_duDo = $SIGMA*1000
     159  vstat -q dPo ; $value_dPo  = $MEAN*1000; $sigma_dPo  = $SIGMA*1000
     160
     161  # echo "Nfit"
     162  vstat -q Nfit; $mean_Nfit = $MEAN
     163
     164  # echo chisq
     165  vstat -q chisq; $mean_chisq = $MEAN
     166
     167  fprintf "%5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %5.2f %5.2f %5.2f %5.2f %5.2f : %6.2f %5.2f" $value_dRf $value_dDf $value_duR $value_duD $value_dP $value_dRo $value_dDo $value_duRo $value_duDo $value_dPo $sigma_dRo $sigma_dDo $sigma_duRo $sigma_duDo $sigma_dPo $mean_Nfit $mean_chisq
     168
     169end
     170
    100171# determine the parallax factor (time) (ra) (dec) (parfR) (parfD) are vectors,
    101172macro getParFactor_vectors
Note: See TracChangeset for help on using the changeset viewer.