I have a slightly better version: fewer copies of the data, and
the structure of the code will allow replacing the linear
interpolation with a cubic spline or the like as needed:
NB. Resampling
NB. y is (x values),:(y values)
NB. x is new x values
NB. Result is new y values
resamp =: 4 : 0
NB. Intervals are numbered by the index of the sample that
NB. ends the interval. So, interval 0 is before the first sample
NB. and interval {:$y is after the last. We calculate the
NB. interval number for each x and then, if it is one of those
NB. off-the-end intervals, adjust to the nearest interior interval.
NB. This means we extrapolate out-of-range values using the slope
NB. of the outermost intervals.
ix =. 1 >. (<:{:$y) <. (0{y) I. x
NB. Calculate the interpolating polynomial for each interval.
NB. Here we use linear interpolation, so the polynomial is (y value),(dy/dx)
NB. Create a polynomial for the first interval (off the beginning),
NB. using the slope of the first internal interval
intpoly =. (1 { y) ,. (,~ {.) %~/ 2 -/\"1 y
NB. The value to return is the interpolating polynomial, evaluated
NB. given the distance between the desired value and the origin point
NB. (i. e. right endpoint) of the interval
(ix { intpoly) p. ((<0;ix) { y) -~ x
)
> -----Original Message-----
> From: [EMAIL PROTECTED]
> [mailto:[EMAIL PROTECTED] On Behalf Of Henry Rich
> Sent: Sunday, June 15, 2008 5:12 PM
> To: 'Programming forum'
> Subject: RE: [Jprogramming] fast 1-d linear interpolation
>
> Here is a verb for resampling. Tested a little bit.
> It uses unboxed right operand because I didn't see that
> boxing was necessary & it usually slows things down.
>
> NB. Resampling
> NB. y is (x values),:(y values)
> NB. x is new x values
> NB. Result is new y values
> resamp =: 4 : 0
> NB. Intervals are numbered by the index of the sample that
> NB. ends the interval. So, interval 0 is before the first sample
> NB. and interval {:$y is after the last. We calculate the
> NB. interval number for each x and then, if it is one of those
> NB. off-the-end intervals, adjust to the nearest interior interval.
> NB. This means we extrapolate out-of-range values using the slope
> NB. of the outermost intervals.
> ix =. 1 >. (<:{:$y) <. (0{y) I. x
> NB. Calculate the slope dy/dx for each interval. Prepend a
> NB. dummy value to account for the fact that the first interval
> NB. is actually interval number 1
> intslope =. 0 , %~/ 2 -/\"1 y
> NB. The value to return is the y at the end of the interval,
> NB. adjusted by the distance of the sampling x from the end of the
> NB. interval
> (ix { 1 { y) + (ix { intslope) * (ix { 0 { y) -~ x
> )
>
> Henry Rich
>
> > -----Original Message-----
> > From: [EMAIL PROTECTED]
> > [mailto:[EMAIL PROTECTED] On Behalf Of sean
> > Sent: Thursday, June 12, 2008 7:34 AM
> > To: Programming forum
> > Subject: [Jprogramming] fast 1-d linear interpolation
> >
> > I'd be interested to see if anyone can suggest improvements to the
> > following code snippets...
> >
> > I have been doing some calculations recently involving
> > interpolations on
> > vectors containing several thousand points. A simple
> > equivalent of the
> > problem might involve a vector of X values and a
> > corresponding vector of
> > Y values: for example,
> >
> > X =: i. 10000
> > Y =: X^1.1
> >
> > I then need to scale the x-values by a small amount and
> interpolate on
> > some of these interior points XI. For example,
> >
> > XI =: 1.1 * i.8000
> >
> > So I want to find YI that corresponds to linearly
> > interpolated values of
> > Y evaluated at XI.
> >
> > I initially used the algorithm posted to this forum
> previously by R.E.
> > Boss, and shown below.
> >
> > interp1=: 4 : 0 NB. From R.E. Boss,
> > NB. eg XI interp1 X;Y, where XI are the interpolation points
> > NB. and X and Y are the original data
> > slopes=. %~/2-~/\&>y
> > index=. <1 i:~"1 x>:/__,}._,~}:>{.y
> > 't0 t1 t2'=. index {&> slopes;y
> > t2+t0*x-t1
> > )
> >
> > This is very elegant and accounts for extrapolation. But on
> > the example
> > matrix here it takes about 5 seconds to run on my computer. I would
> > like it to be faster than this. I will have to call this routine
> > several times as it's part of a least-squares fitting routine.
> >
> > I can arrange XI before interpolating so that it is always
> within the
> > limits of X, removing the need to cater for extrapolation.
> > Also, both X
> > and XI have monotonically increasing values. Looking at
> interp1, most
> > of the time is spent evaluating the x >:/ part of the code.
> >
> > To speed it up, I noted that once the index of X immediately
> > lower than
> > the current XI had been determined, searches of X for subsequent
> > elements of XI don't need to keep checking the elements of
> X prior to
> > that index. Using loops, I wrote the following slight
> modification of
> > the code to avoid these extra checks:
> >
> > indices =: 4 : 0
> > list =. 0
> > for_XI. x do.
> > for_X. ({:list)}.y do.
> > if. (X > XI) do. list=.list,(({: list) + <: X_index) break. end.
> > end.
> > end.
> > < }. list
> > )
> >
> > interp2 =: 4 : 0
> > slopes=. %~/2-~/\&>y
> > 't0 t1 t2'=: (x indices 0{ > y) {&> slopes;y
> > t2+t0*x-t1
> > )
> >
> > This works OK for my data, and the execution time (around 0.3
> > seconds on
> > my computer for the data above) is acceptable for my needs.
> >
> > So is there a way of doing something similar to the indices
> verb that
> > does not involve using for. loops? Or, indeed, is there a
> > neater way of
> > doing this using for. loops? I'd like to see how others
> > might tackle this.
> >
> > Thanks,
> >
> > Sean O'Byrne
> >
> ----------------------------------------------------------------------
> > For information about J forums see
> > http://www.jsoftware.com/forums.htm
>
> ----------------------------------------------------------------------
> For information about J forums see
> http://www.jsoftware.com/forums.htm
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
