Merge pull request #442 from TorkelE/ReactionSystem_2_SteadyStateProblem

Enable creation of Steady State Problems from Reaction Systems.

Author: ChrisRackauckas

Project.toml

@@ -51,8 +51,9 @@ julia = "1.2"
 
 [extras]
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["OrdinaryDiffEq", "Test", "StochasticDiffEq"]
+test = ["OrdinaryDiffEq", "SteadyStateDiffEq", "Test", "StochasticDiffEq"]

src/ModelingToolkit.jl

@@ -112,7 +112,7 @@ include("build_function.jl")
 export ODESystem, ODEFunction
 export SDESystem, SDEFunction
 export JumpSystem
-export ODEProblem, SDEProblem, NonlinearProblem, OptimizationProblem
+export ODEProblem, SDEProblem, NonlinearProblem, OptimizationProblem, SteadyStateProblem
 export JumpProblem, DiscreteProblem
 export NonlinearSystem, OptimizationSystem
 export ode_order_lowering

src/systems/diffeqs/abstractodesystem.jl

@@ -223,3 +223,39 @@ function DiffEqBase.ODEProblem{iip}(sys::AbstractODESystem,u0map,tspan,
                         sparse=sparse)
     ODEProblem{iip}(f,u0,tspan,p;kwargs...)
 end
+
+
+### Enables Steady State Problems ###
+function DiffEqBase.SteadyStateProblem(sys::AbstractODESystem, args...; kwargs...)
+    SteadyStateProblem{true}(sys, args...; kwargs...)
+end
+
+"""
+```julia
+function DiffEqBase.SteadyStateProblem(sys::AbstractODESystem,u0map,tspan,
+                                    parammap=DiffEqBase.NullParameters();
+                                    version = nothing, tgrad=false,
+                                    jac = false, Wfact = false,
+                                    checkbounds = false, sparse = false,
+                                    linenumbers = true, parallel=SerialForm(),
+                                    kwargs...) where iip
+```
+Generates an SteadyStateProblem from an ODESystem and allows for automatically
+symbolically calculating numerical enhancements.
+"""
+function DiffEqBase.SteadyStateProblem(sys::AbstractODESystem,u0map,
+                                    parammap=DiffEqBase.NullParameters();
+                                    version = nothing, tgrad=false,
+                                    jac = false, Wfact = false,
+                                    checkbounds = false, sparse = false,
+                                    linenumbers = true, parallel=SerialForm(),
+                                    kwargs...) where iip
+    dvs = states(sys)
+    ps = parameters(sys)
+    u0 = varmap_to_vars(u0map,dvs)
+    p = varmap_to_vars(parammap,ps)
+    f = ODEFunction(sys,dvs,ps,u0;tgrad=tgrad,jac=jac,Wfact=Wfact,checkbounds=checkbounds,
+                        linenumbers=linenumbers,parallel=parallel,
+                        sparse=sparse)
+    SteadyStateProblem(f,u0,p;kwargs...)
+end

src/systems/reaction/reactionsystem.jl

@@ -236,7 +236,7 @@ explicitly on the independent variable (usually time).
 function ismassaction(rx, rs; rxvars = get_variables(rx.rate),
                               haveivdep = any(var -> isequal(rs.iv,convert(Variable,var)), rxvars))
     # if no dependencies must be zero order
-    if isempty(rxvars) 
+    if isempty(rxvars)
         return true
     else
         return !(haveivdep || rx.only_use_rate || any(convert(Variable,rxv) in states(rs) for rxv in rxvars))
@@ -246,7 +246,7 @@ end
 function assemble_jumps(rs)
     eqs = Vector{Union{ConstantRateJump, MassActionJump, VariableRateJump}}()
 
-    for rx in equations(rs)        
+    for rx in equations(rs)
         rxvars    = (rx.rate isa Operation) ? get_variables(rx.rate) : Operation[]
         haveivdep = any(var -> isequal(rs.iv,convert(Variable,var)), rxvars)
         if ismassaction(rx, rs; rxvars=rxvars, haveivdep=haveivdep)
@@ -366,3 +366,10 @@ end
 function DiffEqJump.JumpProblem(rs::ReactionSystem, prob, aggregator, args...; kwargs...)
     return JumpProblem(convert(JumpSystem,rs), prob, aggregator, args...; kwargs...)
 end
+
+# SteadyStateProblem from AbstractReactionNetwork
+function DiffEqBase.SteadyStateProblem(rs::ReactionSystem, u0::Union{AbstractArray, Number}, p, args...; kwargs...)
+    #u0 = typeof(u0) <: Array{<:Pair} ? u0 : Pair.(rs.states,u0)
+    #p = typeof(p) <: Array{<:Pair} ? p : Pair.(rs.ps,p)
+    return SteadyStateProblem(ODEFunction(convert(ODESystem,rs)),u0,p, args...; kwargs...)
+end

test/runtests.jl

@@ -8,6 +8,7 @@ using SafeTestsets, Test
 @safetestset "Build Function Test" begin include("build_function.jl") end
 @safetestset "ODESystem Test" begin include("odesystem.jl") end
 @safetestset "Mass Matrix Test" begin include("mass_matrix.jl") end
+@safetestset "SteadyStateSystem Test" begin include("steadystatesystems.jl") end
 @safetestset "SDESystem Test" begin include("sdesystem.jl") end
 @safetestset "NonlinearSystem Test" begin include("nonlinearsystem.jl") end
 @safetestset "OptimizationSystem Test" begin include("optimizationsystem.jl") end

test/steadystatesystems.jl

@@ -0,0 +1,17 @@
+using ModelingToolkit
+using SteadyStateDiffEq
+using Test
+
+@parameters t r
+@variables x(t)
+@derivatives D'~t
+eqs = [D(x) ~ x^2-r]
+de = ODESystem(eqs)
+
+for factor in [1e-1, 1e0, 1e10], u0_p in [(2.34,2.676),(22.34,1.632),(.3,15.676),(0.3,0.006)]
+    u0 = [x => factor*u0_p[1]]
+    p = [r => factor*u0_p[2]]
+    ss_prob = SteadyStateProblem(de,u0,p)
+    sol = solve(ss_prob,SSRootfind()).u[1]
+    @test abs(sol^2 - factor*u0_p[2]) < 1e-8
+end