Merge branch 'master' into speed

Author: ChrisRackauckas

docs/src/examples/remake.md

@@ -51,11 +51,8 @@ function loss(x, p)
     odeprob = p[1] # ODEProblem stored as parameters to avoid using global variables
     ps = parameter_values(odeprob) # obtain the parameter object from the problem
     ps = replace(Tunable(), ps, x) # create a copy with the values passed to the loss function
-    T = eltype(x)
-    # we also have to convert the `u0` vector
-    u0 = T.(state_values(odeprob))
     # remake the problem, passing in our new parameter object
-    newprob = remake(odeprob; u0 = u0, p = ps)
+    newprob = remake(odeprob; p = ps)
     timesteps = p[2]
     sol = solve(newprob, AutoTsit5(Rosenbrock23()); saveat = timesteps)
     truth = p[3]

src/systems/abstractsystem.jl

@@ -859,6 +859,11 @@ function namespace_defaults(sys)
     for (k, v) in pairs(defs))
 end
 
+function namespace_guesses(sys)
+    guess = guesses(sys)
+    Dict(unknowns(sys, k) => namespace_expr(v, sys) for (k, v) in guess)
+end
+
 function namespace_equations(sys::AbstractSystem, ivs = independent_variables(sys))
     eqs = equations(sys)
     isempty(eqs) && return Equation[]
@@ -968,7 +973,13 @@ function full_parameters(sys::AbstractSystem)
 end
 
 function guesses(sys::AbstractSystem)
-    get_guesses(sys)
+    guess = get_guesses(sys)
+    systems = get_systems(sys)
+    isempty(systems) && return guess
+    for subsys in systems
+        guess = merge(guess, namespace_guesses(subsys))
+    end
+    return guess
 end
 
 # required in `src/connectors.jl:437`

src/systems/diffeqs/odesystem.jl

@@ -354,6 +354,7 @@ function flatten(sys::ODESystem, noeqs = false)
             get_iv(sys),
             unknowns(sys),
             parameters(sys),
+            guesses = guesses(sys),
             observed = observed(sys),
             continuous_events = continuous_events(sys),
             discrete_events = discrete_events(sys),

src/systems/index_cache.jl

@@ -17,7 +17,8 @@ struct ParameterIndex{P, I}
 end
 
 const ParamIndexMap = Dict{Union{Symbol, BasicSymbolic}, Tuple{Int, Int}}
-const UnknownIndexMap = Dict{Union{Symbol, BasicSymbolic}, Union{Int, UnitRange{Int}}}
+const UnknownIndexMap = Dict{
+    Union{Symbol, BasicSymbolic}, Union{Int, UnitRange{Int}, AbstractArray{Int}}}
 
 struct IndexCache
     unknown_idx::UnknownIndexMap
@@ -40,17 +41,32 @@ function IndexCache(sys::AbstractSystem)
         for sym in unks
             usym = unwrap(sym)
             sym_idx = if Symbolics.isarraysymbolic(sym)
-                idx:(idx + length(sym) - 1)
+                reshape(idx:(idx + length(sym) - 1), size(sym))
             else
                 idx
             end
             unk_idxs[usym] = sym_idx
 
-            if hasname(sym)
+            if hasname(sym) && (!istree(sym) || operation(sym) !== getindex)
                 unk_idxs[getname(usym)] = sym_idx
             end
             idx += length(sym)
         end
+        for sym in unks
+            usym = unwrap(sym)
+            istree(sym) && operation(sym) === getindex || continue
+            arrsym = arguments(sym)[1]
+            all(haskey(unk_idxs, arrsym[i]) for i in eachindex(arrsym)) || continue
+
+            idxs = [unk_idxs[arrsym[i]] for i in eachindex(arrsym)]
+            if idxs == idxs[begin]:idxs[end]
+                idxs = reshape(idxs[begin]:idxs[end], size(idxs))
+            end
+            unk_idxs[arrsym] = idxs
+            if hasname(arrsym)
+                unk_idxs[getname(arrsym)] = idxs
+            end
+        end
     end
 
     disc_buffers = Dict{Any, Set{BasicSymbolic}}()
@@ -124,7 +140,7 @@ function IndexCache(sys::AbstractSystem)
             for (j, p) in enumerate(buf)
                 idxs[p] = (i, j)
                 idxs[default_toterm(p)] = (i, j)
-                if hasname(p)
+                if hasname(p) && (!istree(p) || operation(p) !== getindex)
                     idxs[getname(p)] = (i, j)
                     idxs[getname(default_toterm(p))] = (i, j)
                 end

src/systems/parameter_buffer.jl

