Add groebner_learn / groebner_apply!

experimental/SI/fujita.jl

@@ -0,0 +1,84 @@
+using StructuralIdentifiability
+using AbstractAlgebra, BenchmarkTools, JLD2
+import Nemo, Profile
+
+macro myprof(ex)
+    :(
+        (VSCodeServer.Profile).clear();
+        Profile.init(n=10^7, delay=0.0000001);
+        Profile.@profile $ex;
+        VSCodeServer.view_profile(;)
+    )
+end
+
+ode = @ODEmodel(
+    EGFR'(t) =
+        EGFR_turnover * pro_EGFR(t) + EGF_EGFR(t) * reaction_1_k2 -
+        EGFR(t) * EGFR_turnover - EGF_EGFR(t) * reaction_1_k1,
+    pEGFR'(t) =
+        EGF_EGFR(t) * reaction_9_k1 - pEGFR(t) * reaction_4_k1 +
+        pEGFR_Akt(t) * reaction_2_k2 +
+        pEGFR_Akt(t) * reaction_3_k1 - Akt(t) * pEGFR(t) * reaction_2_k1,
+    pEGFR_Akt'(t) =
+        Akt(t) * pEGFR(t) * reaction_2_k1 - pEGFR_Akt(t) * reaction_3_k1 -
+        pEGFR_Akt(t) * reaction_2_k2,
+    Akt'(t) =
+        pAkt(t) * reaction_7_k1 + pEGFR_Akt(t) * reaction_2_k2 -
+        Akt(t) * pEGFR(t) * reaction_2_k1,
+    pAkt'(t) =
+        pAkt_S6(t) * reaction_5_k2 - pAkt(t) * reaction_7_k1 +
+        pAkt_S6(t) * reaction_6_k1 +
+        pEGFR_Akt(t) * reaction_3_k1 - S6(t) * pAkt(t) * reaction_5_k1,
+    S6'(t) =
+        pAkt_S6(t) * reaction_5_k2 + pS6(t) * reaction_8_k1 -
+        S6(t) * pAkt(t) * reaction_5_k1,
+    pAkt_S6'(t) =
+        S6(t) * pAkt(t) * reaction_5_k1 - pAkt_S6(t) * reaction_6_k1 -
+        pAkt_S6(t) * reaction_5_k2,
+    pS6'(t) = pAkt_S6(t) * reaction_6_k1 - pS6(t) * reaction_8_k1,
+    EGF_EGFR'(t) =
+        EGF_EGFR(t) * reaction_1_k1 - EGF_EGFR(t) * reaction_9_k1 -
+        EGF_EGFR(t) * reaction_1_k2,
+    y1(t) = a1 * (pEGFR(t) + pEGFR_Akt(t)),
+    y2(t) = a2 * (pAkt(t) + pAkt_S6(t)),
+    y3(t) = a3 * pS6(t)
+)
+
+G = StructuralIdentifiability.ideal_generators(ode);
+
+par = parent(G[1])
+FF = Nemo.GF(2^31 - 1)
+point = Nemo.QQ.(rand(1:100, length(AbstractAlgebra.gens(base_ring(base_ring(par))))))
+G_zz = map(
+    f -> map_coefficients(
+        c -> evaluate(numerator(c), point) // evaluate(denominator(c), point),
+        f
+    ),
+    G
+);
+G_zp = map(
+    f -> map_coefficients(c -> FF(BigInt(numerator(c))) // FF(BigInt(denominator(c))), f),
+    G_zz
+);
+typeof(parent(G_zp[1]))
+typeof(G_zp[1])
+
+Groebner.logging_enabled() = false
+Groebner.invariants_enabled() = false
+
+gb = Groebner.groebner(G_zp);
+graph, gb_1 = Groebner.groebner_learn(G_zp);
+graph
+graph.matrix_infos
+flag, gb_2 = Groebner.groebner_apply!(graph, G_zp)
+@assert flag && (gb == gb_1 == gb_2)
+
+@btime Groebner.groebner($G_zp);
+@btime Groebner.groebner_apply!($graph, $G_zp);
+
+@myprof begin
+    for i in 1:1000
+        Groebner.groebner_apply!(graph, G_zp)
+    end
+end
+

src/Groebner.jl

@@ -20,7 +20,7 @@ invariants_enabled() = true
 Specifies if custom logging is enabled.
 If `false`, then all logging is disabled, and entails no runtime overhead.
 
-See also `@log` in `src/utils/logging.jl`.
+See also `@log` in `src/uti0-w2q3=-asqqw2q=[-edrfls/logging.jl`.
 """
 logging_enabled() = true
 
@@ -123,7 +123,7 @@ include("interface.jl")
 using SnoopPrecompile
 include("precompile.jl")
 
-export groebner, groebner_learn, groebner_apply
+export groebner, groebner_learn, groebner_apply!
 export isgroebner
 export normalform
 

src/f4/f4.jl

@@ -66,6 +66,52 @@ function initialize_structs(
     basis, pairset, hashtable
 end
 
+function initialize_structs_learn(
+    ring::PolyRing,
+    monoms::Vector{Vector{M}},
+    coeffs::Vector{Vector{C}},
+    params::AlgorithmParameters;
+    normalize=true
+) where {M <: Monom, C <: Coeff}
+    @log level = 3 "Initializing structs.."
+
+    tablesize = select_hashtable_size(ring, monoms)
+    @log level = 3 "Initial hashtable size is $tablesize"
+
+    # Basis for storing basis elements,
+    # Pairset for storing critical pairs of basis elements,
+    # Hashtable for hashing monomials stored in the basis
+    basis = initialize_basis(ring, length(monoms), C)
+    pairset = initialize_pairset(entrytype(M))
+    hashtable = initialize_hashtable(ring, params.rng, M, tablesize)
+
+    # Filling the basis and hashtable with the given inputs
+    fill_data!(basis, hashtable, monoms, coeffs)
+    fill_divmask!(hashtable)
+
+    @log level = 4 "Hashtable:"
+
+    @log level = 4 "Sorting input polynomials by their leading terms in non-decreasing order"
+    permutation = sort_polys_by_lead_increasing!(basis, hashtable)
+
+    # Divide each polynomial by the leading coefficient
+    if normalize
+        @log level = 4 "Normalizing input polynomials"
+        normalize_basis!(ring, basis)
+    end
+
+    @log level = 5 "Initializing computation graph"
+    graph = initialize_computation_graph_f4(
+        ring,
+        deepcopy_basis(basis),
+        basis,
+        hashtable,
+        permutation
+    )
+
+    graph, basis, pairset, hashtable
+end
+
 # Initializes Basis and MonomialHashtable structures,
 # fills input data from exponents and coeffs
 #
@@ -228,20 +274,23 @@ end
 
 # Given a `basis` object that stores some groebner basis
 # performs basis interreduction and writes the result to `basis` inplace
+# TODO: f4_reducegb!
 function reducegb_f4!(
     ring::PolyRing,
     basis::Basis,
     matrix::MacaulayMatrix,
     ht::MonomialHashtable{M},
     symbol_ht::MonomialHashtable{M}
 ) where {M}
+    @log level = 100000 "Entering autoreduction" basis
+
     etmp = construct_const_monom(M, ht.nvars)
     # etmp is now set to zero, and has zero hash
 
     reinitialize_matrix!(matrix, basis.ndivmasks)
     uprows = matrix.uprows
 
-    # add all non redundant elements from basis
+    # add all non redundant elements from the basis
     # as matrix upper rows
     @inbounds for i in 1:(basis.ndivmasks) #
         matrix.nrows += 1
@@ -646,6 +695,29 @@ function basis_well_formed(key, ring, basis, hashtable)
     if key in (:input_f4!, :input_f4_learn!, :input_f4_apply!)
         (isempty(basis.monoms) || isempty(basis.coeffs)) && return false
         (basis.size == 0 || basis.ntotal == 0) && return false
+    elseif key in (:output_f4!, :output_f4_learn!, :output_f4_apply!)
+        # TODO: also check: 
+        #   are sorted: sort_polys_by_lead_increasing!(basis, ht, ord=ord)
+        basis.ndivmasks ==
+        length(basis.coeffs) ==
+        length(basis.monoms) ==
+        length(basis.divmasks) ==
+        length(basis.nonredundant) ==
+        length(basis.isredundant) || return false
+        basis.nonredundant == collect(1:(basis.ndivmasks)) || return false
+        any(!iszero, basis.isredundant) && return false
+        any(c -> !isone(c[1]), basis.coeffs) && return false
+    else
+        return false
+    end
+    for i in 1:length(basis.coeffs)
+        if !isassigned(basis.coeffs, i)
+            if isassigned(basis.monoms, i)
+                return false
+            end
+        else
+            length(basis.coeffs[i]) != length(basis.monoms[i]) && return false
+        end
     end
     true
 end
@@ -788,7 +860,7 @@ function f4!(
     standardize_basis!(ring, basis, hashtable, hashtable.ord)
 
     # @invariant hashtable_well_formed(:output_f4!, ring, hashtable)
-    # @invariant basis_well_formed(:output_f4!, ring, basis, hashtable)
+    @invariant basis_well_formed(:output_f4!, ring, basis, hashtable)
 
     nothing
 end
@@ -824,7 +896,7 @@ function f4_normalform!(
     ring::PolyRing,
     basis::Basis{C},
     tobereduced::Basis{C},
-    ht::MonomialHashtable,
+    ht::MonomialHashtable
 ) where {C <: Coeff}
     matrix = initialize_matrix(ring, C)
     symbol_ht = initialize_secondary_hashtable(ht)

src/f4/graph.jl

@@ -1,15 +1,23 @@
 
 mutable struct ComputationGraphF4{C, M, Ord}
+    ring::PolyRing{Ord}
     input_basis::Basis{C}
+    buf_basis::Basis{C}
     gb_basis::Basis{C}
     hashtable::MonomialHashtable{M, Ord}
+    input_permutation::Vector{Int}
     # The number of columns,
     # The number of upper rows,
     # The number of lower rows
     matrix_infos::Vector{NamedTuple{(:nup, :nlow, :ncols), Tuple{Int64, Int64, Int64}}}
     matrix_nonzeroed_rows::Vector{Vector{Int}}
     matrix_upper_rows::Vector{Tuple{Vector{Int}, Vector{MonomIdx}}}
     matrix_lower_rows::Vector{Tuple{Vector{Int}, Vector{MonomIdx}}}
+    output_nonredundant_indices::Vector{Int}
+    nonredundant_indices_before_reduce::Vector{Int}
+    output_sort_indices::Vector{Int}
+    # matrix_autored_lower_rows::Vector{Tuple{Vector{Int}, Vector{MonomIdx}}}
+    # matrix_autored_upper_rows::Vector{Tuple{Vector{Int}, Vector{MonomIdx}}}
     # F4 iteration number --> index in the basis
     # matrix_tobereduced_rows::Vector{Vector{Int}}
     # matrix_tobereduced_mult::Vector{Vector{MonomIdx}}
@@ -19,22 +27,50 @@ mutable struct ComputationGraphF4{C, M, Ord}
     # matrix_columns::Vector{Vector{Int}}
 end
 
-function initialize_computation_graph_f4(input_basis, gb_basis, hashtable)
+function initialize_computation_graph_f4(ring, input_basis, gb_basis, hashtable, permutation)
     ComputationGraphF4(
+        ring,
         input_basis,
+        deepcopy_basis(gb_basis),
         gb_basis,
         hashtable,
+        permutation,
         Vector{NamedTuple{(:nup, :nlow, :ncols), Tuple{Int64, Int64, Int64}}}(),
         Vector{Vector{Int}}(),
         Vector{Tuple{Vector{Int}, Vector{MonomIdx}}}(),
         Vector{Tuple{Vector{Int}, Vector{MonomIdx}}}(),
+        Vector{Int}(),
+        Vector{Int}(),
+        Vector{Int}(),
+        # Vector{Tuple{Vector{Int}, Vector{MonomIdx}}}(),
+        # Vector{Tuple{Vector{Int}, Vector{MonomIdx}}}(),
     )
 end
 
+function finalize_graph!(graph::ComputationGraphF4) 
+    # TODO: trim sizes
+    graph.buf_basis = deepcopy_basis(graph.gb_basis)
+    graph.buf_basis.ndivmasks = graph.input_basis.ndivmasks
+    graph.buf_basis.nprocessed = graph.input_basis.nprocessed
+    graph.buf_basis.ntotal = graph.input_basis.ntotal
+    # TODO: set size to the allocated size
+    # done.
+    # graph.buf_basis.size = graph.input_basis.size
+    nothing
+end
+
 function Base.show(io::IO, ::MIME"text/plain", graph::ComputationGraphF4)
-    println(io, "Computation graph of F4.")
-    total_matrix_rows = sum(x -> x.nlow, graph.matrix_infos)
-    # total_redundant_rows = sum(length, graph.matrix_zeroed_rows)
-    println(io, "Matrix rows in total: $(total_matrix_rows)")
-    # print("Nonzero matrix rows: $(total_redundant_rows) ($(round(total_redundant_rows / total_matrix_rows * 100, digits=2)) %)")
+    sz = Base.summarysize(graph) >> 10
+    println(io, "Computation graph of F4 ($sz KiB)")
+    println(io, "Recorded in TODO with parameters: TODO")
+    total_iterations = length(graph.matrix_infos)
+    total_matrix_low_rows = sum(x -> x.nlow, graph.matrix_infos)
+    total_matrix_up_rows = sum(x -> x.nup, graph.matrix_infos)
+    # TODO
+    total_matrix_up_rows_useful = sum(length ∘ first, graph.matrix_upper_rows)
+    total_matrix_low_rows_useful = sum(length ∘ first, graph.matrix_lower_rows)
+    println(io, """
+    Total iterations: $(total_iterations)
+    Total upper matrix rows: $(total_matrix_up_rows) ($(round(total_matrix_up_rows_useful / total_matrix_up_rows * 100, digits=2)) % are useful)
+    Total lower matrix rows: $(total_matrix_low_rows) ($(round(total_matrix_low_rows_useful / total_matrix_low_rows * 100, digits=2)) % are useful)""")
 end