# Taking Part of an InterpolatingFunction

I want a function to extract part of a one-dimensional InterpolatingFunction, similar to how Take works on lists. Based on the code in this question by polyglot, I put together the following function:

InterpolatingFunctionTake[if_InterpolatingFunction, {xmin_?NumericQ, xmax_?NumericQ}] :=
Interpolation[Join[
{{xmin, if[xmin]}},
Select[Transpose[{if["Coordinates"][[1]], if["ValuesOnGrid"]}], xmin < #[[1]] < xmax &],
{{xmax, if[xmax]}}
]];


It generally works OK, but has two problems that I'd like to fix.

Problem 1 - It can be slow

When the InterpolatingFunction is big, it can take too long to run. For example,

if = n /. NDSolve[{n'[t] == n[t] (1 - n[t - 2]), n[0] == 0.1}, n, {t, 0, 20000}, MaxSteps -> ∞][[1]];

InterpolatingFunctionTake[if, {0.5, 1.5}]


takes almost 1 second for the InterpolatingFunctionTake.

Can this be sped up?

Problem 2 - It fails on discontinuous functions

Suppose there's a discontinuity in the function due to an impulse, modeled with a WhenEvent in NDSolve. For example,

if = n /. NDSolve[{n'[t] == 1, n[0] == 0, WhenEvent[Mod[t, 1], n[t] -> 0]}, n, {t, 0, 2}][[1]];

InterpolatingFunctionTake[if, {0.5, 1.5}]


gives the error "Interpolation::inddp: The point 1. in dimension 1 is duplicated." Simply deleting one of the duplicate points is no good because it completely messes up the interpolation.

This very issue was explored in depth in this answer by Michael E2. His solution involves manually building up a new InterpolatingFunction based on a list of {x,y,y'} triplets. It seems that once we're already that deep into the innards of the new InterpolatingFunction, we might as well work directly on the parts of the original InterpolatingFunction without making a list of {x,y,y'} triplets (particularly since there doesn't seem to be anything like if["DerivativesOnGrid"] to get the y' values out easily).

Is there a good way to directly alter the internals of the InterpolatingFunction? (This seems slightly dangerous, but maybe the only way to achieve this?)

• if["DerivativesOnGrid"] would be equivalent to if'["ValuesOnGrid"]. – Michael E2 Jul 20 '17 at 0:35
• There are many kinds of InterpolatingFunction. Do you only want to deal with the kind that has the DeveloperPackedArrayForm kind of interpolation data? – Michael E2 Jul 20 '17 at 0:39
• @MichaelE2 I'm only interested slicing up the 1D InterpolatingFunctions that NDSolve returns when solving time-dependent differential equations. – Chris K Jul 20 '17 at 0:53
• In fact for me (V11.1.1), your example has if["InterpolationMethod"] --> "Local Taylor series", and it is not of the PackedArrayForm you show. – Michael E2 Jul 20 '17 at 0:53
• @MichaelE2 Good catch - the first example's InterpolationMethod is "Local Taylor series" (and has no PackedArray in if[[4]]) and the second example's InterpolationMethod is "Hermite". So attempting InterpolatingFunction surgery might be too risky ;) – Chris K Jul 20 '17 at 0:58

This handles the "Local Taylor series", "Chebyshev" and both the packed and unpacked "Hermite" interpolation method. One thing I would recommend, and implemented below, is not to worry about cutting off the interpolation exactly at tmin or tmax. Rather, take the smallest interval containing tmin and tmax with endpoints on the interpolation grid. Actually the code below might include one more grid point at each end than is necessary. This adds very little overhead. On the other hand, dealing with dividing a subinterval with these methods might be a much bigger pain than it's worth. In any case, the function in the "Code dump" below preserves these (more accurate!) interpolation methods. It seems a shame to throw away the extra accuracy. One might do it if speed is more important than accuracy.

OP's first example:

if1 = n /. NDSolve[{n'[t] == n[t] (1 - n[t - 2]), n[0] == 0.1},
n, {t, 0, 20000}, MaxSteps -> ∞][[1]];
if1["InterpolationMethod"]
(*  "Local Taylor series"  *)

ifTake[if1, {0.5, 1.5}] // AbsoluteTiming
(*  {0.028446, InterpolatingFunction[{{0.478428, 1.53819}}, <>]}  *)

ifTake[if1, {1000.5, 1010.5}] // AbsoluteTiming
(*  {0.03142, InterpolatingFunction[{{1000.46, 1010.55}}, <>]}  *)


Chebyshev example:

if3 = NDSolveValue[{y''[x] + y[x] == 0, y[0] == 1, y'[0] == 0},
y, {x, 0, 100}, InterpolationOrder -> All, Method -> "Extrapolation"];
if3["InterpolationMethod"]
(*  "Chebyshev"  *)

ifTake[if3, {10.5, 11.5}]
(*  InterpolatingFunction[{{9.03213, 12.7461}}, <>]  *)


Code dump

InterpolatingFunction parts. Execute

ifnPart["Properties"] (* OR *)
DownValues[ifnPart][[;; Length@ifnPart["Properties"] + 1]]


to see a list of the part names. Some are valid only for the "Hermite" method, some for "Chebyshev" or "Local Taylor Series". The function of some of the parts are as yet unknown to me. It as complete as my current knowledge, and more than is needed. But it seemed worth sharing.

ClearAll[ifnPart];
ifnPart["Domain"] = Sequence[1];   (*bounding box for domain*)
ifnPart["X1"] = Sequence[1, 1];    (* lower bound for first coordinate *)
ifnPart["X2"] = Sequence[1, 2];    (* upper bound for first coordinate *)
ifnPart["Version"] = Sequence[2, 1];
ifnPart["Flags"] = Sequence[2, 2]; (*flags indicating properties:
bit field positions - inferred, perhaps mistaken
$extrapolation=0; whether to warn about extrapolation$fullArrayBit=1; interpolation data is a full array (not ragged)
$packed=2; packed array form (???)$repeatedBit=4; whether repeated abscissae are permitted*)
ifnPart["DerivativeOrder"] = Sequence[2, 3]; (*max derivative order*)
ifnPart["NGrid"] = Sequence[2, 4];  (*number of points in each coordinate grid*)
ifnPart["InterpolationOrder"] = Sequence[2, 5]; (*interpolation order*)
ifnPart["Derivative"] = Sequence[2, 6];  (*derivative to evaluate:0-->f[x], 1-->f'[x],...*)
ifnPart["Periodic"] = Sequence[2, 7];
(*ifnPart["??"]=Sequence[2,8];*)
(*ifnPart["??"]=Sequence[2,9];*)
ifnPart["ExtrapolationHandler"] = Sequence[2, 10];
(*ifnPart["??"]=Sequence[2,11];*)
(*ifnPart["??"]=Sequence[2,12];*)
(*ifnPart["??"]=Sequence[2,13];*)
ifnPart["Coordinates"] = Sequence[3];  (*list of lists, abscissae of interpolation grid*)
ifnPart["InterpolationData"] = Sequence[4]; (*interpolation data (values or coefficients)*)
ifnPart["Offsets"] = Sequence[4, 2]; (*offsets in function/derivative array (PackedArrayForm)*)
ifnPart["FlatData"] = Sequence[4, 3]; (*flattened function/derivative values (PackedArrayForm)*)
ifnPart["InterpolationStructure"] = Sequence[5]; (*{Automatic}, or dense output interpolation structure:
list of types for each unit/subinterval*)
ifnPart["UnitIndices"] = Sequence[5, 1, 1]; (*dense output:
Indices (to grid) for corresponding coefficients*)
ifnPart["UnitTypes"] = Sequence[5, 1, 2];   (*dense output types:
Automatic | NDSolveCubicHermite | NDSolveLocalSeries | ChebyshevT*)

ifnPart["Properties"] =
Cases[DownValues[ifnPart], Verbatim[ifnPart][prop_] :> prop, Infinity];

ifnPart["ValidPartQ", "Chebyshev" | "Local Taylor Series", "UnitIndices" | "UnitTypes", _] := True;
ifnPart["ValidPartQ", _, "UnitIndices" | "UnitTypes", _] := False;
ifnPart["ValidPartQ", "Hermite", "Offsets" | "FlatData", DeveloperPackedArrayForm] := True;
ifnPart["ValidPartQ", _, "Offsets" | "FlatData", _] := False;
ifnPart["ValidPartQ", method_String, part_String, _] /;
MemberQ[method, "Chebyshev" | "Local Taylor Series" | "Hermite"] &&
MemberQ[part, ifnPart["Properties"]] := True;
ifnPart["ValidPartQ", _, _, _] := False;
ifnPart[if_InterpolatingFunction, part_String] /;
ifnPart["ValidPartQ", if["InterpolationMethod"], part, if[[4, 1]]] :=
if~Part~ifnPart[part];


Taking part of an InterpolatingFunction:

ClearAll[ifTake];
dupeLast[list_] := Append[list, Last@list];
iDataTake["Local Taylor series" | "Chebyshev", data_, span_] := Join[
{data[[First@span, 1 ;; 2]]}, data[[First@span + 1 ;; Last@span]]
];
iDataTake["Hermite", data : {DeveloperPackedArrayForm, _, _},
span : {s1_, s2_}] := ReplacePart[
data,
{Rest@{ifnPart["Offsets"]} ->
data[[2, s1 ;; s2 + 1]] - data[[2, s1]],
Rest@{ifnPart["FlatData"]} ->
data[[3, data[[2, s1]] + 1 ;; data[[2, s2 + 1]] ]]}
];
iDataTake["Hermite", data : {__List}, span_] := data[[Span @@ span]];
iStructureTake["Local Taylor series" | "Chebyshev", structure_,
span_] := ReplacePart[structure,
{Rest@{ifnPart["UnitIndices"]} ->
Join[
{{1}},
1 + structure[[##2 &@ifnPart["UnitIndices"], First@span + 1 ;; Last@span]] -
structure[[##2 &@ifnPart["UnitIndices"], First@span, -1]] //
dupeLast
],
Rest@{ifnPart["UnitTypes"]} ->
Join[
{Automatic},
structure[[##2 &@ifnPart["UnitTypes"],
First@span + 1 ;; Last@span]] // dupeLast
]}
];
iStructureTake["Hermite", structure_, span_] := structure;
ifTake[if_InterpolatingFunction, {tmin_?NumericQ, tmax_?NumericQ}] /;
Length@if["Domain"] == 1 :=
Module[{coords, newif = Hold @@ if, span, method},
method = if["InterpolationMethod"];
coords = First@if["Coordinates"];
span = Clip[
SparseArray[UnitStep[coords - tmin] UnitStep[tmax - coords]][
"AdjacencyLists"][[{1, -1}]] + {-1, 1}, {1, Length@coords}];

newif[[ifnPart["Domain"]]] =
{coords[[span]]};
(*newif[[ifnPart[if,"Flags"]]]=??;
newif[[ifnPart[if, "DerivativeOrder"]}]]=??;*) (* not needed?? *)
newif[[ifnPart["NGrid"]]] =
1 + Differences@span;
newif[[ifnPart["Coordinates"]]] =
DeveloperToPackedArray@{coords[[Span @@ span]]};
newif[[ifnPart["InterpolationData"]]] =
iDataTake[method, if[[ifnPart["InterpolationData"]]], span];
newif[[ifnPart["InterpolationStructure"]]] =
iStructureTake[method, if[[ifnPart["InterpolationStructure"]]], span];

InterpolatingFunction @@ newif
];

• Thanks! Looks like all of the InterpolatingFunctions I commonly get as solutions to first order initial value problems in NDSolve have InterpolatingMethod Hermite. if[[4, 3, ;; ;; 2]] is a faster equivalent to if["ValuesOnGrid"] on the examples I've looked at, but I've already been surprised by the variety of possible structures inside InterpolatingFunctions, so it needs more testing tomorrow. – Chris K Aug 2 '17 at 3:51
• Here's an example that generates a Hermite InterpolatingMethod: sol = NDSolve[{n'[t] == (Sin[2 \[Pi] t] + 1) n[t] - n[t]^2, n[0] == 0.01}, n, {t, 0, 100000}][[1]]. Timing[g1 = (n /. sol)["ValuesOnGrid"];] is 1 second, but Timing[g2 = (n /. sol)[[4, 3, ;; ;; 2]];] is 0.017s, and g1==g2. – Chris K Aug 2 '17 at 3:54

Most (if not all) of my actual InterpolatingFunctions use InterpolationMethod Hermite, so I tried tweaking @Mr.Wizard's and @MichaelE2's answers to deal with them, while still addressing discontinuous functions too. In the end, @MichaelE2's surgical reconstruction approach was easier to get working. Here's my attempt:

ifTake[if_InterpolatingFunction, {tmin_?NumericQ, tmax_?NumericQ}]
/; if["InterpolationMethod"] == "Hermite" :=
Module[{coords, newif = Hold @@ if, span},
coords = First@if["Coordinates"];
span = Clip[SparseArray[UnitStep[coords - tmin] UnitStep[tmax - coords]]["AdjacencyLists"][[{1, -1}]] + {-1, 1}, {1, Length@coords}];
newif[[1]] = {{tmin, tmax}};
newif[[2, 4]] = 1 + Differences@span;
newif[[3]] = DeveloperToPackedArray@{coords[[Span @@ span]]};
newif[[4, 2]] = if[[4, 2, ;; (span[[2]] - span[[1]]) + 2]];
newif[[4, 3]] = if[[4, 3, 2 span[[1]] - 1 ;; 2 span[[2]]]];
InterpolatingFunction @@ newif
];


Here's an example that was too slow with my original, naive InterpolatingFunctionTake:

if = n /. NDSolve[{n'[t] == (Sin[2 π t] + 1) n[t] - n[t]^2, n[0] == 0.01},
n, {t, 0, 10^5}][[1]];
if["InterpolationMethod"]
(* Hermite *)

InterpolatingFunctionTake[if, {99000, 100000}] // Timing
(* 5.74 sec *)

ifTake[if, {99900, 100000}] // Timing
(* 0.163 sec *)


It also works on my discontinuous example 2 in the original post. In both cases there is zero difference between the original and trimmed InterpolatingFunction.

Some notes:

1. @MichaelE2's trick of not worrying about the precise endpoints is great, but I actually also require the length to be exact. I used newif[[1]] = {{tmin, tmax}} to directly set the domain, which seems to work.
2. I'm assuming that the InterpolatingFunction stores only values and first derivatives. That seems to be the case for all the ones I've seen, but maybe isn't general.
3. Messing with InterpolatingFunction internals still seems super sketchy and liable to break with new or old Mathematica versions.
4. This took a lot of trial & error, so there still might be issues!
• +1. I meant to come back to this. I was trying to think of an approach that would incorporate all the methods elegantly, but other work has delayed it. I think this method will work only on first-order ODE solutions. Try it on NDSolve[{y''[x] + y[x] == 0, y[0] == 1, y'[0] == 1}, y, {x, 0, 10}]. The PackedArrayForm can have arbitrary index differences. (Then there is the unpacked form produced by Interpolation but not apparently by NDSolve. One can convert a packed IF to an unpacked one with if["Unpack"].) – Michael E2 Aug 6 '17 at 18:56
• @MichaelE2 Yes, that's exactly the kind of thing I was worried about. Cheers for the amazing solution you provided in your answer. – Chris K Aug 7 '17 at 14:58

With help from a few clues in @MichaelE2's comments, I think I've solved Problem 2 (discontinuous InterpolatingFunctions). First, define discontInterpolation from this answer. Then define

discontInterpolatingFunctionTake[if_InterpolatingFunction, {xmin_?NumericQ, xmax_?NumericQ}] :=
discontInterpolation[Join[
{{{xmin}, if[xmin], if'[xmin]}},
Select[
Transpose[{{#} & /@ if["Coordinates"][[1]], if["ValuesOnGrid"], if'["ValuesOnGrid"]}],
xmin < #[[1, 1]] < xmax &],
{{{xmax}, if[xmax], if'[xmax]}}]];


It works perfectly on my problem 2:

if = n /. NDSolve[{n'[t] == 1, n[0] == 0, WhenEvent[Mod[t, 1], n[t] -> 0]}, n, {t, 0, 2}][[1]]
ift2 = discontInterpolatingFunctionTake[if, {0.5, 1.5}];

Show[
Plot[if[t], {t, 0, 2}],
Plot[ift2[t], {t, 0.5, 1.5}, PlotStyle -> Pink]
]


The only rub is this solution exacerbates my problem 1 (too slow on large InterpolatingFunctions). On that example,

if = n /.
NDSolve[{n'[t] == n[t] (1 - n[t - 2]), n[0] == 0.1}, n, {t, 0, 20000}, MaxSteps -> ∞][[1]];

discontInterpolatingFunctionTake[if, {0.5, 1.5}]//Timing
(* ~1.5 sec *)


Based on another tip from MichaelE2 in chat, it seems that if["ValuesOnGrid"] is responsible for the slowness. Each of the two ValuesOnGrids take ~0.6 sec and the Select[Transpose[...]] bit takes ~0.3 sec.

Maybe manual surgery on the InterpolatingFunction internals could be a way to speed it up, but those internals are just too complex for me to attempt it myself.

• I see that you already knew about the cost of "ValuesOnGrid" and Select but I hope you find my answer of value nevertheless. – Mr.Wizard Aug 1 '17 at 17:11

A couple of observations.

1. Method "ValuesOnGrid" is very slow, and I don't know why. At least for the example given it can be replaced with a simple Part extraction.

2. Select as used is inefficient; a binary search would be better; see: Finding all elements within a certain range in a sorted list

Applying those methods:

Quiet[Needs["Combinatorica"]]

ifnTake1[
if_InterpolatingFunction,
{xmin_?NumericQ, xmax_?NumericQ}
] :=
Module[{c, v, bs},
bs = CombinatoricaBinarySearch;
c = if[[3, 1]];         (* warning: magic numbers *)
v = if[[4, All, 1]];    (* warning: magic numbers *)
Join[
{{xmin, if[xmin]}},
{c, v}\[Transpose][[⌈bs[c, xmin]⌉ ;; ⌊bs[c, xmax]⌋]],
{{xmax, if[xmax]}}
] // Interpolation
]


Timings:

fn1 = InterpolatingFunctionTake[if, {0.5, 1.5}]; // RepeatedTiming
fn2 = ifnTake1[if, {0.5, 1.5}];                  // RepeatedTiming
fn1 === fn2

{0.792, Null}

{0.0318, Null}

True


Caveat: I believe that other InterpolatingFunction expressions have a different internal data structure so this is not a general solution as written. However it would probably be possible to detect these different formats and extract the needed parts accordingly.

• Thanks for the start. Unfortunately I need a general approach that works on any type of InterpolatingFunction, which is beyond me right now! – Chris K Aug 1 '17 at 18:05
• More generally, a great solution would speed up Method "ValuesOnGrid" for all purposes. – Chris K Aug 1 '17 at 18:16
• @ChrisK If you provide me with a series of InterpolatingFunction expressions to use as a target I shall attempt that. – Mr.Wizard Aug 1 '17 at 19:49
• @ChrisK, Mr. Wizard: The "ValuesOnGrid" is slow on a packed "Hermite" interpolation because it stores a potentially ragged array of function, derivative and higher derivative values as a flattened, packed array of those values and a flattened, packed array of indices (suitable for InternalPartitionRagged IIRC). Extracting the values is done in a way that unpacks the values and apparently other arrays (use On["Packing"]). – Michael E2 Aug 2 '17 at 11:58