Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I want a linear interpolation from the following example list:

list={{0.0005023, 22.24}, {0.01457, 21.47}, {0.04922, 19.79}, 
      {0.07484, 18.7}, {0.104, 17.55}, {0.1331, 16.52}, {0.1632, 15.49},
      {0.1888, 14.52}, {0.2215, 13.31}, {0.2506, 12.16}, {0.3024, 10.01}, 
      {0.3435, 8.304}, {0.3943, 6.036}, {0.4098, 5.329}, {0.4726, 2.384}};

The easiest way is to use:

Interpolation[list, InterpolationOrder -> 1]

but my list will be changing a lot, and the InterpolatingFunction takes a lot of time to build:

Timing[
   Table[Interpolation[list, InterpolationOrder -> 1][q], {q, 
     0.0006, 0.4, 0.00001}];]

is 10x slower than:

test=Interpolation[list, InterpolationOrder -> 1];
Timing[Table[test[q], {q,0.0006, 0.4, 0.00001}];]

How can I remove the overhead?

EDIT (following JxB comment)

This compiled version is 5 times faster than the original version, but I don't think Partition is compiling (it appears between all the Lists when I use FullForm); and there's also a CopyTensor that doesn't look good:

Compile[{{list, _Real, 2}, {value, _Real, 0}},
 Module[{temp},
  temp = Select[
     Partition[list, 2, 1], #[[1, 1]] <= value && #[[2, 1]] > value &][[1]
   ];
  temp[[1, 2]] +
   (value - temp[[1, 1]])/(temp[[2, 1]] - temp[[1, 1]])*(temp[[2, 2]] - temp[[1, 2]])
  ]
 ]

