Taking Part of an InterpolatingFunction
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 | NDSolve`CubicHermite | NDSolve`LocalSeries | 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", Developer`PackedArrayForm] := 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 : {Developer`PackedArrayForm, _, _},
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"]]] =
Developer`ToPackedArray@{coords[[Span @@ span]]};
newif[[ifnPart["InterpolationData"]]] =
iDataTake[method, if[[ifnPart["InterpolationData"]]], span];
newif[[ifnPart["InterpolationStructure"]]] =
iStructureTake[method, if[[ifnPart["InterpolationStructure"]]], span];
InterpolatingFunction @@ newif
];
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]] = Developer`ToPackedArray@{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:
- @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. - I'm assuming that the
InterpolatingFunctionstores only values and first derivatives. That seems to be the case for all the ones I've seen, but maybe isn't general. - Messing with
InterpolatingFunctioninternals still seems super sketchy and liable to break with new or old Mathematica versions. - This took a lot of trial & error, so there still might be issues!
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.