SciML / ModelingToolkit.jl

An acausal modeling framework for automatically parallelized scientific machine learning (SciML) in Julia. A computer algebra system for integrated symbolics for physics-informed machine learning and automated transformations of differential equations
https://mtk.sciml.ai/dev/
Other
1.41k stars 204 forks source link

Regression: Problem with indexing of JumpProblem #2838

Open TorkelE opened 3 months ago

TorkelE commented 3 months ago

Also one case which did not work still does.

using ModelingToolkit, JumpProcesses

@parameters p d
@variables t X(t)
rate1   = p
rate2   = X*d
affect1 = [X ~ X + 1]
affect2 = [X ~ X - 1]
j1 = ConstantRateJump(rate1, affect1)
j2 = ConstantRateJump(rate2, affect2)

# Works.
@mtkbuild js = JumpSystem([j1, j2], t, [X], [p,d])
dprob = DiscreteProblem(js, [X => 15], (0.0, 10.), [p => 2.0, d => 0.5])
jprob = JumpProblem(js, dprob, Direct())
sol = solve(jprob, SSAStepper())

jprob[X]
jprob[[X]] # Error (regression)
jprob[(X,)] # Error, did error previously.
AayushSabharwal commented 3 months ago

Yeah, JumpProcesses should not have a getindex method for JumpProblem. Those tests are marked as broken

TorkelE commented 3 months ago

Glad to hear that it is under control. I am a bit confused though, was this feature intentionally withdrawn for now? Maybe more importantly, it is intended to re introduce this type of indexing?

AayushSabharwal commented 3 months ago

It wasn't intentionally withdrawn, just that it was in a way inevitable. It's a problem sort of unique to Julia. There are two ways out of this:

  1. JumpProcesses removes the getindex method
  2. The getindex method in SciMLBase is redefined to apply to apply to all of these types except AbstractJumpProblem:
    
    julia> subtypes(SciMLBase.AbstractSciMLProblem)
    5-element Vector{Any}:
    SciMLBase.AbstractDEProblem
    SciMLBase.AbstractEnsembleProblem
    SciMLBase.AbstractIntegralProblem
    SciMLBase.AbstractLinearProblem
    SciMLBase.AbstractOptimizationProblem

julia> subtypes(SciMLBase.AbstractDEProblem) 9-element Vector{Any}: SciMLBase.AbstractDAEProblem SciMLBase.AbstractDDEProblem SciMLBase.AbstractJumpProblem SciMLBase.AbstractNoiseProblem SciMLBase.AbstractNonlinearProblem SciMLBase.AbstractODEProblem SciMLBase.AbstractPDEProblem SciMLBase.AbstractRODEProblem SciMLBase.AbstractSDDEProblem

AayushSabharwal commented 3 months ago

This sort of thing crops up a lot in Julia, where you want a subtype to just wrap another type and yet forwarding the method can lead to unintentional ambiguities down the line, despite the fact that no semantic guarantees were broken

TorkelE commented 3 months ago

got it, makes sense. thanks again for helping improve the indexing test :)

ChrisRackauckas commented 3 months ago

The getindex method in SciMLBase is redefined to apply to apply to all of these types except AbstractJumpProblem:

We might want to do something of the sort. In particular, some AbstractConcreteSciMLProblem or something, and then better define the interface on that, where everything has a u0, etc. This would be a good part of making our interface definitions more exact.

AayushSabharwal commented 3 months ago

Changing the type hierarchy is breaking, though?

ChrisRackauckas commented 3 months ago

No, we'd keep AbstractSciMLProblem the same, we'd just add something below it that makes more assumptions and update the dispatches to those that require those assumptions. AbstractConcreteSciMLProblem <: AbstractSciMLProblem is assumed.

AayushSabharwal commented 3 months ago

If AbstractDEProblem <: AbstractConcreteSciMLProblem <: AbstractSciMLProblem then AbstractJumpProblem would have to stop subtyping AbstractDEProblem. If AbstractConcreteSciMLProblem <: AbstractDEProblem <: AbstractSciMLProblem then it doesn't really make sense because there are concrete problems which are not DE problems. The concept of concreteness and the categorization we have here are orthogonal

ChrisRackauckas commented 3 months ago

If AbstractDEProblem <: AbstractConcreteSciMLProblem <: AbstractSciMLProblem then AbstractJumpProblem would have to stop subtyping AbstractDEProblem.

Hmm... we'd need to hash out the whole plan, but I think it would be good to spend the time to do this and then actually define what these abstract types stand for.