Take some data, return a pure piecewise function

Question

Given a list of {x,y} datapoints, return a pure function f (from the reals to the reals) such that f[x]==y for every {x,y} in the data. If x is not one of the x-values then return the y-value for the previous point (the one with x-value less than x). If the function gets a value less than that of the first x-value in the data -- i.e., there is no previous point -- then return 0.

For example, given data {{1,20}, {2,10}}, return a pure function that looks like this:

Graph of the function given {{1,20},{2,10}} http://yootles.com/outbox/so/piecewise.png

I wrote something using Function and Piecewise that I'll include as an answer but it seems like it might be inefficient, especially for a large list of points. [UPDATE: My answer may actually be decent now. I'll probably go with it if no one has better ideas.]

To be clear, we're looking for function that takes a single argument -- a list of pairs of numbers -- and returns a pure function. That pure function should take a number and return a number.

Answer 1

Hand-Coded Binary Search

If one is willing to sacrifice conciseness for performance, then an imperative binary search approach performs well:

stepifyWithBinarySearch[data_] :=
  With[{sortedData = SortBy[data, First], len = Length @ data}
  , Module[{min = 1, max = len, i, x, list = sortedData}
    , While[min <= max
      , i = Floor[(min + max) / 2]
      ; x = list[[i, 1]]
      ; Which[
          x == #, min = max = i; Break[]
        , x < #, min = i + 1
        , True, max = i - 1
        ]
      ]
    ; If[0 == max, 0, list[[max, 2]]]
    ]&
  ]

Equipped with some test scaffolding...

test[s_, count_] :=
  Module[{data, f}
  , data = Table[{n, n^2}, {n, count}]
  ; f = s[data]
  ; Timing[Plot[f[x], {x, -5, count + 5}]]
]

...we can test and time various solutions:

test[stepifyWithBinarySearch, 10]

On my machine, the following timings are obtained:

test[stepify (*version 1*), 100000]      57.034 s
test[stepify (*version 2*), 100000]      40.903 s
test[stepifyWithBinarySearch, 100000]     2.902 s

I expect that further performance gains could be obtained by compiling the various functions, but I'll leave that as an exercise for the reader.

Better Still: Precomputed Interpolation (response To dreeves' comment)

It is baffling that a hand-coded, uncompiled binary search would beat a Mathematica built-in function. It is perhaps not so surprising for Piecewise since, barring optimizations, it is really just a glorified IF-THEN-ELSEIF chain testing expressions of arbitrary complexity. However, one would expect Interpolation to fare much better since it is essentially purpose-built for this task.

The good news is that Interpolation does provide a very fast solution, provided one arranges to compute the interpolation only once:

stepifyWithInterpolation[data_] :=
  With[{f=Interpolation[
            {-1,1}*#& /@ Join[{{-9^99,0}}, data, {{9^99, data[[-1,2]]}}]
            , InterpolationOrder->0 ]}
    , f[-#]&
  ]

This is blindingly fast, requiring only 0.016 seconds on my machine to execute test[stepifyWithInterpolation, 100000].

Answer 2

You could also do this with Interpolation (with InterpolationOrder->0) but that interpolates by using the value of the next point instead of the previous. But then I realized you can reverse that with a simple double-negation trick:

stepify[data_] := Function[x,
  Interpolation[{-1,1}*#& /@ Join[{{-9^99,0}}, data, {{9^99, data[[-1,2]]}}],
                InterpolationOrder->0][-x]]