Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Composing SDESystems does not work #2085

Closed
hstrey opened this issue Feb 15, 2023 · 13 comments
Closed

Composing SDESystems does not work #2085

hstrey opened this issue Feb 15, 2023 · 13 comments

Comments

@hstrey
Copy link
Contributor

hstrey commented Feb 15, 2023

I am trying to compose SDESystems, and when I do, Julia goes into an infinite loop after I execute compose. Am I doing it the right way? (see below for a slightly better way of attempting to compose SDEs).

using ModelingToolkit

@variables t
D = Differential(t)

function van_der_pol(;name, θ=1.0,ϕ=0.1)
    params  = @parameters θ=θ ϕ=ϕ
    sts = @variables x(t)=1.0 y(t)=1.0 jcn(t)=0.0

    eqs = [D(x) ~ y + jcn,
           D(y) ~ θ*(1-x^2)*y - x]

    noiseeqs = [ϕ,ϕ]

    return SDESystem(eqs,noiseeqs,t,sts,params; name=name)
end

@named VP1 = van_der_pol()
@named VP2 = van_der_pol()

eqs = [VP1.jcn ~ VP2.x,
        VP2.jcn ~ VP1.x]

sys = [VP1,VP2]

@named coupledVP = compose(eqs,sys)
@isaacsas
Copy link
Member

Don't eqs need to be in a system to compose?

@hstrey
Copy link
Contributor Author

hstrey commented Feb 15, 2023

For ODESystems, that is the way we use compose and it works for ODESystems. @ChrisRackauckas told me that this should work also for SDEs, but maybe I am missing a detail or two.

@hstrey
Copy link
Contributor Author

hstrey commented Feb 15, 2023

@isaacsas you are correct. For ODESystems the correct compose is: @named coupledVP = compose(ODESystem(eqs;name=:connected),sys)
but there is no equivalent way for SDESystem.

@isaacsas
Copy link
Member

Yeah, SDESystems are missing both compose and flattening when generating code for DifferentialEquations. Without the latter you might be able to compose systems but you wouldn't be able to actually generate a correct SDEProblem from the system. So there is some work needed to get this together I think.

@hstrey
Copy link
Contributor Author

hstrey commented Feb 15, 2023

I just tried: @named coupledVP = compose(SDESystem(eqs,[],t,[],[];name=:connected),sys)
that creates the correct set of equations, but structural_simplify fails

@isaacsas
Copy link
Member

Does it actually get the correct noise equations too though?

@hstrey
Copy link
Contributor Author

hstrey commented Feb 15, 2023

here is the updated code:

using ModelingToolkit

@variables t
D = Differential(t)

function van_der_pol(;name, θ=1.0,ϕ=0.1)
    params  = @parameters θ=θ ϕ=ϕ
    sts = @variables x(t)=1.0 y(t)=1.0 jcn(t)=0.0

    eqs = [D(x) ~ y + jcn,
           D(y) ~ θ*(1-x^2)*y - x]

    noiseeqs = [ϕ,ϕ]

    return SDESystem(eqs,noiseeqs,t,sts,params; name=name)
end

@named VP1 = van_der_pol()
@named VP2 = van_der_pol()

eqs = [VP1.jcn ~ VP2.x,
        VP2.jcn ~ VP1.x]

sys = [VP1,VP2]

@named coupledVP = compose(SDESystem(eqs,[],t,[],[];name=:connected),sys)
coupledVPs = structural_simplify(coupledVP)

structural_simplify fails with:

