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
