IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 42967


Ignore:
Timestamp:
Mar 9, 2026, 11:26:15 AM (2 months ago)
Author:
eugene
Message:

add spline options to extrapolate or saturate at end points; add vweave function to interpolate a vector with varying bin sizes

Location:
trunk/Ohana/src/opihi
Files:
6 edited
1 copied

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/opihi/cmd.data/Makefile

    r42821 r42967  
    185185$(SRC)/vbin.$(ARCH).o              \
    186186$(SRC)/vgroup.$(ARCH).o            \
     187$(SRC)/vweave.$(ARCH).o            \
    187188$(SRC)/vclip.$(ARCH).o             \
    188189$(SRC)/vgauss.$(ARCH).o            \

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

    r42821 r42967  
    170170int vbin             PROTO((int, char **));
    171171int vgroup           PROTO((int, char **));
     172int vweave           PROTO((int, char **));
    172173int vclip            PROTO((int, char **));
    173174int vect_select      PROTO((int, char **));
    … …  
    386387  {1, "vbin",         vbin,             "rebin vector data by a factor of N"},
    387388  {1, "vgroup",       vgroup,           "group y vector into bins defined by x vector values"},
     389  {1, "vweave",       vweave,           "drizzle input vectors to output vector with correct fractional bin weighting"},
    388390  {1, "vclip",        vclip,            "clip values in a vector to be within a range"},
    389391  {1, "vectors",      list_vectors,     "list vectors"},

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

    r42332 r42967  
    6262int spline_apply (int argc, char **argv) {
    6363
    64   int i;
     64  int i, N;
    6565  Vector *xvec, *yvec;
     66
     67  spline_contruct_mode (SPLINE_SATURATE);
     68  if ((N = get_argument (argc, argv, "-extrapolate"))) {
     69    spline_contruct_mode (SPLINE_EXTRAPOLATE);
     70    remove_argument (N, &argc, argv);
     71  }
    6672
    6773  if (argc != 4) {

  • trunk/Ohana/src/opihi/cmd.data/test/spline.sh

    r41341 r42967  
    11
    22macro test1
    3 
    43  for i 0 100
    54    create x 0 50
    … …  
    109  end
    1110end
     11
     12macro mytest.saturate
     13  vlist x 400 500  600  700  800  900
     14  vlist y  15 425  610  700  745  775
     15
     16  spline create t1 x y
     17
     18  create X 300 1000 5
     19  spline apply t1 X Y
     20
     21  lim X Y; clear; box; plot x y; plot -c red -x line X Y
     22end
     23 
     24macro mytest.extrapolate
     25  vlist x 400 500  600  700  800  900
     26  vlist y  15 425  610  700  745  775
     27
     28  spline create t1 x y
     29
     30  create X 300 1000 5
     31  spline apply t1 X Y -extrapolate
     32
     33  lim X Y; clear; box; plot x y; plot -c red -x line X Y
     34end
     35 

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

    r42946 r42967  
    66 
    77  int i, j, N;
    8   Vector *xin, *yin, *xout, *yout;
     8  Vector *xin, *yin, *yout;
     9  Vector *xminV = NULL, *xmaxV = NULL, *xout = NULL;
    910  opihi_flt xmin, xmax, Xmin, Xmax;
    1011
    … …  
    2223  }
    2324
    24   if (argc != 5) {
     25  if ((argc != 5) && (argc != 6)) {
    2526    gprint (GP_ERR, "USAGE: vweave <xin> <yin> <xout> <yout>\n");
     27    gprint (GP_ERR, "   OR: vweave <xin> <yin> <xmin> <xmax> <yout>\n");
    2628    gprint (GP_ERR, " fully interpolate input bin values into output bins, using fractional input bin contributions\n");
    2729    gprint (GP_ERR, " by default, yout has the sum of the associated input values\n");
    2830    gprint (GP_ERR, " use -count to find the fractional input bin contributions\n");
    29     gprint (GP_ERR, " output bin boundaries defined as midpoint of xout values\n");
     31    gprint (GP_ERR, " with a single <xout> vector, output bin boundaries defined as midpoint of xout values\n");
     32    gprint (GP_ERR, " with <xmin> and <xmax>, the output bin boundaries are explicitly provided\n");
    3033    gprint (GP_ERR, " use -binsize to specify a fixed bin width (<xout> will define the bin center)\n");
    3134    gprint (GP_ERR, " input values are treated as uniformly spread between midpoints of xin values\n");
    … …  
    3841    return FALSE;
    3942  }
    40   if ((xout = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
    41   if (xout->Nelements < 2) {
    42     gprint (GP_ERR, "output vector of length < 2 is invalid\n");
    43     return FALSE;
     43  if ((yin  = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
     44
     45  int Nout = -1;
     46  if (argc == 5) {
     47    if ((xout = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
     48    if (xout->Nelements < 2) {
     49      gprint (GP_ERR, "output vector of length < 2 is invalid\n");
     50      return FALSE;
     51    }
     52    Nout = xout->Nelements;
     53    if ((yout = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
     54  } else {
     55    if (!isnan(binsize)) {
     56      gprint (GP_ERR, "-binsize is incompatible with <xmin>, <xmax>\n");
     57      return FALSE;
     58    }
     59    if ((xminV = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
     60    if (xminV->Nelements < 2) {
     61      gprint (GP_ERR, "output vector of length < 2 is invalid\n");
     62      return FALSE;
     63    }
     64    Nout = xminV->Nelements;
     65    if ((xmaxV = SelectVector (argv[4], OLDVECTOR, TRUE)) == NULL) return (FALSE);
     66    if (xmaxV->Nelements != Nout) {
     67      gprint (GP_ERR, "<xmin> and <xmax> lengths must match\n");
     68      return FALSE;
     69    }
     70    if ((yout = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) return (FALSE);
    4471  }
    4572
    46   if ((yin  = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
    47   if ((yout = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
    48 
    4973  // re-binning creates a float vector
    50   ResetVector (yout, OPIHI_FLT, xout[0].Nelements);
     74  ResetVector (yout, OPIHI_FLT, Nout);
    5175
    5276  // XXX this algorithm is not optimized: full scale of input vector for each output bin
    5377  // could speed up with sorting and pre-defined boundaries
    54   for (i = 0; i < xout[0].Nelements; i++) {
     78  for (i = 0; i < Nout; i++) {
    5579    if (!isnan(binsize)) {
    5680      xmin = xout[0].elements.Flt[i] - 0.5*binsize;
    5781      xmax = xout[0].elements.Flt[i] + 0.5*binsize;
    5882    } else {
    59       if (i > 0) {
    60         xmin = 0.5*(xout[0].elements.Flt[i] + xout[0].elements.Flt[i-1]);
     83      if (xout) {
     84        if (i > 0) {
     85          xmin = 0.5*(xout[0].elements.Flt[i] + xout[0].elements.Flt[i-1]);
     86        } else {
     87          xmin = 1.5*xout[0].elements.Flt[i] - 0.5*xout[0].elements.Flt[i+1];
     88        }
     89        if (i < xout[0].Nelements - 1) {
     90          xmax = 0.5*(xout[0].elements.Flt[i] + xout[0].elements.Flt[i+1]);
     91        } else {
     92          xmax = 1.5*xout[0].elements.Flt[i] - 0.5*xout[0].elements.Flt[i-1];
     93        }
    6194      } else {
    62         xmin = 1.5*xout[0].elements.Flt[i] + 0.5*xout[0].elements.Flt[i+1];
    63       }
    64       if (i < xout[0].Nelements - 1) {
    65         xmax = 0.5*(xout[0].elements.Flt[i] + xout[0].elements.Flt[i+1]);
    66       } else {
    67         xmax = 1.5*xout[0].elements.Flt[i] - 0.5*xout[0].elements.Flt[i-1];
     95        xmin = xminV[0].elements.Flt[i];
     96        xmax = xmaxV[0].elements.Flt[i];
    6897      }
    6998    }

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

    r42456 r42967  
    3333  char **pageIDs;
    3434} Book;
     35
     36typedef enum {SPLINE_SATURATE, SPLINE_EXTRAPOLATE} SplineMode;
    3537
    3638// the interpolating spline has valu
    … …  
    138140void spline_construct_dbl (opihi_flt *x, opihi_flt *y, int N, opihi_flt *y2, opihi_flt dyLower, opihi_flt dyUpper);
    139141opihi_flt spline_apply_dbl (opihi_flt *x, opihi_flt *y, opihi_flt *y2, int N, opihi_flt X);
     142int spline_contruct_mode (SplineMode mode);
    140143
    141144/* in svdcmp.c */

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

    r42821 r42967  
    11# include "data.h"
     2
     3static SplineMode MODE = SPLINE_SATURATE; // SATURATE or EXTRAPOLATE?
     4
     5int spline_contruct_mode (SplineMode mode) {
     6 
     7  switch (mode) {
     8  case SPLINE_SATURATE:
     9  case SPLINE_EXTRAPOLATE:
     10    MODE = mode;
     11    break;
     12  default:
     13    return FALSE;
     14  }
     15  return TRUE;
     16}
    217
    318/* construct the natural spline for x, y in y2 */
    … …  
    3651 
    3752  // saturate correction at high and low ends
    38   if (X < x[0]) return y[0];
    39   if (X > x[N-1]) return y[N-1];
     53  if (MODE == SPLINE_SATURATE) {
     54    if (X < x[0]) return y[0];
     55    if (X > x[N-1]) return y[N-1];
     56  }
    4057
    4158  /* find correct element in array (x must be sorted) */
    … …  
    113130  opihi_flt dx, a, b, value;
    114131 
    115   // linear extrapolation past endpoints
    116   if (X < x[0]) {
    117     hi = 1;
    118     lo = 0;
    119     goto evaluate;
    120     // alternative: saturate correction at high and low ends
    121     // return y[0];
    122   }
    123   if (X > x[N-1]) {
    124     hi = N - 1;
    125     lo = N - 2;
    126     goto evaluate;
    127     // alternative: saturate correction at high and low ends
    128     // return y[N-1];
     132  // saturate correction at high and low ends
     133  if (MODE == SPLINE_SATURATE) {
     134    if (X < x[0]) return y[0];
     135    if (X > x[N-1]) return y[N-1];
    129136  }
    130137
    … …  
    140147    }
    141148  }
    142 
    143 evaluate:
    144149
    145150  /* error condition: duplicate abssisca */
Note: See TracChangeset for help on using the changeset viewer.