diff --git a/benchmarking/IdentifiableFunctions/results/Transfection_4State/logs_with_states b/benchmarking/IdentifiableFunctions/results/Transfection_4State/logs_with_states new file mode 100644 index 000000000..e69de29bb diff --git a/benchmarking/IdentifiableFunctions/results/Treatment_io/logs_with_states b/benchmarking/IdentifiableFunctions/results/Treatment_io/logs_with_states new file mode 100644 index 000000000..e69de29bb diff --git a/benchmarking/IdentifiableFunctions/results/cLV1 (2o)/logs_with_states b/benchmarking/IdentifiableFunctions/results/cLV1 (2o)/logs_with_states new file mode 100644 index 000000000..e69de29bb diff --git a/benchmarking/IdentifiableFunctions/results/cholera/logs_with_states b/benchmarking/IdentifiableFunctions/results/cholera/logs_with_states new file mode 100644 index 000000000..b428fd32c --- /dev/null +++ b/benchmarking/IdentifiableFunctions/results/cholera/logs_with_states @@ -0,0 +1,13 @@ +┌ Info: +└ FUNCTION_NAME = "find_identifiable_functions" +┌ Info: +└ PROBLEM_NAME = "cholera" +┌ Info: +└ KWARGS = (with_states = true,) +┌ Info: +└ GLOBAL_ID = :with_states +Internal error: encountered unexpected error in runtime: +ReadOnlyMemoryError() + +Please submit a bug report with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks. +Exception: UNKNOWN at 0x7ffedb2acf19 -- \ No newline at end of file diff --git a/benchmarking/IdentifiableFunctions/results/generalizedLoktaVolterra (1o)/logs_with_states b/benchmarking/IdentifiableFunctions/results/generalizedLoktaVolterra (1o)/logs_with_states new file mode 100644 index 000000000..b40eace03 --- /dev/null +++ b/benchmarking/IdentifiableFunctions/results/generalizedLoktaVolterra (1o)/logs_with_states @@ -0,0 +1,13 @@ +┌ Info: +└ FUNCTION_NAME = "find_identifiable_functions" +┌ Info: +└ PROBLEM_NAME = "generalizedLoktaVolterra (1o)" +┌ Info: +└ KWARGS = (with_states = true,) +┌ Info: +└ GLOBAL_ID = :with_states +Internal error: encountered unexpected error in runtime: +ReadOnlyMemoryError() + +Please submit a bug report with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks. +Exception: UNKNOWN at 0x7ffedb2acf19 -- \ No newline at end of file diff --git a/benchmarking/IdentifiableFunctions/results/generalizedLoktaVolterra (2o)/logs_with_states b/benchmarking/IdentifiableFunctions/results/generalizedLoktaVolterra (2o)/logs_with_states new file mode 100644 index 000000000..775967bc5 --- /dev/null +++ b/benchmarking/IdentifiableFunctions/results/generalizedLoktaVolterra (2o)/logs_with_states @@ -0,0 +1,13 @@ +┌ Info: +└ FUNCTION_NAME = "find_identifiable_functions" +┌ Info: +└ PROBLEM_NAME = "generalizedLoktaVolterra (2o)" +┌ Info: +└ KWARGS = (with_states = true,) +┌ Info: +└ GLOBAL_ID = :with_states +Internal error: encountered unexpected error in runtime: +ReadOnlyMemoryError() + +Please submit a bug report with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks. +Exception: UNKNOWN at 0x7ffedb2acf19 -- \ No newline at end of file diff --git a/benchmarking/IdentifiableFunctions/results/p53/logs_with_states b/benchmarking/IdentifiableFunctions/results/p53/logs_with_states new file mode 100644 index 000000000..e8ecfc6da --- /dev/null +++ b/benchmarking/IdentifiableFunctions/results/p53/logs_with_states @@ -0,0 +1,13 @@ +┌ Info: +└ FUNCTION_NAME = "find_identifiable_functions" +┌ Info: +└ PROBLEM_NAME = "p53" +┌ Info: +└ KWARGS = (with_states = true,) +┌ Info: +└ GLOBAL_ID = :with_states +Internal error: encountered unexpected error in runtime: +ReadOnlyMemoryError() + +Please submit a bug report with steps to reproduce this fault, and any error messages that follow (in their entirety). Thanks. +Exception: UNKNOWN at 0x7ffedb2acf19 -- \ No newline at end of file diff --git a/src/RationalFunctionFields/IdealMQS.jl b/src/RationalFunctionFields/IdealMQS.jl index ba936eb77..bd4d7aabd 100644 --- a/src/RationalFunctionFields/IdealMQS.jl +++ b/src/RationalFunctionFields/IdealMQS.jl @@ -55,20 +55,22 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal function IdealMQS( funcs_den_nums::Vector{Vector{PolyQQ}}; sat_varname = "t", + sat_var_position = :first, ordering = :degrevlex, ) where {PolyQQ} # We are given polynomials of form # [[f1, f2, f3, ...], [g1, g2, g3, ...], ...] # We prepare and store them to construct ideal specializations later @assert !isempty(funcs_den_nums) + @assert sat_var_position in (:first, :last) ordering !== :degrevlex && (@warn "Ordering is not degrevlex but $ordering") ring = parent(first(first(funcs_den_nums))) @debug "Constructing the MQS ideal in $ring" K, n = base_ring(ring), nvars(ring) @debug "Finding pivot polynomials" # In the component f1,f2,... find the polynomial with the minimal total - # degree. Such element will serve as a normalizing term for the - # component + # degree and length. Such element will serve as a normalizing term for + # the component funcs_den_nums = map(plist -> filter(!iszero, plist), funcs_den_nums) # funcs_den_nums = filter(plist -> length(plist) > 1, funcs_den_nums) @assert !isempty(funcs_den_nums) "All elements of the ideal are zero" @@ -90,16 +92,25 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal existing_varnames = map(String, symbols(ring)) ystrs = ["y$i" for i in 1:length(existing_varnames)] @assert !(sat_varname in ystrs) "The name of the saturation variable collided with a primary variable" - sat_var_index = length(ystrs) + 1 - varnames = push!(ystrs, sat_varname) + sat_var_index = if sat_var_position === :first + 1 + else + @assert sat_var_position === :last + length(ystrs) + 1 + end + varnames = append_at_index(ystrs, sat_var_index, sat_varname) @debug "Saturating variable is $sat_varname, index is $sat_var_index" flush(stdout) R_sat, v_sat = Nemo.PolynomialRing(K, varnames, ordering = ordering) # Saturation - @assert sat_var_index == length(v_sat) t_sat = v_sat[sat_var_index] den_lcm_orig = den_lcm - den_lcm = parent_ring_change(den_lcm, R_sat, matching = :byindex) + den_lcm = parent_ring_change( + den_lcm, + R_sat, + matching = :byindex, + shift = Int(sat_var_index == 1), + ) den_lcm_sat = parent_ring_change(den_lcm, R_sat) sat_qq = den_lcm_sat * t_sat - 1 # We construct the array of numerators nums_qq and the array of @@ -114,13 +125,23 @@ mutable struct IdealMQS{T} <: AbstractBlackboxIdeal for i in 1:length(funcs_den_nums) plist = funcs_den_nums[i] den = plist[pivots_indices[i]] - den = parent_ring_change(den, R_sat, matching = :byindex) + den = parent_ring_change( + den, + R_sat, + matching = :byindex, + shift = Int(sat_var_index == 1), + ) push!(dens_qq, den) push!(dens_indices, (length(nums_qq) + 1, length(nums_qq) + length(plist) - 1)) for j in 1:length(plist) j == pivots_indices[i] && continue num = plist[j] - num = parent_ring_change(num, R_sat, matching = :byindex) + num = parent_ring_change( + num, + R_sat, + matching = :byindex, + shift = Int(sat_var_index == 1), + ) push!(nums_qq, num) end end @@ -170,6 +191,7 @@ end end function fractionfree_generators_raw(mqs::IdealMQS) + # TODO: this assumes mqs.sat_var_index is last, and thus is broken ring_params = ParamPunPam.parent_params(mqs) K = base_ring(ring_params) varnames = map(string, Nemo.symbols(ring_params)) @@ -243,15 +265,13 @@ function ParamPunPam.specialize_mod_p( @assert length(point) == nvars(ParamPunPam.parent_params(mqs)) # +1 actual variable because of the saturation! @assert length(point) + 1 == nvars(parent(nums_gf[1])) - # NOTE: Assuming the saturating variable is the last one - @assert mqs.sat_var_index == length(point) + 1 - point_sat = vcat(point, one(K_1)) + point_sat = append_at_index(point, mqs.sat_var_index, one(K_1)) 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) @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("Ideal: $mqs\nPoint: $point") + iszero(den_spec) && __throw_unlucky_evaluation("Point: $point") span = dens_indices[i] for j in span[1]:span[2] num, num_spec = nums_gf[j], nums_gf_spec[j] @@ -273,15 +293,13 @@ function specialize(mqs::IdealMQS, point::Vector{Nemo.fmpq}; saturated = true) @assert length(point) == nvars(ParamPunPam.parent_params(mqs)) # +1 actual variable because of the saturation! @assert length(point) + 1 == nvars(parent(nums_qq[1])) - # NOTE: Assuming the saturating variable is the last one - @assert mqs.sat_var_index == length(point) + 1 - point_sat = vcat(point, one(K)) + point_sat = append_at_index(point, mqs.sat_var_index, one(K)) 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) @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("Ideal: $mqs\nPoint: $point") + iszero(den_spec) && __throw_unlucky_evaluation("Point: $point") span = dens_indices[i] for j in span[1]:span[2] num, num_spec = nums_qq[j], nums_qq_spec[j] diff --git a/src/RationalFunctionFields/RationalFunctionField.jl b/src/RationalFunctionFields/RationalFunctionField.jl index fa17964be..be4f7203a 100644 --- a/src/RationalFunctionFields/RationalFunctionField.jl +++ b/src/RationalFunctionFields/RationalFunctionField.jl @@ -493,7 +493,7 @@ function simplified_generating_set( check_variables = false, # almost always slows down and thus turned off rational_interpolator = :VanDerHoevenLecerf, ) - @info "Simplifying identifiable functions" + @info "Simplifying generating set" _runtime_logger[:id_groebner_time] = 0.0 _runtime_logger[:id_calls_to_gb] = 0 _runtime_logger[:id_inclusion_check_mod_p] = 0.0 @@ -502,8 +502,7 @@ function simplified_generating_set( _runtime_logger[:id_normalforms_time] = 0.0 _runtime_logger[:id_ranking] = 0 - # Checking identifiability of particular variables and adding them to the - # field + # Checking membership of particular variables and adding them to the field if check_variables vars = gens(poly_ring(rff)) containment = field_contains(rff, vars, (1.0 + p) / 2) diff --git a/src/RationalFunctionFields/normalforms.jl b/src/RationalFunctionFields/normalforms.jl index 219c6b0b6..14d7a00df 100644 --- a/src/RationalFunctionFields/normalforms.jl +++ b/src/RationalFunctionFields/normalforms.jl @@ -54,15 +54,14 @@ function local_normal_forms( ) @assert !isempty(point) @assert parent(first(point)) == finite_field - point_ff_ext = vcat(point, one(finite_field)) + point_ff_ext = append_at_index(point, mqs.sat_var_index, one(finite_field)) 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) - @assert mqs.sat_var_index == length(xs_ff) - xs_ff = xs_ff[1:(end - 1)] + xs_ff = cut_at_index(xs_ff, mqs.sat_var_index) pivot_vectors = map(f -> exponent_vector(f, 1), xs_ff) @debug """ variables in the finite field: $(xs_ff) @@ -358,7 +357,7 @@ 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 = vcat(point, one(finite_field)) + point_ext = append_at_index(point, mqs.sat_var_index, one(finite_field)) for i in 1:length(complete_intersection_relations_ff) relation = complete_intersection_relations_ff[i] relation_mqs = relation - evaluate(relation, point_ext) @@ -409,7 +408,10 @@ function linear_relations_between_normal_forms( modulo: $finite_field""" throw(ErrorException("Rational reconstruction failed.")) end - relation_qq_param = evaluate(relation_qq, vcat(xs_param, one(ring))) + relation_qq_param = evaluate( + relation_qq, + append_at_index(xs_param, mqs.sat_var_index, one(ring_param)), + ) relations_qq[i] = relation_qq_param // one(relation_qq_param) end relations_qq diff --git a/src/util.jl b/src/util.jl index 474c73d32..af9d95c2d 100644 --- a/src/util.jl +++ b/src/util.jl @@ -10,6 +10,27 @@ end # ------------------------------------------------------------------------------ +function append_at_index(vec::Vector{T}, idx::Integer, val::T) where {T} + # NOTE: could also just use insert! + @assert (idx == 1) || (idx == length(vec) + 1) + if idx == 1 + vcat(val, vec) + else + vcat(vec, val) + end +end + +function cut_at_index(vec::Vector{T}, idx::Integer) where {T} + @assert (idx == 1) || (idx == length(vec)) + if idx == 1 + vec[2:end] + else + vec[1:(end - 1)] + end +end + +# ------------------------------------------------------------------------------ + """ eval_at_dict(f, d) @@ -209,7 +230,12 @@ Input Output: - a polynomial in `new_ring` “equal” to `poly` """ -function parent_ring_change(poly::MPolyElem, new_ring::MPolyRing; matching = :byname) +function parent_ring_change( + poly::MPolyElem, + new_ring::MPolyRing; + matching = :byname, + shift = 0, +) old_ring = parent(poly) # Construct a mapping for the variable indices. # Zero indicates no image of the old variable in the new ring @@ -223,7 +249,7 @@ function parent_ring_change(poly::MPolyElem, new_ring::MPolyRing; matching = :by var_mapping[i] = found end elseif matching === :byindex - var_mapping[1:nvars(new_ring)] .= 1:nvars(new_ring) + var_mapping[1:(nvars(new_ring) - shift)] .= (1 + shift):nvars(new_ring) else throw(Base.ArgumentError("Unknown matching type: $matching")) end @@ -393,7 +419,7 @@ end function str_to_var(s::String, ring::MPolyRing) ind = findfirst(v -> (string(v) == s), symbols(ring)) - if ind == nothing + if ind === nothing throw(Base.KeyError("Variable $s is not found in ring $ring")) end return gens(ring)[ind] @@ -451,6 +477,7 @@ For a variable `v`, returns a variable in `ring` with the same name """ function switch_ring(v::MPolyElem, ring::MPolyRing) ind = findfirst(vv -> vv == v, gens(parent(v))) + return str_to_var(string(symbols(parent(v))[ind]), ring) end