Taking Part of an InterpolatingFunction

Question

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

Answer 1

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)
      $machPrecBit=2; whether the data (f,f',...) is MachinePrecision
      $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
   ];

Answer 2

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:

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!

Answer 3

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.

Answer 4

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 = Combinatorica`BinarySearch;
   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.