The line
index=. <1 i:~"1 x>:/__,}._,~}:>{.y
creates a 2-dimensional array, shape (#y),(#y) which
no doubt is what is slow.
What you want is something like
index =. ({.y) I.&:> x
but not exactly, as you can see by experimenting. I think if
you change the final computation you can make the faster
index computation serve.
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
