Add tests & fixes for time derivatives

Author: serenity4

src/analysis/lattice.jl

@@ -139,11 +139,13 @@ struct Incidence
         if is_non_incidence_type(type)
             throw(DomainError(type, "Invalid type for Incidence"))
         end
-        if !isa(row, IncidenceVector)
+        if !isa(row, SparseVector)
             vec, row = row, _zero_row()
             for (i, val) in enumerate(vec)
                 row[i] = val
             end
+        else
+            row = convert(IncidenceVector, row)
         end
         time = row[1]
         if in(time, (linear_time_dependent, linear_time_and_state_dependent))

src/analysis/refiner.jl

@@ -79,8 +79,8 @@ end
 function structural_inc_ddt(var_to_diff::DiffGraph, varclassification::Union{Vector{VarEqClassification}, Nothing}, varkinds::Union{Vector{Intrinsics.VarKind}, Nothing}, inc::Union{Incidence, Const})
     isa(inc, Const) && return Const(zero(inc.val))
     r = _zero_row()
-    function get_or_make_diff(v_offset::Int)
-        v = v_offset - 1
+    function get_or_make_diff(i::Int)
+        v = i - 1
         var_to_diff[v] !== nothing && return var_to_diff[v] + 1
         dv = add_vertex!(var_to_diff)
         if varclassification !== nothing
