diff --git a/Manifest.toml b/Manifest.toml
index 4f0a7cc..2bb45ff 100644
--- a/Manifest.toml
+++ b/Manifest.toml
@@ -2,7 +2,7 @@
 
 julia_version = "1.13.0-DEV"
 manifest_format = "2.0"
-project_hash = "db4caed2be78d149ffc0150f3ab04e62547353c7"
+project_hash = "d2c28a8e33664424dc750db4dae46c782768f682"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "e2478490447631aedba0823d4d7a80b2cc8cdb32"
@@ -361,11 +361,11 @@ version = "4.1.1"
 
 [[deps.Cthulhu]]
 deps = ["CodeTracking", "FoldingTrees", "InteractiveUtils", "JuliaSyntax", "PrecompileTools", "Preferences", "REPL", "TypedSyntax", "UUIDs", "Unicode", "WidthLimitedIO"]
-git-tree-sha1 = "e4dce1eefee8a3375585bce3d5c936eb797c219f"
+git-tree-sha1 = "c1e4aefb264d17e30613fad7d32779cddd019af8"
 repo-rev = "master"
 repo-url = "https://github.com/JuliaDebug/Cthulhu.jl.git"
 uuid = "f68482b8-f384-11e8-15f7-abe071a5a75f"
-version = "2.17.2"
+version = "2.17.4"
 weakdeps = ["Compiler"]
 
     [deps.Cthulhu.extensions]
