WIP: Output saturation

benchmarking/benchmarks.jl

@@ -1279,6 +1279,18 @@ benchmarks = Dict(
             y(t) = D(t)
         )
     ),
+    # equations (1) - (5) from https://doi.org/10.3389/fimmu.2020.01376
+    :TranAlRadhawi => Dict(
+        :name => "Chemotherapy model",
+        :ode => @ODEmodel(
+            C'(t) = u(t) - k1 * C(t) / (k2 + C(t)),
+            T'(t) = ka * T(t) - kb * C(t) * T(t) / (kc + T(t)) - kd * T(t) * I(t),
+            I'(t) = X(t) - ke * T(t) * I(t) - kf * C(t) * I(t) - kg * I(t) * Y(t) - kh * I(t),
+            X'(t) = C(t) / (1 + C(t) * ki_inv) - kj * X(t) - kk * X(t) * Y(t),
+            Y'(t) = I(t) / (1 + C(t) * kl_inv) - km * Y(t) * C(t),
+            y(t) = T(t)
+        )
+    )
 )
 
 # the NFkB example

src/StructuralIdentifiability.jl

@@ -75,6 +75,7 @@ include("elimination.jl")
 include("primality_check.jl")
 include("io_equation.jl")
 include("states.jl")
+include("output_saturation.jl")
 include("global_identifiability.jl")
 include("identifiable_functions.jl")
 include("parametrizations.jl")

src/output_saturation.jl

@@ -0,0 +1,93 @@
+"""
+    saturate_outputs(ode, orders, max_deg = 5)
+
+Adds extra outputs extracted from the Lie derivatives of the existing outputs up to the prescribed order.
+Adds only the outputs which are algebraically independent over the existing outputs, parameters, and inputs
+(and, thus, may be useful in the differential elimination process).
+
+Inputs:
+- `ode` - ode model
+- `orders` - dictionary from outputs to the order of the Lie derivative to consider
+- `max_deg` - the maximal degree (as a rational function) of possible new outputs to consider
+
+Output: a new ode model with inputs inserted
+"""
+function saturate_outputs(
+        ode::ODE{P},
+        orders::Dict{P, Int},
+        max_deg = 5,
+    ) where {P <: MPolyRingElem}
+    lie_ders = lie_derivatives_up_to(ode, Dict(ode.y_equations[y] => ord for (y, ord) in orders))
+    @info "Degrees: $([(total_degree(numerator(f)), total_degree(denominator(f))) for f in lie_ders])"
+
+    current_y = RationalFunctionField(vcat(collect(values(ode.y_equations)), ode.u_vars, ode.parameters))
+    lie_ders = filter(f -> total_degree_frac(f) <= max_deg, lie_ders)
+    sort!(lie_ders, lt = RationalFunctionFields.rational_function_cmp)
+    for f in lie_ders
+        if !first(RationalFunctionFields.check_algebraicity_modp(current_y, [f]))
+            current_y = RationalFunctionField([generators(current_y)..., f])
+        end
+    end
+    new_outputs = generators(current_y)[(length(ode.y_vars) + length(ode.u_vars) + length(ode.parameters) + 1):end]
+
+    idx = 1
+    old_y_names = map(var_to_str, ode.y_vars)
+    new_ys = Dict{String, Generic.FracFieldElem{P}}()
+    @info "New outputs $new_outputs"
+    for new_y in new_outputs
+        while "y_aux_$idx(t)" in old_y_names
+            idx += 1
+        end
+        new_ys["y_aux_$idx(t)"] = new_y
+        idx += 1
+    end
+    @info new_ys
+
+    return add_outputs(ode, new_ys)
+end
+
+#------------------------------------------------------------------------------
+"""
+    propose_orders(ode)
+
+For a given ode, proposes (using a heuristic) the orders for the outputs to consider in the saturation process.
+The order is defined a the highest order of differentiation at which a new output may occur plus two.
+
+Inputs:
+- `ode` - ode model
+
+Output: dictionary from outputs to orders
+"""
+function propose_orders(ode)
+    x_with_u = filter(x -> length(intersect(vars(ode.x_equations[x]), ode.u_vars)) > 0, ode.x_vars)
+
+    graph = Dict(x => Set(intersect(vars(ode.x_equations[x]), ode.x_vars)) for x in ode.x_vars)
+
+    result = Dict(y => 0 for y in ode.y_vars)
+    for y in ode.y_vars
+        current_x = Set(intersect(vars(ode.y_equations[y]), ode.x_vars))
+        reached_u = length(intersect(current_x, x_with_u))
+        if reached_u > 0
+            result[y] = 2
+        end
+        step = 1
+        while true
+            new_x = copy(current_x)
+            for x in current_x
+                new_x = union(new_x, graph[x])
+            end
+            if current_x == new_x
+                break
+            end
+            current_x = new_x
+            new_reached = length(intersect(current_x, x_with_u))
+            if new_reached > reached_u
+                result[y] = step + 2
+            end
+            reached_u = new_reached
+            step += 1
+        end
+    end
+
+    return result
+end

