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
