|
| 1 | +include(joinpath(@__DIR__, "..", "shared", "test_setup.jl")) |
| 2 | + |
| 3 | +if GROUP == "All" || GROUP == "Core" |
| 4 | + @testset "Computing variations around a sequence solution" begin |
| 5 | + # Computing sensitivities directly be explicitly writing down Lie derivatives |
| 6 | + function use_lie_derivatives( |
| 7 | + dds::StructuralIdentifiability.DDS{P}, |
| 8 | + params::Dict{P, T}, |
| 9 | + ic::Dict{P, T}, |
| 10 | + input_values::Dict{P, Array{T, 1}}, |
| 11 | + num_terms::Int, |
| 12 | + ) where {T <: Generic.FieldElem, P <: MPolyRingElem{T}} |
| 13 | + newvars = [var_to_str(v) for v in gens(parent(dds))] |
| 14 | + append!( |
| 15 | + newvars, |
| 16 | + [var_to_str(v) * "$i" for v in inputs(dds) for i in 1:num_terms], |
| 17 | + ) |
| 18 | + R, _ = StructuralIdentifiability.Nemo.polynomial_ring( |
| 19 | + base_ring(parent(dds)), |
| 20 | + newvars, |
| 21 | + ) |
| 22 | + explicit_sol = merge( |
| 23 | + Dict( |
| 24 | + parent_ring_change(x, R) => Vector{Any}([parent_ring_change(x, R)]) |
| 25 | + for (x, eq) in x_equations(dds) |
| 26 | + ), |
| 27 | + Dict( |
| 28 | + parent_ring_change(y, R) => |
| 29 | + Vector{Any}([parent_ring_change(eq, R)]) for |
| 30 | + (y, eq) in y_equations(dds) |
| 31 | + ), |
| 32 | + ) |
| 33 | + time_step = merge( |
| 34 | + Dict( |
| 35 | + parent_ring_change(x, R) => parent_ring_change(eq, R) for |
| 36 | + (x, eq) in x_equations(dds) |
| 37 | + ), |
| 38 | + Dict( |
| 39 | + parent_ring_change(p, R) => parent_ring_change(p, R) for |
| 40 | + p in StructuralIdentifiability.parameters(dds) |
| 41 | + ), |
| 42 | + Dict( |
| 43 | + parent_ring_change(u, R) => str_to_var(var_to_str(u) * "1", R) for |
| 44 | + u in inputs(dds) |
| 45 | + ), |
| 46 | + Dict( |
| 47 | + str_to_var(s * "$i", R) => str_to_var(s * "$(i + 1)", R) for |
| 48 | + s in map(var_to_str, inputs(dds)) for i in 1:(num_terms - 1) |
| 49 | + ), |
| 50 | + ) |
| 51 | + eval_dict = merge( |
| 52 | + Dict(parent_ring_change(p, R) => v for (p, v) in params), |
| 53 | + Dict(parent_ring_change(x, R) => v for (x, v) in ic), |
| 54 | + Dict(parent_ring_change(u, R) => val[1] for (u, val) in input_values), |
| 55 | + Dict( |
| 56 | + str_to_var(var_to_str(u) * "$i", R) => input_values[u][i + 1] for |
| 57 | + u in inputs(dds) for i in 1:(num_terms - 1) |
| 58 | + ), |
| 59 | + ) |
| 60 | + generalized_parameters = [ |
| 61 | + parent_ring_change(p, R) for |
| 62 | + p in vcat(x_vars(dds), StructuralIdentifiability.parameters(dds)) |
| 63 | + ] |
| 64 | + for i in 2:num_terms |
| 65 | + for (k, v) in explicit_sol |
| 66 | + push!(explicit_sol[k], eval_at_dict(v[end], time_step)) |
| 67 | + end |
| 68 | + end |
| 69 | + part_diffs = Dict( |
| 70 | + (f, p) => [] for f in keys(explicit_sol) for p in generalized_parameters |
| 71 | + ) |
| 72 | + for i in 1:num_terms |
| 73 | + for (f, ders) in explicit_sol |
| 74 | + for p in generalized_parameters |
| 75 | + push!( |
| 76 | + part_diffs[(f, p)], |
| 77 | + eval_at_dict(derivative(ders[i], p), eval_dict), |
| 78 | + ) |
| 79 | + end |
| 80 | + end |
| 81 | + end |
| 82 | + return Dict( |
| 83 | + ( |
| 84 | + parent_ring_change(k[1], parent(dds)), |
| 85 | + parent_ring_change(k[2], parent(dds)), |
| 86 | + ) => res for (k, res) in part_diffs |
| 87 | + ) |
| 88 | + end |
| 89 | + |
| 90 | + locQQ = StructuralIdentifiability.Nemo.QQ |
| 91 | + |
| 92 | + dds = @DDSmodel(a(t + 1) = a(t)^2 + b, y(t) = 1 / (a(t) * c(t))) |
| 93 | + params = Dict(b => locQQ(1)) |
| 94 | + ic = Dict(a => locQQ(2)) |
| 95 | + input_values = Dict(c => [locQQ(1), locQQ(-2), locQQ(3), locQQ(-4), locQQ(5)]) |
| 96 | + seq_sol, diff_sol = |
| 97 | + differentiate_sequence_solution(dds, params, ic, input_values, 4) |
| 98 | + diff_y = differentiate_sequence_output(dds, params, ic, input_values, 4) |
| 99 | + lie_ders_sol = use_lie_derivatives(dds, params, ic, input_values, 4) |
| 100 | + merged = merge(diff_sol, diff_y) |
| 101 | + @test merged == lie_ders_sol |
| 102 | + |
| 103 | + dds = @DDSmodel( |
| 104 | + a(t + 1) = |
| 105 | + (23 * k1 * a(t) - 7 * b(t)^3) // (a(t)^2 + b(t)^2) - c(t)^3 * k1 * b(t), |
| 106 | + b(t + 1) = a(t) + 17 * (b(t) - c(t))^2 + 1 // (a(t) + b(t) - k2), |
| 107 | + y(t) = (a(t) + b(t) - c(t)) // (a(t)^2 + k2) |
| 108 | + ) |
| 109 | + params = Dict(k1 => locQQ(1), k2 => locQQ(2)) |
| 110 | + ic = Dict(a => locQQ(3), b => locQQ(-4)) |
| 111 | + input_values = Dict(c => [locQQ(5), locQQ(-6), locQQ(7), locQQ(-8)]) |
| 112 | + seq_sol, diff_sol = |
| 113 | + differentiate_sequence_solution(dds, params, ic, input_values, 2) |
| 114 | + diff_y = differentiate_sequence_output(dds, params, ic, input_values, 2) |
| 115 | + lie_ders_sol = use_lie_derivatives(dds, params, ic, input_values, 2) |
| 116 | + merged = merge(diff_sol, diff_y) |
| 117 | + @test merged == lie_ders_sol |
| 118 | + end |
| 119 | +end |
0 commit comments