Merge pull request #222 from SciML/moving_saturation

Abstract the saturation variable index

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 -- 

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 -- 

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 -- 

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 -- 

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]

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)

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

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