fix: detect observed variables and dependent parameters dependent on discrete parameters

[AayushSabharwal]

Close SciML/SciMLBase.jl#786, #3199

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

Add any other context about the problem here.

[ChrisRackauckas]

tests fail

[AayushSabharwal]

It's 1.11 + #3094 causing the failures. Though the docstring error is something I haven't been able to reproduce on 1.11 either, that's just the only thing that changed between CI runs

[hexaeder]

I think this MWE I provided in the issue is not a sufficient testcase. I'd propose the following

    @variables x1(t)=0 x2(t)=0 obs1(t) obs2(t)
    @parameters c1(t)=1 c2=1
    @mtkbuild sys = ODESystem(
        [D(x1) ~ c1,
        D(x2) ~ c2,
        obs1 ~ x1*c1,
        obs2 ~ x2*c2], t; discrete_events = [[1.0] => [c1 ~ 0, c2 ~ 0]])
    prob = ODEProblem(sys, [x1=>0, x2=>0], (0.0, 2))
    sol = solve(prob, Tsit5())

    # tests that should pass (?)
    @test sol([0,2], idxs=c1) == [1.0, 0.0]
    @test sol([0,0.9,1.1,2], idxs=obs1) ≈ [0, 0.9, 0, 0]
    @test sol[obs1] == sol(sol.t, idxs=obs1) # errors because of mixed timeseries

    # the following tests check an observable which depends on a parameter which is not declared time dependent, which is done in the docs on discrete events
    # i don't know how this should be handled. Personally, as a user i'd expect all parameters to be discrete timeseries implicitly
    # Depending on your API design, those failures might be by design and don't need fixing.
    @test sol([0,2], idxs=c2) == [1.0, 0.0]
    @test sol([0,0.9,1.1,2], idxs=obs2) ≈ [0, 0.9, 0, 0]
    @test sol[obs2] == sol(sol.t, idxs=obs2)

[ChrisRackauckas]

Add this?

[AayushSabharwal]

This won't work with the existing SII infrastructure, because parameters changed in callbacks are piecewise continuous (green in the diagram in #3106 (comment)) but SII always assumes discrete variables are clocked (red). It explicitly disallows indexing variables like sol[obs1] because it is a mix of discrete and continuous, and SII doesn't know which timeseries to return. As a workaround, the interpolation syntax still works, so sol(sol.t, idxs=obs1) returns the required values.

[AayushSabharwal]

Supporting this is part of the vision for SII@1.0

[ChrisRackauckas]

okay.

[ChrisRackauckas]

Rebase this for LTS tests?

[AayushSabharwal]

I'll quote something I wrote earlier so that everyone tracking this PR is aware of the state of our discrete variable support:

There are effectively 3 types of variables in a simulation, illustrated in this diagram from the FMI spec. Unknowns are the blue “continuous” variables. The variables in your example are the green type, which FMI terms as “discrete”. The red variables are “clocked” variables. Blue variables are integrated. Green variables change in callbacks and hold their value in between triggers. Red variables only exist at clock ticks. They must be held (hold operator) to use in a continuous domain, or a clockchange to use in a different clock.

The discrete variable/indexing functionality in MTK is specifically geared around the red variables since those were all I was aware of at the time of writing it. As such, they explicitly disallow mixing with other colors, and with other types of clocks. MTK does not have the compiler infrastructure to handle them (it technically can, but the handling is far from robust and has some correctness issues). Time-dependent parameters as in your example are lumped together with the red variables, and thus the indexing you want to do won’t work even with the PR. It will complain that you’re mixing discrete and continuous variables. There is a plan to build the indexing infrastructure to allow for this differentiation in SymbolicIndexingInterface.jl, but unfortunately I have too many other tasks to tackle it. The refactor is very significant and time consuming. It is also something we need to get right because changing it after the fact is very difficult and likely breaking.

You can work around this, however. Your current approach works for red variables. For green variables, mark them as irreducible with [irreducible = true], and provide D(x) ~ 0 as the equation. This should allow you to change the variable in a callback, while having it hold a constant value.

Also worth noting is @hexaeder's example in #3106 (comment). sol[obs1] won't work, but as a workaround sol(sol.t; idxs = obs1) will work since interpolation doesn't have as strict checks.