@@ -94,21 +94,20 @@ function structural_inc_ddt(var_to_diff::DiffGraph, varclassification::Union{Vec
     end
     base = isa(inc.typ, Const) ? Const(zero(inc.typ.val)) : inc.typ
     indices = rowvals(inc.row)
-    for (v_offset, coeff) in zip(indices, nonzeros(inc.row))
-        if v_offset == 1
-            # t
+    for (i, coeff) in zip(indices, nonzeros(inc.row))
+        if i == 1 # time
             if isa(coeff, Float64) # known constant coefficient
-                # Do not set r[v_offset]; d/dt t = 1
+                # Do not set r[i]; d/dt t = 1
                 if isa(base, Const)
                     base = Const(base.val + coeff)
                 end
             elseif coeff.nonlinear
-                r[v_offset] = nonlinear
+                r[i] = nonlinear
             else
                 @assert !coeff.time_dependent # should be nonlinear if time-dependent
                 if coeff.state_dependent # e.g. u₁ * t
                     # State dependence will not be eliminated because of the chain rule.
-                    r[v_offset] = coeff
+                    r[i] = coeff
                 else # unknown constant coefficient
                     if isa(base, Const)
                         # We are adding an unknown but constant value to the
@@ -121,16 +120,30 @@ function structural_inc_ddt(var_to_diff::DiffGraph, varclassification::Union{Vec
         end
         if isa(coeff, Float64)
             # Linear with a known constant coefficient, just add to the derivative
-            r[get_or_make_diff(v_offset)] += coeff
+            r[get_or_make_diff(i)] += coeff
         elseif !coeff.state_dependent && !coeff.time_dependent
             # Linear with an unknown constant coefficient.
-            r[get_or_make_diff(v_offset)] = coeff
+            r[get_or_make_diff(i)] = coeff
         elseif coeff.nonlinear
-            r[v_offset] = nonlinear
-            r[get_or_make_diff(v_offset)] = nonlinear
+            r[i] = nonlinear
+            r[get_or_make_diff(i)] = nonlinear
+            # derivative may yield constant terms (only if time-dependent,
+            # but we conservatively assume nonlinear not to be time-independent)
+            base = widenconst(base)
         else # time- or state-dependent linear coefficient
-            r[v_offset] = coeff
-            r[get_or_make_diff(v_offset)] = coeff
+            if !coeff.state_dependent && coeff.time_dependent && length(nonzeros(inc.row)) == 2 && inc.row[1] ≠ nonlinear
+                # `f(∝t, ∝ₜuᵢ)` is of the form `(a + bt)(c + duᵢ)`, we can simplify to `(a + bt)t∂ₜuᵢ + b(c + duᵢ)`,
+                # yielding `∝ₜ∂uᵢ + bc + ∝uᵢ`. (note that in this case we'll also need to widen `base`, as done further below)
+                r[i] = linear
+                r[get_or_make_diff(i)] = linear_time_dependent
+            else
+                r[i] = coeff
+                r[get_or_make_diff(i)] = coeff
+            end
+            if coeff.time_dependent
+                # The product rule with a time factor may yield constant terms.
+                base = widenconst(base)
+            end
         end
     end
     return Incidence(base, r)

test/incidence.jl

@@ -23,7 +23,7 @@ postwalk(f, ex) = walk(ex, x -> postwalk(f, x), f)
       x + u₃
     end
     @infer_incidence u₁ + u₂ true    # a second argument of `true` makes it return the IR, not the Incidence
-    @infer_incidence exp(u₂) + u₁    # will create u₂ as the second continuous variable
+    @infer_incidence 1/u₂ + u₁       # will create u₂ as the second continuous variable
 
 Return the `Incidence` object after inferring the structure of the provided code,
 substituting any variables starting by 'u'. Variables are created as `continuous()`
@@ -174,6 +174,11 @@ dependencies(row) = sort(rowvals(row) .=> nonzeros(row), by = first)
     @test dependencies(incidence.row) == [2 => 1]
     @test incidence == incidence"u₁"
 
+    incidence = @infer_incidence 3.0 *ᵢ u₁
+    @test incidence.typ === Const(0.0)
+    @test dependencies(incidence.row) == [2 => 3.0]
+    @test incidence == incidence"3.0u₁"
+
     incidence = @infer_incidence u₁ +ᵢ u₂
     @test incidence.typ === Const(0.0)
     @test dependencies(incidence.row) == [2 => 1, 3 => 1]
@@ -274,26 +279,66 @@ dependencies(row) = sort(rowvals(row) .=> nonzeros(row), by = first)
     @test dependencies(incidence.row) == [2 => nonlinear, 3 => linear_state_dependent, 4 => linear_state_dependent]
     @test incidence == incidence"f(u₁, ∝ₛu₂, ∝ₛu₃)"
 
-    incidence = @infer_incidence exp(u₁)
+    incidence = @infer_incidence 1/u₁
     @test dependencies(incidence.row) == [2 => nonlinear]
     @test incidence == incidence"a + f(u₁)"
 
-    incidence = @infer_incidence t * exp(u₁)
+    incidence = @infer_incidence t * (1/u₁)
     @test dependencies(incidence.row) == [1 => linear_state_dependent, 2 => nonlinear]
     @test incidence == incidence"a + f(∝ₛt, u₁)"
 
-    incidence = @infer_incidence u₁ * exp(t)
+    incidence = @infer_incidence u₁ * (1/t)
     @test dependencies(incidence.row) == [1 => nonlinear, 2 => linear_time_dependent]
     @test incidence == incidence"a + f(t, ∝ₜu₁)"
 
-    incidence = @infer_incidence u₁ * exp(t + u₂)
+    incidence = @infer_incidence u₁ * (1/(t + u₂))
     @test dependencies(incidence.row) == [1 => nonlinear, 2 => linear_time_and_state_dependent, 3 => nonlinear]
     @test incidence == incidence"a + f(t, ∝ₜₛu₁, u₂)"
 
-    incidence = @infer_incidence atan(u₁, u₂)
+    incidence = @infer_incidence 1/(u₁ * u₂)
     @test dependencies(incidence.row) == [2 => nonlinear, 3 => nonlinear]
     @test incidence == incidence"a + f(u₁, u₂)"
   end
-end;
+
+  @testset "Time derivatives" begin
+    incidence = @infer_incidence ddt(3.0 *ᵢ t)
+    @test incidence == incidence"3.0"
+
+    incidence = @infer_incidence ddt(3.0 *ᵢ t +ᵢ 5.0)
+    @test incidence == incidence"3.0"
+
+    incidence = @infer_incidence ddt(3.0 * t)
+    @test incidence == incidence"a"
+
+    incidence = @infer_incidence ddt(u₁)
+    @test incidence == incidence"u₂"
+
+    incidence = @infer_incidence ddt(1.0 +ᵢ u₁)
+    @test incidence == incidence"u₂"
+
+    incidence = @infer_incidence u₁ +ᵢ ddt(u₁)
+    @test incidence == incidence"u₁ + u₂"
+
+    incidence = @infer_incidence ddt(u₁ *ᵢ u₂)
+    @test incidence == incidence"f(∝ₛu₁, ∝ₛu₂, ∝ₛu₃, ∝ₛu₄)"
+
+    incidence = @infer_incidence ddt(u₁ *ᵢ t)
+    @test incidence == incidence"a + ∝u₁ + f(∝ₛt, ∝ₜu₂)"
+
+    incidence = @infer_incidence ddt(u₁ *ᵢ u₁)
+    @test incidence == incidence"a + f(u₁, u₂)"
+    # Note that the constant term may be removed if we
+    # model nonlinear time-independent incidences.
+
+    incidence = @infer_incidence ddt(1/u₁)
+    @test incidence == incidence"a + f(u₁, u₂)"
+
+    incidence = @infer_incidence ddt((2.0 +ᵢ u₁) *ᵢ (3.0 +ᵢ u₂))
+    @test incidence == incidence"f(∝ₛu₁, ∝ₛu₂, ∝ₛu₃, ∝ₛu₄)"
+
+    incidence = @infer_incidence ddt((2.0 +ᵢ u₁) *ᵢ (3.0 +ᵢ t))
+    @test incidence == incidence"a + ∝u₁ + f(∝ₛt, ∝ₜu₂)"
+  end
+end
 
 end