8
$\begingroup$

Bug introduced in 10.1 or earlier and fixed in 10.4.0


Using Indexed on an InterpolatingFunction leads to an $IterationLimit::itlim error.

{foo} = NDSolve[{y'[t] == y[t], y[0] == {1, 2}}, y, {t, 0, 1}];

sol = Table[Indexed[y[t], i], {i, 2}];
sol /. foo

\$IterationLimit::itlim: Iteration limit of 4096 exceeded. >>
\$IterationLimit::itlim: Iteration limit of 4096 exceeded. >>

Mathematica graphics

It works find if t has a numeric value (since the InterpolatingFunction will evaluate to a List before Indexed is evaluated):

Block[{t = 0.5}, sol /. foo]
(*  {1.64872, 3.29744}  *)

This came up in trying to Plot the components of the solutions, which combines both of the above.

Plot[sol /. foo // Evaluate, {t, 0, 1}]
(*  error messages and a nice two-tone plot  *)

Is this a bug, or am I using Indexed wrong?

Side notes:

(1) Here is a simpler example, easier for tracing (Block[{$IterationLimit = 20}, Trace[...]]):

    foo = y -> Interpolation[{{{0}, {1, 2}}, {{1}, {3, 4}}}, InterpolationOrder -> 1];

or, in essence, indexing just one part,

if = Interpolation[{{{0}, {1, 2}}, {{1}, {3, 4}}}, InterpolationOrder -> 1];
Indexed[if[t], 2]

(2) Workarounds may be found in the answers to Prevent Part[] from trying to extract parts of symbolic expressions

Confirmed by WRI.

$\endgroup$
  • $\begingroup$ Does anyone agree that this question should be tagged a bug? $\endgroup$ – Michael E2 Jan 10 '16 at 14:47
4
$\begingroup$

Indexed is hardly well-documented, so it is hard to say what is or what is not a bug. However, I will argue that what you encountered is not a bug. The basis of my argument is that what documentation there is indicates the first argument given to Indexed should be an expression that is a valid first argument for Part. An interpolation function does not meet this condition.

f = NDSolve[{y'[t] == y[t], y[0] == {1, 2}}, y, {t, 0, 1}][[1, 1, 2]]
Table[Part[f[t], i], {i, 2}]

Part::partw: Part 2 of InterpolatingFunction[{{0.,1.}},{5,3,1,{26},{4},0,0,0,0,Automatic,{},{},False},{{<<1>>}},{{{1.,2.},{1.,2.}},<<24>>,{{2.71828,5.43656},{2.71828,<<19>>}}},{Automatic}][t] does not exist. >>

I am arguing that one can only use Indexed safely in situations where

 Indexed`[<expression>, <index>] /. Indexed -> Part

produces a valid result.

$\endgroup$
  • $\begingroup$ Clear[x]; Indexed[x[t], 2] does not produce an error, even though Part would. It has only happened for me when the head is InterpolatingFunction, although there are lots of heads I haven't checked. $\endgroup$ – Michael E2 Sep 5 '15 at 22:12
  • $\begingroup$ Thanks for your ideas. My first thought is that you are onto something, but my feeling is that the connection to Part is stated too strongly. I thought perhaps Indexed evaluated effectively to Part on any expression that has parts, but the example of my comment shows that it not quite right. (Of course, the docs mention it working on List but don't say that it's the only head it will evaluate on.) $\endgroup$ – Michael E2 Sep 6 '15 at 0:14
  • $\begingroup$ @MichaelE2. I accept that I may have overstated the constraint. It's hard to interpret the sparse documentation. I feel pretty confident about "if it works with Part, it will work with Indexed. I am quite uncertain about what qualifies as a necessary condition. $\endgroup$ – m_goldberg Sep 6 '15 at 2:29
3
$\begingroup$

[Ref: Please see The clearest way to represent Mathematica's evaluation sequence for WReach's traceView2 function, which I will use to trace evaluation.]

I'm a little more confused than I like to be when I answer. Perhaps sharing what I think I know will prompt others who better understand evaluation in Mathematica to contribute.

One of the mysterious things about Mathematica is how it determines when to evaluate an expression as it moves through the digestion process of evaluation. No doubt it's perfectly clear to those who know the internal workings. What I mean is that, putting held expressions and so forth aside, once Mathematica has evaluated an InterpolatingFunction or an Integrate that results in the same (or even a different) expression, it does not need to reevaluate it -- unless something changes which either does or might lead to a different result. There must be something like a been-there-done-that flag associated with the expression, but however it is done, I have to hope and treat it like it does not matter, because I don't how it's done.

Working answer: What seems to be happening is that inside Indexed something happens that suggests to Mathematica that something changed that might lead to a different result for the InterpolatingFunction. The InterpolatingFunction is reevaluated and that causes Indexed to be reevaluated. This leads to an infinite loop, since inside Indexed, once again the InterpolatingFunction is reevaluated, etc.

If you wish to explore this, here is a minimal example:

if = Interpolation[Range@2, InterpolationOrder -> 1];
Indexed[if, 2]

Unlike the example in the OP, using Indexed this way does not make sense; however, it does give the same $IterationLimit::itlim behavior as the original with a minimal-sized trace.

You might first want to see the evaluation process that constructs an InterpolatingFunction. First Interpolation constructs an InterpolatingFunction expression: The standard evaluation sequence, head Interpolation, first argument Range@2, etc. results in

InterpolatingFunction[{{1, 2}}, {5, 3, 0, {2}, {2}, 0, 0, 0, 0, 
  Automatic, {}, {}, False}, {{1, 2}}, {{1}, {2}}, {Automatic}]

which is itself evaluated in six steps.

traceView2[Interpolation[Range@2, InterpolationOrder -> 1]]

Mathematica graphics

The final result has the same input form, but is formatted(!). I don't know for sure why it is formatted, but it seems consistent that it is formatted when Mathematica is done evaluating the expression. I have assumed this and used it as a key in interpreting the trace of the evaluation.

The evaluation of InterpolatingFunction[..] also follows the standard evaluation sequence and results in an expression with the same input form, although now that it is evaluated, it will be formatted when returned as seen above.

Mathematica graphics

Another trick to help read the trace output is to replace all InterpolationFunction[..] expressions by "IFN". You lose the formatting, but the output is smaller and easier to deal with. Since the expressions in a trace are wrapped in HoldForm, this replacement does not cause reevaluation of the trace output.

Here we can see the start of the infinite loop described in the working answer. After if is evaluated, Indexed[IFN, 2] evaluates to Indexed[IFN, {2}], which begins the loop. Note that the InterpolatingFunction is reevaluated (the 7-step block, not expanded). So after the head Indexed and the arguments IFN and {2}, Indexed[IFN, {2}] is evaluated. Instead of evaluation stopping at this point, as it should, the evaluation process is repeated.

Block[{$IterationLimit = 20},
  traceView2[Indexed[if, 2]]
  ] /. i_InterpolatingFunction :> "IFN"

$IterationLimit::itlim: Iteration limit of 20 exceeded. >>

Mathematica graphics

I haven't found another function with this problem, so far. Only 5000 more to check. :)

$\endgroup$
  • 1
    $\begingroup$ This was confirmed to be the bug by WRI. $\endgroup$ – Michael E2 Jan 19 '16 at 1:00
2
$\begingroup$

This is an extended comment, not an answer.

After running into problems similar to yours with evaluations that returned rules involving interpolation functions, I decided to cut the Gordian knot and simply extract the interpolation function(s) for use in any further computations. I have never regretted adopting this strategy and it certainly works well for this problem.

f = NDSolve[{y'[t] == y[t], y[0] == {1, 2}}, y, {t, 0, 1}][[1, 1, 2]];
Plot[f[t], {t, 0., 1.}]

plot

Voila, no fuss, no muss, no bother (old tv commercial catch phrase).

Update

Addressing an issue raised in a comment.

{f1, f2} = Table[With[{i = i}, f[#][[i]] &], {i, 2}];
Plot[{f1[t], f2[t]}, {t, 0., 1.}]

color-plot

$\endgroup$
  • $\begingroup$ Can you figure out how to style each curve separately? $\endgroup$ – Michael E2 Sep 5 '15 at 22:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.