Add fallback incidence

Fixes #26, although somewhat unsatisfying, because we should be able to
get this precise by improving irinterp.

Author: Keno

Project.toml

@@ -43,14 +43,18 @@ Tracy = "e689c965-62c8-4b79-b2c5-8359227902fd"
 [weakdeps]
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 
+[sources]
+DifferentiationInterface = {rev = "main", subdir = "DifferentiationInterface", url = "https://github.com/Keno/DifferentiationInterface.jl#main"}
+SimpleNonlinearSolve = {rev = "master", subdir = "lib/SimpleNonlinearSolve", url = "https://github.com/SciML/NonlinearSolve.jl.git"}
+StateSelection = {rev = "main", url = "https://github.com/JuliaComputing/StateSelection.jl.git"}
+
 [compat]
 Accessors = "0.1.36"
 AutoHashEquals = "2.2.0"
 CentralizedCaches = "1.1.0"
 ChainRules = "1.50"
 ChainRulesCore = "1.20"
 Compiler = ">= 0.0.1"
-Cthulhu = "2.10.1"
 DiffEqBase = "6.149.2"
 DifferentiationInterface = "0.6.52"
 Diffractor = "0.2.7"

src/analysis/structural.jl

@@ -350,15 +350,14 @@ function process_ipo_return!(𝕃, ultimate_rt::Incidence, eqclassification, var
     new_eq_row = _zero_row()
     for (v_offset, coeff) in zip(rowvals(ultimate_rt.row), nonzeros(ultimate_rt.row))
         v = v_offset - 1
-        if varclassification[v] != External && coeff == nonlinear
+        if v != 0 && varclassification[v] != External && coeff == nonlinear
             get_nonlinrepl()
             new_eq_row[v_offset] = nonlinear
         else
             new_row[v_offset] = coeff
-            while true
+            while v != 0 && v !== nothing
                 varclassification[v] = External
                 v = invview(var_to_diff)[v]
-                v === nothing && break
             end
         end
     end
@@ -426,6 +425,17 @@ function add_internal_equations_to_structure!(refiner::StructuralRefiner, eqkind
     return true
 end
 
+function process_ipo_return!(𝕃, ultimate_rt::Type, eqclassification, varclassification, var_to_diff, total_incidence, eq_callee_mapping)
+    # If we don't have any internal variables (in which case we might have to to do a more aggressive rewrite), strengthen the incidence
+    # by demoting to full incidence over the argument variables. Incidence is not allowed to propagate through global mutable state, so
+    # the incidence of the return type is bounded by the incidence of the arguments in this case.
+    if !all(==(External), varclassification)
+        return ultimate_rt, 0
+    end
+    # TODO: Keep track of whether we have any time dependence?
+    return Incidence(ultimate_rt, IncidenceVector(MAX_EQS, Int[1:length(varclassification)+1;], Union{Float64, NonLinear}[nonlinear for _ in 1:length(varclassification)+1])), 0
+end
+
 function process_ipo_return!(𝕃, ultimate_rt::Eq, eqclassification, args...)
     eqclassification[ultimate_rt.id] = External
     return ultimate_rt, 0

test/regression.jl

@@ -72,4 +72,11 @@ end
 sol = solve(DAECProblem(tfb5, (1,) .=> 1.), IDA())
 @test all(map((x,y)->isapprox(x[], y, atol=1e-2), sol.u[:, 1], (sol.t .+ 1) .* log.(sol.t .+ 1) .+ 1 .- sol.t))
 
+function tfb6()
+    x = continuous()
+    always!(ddt(x) - x^0.5)
+end
+sol = solve(DAECProblem(tfb6, (1,) .=> 1.), IDA())
+@test all(map((x,y)->isapprox(x[], y, atol=1e-2), sol.u[:, 1], 1/4 .* (sol.t .+ 2).^2))
+
 end