# How can the behavior of InterpolationOrder->0 be controlled?

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}]


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

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


How can I make Interpolation use the ListPlot convention?

-
Great question +1. I should have asked this myself. –  Mr.Wizard Aug 8 '13 at 4:24

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}]


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:

The result on lists of pairs is similar:

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

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]


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]


-
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

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}]


This is a bit clumsy.

-
"A bit clumsy," but very fast. –  Mr.Wizard Aug 9 '13 at 0:53

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]


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]


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]]]
]

-
+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

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]
]


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

-
Ouch! that's fast. +1 –  Michael E2 Aug 9 '13 at 12:26