@@ -402,6 +402,15 @@ function SymbolicIndexingInterface.remake_buffer(sys, oldbuf::MTKParameters, val
     return newbuf
 end
 
+function DiffEqBase.anyeltypedual(
+        p::MTKParameters, ::Type{Val{counter}} = Val{0}) where {counter}
+    DiffEqBase.anyeltypedual(p.tunable)
+end
+function DiffEqBase.anyeltypedual(p::Type{<:MTKParameters{T}},
+        ::Type{Val{counter}} = Val{0}) where {counter} where {T}
+    DiffEqBase.__anyeltypedual(T)
+end
+
 _subarrays(v::AbstractVector) = isempty(v) ? () : (v,)
 _subarrays(v::ArrayPartition) = v.x
 _subarrays(v::Tuple) = v

src/variables.jl

@@ -456,7 +456,7 @@ end
 ## Guess ======================================================================
 struct VariableGuess end
 Symbolics.option_to_metadata_type(::Val{:guess}) = VariableGuess
-getguess(x::Num) = getguess(Symbolics.unwrap(x))
+getguess(x::Union{Num, Symbolics.Arr}) = getguess(Symbolics.unwrap(x))
 
 """
     getguess(x)
@@ -469,8 +469,6 @@ Create variables with a guess like this
 ```
 """
 function getguess(x)
-    p = Symbolics.getparent(x, nothing)
-    p === nothing || (x = p)
     Symbolics.getmetadata(x, VariableGuess, nothing)
 end
 

test/mtkparameters.jl

@@ -2,6 +2,7 @@ using ModelingToolkit
 using ModelingToolkit: t_nounits as t, D_nounits as D, MTKParameters
 using SymbolicIndexingInterface
 using SciMLStructures: SciMLStructures, canonicalize, Tunable, Discrete, Constants
+using OrdinaryDiffEq
 using ForwardDiff
 using JET
 
@@ -206,3 +207,24 @@ end
             portion, ps, ones(length(buffer)))
     end
 end
+
+# Issue#2642
+@parameters α β γ δ
+@variables x(t) y(t)
+eqs = [D(x) ~ (α - β * y) * x
+       D(y) ~ (δ * x - γ) * y]
+@mtkbuild odesys = ODESystem(eqs, t)
+odeprob = ODEProblem(
+    odesys, [x => 1.0, y => 1.0], (0.0, 10.0), [α => 1.5, β => 1.0, γ => 3.0, δ => 1.0])
+tunables, _... = canonicalize(Tunable(), odeprob.p)
+@test tunables isa AbstractVector{Float64}
+
+function loss(x)
+    ps = odeprob.p
+    newps = SciMLStructures.replace(Tunable(), ps, x)
+    newprob = remake(odeprob, p = newps)
+    sol = solve(newprob, Tsit5())
+    return sum(sol)
+end
+
+@test_nowarn ForwardDiff.gradient(loss, collect(tunables))

test/odesystem.jl

@@ -1120,3 +1120,42 @@ tearing_state = TearingState(expand_connections(sys))
 ts_vars = tearing_state.fullvars
 orig_vars = unknowns(sys)
 @test isempty(setdiff(ts_vars, orig_vars))
+
+# Guesses in hierarchical systems
+@variables x(t) y(t)
+@named sys = ODESystem(Equation[], t, [x], []; guesses = [x => 1.0])
+@named outer = ODESystem(
+    [D(y) ~ sys.x + t, 0 ~ t + y - sys.x * y], t, [y], []; systems = [sys])
+@test ModelingToolkit.guesses(outer)[sys.x] == 1.0
+outer = structural_simplify(outer)
+@test ModelingToolkit.get_guesses(outer)[sys.x] == 1.0
+prob = ODEProblem(outer, [outer.y => 2.0], (0.0, 10.0))
+int = init(prob, Rodas4())
+@test int[outer.sys.x] == 1.0
+
+# Ensure indexes of array symbolics are cached appropriately
+@variables x(t)[1:2]
+@named sys = ODESystem(Equation[], t, [x], [])
+sys1 = complete(sys)
+@named sys = ODESystem(Equation[], t, [x...], [])
+sys2 = complete(sys)
+for sys in [sys1, sys2]
+    for (sym, idx) in [(x, 1:2), (x[1], 1), (x[2], 2)]
+        @test is_variable(sys, sym)
+        @test variable_index(sys, sym) == idx
+    end
+end
+
+@variables x(t)[1:2, 1:2]
+@named sys = ODESystem(Equation[], t, [x], [])
+sys1 = complete(sys)
+@named sys = ODESystem(Equation[], t, [x...], [])
+sys2 = complete(sys)
+for sys in [sys1, sys2]
+    @test is_variable(sys, x)
+    @test variable_index(sys, x) == [1 3; 2 4]
+    for i in eachindex(x)
+        @test is_variable(sys, x[i])
+        @test variable_index(sys, x[i]) == variable_index(sys, x)[i]
+    end
+end

test/test_variable_metadata.jl

@@ -16,6 +16,18 @@ using ModelingToolkit
 @test hasguess(y) === true
 @test ModelingToolkit.dump_variable_metadata(y).guess == 0
 
+# Issue#2653
+@variables y[1:3] [guess = ones(3)]
+@test getguess(y) == ones(3)
+@test hasguess(y) === true
+@test ModelingToolkit.dump_variable_metadata(y).guess == ones(3)
+
+for i in 1:3
+    @test getguess(y[i]) == 1.0
+    @test hasguess(y[i]) === true
+    @test ModelingToolkit.dump_variable_metadata(y[i]).guess == 1.0
+end
+
 @variables y
 @test hasguess(y) === false
 @test !haskey(ModelingToolkit.dump_variable_metadata(y), :guess)