Merge pull request #2391 from AayushSabharwal/as/sii-multivariate

refactor: handle NullParameters, fix independent_variable_symbols for multivariate systems

Author: ChrisRackauckas

src/systems/abstractsystem.jl

@@ -185,75 +185,77 @@ end
 #Treat the result as a vector of symbols always
 function SymbolicIndexingInterface.is_variable(sys::AbstractSystem, sym)
     if unwrap(sym) isa Int    # [x, 1] coerces 1 to a Num
-        return unwrap(sym) in 1:length(unknown_states(sys))
+        return unwrap(sym) in 1:length(variable_symbols(sys))
     end
-    return any(isequal(sym), unknown_states(sys)) ||
+    return any(isequal(sym), variable_symbols(sys)) ||
            hasname(sym) && is_variable(sys, getname(sym))
 end
 
 function SymbolicIndexingInterface.is_variable(sys::AbstractSystem, sym::Symbol)
-    return any(isequal(sym), getname.(unknown_states(sys))) ||
+    return any(isequal(sym), getname.(variable_symbols(sys))) ||
            count('₊', string(sym)) == 1 &&
-           count(isequal(sym), Symbol.(sys.name, :₊, getname.(unknown_states(sys)))) == 1
+           count(isequal(sym), Symbol.(sys.name, :₊, getname.(variable_symbols(sys)))) == 1
 end
 
 function SymbolicIndexingInterface.variable_index(sys::AbstractSystem, sym)
     if unwrap(sym) isa Int
         return unwrap(sym)
     end
-    idx = findfirst(isequal(sym), unknown_states(sys))
+    idx = findfirst(isequal(sym), variable_symbols(sys))
     if idx === nothing && hasname(sym)
         idx = variable_index(sys, getname(sym))
     end
     return idx
 end
 
 function SymbolicIndexingInterface.variable_index(sys::AbstractSystem, sym::Symbol)
-    idx = findfirst(isequal(sym), getname.(unknown_states(sys)))
+    idx = findfirst(isequal(sym), getname.(variable_symbols(sys)))
     if idx !== nothing
         return idx
     elseif count('₊', string(sym)) == 1
-        return findfirst(isequal(sym), Symbol.(sys.name, :₊, getname.(unknown_states(sys))))
+        return findfirst(isequal(sym), Symbol.(sys.name, :₊, getname.(variable_symbols(sys))))
     end
     return nothing
 end
 
+SymbolicIndexingInterface.variable_symbols(sys::AbstractMultivariateSystem) = sys.dvs
+
 function SymbolicIndexingInterface.variable_symbols(sys::AbstractSystem)
     return unknown_states(sys)
 end
 
 function SymbolicIndexingInterface.is_parameter(sys::AbstractSystem, sym)
     if unwrap(sym) isa Int
-        return unwrap(sym) in 1:length(parameters(sys))
+        return unwrap(sym) in 1:length(parameter_symbols(sys))
     end
 
-    return any(isequal(sym), parameters(sys)) ||
+    return any(isequal(sym), parameter_symbols(sys)) ||
            hasname(sym) && is_parameter(sys, getname(sym))
 end
 
 function SymbolicIndexingInterface.is_parameter(sys::AbstractSystem, sym::Symbol)
-    return any(isequal(sym), getname.(parameters(sys))) ||
+    return any(isequal(sym), getname.(parameter_symbols(sys))) ||
            count('₊', string(sym)) == 1 &&
-           count(isequal(sym), Symbol.(sys.name, :₊, getname.(parameters(sys)))) == 1
+           count(isequal(sym), Symbol.(sys.name, :₊, getname.(parameter_symbols(sys)))) == 1
 end
 
 function SymbolicIndexingInterface.parameter_index(sys::AbstractSystem, sym)
     if unwrap(sym) isa Int
         return unwrap(sym)
     end
-    idx = findfirst(isequal(sym), parameters(sys))
+    idx = findfirst(isequal(sym), parameter_symbols(sys))
     if idx === nothing && hasname(sym)
         idx = parameter_index(sys, getname(sym))
     end
     return idx
 end
 
 function SymbolicIndexingInterface.parameter_index(sys::AbstractSystem, sym::Symbol)
-    idx = findfirst(isequal(sym), getname.(parameters(sys)))
+    idx = findfirst(isequal(sym), getname.(parameter_symbols(sys)))
     if idx !== nothing
         return idx
     elseif count('₊', string(sym)) == 1
-        return findfirst(isequal(sym), Symbol.(sys.name, :₊, getname.(parameters(sys))))
+        return findfirst(isequal(sym), Symbol.(sys.name, :₊, getname.(parameter_symbols(sys))))
     end
     return nothing
 end
