IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 27, 2015, 4:49:06 PM (11 years ago)
Author:
eugene
Message:

extensive work on relphot, relastro, uniphot, dvomerge aiming to the construction and calibration of PV3

Location:
trunk/Ohana
Files:
18 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana

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

    r37807 r38986  
    77
    88typedef struct {
     9  double *X;
     10  double *Y;
     11  double *t;
     12  double *pX;
     13  double *pY;
     14  double *dX;
     15  double *dY;
     16  int *index;
     17  int Npts;
     18} PlxFitData;
     19
     20typedef struct {
    921  double Ro, dRo;
    1022  double Do, dDo;
    … …  
    1729  double chisq;
    1830  int Nfit;
     31  int getChisq;
    1932} PlxFit;
     33
     34int VectorRobustStats (Vector *vector, double *median, double *sigma);
     35double VectorFractionInterpolate (double *values, float fraction, int Npts);
     36
     37int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal);
     38int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean);
     39int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE);
     40
     41int PlxFitDataAlloc (PlxFitData *data, int N);
     42void PlxFitDataFree (PlxFitData *data);
     43int PlxBootstrapResample (PlxFitData *src, PlxFitData *tgt);
    2044
    2145int FitPMandPar (PlxFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts, int VERBOSE);
    … …  
    4771  }
    4872
     73  int Noutlier = 0;
     74  float dPsigMax = FLT_MAX;
     75  if ((N = get_argument (argc, argv, "-outlier-tests"))) {
     76    remove_argument (N, &argc, argv);
     77    Noutlier = atoi(argv[N]);
     78    remove_argument (N, &argc, argv);
     79    dPsigMax = atof(argv[N]);
     80    remove_argument (N, &argc, argv);
     81  }
     82
     83  int Nresample = 0;
     84  if ((N = get_argument (argc, argv, "-bootstrap-resample"))) {
     85    remove_argument (N, &argc, argv);
     86    Nresample = atoi(argv[N]);
     87    remove_argument (N, &argc, argv);
     88  }
     89
     90  Vector *dPvec = NULL;
     91  if ((N = get_argument (argc, argv, "-dPsig"))) {
     92    if (!Noutlier) { gprint (GP_ERR, "-dPsig requires -outlier-tests to be non-zero\n"); return FALSE; }
     93    remove_argument (N, &argc, argv);
     94    if (!(dPvec = SelectVector (argv[N], ANYVECTOR, TRUE))) return FALSE;
     95    remove_argument (N, &argc, argv);
     96  }
     97
    4998  if (argc != 6) {
    50     gprint (GP_ERR, "USAGE: fitplx (ra) (dR) (dec) (dD) (mjd) [-mask mask]\n");
    51     // what about the errors?
     99    gprint (GP_ERR, "USAGE: fitplx (ra) (dR) (dec) (dD) (mjd) [-mask mask] [-v] [-vv]\n");
     100    gprint (GP_ERR, "  -outlier-tests Nsamples dPsigMax : run Nsample bootstrap-resamples to define the path deviations and reject based on dPsigMax\n");
     101    gprint (GP_ERR, "  -dPsig vec : save path deviations in vec\n");
    52102    return (FALSE);
    53103  }
    … …  
    77127  }
    78128
    79   N = tvec->Nelements; // XXX check other lengths
    80 
    81   // find mean values to remove
    82   double Npts = 0;
    83   double Tmean = 0;
    84   double Rmean = 0;
    85   double Dmean = 0;
    86   double Tmin = +1000000;
    87   double Tmax = -1000000;
    88   for (i = 0; i < N; i++) {
    89     if (mask && !mask[i]) continue;
    90     Rmean += R[i];
    91     Dmean += D[i];
    92     Tmean += T[i];
    93     Tmin = MIN(Tmin, T[i]);
    94     Tmax = MAX(Tmax, T[i]);
    95     Npts += 1.0;
    96   }
    97   Rmean /= Npts;
    98   Dmean /= Npts;
    99   Tmean /= Npts;
    100 
    101   float Trange = Tmax - Tmin;
    102   // fprintf (stderr, "R,D : %f,%f, T: %f, Trange: %f, Tmin: %f, Tmax: %f\n", Rmean, Dmean, Tmean, Trange, Tmin, Tmax);
     129  // Ntotal : all points supplied by user
     130  // Nsubset : unmasked points
     131  int Ntotal = tvec->Nelements; // XXX check other lengths
     132  if (dPvec) ResetVector (dPvec, OPIHI_FLT, Ntotal);
     133
     134  double Rmean, Dmean, Tmean;
     135  PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
    103136
    104137  /* project coordinates to a plane centered on the object with units of arcsec */
    … …  
    109142  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
    110143
    111   double *X, *Y, *t, *pX, *pY, *dX, *dY;
    112   ALLOCATE (X, double, N);
    113   ALLOCATE (Y, double, N);
    114   ALLOCATE (dX, double, N);
    115   ALLOCATE (dY, double, N);
    116   ALLOCATE (t, double, N);
    117   ALLOCATE (pX, double, N);
    118   ALLOCATE (pY, double, N);
    119 
    120   float pXmin = +2.0;
    121   float pXmax = -2.0;
    122   float pYmin = +2.0;
    123   float pYmax = -2.0;
    124 
    125   int n = 0;
    126   for (i = 0; i < N; i++) {
    127     if (mask && !mask[i]) continue;
    128     RD_to_XY (&X[n], &Y[n], R[i], D[i], &coords);
    129     dX[n] = dR[i];
    130     dY[n] = dD[i];
    131     t[n] = (T[i] - Tmean) / 365.25;
    132     ParFactor (&pX[n], &pY[n], R[i], D[i], T[i]);
    133     pXmin = MIN (pXmin, pX[n]);
    134     pXmax = MAX (pXmax, pX[n]);
    135     pYmin = MIN (pYmin, pY[n]);
    136     pYmax = MAX (pYmax, pY[n]);
    137     n++;
    138   }
    139   float dXRange = pXmax - pXmin;
    140   float dYRange = pYmax - pYmin;
    141   float parRange = hypot (dXRange, dYRange);
    142        
    143   // fprintf (stderr, "par factor range: %f\n", parRange);
    144 
    145   PlxFit fit;
    146   if (!FitPMandPar (&fit, X, dX, Y, dY, t, pX, pY, n, VERBOSE)) {
     144  PlxFitData fitdata;
     145  PlxFitDataAlloc (&fitdata, Ntotal);
     146  PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
     147
     148  PlxFit fit; memset (&fit, 0, sizeof(PlxFit));
     149
     150  // determine dPsig for detections based on Noutlier attempts
     151  if (Noutlier) {
     152    int clipRetry = TRUE;
     153    for (i = 0; clipRetry && (i < 3); i++) {
     154      clipRetry = !PlxOutlierClip (&fitdata, mask, Noutlier, dPsigMax, dPvec, VERBOSE);
     155
     156      // using the new mask values, reset fitdata
     157      PlxSetMeanEpoch (R, D, T, &Rmean, &Dmean, &Tmean, mask, Ntotal);
     158      PlxSetEpochPosition (&fitdata, R, D, dR, dD, T, mask, Ntotal, &coords, Tmean);
     159      if (VERBOSE) fprintf (stderr, "keep %d of %d\n", fitdata.Npts, Ntotal);
     160    }
     161  }
     162
     163  for (i = 0; (VERBOSE == 2) && (i < fitdata.Npts); i++) {
     164    int n = fitdata.index[i];
     165    int maskValue = mask ? mask[n] : 0;
     166    fprintf (stderr, "%f %f : %f %d : %f %f %f\n", R[n], D[n], T[n], maskValue, fitdata.t[i], fitdata.X[i], fitdata.Y[i]);
     167  }
     168
     169  fit.getChisq = TRUE;
     170  if (!FitPMandPar (&fit,
     171                    fitdata.X, fitdata.dX,
     172                    fitdata.Y, fitdata.dY,
     173                    fitdata.t, fitdata.pX, fitdata.pY, fitdata.Npts, VERBOSE)) {
    147174    return FALSE;
    148175  }
     176
     177  if (Nresample){
     178    PlxFitData sample;
     179    PlxFitDataAlloc (&sample, fitdata.Npts);
     180
     181    PlxFit *testfit = NULL;
     182    ALLOCATE (testfit, PlxFit, Nresample);
     183
     184    int Ngood = 0;
     185    for (i = 0; i < Nresample; i++) {
     186      PlxBootstrapResample (&fitdata, &sample);
     187     
     188      if (i % 100000 == 99999) fprintf (stderr, ".");
     189
     190      // fit the sample
     191      testfit[Ngood].getChisq = FALSE;
     192      if (!FitPMandPar (&testfit[Ngood],
     193                        sample.X, sample.dX,
     194                        sample.Y, sample.dY, sample.t,
     195                        sample.pX, sample.pY, sample.Npts, VERBOSE)) continue;
     196      Ngood ++;
     197    }
     198
     199    Vector *pvec, *uRvec, *uDvec, *Rvec, *Dvec;
     200
     201    // save the Nresample histograms
     202    if ((pvec  = SelectVector ("plxVector", ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "plxVector");
     203    if ((uRvec = SelectVector ("uRVector",  ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "uDVector");
     204    if ((uDvec = SelectVector ("uDVector",  ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "uRVector");
     205    if ((Rvec  = SelectVector ("RoVector",  ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "RoVector");
     206    if ((Dvec  = SelectVector ("DoVector",  ANYVECTOR, TRUE)) == NULL) ESCAPE ("missing vector %s\n", "DoVector");
     207   
     208    ResetVector ( pvec, OPIHI_FLT, Ngood);
     209    ResetVector (uRvec, OPIHI_FLT, Ngood);
     210    ResetVector (uDvec, OPIHI_FLT, Ngood);
     211    ResetVector ( Rvec, OPIHI_FLT, Ngood);
     212    ResetVector ( Dvec, OPIHI_FLT, Ngood);
     213   
     214    for (i = 0; i < Ngood; i++) {
     215      pvec->elements.Flt[i]  = testfit[i].p;
     216      uRvec->elements.Flt[i] = testfit[i].uR;
     217      uDvec->elements.Flt[i] = testfit[i].uD;
     218      Rvec->elements.Flt[i]  = testfit[i].Ro;
     219      Dvec->elements.Flt[i]  = testfit[i].Do;
     220    }
     221
     222    // now calculate median and sigma for each vector
     223    VectorRobustStats (pvec,  &fit.p,  &fit.dp);
     224    VectorRobustStats (uRvec, &fit.uR, &fit.duR);
     225    VectorRobustStats (uDvec, &fit.uD, &fit.duD);
     226    VectorRobustStats (Rvec,  &fit.Ro, &fit.dRo);
     227    VectorRobustStats (Dvec,  &fit.Do, &fit.dDo);
     228  }
     229
     230  // fprintf (stderr, "%f +/- %f | %f %f\n", fit.p, fit.dp, fit.uR, fit.uD);
     231
     232/*
     233  FILE *f = fopen ("test.pf.dat", "w");
     234  for (i = 0; i < Ntotal; i++) {
     235    double Xf = fit.Ro + fit.uR*fitdata.t[i] + fit.p*fitdata.pX[i];
     236    double Yf = fit.Do + fit.uD*fitdata.t[i] + fit.p*fitdata.pY[i];
     237    fprintf (f, "%f : %f %f : %f %f : %f : %f %f : %f %f\n", T[i], R[i], D[i], Xf, Yf, fitdata.t[i], fitdata.X[i], fitdata.Y[i], fitdata.pX[i], fitdata.pY[i]);
     238  }
     239  fclose (f);
     240*/
    149241
    150242  // fprintf (stderr, "Roff, Doff: %f, %f; dRo, dDo: %f, %f\n", fit.Ro, fit.Do, fit.dRo, fit.dDo);
    … …  
    152244  XY_to_RD (&Rmean, &Dmean, fit.Ro, fit.Do, &coords);
    153245  if (VERBOSE) {
    154     fprintf (stderr, "Ro, Do: %f, %f +/- %f, %f\n", Rmean, Dmean, fit.dRo, fit.dDo);
     246    fprintf (stderr, "Ro, Do: %f, %f +/- %f, %f (%f, %f)\n", Rmean, Dmean, fit.dRo, fit.dDo, fit.Ro, fit.Do);
    155247    fprintf (stderr, "uR, uD: %f, %f; duR, duD: %f, %f\n", fit.uR, fit.uD, fit.duR, fit.duD);
    156248    fprintf (stderr, "par: %f +/- %f\n", fit.p, fit.dp);
    … …  
    164256  set_variable ("uR",   fit.uR);
    165257  set_variable ("uD",   fit.uD);
    166   set_variable ("duR",   fit.duR);
    167   set_variable ("duD",   fit.duD);
     258  set_variable ("duR",  fit.duR);
     259  set_variable ("duD",  fit.duD);
    168260  set_variable ("plx",  fit.p);
    169261  set_variable ("dplx", fit.dp);
    170262 
    171263  set_variable ("Tmean",  Tmean);
    172   set_variable ("Trange", Trange);
    173   set_variable ("Prange", parRange);
    174264
    175265  set_variable ("chisq", fit.chisq);
    … …  
    293383  fit[0].dp  = sqrt(A[4][4]);
    294384 
    295   // add up the chi square for the fit
    296   chisq = 0.0;
    297   for (i = 0; i < Npts; i++) {
    298     Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].p*pR[i];
    299     Yf = fit[0].Do + fit[0].uD*T[i] + fit[0].p*pD[i];
    300     wx = (fabs(dX[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dX[i]);
    301     wy = (fabs(dY[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dY[i]);
    302     chisq += SQ(X[i] - Xf) * wx;
    303     chisq += SQ(Y[i] - Yf) * wy;
    304     // if (VERBOSE) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq);
    305   }
     385  // (optionally) add up the chi square for the fit
     386  if (fit->getChisq) {
     387    chisq = 0.0;
     388    for (i = 0; i < Npts; i++) {
     389      Xf = fit[0].Ro + fit[0].uR*T[i] + fit[0].p*pR[i];
     390      Yf = fit[0].Do + fit[0].uD*T[i] + fit[0].p*pD[i];
     391      wx = (fabs(dX[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dX[i]);
     392      wy = (fabs(dY[i]) < 0.0001) ? 1.0 : 1.0 / SQ(dY[i]);
     393      chisq += SQ(X[i] - Xf) * wx;
     394      chisq += SQ(Y[i] - Yf) * wy;
     395      // if (VERBOSE) fprintf (stderr, "chisq contrib : %f %f : %f %f : %f %f : %f %f : %f\n", Xf, Yf, X[i] - Xf, Y[i] - Yf, dX[i], dY[i], (X[i] - Xf) / dX[i], (Y[i] - Yf) / dY[i], chisq);
     396    }
     397    // the reduced chisq is divided by (Ndof = 2*Npts - 5)
     398    fit[0].chisq = chisq / (2.0*Npts - 5.0);
     399  }
     400 
    306401  fit[0].Nfit = Npts;
    307 
    308   // the reduced chisq is divided by (Ndof = 2*Npts - 5)
    309   fit[0].chisq = chisq / (2.0*Npts - 5.0);
    310402  return (TRUE);
    311403}
    … …  
    349441  return TRUE;
    350442}
     443
     444// allocate arrays but not the container
     445int PlxFitDataAlloc (PlxFitData *data, int N) {
     446
     447  data->Npts = N;
     448  ALLOCATE (data->X, double, N);
     449  ALLOCATE (data->Y, double, N);
     450  ALLOCATE (data->dX, double, N);
     451  ALLOCATE (data->dY, double, N);
     452  ALLOCATE (data->t, double, N);
     453  ALLOCATE (data->pX, double, N);
     454  ALLOCATE (data->pY, double, N);
     455  ALLOCATE (data->index, int, N);
     456  return TRUE;
     457}
     458
     459void PlxFitDataFree (PlxFitData *data) {
     460  FREE (data->X);
     461  FREE (data->Y);
     462  FREE (data->dX);
     463  FREE (data->dY);
     464  FREE (data->t);
     465  FREE (data->pX);
     466  FREE (data->pY);
     467  FREE (data->index);
     468}
     469
     470int PlxBootstrapResample (PlxFitData *src, PlxFitData *tgt) {
     471  int i;
     472  tgt->Npts = src->Npts;
     473  for (i = 0; i < src->Npts; i++) {
     474    int N = tgt->Npts * drand48();
     475    // int N = i;
     476    tgt->X [i] = src->X [N];
     477    tgt->Y [i] = src->Y [N];
     478    tgt->dX[i] = src->dX[N];
     479    tgt->dY[i] = src->dY[N];
     480    tgt->t [i] = src->t [N];
     481    tgt->pX[i] = src->pX[N];
     482    tgt->pY[i] = src->pY[N];
     483  }
     484  return TRUE;
     485}
     486
     487int PlxSetMeanEpoch (double *R, double *D, double *T, double *Rmean, double *Dmean, double *Tmean, int *mask, int Ntotal) {
     488
     489  int i;
     490
     491  // find mean values to remove
     492  double Nmean = 0;
     493  *Tmean = 0;
     494  *Rmean = 0;
     495  *Dmean = 0;
     496  double Tmin = +1000000;
     497  double Tmax = -1000000;
     498  for (i = 0; i < Ntotal; i++) {
     499    if (mask && !mask[i]) continue;
     500    *Rmean += R[i];
     501    *Dmean += D[i];
     502    *Tmean += T[i];
     503    Tmin = MIN(Tmin, T[i]);
     504    Tmax = MAX(Tmax, T[i]);
     505    Nmean += 1.0;
     506  }
     507  *Rmean /= Nmean;
     508  *Dmean /= Nmean;
     509  *Tmean /= Nmean;
     510 
     511  double Trange = Tmax - Tmin;
     512
     513  // fprintf (stderr, "R,D : %f,%f, T: %f, Trange: %f, Tmin: %f, Tmax: %f\n", *Rmean, *Dmean, *Tmean, Trange, Tmin, Tmax);
     514
     515  set_variable ("Trange", Trange);
     516  return TRUE;
     517}
     518
     519// generate the fit values (projected X,Y; parallax factors;
     520int PlxSetEpochPosition (PlxFitData *fitdata, double *R, double *D, double *dR, double *dD, double *T, int *mask, int Ntotal, Coords *coords, double Tmean) {
     521
     522  int i;
     523
     524  float pXmin = +2.0;
     525  float pXmax = -2.0;
     526  float pYmin = +2.0;
     527  float pYmax = -2.0;
     528
     529  int Nsubset = 0;
     530  for (i = 0; i < Ntotal; i++) {
     531    if (mask && !mask[i]) continue;
     532    RD_to_XY (&fitdata->X[Nsubset], &fitdata->Y[Nsubset], R[i], D[i], coords);
     533    fitdata->dX[Nsubset] = dR[i];
     534    fitdata->dY[Nsubset] = dD[i];
     535    fitdata->t[Nsubset] = (T[i] - Tmean) / 365.25;
     536    ParFactor (&fitdata->pX[Nsubset], &fitdata->pY[Nsubset], R[i], D[i], T[i]);
     537    pXmin = MIN (pXmin, fitdata->pX[Nsubset]);
     538    pXmax = MAX (pXmax, fitdata->pX[Nsubset]);
     539    pYmin = MIN (pYmin, fitdata->pY[Nsubset]);
     540    pYmax = MAX (pYmax, fitdata->pY[Nsubset]);
     541    fitdata->index[Nsubset] = i;
     542    Nsubset++;
     543  }
     544  fitdata->Npts = Nsubset;
     545  float dXRange = pXmax - pXmin;
     546  float dYRange = pYmax - pYmin;
     547  float parRange = hypot (dXRange, dYRange);
     548
     549  set_variable ("Prange", parRange);
     550  // fprintf (stderr, "par factor range: %f\n", parRange);
     551
     552  return TRUE;
     553}
     554
     555/* Outlier clipping based on bootstrap-resampling tests of the plx path
     556 * generate Noutlier resampled datasets
     557 * fit the Noutlier plx paths
     558 * determine and save the distribution of dXsig and dYsig for each point
     559 * sort the resulting distributions and find dPsig (median point) for each measurement
     560 * find the 90% point of dPsig : if > dPsigMax, only clip the 10% most deviant points
     561 * set the dPvec values if desired
     562 * -- mask is modified, dPvec values are set
     563 * -- fitdata is unchanged
     564 */
     565
     566# define MAX_REJECT 0.1
     567
     568int PlxOutlierClip (PlxFitData *fitdata, int *mask, int Noutlier, float dPsigMax, Vector *dPvec, int VERBOSE) {
     569
     570  int i, n;
     571
     572  PlxFit testfit;
     573  testfit.getChisq = FALSE;
     574
     575  PlxFitData sample;
     576  PlxFitDataAlloc (&sample, fitdata->Npts);
     577
     578  double **dXsig, **dYsig;
     579  ALLOCATE (dXsig, double *, fitdata->Npts);
     580  ALLOCATE (dYsig, double *, fitdata->Npts);
     581  for (i = 0; i < fitdata->Npts; i++) {
     582    ALLOCATE (dXsig[i], double, Noutlier);
     583    ALLOCATE (dYsig[i], double, Noutlier);
     584  }
     585
     586  int Nsamples = 0;
     587  for (n = 0; n < Noutlier; n++) {
     588    // bootstrap resample (fitdata -> sample)
     589    PlxBootstrapResample (fitdata, &sample);
     590     
     591    if (n % 100000 == 99999) fprintf (stderr, ".");
     592
     593    // fit the sample
     594    if (!FitPMandPar (&testfit,
     595                      sample.X, sample.dX,
     596                      sample.Y, sample.dY, sample.t,
     597                      sample.pX, sample.pY, sample.Npts, VERBOSE)) continue;
     598
     599    // fprintf (stderr, "%f +/- %f | %f %f\n", testfit.p, testfit.dp, testfit.uR, testfit.uD);
     600
     601    // find the distances to the path
     602    for (i = 0; i < fitdata->Npts; i++) {
     603      double Xf = testfit.Ro + testfit.uR*fitdata->t[i] + testfit.p*fitdata->pX[i];
     604      double Yf = testfit.Do + testfit.uD*fitdata->t[i] + testfit.p*fitdata->pY[i];
     605      dXsig[i][Nsamples] = fabs(fitdata->X[i] - Xf) / fitdata->dX[i];
     606      dYsig[i][Nsamples] = fabs(fitdata->Y[i] - Yf) / fitdata->dY[i];
     607      // fprintf (stderr, "%f : %f %f : %f %f : %f %f : %f %f %f\n", T[i], Xf, Yf, fitdata->X[i], fitdata->Y[i], fitdata->dX[i], fitdata->dY[i], fitdata->t[i], fitdata->pX[i], fitdata->pY[i]);
     608    }
     609    Nsamples ++;
     610  }
     611
     612  double *dPsig;
     613  ALLOCATE (dPsig, double, fitdata->Npts);
     614   
     615  for (i = 0; i < fitdata->Npts; i++) {
     616    dsort (dXsig[i], Nsamples);
     617    dsort (dYsig[i], Nsamples);
     618
     619    // choose the median values
     620    double dXsigMedian, dYsigMedian;
     621    if (Nsamples % 2) {
     622      int Ncenter = Nsamples / 2;
     623      dXsigMedian = dXsig[i][Ncenter];
     624      dYsigMedian = dYsig[i][Ncenter];
     625    } else {
     626      int Ncenter = Nsamples / 2 - 1;
     627      dXsigMedian = 0.5*(dXsig[i][Ncenter] + dXsig[i][Ncenter + 1]);
     628      dYsigMedian = 0.5*(dYsig[i][Ncenter] + dYsig[i][Ncenter + 1]);
     629    }
     630    // XXX replace with hypotenuse?
     631    dPsig[i] = 0.5*(dXsigMedian + dYsigMedian);
     632    // fprintf (stderr, "%d %10.6f %10.6f %10.6f  %f %f : %f\n", i, R[i], D[i], T[i], dXsig[i][Ncenter], dYsig[i][Ncenter], dPsig[i]);
     633  }
     634
     635  // make a copy of dPsig[] and check if > 10% are > dPsigMax
     636  double *dPsigSort;
     637  ALLOCATE (dPsigSort, double, fitdata->Npts);
     638  for (i = 0; i < fitdata->Npts; i++) {
     639    dPsigSort[i] = dPsig[i];
     640  }
     641  dsort (dPsigSort, fitdata->Npts);
     642  int Nmax = (1.0 - MAX_REJECT)*fitdata->Npts;
     643
     644  int completeClip = TRUE;
     645  if (dPsigSort[Nmax] > dPsigMax) {
     646    if (VERBOSE) fprintf (stderr, "too many outliers: %f at 90\n", dPsigSort[Nmax]);
     647    dPsigMax = dPsigSort[Nmax];
     648    completeClip = FALSE;
     649  }
     650
     651  for (i = 0; i < fitdata->Npts; i++) {
     652    if (dPsig[i] < dPsigMax) continue;
     653    int n = fitdata->index[i];
     654    // fprintf (stderr, "clip %d: %f : %f\n", i, fitdata->t[i], dPsig[i]);
     655    mask[n] = 0; // mask these points
     656  }
     657
     658  // only set dPvec if we have completed the clipping?
     659  if (dPvec) {
     660    for (i = 0; i < dPvec->Nelements; i++) {
     661      dPvec->elements.Flt[i] = NAN;
     662    }
     663    for (i = 0; i < fitdata->Npts; i++) {
     664      int n = fitdata->index[i];
     665      dPvec->elements.Flt[n] = dPsig[i];
     666    }
     667  }
     668
     669  free (dPsig);
     670  free (dPsigSort);
     671 
     672  for (i = 0; i < fitdata->Npts; i++) {
     673    free (dXsig[i]);
     674    free (dYsig[i]);
     675  }
     676  free (dXsig);
     677  free (dYsig);
     678
     679  return completeClip;
     680}
     681
     682int VectorRobustStats (Vector *vector, double *median, double *sigma) {
     683
     684  // warn if vector->Nelements > 1000? 10000?)
     685  // warn if vector is not float
     686
     687  // we need to copy the vector to avoid changing the sort order
     688  double *values = NULL;
     689  ALLOCATE (values, double, vector->Nelements);
     690
     691  int i;
     692  int Npts = 0;
     693  for (i = 0; i < vector->Nelements; i++) {
     694    if (!isfinite(vector->elements.Flt[i])) continue;
     695    values[Npts] = vector->elements.Flt[i];
     696    Npts++;
     697  }
     698
     699  dsort (values, Npts);
     700
     701  if (Npts % 2) {
     702    int Ncenter = Npts / 2;
     703    *median = values[Ncenter];
     704  } else {
     705    int Ncenter = Npts / 2 - 1;
     706    *median = 0.5*(values[Ncenter] + values[Ncenter + 1]);
     707  }
     708
     709  double Slo = VectorFractionInterpolate (values, 0.158655, Npts);
     710  double Shi = VectorFractionInterpolate (values, 0.841345, Npts);
     711
     712  *sigma = (Shi - Slo) / 2.0;
     713
     714  return TRUE;
     715}
     716
     717double VectorFractionInterpolate (double *values, float fraction, int Npts) {
     718
     719  float F = fraction * Npts;
     720  int   N = fraction * Npts;
     721
     722  if (N < 0        ) return NAN;
     723  if (N >= Npts - 2) return NAN;
     724
     725  // interpolate between N,N+1
     726   
     727  double S = (F - N) * (values[N+1] - values[N]) + values[N];
     728  return S;
     729}

  • trunk/Ohana/src/opihi/cmd.basic/date.c

    r27435 r38986  
    33int date (int argc, char **argv) {
    44 
    5   int N, SECONDS, REFTIME;
     5  int N, SECONDS;
     6  double REFTIME;
    67  struct timeval now;
    78  char *tstring = NULL;
    … …  
    1617  }
    1718
    18   REFTIME = 0;
     19  REFTIME = 0.0;
    1920  if ((N = get_argument (argc, argv, "-reftime"))) {
    2021    remove_argument (N, &argc, argv);
    21     REFTIME = atoi (argv[N]);
     22    REFTIME = atof (argv[N]);
    2223    remove_argument (N, &argc, argv);
    2324  }
    … …  
    3738  gettimeofday (&now, NULL);
    3839  if (SECONDS) {
     40    double nowSec = now.tv_sec + 1e-6*now.tv_usec;
    3941    if (varName) {
    40       set_int_variable (varName, now.tv_sec - REFTIME);
     42      set_variable (varName, nowSec - REFTIME);
    4143    } else {
    42       gprint (GP_ERR, "%d\n", (int) now.tv_sec - REFTIME);
     44      gprint (GP_ERR, "%.12g\n", nowSec - REFTIME);
    4345    }
    4446  } else {

  • trunk/Ohana/src/opihi/cmd.data/box.c

    r34749 r38986  
    4848    if (strlen (graphmode.ticks) != 4) { goto usage; }
    4949    for (i = 0; i < strlen (graphmode.ticks); i++) {
    50       if ((graphmode.ticks[i] != '0') && (graphmode.ticks[i] != '1') && (graphmode.ticks[i] != '2')) { goto usage; }
     50      if ((graphmode.ticks[i] != '0') && (graphmode.ticks[i] != '1') && (graphmode.ticks[i] != '2') && (graphmode.ticks[i] != '3')) { goto usage; }
    5151    }
    5252  }
    … …  
    139139    remove_argument (N, &argc, argv);
    140140    graphmode.padYp = atof(argv[N]);
     141    remove_argument (N, &argc, argv);
     142  }
     143
     144  if ((N = get_argument (argc, argv, "-fminor"))) {
     145    remove_argument (N, &argc, argv);
     146    graphmode.fMinorXm = atof(argv[N]);
     147    graphmode.fMinorXp = atof(argv[N]);
     148    graphmode.fMinorYm = atof(argv[N]);
     149    graphmode.fMinorYp = atof(argv[N]);
     150    remove_argument (N, &argc, argv);
     151  }
     152  if ((N = get_argument (argc, argv, "-xfminor"))) {
     153    remove_argument (N, &argc, argv);
     154    graphmode.fMinorXm = atof(argv[N]);
     155    remove_argument (N, &argc, argv);
     156  }
     157  if ((N = get_argument (argc, argv, "+xfminor"))) {
     158    remove_argument (N, &argc, argv);
     159    graphmode.fMinorXp = atof(argv[N]);
     160    remove_argument (N, &argc, argv);
     161  }
     162  if ((N = get_argument (argc, argv, "-yfminor"))) {
     163    remove_argument (N, &argc, argv);
     164    graphmode.fMinorYm = atof(argv[N]);
     165    remove_argument (N, &argc, argv);
     166  }
     167  if ((N = get_argument (argc, argv, "+yfminor"))) {
     168    remove_argument (N, &argc, argv);
     169    graphmode.fMinorYp = atof(argv[N]);
     170    remove_argument (N, &argc, argv);
     171  }
     172
     173  if ((N = get_argument (argc, argv, "-flabel"))) {
     174    remove_argument (N, &argc, argv);
     175    graphmode.fLabelRangeXm = atof(argv[N]);
     176    graphmode.fLabelRangeXp = atof(argv[N]);
     177    graphmode.fLabelRangeYm = atof(argv[N]);
     178    graphmode.fLabelRangeYp = atof(argv[N]);
     179    remove_argument (N, &argc, argv);
     180  }
     181  if ((N = get_argument (argc, argv, "-xflabel"))) {
     182    remove_argument (N, &argc, argv);
     183    graphmode.fLabelRangeXm = atof(argv[N]);
     184    remove_argument (N, &argc, argv);
     185  }
     186  if ((N = get_argument (argc, argv, "+xflabel"))) {
     187    remove_argument (N, &argc, argv);
     188    graphmode.fLabelRangeXp = atof(argv[N]);
     189    remove_argument (N, &argc, argv);
     190  }
     191  if ((N = get_argument (argc, argv, "-yflabel"))) {
     192    remove_argument (N, &argc, argv);
     193    graphmode.fLabelRangeYm = atof(argv[N]);
     194    remove_argument (N, &argc, argv);
     195  }
     196  if ((N = get_argument (argc, argv, "+yflabel"))) {
     197    remove_argument (N, &argc, argv);
     198    graphmode.fLabelRangeYp = atof(argv[N]);
    141199    remove_argument (N, &argc, argv);
    142200  }
    … …  
    174232  gprint (GP_ERR, "         alternatively, set each axis independently with:\n");
    175233  gprint (GP_ERR, "        -xpad, -ypad, +xpad, +ypad\n");
     234  gprint (GP_ERR, "  \n");
     235  gprint (GP_ERR, "  -fminor : set the number of minor ticks per major (all axes)\n");
     236  gprint (GP_ERR, "         alternatively, set each axis independently with:\n");
     237  gprint (GP_ERR, "        -xfminor, -yfminor, +xfminor, +yfminor\n");
     238  gprint (GP_ERR, "  \n");
     239  gprint (GP_ERR, "  -flabel : set the fraction of axis over which major ticks have labels (all axes)\n");
     240  gprint (GP_ERR, "         alternatively, set each axis independently with:\n");
     241  gprint (GP_ERR, "        -xflabel, -yflabel, +xflabel, +yflabel\n");
    176242
    177243  return (FALSE);

  • trunk/Ohana/src/opihi/cmd.data/device.c

    r14590 r38986  
    66  char *name;;
    77  /* set / get current graphics device */
     8
     9  int QUIET = FALSE;
     10  if ((N = get_argument (argc, argv, "-q"))) {
     11    remove_argument (N, &argc, argv);
     12    QUIET = TRUE;
     13  }
    814
    915  name = NULL;
    … …  
    2329    if (!GetGraph (NULL, &kapa, name)) return (FALSE);
    2430  }
    25   gprint (GP_ERR, "kapa %s\n", name);
     31  if (!QUIET) gprint (GP_ERR, "kapa %s\n", name);
    2632
    2733  return (TRUE);

  • trunk/Ohana/src/opihi/cmd.data/read_vectors.c

    r38553 r38986  
    511511    if (!gfits_fread_ftable_data (f, &table, padIfShort)) ESCAPE ("error reading table for extension %d\n", Nextend);
    512512  } else {
    513     if (!gfits_fread_ftable_range (f, padIfShort, &table, start, Nrows)) ESCAPE ("error reading table for extension %d\n", Nextend);
     513    // arg3 (FALSE) : seek to this segment start
     514    if (!gfits_fread_ftable_range (f, padIfShort, FALSE, &table, start, Nrows)) ESCAPE ("error reading table for extension %d\n", Nextend);
    514515  }
    515516

  • trunk/Ohana/src/opihi/cmd.data/uniq.c

    r20936 r38986  
    11# include "data.h"
     2// NOTE: if there are only a few uniq values, the old algorithm is not bad. 
     3// for 10000 uniq values, 30M points take ~20sec in the new algorithm,
     4// 3M points takess 45 sec in the old method.
    25
    36int uniq (int argc, char **argv) {
    47 
    5   int Nnew, i, j, found;
     8  int Nnew, i, N;
    69  Vector *ivec, *ovec;
    710
     11  Vector *cvec = NULL;
     12  if ((N = get_argument (argc, argv, "-c"))) {
     13    remove_argument (N, &argc, argv);
     14    if ((cvec = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) {
     15      gprint (GP_ERR, "invalid vector %s\n", argv[N]);
     16      return FALSE;
     17    }
     18    remove_argument (N, &argc, argv);
     19  }
     20
    821  if (argc != 3) {
    9     gprint (GP_ERR, "USAGE: uniq (in) (out)\n");
     22    gprint (GP_ERR, "USAGE: uniq (in) (out) -c count\n");
    1023    return (FALSE);
    1124  }
    … …  
    1629  /* allocate the maximum possible needed */
    1730  ResetVector (ovec, ivec->type, ivec->Nelements);
     31  if (cvec) {
     32    ResetVector (cvec, OPIHI_INT, ivec->Nelements);
     33  }
    1834
    1935  Nnew = 0;
    2036
    2137  if (ivec->type == OPIHI_FLT) {
    22     opihi_flt *v1 = ivec[0].elements.Flt;
    23     for (i = 0; i < ivec[0].Nelements; i++, v1++) {
    24       opihi_flt *v2 = ovec[0].elements.Flt;
    25       found = FALSE;
    26       for (j = 0; !found && (j < Nnew); j++, v2++) {
    27         if (*v1 == *v2) found = TRUE;
     38    // copy the input data to a temporary array to avoid damaging it with sort
     39    opihi_flt *indata = NULL;
     40    ALLOCATE (indata, opihi_flt, ivec[0].Nelements);
     41    memcpy (indata, ivec->elements.Flt, ivec[0].Nelements*sizeof(opihi_flt));
     42
     43    dsort (indata, ivec->Nelements);
     44
     45    Nnew = 0;
     46    opihi_flt *vtgt = ovec[0].elements.Flt;
     47
     48    opihi_flt *vsrc = indata;
     49
     50    for (i = 0; i < ivec->Nelements; Nnew++) {
     51      vtgt[Nnew] = *vsrc;
     52      int Ndup = 0;
     53      opihi_flt lastValue = *vsrc;
     54      while ((i < ivec->Nelements) && (*vsrc == lastValue)) {
     55        i++;
     56        vsrc ++;
     57        Ndup ++;
    2858      }
    29       if (!found) {
    30         ovec[0].elements.Flt[Nnew] = *v1;
    31         Nnew ++;
     59      if (cvec) {
     60        cvec->elements.Int[Nnew] = Ndup;
    3261      }
    3362    }
     63    free (indata);
    3464  } else {
    35     opihi_int *v1 = ivec[0].elements.Int;
    36     for (i = 0; i < ivec[0].Nelements; i++, v1++) {
    37       opihi_int *v2 = ovec[0].elements.Int;
    38       found = FALSE;
    39       for (j = 0; !found && (j < Nnew); j++, v2++) {
    40         if (*v1 == *v2) found = TRUE;
     65    // copy the input data to a temporary array to avoid damaging it with sort
     66    opihi_int *indata = NULL;
     67    ALLOCATE (indata, opihi_int, ivec[0].Nelements);
     68    memcpy (indata, ivec->elements.Int, ivec[0].Nelements*sizeof(opihi_int));
     69
     70    isort (indata, ivec->Nelements);
     71
     72    Nnew = 0;
     73    opihi_int *vtgt = ovec[0].elements.Int;
     74
     75    opihi_int *vsrc = indata;
     76
     77    for (i = 0; i < ivec->Nelements; Nnew++) {
     78      vtgt[Nnew] = *vsrc;
     79      int Ndup = 0;
     80      opihi_int lastValue = *vsrc;
     81      while ((i < ivec->Nelements) && (*vsrc == lastValue)) {
     82        i++;
     83        vsrc ++;
     84        Ndup ++;
    4185      }
    42       if (!found) {
    43         ovec[0].elements.Int[Nnew] = *v1;
    44         Nnew ++;
     86      if (cvec) {
     87        cvec->elements.Int[Nnew] = Ndup;
    4588      }
    4689    }
     90    free (indata);
    4791  }
    4892
    4993  // free up extra memory
    5094  ResetVector (ovec, ivec->type, Nnew);
     95  if (cvec) ResetVector (cvec, OPIHI_INT, Nnew);
    5196
    5297  return (TRUE);

  • trunk/Ohana/src/opihi/dvo/mextract.c

    r38471 r38986  
    213213  }
    214214
    215   // int needLensing = dbFieldNeedLensing (fields, Nfields);
     215  int needLensing = dbFieldNeedLensing (fields, Nfields);
    216216  int needStarpar = dbFieldNeedStarpar (fields, Nfields, FALSE);
    217217
    … …  
    234234    catalog.filename = HOST_ID ? hostfile : skylist[0].filename[i];
    235235    catalog.catflags = DVO_LOAD_AVERAGE | DVO_LOAD_MEASURE | DVO_LOAD_SECFILT;
    236     // catalog.catflags |= needLensing ? DVO_LOAD_LENSING : DVO_SKIP_LENSING;
     236    catalog.catflags |= needLensing ? DVO_LOAD_LENSING : DVO_SKIP_LENSING;
    237237    catalog.catflags |= needStarpar ? DVO_LOAD_STARPAR : DVO_SKIP_STARPAR;
    238238    catalog.Nsecfilt = Nsecfilt;
    … …  
    274274        StarPar *starpar = needStarpar ? &catalog.starpar[Nstarpar] : NULL;
    275275
    276         // int Nlensing = average->lensobjOffset;
    277         // Lensobj *lensobj = needLensobj ? &catalog.lensobj[m] : NULL;
     276        int Nlensing = average->lensingOffset;
     277        Lensing *lensing = needLensing ? &catalog.lensing[Nlensing] : NULL;
    278278
    279279        int Nsec = j*Nsecfilt;
    … …  
    281281
    282282        for (n = 0; n < Nfields; n++) {
    283           values[n] = dbExtractMeasures (average, secfilt, &catalog.measure[m], NULL, starpar, &fields[n]);
     283          values[n] = dbExtractMeasures (average, secfilt, &catalog.measure[m], lensing, starpar, &fields[n]);
    284284        }
    285285        // fprintf (stderr, "object: ave: %f, cat: %f, averef %d\n", fields[n].name, values[2], values[3], catalog.measure[m].averef);

  • trunk/Ohana/src/opihi/lib.data/gaussian.c

    r3894 r38986  
    2020 
    2121  int i;
    22   long A, B;
    2322  double val, x, dx, dx1, dx2, dx3, df;
    2423  double mean, sigma;
    … …  
    2726  if (Ngaussint == Nbin) return;
    2827
    29   A = time(NULL);
    30   for (B = 0; A == time(NULL); B++);
    31   srand48(B);
     28  long A = time(NULL);
     29  srand48(A);
    3230 
    3331  Ngaussint = Nbin;

  • trunk/Ohana/src/opihi/lib.shell/ConfigInit.c

    r16900 r38986  
    3838  }
    3939
     40  // XXX this is a bit dangerous : answer may have an arbitrary
     41  // length, but we do not know the available space in 'ptr'
     42  // we really should be passing in the buffer size and limiting the copy
     43
     44
    4045  if (!strcmp (mode, "%s"))  strcpy ((char *) ptr, answer);
    4146  if (!strcmp (mode, "%d"))  *(int *) ptr       = atoi (answer);

  • trunk/Ohana/src/opihi/lib.shell/SocketOps.c

    r32632 r38986  
    2323  status = gethostname (myHostname, HOST_NAME_MAX);
    2424
    25   fprintf (stderr, "target host: %s, real host: %s\n", hostname, myHostname);
     25  if (strcmp (hostname, myHostname)) {
     26    fprintf (stderr, "target host: %s, real host: %s\n", hostname, myHostname);
     27    fprintf (stderr, "please run on the correct host\n");
     28    exit (2);
     29  }
    2630
    2731  GetPortRange (&start, &stop, portinfo);
    … …  
    197201
    198202  host = gethostbyname (hostname);
     203  if (!host) {
     204    gprint (GP_ERR, "cannot connect to pantasks server %s\n", hostname);
     205    exit (3);
     206  }
     207
    199208  bzero (hostip, 80);
    200209  for (i = 0; i < host[0].h_length; i++) {

  • trunk/Ohana/src/opihi/lib.shell/gprint.c

    r27592 r38986  
    255255  if (file == NULL) {
    256256    // XXX this is a problem: we are leaving open the old file
    257     fprintf (stderr, "cannot open file %s\n", stream[0].name);
     257    fprintf (stderr, "gprint cannot open file %s\n", stream[0].name);
    258258    free (stream[0].name);
    259259    file = (dest == GP_LOG) ? stdout : stderr;

  • trunk/Ohana/src/opihi/pantasks

  • trunk/Ohana/src/opihi/pantasks/CheckJobs.c

    r27614 r38986  
    11# include "pantasks.h"
     2static int Ncheck = 0;
    23
    34float CheckJobs () {
    … …  
    1213  float time_running, next_timeout;
    1314
    14   // int Ncheck;
    15   // Ncheck = 0;
     15  Ncheck ++;
    1616
    1717  // actual maximum delay is controlled in job_threads.c
    … …  
    2121  /** test all jobs: ready to test?  finished? **/
    2222  while ((job = NextJob ()) != NULL) {
    23     // Ncheck ++;
    2423
    2524    task = job[0].task;
    … …  
    220219    SetTaskTimer (&job[0].last);
    221220  }
    222   // fprintf (stderr, "check %d jobs\n", Ncheck);
     221
    223222  JobTaskUnlock();
    224223  return (next_timeout);

  • trunk/Ohana/src/opihi/pantasks/CheckTasks.c

    r31666 r38986  
    11# include "pantasks.h"
     2static int Ncheck = 0;
    23
    34float CheckTasks () {
    … …  
    78  int status;
    89  float time_running, next_timeout, fuzz;
     10
     11  Ncheck ++;
    912
    1013  // actual maximum delay is controlled in job_threads.c

  • trunk/Ohana/src/opihi/pantasks/TaskOps.c

    r36623 r38986  
    2222void FreeTasks () {
    2323  int i;
    24   for (i = 0; i < Ntasks; i++) {
     24  int ntasks = Ntasks;
     25  Ntasks = 0;
     26  for (i = 0; i < ntasks; i++) {
    2527    FreeTask (tasks[i]);
    2628  }

  • trunk/Ohana/src/opihi/pantasks/pantasks_server.c.in

    r32632 r38986  
    1919 
    2020  char hostname[256], portinfo[256];
    21   char log_stdout[1024], log_stderr[1024];
     21  char log_stdout[1024], log_stderr[1024], tmpname[1024];
    2222  pthread_t JobsAndTasksThread;
    2323  pthread_t clientsThread;
    … …  
    3434  stdin = freopen ("/dev/zero", "r", stdin);
    3535
    36   // this block redirects the actual stderr, stdout steams to the output files
    37   if (VarConfig ("PANTASKS_SERVER_STDOUT", "%s", log_stdout) != NULL) {
     36  // this block redirects the actual stderr, stdout streams to the output files
     37  if (VarConfig ("PANTASKS_SERVER_STDOUT", "%s", tmpname) != NULL) {
     38    strcpy (log_stdout, tmpname);
    3839    if (strcmp (log_stdout, "stdout")) {
    3940      stdout = freopen (log_stdout, "a", stdout);
    … …  
    4445    }
    4546  }
    46   if (VarConfig ("PANTASKS_SERVER_STDERR", "%s", log_stderr) != NULL) {
     47  if (VarConfig ("PANTASKS_SERVER_STDERR", "%s", tmpname) != NULL) {
     48    strcpy (log_stderr, tmpname);
    4749    if (strcmp (log_stderr, "stderr")) {
    4850      stderr = freopen (log_stderr, "a", stderr);

  • trunk/Ohana/src/opihi/pantasks/test/local.sh

    r10647 r38986  
    11
    22## a basic test of memory allocation
     3if (not($?VERBOSE)) set VERBOSE = 0
    34
    45exec rm -f tmp.txt
    … …  
    3839 $startmem = $word:1
    3940end
     41
    4042macro memcheck
    4143 list word -x "ps -p $PID -o rss"
Note: See TracChangeset for help on using the changeset viewer.