src/primality_check.jl

@@ -6,7 +6,6 @@ function check_primality_zerodim(J::Array{QQMPolyRingElem, 1})
     dim = length(basis)
     S = Nemo.matrix_space(Nemo.QQ, dim, dim)
     matrices = []
-    @debug "$J $basis"
     @debug "Dim is $dim"
     for v in gens(parent(first(J)))
         M = zero(S)
@@ -17,15 +16,31 @@ function check_primality_zerodim(J::Array{QQMPolyRingElem, 1})
             end
         end
         push!(matrices, M)
-        @debug "Multiplication by $v: $M"
     end
     generic_multiplication = sum(Nemo.QQ(rand(1:100)) * M for M in matrices)
-    @debug generic_multiplication
+    @debug "Generic multiplication matrix computed"
+    @debug "Trying reductions over primed first"
+    NUM_PRIMES = 10
+    p = 2^31 - 1
+    for _ in 1:NUM_PRIMES
+        @debug "Prime is $p"
+        F = Nemo.GF(p)
+        S_F = Nemo.matrix_space(F, dim, dim)
+        generic_multiplication_modp = S_F([generic_multiplication[i, j] for i in 1:dim for j in 1:dim])
+        R, t = Nemo.polynomial_ring(F, "t")
+        chpoly = Nemo.charpoly(R, generic_multiplication_modp)
+        @debug "Charpoly computed"
+        if Nemo.is_irreducible(chpoly)
+            return true
+        end
+        p = Primes.nextprime(p + 1)
+    end
 
+    @debug "No conclusion modulo primes, computing over Q"
     R, t = Nemo.polynomial_ring(Nemo.QQ, "t")
-    @debug "$(Nemo.charpoly(R, generic_multiplication))"
+    chpoly = Nemo.charpoly(R, generic_multiplication)
 
-    return Nemo.is_irreducible(Nemo.charpoly(R, generic_multiplication))
+    return Nemo.is_irreducible(chpoly)
 end
 
 #------------------------------------------------------------------------------

src/states.jl

