diff --git a/benchmarking/IdentifiableFunctions/experiments.jl b/benchmarking/IdentifiableFunctions/experiments.jl index 72dada132..a12876571 100644 --- a/benchmarking/IdentifiableFunctions/experiments.jl +++ b/benchmarking/IdentifiableFunctions/experiments.jl @@ -988,14 +988,15 @@ begin ) using Nemo, Logging - using JuliaInterpreter Groebner = StructuralIdentifiability.Groebner # ParamPunPam = StructuralIdentifiability.ParamPunPam Base.global_logger(ConsoleLogger(Logging.Info)) end -dennums, ring = - StructuralIdentifiability.initial_identifiable_functions(Pivastatin, p = 0.99); +res = StructuralIdentifiability.find_identifiable_functions(CGV1990); + +res = StructuralIdentifiability.assess_identifiability(CGV1990); + fracs = StructuralIdentifiability.dennums_to_fractions(dennums); rff = StructuralIdentifiability.RationalFunctionField(dennums); diff --git a/src/RationalFunctionFields/IdealMQS.jl b/src/RationalFunctionFields/IdealMQS.jl index bd4d7aabd..4864ffd88 100644 --- a/src/RationalFunctionFields/IdealMQS.jl +++ b/src/RationalFunctionFields/IdealMQS.jl @@ -25,14 +25,16 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal # NB: this lists may have different length nums_qq::Vector{T} dens_qq::Vector{T} - sat_qq::T + sat_qq::Vector{T} # dens_indices[i] is a pair of form (from, to). # Denominator as index i corresponds to numerators at indices [from..to] dens_indices::Vector{Tuple{Int, Int}} # Indices of pivot elements for each component of funcs_den_nums pivots_indices::Vector{Int} den_lcm::T + # Saturation sat_var_index::Int + sat_var_count::Int # Numerators and denominators over GF. # We cache them and maintain a map: # a finite field --> an image over this finite field @@ -46,16 +48,26 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal """ IdealMQS(funcs_den_nums::Vector{Vector}) - Given an array of polynomials, that is generators of form + Given an array of polynomials, that is generators of the form [[f1, f2, f3, ...], [g1, g2, g3, ...], ...] - constructs an MQS ideal. + constructs the corresponding MQS ideal. + + ## Options + + This constructor takes the following optional arguments: + - `sat_varname`: a string, name of the saturation variable. + - `sat_var_position`: a symbol, position of the saturation variable in the + monomial ordering. Possible options are `:first` and `:last`. + - `ordering`: a symbol, monomial ordering. + - `sat_factorize`: a bool, if the saturating `Q(y)` should be factored. """ function IdealMQS( funcs_den_nums::Vector{Vector{PolyQQ}}; sat_varname = "t", sat_var_position = :first, + sat_factorize = true, ordering = :degrevlex, ) where {PolyQQ} # We are given polynomials of form @@ -102,7 +114,6 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal @debug "Saturating variable is $sat_varname, index is $sat_var_index" flush(stdout) R_sat, v_sat = Nemo.PolynomialRing(K, varnames, ordering = ordering) - # Saturation t_sat = v_sat[sat_var_index] den_lcm_orig = den_lcm den_lcm = parent_ring_change( @@ -111,8 +122,32 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal matching = :byindex, shift = Int(sat_var_index == 1), ) + # Construct the saturating polynomial Q(x) * t - 1. + # If `sat_factorize` is true, then factorize Q(x). den_lcm_sat = parent_ring_change(den_lcm, R_sat) - sat_qq = den_lcm_sat * t_sat - 1 + den_lcm_sat = divexact(den_lcm_sat, leading_coefficient(den_lcm_sat)) + @assert isone(leading_coefficient(den_lcm_sat)) + sat_qq = Vector{PolyQQ}() + # TODO: in addition to these cases, also do not factor the saturation + # polynomial when the ideal contains a single polynomial plus saturaton + # (this case is the most frequent one) + if !sat_factorize || isone(den_lcm_sat) + push!(sat_qq, den_lcm_sat * t_sat - 1) + else + fact_qq = Nemo.factor(den_lcm_sat) + # Potentially extend the ring with more variables for saturation + sat_varnames = map(i -> "$sat_varname$i", 1:length(fact_qq)) + varnames = append_at_index(ystrs, sat_var_index, sat_varnames) + @debug "Saturating variables are $sat_varnames, index is $sat_var_index" + R_sat, v_sat = Nemo.PolynomialRing(K, varnames, ordering = ordering) + t_sat_many = v_sat[sat_var_index:(sat_var_index + length(sat_varnames) - 1)] + for ((sat_qq_factor, deg), sat_var_i) in zip(fact_qq, t_sat_many) + sat_qq_factor = parent_ring_change(sat_qq_factor, R_sat) + push!(sat_qq, sat_qq_factor^deg * sat_var_i - 1) + end + end + sat_var_count = length(sat_qq) + @info "Degrees of saturating polynomials are $(map(total_degree, sat_qq))" # We construct the array of numerators nums_qq and the array of # denominators dens_qq. Since there are usually more unique numerators # than denominators, we condense the array dens_qq and produce an @@ -129,7 +164,7 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal den, R_sat, matching = :byindex, - shift = Int(sat_var_index == 1), + shift = sat_var_count * (sat_var_index == 1), ) push!(dens_qq, den) push!(dens_indices, (length(nums_qq) + 1, length(nums_qq) + length(plist) - 1)) @@ -140,16 +175,17 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal num, R_sat, matching = :byindex, - shift = Int(sat_var_index == 1), + shift = sat_var_count * (sat_var_index == 1), ) push!(nums_qq, num) end end parent_ring_param, _ = PolynomialRing(ring, varnames, ordering = ordering) - @debug "Constructed MQS ideal in $R_sat with $(length(nums_qq) + 1) elements" + @debug "Constructed MQS ideal in $R_sat with $(length(nums_qq)) + $(length(sat_qq)) elements" flush(stdout) @assert length(pivots_indices) == length(dens_indices) == length(dens_qq) @assert length(pivots_indices) == length(funcs_den_nums) + @assert length(sat_qq) == sat_var_count && length(sat_qq) > 0 new{elem_type(R_sat)}( funcs_den_nums, @@ -162,6 +198,7 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal pivots_indices, den_lcm, sat_var_index, + sat_var_count, Dict(), Dict(), Dict(), @@ -170,7 +207,7 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal end end -Base.length(ideal::IdealMQS) = length(ideal.nums_qq) + 1 +Base.length(ideal::IdealMQS) = length(ideal.nums_qq) + length(ideal.sat_qq) AbstractAlgebra.base_ring(ideal::IdealMQS) = base_ring(ideal.nums_qq[1]) AbstractAlgebra.parent(ideal::IdealMQS) = ideal.parent_ring_param ParamPunPam.parent_params(ideal::IdealMQS) = base_ring(ideal.parent_ring_param) @@ -191,7 +228,8 @@ end end function fractionfree_generators_raw(mqs::IdealMQS) - # TODO: this assumes mqs.sat_var_index is last, and thus is broken + # TODO: this assumes mqs.sat_var_index is last, and thus is broken. + # TODO: well, now this is totally broken. ring_params = ParamPunPam.parent_params(mqs) K = base_ring(ring_params) varnames = map(string, Nemo.symbols(ring_params)) @@ -239,8 +277,8 @@ function ParamPunPam.reduce_mod_p!( return nothing end nums_qq, dens_qq, sat_qq = mqs.nums_qq, mqs.dens_qq, mqs.sat_qq - sat_gf = map_coefficients(c -> ff(c), sat_qq) - ring_ff = parent(sat_gf) + sat_gf = map(poly -> map_coefficients(c -> ff(c), poly), sat_qq) + ring_ff = parent(first(sat_gf)) nums_gf = map(poly -> map_coefficients(c -> ff(c), poly, parent = ring_ff), nums_qq) dens_gf = map(poly -> map_coefficients(c -> ff(c), poly, parent = ring_ff), dens_qq) mqs.cached_nums_gf[ff] = nums_gf @@ -263,12 +301,13 @@ function ParamPunPam.specialize_mod_p( K_2 = base_ring(nums_gf[1]) @assert K_1 == K_2 @assert length(point) == nvars(ParamPunPam.parent_params(mqs)) - # +1 actual variable because of the saturation! - @assert length(point) + 1 == nvars(parent(nums_gf[1])) - point_sat = append_at_index(point, mqs.sat_var_index, one(K_1)) + # +k actual variables because of the saturation! + @assert length(point) + mqs.sat_var_count == nvars(parent(nums_gf[1])) + point_sat = + append_at_index(point, mqs.sat_var_index, map(_ -> one(K_1), 1:(mqs.sat_var_count))) nums_gf_spec = map(num -> evaluate(num, point_sat), nums_gf) dens_gf_spec = map(den -> evaluate(den, point_sat), dens_gf) - polys = Vector{typeof(sat_gf)}(undef, length(nums_gf_spec) + 1) + polys = Vector{eltype(sat_gf)}(undef, length(nums_gf_spec) + length(sat_gf)) @inbounds for i in 1:length(dens_gf_spec) den, den_spec = dens_gf[i], dens_gf_spec[i] iszero(den_spec) && __throw_unlucky_evaluation("Point: $point") @@ -278,9 +317,12 @@ function ParamPunPam.specialize_mod_p( polys[j] = num * den_spec - den * num_spec end end - polys[end] = sat_gf - if !saturated - resize!(polys, length(polys) - 1) + if saturated + for i in 1:length(sat_gf) + polys[length(nums_gf_spec) + i] = sat_gf[i] + end + else + resize!(polys, length(nums_gf_spec)) end return polys end @@ -291,12 +333,13 @@ function specialize(mqs::IdealMQS, point::Vector{Nemo.fmpq}; saturated = true) dens_indices = mqs.dens_indices K = base_ring(mqs) @assert length(point) == nvars(ParamPunPam.parent_params(mqs)) - # +1 actual variable because of the saturation! - @assert length(point) + 1 == nvars(parent(nums_qq[1])) - point_sat = append_at_index(point, mqs.sat_var_index, one(K)) + # +k actual variables because of the saturation! + @assert length(point) + mqs.sat_var_count == nvars(parent(nums_qq[1])) + point_sat = + append_at_index(point, mqs.sat_var_index, map(_ -> one(K), 1:(mqs.sat_var_count))) nums_qq_spec = map(num -> evaluate(num, point_sat), nums_qq) dens_qq_spec = map(den -> evaluate(den, point_sat), dens_qq) - polys = Vector{typeof(sat_qq)}(undef, length(nums_qq_spec) + 1) + polys = Vector{eltype(sat_qq)}(undef, length(nums_qq_spec) + length(sat_qq)) @inbounds for i in 1:length(dens_qq_spec) den, den_spec = dens_qq[i], dens_qq_spec[i] iszero(den_spec) && __throw_unlucky_evaluation("Point: $point") @@ -306,9 +349,12 @@ function specialize(mqs::IdealMQS, point::Vector{Nemo.fmpq}; saturated = true) polys[j] = num * den_spec - den * num_spec end end - polys[end] = sat_qq - if !saturated - resize!(polys, length(polys) - 1) + if saturated + for i in 1:length(sat_qq) + polys[length(nums_qq_spec) + i] = sat_qq[i] + end + else + resize!(polys, length(nums_qq_spec)) end return polys end diff --git a/src/RationalFunctionFields/normalforms.jl b/src/RationalFunctionFields/normalforms.jl index 02aafcbbd..f2bc58332 100644 --- a/src/RationalFunctionFields/normalforms.jl +++ b/src/RationalFunctionFields/normalforms.jl @@ -54,14 +54,20 @@ function local_normal_forms( ) @assert !isempty(point) @assert parent(first(point)) == finite_field - point_ff_ext = append_at_index(point, mqs.sat_var_index, one(finite_field)) + # TODO: it is possible that the saturation details should be hidden in the + # MQS ideal and not visible to the outside code + point_ff_ext = append_at_index( + point, + mqs.sat_var_index, + map(_ -> one(finite_field), 1:(mqs.sat_var_count)), + ) gens_ff_spec = specialize_mod_p(mqs, point) gb_ff_spec = Groebner.groebner(gens_ff_spec) ring_ff = parent(gb_ff_spec[1]) xs_ff = gens(ring_ff) normal_forms_ff = Vector{elem_type(ring_ff)}(undef, 0) monoms_ff = Vector{elem_type(ring_ff)}(undef, 0) - xs_ff = cut_at_index(xs_ff, mqs.sat_var_index) + xs_ff = cut_at_index(xs_ff, mqs.sat_var_index, mqs.sat_var_count) pivot_vectors = map(f -> exponent_vector(f, 1), xs_ff) @debug """ variables in the finite field: $(xs_ff) @@ -356,7 +362,11 @@ function linear_relations_between_normal_forms( n_relations_ff = length(complete_intersection_relations_ff) # Filter out some relations from the complete intersection zeroed_relations_inds = Vector{Int}() - point_ext = append_at_index(point, mqs.sat_var_index, one(finite_field)) + point_ext = append_at_index( + point, + mqs.sat_var_index, + map(_ -> one(finite_field), 1:(mqs.sat_var_count)), + ) for i in 1:length(complete_intersection_relations_ff) relation = complete_intersection_relations_ff[i] relation_mqs = relation - evaluate(relation, point_ext) @@ -408,7 +418,11 @@ function linear_relations_between_normal_forms( end relation_qq_param = evaluate( relation_qq, - append_at_index(xs_param, mqs.sat_var_index, one(ring_param)), + append_at_index( + xs_param, + mqs.sat_var_index, + map(_ -> one(ring_param), 1:(mqs.sat_var_count)), + ), ) relations_qq[i] = relation_qq_param // one(relation_qq_param) end diff --git a/src/util.jl b/src/util.jl index af9d95c2d..6c3d47063 100644 --- a/src/util.jl +++ b/src/util.jl @@ -10,7 +10,7 @@ end # ------------------------------------------------------------------------------ -function append_at_index(vec::Vector{T}, idx::Integer, val::T) where {T} +function append_at_index(vec::Vector{T}, idx::Integer, val) where {T} # NOTE: could also just use insert! @assert (idx == 1) || (idx == length(vec) + 1) if idx == 1 @@ -20,12 +20,12 @@ function append_at_index(vec::Vector{T}, idx::Integer, val::T) where {T} end end -function cut_at_index(vec::Vector{T}, idx::Integer) where {T} +function cut_at_index(vec::Vector{T}, idx::Integer, count::Integer) where {T} @assert (idx == 1) || (idx == length(vec)) if idx == 1 - vec[2:end] + vec[(count + 1):end] else - vec[1:(end - 1)] + vec[1:(end - count)] end end diff --git a/test/identifiable_functions.jl b/test/identifiable_functions.jl index 6adcd51ca..dab6aebc4 100644 --- a/test/identifiable_functions.jl +++ b/test/identifiable_functions.jl @@ -935,7 +935,7 @@ ode = StructuralIdentifiability.@ODEmodel( ident_funcs = [k3, k1 // k7, k5 // k2, k6 // k2, k7 * EpoR_A] push!(test_cases, (ode = ode, ident_funcs = ident_funcs)) -ode = @ODEmodel(x1'(t) = x1, x2'(t) = x2, y(t) = x1 + x2(t)) +ode = StructuralIdentifiability.@ODEmodel(x1'(t) = x1, x2'(t) = x2, y(t) = x1 + x2(t)) ident_funcs = [x1 + x2] push!(test_cases, (ode = ode, ident_funcs = ident_funcs, with_states = true))