Cannot index/call solution by pure parameters

[hersle]

The last line in

using ModelingToolkit, OrdinaryDiffEq, Test
using ModelingToolkit: t_nounits as t, D_nounits as D
vars = @variables x(t)
pars = @parameters x0
@named S = ODESystem([D(x) ~ 0.0], t, vars, pars; defaults = [x => x0])
S = structural_simplify(S)
prob = ODEProblem(S, [], (0.0, 1.0), [x0 => 1.0])
sol = solve(prob)
sol[x][begin] # works
sol[x0*2] # works
sol[x0] # fails

fails with

ERROR: Indexing with parameters is deprecated. Use `getp(A, x0)` for parameter indexing.

I think that if indexing by x0*2 works, then so should x0.

[AayushSabharwal]

Yeah we explicitly disallowed indexing with parameters directly, preferring sol.ps[param] instead. I forgot why we did this 😅. However, the check isn't smart enough to figure out 2x0 is a parameter expression and should also be disallowed

[ChrisRackauckas]

It's required for type inference.

[hersle]

Could sol[param] just return a (constant) time series of parameter values, so it is consistent with sol[unknown] and sol[observed]? Then the output size of sol[expr] would be the same for all expr.

[AayushSabharwal]

Could sol[param] just return a (constant) time series of parameter values

That's breaking

Then the output size of sol[expr] would be the same for all expr.

Still not necessarily, since a solution can have multiple timeseries (hybrid systems)

[hersle]

Hmm I see, that makes sense, I suppose this is more subtle than I thought. So only for unknowns/observables there is no ambiguity in which timestream to select? In that case I suppose all pure parameter expressions should indeed be disallowed.

What about when the solution is called/interpolated? Then the timeseries is user-specified, and a similar issue arises:

ts = 0.0:0.1:1.0
sol(ts, idxs=S.x) # works
sol(ts, idxs=S.x*2) # works
sol(ts, idxs=S.x0) # fails
sol(ts, idxs=S.x0*2) # fails

In all cases I would expect the output to have the same size as ts, even if idxs contains only parameters.

[AayushSabharwal]

So only for unknowns/observables there is no ambiguity in which timestream to select?

For unknowns, no. For observed, MTK keeps track of which timeseries each observed variable belongs to. Any observed expression can only be associated with one timeseries. So e.g. if x is a continuous variable, and y is a discrete variable, sol[x * y] is invalid. HOWEVER, interpolation intentionally throws all this out the window and just returns the value at the given time(s). Discrete values are assumed to hold the last value they had.

This discrepancy is due to the fact that there are 3 "types of timeseries".

1. Continuous. The unknowns of the system.
2. Clocked. These variables change on a (possibly periodic, but not necessarily) clock. They are only defined at the time that their clock ticks. Mixing two such variables is not allowed. This is the design philosophy behind sol[x * y] being disallowed.
3. Piecewise continuous. These variables changes at discrete intervals, but are otherwise continuous and hold their last value until it is changed again. SII wasn't written with this in mind, and the interpolation syntax is written to allow this as a workaround.

The eventual goal is for SII to recognize these types of timeseries and handle them appropriately. It's an annoyingly difficult rewrite, though.