24
$\begingroup$

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 10× 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.)

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ Perhaps Timing[Table[ Evaluate@Interpolation[list, InterpolationOrder -> 1][q], {q, 0.0006, 0.4, 0.00001}];] ? $\endgroup$
    – Dr. belisarius
    May 2, 2012 at 12:46
  • 1
    $\begingroup$ Since you are using linear interpolation, it might be straightforward to build your own compiled version of an interpolating function. $\endgroup$
    – JxB
    May 2, 2012 at 12:53
  • $\begingroup$ Do your x-value remain the same from list to list? $\endgroup$
    – rcollyer
    May 2, 2012 at 14:50
  • $\begingroup$ @rcollyer my x-value also changes $\endgroup$
    – P. Fonseca
    May 2, 2012 at 16:21
  • $\begingroup$ 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. $\endgroup$
    – P. Fonseca
    May 2, 2012 at 16:30

4 Answers 4

Reset to default
21
$\begingroup$

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}

Improve this answer
$\endgroup$
18
  • $\begingroup$ I just realized the OP specifically said "I don't want to compile to C." What timings do you get on v8 using WVM? $\endgroup$
    – Mr.Wizard
    May 2, 2012 at 19:18
  • $\begingroup$ @Mr.Wizard It's about twice slower. The difference is not so dramatic, because the complexity itself is only logarithmic. $\endgroup$
    – Leonid Shifrin
    May 2, 2012 at 19:48
  • 1
    $\begingroup$ @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. $\endgroup$
    – Leonid Shifrin
    May 2, 2012 at 21:03
  • 1
    $\begingroup$ @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. $\endgroup$
    – rcollyer
    May 3, 2012 at 1:23
  • 1
    $\begingroup$ @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. $\endgroup$
    – Leonid Shifrin
    May 3, 2012 at 9:47
7
$\begingroup$

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
Improve this answer
$\endgroup$
3
  • $\begingroup$ I was fighting with inlining in Compile (and lost), and was beaten to it by you, as a result (so lost again). +1. $\endgroup$
    – Leonid Shifrin
    May 2, 2012 at 18:57
  • $\begingroup$ @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. :-) $\endgroup$
    – Mr.Wizard
    May 2, 2012 at 19:11
  • 2
    $\begingroup$ It seems you don't need to load Combinatorica anymore; you can use GeometricFunctions`BinarySearch[] instead. $\endgroup$
    – J. M.'s persistent exhaustion ♦
    May 10, 2013 at 14:03
6
$\begingroup$

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 = Min[Max[2, 1 + Quotient[x - First[xa],
                                     (Last[xa] - First[xa])/(n - 1)]], n - 1];

         If[xa[[l]] <= x,
            r = l + 1;
            While[r < n && xa[[r]] <= x,
                  l = r; k *= 2; r = Min[l + k, n]],
            {l, r} = {l - 1, l};
            While[1 < l && 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]];

         ({xa[[r]] - x, x - xa[[l]]}/(xa[[r]] - xa[[l]])).ya[[{l, r}]]],
         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
Improve this answer
$\endgroup$
4
$\begingroup$

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} *)
Improve this answer
$\endgroup$
3
  • $\begingroup$ 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…)? $\endgroup$
    – P. Fonseca
    May 2, 2012 at 16:37
  • $\begingroup$ @P.Fonseca Do you mean remove the previous definitions if they exceed a certain memory usage ? If so, sorry I don't. $\endgroup$
    – b.gates.you.know.what
    May 2, 2012 at 18:38
  • $\begingroup$ Yes. If the total myTest definitions exceed a certain quantity or memory usage $\endgroup$
    – P. Fonseca
    May 2, 2012 at 18:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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