This is what I've used for linear interpolation in large arrays. I
tried it on the example below and it took ~0.002 sec on my machine.
---------------------------------------------------------------------
NB. linear interpolation: (x;y) lintrp x0 --> y(x0)
NB. x & y boxed, x0 may be a vector, extrapolate off ends.
lintrp=: 4 : 0
'xx yy'=. x
g0=. (}: yy)% d=. DEL xx
g1=. (}. yy)% d
u=. (g0* }.xx) - g1* }:xx
n=. xx IdX y
(n{u)+ y*n{g1-g0
)
NB. x IdX x0 --> index where x0 fits in x-array, e.g.:
NB. 0 1 2 3 4 5 6 <- x values
NB. | | | | | | |
NB. <-- 0 | 1 | 2 | 3 | 4 | 5 ---> IdX values
IdX=: ] I.~ [: }: [: }. [
DEL =: }. - }:
------------------------------------------------------------------------
Patrick
PS. I've sent this to Sean, but had trouble posting to the newsgroup
last week as my "From:" name had changed.
On Thu, 12 Jun 2008, sean wrote:
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
