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.

In ListPlot and Interpolation, the option InterpolationOrder->0 uses two different conventions.

Interpolation gives the y-value from the nearest greater x-value.

Plot[Quiet@Interpolation[{{0,0},{1,1},{2,2}},InterpolationOrder->0][x],{x,-1,3}]

enter image description here

ListLinePlot gives the y-value from the nearest smaller x-value.

ListLinePlot[{{0,0},{1,1},{2,2}},InterpolationOrder->0,PlotRange->{{-1,3},All}]

enter image description here

How can I make Interpolation use the ListPlot convention?

share|improve this question
1  
Great question +1. I should have asked this myself. –  Mr.Wizard Aug 8 '13 at 4:24
add comment

4 Answers 4

The question made me wonder about zero-order interpolations. It's hardly clear which is the right way. When I tried to figure out why ListLinePlot would use a different Interpolation, I noticed it didn't seem to use an Interpolation for orders 0 or 1 at all, but did it the simple way which you might use by hand: connect the dots. This was probably done for speed, and it is an important speed-up for linear interpolations. For a zero-order plot, there are at least two choices, "up/down and over" or "over and up/down" -- like walking around a block. Somehow the implementation for ListLinePlot chooses the way that does not agree with the choice of Interpolation. I tend to favor consistency, but this way there is a convenient way to get both plots.

It turns out there are three common zero-order interpolations. Of the two we see in the question, one sort of take the floor of the first coordinate, the datum whose first coordinate is nearest but less than or equal to an input x; the second, takes the ceiling with respect to the data in a similar way. The third method chooses the datum whose first coordinate is absolutely nearest x, like Voronoi cells. Indeed, ListPlot3D seems to use this method.

This third way led me to think that Nearest could be used to make step functions for all three ways that are faster than Piecewise on large data sets. They are called stepFloor (like Interpolation), stepCeiling (like ListLinePlot) and stepNearest and each returns a function.

Plot[stepFloor[steps][x],
     stepCeiling[steps][x]},
 {x, -1, 11},
 PlotStyle -> {
     Directive[Thick, Dashed],
     Thickness[Medium]}, 
 Epilog -> {Red, PointSize[Large], 
   Point@Transpose@{Range@Length@steps, steps}}, PlotRange -> {0, 10}]

Plot[stepNearest[steps][x],
 {x, -1, 11}, 
 PlotStyle -> {Thickness[Medium]}, 
 Epilog -> {Red, PointSize[Large], 
   Point@Transpose@{Range@Length@steps, steps}}, PlotRange -> {0, 10}]

Mathematica graphics Mathematica graphics

Functions

The three are implemented via the general stepFunction[stepFn_, steps_] that evaluates when applied to a numeric argument x (that is, the PatternTest x_?NumericQ prevents evaluation on symbolic argument x and avoids error messages from NearestFunction). For whatever reason, I wanted to make a decent interface, and I was led to incorporating intermediate, auxiliary functions. Perhaps there's a shorter way, but it takes a long time to write something short.

One of the helper functions is the intermediate stepFunction[nearest, steps, {test, next}], which can summarized as follows: it takes the nearest value in steps, and if test returns True, then it returns the next value instead. The argument nearest is a NearestFunction, and test and next are functions as well. The function nearest returns the index of the nearest datum and the desired one will be this one or the one next to it returned by next. In stepCeiling for instance, test is Less and next is #-1 &, or the previous element in steps, In other words, test indicates that nearest did not return the desired value and that it should be adjusted by next. (I waver whether the sense of test should be reversed -- a matter of taste, perhaps.)

There is also an intermediate step to sort the data steps, except in stepNearest where it is unnecessary.

ClearAll[stepFunction, stepCeiling, stepFloor, stepNearest];

