Cannot reconstruct parameter-timeseries dependend observables anymore

[hexaeder]

Describe the example

The recently merged #779 broke some things.

Specificially:

• interpolation for discrete values behaves differently wether you interpolate for a vector of time points or for single timepoints
• observables which depend on changing parameters won't be resolved correctly

Minimal Reproducible Example 👇

using ModelingToolkit
using OrdinaryDiffEqTsit5
t = ModelingToolkit.t_nounits
D = ModelingToolkit.D_nounits
@variables x(t) obs(t)
@parameters c(t)
@mtkbuild sys = ODESystem(
    [D(x) ~ c * cos(x), obs ~ c], t, [x], [c]; discrete_events = [1.0 => [c ~ c + 1]])
prob = ODEProblem(sys, [x => 0.0], (0.0, 2pi), [c => 1.0])
sol = solve(prob, Tsit5())

Different behavior in extrapolation between single value and timeseries

sol(1:9999, idxs=sys.c) # works
sol(sol.t[end], idxs=sys.c) # errors

I think it is reasonable to extrapolate discrete time series (just hold the last value forever). That works as expected if one resolves c for an entire vector of time points. Won't do it if you resolve it for a single timepoint though.

ERROR: Solution interpolation cannot extrapolate past the final timepoint. Either solve on a longer timespan or use the local extrapolation from the integrator interface.

Now that I think about it, I am not sure wether this behavior might have been like that forever and I just never noticed though.

Wrong resolution of parameter-dependent observables

plot(sol; idxs = obs)

(you'd expect some staircase behavior)
Since #779, the observable function gets called without a time dependent parameters, because is_discrete_expression(sol, obs) = false and thus this branch is never hit

SciMLBase.jl/src/solutions/ode_solutions.jl

Lines 278 to 287 in 012ff4e

	if is_parameter_timeseries(sol) == Timeseries() && is_discrete_expression(sol, idxs)
	discs::ParameterTimeseriesCollection = RecursiveArrayTools.get_discretes(sol)
	ps = parameter_values(discs)
	for ts_idx in eachindex(discs)
	partition = discs[ts_idx]
	interp_val = ConstantInterpolation(partition.t, partition.u)(
	t, nothing, deriv, nothing, continuity)
	ps = with_updated_parameter_timeseries_values(sol, ps, ts_idx => interp_val)
	end
	end

[AayushSabharwal]

Now that I think about it, I am not sure wether this behavior might have been like that forever and I just never noticed though.

I'm fairly certain this was always the behavior

As for better resolution of observed variable dependencies, that shouldn't be too hard