@@ -67,6 +67,28 @@ function lie_derivative(f::P, ode::ODE{<:P}) where {P <: MPolyRingElem{<:FieldEl
     return lie_derivative(f // 1, ode)
 end
 
+#------------------------------------------------------------------------------
+
+function lie_derivatives_up_to(ode::ODE, orders::Dict{P, Int}) where {P <: AbstractAlgebra.RingElem}
+    result = Array{Generic.FracFieldElem, 1}()
+    for (f, ord) in orders
+        curr = extract_coefficients_ratfunc(f, ode.u_vars)
+        for _ in 0:ord
+            filter!(!is_rational_func_const, curr)
+            append!(result, curr)
+            curr = reduce(
+                vcat,
+                [
+                    extract_coefficients_ratfunc(lie_derivative(g, ode), ode.u_vars) for
+                        g in curr
+                ],
+                init = empty(curr),
+            )
+        end
+    end
+    return result
+end
+
 #------------------------------------------------------------------------------
 """
     states_generators(ode, io_equations)
@@ -95,22 +117,5 @@ identifiable functions of parameters only
         end
     end
 
-    result = Array{Generic.FracFieldElem{P}, 1}()
-    for (y, ord) in y_to_ord
-        curr = extract_coefficients_ratfunc(ode.y_equations[y], ode.u_vars)
-        for _ in 0:ord
-            filter!(!is_rational_func_const, curr)
-            append!(result, curr)
-            curr = reduce(
-                vcat,
-                [
-                    extract_coefficients_ratfunc(lie_derivative(f, ode), ode.u_vars) for
-                        f in curr
-                ],
-                init = empty(curr),
-            )
-        end
-    end
-
-    return result
+    return lie_derivatives_up_to(ode, Dict(ode.y_equations[y] => ord for (y, ord) in y_to_ord))
 end

src/submodels.jl

@@ -53,6 +53,7 @@ function traverse_outputs(
     raw_models = Dict{QQMPolyRingElem, Set{QQMPolyRingElem}}()
     for y in ys
         model = dfs(graph, y, Set{QQMPolyRingElem}())
+        delete!(model, y)
         raw_models[y] = model
     end
     return raw_models
@@ -79,11 +80,14 @@ end
 
 # ------------------------------------------------------------------------------
 
-function search_add_unions(submodels::Array{Set{QQMPolyRingElem}, 1})
+function search_add_unions(submodels::Set{Set{QQMPolyRingElem}})
     result = Array{Set{QQMPolyRingElem}, 1}([Set{QQMPolyRingElem}()])
     for model in submodels
         for index in 1:length(result)
-            push!(result, union(result[index], model))
+            uni = union(result[index], model)
+            if !(uni in result)
+                push!(result, uni)
+            end
         end
     end
     return result
@@ -200,8 +204,9 @@ function find_submodels(ode::ODE{P}) where {P <: MPolyRingElem}
     ys = ode.y_vars
     xs = ode.x_vars
     raw_models = traverse_outputs(graph, ys)
-    input_unions = [raw_models[y] for y in ys]
-    unions = (search_add_unions(input_unions))
+    input_unions = Set([raw_models[y] for y in ys])
+    unions = search_add_unions(input_unions)
+    @debug "We have $(length(unions)) candidate submodels to check"
     saturate_ys(unions, ys, graph, xs)
     result = filter_max_empty(ode, union(unions))
     back2ode = submodel2ode(ode, result)

test/check_primality_zerodim.jl

@@ -9,5 +9,8 @@ if GROUP == "All" || GROUP == "Core"
         @test check_primality_zerodim([x^3 - 5, y - 1]) == true
 
         @test check_primality_zerodim([x^2 + 1, y^3 - 3 * x + x + 5]) == true
+
+        # not prime over any modulous but prime over Q
+        @test check_primality_zerodim([x, y^4 + 1]) == true
     end
 end

test/runtests.jl

@@ -62,7 +62,9 @@ using StructuralIdentifiability:
     x_equations,
     y_equations,
     inputs,
-    quotient_basis
+    quotient_basis,
+    propose_orders,
+    saturate_outputs
 
 const GROUP = get(ENV, "GROUP", "All")
 

test/y_saturation.jl

@@ -0,0 +1,72 @@
+if GROUP == "All" || GROUP == "Core"
+    @testset "Output saturation" begin
+        # full nfkb model from benchmarks
+        nfkb = @ODEmodel(
+            x1'(t) = k_prod - k_deg * x1(t) - k1 * x1(t) * u(t),
+            x2'(t) =
+                -k3 * x2(t) - k_deg * x2(t) - a2 * x2(t) * x10(t) + t1 * x4(t) -
+                a3 * x2(t) * x13(t) +
+                t2 * x5(t) +
+                (k1 * x1(t) - k2 * x2(t) * x8(t)) * u(t),
+            x3'(t) = k3 * x2(t) - k_deg * x3(t) + k2 * x2(t) * x8(t) * u(t),
+            x4'(t) = a2 * x2(t) * x10(t) - t1 * x4(t),
+            x5'(t) = a3 * x2(t) * x13(t) - t2 * x5(t),
+            x6'(t) = c6a * x13(t) - a1 * x6(t) * x10(t) + t2 * x5(t) - i1 * x6(t),
+            x7'(t) = i1 * kv * x6(t) - a1 * x11(t) * x7(t),
+            x8'(t) = c4 * x9(t) - c5 * x8(t),
+            x9'(t) = c2 - c1 * x7(t) - c3 * x9(t),
+            x10'(t) =
+                -a2 * x2(t) * x10(t) - a1 * x10(t) * x6(t) + c4a * x12(t) - c5a * x10(t) -
+                i1a * x10(t) + e1a * x11(t),
+            x11'(t) = -a1 * x11(t) * x7(t) + i1a * kv * x10(t) - e1a * kv * x11(t),
+            x12'(t) = c2a + c1a * x7(t) - c3a * x12(t),
+            x13'(t) = a1 * x10(t) * x6(t) - c6a * x13(t) - a3 * x2(t) * x13(t) + e2a * x14(t),
+            x14'(t) = a1 * x11(t) * x7(t) - e2a * kv * x14(t),
+            x15'(t) = c2c + c1c * x7(t) - c3c * x15(t),
+            y1(t) = x7(t),
+            y2(t) = x10(t) + x13(t),
+            y3(t) = x9(t),
+            y4(t) = x1(t) + x2(t) + x3(t),
+            y5(t) = x2(t),
+            y6(t) = x12(t)
+        )
+
+        nfkb_orders = propose_orders(nfkb)
+        new_nfkb = saturate_outputs(nfkb, nfkb_orders)
+        @test length(new_nfkb.y_vars) == 14
+        nfkb_time = @elapsed nfkb_ident_results = assess_identifiability(new_nfkb)
+        @test nfkb_time < 500
+        @test count(x -> x == :globally, values(nfkb_ident_results)) == 37
+        @test count(x -> x == :nonidentifiable, values(nfkb_ident_results)) == 7
+
+        # TumorPilis model from the benchmarks
+        tum_pil = @ODEmodel(
+            T'(t) =
+                a * T(t) * (1 - b * T(t)) - c1 * N(t) * T(t) - D(t) * T(t) -
+                KT * M(t) * T(t), #tumor cells
+            L'(t) =
+                -m * L(t) - q * L(t) * T(t) - ucte * L(t)^2 +
+                r2 * C(t) * T(t) +
+                pI * L(t) * I(t) / (gI + I(t)) +
+                u1(t) - KL * M(t) * L(t), # tumor-specific effector cells, T-celss
+            N'(t) =
+                alpha1 - f * N(t) + g * T(t) * N(t) / (h + T(t)) - p * N(t) * T(t) - KN * M(t) * N(t), # non-specific effector cells, NK cells
+            C'(t) = alpha2 - beta * C(t) - KC * M(t) * C(t), #circulating lymphocytes
+            I'(t) =
+                pt * T(t) * L(t) / (gt + T(t)) + w * L(t) * I(t) - muI * I(t) + u2(t), # IL-2, VI = u2 aplicación directa, terapia de IL2
+            M'(t) = -gamma * M(t) + u1(t), #chemotherapy drug, terapia/aplicación de quimio, u1 = VM
+            y1(t) = L(t),
+            y2(t) = N(t),
+            y3(t) = M(t)
+        )
+
+        tum_pil_orders = propose_orders(tum_pil)
+        new_tum_pil = saturate_outputs(tum_pil, tum_pil_orders)
+        @test length(new_tum_pil.y_vars) == 6
+        tum_pil_time = @elapsed tum_pil_results = assess_identifiability(new_tum_pil)
+        @test tum_pil_time < 500
+        @test count(x -> x == :globally, values(tum_pil_results)) == 22
+        @test count(x -> x == :nonidentifiable, values(tum_pil_results)) == 9
+
+    end
+end