stepFunction[stepFn_, steps_][x_?NumericQ] := stepFn[steps, x];
stepFunction[nearest_NearestFunction, steps_, {test_, next_}] := 
  stepFunction[
   With[{nearIDX = First@nearest@#2},
     If[test[#2, #1[[nearIDX, 1]]],
      #1[[Clip[next@nearIDX, {1, Length@#1}], 2]],
      #1[[nearIDX, 2]]]] &,
   steps
   ];
Format[stepFunction[args__]] := Shallow[HoldForm[stepFunction][args], {3, 2}];

stepFloor[steps_?VectorQ] := stepFloor[Transpose @ {N@Range@Length@steps, steps}, True];
stepFloor[steps_?MatrixQ] := stepFloor[Sort @ steps, True];
stepFloor[steps_?MatrixQ, True] := (* True => sorted *)
  stepFunction[Nearest[First /@ steps -> Automatic], steps, {Greater, # + 1 &}];

stepCeiling[steps_?VectorQ] := stepCeiling[Transpose @ {N@Range@Length@steps, steps}, True];
stepCeiling[steps_?MatrixQ] := stepCeiling[Sort @ steps, True];
stepCeiling[steps_?MatrixQ, True] := (* True => sorted *)
  stepFunction[Nearest[First /@ steps -> Automatic], steps, {Less, # - 1 &}];

stepNearest[steps_?VectorQ] := stepNearest[Transpose @ {N@Range@Length@steps, steps}];
stepNearest[steps_?MatrixQ] := With[{nf = Nearest[Rule @@ Transpose@steps]},
  stepFunction[First @ nf @ #2 &, {}]];

Efficiency

Time: As for speed, Mr.Wizard's Piecewise solution interp0 starts out about a bit faster than the corresponding stepCeiling, but as the length of the data grows, stepCeiling quickly becomes markedly faster. On the other hand, Plot recognizes Piecewise and the advantage in plotting interp0 lasts up until the length of the data exceeds a couple hundred. Plot generally produces a better-looking graph with interp0 than with stepCeiling.

In this table, the time to construct the function once plus the time to evaluate each function 1000 times on data of length 10^i:

Mathematica graphics

The result on lists of pairs is similar:

Mathematica graphics

Space: Here are the ByteCounts of the two functions and Interpolation on data of size 10^i, as well as the ratios between interp0 and stepCeiling (i/s) and between Interpolation and stepCeiling (I/s).

Mathematica graphics

Timing code

SeedRandom[2];
times = Table[
  dat = RandomReal[999, 10^i];
  {SeedRandom[i]; 
   First @ AbsoluteTiming[
      Module[{f = Function[x, Evaluate@interp0[dat, x]]}, Do[f@RandomReal[999], {1000}]]],
   SeedRandom[i]; 
   First @ AbsoluteTiming[
      Module[{f = stepCeiling[dat]}, Do[f@RandomReal[999], {1000}]]]},
   {i, 5}]

Plotting Manhattan Walks

While not precisely part of the OP's question, the difference between the plots of the methods of interpolation was raised. It seems that ListLinePlot with InterpolationOrder -> 0 simply constructs a "Manhattan path" through the points, that is, an alternating sequence of horizontal and vertical line segments. The function manhattanWalk below constructs such a path; the form manhattanWalk[list, Horizontal] begins with a Horizontal segment and matches what ListLinePlot does.

manhattanWalk[pts_, Left | Right | Horizontal | ___] := 
  Most @ Flatten[
    Transpose[{{#1, RotateLeft[#1]}, {#2, #2}} & @@ Transpose[pts], {1, 3, 2}], {{2, 3}}];
manhattanWalk[pts_, Up | Down | Vertical] := 
  Most @ Flatten[
    Transpose[{{#1, #1}, {#2, RotateLeft@#2}} & @@ Transpose[pts], {1, 3, 2}], {{2, 3}}];

We can see that LineLinePlot produces an equivalent plot to manhattanWalk (aside from styling), even when the points are not sorted.

SeedRandom[2];
data = RandomReal[10, {5, 2}];

Graphics[{Opacity[0.6], Line @ manhattanWalk[data], Red, 
  PointSize@Medium, Point[data], Blue, PointSize@Large, Point @ First @ data},
  ImageSize -> 250, Frame -> True]

Show[
 ListLinePlot[data, InterpolationOrder -> 0, 
   AspectRatio -> Automatic, Frame -> True, Axes -> False, ImageSize -> 250],
 PlotRange -> Automatic]

Mathematica graphics Mathematica graphics

This is not the same as Interpolation at all, which sorts the data (and does a manhattanWalk beginning with vertical move:

Module[{f = Interpolation[data, InterpolationOrder -> 0]},
 Plot[f[x], Evaluate[(First@f["Domain"] + {-0.01, 0})~Prepend~x], Frame -> True]]

Graphics[{Opacity[0.6], Line @ manhattanWalk[Sort @ data, Vertical],
  Red, PointSize@Medium, Point[data], Blue, PointSize@Large, Point @ First @ Sort @ data},
 ImageSize -> 250, Frame -> True, AspectRatio -> 1/GoldenRatio]

Mathematica graphics Mathematica graphics

share|improve this answer
    
I've used Nearest myself for the centered version, and I tried to use it for this problem as well but I failed. I don't have time now but I'll be very interested to see how you pulled this off. –  Mr.Wizard Aug 8 '13 at 22:02
    
Michael, what do you think of my new answer? –  Mr.Wizard Aug 9 '13 at 12:08
add comment
up vote 9 down vote accepted

You can negate the x-values in the input data and then negate the trial x-value.

Plot[Quiet@Interpolation[{{0,0},{1,1},{2,2}}.DiagonalMatrix[{-1,1}],InterpolationOrder->0][-x],{x,-1,3}]

enter image description here

This is a bit clumsy.

share|improve this answer
1  
"A bit clumsy," but very fast. –  Mr.Wizard Aug 9 '13 at 0:53
add comment

Since I can't think of anything else here is a method using a Piecewise object. For a limited number of points its quite fast but it slows down as the number increases. It's not a replacement for InterpolatingFunction but it can be handy at times.

SeedRandom[2]
dat = RandomInteger[999, {35, 2}];

ListLinePlot[SortBy[dat, {First}], InterpolationOrder -> 0]
Plot[interp0[dat, x], {x, 0, 999}, Exclusions -> None, Evaluated -> True]

enter image description here enter image description here

Also works with one-dimensional data:

SeedRandom[2]
dat = RandomReal[99, 35];

ListLinePlot[dat, InterpolationOrder -> 0]
Plot[interp0[dat, x], {x, 0, 36}, Exclusions -> None, Evaluated -> True]

enter image description here enter image description here

The optional second argument specifies the variable to use; it defaults to Formal x.

Code:

interp0[dat_?VectorQ, sym_: \[FormalX]] :=
 {Range[2, Length@dat], Most@dat}\[Transpose] //
  Piecewise[{#2, sym < #} & @@@ #, Last @ dat] &  

interp0[dat_, sym_: \[FormalX]] :=
 Module[{d},
  d = SortBy[dat, {First}];
  d[[All, 1]] = RotateLeft[d[[All, 1]]];
  Piecewise[{#2, sym < #} & @@@ d, d[[-1, 2]]]
 ]
share|improve this answer
1  
+1 Nice and simple. One of the nice things about Piecewise is that, if desired, Plot won't connect the steps when the Exclusions option is omitted. –  Michael E2 Aug 8 '13 at 20:36
add comment

It has occurred to me that I was thinking about this problem all wrong:
We merely need a binary search within the index list.

SeedRandom[2]
dat = RandomInteger[999, {35, 2}];

ListLinePlot[SortBy[dat, {First}], InterpolationOrder -> 0, Epilog -> Point@dat]

Needs["Combinatorica`"]
Module[{d1, d2},
 {d1, d2} = SortBy[dat, {First}]\[Transpose];
 Plot[d2[[Floor @ BinarySearch[d1, x]]], {x, 0, 999}, Epilog -> Point@dat]
]

enter image description here enter image description here

A binary search can also be quite efficient. I will use Leonid's bsearchMin code, compiled:

bsInterp0[dat_] :=
  Module[{d1, d2, cf},
   {d1, d2} = SortBy[dat, {First}]\[Transpose];
   With[{list = d1},
    cf = Compile[{elem},
      Module[{n0 = 1, n1 = Length[list], m = 0},
       While[n0 <= n1, m = Floor[(n0 + n1)/2];
        If[list[[m]] == elem, While[list[[m]] == elem, m++];
         Return[m - 1]];
        If[list[[m]] < elem, n0 = m + 1, n1 = m - 1]];
       If[list[[m]] < elem, m, m - 1]]
    ]
   ];
   d2[[cf[#]]] &
  ];

Timings:

SeedRandom[2]
dat = RandomReal[999, {500000, 2}];

if = Interpolation[dat, InterpolationOrder -> 0];
bsf = bsInterp0[dat];
scf = stepCeiling[dat];  (* Michael's function *)

if  /@ Range[0, 999, 0.01] // Timing // First
bsf /@ Range[0, 999, 0.01] // Timing // First

0.094

0.514

On a set with half a million elements it is ~5.5X slower than InterpolatingFunction.

Compared to Michael's stepCeiling with the same data it is ~45X faster:

bsf /@ Range[0, 999, 0.05] // Timing // First
scf /@ Range[0, 999, 0.05] // Timing // First

0.1062

4.758

share|improve this answer
    
Ouch! that's fast. +1 –  Michael E2 Aug 9 '13 at 12:26
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.