Merge branch 'master' into dps/set_operations

Author: ChrisRackauckas

Project.toml

@@ -1,7 +1,7 @@
 name = "ModelingToolkit"
 uuid = "961ee093-0014-501f-94e3-6117800e7a78"
 authors = ["Chris Rackauckas <[email protected]>"]
-version = "4.3.4"
+version = "4.4.1"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"

src/extra_functions.jl

@@ -44,4 +44,34 @@ ModelingToolkit.@register Distributions.Normal(mu,sigma)
 @register ∩(x, y)
 @register ∨(x, y)
 @register ∧(x, y)
-@register ⊆(x, y)
+@register ⊆(x, y)
+
+function LinearAlgebra.det(A::AbstractMatrix{<:Num}; laplace=true)
+    if laplace
+        n = LinearAlgebra.checksquare(A)
+        if n == 1 
+            return A[1, 1]
+        elseif n == 2
+            return A[1, 1] * A[2, 2] - A[1, 2] * A[2, 1]
+        else
+            temp = 0
+            # Laplace expansion along the first column
+            M′ = A[:, 2:end]
+            for i in axes(A, 1)
+                M = M′[(1:n) .!= i, :]
+                d′ = A[i, 1] * det(M)
+                if iseven(i)
+                    temp = iszero(temp) ? d′ : temp - d′
+                else
+                    temp = iszero(temp) ? d′ : temp + d′
+                end
+            end
+        end
+        return temp
+    else
+        if istriu(A) || istril(A)
+            return det(UpperTriangular(A))
+        end
+        return det(lu(A; check = false))
+    end
+end

src/systems/jumps/jumpsystem.jl

@@ -185,17 +185,8 @@ dprob = DiscreteProblem(js, u₀map, tspan, parammap)
 """
 function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Tuple,
                                     parammap=DiffEqBase.NullParameters(); kwargs...)
-
-    (u0map isa AbstractVector) || error("For DiscreteProblems u0map must be an AbstractVector.")
-    u0d = Dict( value(u[1]) => u[2] for u in u0map)
-    u0 = [u0d[u] for u in states(sys)]
-    if parammap != DiffEqBase.NullParameters()
-        (parammap isa AbstractVector) || error("For DiscreteProblems parammap must be an AbstractVector.")
-        pd  = Dict( value(u[1]) => u[2] for u in parammap)
-        p  = [pd[u] for u in parameters(sys)]
-    else
-        p = parammap
-    end
+    u0 = varmap_to_vars(u0map, states(sys))
+    p  = varmap_to_vars(parammap, parameters(sys))
     f  = DiffEqBase.DISCRETE_INPLACE_DEFAULT
     df = DiscreteFunction{true,true}(f, syms=Symbol.(states(sys)))
     DiscreteProblem(df, u0, tspan, p; kwargs...)
@@ -207,7 +198,7 @@ function DiffEqBase.DiscreteProblemExpr(sys::JumpSystem, u0map, tspan,
                                     parammap=DiffEqBase.NullParameters; kwargs...)
 ```
 
-Generates a black DiscreteProblem for a JumpSystem to utilize as its
+Generates a blank DiscreteProblem for a JumpSystem to utilize as its
 solving `prob.prob`. This is used in the case where there are no ODEs
 and no SDEs associated with the system.
 

src/variables.jl

@@ -227,6 +227,22 @@ function varmap_to_vars(varmap::AbstractArray{Pair{T,S}},varlist) where {T,S}
     # Does things like MArray->SArray
     ArrayInterface.restructure(varmap,out)
 end
+
+function varmap_to_vars(varmap::NTuple{N,Pair{T,S}},varlist) where {N,T,S}
+    out = MArray{Tuple{N},S}(undef)
+    for (ivar, ival) in varmap
+        j = findfirst(isequal(ivar),varlist)
+        if isnothing(j)
+            throw(ArgumentError("Value $(ivar) provided in map not found in $(varlist)"))
+        elseif j > length(varmap)
+            throw(ArgumentError("Missing value in $(varmap), need $(varlist)"))
+        end
+        out[j] = ival
+    end
+    out.data
+end
+
 varmap_to_vars(varmap::AbstractArray,varlist) = varmap
+varmap_to_vars(varmap::Tuple,varlist) = varmap
 varmap_to_vars(varmap::DiffEqBase.NullParameters,varlist) = varmap
 varmap_to_vars(varmap::Nothing,varlist) = varmap

test/odesystem.jl

@@ -192,8 +192,14 @@ u0 = [y₁ => 1.0,
 p  = [k₁ => 0.04,
       k₂ => 3e7,
       k₃ => 1e4]
+p2 = (k₁ => 0.04,
+      k₂ => 3e7,
+      k₃ => 1e4)
 tspan = (0.0,100000.0)
 prob1 = ODEProblem(sys,u0,tspan,p)
+prob12 = ODEProblem(sys,u0,tspan,[0.04,3e7,1e4])
+prob13 = ODEProblem(sys,u0,tspan,(0.04,3e7,1e4))
+prob14 = ODEProblem(sys,u0,tspan,p2)
 prob2 = ODEProblem(sys,u0,tspan,p,jac=true)
 prob3 = ODEProblem(sys,u0,tspan,p,jac=true,sparse=true)
 for (prob, atol) in [(prob1, 1e-12), (prob2, 1e-12), (prob3, 1e-12)]

test/operation_overloads.jl

@@ -19,7 +19,7 @@ aa = a; # old a
 
 # test some matrix operations don't throw errors
 X = [0 b c; d e f; g h i]
-@test iszero(simplify(det(X) - ((b * f * g) + (c * d * h) - (b * d * i) - (c * e * g)), polynorm=true))
+@test iszero(simplify(det(X) - ((d * ((b * i) - (c * h))) + (g * ((b * f) - (c * e))))))
 F = lu(X)
 @test F.p == [2, 1, 3]
 R = simplify.(F.L * F.U - X[F.p, :], polynorm=true)
@@ -85,3 +85,27 @@ M \ reshape(b,2,1)
 M = [1 a; 0 2]
 M \ b
 M \ [1, 2]
+
+# test det
+@variables X[1:4,1:4]
+d1 = det(X, laplace=true)
+d2 = det(X, laplace=false)
+_det1 = eval(build_function(d1, X))
+_det2 = eval(build_function(d2, X))
+A = [1 1 1 1
+     1 0 1 1
+     1 1 0 1
+     1 1 1 0]
+@test _det1(A) == -1
+@test _det2(A) == -1
+
+@variables X[1:3,1:3]
+d1 = det(X, laplace=true)
+d2 = det(X, laplace=false)
+_det1 = eval(build_function(d1, X))
+_det2 = eval(build_function(d2, X))
+A = [1 1 1
+     1 0 1
+     1 1 1]
+@test _det1(A) == 0
+@test _det2(A) == 0