IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39640


Ignore:
Timestamp:
Jul 28, 2016, 2:14:52 PM (10 years ago)
Author:
eugene
Message:

new api to set vector value en-masse; fix the fit1d_irls mean calculation

Location:
trunk/Ohana/src/opihi
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r39638 r39640  
    22
    33// These should probably be tunable:
    4 # define MAX_ITERATIONS 10
    54# define FIT_TOLERANCE 1e-4
    65# define FLT_TOLERANCE 1e-6
    76# define WEIGHT_THRESHOLD 0.3
    8 # define NBOOTSTRAP 100
    97
    108int fit_least_squares (double *s, double **b, double **c, int nterm, int mterm, opihi_flt *x, opihi_flt *y, opihi_flt *dy, opihi_flt *wt, char *mask, int Npts);
    … …  
    3230  }
    3331
     32  Vector *wtVec = NULL;
     33  if ((N = get_argument (argc, argv, "-wt"))) {
     34    remove_argument (N, &argc, argv);
     35    if (!(wtVec = SelectVector (argv[N], ANYVECTOR, TRUE))) return FALSE;
     36    remove_argument (N, &argc, argv);
     37  }
     38
     39  Vector *yfitVec = NULL;
     40  if ((N = get_argument (argc, argv, "-yfit"))) {
     41    remove_argument (N, &argc, argv);
     42    if (!(yfitVec = SelectVector (argv[N], ANYVECTOR, TRUE))) return FALSE;
     43    remove_argument (N, &argc, argv);
     44  }
     45  Vector *mkVec = NULL;
     46  if ((N = get_argument (argc, argv, "-mask"))) {
     47    remove_argument (N, &argc, argv);
     48    if (!(mkVec = SelectVector (argv[N], ANYVECTOR, TRUE))) return FALSE;
     49    remove_argument (N, &argc, argv);
     50  }
     51
     52  int NBOOTSTRAP = 100;
     53  if ((N = get_argument (argc, argv, "-Nboot"))) {
     54    remove_argument (N, &argc, argv);
     55    NBOOTSTRAP = atoi(argv[N]);
     56    remove_argument (N, &argc, argv);
     57  }
     58  int MaxIterations = 10;
     59  if ((N = get_argument (argc, argv, "-max-iterations"))) {
     60    remove_argument (N, &argc, argv);
     61    MaxIterations = atof(argv[N]);
     62    remove_argument (N, &argc, argv);
     63  }
     64  float WeightThreshold = 0.3;
     65  if ((N = get_argument (argc, argv, "-weight-threshold"))) {
     66    remove_argument (N, &argc, argv);
     67    WeightThreshold = atof(argv[N]);
     68    remove_argument (N, &argc, argv);
     69  }
     70  int UseMean = TRUE;
     71  if ((N = get_argument (argc, argv, "-use-median"))) {
     72    remove_argument (N, &argc, argv);
     73    UseMean = FALSE;
     74  }
     75
    3476  if ((argc != 4) || !dyvec) {
    3577    // require a y error value
    36     gprint (GP_ERR, "USAGE: fit1d x y order -dy wt [-quiet/-q] [-clip Nsigma Niter]\n");
     78    gprint (GP_ERR, "USAGE: fit1d x y order -dy wt [-quiet/-q] [-wt wtOut] [-yfit yfit] [-mask maskOut]\n");
    3779    return (FALSE);
    3880  }
    … …  
    5395  }
    5496
     97  if (mkVec) ResetVector (mkVec, OPIHI_INT, xvec[0].Nelements);
     98
    5599  /* nterm is number of polynomial terms, starting at x^0 */
    56100  order = atof (argv[3]);
    … …  
    87131
    88132  int converged = FALSE;
    89   for (int iterations = 0; !converged && (iterations < MAX_ITERATIONS); iterations++) {
     133  for (int iterations = 0; !converged && (iterations < MaxIterations); iterations++) {
    90134
    91135    for (i = 0; i < xvec[0].Nelements; i++) {
    … …  
    116160  }
    117161     
     162  /* XXX TEST
     163  for (i = 0; i < nterm; i++) {
     164    snprintf (name, 64, "C%d", i);
     165    set_variable (name, b[i][0]);
     166    if (!Quiet) gprint (GP_ERR, "%f x^%d ", b[i][0], i);
     167    snprintf (name, 64, "dC%d", i);
     168    set_variable (name, sqrt(c[i][i]));
     169    if (!Quiet) gprint (GP_ERR, "%f x^%d ", sqrt(c[i][i]), i);
     170  }
     171  if (!Quiet) gprint (GP_ERR, "\n");
     172  return TRUE; */
     173
     174  for (i = 0; i < nterm; i++) {
     175    b_old[i][0] = b[i][0];
     176  }
     177
     178
    118179  // calculate the weight thresholds to mask the bad points:
    119180  double Sum_W = 0.0;
     181  ALLOCATE_PTR (wtlist, double, xvec[0].Nelements);
    120182  for (i = 0; i < xvec[0].Nelements; i++) {
    121183    wt[i] = weight_cauchy (yoff[i] / dy[i]);
     184    wtlist[i] = wt[i];
    122185    Sum_W += wt[i];
    123186  }
    124   double WtThreshold = WEIGHT_THRESHOLD * Sum_W / (1.0 * xvec[0].Nelements);
     187  dsort (wtlist, xvec[0].Nelements);
     188  int midpt = 0.5 * xvec[0].Nelements;
     189  double WtMedian = (xvec[0].Nelements % 2) ? wtlist[midpt] : 0.5*(wtlist[midpt] + wtlist[midpt-1]);
     190  double WtThreshold = UseMean ? WeightThreshold * Sum_W / (1.0 * xvec[0].Nelements) : WeightThreshold * WtMedian;
     191  free (wtlist);
    125192
    126193  // generate the unmasked subset
    … …  
    132199  int Npoints = 0;
    133200  for (i = 0; i < xvec[0].Nelements; i++) {
    134     if (wt[i] < WtThreshold) continue;
     201    if (mkVec) mkVec->elements.Int[i] = 0;
     202    if (wt[i] < WtThreshold) {
     203      if (mkVec) mkVec->elements.Int[i] = 1;
     204      continue;
     205    }
    135206     xunmask[Npoints] =  xvec[0].elements.Flt[i];
    136207     yunmask[Npoints] =  yvec[0].elements.Flt[i];
    … …  
    184255  for (i = 0; i < nterm; i++) {
    185256    snprintf (name, 64, "C%d", i);
    186     set_variable (name, b[i][0]);
    187     if (!Quiet) gprint (GP_ERR, "%f x^%d ", b[i][0], i);
     257    set_variable (name, b_old[i][0]);
     258    if (!Quiet) gprint (GP_ERR, "%f x^%d ", b_old[i][0], i);
    188259  }
    189260  if (!Quiet) gprint (GP_ERR, "\n");
    … …  
    216287  free (s);
    217288  free (yoff);
    218   free (yfit);
    219   free (wt);
     289
     290  if (wtVec) {
     291    SetVectorValues (wtVec, wt, OPIHI_FLT, xvec[0].Nelements);
     292  } else {
     293    free (wt);
     294  }
     295  if (yfitVec) {
     296    SetVectorValues (yfitVec, yfit, OPIHI_FLT, xvec[0].Nelements);
     297  } else {
     298    free (yfit);
     299  }
    220300 
    221301  free (xunmask);

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

    r39360 r39640  
    150150int           CopyVector            PROTO((Vector *out, Vector *in));
    151151int           ResetVector           PROTO((Vector *vec, char type, int Nelements));
     152int           SetVectorValues       PROTO((Vector *vec, void *data, char type, int Nelements));
    152153int           SetVector             PROTO((Vector *vec, char type, int Nelements));
    153154int           CastVector            PROTO((Vector *vec, char type));

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

    r39225 r39640  
    229229  } else {
    230230    REALLOCATE (vec[0].elements.Int, opihi_int, MAX(1, Nelements));
     231    vec[0].type = OPIHI_INT;
     232  }
     233  return TRUE;
     234}
     235
     236int SetVectorValues (Vector *vec, void *data, char type, int Nelements) {
     237
     238  // a vector can only have >= 0 elements
     239  vec[0].Nelements = MAX(Nelements,0);
     240  if (type == OPIHI_FLT) {
     241    free (vec[0].elements.Flt);
     242    vec[0].elements.Flt = (opihi_flt *) data;
     243    vec[0].type = OPIHI_FLT;
     244  } else {
     245    free (vec[0].elements.Int);
     246    vec[0].elements.Int = (opihi_int *) data;
    231247    vec[0].type = OPIHI_INT;
    232248  }
Note: See TracChangeset for help on using the changeset viewer.