@@ -263,7 +265,7 @@ function SymbolicIndexingInterface.parameter_symbols(sys::AbstractSystem)
 end
 
 function SymbolicIndexingInterface.is_independent_variable(sys::AbstractSystem, sym)
-    return any(isequal(sym), independent_variables(sys))
+    return any(isequal(sym), independent_variable_symbols(sys))
 end
 
 function SymbolicIndexingInterface.is_independent_variable(sys::AbstractSystem, sym::Symbol)
@@ -650,6 +652,12 @@ end
 
 function parameters(sys::AbstractSystem)
     ps = get_ps(sys)
+    if ps == SciMLBase.NullParameters()
+        return []
+    end
+    if eltype(ps) <: Pair
+        ps = first.(ps)
+    end
     systems = get_systems(sys)
     unique(isempty(systems) ? ps : [ps; reduce(vcat, namespace_parameters.(systems))])
 end

test/symbolic_indexing_interface.jl

@@ -0,0 +1,77 @@
+using ModelingToolkit, SymbolicIndexingInterface, SciMLBase
+
+@parameters t a b
+@variables x(t) y(t)
+D = Differential(t)
+eqs = [D(x) ~ a * y + t, D(y) ~ b * t]
+@named odesys = ODESystem(eqs, t, [x, y], [a, b])
+
+@test all(is_variable.((odesys,), [x, y, 1, 2, :x, :y]))
+@test all(.!is_variable.((odesys,), [a, b, t, 3, 0, :a, :b]))
+@test variable_index.((odesys,), [x, y, a, b, t, 1, 2, :x, :y, :a, :b]) == [1, 2, nothing, nothing, nothing, 1, 2, 1, 2, nothing, nothing]
+@test isequal(variable_symbols(odesys), [x, y])
+@test all(is_parameter.((odesys,), [a, b, 1, 2, :a, :b]))
+@test all(.!is_parameter.((odesys,), [x, y, t, 3, 0, :x, :y]))
+@test parameter_index.((odesys,), [x, y, a, b, t, 1, 2, :x, :y, :a, :b]) == [nothing, nothing, 1, 2, nothing, 1, 2, nothing, nothing, 1, 2]
+@test isequal(parameter_symbols(odesys), [a, b])
+@test all(is_independent_variable.((odesys,), [t, :t]))
+@test all(.!is_independent_variable.((odesys,), [x, y, a, :x, :y, :a]))
+@test isequal(independent_variable_symbols(odesys), [t])
+@test is_time_dependent(odesys)
+@test constant_structure(odesys)
+
+@variables x y z
+@parameters σ ρ β
+
+eqs = [0 ~ σ*(y-x),
+       0 ~ x*(ρ-z)-y,
+       0 ~ x*y - β*z]
+@named ns = NonlinearSystem(eqs, [x,y,z],[σ,ρ,β])
+
+@test !is_time_dependent(ns)
+
+@parameters x
+@variables t u(..)
+Dxx = Differential(x)^2
+Dtt = Differential(t)^2
+Dt = Differential(t)
+
+#2D PDE
+C=1
+eq  = Dtt(u(t,x)) ~ C^2*Dxx(u(t,x))
+
+# Initial and boundary conditions
+bcs = [u(t,0) ~ 0.,# for all t > 0
+       u(t,1) ~ 0.,# for all t > 0
+       u(0,x) ~ x*(1. - x), #for all 0 < x < 1
+       Dt(u(0,x)) ~ 0. ] #for all  0 < x < 1]
+
+# Space and time domains
+domains = [t ∈ (0.0,1.0),
+           x ∈ (0.0,1.0)]
+
+@named pde_system = PDESystem(eq,bcs,domains,[t,x],[u])
+
+@test pde_system.ps == SciMLBase.NullParameters()
+@test parameter_symbols(pde_system) == []
+
+@parameters t x
+@constants h = 1
+@variables u(..)
+Dt = Differential(t)
+Dxx = Differential(x)^2
+eq = Dt(u(t, x)) ~ h * Dxx(u(t, x))
+bcs = [u(0, x) ~ -h * x * (x - 1) * sin(x),
+    u(t, 0) ~ 0, u(t, 1) ~ 0]
+
+domains = [t ∈ (0.0, 1.0),
+    x ∈ (0.0, 1.0)]
+
+analytic = [u(t, x) ~ -h * x * (x - 1) * sin(x) * exp(-2 * h * t)]
+analytic_function = (ps, t, x) -> -ps[1] * x * (x - 1) * sin(x) * exp(-2 * ps[1] * t)
+
+@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u], [h => 1], analytic = analytic)
+
+@test isequal(pdesys.ps, [h => 1])
+@test isequal(parameter_symbols(pdesys), [h])
+@test isequal(parameters(pdesys), [h])