ERROR: BoundsError: attempt to access 10-element Vector{Vector{Int64}} at index [11]
Stacktrace:
  [1] getindex
    @ ./array.jl:924 [inlined]
  [2] 𝑠neighbors (repeats 2 times)
    @ ~/.julia/packages/ModelingToolkit/QxF9d/src/bipartite_graph.jl:345 [inlined]
  [3] find_eq_solvables!(state::TearingState{SDESystem}, ieq::Int64, to_rm::Vector{Int64}, coeffs::Vector{Int64}; may_be_zero::Bool, allow_symbolic::Bool, allow_parameter::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/QxF9d/src/structural_transformation/utils.jl:185
  [4] find_eq_solvables!(state::TearingState{SDESystem}, ieq::Int64, to_rm::Vector{Int64}, coeffs::Vector{Int64})
    @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/QxF9d/src/structural_transformation/utils.jl:175
  [5] linear_subsys_adjmat!(state::TearingState{SDESystem}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit.StructuralTransformations ~/.julia/packages/ModelingToolkit/QxF9d/src/structural_transformation/utils.jl:257
  [6] linear_subsys_adjmat!
    @ ~/.julia/packages/ModelingToolkit/QxF9d/src/structural_transformation/utils.jl:245 [inlined]
  [7] alias_eliminate_graph!(state::TearingState{SDESystem}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QxF9d/src/systems/alias_elimination.jl:7
  [8] alias_eliminate_graph!
    @ ~/.julia/packages/ModelingToolkit/QxF9d/src/systems/alias_elimination.jl:6 [inlined]
  [9] alias_elimination!(state::TearingState{SDESystem}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QxF9d/src/systems/alias_elimination.jl:52
 [10] alias_elimination!(state::TearingState{SDESystem})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QxF9d/src/systems/alias_elimination.jl:48
 [11] _structural_simplify!(state::TearingState{SDESystem}, io::Nothing; simplify::Bool, check_consistency::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit.SystemStructures ~/.julia/packages/ModelingToolkit/QxF9d/src/systems/systemstructure.jl:538
 [12] #structural_simplify!#23
    @ ~/.julia/packages/ModelingToolkit/QxF9d/src/systems/systemstructure.jl:523 [inlined]
 [13] structural_simplify(sys::SDESystem, io::Nothing; simplify::Bool, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/QxF9d/src/systems/systems.jl:39
 [14] structural_simplify (repeats 2 times)
    @ ~/.julia/packages/ModelingToolkit/QxF9d/src/systems/systems.jl:19 [inlined]
 [15] top-level scope
    @ ~/Documents/programming/Modeling-SDEs/coupled_SDEs.jl:27

@hstrey
Copy link
Contributor Author

hstrey commented Feb 15, 2023

Does it actually get the correct noise equations too though?

no, the noise equations are missing.

julia> coupledVP.noiseeqs
Any[][1:0]

@isaacsas
Copy link
Member

Right, so it needs a custom compose too in addition to flattening code in generating the needed DiffEqs functions.

@hstrey hstrey changed the title Composing SDESystems goes into an infinite loop Composing SDESystems does not work Feb 16, 2023
@hstrey
Copy link
Contributor Author

hstrey commented Feb 27, 2023

I just noticed this pull request: #1642
@YingboMa, is this still current?

@YingboMa
Copy link
Member

This works:

using ModelingToolkit

@variables t
D = Differential(t)

function van_der_pol(;name, θ=1.0=0.1)
    params  = @parameters θ=θ
    sts = @variables x(t)=1.0 y(t)=1.0 jcn(t)=0.0
    @brownian ϕ

    eqs = [D(x) ~ y + jcn + ϕ,
           D(y) ~ θ*(1-x^2)*y - x + ϕ]

    return System(eqs, t; name=name)
end

@named VP1 = van_der_pol()
@named VP2 = van_der_pol()

eqs = [VP1.jcn ~ VP2.x,
        VP2.jcn ~ VP1.x]

sys = [VP1,VP2]

@named coupledVP = compose(System(eqs,t;name=:connected),sys)
coupledVPs = structural_simplify(coupledVP)

@hstrey
Copy link
Contributor Author

hstrey commented Mar 31, 2023

@YingboMa thanks so much

@arnold-pdev
Copy link

arnold-pdev commented Nov 29, 2023

This works:

using ModelingToolkit

@variables t
D = Differential(t)

function van_der_pol(;name, θ=1.0=0.1)
    params  = @parameters θ=θ
    sts = @variables x(t)=1.0 y(t)=1.0 jcn(t)=0.0
    @brownian ϕ

    eqs = [D(x) ~ y + jcn + ϕ,
           D(y) ~ θ*(1-x^2)*y - x + ϕ]

    return System(eqs, t; name=name)
end

@named VP1 = van_der_pol()
@named VP2 = van_der_pol()

eqs = [VP1.jcn ~ VP2.x,
        VP2.jcn ~ VP1.x]

sys = [VP1,VP2]

@named coupledVP = compose(System(eqs,t;name=:connected),sys)
coupledVPs = structural_simplify(coupledVP)

Hi Yingbo-

Could your provide an example of a construction of a (ODE? SDE?)Problem using this System? I am running into the issue that the number of equations and variables don't match.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants