Clean-up comments & tests.

Author: Lucas Morton

src/systems/validation.jl

@@ -5,19 +5,33 @@ struct ValidationError <: Exception
     message::String
 end
 
+"Throw exception on invalid unit types, otherwise return argument."
 function screen_unit(result)
     result isa Unitful.Unitlike || throw(ValidationError("Unit must be a subtype of Unitful.Unitlike, not $(typeof(result))."))
     result isa Unitful.ScalarUnits || throw(ValidationError("Non-scalar units such as $result are not supported. Use a scalar unit instead."))
     result == u"°" && throw(ValidationError("Degrees are not supported. Use radians instead."))
     result
 end
+
+"""Test unit equivalence.
+
+Example of implemented behavior:
+```julia
+using ModelingToolkit, Unitful
+MT = ModelingToolkit
+@parameters γ P [unit = u"MW"] E [unit = u"kJ"] τ [unit = u"ms"]
+@test MT.equivalent(u"MW" ,u"kJ/ms") # Understands prefixes
+@test !MT.equivalent(u"m", u"cm") # Units must be same magnitude
+@test MT.equivalent(MT.get_unit(P^γ), MT.get_unit((E/τ)^γ)) # Handles symbolic exponents
+```
+"""
 equivalent(x,y) = isequal(1*x,1*y)
 unitless = Unitful.unit(1)
 
 #For dispatching get_unit
 Literal = Union{Sym,Symbolics.ArrayOp,Symbolics.Arr,Symbolics.CallWithMetadata}
 Conditional = Union{typeof(ifelse),typeof(IfElse.ifelse)}
-Comparison = Union{typeof(Base.:>), typeof(Base.:<), typeof(==)}
+Comparison = Union{typeof.([==, !=, ≠, <, <=, ≤, >, >=, ≥])...}
 
 "Find the unit of a symbolic item."
 get_unit(x::Real) = unitless
@@ -26,9 +40,9 @@ get_unit(x::AbstractArray) = map(get_unit,x)
 get_unit(x::Num) = get_unit(value(x))
 get_unit(x::Literal) = screen_unit(getmetadata(x,VariableUnit, unitless))
 get_unit(op::Differential, args) = get_unit(args[1]) / get_unit(op.x)
-get_unit(op::Difference, args) =   get_unit(args[1]) / get_unit(op.t) #why are these not identical?!?
+get_unit(op::Difference, args) =   get_unit(args[1]) / get_unit(op.t)
 get_unit(op::typeof(getindex),args) = get_unit(args[1]) 
-function get_unit(op,args) #Fallback
+function get_unit(op,args) # Fallback
     result = op(1 .* get_unit.(args)...)
     try 
         unit(result)
@@ -87,24 +101,22 @@ end
 
 function get_unit(op::Comparison, args)
     terms = get_unit.(args)
-    equivalent(terms[1], terms[2]) || throw(ValidationError(", in comparison $x, units [$(terms[1])] and [$(terms[2])] do not match."))
+    equivalent(terms[1], terms[2]) || throw(ValidationError(", in comparison $op, units [$(terms[1])] and [$(terms[2])] do not match."))
     return unitless
 end
 
 function get_unit(x::Symbolic) 
     if SymbolicUtils.istree(x)
         op = operation(x)
-        if op isa Sym # Not a real function call, just a dependent variable. Unit is on the Sym.
-            return screen_unit(getmetadata(x, VariableUnit, unitless)) 
-        elseif op isa Term && !(operation(op) isa Term) #
-            gp = getmetadata(x,Symbolics.GetindexParent,nothing)
+        if op isa Sym || (op isa Term && operation(op) isa Term) # Dependent variables, not function calls
+            return screen_unit(getmetadata(x, VariableUnit, unitless)) # Like x(t) or x[i]
+        elseif op isa Term && !(operation(op) isa Term) 
+            gp = getmetadata(x,Symbolics.GetindexParent,nothing) # Like x[1](t)
             return screen_unit(getmetadata(gp, VariableUnit, unitless))
-        elseif op isa Term
-            return screen_unit(getmetadata(x, VariableUnit, unitless))
-        end        
+        end  # Actual function calls:
         args = arguments(x)
         return get_unit(op, args)
-    else #This function should only be reached by Terms, for which `istree` is true
+    else # This function should only be reached by Terms, for which `istree` is true
         throw(ArgumentError("Unsupported value $x."))
     end
 end
@@ -119,7 +131,7 @@ function safe_get_unit(term, info)
             @warn("$info: $(err.x) and $(err.y) are not dimensionally compatible.")
         elseif err isa ValidationError
             @warn(info*err.message)
-        elseif err isa MethodError #Warning: Unable to get unit for operation x[1] with arguments SymbolicUtils.Sym{Real, Base.ImmutableDict{DataType, Any}}[t].
+        elseif err isa MethodError
             @warn("$info: no method matching $(err.f) for arguments $(typeof.(err.args)).")
         else
             rethrow()

test/units.jl

@@ -1,65 +1,59 @@
 using ModelingToolkit, Unitful, OrdinaryDiffEq, DiffEqJump, IfElse
 using Test
 MT = ModelingToolkit
-@parameters τ [unit = u"ms"]
+@parameters τ [unit = u"ms"] γ
 @variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
 D = Differential(t)
 
+#This is how equivalent works:
+@test MT.equivalent(u"MW" ,u"kJ/ms")
+@test !MT.equivalent(u"m", u"cm")
+@test MT.equivalent(MT.get_unit(P^γ), MT.get_unit((E/τ)^γ))
+
+# Basic access
 @test MT.get_unit(t) == u"ms"
 @test MT.get_unit(E) == u"kJ"
 @test MT.get_unit(τ) == u"ms"
+@test MT.get_unit(γ) == MT.unitless
+@test MT.get_unit(0.5) == MT.unitless
 
+# Prohibited unit types
 @parameters β [unit = u"°"] α [unit = u"°C"] γ [unit = 1u"s"]
 @test_throws MT.ValidationError MT.get_unit(β)
 @test_throws MT.ValidationError MT.get_unit(α)
 @test_throws MT.ValidationError MT.get_unit(γ)
 
-unitless = Unitful.unit(1)
-@test MT.get_unit(0.5) == unitless
-@test MT.get_unit(t) == u"ms"
-@test MT.get_unit(P) == u"MW"
-@test MT.get_unit(τ) == u"ms"
-
+# Non-trivial equivalence & operators
 @test MT.get_unit(τ^-1) == u"ms^-1"
 @test MT.equivalent(MT.get_unit(D(E)),u"MW")
 @test MT.equivalent(MT.get_unit(E/τ), u"MW")
 @test MT.get_unit(2*P) == u"MW"
-@test MT.get_unit(t/τ) == unitless
+@test MT.get_unit(t/τ) == MT.unitless
 @test MT.equivalent(MT.get_unit(P - E/τ),u"MW")
 @test MT.equivalent(MT.get_unit(D(D(E))),u"MW/ms")
-
-@test MT.get_unit(1.0^(t/τ)) == unitless
-@test MT.get_unit(exp(t/τ)) == unitless
-@test MT.get_unit(sin(t/τ)) == unitless
-@test MT.get_unit(sin(1u"rad")) == unitless
+@test MT.get_unit(IfElse.ifelse(t>t,P,E/τ)) == u"MW"
+@test MT.get_unit(1.0^(t/τ)) == MT.unitless
+@test MT.get_unit(exp(t/τ)) == MT.unitless
+@test MT.get_unit(sin(t/τ)) == MT.unitless
+@test MT.get_unit(sin(1u"rad")) == MT.unitless
 @test MT.get_unit(t^2) == u"ms^2"
 
-@test !MT.validate(E^1.5 ~ E^(t/τ))
-@test MT.validate(E^(t/τ) ~ E^(t/τ))
-
 eqs = [D(E) ~ P - E/τ
         0 ~ P]
 @test MT.validate(eqs)
 @named sys = ODESystem(eqs)
-@named sys = ODESystem(eqs, t, [P, E], [τ])
 
 @test !MT.validate(D(D(E)) ~ P)
 @test !MT.validate(0 ~ P + E*τ)
 
-#Unit-free
-@variables x y z u
-@parameters σ ρ β
-eqs = [0 ~ σ*(y - x)]
-@test MT.validate(eqs)
-
-##Array variables
+# Array variables
 @variables t [unit = u"s"] x[1:3](t) [unit = u"m"]
 @parameters v[1:3] = [1,2,3] [unit = u"m/s"]
 D = Differential(t)
 eqs = D.(x) .~ v
 ODESystem(eqs,name=:sys)
 
-#Difference equation with units
+# Difference equation
 @parameters t [unit = u"s"] a [unit = u"s"^-1]
 @variables x(t) [unit = u"kg"]
 δ = Differential(t)
@@ -69,18 +63,6 @@ eqs = [
 ]
 de = ODESystem(eqs, t, [x], [a],name=:sys)
 
-
-@parameters t
-@variables y[1:3](t)
-@parameters k[1:3]
-D = Differential(t)
-
-eqs = [D(y[1]) ~ -k[1]*y[1] + k[3]*y[2]*y[3],
-       D(y[2]) ~  k[1]*y[1] - k[3]*y[2]*y[3] - k[2]*y[2]^2,
-       0 ~  y[1] + y[2] + y[3] - 1]
-
-@named sys = ODESystem(eqs,t,y,k)
-
 # Nonlinear system
 @parameters a [unit = u"kg"^-1]
 @variables x [unit = u"kg"]
@@ -99,16 +81,18 @@ eqs = [D(E) ~ P - E/τ
 noiseeqs = [0.1u"MW",
             0.1u"MW"]
 @named sys = SDESystem(eqs, noiseeqs, t, [P, E], [τ, Q])
+
 # With noise matrix
 noiseeqs = [0.1u"MW" 0.1u"MW"
             0.1u"MW" 0.1u"MW"]
 @named sys = SDESystem(eqs,noiseeqs, t, [P, E], [τ, Q])
 
+# Invalid noise matrix 
 noiseeqs = [0.1u"MW" 0.1u"MW"
             0.1u"MW" 0.1u"s"]
 @test !MT.validate(eqs,noiseeqs)
 
-#Test non-trivial simplifications
+# Non-trivial simplifications
 @variables t [unit = u"s"] V(t) [unit = u"m"^3] L(t) [unit = u"m"]
 @parameters v [unit = u"m/s"] r [unit =u"m"^3/u"s"]
 D = Differential(t)
@@ -166,12 +150,3 @@ maj1 = MassActionJump(2.0, [0 => 1], [S => 1])
 maj2 = MassActionJump(γ, [S => 1], [S => -1])
 @named js4  = JumpSystem([maj1, maj2], t, [S], [β, γ])
 
-#Test comparisons
-@parameters t
-vars = @variables x(t)
-D = Differential(t)
-eqs = [
-    D(x) ~ IfElse.ifelse(t>0.1,2,1)
-]
-@named sys = ODESystem(eqs, t, vars, [])
-