dense_idxs

[devmotion]

As discussed in SciML/DelayDiffEq.jl#68 (comment), it would be nice to be able to specify a subset of indices for which a dense solution is required which can be different from save_idxs (in particular thinking about higher-dimensional DDEs, for which a solution of all components is required at certain time points without dense interpolation but integration requires a dense interpolation of a small subset of indices). The proposed idea was basically, as far as I understood, to introduce an option which could be named dense_idxs or similar and be set to save_idxs by default that would be used for copying of k in https://github.com/JuliaDiffEq/OrdinaryDiffEq.jl/blob/master/src/integrators/integrator_utils.jl#L81-L85.

However, this would break the interpolation of the resulting solution in interesting cases mentioned above in which save_idxs != dense_idxs. Hence I thought about creating a new CompositeInterpolation (the name can be discussed, of course) that would use a suitable interpolation for dense_idxs and the default linear or Hermite interpolation for save_idxs that are not in dense_idxs. Does this seem like a good idea? Is there a better approach? In particular, I was wondering how one could efficiently compute interpolations if the provided idxs are spread over different indices in dense_idxs and other save_idxs - but maybe there is no better way than just computing the interpolations for both subsets and then combining them?

[ChrisRackauckas]

Saveat happens during the stepping time, so it doesn't require the k's being saved at all. So it can always use the full interpolation no matter what. dense_idxs could just cut down to a smaller set of k to do the same interpolation. I don't think there's a need to fall back to linear or Hermite here.

[devmotion]

I guess I did not explain the problem I see well enough, or maybe it's not a problem actually. So I try to illustrate what I mean.

I added the changes you suggested (or at least the changes I assumed you suggested) to a "dense_idxs" branch in my fork. However, as I expected, these changes yield incorrectly interpolated values even for dense_idxs:

using OrdinaryDiffEq, DiffEqProblemLibrary

sol_dense = solve(prob_ode_2Dlinear, Tsit5())
sol_nodense = solve(prob_ode_2Dlinear, Tsit5(); dense = false)
sol_denseidxs = solve(prob_ode_2Dlinear, Tsit5(); dense_idxs = [2, 4])

The following plots visualize the resulting interpolation:

using Plots
plotly()

plot(sol_dense)

plot(sol_nodense)

plot(sol_denseidxs)

This can also checked numerically:

sol_dense(0.6, idxs = [2,4]) # [0.348859, 0.123171]
sol_nodense(0.6, idxs = [2,4]) # [0.352189, 0.124346]
sol_denseidxs(0.6, idxs = [2,4]) # [0.297762, 0.0952841]

sol_dense(0.6, idxs=1) # 0.4768310628469411
sol_nodense(0.6, idxs=1) # 0.4813824589230935
sol_denseidxs(0.6, idxs=1) #0.4478574621174054

I assume, the main problem is that the current interpolation implementation cannot handle different indices for u and k, which results in completely messed up values. In my opinion, the expected and desired output would be sol_denseidxs(0.6, idxs = [2,4]) = sol_dense(0.6, idxs = [2,4]) and sol_denseidxs(0.6, idxs = 1) = sol_nodense(0.6, idxs = 1).

[ChrisRackauckas]

Yes, it would need to index u as well at the dense indices.