@@ -420,9 +420,9 @@ version = "1.11.0"
 
 [[deps.DiffEqBase]]
 deps = ["ArrayInterface", "ConcreteStructs", "DataStructures", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "FastClosures", "FastPower", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "Parameters", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "Setfield", "Static", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface", "TruncatedStacktraces"]
-git-tree-sha1 = "6a0f6f17adbeb77688220cfb3e665c8f6e2a60ea"
+git-tree-sha1 = "1bcd3a5c585c477e5d0595937ea7b5adcda6c621"
 uuid = "2b5f629d-d688-5b77-993f-72d75c75574e"
-version = "6.172.0"
+version = "6.174.0"
 
     [deps.DiffEqBase.extensions]
     DiffEqBaseCUDAExt = "CUDA"
@@ -458,9 +458,9 @@ version = "6.172.0"
 
 [[deps.DiffEqCallbacks]]
 deps = ["ConcreteStructs", "DataStructures", "DiffEqBase", "DifferentiationInterface", "Functors", "LinearAlgebra", "Markdown", "RecipesBase", "RecursiveArrayTools", "SciMLBase", "StaticArraysCore"]
-git-tree-sha1 = "f98c17df6b2f3ac7d4a7c9b33c161b85c9b496f0"
+git-tree-sha1 = "76292e889472e810d40a844b714743c0ffb1c53b"
 uuid = "459566f4-90b8-5000-8ac3-15dfb0a30def"
-version = "4.4.1"
+version = "4.6.0"
 
 [[deps.DiffEqNoiseProcess]]
 deps = ["DiffEqBase", "Distributions", "GPUArraysCore", "LinearAlgebra", "Markdown", "Optim", "PoissonRandom", "QuadGK", "Random", "Random123", "RandomNumbers", "RecipesBase", "RecursiveArrayTools", "ResettableStacks", "SciMLBase", "StaticArraysCore", "Statistics"]
@@ -489,7 +489,7 @@ deps = ["ADTypes", "LinearAlgebra"]
 git-tree-sha1 = "ab33044805c0c6a3d85c9e75f35c3f7d765235f2"
 repo-rev = "main"
 repo-subdir = "DifferentiationInterface"
-repo-url = "https://github.com/Keno/DifferentiationInterface.jl#main"
+repo-url = "https://github.com/Keno/DifferentiationInterface.jl"
 uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 version = "0.6.52"
 
@@ -539,7 +539,7 @@ version = "0.6.52"
 
 [[deps.Diffractor]]
 deps = ["AbstractDifferentiation", "ChainRules", "ChainRulesCore", "Combinatorics", "Compiler", "Cthulhu", "InteractiveUtils", "OffsetArrays", "PrecompileTools", "StaticArrays", "StructArrays"]
-git-tree-sha1 = "5d0d60c0fc35d36bd6f60eaec296e9380aeeed18"
+git-tree-sha1 = "4b30257efc445d38523a7e38d74f54206fe42c80"
 repo-rev = "main"
 repo-url = "https://github.com/JuliaDiff/Diffractor.jl.git"
 uuid = "9f5e2b26-1114-432f-b630-d3fe2085c51c"
@@ -563,9 +563,9 @@ version = "1.11.0"
 
 [[deps.Distributions]]
 deps = ["AliasTables", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SpecialFunctions", "Statistics", "StatsAPI", "StatsBase", "StatsFuns"]
-git-tree-sha1 = "6d8b535fd38293bc54b88455465a1386f8ac1c3c"
+git-tree-sha1 = "3e6d038b77f22791b8e3472b7c633acea1ecac06"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
-version = "0.25.119"
+version = "0.25.120"
 
     [deps.Distributions.extensions]
     DistributionsChainRulesCoreExt = "ChainRulesCore"
@@ -629,9 +629,9 @@ uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 version = "1.0.5"
 
 [[deps.EnzymeCore]]
-git-tree-sha1 = "0cdb7af5c39e92d78a0ee8d0a447d32f7593137e"
+git-tree-sha1 = "1eb59f40a772d0fbd4cb75e00b3fa7f5f79c975a"
 uuid = "f151be2c-9106-41f4-ab19-57ee4f262869"
-version = "0.8.8"
+version = "0.8.9"
 weakdeps = ["Adapt"]
 
     [deps.EnzymeCore.extensions]
@@ -1083,9 +1083,9 @@ weakdeps = ["ChainRulesCore", "SparseArrays", "Statistics"]
 
 [[deps.LinearSolve]]
 deps = ["ArrayInterface", "ChainRulesCore", "ConcreteStructs", "DocStringExtensions", "EnumX", "GPUArraysCore", "InteractiveUtils", "Krylov", "LazyArrays", "Libdl", "LinearAlgebra", "MKL_jll", "Markdown", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "Setfield", "StaticArraysCore", "UnPack"]
-git-tree-sha1 = "c2685cb9d01923f0e63155149c390504e72a8fcc"
+git-tree-sha1 = "44ead334449fcbb8469309582e56ea4703e9185f"
 uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-version = "3.14.0"
+version = "3.15.0"
 
     [deps.LinearSolve.extensions]
     LinearSolveBandedMatricesExt = "BandedMatrices"
@@ -1171,9 +1171,9 @@ version = "1.11.0"
 
 [[deps.MatrixEquations]]
 deps = ["LinearAlgebra", "LinearMaps"]
-git-tree-sha1 = "928b95ce3845485bfc4571dc80df6b4d6ddcbb65"
+git-tree-sha1 = "7363bf622892e72b8281323266ef94f6c238144b"
 uuid = "99c1a7ee-ab34-5fd5-8076-27c950a045f4"
-version = "2.4.5"
+version = "2.5.0"
 
 [[deps.MatrixPencils]]
 deps = ["LinearAlgebra", "Polynomials", "Random"]
@@ -1231,7 +1231,7 @@ version = "0.3.5"
 
 [[deps.MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
-version = "2025.2.25"
+version = "2025.5.20"
 
 [[deps.MuladdMacro]]
 git-tree-sha1 = "cac9cc5499c25554cba55cd3c30543cff5ca4fab"
@@ -1292,9 +1292,9 @@ version = "1.3.0"
 
 [[deps.NonlinearSolve]]
 deps = ["ADTypes", "ArrayInterface", "BracketingNonlinearSolve", "CommonSolve", "ConcreteStructs", "DiffEqBase", "DifferentiationInterface", "FastClosures", "FiniteDiff", "ForwardDiff", "LineSearch", "LinearAlgebra", "LinearSolve", "NonlinearSolveBase", "NonlinearSolveFirstOrder", "NonlinearSolveQuasiNewton", "NonlinearSolveSpectralMethods", "PrecompileTools", "Preferences", "Reexport", "SciMLBase", "SimpleNonlinearSolve", "SparseArrays", "SparseMatrixColorings", "StaticArraysCore", "SymbolicIndexingInterface"]
-git-tree-sha1 = "7fd96e0e6585063a7193007349799155ba5a069f"
+git-tree-sha1 = "aeb6fb02e63b4d4f90337ed90ce54ceb4c0efe77"
 uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-version = "4.8.0"
+version = "4.9.0"
 
     [deps.NonlinearSolve.extensions]
     NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt"
@@ -1431,15 +1431,15 @@ version = "0.4.6"
     Reactant = "3c362404-f566-11ee-1572-e11a4b42c853"
 
 [[deps.OrderedCollections]]
-git-tree-sha1 = "cc4054e898b852042d7b503313f7ad03de99c3dd"
+git-tree-sha1 = "05868e21324cede2207c6f0f466b4bfef6d5e7ee"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.8.0"
+version = "1.8.1"
 
 [[deps.OrdinaryDiffEq]]
 deps = ["ADTypes", "Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "EnumX", "ExponentialUtilities", "FastBroadcast", "FastClosures", "FillArrays", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "InteractiveUtils", "LineSearches", "LinearAlgebra", "LinearSolve", "Logging", "MacroTools", "MuladdMacro", "NonlinearSolve", "OrdinaryDiffEqAdamsBashforthMoulton", "OrdinaryDiffEqBDF", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqExplicitRK", "OrdinaryDiffEqExponentialRK", "OrdinaryDiffEqExtrapolation", "OrdinaryDiffEqFIRK", "OrdinaryDiffEqFeagin", "OrdinaryDiffEqFunctionMap", "OrdinaryDiffEqHighOrderRK", "OrdinaryDiffEqIMEXMultistep", "OrdinaryDiffEqLinear", "OrdinaryDiffEqLowOrderRK", "OrdinaryDiffEqLowStorageRK", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqNordsieck", "OrdinaryDiffEqPDIRK", "OrdinaryDiffEqPRK", "OrdinaryDiffEqQPRK", "OrdinaryDiffEqRKN", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqSSPRK", "OrdinaryDiffEqStabilizedIRK", "OrdinaryDiffEqStabilizedRK", "OrdinaryDiffEqSymplecticRK", "OrdinaryDiffEqTsit5", "OrdinaryDiffEqVerner", "Polyester", "PreallocationTools", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "SimpleNonlinearSolve", "SimpleUnPack", "SparseArrays", "Static", "StaticArrayInterface", "StaticArrays", "TruncatedStacktraces"]
-git-tree-sha1 = "dfae5ed215f5949f52de22b92e2e42ea3e3e652d"
+git-tree-sha1 = "56d5500e9970f0112a4e1ab6474d6fedde61ef64"
 uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-version = "6.96.0"
+version = "6.97.0"
 
 [[deps.OrdinaryDiffEqAdamsBashforthMoulton]]
 deps = ["DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqLowOrderRK", "Polyester", "RecursiveArrayTools", "Reexport", "Static"]
@@ -1965,12 +1965,12 @@ version = "1.11.0"
 
 [[deps.SimpleNonlinearSolve]]
 deps = ["ADTypes", "ArrayInterface", "BracketingNonlinearSolve", "CommonSolve", "ConcreteStructs", "DifferentiationInterface", "FastClosures", "FiniteDiff", "ForwardDiff", "LineSearch", "LinearAlgebra", "MaybeInplace", "NonlinearSolveBase", "PrecompileTools", "Reexport", "SciMLBase", "Setfield", "StaticArraysCore"]
-git-tree-sha1 = "068c16a16834c1483c299b0e27e901599439570d"
+git-tree-sha1 = "7aaa5fe4617271b64fce0466d187f2a72edbd81a"
 repo-rev = "master"
 repo-subdir = "lib/SimpleNonlinearSolve"
 repo-url = "https://github.com/SciML/NonlinearSolve.jl.git"
 uuid = "727e6d20-b764-4bd8-a329-72de5adea6c7"
-version = "2.4.0"
+version = "2.5.0"
 weakdeps = ["ChainRulesCore", "DiffEqBase", "ReverseDiff", "Tracker"]
 
     [deps.SimpleNonlinearSolve.extensions]
@@ -2012,10 +2012,10 @@ uuid = "dc90abb0-5640-4711-901d-7e5b23a2fada"
 version = "0.1.2"
 
 [[deps.SparseMatrixColorings]]
-deps = ["ADTypes", "DocStringExtensions", "LinearAlgebra", "Random", "SparseArrays"]
-git-tree-sha1 = "76e9564f0de0d1d7a46095e758ae13ceba680cfb"
+deps = ["ADTypes", "DocStringExtensions", "LinearAlgebra", "PrecompileTools", "Random", "SparseArrays"]
+git-tree-sha1 = "ab958b4fec46d1f1d057bb8e2a99bfdb90744646"
 uuid = "0a514795-09f3-496d-8182-132a7b665d35"
-version = "0.4.19"
+version = "0.4.20"
 
     [deps.SparseMatrixColorings.extensions]
     SparseMatrixColoringsCliqueTreesExt = "CliqueTrees"
@@ -2194,9 +2194,9 @@ version = "3.27.0"
 
 [[deps.Symbolics]]
 deps = ["ADTypes", "ArrayInterface", "Bijections", "CommonWorldInvalidations", "ConstructionBase", "DataStructures", "DiffRules", "Distributions", "DocStringExtensions", "DomainSets", "DynamicPolynomials", "LaTeXStrings", "Latexify", "Libdl", "LinearAlgebra", "LogExpFunctions", "MacroTools", "Markdown", "NaNMath", "OffsetArrays", "PrecompileTools", "Primes", "RecipesBase", "Reexport", "RuntimeGeneratedFunctions", "SciMLBase", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArraysCore", "SymbolicIndexingInterface", "SymbolicLimits", "SymbolicUtils", "TermInterface"]
-git-tree-sha1 = "c8bc71edb71a1dcbb241a75d5f57f40804493678"
+git-tree-sha1 = "e14834f421edaa8a30493f7864dfc8582855bb3c"
 uuid = "0c5d862f-8b57-4792-8d23-62f2024744c7"
-version = "6.39.1"
+version = "6.40.0"
 
     [deps.Symbolics.extensions]
     SymbolicsForwardDiffExt = "ForwardDiff"
@@ -2258,15 +2258,15 @@ version = "1.0.0"
 
 [[deps.ThreadingUtilities]]
 deps = ["ManualMemory"]
-git-tree-sha1 = "18ad3613e129312fe67789a71720c3747e598a61"
+git-tree-sha1 = "2d529b6b22791f3e22e7ec5c60b9016e78f5f6bf"
 uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5"
-version = "0.5.3"
+version = "0.5.4"
 
 [[deps.TimerOutputs]]
 deps = ["ExprTools", "Printf"]
-git-tree-sha1 = "f57facfd1be61c42321765d3551b3df50f7e09f6"
+git-tree-sha1 = "3748bd928e68c7c346b52125cf41fff0de6937d0"
 uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
-version = "0.5.28"
+version = "0.5.29"
 
     [deps.TimerOutputs.extensions]
     FlameGraphsExt = "FlameGraphs"
diff --git a/Project.toml b/Project.toml
index 7d42d18..c1d8752 100644
--- a/Project.toml
+++ b/Project.toml
@@ -47,7 +47,7 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 [sources]
 Compiler = {rev = "master", url = "https://github.com/JuliaLang/BaseCompiler.jl.git"}
 Cthulhu = {rev = "master", url = "https://github.com/JuliaDebug/Cthulhu.jl.git"}
-DifferentiationInterface = {rev = "main", subdir = "DifferentiationInterface", url = "https://github.com/Keno/DifferentiationInterface.jl#main"}
+DifferentiationInterface = {rev = "main", subdir = "DifferentiationInterface", url = "https://github.com/Keno/DifferentiationInterface.jl"}
 Diffractor = {rev = "main", url = "https://github.com/JuliaDiff/Diffractor.jl.git"}
 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"}
diff --git a/benchmark/Project.toml b/benchmark/Project.toml
index aeb9d1b..3c43db1 100644
--- a/benchmark/Project.toml
+++ b/benchmark/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 Compiler = "807dbc54-b67e-4c79-8afb-eafe4df6f2e1"
 DAECompiler = "32805668-c3d0-42c2-aafd-0d0a9857a104"
+Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 XSteam = "95ff35a0-be81-11e9-2ca3-5b4e338e8476"
diff --git a/src/analysis/cache.jl b/src/analysis/cache.jl
index 8f007bc..5395609 100644
--- a/src/analysis/cache.jl
+++ b/src/analysis/cache.jl
@@ -60,7 +60,7 @@ function add_equation_row!(graph, solvable_graph, ieq::Int, inc::Incidence)
     for (v, coeff) in zip(rowvals(inc.row), nonzeros(inc.row))
         v == 1 && continue
         add_edge!(graph, ieq, v-1)
-        coeff !== nonlinear && add_edge!(solvable_graph, ieq, v-1)
+        isa(coeff, Float64) && add_edge!(solvable_graph, ieq, v-1)
     end
 end
 add_equation_row!(graph, solvable_graph, ieq::Int, c::Const) = nothing
diff --git a/src/analysis/ipoincidence.jl b/src/analysis/ipoincidence.jl
index 41f7ee2..bd9be52 100644
--- a/src/analysis/ipoincidence.jl
+++ b/src/analysis/ipoincidence.jl
@@ -45,49 +45,99 @@ end
 apply_linear_incidence(𝕃, ret::Type, caller::CallerMappingState, mapping::CalleeMapping) = ret
 apply_linear_incidence(𝕃, ret::Const, caller::CallerMappingState, mapping::CalleeMapping) = ret
 function apply_linear_incidence(𝕃, ret::Incidence, caller::Union{CallerMappingState, Nothing}, mapping::CalleeMapping)
-    coeffs = mapping.var_coeffs
-
-    const_val = ret.typ
-    new_row = _zero_row()
-
-    for (v_offset, coeff) in zip(rowvals(ret.row), nonzeros(ret.row))
-        v = v_offset - 1
+    # Substitute variables returned by the callee with the incidence defined by the caller.
+    # The composition will be additive in the constant terms, and multiplicative for linear coefficients.
+    caller_variables = mapping.var_coeffs
+
+    typ = ret.typ
+    row = _zero_row()
+
+    used_caller_variables = Int[]
+    for i in rowvals(ret.row)
+        i == 1 && continue # skip time
+        v = i - 1
+        if !isassigned(caller_variables, v)
+            compute_missing_coeff!(caller_variables, caller::CallerMappingState, v)
+        end
+        substitution = caller_variables[i - 1]
+        isa(substitution, Incidence) || continue
+        for j in rowvals(substitution.row)
+            push!(used_caller_variables, j)
+        end
+    end
+    did_use_time = in(1, used_caller_variables)
 
+    for (i, coeff) in zip(rowvals(ret.row), nonzeros(ret.row))
         # Time dependence persists as itself
-        if v == 0
-            new_row[v_offset] += coeff
+        if i == 1
+            row[i] = coeff
             continue
         end
 
-        if !isassigned(coeffs, v)
-            @assert caller !== nothing
-            compute_missing_coeff!(coeffs, caller, v)
-        end
-
-        replacement = coeffs[v]
-        if isa(replacement, Incidence)
-            new_row .+= replacement.row .* coeff
-        else
-            if isa(replacement, Const)
-                if isa(const_val, Const)
-                    new_const_val = const_val.val + replacement.val * coeff
-                    if isa(new_const_val, Float64)
-                        const_val = Const(new_const_val)
-                    else
-                        const_val = widenconst(const_val)
-                    end
+        v = i - 1
+        substitution = caller_variables[v]
+        if isa(substitution, Incidence)
+            # Distribute the coefficient to all terms.
+            # Because the coefficient is expressed in the reference of the callee,
+            # state dependence must be processed carefully.
+            typ = compose_additive_term(typ, substitution.typ, coeff)
+            for (j, substitute) in zip(rowvals(substitution.row), nonzeros(substitution.row))
+                row[j] === nonlinear && continue # no more information to be gained
+                if substitute === nonlinear || coeff === nonlinear
+                    row[j] = nonlinear
+                elseif isa(coeff, Float64)
+                    row[j] += coeff * substitute
                 else
-                    const_val = widenconst(const_val)
+                    time_dependent = coeff.time_dependent
+                    state_dependent = false
+                    if isa(substitute, Linearity)
+                        time_dependent |= substitute.time_dependent
+                        state_dependent |= substitute.state_dependent
+                    end
+                    if coeff.state_dependent
+                        if coeff.time_dependent && did_use_time
+                            # The term is at least bilinear in another state, and this state
+                            # from the callee may alias time from the caller, so we must mark
+                            # time as nonlinear.
+                            row[1] = nonlinear
+                        end
+                        if count(==(j), used_caller_variables) > 1
+                            # The term is at least bilinear in another state, but we don't
+                            # know which state, so we must fall back to nonlinear.
+                            row[j] = nonlinear
+                            continue
+                        end
+                        # We'll only be state-dependent if variables from the callee
+                        # map to at least one other variable than `j`.
+                        if j == 1 # time
+                            state_dependent |= length(used_caller_variables) > 1
+                        else # state
+                            state_dependent |= length(used_caller_variables) - did_use_time > 1
+                        end
+                        # If another state may contain time, we may be time-dependent too.
+                        time_dependent |= did_use_time
+                    end
+                    j == 1 && (time_dependent = false)
+                    row[j] += Linearity(; nonlinear = false, state_dependent, time_dependent)
                 end
-            else
-                # The replacement has some unknown type - we need to widen
-                # all the way here.
-                return widenconst(const_val)
             end
+        elseif isa(substitution, Const)
+            typ = compose_additive_term(typ, substitution, coeff)
+        else
+            return widenconst(typ) # unknown lattice element, we should widen
         end
     end
 
-    return Incidence(const_val, new_row)
+    return Incidence(typ, row)
+end
+
+function compose_additive_term(@nospecialize(a), @nospecialize(b), coeff)
+    isa(a, Const) || return widenconst(a)
+    isa(b, Const) || return widenconst(a)
+    isa(coeff, Linearity) && return b.val == 0 ? a : widenconst(a)
+    val = a.val + b.val * coeff
+    isa(val, Float64) || return widenconst(a)
+    return Const(val)
 end
 
 function apply_linear_incidence(𝕃, ret::Eq, caller::CallerMappingState, mapping::CalleeMapping)
diff --git a/src/analysis/lattice.jl b/src/analysis/lattice.jl
index f42b8b0..ffbc737 100644
--- a/src/analysis/lattice.jl
+++ b/src/analysis/lattice.jl
@@ -39,56 +39,89 @@ Compiler.widenlattice(::EqStructureLattice) = Compiler.ConstsLattice()
 Compiler.is_valid_lattice_norec(::EqStructureLattice, @nospecialize(v)) = isa(v, Incidence) || isa(v, Eq) || isa(v, PartialScope) || isa(v, PartialKeyValue)
 Compiler.has_extended_unionsplit(::EqStructureLattice) = true
 
-############################## NonLinear #######################################
+############################## Linearity #######################################
 
 """
-    struct NonLinear
-
-This singleton number is similar to `missing` in that arithmatic with it is
-saturating. When used as a coefficient in the Incidence lattice, it indicates
-that the corresponding variable does not have a (constant) linear coefficient.
-This may either mean that the variable in question has a non-constant linear
-coefficient or that the variable is used non-linearly. We do not currently
-distinguish the two situations.
+    Base.@kwdef struct Linearity
+        time_dependent::Bool = true
+        state_dependent::Bool = true
+        nonlinear::Bool = true
+    end
+
+Together with `Float64` values in `Incidence`, this struct expresses linearity information:
+- Linear, coefficient is a known constant (a `Float64` value).
+- Linear, coefficient is an unknown constant (`linear`).
+- Linear, coefficient is unknown with a known state or time dependence (`linear_state_dependent`, `linear_time_dependent`).
+- Linear, coefficient is unknown with both state and time dependence (`linear_time_and_state_dependent`).
+- Nonlinear (`nonlinear` - always taints state and time dependence currently).
+
+A known constant-coefficient contains the highest level of information. `nonlinear` contains the lowest
+(and most conservative) level of information, and should be assumed used by default or if uncertain.
 """
-struct NonLinear; end
-Base.iszero(::NonLinear) = false
-Base.zero(::Type{Union{NonLinear, Float64}}) = 0.
+Base.@kwdef struct Linearity
+    time_dependent::Bool = true
+    state_dependent::Bool = true
+    nonlinear::Bool = true
+    function Linearity(time_dependent, state_dependent, nonlinear)
+        if nonlinear && (!time_dependent || !state_dependent)
+            throw(ArgumentError("Modeling of state or time independence is not supported for nonlinearities"))
+        end
+        new(time_dependent, state_dependent, nonlinear)
+    end
+end
+
+"The variable is used linearly, with an unknown constant."
+const linear = Linearity(time_dependent = false, state_dependent = false, nonlinear = false)
+"The variable is used linearly, with an unknown constant that may depend on time."
+const linear_time_dependent = Linearity(state_dependent = false, nonlinear = false)
+"The variable is used linearly, with an unknown constant that may depend on time and on other states."
+const linear_state_dependent = Linearity(time_dependent = false, nonlinear = false)
+"The variable is used linearly, with an unknown constant that may depend on time and on other states."
+const linear_time_and_state_dependent = Linearity(nonlinear = false)
+"The variable is used nonlinearly, with a possible dependence on time and other states."
+const nonlinear = Linearity()
+
+join_linearity(a::Linearity, b::Real) = a
+join_linearity(a::Real, b::Linearity) = b
+join_linearity(a::Real, b::Real) = a == b ? a : linear
+function join_linearity(a::Linearity, b::Linearity)
+    (a.nonlinear | b.nonlinear) && return nonlinear
+    return Linearity(; time_dependent = a.time_dependent | b.time_dependent, state_dependent = a.state_dependent | b.state_dependent, nonlinear = false)
+end
+
+Base.iszero(::Linearity) = false
+Base.zero(::Type{Union{Linearity, Float64}}) = 0.
 for f in (:+, :-)
     @eval begin
-        Base.$f(a::Real, b::NonLinear) = b
-        Base.$f(a::NonLinear, b::Real) = a
-        Base.$f(a::NonLinear, b::NonLinear) = nonlinear
-        Base.$f(::NonLinear) = nonlinear
-    end
-end
-
-Base.:(*)(a::Real, b::NonLinear) = iszero(a) ? a : b
-Base.:(*)(a::NonLinear, b::Real) = iszero(b) ? b : a
-Base.:(*)(a::NonLinear, b::NonLinear) = nonlinear
-Base.div(a::Real, b::NonLinear) = iszero(a) ? a : b
-Base.div(a::NonLinear, b::Real) = a
-Base.div(a::NonLinear, b::NonLinear) = a
-Base.:(/)(a::Real, b::NonLinear) = iszero(a) ? a : b
-Base.:(/)(a::NonLinear, b::Real) = a
-Base.:(/)(a::NonLinear, b::NonLinear) = a
-Base.rem(a::Real, b::NonLinear) = b
-Base.rem(a::NonLinear, b::Real) = a
-Base.rem(a::NonLinear, b::NonLinear) = a
-Base.abs(a::NonLinear) = nonlinear
-Base.isone(a::NonLinear) = false
-
-const nonlinear = NonLinear.instance
-Base.Broadcast.broadcastable(::NonLinear) = Ref(nonlinear)
+        Base.$f(a::Real, b::Linearity) = b
+        Base.$f(a::Linearity, b::Real) = a
+        Base.$f(a::Linearity, b::Linearity) = join_linearity(a, b)
+        Base.$f(a::Linearity) = a
+    end
+end
+
+Base.:(*)(a::Real, b::Linearity) = iszero(a) ? a : b
+Base.:(*)(a::Linearity, b::Real) = iszero(b) ? b : a
+Base.div(a::Real, b::Linearity) = iszero(a) ? a : nonlinear
+Base.div(a::Linearity, b::Real) = a
+Base.div(a::Linearity, b::Linearity) = nonlinear
+Base.:(/)(a::Real, b::Linearity) = iszero(a) ? a : nonlinear
+Base.:(/)(a::Linearity, b::Real) = a
+Base.:(/)(a::Linearity, b::Linearity) = nonlinear
+Base.rem(a::Real, b::Linearity) = b
+Base.rem(a::Linearity, b::Real) = a
+Base.rem(a::Linearity, b::Linearity) = a
+Base.abs(a::Linearity) = nonlinear
+Base.isone(a::Linearity) = false
+
+Base.Broadcast.broadcastable(x::Linearity) = Ref(x)
 
 ############################## Incidence #######################################
 # TODO: Just use Infinities.jl here?
 const MAX_EQS = div(typemax(Int), 2)
 
-# For now, we only track exact, integer linearities, because that's what
-# MTK can handle, so `nonlinear` includes linear operations with floating
-# point values.
-const IncidenceVector = SparseVector{Union{Float64, NonLinear}, Int}
+const IncidenceValue = Union{Float64, Linearity}
+const IncidenceVector = SparseVector{IncidenceValue, Int}
 
 is_non_incidence_type(@nospecialize(type)) = type === Union{} || Base.issingletontype(type)
 
@@ -102,6 +135,9 @@ elements are defined by subset inclusion. Note that in particular this implies
 that plain `T` lattice elements have unknown incidence and `Const` lattice elements
 have no incidence. A lattice element of type `T` that is known to be state-independent
 would have lattice element `Incidence(T, {})`.
+
+For convenience, you may want to use the `incidence"..."` string macro to construct an
+[`Incidence`](@ref) from its printed representation.
 """
 struct Incidence
     # Considered additive to `row`. In particular, if the `typ` is Float64,
@@ -113,6 +149,31 @@ struct Incidence
         if is_non_incidence_type(type)
             throw(DomainError(type, "Invalid type for Incidence"))
         end
+        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))
+            throw(ArgumentError("Time incidence cannot be both linear and time-dependent, otherwise it would be nonlinear"))
+        end
+        for (i, coeff) in zip(rowvals(row), nonzeros(row))
+            isa(coeff, Linearity) || continue
+            coeff.nonlinear && continue
+            if coeff.time_dependent && !in(1, rowvals(row))
+                throw(ArgumentError("Time-dependent incidence annotation for $(subscript_state(i)) is inconsistent with an absence of time incidence"))
+            end
+            if coeff.state_dependent && !any(x -> x != 1, rowvals(row))
+                throw(ArgumentError("State-dependent incidence annotation for $(subscript_state(i)) is inconsistent with an absence of state incidence"))
+            end
+            if i > 1 && coeff.time_dependent && (isa(time, Float64) || !time.state_dependent)
+                throw(ArgumentError("Time-dependent state incidence for $(subscript_state(i)) is inconsistent with an absence of state dependence for time"))
+            end
+        end
         return new(type, row)
     end
 end
@@ -134,47 +195,185 @@ end
 function Base.show(io::IO, inc::Incidence)
     print(io, "Incidence(")
     first = true
+    print_plusminus(io, minus=false) = if first
+        first = false
+        minus && print(io, "-")
+    else
+        print(io, minus ? " - " : " + ")
+    end
     if isa(inc.typ, Const) && isa(inc.typ.val, Float64)
         if !iszero(inc.typ.val)
             print(io, inc.typ.val)
             first = false
         end
-    else inc.typ !== Float64
-        print(io, inc.typ, ", ")
-    end
-    print_plusminus(io, minus=false) = if first
-            first = false
-            minus && print(io, "-")
+    else
+        if inc.typ === Float64
+            print(io, 'a')
+            !isempty(rowvals(inc.row)) && print_plusminus(io)
         else
-            print(io, minus ? " - " : " + ")
+            print(io, inc.typ)
+            !isempty(rowvals(inc.row)) && print(io, ", ")
         end
+    end
+    time = inc.row[1]
+    is_grouped(v, i) = isa(v, Linearity) && (v.state_dependent || (v.time_dependent || i == 1) && in(time, (linear_state_dependent, nonlinear)))
+    function propto(linearity::Linearity)
+        str = "∝"
+        linearity.time_dependent && (str *= 'ₜ')
+        linearity.state_dependent && (str *= 'ₛ')
+        return str
+    end
     for (i, v) in zip(rowvals(inc.row), nonzeros(inc.row))
-        v !== nonlinear || continue
-        print_plusminus(io, v < 0)
-        if abs(v) != 1 || i == 1
-            print(io, abs(v))
+        is_grouped(v, i) && continue
+        if isa(v, Float64)
+            print_plusminus(io, v < 0)
+            if abs(v) != 1
+                print(io, abs(v))
+            end
+        else
+            !first && print(io, " + ")
+            first = false
+            if is_grouped(inc.row[1], 1) && v.time_dependent
+                print(io, propto(inc.row[1]::Linearity), 't', " * ")
+            else # unknown constant coefficient
+                print(io, propto(v))
+            end
         end
         print(io, subscript_state(i))
     end
-    first_nonlinear = true
-    for (i, v) in zip(rowvals(inc.row), nonzeros(inc.row))
-        v === nonlinear || continue
-        if first_nonlinear
-            print_plusminus(io)
-            first_nonlinear = false
-            print(io, "f(")
+    first_grouped = true
+    if any(i -> is_grouped(inc.row[i], i), rowvals(inc.row))
+        for (i, v) in zip(rowvals(inc.row), nonzeros(inc.row))
+            !is_grouped(v, i) && continue
+            if first_grouped
+                print_plusminus(io)
+                first_grouped = false
+                print(io, "f(")
+            else
+                print(io, ", ")
+            end
+            !v.nonlinear && print(io, propto(v))
+            print(io, subscript_state(i))
+        end
+        if !first_grouped
+            print(io, ")")
+        end
+    end
+    print(io, ")")
+end
+
+"""
+    incidence"a + f(∝ₛt, u₁)"
+    incidence"Incidence(a + f(∝ₛt, u₁))" # you may copy-paste straight from its printed output
+
+Construct an [`Incidence`](@ref) from its printed representation.
+"""
+macro incidence_str(str) generate_incidence(str) end
+
+function generate_incidence(str::String)
+    if startswith(str, "Incidence(") && endswith(str, ')')
+        # Support `incidence"Incidence(...)"` so the user doesn't have to
+        # manually remove the `Incidence` call when copy-pasting.
+        str = str[11:(end - 1)]
+    end
+    str = replace(str, '∝' => '~')
+    ex = Meta.parse(str)
+    generate_incidence(ex)
+end
+
+function generate_incidence(ex)
+    T = nothing
+    if isexpr(ex, :tuple, 2)
+        T, ex = ex.args[1], ex.args[2]
+    end
+    generate_incidence(T, ex)
+end
+
+function generate_incidence(T, ex)
+    if isexpr(ex, :call) && ex.args[1] === :+
+        terms = ex.args[2:end]
+    else
+        terms = Any[ex]
+    end
+    pairs = Dict{Int,Any}()
+    for term in terms
+        if term === :a
+            T === nothing || throw(ArgumentError("The incidence type must not be provided if a constant `Float64` term is already present"))
+            T = Float64
+            continue
+        elseif isa(term, Float64)
+            T === nothing || throw(ArgumentError("The incidence type must not be provided if a literal `Float64` term is already present"))
+            T = Const(term)
+            continue
+        end
+
+        @assert isa(term, Symbol) || isexpr(term, :call)
+
+        ispropto(x) = isexpr(x, :call, 2) && startswith(string(x.args[1]), '~')
+
+        if isa(term, Symbol)
+            i = parse_variable(string(term))
+            pairs[i] = 1.0
+        elseif isexpr(term, :call, 3) && term.args[1] === :*
+            factor = parse(Float64, string(term.args[2]))
+            i = parse_variable(string(term.args[3]))
+            pairs[i] = factor
+        elseif ispropto(term)
+            coefficient, variable = separate_coefficient_and_variable(term)
+            coefficient = parse_coefficient(coefficient)
+            i = parse_variable(variable)
+            pairs[i] = coefficient
+        elseif isexpr(term, :call) && term.args[1] === :f
+            for argument in @view term.args[2:end]
+                if ispropto(argument)
+                    coefficient, variable = separate_coefficient_and_variable(argument)
+                    coefficient = parse_coefficient(coefficient)
+                    i = parse_variable(variable)
+                else
+                    i = parse_variable(string(argument))
+                    coefficient = nonlinear
+                end
+                pairs[i] = coefficient
+            end
         else
-            print(io, ", ")
+            throw(ArgumentError("Unrecognized call to function '$(term.args[1])'"))
         end
-        print(io, i == 1 ? "t" : subscript_state(i))
     end
-    if !first_nonlinear
-        print(io, ")")
+    values = IncidenceValue[]
+    for i in 1:maximum(keys(pairs); init = 0)
+        val = get(pairs, i, 0.0)
+        isa(val, Pair) && (val = val.second)
+        push!(values, val)
     end
-    print(io, ")")
+    T = something(T, Const(0.0))
+    :(Incidence($T, IncidenceValue[$(values...)]))
 end
 
-_zero_row() = IncidenceVector(MAX_EQS, Int[], Union{Float64, NonLinear}[])
+function separate_coefficient_and_variable(term::Expr)
+    str = string(term)
+    i = findfirst(in(('t', 'u')), str)::Int
+    @view(str[1:prevind(str, i)]), @view(str[i:end])
+end
+
+function parse_coefficient(coefficient::AbstractString)
+    matched = match(r"^~(ₜ)?(ₛ)?$", coefficient)
+    @assert matched !== nothing
+    time_dependent = matched.captures[1] !== nothing
+    state_dependent = matched.captures[2] !== nothing
+    return Linearity(; time_dependent, state_dependent, nonlinear = false)
+end
+
+function parse_variable(term)
+    term == "t" && return 1
+    matched = match(r"^u([₀₁₂₃₄₅₆₇₈₉]+)$", term)
+    @assert matched !== nothing
+    capture = matched[1]
+    return 1 + parse(Int, map(subscript_to_number, capture))
+end
+
+subscript_to_number(char) = Char(48 + (UInt32(char) - 8320))
+
+_zero_row() = IncidenceVector(MAX_EQS, Int[], IncidenceValue[])
 const _ZERO_ROW = _zero_row()
 const _ZERO_CONST = Const(0.0)
 Base.zero(::Type{Incidence}) = Incidence(_ZERO_CONST, _zero_row())
@@ -514,11 +713,7 @@ function Compiler.tmerge(🥬::EqStructureLattice, @nospecialize(a), @nospeciali
             merged_typ = Compiler.tmerge(Compiler.widenlattice(🥬), a.typ, b.typ)
             row = _zero_row()
             for i in union(rowvals(a.row), rowvals(b.row))
-                if a.row[i] == b.row[i]
-                    row[i] = a.row[i]
-                else
-                    row[i] = nonlinear
-                end
+                row[i] = join_linearity(a.row[i], b.row[i])
             end
             return Incidence(merged_typ, row)
         elseif isa(b, Const)
@@ -590,17 +785,14 @@ function _aggressive_incidence_join(@nospecialize(rt), argtypes::Vector{Any})
             for (i, v) in zip(rowvals(a.row), nonzeros(a.row))
                 # as long as they are equal then it is correct for both so nothing to do
                 if inci.row[i] != v
-                    # Otherwise it can't be either but must allow both. We would ideally represent this as
-                    # `LinearUnion{rr[i], v}` or `Linear`, but we don't have lattice elements like that
-                    # `NonLinear` is our more general representation
-                    inci.row[i] = nonlinear
+                    # Otherwise it can't be either but must allow both.
+                    inci.row[i] = join_linearity(inci.row[i], v)
                 end
             end
             # and the the other way: catch places that `rr` is nonzero but `aa` is zero
             for (i, v) in zip(rowvals(inci.row), nonzeros(inci.row))
                 if a.row[i] != v
-                    # mix of nonlinear and linear, or again: a mix of two different linear coefficients
-                    inci.row[i] = nonlinear
+                    inci.row[i] = join_linearity(a.row[i], v)
                 end
             end
         end
diff --git a/src/analysis/refiner.jl b/src/analysis/refiner.jl
index fab10f4..aaf6d76 100644
--- a/src/analysis/refiner.jl
+++ b/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
@@ -93,28 +93,57 @@ function structural_inc_ddt(var_to_diff::DiffGraph, varclassification::Union{Vec
         return dv + 1
     end
     base = isa(inc.typ, Const) ? Const(zero(inc.typ.val)) : inc.typ
-    for (v_offset, coeff) in zip(rowvals(inc.row), nonzeros(inc.row))
-        if v_offset == 1
-            # t
-            if isa(base, Const)
-                if isa(coeff, Float64)
+    indices = rowvals(inc.row)
+    for (i, coeff) in zip(indices, nonzeros(inc.row))
+        if i == 1 # time
+            if isa(coeff, Float64) # known constant coefficient
+                # Do not set r[i]; d/dt t = 1
+                if isa(base, Const)
                     base = Const(base.val + coeff)
-                    # Do not set r[v_offset]; d/dt t = 1
-                else
-                    r[v_offset] = nonlinear
                 end
-            elseif !isa(coeff, Const)
-                r[v_offset] = nonlinear
+            elseif coeff.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[i] = coeff
+                else # unknown constant coefficient
+                    if isa(base, Const)
+                        # We are adding an unknown but constant value to the
+                        # result, additive `Const` information is no longer accurate.
+                        base = widenconst(base)
+                    end
+                end
             end
             continue
         end
         if isa(coeff, Float64)
-            # Linear, just add to the derivative
-            r[get_or_make_diff(v_offset)] += coeff
-        else
-            # nonlinear
-            r[v_offset] = nonlinear
-            r[get_or_make_diff(v_offset)] = nonlinear
+            # Linear with a known constant coefficient, just add to the derivative
+            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(i)] = coeff
+        elseif coeff.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
+            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)
@@ -179,8 +208,30 @@ function tfunc(::Union{Val{Core.Intrinsics.mul_float}, Val{Core.Intrinsics.mul_i
         return Incidence(builtin_math_tfunc(Core.Intrinsics.mul_float, a.typ, b), a.row * b.val)
     end
     rrow = _zero_row()
-    for i in Iterators.flatten((rowvals(a.row), rowvals(b.row)))
-        rrow[i] = nonlinear
+    ia = rowvals(a.row)
+    ib = rowvals(b.row)
+    states_a = filter(≠(1), ia)
+    states_b = filter(≠(1), ib)
+    time_dependent_a = in(1, ia)
+    time_dependent_b = in(1, ib)
+    for i in Iterators.flatten((ia, ib))
+        x = a.row[i]
+        y = b.row[i]
+        if x ≠ 0 && y ≠ 0
+            rrow[i] = nonlinear # uᵢ²
+            continue
+        end
+        if x == 0
+            val = y
+            state_dependent = !isempty(states_a)
+            time_dependent = time_dependent_a
+        else; @assert y == 0
+            val = x
+            state_dependent = !isempty(states_b)
+            time_dependent = time_dependent_b
+        end
+        val = join_linearity(val, Linearity(; time_dependent, state_dependent, nonlinear = false))
+        rrow[i] = val
     end
     return Incidence(builtin_math_tfunc(Core.Intrinsics.mul_float, a.typ, b.typ), rrow)
 end
@@ -217,7 +268,9 @@ function tfunc(::Val{Core.Intrinsics.div_float}, @nospecialize(a::Union{Const, T
     if isa(a, Incidence)
         rrow .+= a.row
     end
-    rrow .*= nonlinear
+    for i in rowvals(rrow)
+        rrow[i] = nonlinear
+    end
     return Incidence(builtin_math_tfunc(Core.Intrinsics.div_float, isa(a, Incidence) ? a.typ : a, widenconst(b.typ)), rrow)
 end
 
diff --git a/src/analysis/structural.jl b/src/analysis/structural.jl
index 5c04af6..ba187c3 100644
--- a/src/analysis/structural.jl
+++ b/src/analysis/structural.jl
@@ -139,7 +139,7 @@ function _structural_analysis!(ci::CodeInstance, world::UInt)
             if isa(newT, Const) && newT.val === Intrinsics.ScopeIdentity
                 # Allocate the identity now. After inlining, we're guaranteed that
                 # every Expr(:new) uniquely corresponds to a scope identity, so this
-                # is legal here (bug not before)
+                # is legal here (but not before)
                 compact[SSAValue(i)][:stmt] = Intrinsics.ScopeIdentity()
                 compact[SSAValue(i)][:flag] |= Compiler.IR_FLAG_REFINED
             end
@@ -186,6 +186,8 @@ function _structural_analysis!(ci::CodeInstance, world::UInt)
     end
 
     # Go through each variable we previously identified and record the (post-propagation scope and kind)
+    # For this we inspect the `kind` and `scope` arguments from `variable(kind, scope)` invokes, then
+    # we turn them into `nothing` when we are done.
     for (ivar, ssa) in enumerate(varssa)
         isa(ssa, SSAValue) || continue
         inst = ir[ssa][:inst]::Union{Nothing, Expr}
@@ -349,7 +351,6 @@ function _structural_analysis!(ci::CodeInstance, world::UInt)
     end
     ir = Compiler.finish(compact)
 
-    nimplicitoutpairs = 0
     var_to_diff = StateSelection.complete(var_to_diff)
     ultimate_rt, nimplicitoutpairs = process_ipo_return!(Compiler.typeinf_lattice(refiner), ultimate_rt, eqclassification, eqkinds, varclassification, varkinds,
         var_to_diff, total_incidence, eq_callee_mapping)
@@ -386,9 +387,10 @@ function process_ipo_return!(𝕃, ultimate_rt::Incidence, eqclassification, eqk
     new_eq_row = _zero_row()
     for (v_offset, coeff) in zip(rowvals(ultimate_rt.row), nonzeros(ultimate_rt.row))
         v = v_offset - 1
-        if v != 0 && varclassification[v] != External && coeff == nonlinear
+        if v != 0 && varclassification[v] != External && coeff === nonlinear
+            # nonlinear state variable created within this call tree
             get_nonlinrepl()
-            new_eq_row[v_offset] = nonlinear
+            new_eq_row[v_offset] = nonlinear # we might want to refine this to something linear
         else
             new_row[v_offset] = coeff
             while v != 0 && v !== nothing
@@ -470,7 +472,7 @@ function process_ipo_return!(𝕃, ultimate_rt::Type, eqclassification, eqkinds,
         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
+    return Incidence(ultimate_rt, IncidenceVector(MAX_EQS, Int[1:length(varclassification)+1;], IncidenceValue[nonlinear for _ in 1:length(varclassification)+1])), 0
 end
 
 function process_ipo_return!(𝕃, ultimate_rt::Eq, eqclassification, args...)
diff --git a/src/reflection.jl b/src/reflection.jl
index 94aa4e8..7bb595f 100644
--- a/src/reflection.jl
+++ b/src/reflection.jl
@@ -1,5 +1,5 @@
 export code_structure_by_type, code_structure, @code_structure,
-       code_ad_by_type, code_ad, @code_ad
+       code_ad_by_type, code_ad, @code_ad, @incidence_str
 
 function code_structure(@nospecialize(f), @nospecialize(types = Base.default_tt(f)); kwargs...)
   tt = Base.signature_type(f, types)
@@ -30,6 +30,7 @@ function code_structure_by_type(@nospecialize(tt::Type); world::UInt = Base.tls_
   ci = _code_ad_by_type(tt; world, kwargs...)
   # Perform or lookup DAECompiler specific analysis for this system.
   _result = structural_analysis!(ci, world)
+  isa(_result, UncompilableIPOResult) && throw(_result.error)
   return result ? _result : _result.ir
 end
 
diff --git a/src/transform/autodiff/index_lowering.jl b/src/transform/autodiff/index_lowering.jl
index d7465a9..448068d 100644
--- a/src/transform/autodiff/index_lowering.jl
+++ b/src/transform/autodiff/index_lowering.jl
@@ -72,7 +72,7 @@ function index_lowering_ad!(state::TransformationState, key::TornCacheKey)
     # Mark all non-trivial `ddt()` statements as ones that we should differentiate
     diff_ssas = Pair{SSAValue,Int}[]
     for i = 1:length(ir.stmts)
-        if is_known_invoke(ir.stmts[i][:stmt], ddt, ir) && !is_const_plus_var_linear(argextype(ir.stmts[i][:stmt].args[end], ir))
+        if is_known_invoke(ir.stmts[i][:stmt], ddt, ir) && !is_const_plus_var_known_linear(argextype(ir.stmts[i][:stmt].args[end], ir))
             push!(diff_ssas, SSAValue(i) => 0)
         end
     end
diff --git a/src/transform/state_selection.jl b/src/transform/state_selection.jl
index 99b5c48..8ec1246 100644
--- a/src/transform/state_selection.jl
+++ b/src/transform/state_selection.jl
@@ -14,7 +14,7 @@ function StateSelection.linear_subsys_adjmat!(state::TransformationState)
     linear_equations = Vector{Int}()
     for (i, inc) in enumerate(state.total_incidence)
         isa(inc, Const) && continue
-        any(x->x === nonlinear || !isinteger(x), nonzeros(inc.row)) && continue
+        any(x->isa(x, Linearity) || !isinteger(x), nonzeros(inc.row)) && continue
         (isa(inc.typ, Const) && iszero(inc.typ.val)) || continue
         # Skip any equations involving `t` for now
         # TODO: We may want to adjust our ILS to include equations with `t` in them.
@@ -92,7 +92,7 @@ function ssrm!(state::TransformationState)
         for g in (s.graph, s.solvable_graph)
             StateSelection.set_neighbors!(g, e, ils.row_cols[ei])
         end
-        state.total_incidence[e] = Incidence(Const(0.), IncidenceVector(MAX_EQS, map(x->x+1, ils.row_cols[ei]), Union{Float64, NonLinear}[Float64(x) for x in ils.row_vals[ei]]))
+        state.total_incidence[e] = Incidence(Const(0.), IncidenceVector(MAX_EQS, map(x->x+1, ils.row_cols[ei]), IncidenceValue[Float64(x) for x in ils.row_vals[ei]]))
     end
 end
 
diff --git a/src/transform/tearing/schedule.jl b/src/transform/tearing/schedule.jl
index dd27a41..56ba148 100644
--- a/src/transform/tearing/schedule.jl
+++ b/src/transform/tearing/schedule.jl
@@ -80,7 +80,7 @@ function schedule_incidence!(compact, curval, incT::Incidence, var, line; vars=n
             continue
         end
 
-        coeff == nonlinear && continue
+        isa(coeff, Float64) || continue
 
         if lin_var == 0
             lin_var_ssa = insert_node_here!(compact,
@@ -120,12 +120,12 @@ Base.setindex!(rno::RenameOverlayVector, @nospecialize(val), i::Int) =
     Base.setindex!(rno.overlay, val, i)
 Base.size(rno::RenameOverlayVector) = Base.size(rno.base)
 
-is_var_part_linear(incT::Const) = true
-is_var_part_linear(incT::Incidence) = !any(==(nonlinear), incT.row)
-is_const_plus_var_linear(incT) = is_var_part_linear(incT) && isa(incT.typ, Const)
-is_const_plus_var_linear(incT::Const) = true
+is_var_part_known_linear(incT::Const) = true
+is_var_part_known_linear(incT::Incidence) = all(x -> isa(x, Float64), incT.row)
+is_const_plus_var_known_linear(incT) = is_var_part_known_linear(incT) && isa(incT.typ, Const)
+is_const_plus_var_known_linear(incT::Const) = true
 
-is_fully_state_linear(incT, param_vars) = is_const_plus_var_linear(incT) && is_fully_state_linear(incT.typ, param_vars)
+is_fully_state_linear(incT, param_vars) = is_const_plus_var_known_linear(incT) && is_fully_state_linear(incT.typ, param_vars)
 is_fully_state_linear(incT::Const, param_vars) = iszero(incT.val)
 
 function schedule_nonlinear!(compact, param_vars, var_eq_matching, ir, ordinal, val::Union{SSAValue, Argument}, ssa_rename::AbstractVector{Any}; vars, schedule_missing_var! = nothing)
@@ -181,7 +181,7 @@ function schedule_nonlinear!(compact, param_vars, var_eq_matching, ir, ordinal,
         end
 
         # TODO: SICM
-        if !is_const_plus_var_linear(typ::Incidence)
+        if !is_const_plus_var_known_linear(typ::Incidence)
             this_nonlinear = schedule_nonlinear!(compact, param_vars, var_eq_matching, ir, ordinal, arg, ssa_rename; vars, schedule_missing_var!)
         else
             if @isdefined(result)
@@ -201,7 +201,7 @@ function schedule_nonlinear!(compact, param_vars, var_eq_matching, ir, ordinal,
         return schedule_incidence!(compact, this_nonlinear, typ, -1, inst[:line]; vars, schedule_missing_var!)[1]
     end
 
-    if is_const_plus_var_linear(incT)
+    if is_const_plus_var_known_linear(incT)
         ret = schedule_incidence!(compact, nothing, info.result.extended_rt, -1, inst[:line]; vars=
             [arg === nothing ? 0.0 : arg for arg in args], schedule_missing_var! = var->error((var, incT, args)))[1]
     else
@@ -328,7 +328,7 @@ function compute_eq_schedule(key::TornCacheKey, total_incidence, result, mss::St
                     i in this_callee_eqs && continue # We already scheduled this
                     callee_incidence = callee_info.result.total_incidence[i]
                     incidence = apply_linear_incidence(nothing, callee_incidence, nothing, callee_info.mapping)
-                    if is_const_plus_var_linear(incidence)
+                    if is_const_plus_var_known_linear(incidence)
                         # No non-linear components - skip it
                         push!(previously_scheduled_or_ignored, i)
                         continue
@@ -363,7 +363,7 @@ function compute_eq_schedule(key::TornCacheKey, total_incidence, result, mss::St
         end
 
         function schedule_eq!(eq; force=false)
-            if is_const_plus_var_linear(total_incidence[eq])
+            if is_const_plus_var_known_linear(total_incidence[eq])
                 # This is a linear equation, we can schedule it now
                 return true
             end
@@ -381,7 +381,7 @@ function compute_eq_schedule(key::TornCacheKey, total_incidence, result, mss::St
                 schedule = get!(callee_schedules, ssa, Vector{Pair{BitSet, BitSet}}())
                 callee_info = result.ir[SSAValue(ssa.id)][:info]::MappingInfo
 
-                if is_const_plus_var_linear(callee_info.result.total_incidence[callee_eq])
+                if is_const_plus_var_known_linear(callee_info.result.total_incidence[callee_eq])
                     # This portion of the calle is linear, we can schedule it
                     continue
                 end
@@ -625,7 +625,7 @@ function matching_for_key(state::TransformationState, key::TornCacheKey)
             invview(var_eq_matching)[eq] = WrongEquation()
             continue
         end
-        if is_const_plus_var_linear(total_incidence[eq]) && invview(var_eq_matching)[eq] === unassigned
+        if is_const_plus_var_known_linear(total_incidence[eq]) && invview(var_eq_matching)[eq] === unassigned
             invview(var_eq_matching)[eq] = FullyLinear()
         end
     end
@@ -786,8 +786,8 @@ function tearing_schedule!(state::TransformationState, ci::CodeInstance, key::To
                         push!(callee_explicit_eqs, callee_eq)
                     else
                         var = invview(var_eq_matching)[caller_eq]
-                        (is_const_plus_var_linear(callee_result.total_incidence[callee_eq]) ||
-                         is_const_plus_var_linear(total_incidence[caller_eq])) && continue
+                        (is_const_plus_var_known_linear(callee_result.total_incidence[callee_eq]) ||
+                         is_const_plus_var_known_linear(total_incidence[caller_eq])) && continue
                         @assert var !== unassigned
                         if !any(out->callee_eq in out[2], callee_var_schedule)
                             display(mss)
@@ -846,7 +846,7 @@ function tearing_schedule!(state::TransformationState, ci::CodeInstance, key::To
                     for var in callee_param_vars
                         varmap = info.mapping.var_coeffs[var]
                         nonlin = nothing
-                        if !is_const_plus_var_linear(varmap)
+                        if !is_const_plus_var_known_linear(varmap)
                             nonlin = schedule_nonlinear!(compact, key.param_vars, var_eq_matching, ir, ordinal, stmt.args[1+var], sicm_rename; vars=var_sols)
                         end
                         (argval, _) = schedule_incidence!(compact,
@@ -900,7 +900,7 @@ function tearing_schedule!(state::TransformationState, ci::CodeInstance, key::To
                 result.eq_callee_mapping[baseeq(result, structure, eq)] === nothing &&
                 eqclassification(result, structure, eq) != External &&
                 eqkind(result, structure, eq) == Intrinsics.Always &&
-                !is_const_plus_var_linear(total_incidence[eq])
+                !is_const_plus_var_known_linear(total_incidence[eq])
             push!(eq_orders[end], eq)
         end
     end
@@ -999,7 +999,7 @@ function tearing_schedule!(state::TransformationState, ci::CodeInstance, key::To
                 for var in callee_in_vars
                     nonlin = nothing
                     varmap = info.mapping.var_coeffs[var]
-                    if !is_const_plus_var_linear(varmap)
+                    if !is_const_plus_var_known_linear(varmap)
                         nonlin = schedule_nonlinear!(compact1, key.param_vars, var_eq_matching, ir, ordinal, eqinst[:stmt].args[1+var], ssa_rename; vars=var_sols, schedule_missing_var!)
                     end
                     (argval, _) = schedule_incidence!(compact1,
@@ -1043,7 +1043,7 @@ function tearing_schedule!(state::TransformationState, ci::CodeInstance, key::To
                 var = invview(var_eq_matching)[eq]
 
                 incT = state.total_incidence[eq]
-                anynonlinear = !is_const_plus_var_linear(incT)
+                anynonlinear = !is_const_plus_var_known_linear(incT)
                 nonlinearssa = nothing
                 if anynonlinear
                     if isa(var, Int) && isa(vars[var], SolvedVariable)
@@ -1085,7 +1085,7 @@ function tearing_schedule!(state::TransformationState, ci::CodeInstance, key::To
                 if isa(var, Int)
                     curval = nonlinearssa
                     (curval, thiscoeff) = schedule_incidence!(compact1, curval, incT, var, line; vars=var_sols, schedule_missing_var!)
-                    @assert thiscoeff != nonlinear
+                    @assert isa(thiscoeff, Float64)
                     curval = ir_mul_const!(compact1, line, 1/thiscoeff, curval)
                     var_sols[var] = isa(curval, SSAValue) ? CarriedSSAValue(ordinal, curval.id) : curval
                     insert_solved_var_here!(compact1, var, curval, line)
diff --git a/test/Project.toml b/test/Project.toml
index 486e03d..94f2f44 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -15,6 +15,7 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
 REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
diff --git a/test/incidence.jl b/test/incidence.jl
new file mode 100644
index 0000000..35ebeef
--- /dev/null
+++ b/test/incidence.jl
@@ -0,0 +1,344 @@
+module _Incidence
+
+using Test
+using DAECompiler
+using DAECompiler: Incidence, IncidenceValue, Linearity, refresh, linear, linear_time_dependent, linear_state_dependent, linear_time_and_state_dependent, nonlinear
+using DAECompiler.Intrinsics
+using Compiler: SSAValue, Const
+using SparseArrays: rowvals, nonzeros
+
+const *ᵢ = Core.Intrinsics.mul_float
+const +ᵢ = Core.Intrinsics.add_float
+const -ᵢ = Core.Intrinsics.sub_float
+
+walk(ex::Expr, inner, outer) = outer(Expr(ex.head, map(inner, ex.args)...))
+walk(ex, inner, outer) = outer(ex)
+
+postwalk(f, ex) = walk(ex, x -> postwalk(f, x), f)
+
+"""
+    @infer_incidence u₁ + u₂ + t
+    @infer_incidence begin
+      x = u₁ + u₂
+      x + u₃
+    end
+    @infer_incidence u₁ + u₂ true    # a second argument of `true` makes it return the IR, not the Incidence
+    @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()`
+in lexicographic order, so u₂ will appear before u₁. A special `t` variable may be
+used, which will be created as `sim_time()`.
+
+Note that if variable indices are not contiguous starting from 1, incidence annotations
+will not correspond to the input symbols. `@infer_incidence 1 + u₂` for example will be inferred
+to `Incidence(1.0 + u₁)`.
+"""
+macro infer_incidence(ex, return_ir = false)
+  variables = Symbol[]
+  uses_time = false
+  ex = postwalk(ex) do subex
+    # Force all observed function calls to be `@noinline`,
+    # the caller probably wants to test IPO.
+    Meta.isexpr(subex, :call) && return :(@noinline $subex)
+    isa(subex, Symbol) || return subex
+    uses_time |= subex === :t
+    startswith(string(subex), "u") || return subex
+    !in(subex, variables) && push!(variables, subex)
+    return subex
+  end
+  sort!(variables, by = string)
+  prelude = Expr(:block, [:($(esc(var)) = continuous()) for var in variables]...)
+  uses_time && pushfirst!(prelude.args, :($(esc(:t)) = sim_time()))
+  quote
+    function incidence()
+      $prelude
+      $(esc(ex))
+    end
+    ir = code_structure_by_type(Tuple{typeof(incidence)})
+    $return_ir ? ir : begin
+      ssa = SSAValue(length(ir.stmts) - 1)
+      ir[ssa][:type]::Incidence
+    end
+  end
+end
+
+dependencies(row) = sort(rowvals(row) .=> nonzeros(row), by = first)
+
+@testset "Incidence" begin
+  err = "state or time independence is not supported for nonlinearities"
+  @test_throws err Linearity(state_dependent = false, nonlinear = true)
+  @test_throws err Linearity(time_dependent = false, nonlinear = true)
+
+  @testset "Construction & printing" begin
+    incidence = Incidence(Const(1.0))
+    @test incidence.typ === Const(1.0)
+    @test dependencies(incidence.row) == []
+    @test repr(incidence) == "Incidence(1.0)"
+    @test incidence == incidence"1.0"
+
+    incidence = Incidence(Float64)
+    @test repr(incidence) == "Incidence(a)"
+    @test incidence == incidence"a"
+
+    incidence = Incidence(Float64, IncidenceValue[1.0])
+    @test dependencies(incidence.row) == [1 => 1]
+    @test repr(incidence) == "Incidence(a + t)"
+    @test incidence == incidence"a + t"
+
+    incidence = Incidence(String, IncidenceValue[1.0])
+    @test repr(incidence) == "Incidence(String, t)"
+    @test incidence == incidence"String, t"
+
+    incidence = Incidence(1)
+    @test incidence.typ === Const(0.0)
+    @test dependencies(incidence.row) == [2 => 1]
+    @test repr(incidence) == "Incidence(u₁)"
+    @test incidence == incidence"u₁"
+
+    incidence = Incidence(3)
+    @test dependencies(incidence.row) == [4 => 1]
+    @test repr(incidence) == "Incidence(u₃)"
+    @test incidence == incidence"u₃"
+
+    incidence = Incidence(Const(3.0), IncidenceValue[0.0, 0.0, 2.0, 1.0])
+    @test repr(incidence) == "Incidence(3.0 + 2.0u₂ + u₃)"
+    @test incidence == incidence"3.0 + 2.0u₂ + u₃"
+
+    incidence = Incidence(Const(0.0), IncidenceValue[4.0, 0.0, 2.0])
+    @test repr(incidence) == "Incidence(4.0t + 2.0u₂)"
+    @test incidence == incidence"4.0t + 2.0u₂"
+
+    incidence = Incidence(Const(0.0), IncidenceValue[nonlinear])
+    @test repr(incidence) == "Incidence(f(t))"
+    @test incidence == incidence"f(t)"
+
+    incidence = Incidence(Const(0.0), IncidenceValue[linear])
+    @test repr(incidence) == "Incidence(∝t)"
+    @test incidence == incidence"∝t"
+
+    incidence = Incidence(Const(0.0), IncidenceValue[1.0, nonlinear])
+    @test repr(incidence) == "Incidence(t + f(u₁))"
+    @test incidence == incidence"t + f(u₁)"
+
+    incidence = Incidence(Const(0.0), IncidenceValue[1.0, linear])
+    @test repr(incidence) == "Incidence(t + ∝u₁)"
+    @test incidence == incidence"t + ∝u₁"
+
+    incidence = Incidence(Const(0.0), IncidenceValue[linear, linear, linear])
+    @test repr(incidence) == "Incidence(∝t + ∝u₁ + ∝u₂)"
+    @test incidence == incidence"∝t + ∝u₁ + ∝u₂"
+
+    incidence = Incidence(Const(0.0), IncidenceValue[linear_state_dependent, linear_time_dependent, linear])
+    @test repr(incidence) == "Incidence(∝u₂ + f(∝ₛt, ∝ₜu₁))"
+    @test incidence == incidence"∝u₂ + f(∝ₛt, ∝ₜu₁)"
+    @test incidence == incidence"Incidence(∝u₂ + f(∝ₛt, ∝ₜu₁))"
+
+    incidence = Incidence(Const(0.0), IncidenceValue[linear_state_dependent, linear_time_dependent, nonlinear])
+    @test repr(incidence) == "Incidence(f(∝ₛt, ∝ₜu₁, u₂))"
+    @test incidence == incidence"f(∝ₛt, ∝ₜu₁, u₂)"
+
+    incidence = Incidence(Const(0.0), IncidenceValue[nonlinear, linear_time_dependent, nonlinear])
+    @test repr(incidence) == "Incidence(f(t, ∝ₜu₁, u₂))"
+    @test incidence == incidence"f(t, ∝ₜu₁, u₂)"
+
+    @test_throws "inconsistent with an absence of time incidence" Incidence(Const(0.0), IncidenceValue[0.0, linear_time_dependent])
+    @test_throws "inconsistent with an absence of state incidence" Incidence(Const(0.0), IncidenceValue[linear_state_dependent])
+    @test_throws "absence of state dependence for time" Incidence(Const(0.0), IncidenceValue[1.0, linear_time_dependent])
+    @test_throws "absence of state dependence for time" Incidence(Const(0.0), IncidenceValue[linear, linear_time_dependent])
+  end
+
+  @testset "Basic lattice operations with intrinsic tfuncs" begin
+    incidence = @infer_incidence t
+    @test incidence.typ === Const(0.0)
+    @test dependencies(incidence.row) == [1 => 1]
+    @test incidence == incidence"t"
+
+    incidence = @infer_incidence 5. +ᵢ u₁
+    @test incidence.typ === Const(5.0)
+    @test dependencies(incidence.row) == [2 => 1]
+    @test incidence == incidence"5.0 + u₁"
+
+    incidence = @infer_incidence 5.0 +ᵢ t +ᵢ u₁
+    @test incidence.typ === Const(5.0)
+    @test dependencies(incidence.row) == [1 => 1, 2 => 1]
+    @test incidence == incidence"5.0 + t + u₁"
+
+    incidence = @infer_incidence t *ᵢ u₁
+    @test incidence.typ === Const(0.0)
+    @test dependencies(incidence.row) == [1 => linear_state_dependent, 2 => linear_time_dependent]
+    @test incidence == incidence"f(∝ₛt, ∝ₜu₁)"
+
+    incidence = @infer_incidence u₁
+    @test incidence.typ === Const(0.0)
+    @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]
+    @test incidence == incidence"u₁ + u₂"
+
+    incidence = @infer_incidence u₁ *ᵢ u₂
+    @test incidence.typ === Const(0.0)
+    @test dependencies(incidence.row) == [2 => linear_state_dependent, 3 => linear_state_dependent]
+    @test incidence == incidence"f(∝ₛu₁, ∝ₛu₂)"
+
+    incidence = @infer_incidence (2.0 +ᵢ u₁) *ᵢ (3.0 +ᵢ u₂)
+    @test incidence.typ === Const(6.0)
+    @test dependencies(incidence.row) == [2 => linear_state_dependent, 3 => linear_state_dependent]
+    @test incidence == incidence"6.0 + f(∝ₛu₁, ∝ₛu₂)"
+
+    incidence = @infer_incidence (2.0 +ᵢ u₁) *ᵢ (3.0 +ᵢ u₁ *ᵢ u₂)
+    @test incidence.typ === Const(6.0)
+    @test dependencies(incidence.row) == [2 => nonlinear, 3 => linear_state_dependent]
+    @test incidence == incidence"6.0 + f(u₁, ∝ₛu₂)"
+
+    incidence = @infer_incidence (2.0 +ᵢ u₁) *ᵢ (3.0 +ᵢ u₁ *ᵢ u₂) +ᵢ u₃
+    @test incidence.typ === Const(6.0)
+    @test dependencies(incidence.row) == [2 => nonlinear, 3 => linear_state_dependent, 4 => 1.0]
+    @test incidence == incidence"6.0 + u₃ + f(u₁, ∝ₛu₂)"
+
+    incidence = @infer_incidence (u₁ +ᵢ u₂) *ᵢ t
+    @test incidence.typ === Const(0.0)
+    @test dependencies(incidence.row) == [1 => linear_state_dependent, (2:3 .=> linear_time_dependent)...]
+    @test incidence == incidence"f(∝ₛt, ∝ₜu₁, ∝ₜu₂)"
+
+    incidence = @infer_incidence u₁ *ᵢ u₂ +ᵢ (u₁ +ᵢ t) *ᵢ u₃
+    @test dependencies(incidence.row) == [(1:3 .=> linear_state_dependent)..., 4 => linear_time_and_state_dependent]
+    @test incidence == incidence"f(∝ₛt, ∝ₛu₁, ∝ₛu₂, ∝ₜₛu₃)"
+  end
+
+  @testset "IPO" begin
+    # NOTE: Most of additive terms (`.typ`) can't be precise given the current IPO representation.
+    # We expect `Const(0.0)` in most cases, but widen to `Float64`.
+
+    incidence = @infer_incidence 1.0t * 3.0
+    @test dependencies(incidence.row) == [1 => linear]
+    @test incidence == incidence"a + ∝t"
+
+    incidence = @infer_incidence (1.0 + u₁ +ᵢ u₂) * 1.0
+    @test incidence.typ === Float64
+    @test dependencies(incidence.row) == (2:3 .=> linear_state_dependent)
+    @test incidence == incidence"a + f(∝ₛu₁, ∝ₛu₂)"
+
+    incidence = @infer_incidence (2.0 + u₁) * (3.0 + u₂)
+    @test dependencies(incidence.row) == [2 => linear_state_dependent, 3 => linear_state_dependent]
+    @test incidence == incidence"a + f(∝ₛu₁, ∝ₛu₂)"
+
+    incidence = @infer_incidence 5.0 + u₁
+    @test incidence.typ === Const(5.0)
+    @test dependencies(incidence.row) == [2 => 1]
+    @test incidence == incidence"5.0 + u₁"
+
+    incidence = @infer_incidence u₁ * u₁
+    @test dependencies(incidence.row) == [2 => nonlinear]
+    @test incidence == incidence"f(u₁)"
+
+    incidence = @infer_incidence t * t
+    @test dependencies(incidence.row) == [1 => nonlinear]
+    @test incidence == incidence"f(t)"
+
+    mul3(a, b, c) = a *ᵢ (b *ᵢ c)
+    incidence = @infer_incidence mul3(t, u₁, u₂)
+    @test dependencies(incidence.row) == [1 => linear_state_dependent, (2:3 .=> linear_time_and_state_dependent)...]
+    @test incidence == incidence"f(∝ₛt, ∝ₜₛu₁, ∝ₜₛu₂)"
+
+    incidence = @infer_incidence mul3(t, u₁, u₁)
+    @test dependencies(incidence.row) == [1 => linear_state_dependent, 2 => nonlinear]
+    @test incidence == incidence"f(∝ₛt, u₁)"
+
+    incidence = @infer_incidence mul3(t, u₁, t)
+    # If we knew which state is used for state dependence,
+    # state should be inferred as linear_time_dependent.
+    @test dependencies(incidence.row) == [1 => nonlinear, 2 => linear_time_and_state_dependent]
+    @test incidence == incidence"f(t, ∝ₜₛu₁)"
+
+    incidence = @infer_incidence mul3(u₂, u₁, u₂)
+    @test dependencies(incidence.row) == [2 => linear_state_dependent, 3 => nonlinear]
+    @test incidence == incidence"f(∝ₛu₁, u₂)"
+
+    _muladd(a, b, c) = a +ᵢ b *ᵢ c
+    incidence = @infer_incidence _muladd(u₁, u₁, u₂)
+    # We widen to `nonlinear` because we can't yet infer that `b := u₁` is
+    # not multiplied by `a := u₁`. The solution would be to see that `a`
+    # is linear but state-independent and therefore can't be a factor of `b`.
+    @test dependencies(incidence.row) == [2 => nonlinear, 3 => linear_state_dependent]
+    @test incidence == incidence"f(u₁, ∝ₛu₂)"
+
+    # Here we still wouldn't be able to use the above solution because `a := u₁` is state-dependent.
+    # So `c := u₁` having a state-dependent coefficient might be multiplied by `a` a.k.a itself
+    # which would make it nonlinear, so IPO can only infer `u₁` as nonlinear.
+    _muladd2(a, b, c, d) = d *ᵢ a +ᵢ b *ᵢ c
+    incidence = @infer_incidence _muladd2(u₁, u₂, u₁, u₃)
+    @test dependencies(incidence.row) == [2 => nonlinear, 3 => linear_state_dependent, 4 => linear_state_dependent]
+    @test incidence == incidence"f(u₁, ∝ₛu₂, ∝ₛu₃)"
+
+    incidence = @infer_incidence 1/u₁
+    @test dependencies(incidence.row) == [2 => nonlinear]
+    @test incidence == incidence"a + f(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₁ * (1/t)
+    @test dependencies(incidence.row) == [1 => nonlinear, 2 => linear_time_dependent]
+    @test incidence == incidence"a + f(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 1/(u₁ * u₂)
+    @test dependencies(incidence.row) == [2 => nonlinear, 3 => nonlinear]
+    @test incidence == incidence"a + f(u₁, u₂)"
+  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
diff --git a/test/runtests.jl b/test/runtests.jl
index 2c79a91..3396e8e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,8 +1,13 @@
 include("basic.jl")
 include("reflection.jl")
+include("incidence.jl")
 include("ipo.jl")
 include("ssrm.jl")
 include("regression.jl")
 include("errors.jl")
 include("invalidation.jl")
-include("../benchmark/thermalfluid.jl")
\ No newline at end of file
+
+using Pkg
+Pkg.activate(joinpath(dirname(@__DIR__), "benchmark")) do
+  include("../benchmark/thermalfluid.jl")
+end