Any suggestions? (I don't want to compile to C)

share|improve this question
1  
Perhaps Timing[Table[ Evaluate@Interpolation[list, InterpolationOrder -> 1][q], {q, 0.0006, 0.4, 0.00001}];] ? –  belisarius May 2 '12 at 12:46
1  
Since you are using linear interpolation, it might be straightforward to build your own compiled version of an interpolating function. –  JxB May 2 '12 at 12:53
    
Do your x-value remain the same from list to list? –  rcollyer May 2 '12 at 14:50
    
@rcollyer my x-value also changes –  P. Fonseca May 2 '12 at 16:21
    
Opps... from the answers and comments I see that apparently I didn't explain myself correctly. The Table was there just to raise the timing to readable measurements (I should have used Do...). The calls to "test" will be made almost for one q at a time, and in-between list may change. –  P. Fonseca May 2 '12 at 16:30
add comment

4 Answers

up vote 8 down vote accepted

You can use binary search with Compile. I failed inlining (Compile was complaining endlessly about types mismatch), so I included a binary search directly into Compile-d function. The code for binary search itself corresponds to the bsearchMin function from this answer.

Clear[linterp];
linterp =
   Compile[{{lst, _Real, 2}, {pt, _Real}},
     Module[{pos  = -1 , x = lst[[All, 1]], y = lst[[All, 2]], n0 = 1, 
          n1 = Length[lst], m = 0},
      While[n0 <= n1, m = Floor[(n0 + n1)/2];
        If[x[[m]] == pt,
          While[x[[m]] == pt  && m < Length[lst], m++];
          pos = If[m == Length[lst], m, m - 1];
          Break[];
        ];
        If[x[[m]] < pt, n0 = m + 1, n1 = m - 1]
      ];
      If[pos == -1, pos = If[x[[m]] < pt, m, m - 1]];
      Which[
        pos == 0,
           y[[1]],
        pos == Length[x],
           y[[-1]],
        True,
        y[[pos]] + (y[[pos + 1]] - y[[pos]])/(x[[pos + 1]] - 
              x[[pos]])*(pt - x[[pos]])
      ]],
      CompilationTarget -> "C"];

This is about 20 times faster, on my benchamrks:

AbsoluteTiming[
   Table[Interpolation[list,InterpolationOrder->1][q],{q,0.0006,0.4,0.00001}];
]

{1.453,Null}

AbsoluteTiming[
   Table[linterp[list,q],{q,0.0006,0.4,0.00001}];
]

{0.063,Null}

share|improve this answer
    
I just realized the OP specifically said "I don't want to compile to C." What timings do you get on v8 using WVM? –  Mr.Wizard May 2 '12 at 19:18
    
@Mr.Wizard It's about twice slower. The difference is not so dramatic, because the complexity itself is only logarithmic. –  Leonid Shifrin May 2 '12 at 19:48
1  
@P.Fonseca It is, but it is a hassle. You will need to do this for all platforms (since shared libs are platform-dependent), and you will have to dispatch to a right one depending on the platform. Nothing too complicated though, and sounds like a very good project (since this is a general problem). I think, a general module for serialization of compiled functions would be very useful. The problem would be, that to use it on a single machine, one would need all cross-compilers to other platforms installed. Macs seem best equipped for that, but it would still be some work to set that up. –  Leonid Shifrin May 2 '12 at 21:03
1  
@P.Fonseca if you are precompiling, then sure, I can understand only using WVM. However, if you load CCompilerDriver` the function DefaultCCompiler[] will generate the message CreateLibrary::nocomp and return $Failed if no default compiler was set. Personally, I'd Quiet the message and test for $Failed directly, as there is less chance of the message escaping. Using that info, you can choose your compilation target. –  rcollyer May 3 '12 at 1:23
1  
@rcollyer "that may be counter productive with the OPs purpose" - this is exactly the reason why I did not go there. Inlining a list would make sense when the list doesn't change, which is not the case OP is interested in. –  Leonid Shifrin May 3 '12 at 9:47
show 11 more comments

Combinatorica functions are often not well optimized, so there may very well be a faster binary search algorithm. If that can be found, this might be effective:

Needs["Combinatorica`"]

f[{{a_, b_}, {c_, d_}}][x_] := b + (d - b)/(c - a) (x - a)

list[[Floor@{#, # + 1}]] & @ BinarySearch[list[[All, 1]], 0.33]

f[%][0.33]
8.86436

Check:

Interpolation[list, InterpolationOrder -> 1][0.33]
8.86436
share|improve this answer
    
I was fighting with inlining in Compile (and lost), and was beaten to it by you, as a result (so lost again). +1. –  Leonid Shifrin May 2 '12 at 18:57
    
@Leonid if I had taken the time to look up the faster binary search (which I believe you posted before) you'd have been first. +1 to you as well. :-) –  Mr.Wizard May 2 '12 at 19:11
1  
It seems you don't need to load Combinatorica anymore; you can use GeometricFunctions`BinarySearch[] instead. –  J. M. May 10 '13 at 14:03
add comment

Something with memory ?

myTest[alist_] :=  myTest[alist] = Interpolation[alist, InterpolationOrder -> 1]

Timing[Table[myTest[list][q], {q, 0.0006, 0.4, 0.00001}];]

(* {0.187,Null} *)

test=Interpolation[list,InterpolationOrder->1];
Timing[Table[test[q],{q,0.0006,0.4,0.00001}];]

(* {0.172,Null} *)
share|improve this answer
    
This might work. It builds the InterpolatingFunction just once for each new list. I don't think I will have more than a couple thousands lists during a session, and so, I believe still manageable. Do you see a way to purge the memorization (memory or quantity) if it passes over a certain value (obvious, without too much overhead…)? –  P. Fonseca May 2 '12 at 16:37
    
@P.Fonseca Do you mean remove the previous definitions if they exceed a certain memory usage ? If so, sorry I don't. –  b.gatessucks May 2 '12 at 18:38
    
Yes. If the total myTest definitions exceed a certain quantity or memory usage –  P. Fonseca May 2 '12 at 18:41
add comment

Here is a linear interpolation routine that uses binary search with a few refinements (in particular, the binary search is skipped in the case of equispaced abscissas), as well as a stabilized version of the linear interpolation formula:

lerp = Compile[{{dat, _Real, 2}, {x, _Real}},
  Module[{n = Length[dat], k = 1, l, m, r, xa, ya},
         {xa, ya} = Transpose[dat];
         l = 1 + Quotient[(n - 2) (x - First[xa]), Last[xa] - First[xa]];

         If[xa[[l]] <= x,
            r = l + 1;
            While[xa[[r]] < x,
                  l = r; k *= 2; r = Min[l + k, n]],
            {l, r} = {l - 1, l};
            While[x < xa[[l]],
                  r = l; k *= 2; l = Max[1, r - k]]];

         While[r - l > 1,
               m = Quotient[l + r, 2];
               If[x < xa[[m]], r = m, l = m]];

         {ya[[l]], ya[[r]]}.{xa[[r]] - x, x - xa[[l]]}/(xa[[r]] - xa[[l]])],
         RuntimeOptions -> "Speed"]

Even without the compilation to C, the method is quite fast on my box:

AbsoluteTiming[Table[Interpolation[data, InterpolationOrder -> 1][q],
                     {q, 0.0006, 0.4, 0.00001}];][[1]]
   15.206078

AbsoluteTiming[Table[lerp[data, q], {q, 0.0006, 0.4, 0.00001}];][[1]]
   0.693506
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.