add SparseVector representation for monomials

Project.toml

@@ -5,9 +5,12 @@ version = "1.0.0-DEV"
 
 [deps]
 AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
+Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 GAP = "c863536a-3901-11e9-33e7-d5cd0df7b904"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
 julia = "1.8"

src/DecomposingPolynomialSystems.jl

@@ -8,12 +8,15 @@ export @var, Variable, Expression, System
 # Exports for testing
 # export degree
 
+using SparseArrays: SparseVector, SparseMatrixCSC, spzeros, AbstractSparseVector, findnz
+using Combinatorics: partitions, multiset_permutations, combinations
 using LinearAlgebra: nullspace
+using UnPack: @unpack
 
 include("utils.jl")
-include("monomials.jl")
-include("expression_map.jl")
 include("sampled_system.jl")
+include("monomials.jl")
+# include("expression_map.jl")
 include("scalings.jl")
 # include("interpolation.jl")
 # include("deck_transformations.jl")

src/deck_transformations.jl

@@ -214,7 +214,6 @@ function symmetries_fixing_parameters_graded!(
     tols::Tolerances=Tolerances(),
     logging::Bool=false
 )
-
     n_unknowns, n_sols, _ = size(F.samples.solutions)  # TODO: what if n_sols is huge?
     C = F.deck_permutations
     symmetries = _init_symmetries(length(C), unknowns(F))
@@ -236,9 +235,9 @@ function symmetries_fixing_parameters_graded!(
                     g = gcd(vcat(num_mons, denom_mons))
                     if isone(g) && !only_param_dep(vcat(num_mons, denom_mons), Vector(1:n_unknowns))
                         if isnothing(eval_num_mons)
-                            eval_num_mons = HC.evaluate(num_mons, F.samples)
+                            eval_num_mons = evaluate(num_mons, F.samples)
                         end
-                        eval_denom_mons = HC.evaluate(denom_mons, F.samples)
+                        eval_denom_mons = evaluate(denom_mons, F.samples)
                         for (j, symmetry) in enumerate(symmetries)
                             if ismissing(symmetry[i])
                                 symmetry[i] = _interpolate_deck_function(
@@ -300,7 +299,10 @@ function symmetries_fixing_parameters_dense!(
         sample_system!(F, n_instances)
 
         logging && println("Evaluating monomials...\n")
-        evaluated_mons = HC.evaluate(mons, F.samples)
+        # TODO: pick samples at which to evaluate monomials
+        # TODO: do it according to n_mons (influences n_samples)
+        # and n_param_dep_only_mons (influences n_instances)
+        evaluated_mons = evaluate(mons, F.samples)
 
         for (i, symmetry) in enumerate(symmetries)
             logging && printstyled("Interpolating the ", i, "-th symmetry map...\n"; color=:blue)

src/interpolation.jl

@@ -4,20 +4,14 @@ export rational_function,
 
 using LinearAlgebra: norm, dot
 
-function check_func_type(func_type::String)
-    if (func_type != "polynomial" && func_type != "rational")
-        error("func_type argument must be either \"polynomial\" or \"rational\"")
-    end
-end
-
 function rational_function(
     coeffs::AbstractVector{<:Number},
     num_mons::MonomialVector,
     denom_mons::MonomialVector;
     logging::Bool=false,
     tol::Float64=1e-5
 )
-    n_num_mons, n_denom_mons = length(num_mons.mds), length(denom_mons.mds)
+    n_num_mons, n_denom_mons = length(num_mons), length(denom_mons)
     @assert length(coeffs) == n_num_mons + n_denom_mons
     num, denom = coeffs[1:n_num_mons], coeffs[n_num_mons+1:end]
     if norm(num) < tol
@@ -53,6 +47,12 @@ function polynomial_function(
     return p
 end
 
+function check_func_type(func_type::String)
+    if (func_type != "polynomial" && func_type != "rational")
+        error("func_type argument must be either \"polynomial\" or \"rational\"")
+    end
+end
+
 function reconstruct_function(
     coeffs::AbstractVector{<:Number},
     mons::MonomialVector;
@@ -160,7 +160,7 @@ end
 # supposes each md in mds is a multidegree in both unknowns and parameters
 # TODO: The implementation below (with _) is more efficient (approx 2x),
 # TODO: since it exploits the sparsity of multidegrees. REMOVE THIS METHOD?
-# function HomotopyContinuation.evaluate(mons::MonomialVector, samples::VarietySamples)
+# function HomotopyContinuation.evaluate(mons::MonomialVector, samples::Samples)
 #     solutions = samples.solutions
 #     parameters = samples.parameters
 
@@ -180,38 +180,6 @@ end
 #     return evaluated_mons
 # end
 
-# TODO: consider view for slices
-function HC.evaluate(
-    mons::MonomialVector,
-    samples::VarietySamples;
-    sparse::Bool=false
-)
-    solutions = samples.solutions
-    parameters = samples.parameters
-
-    n_unknowns, n_sols, n_instances = size(solutions)
-    mds = mons.mds
-    n_mds = length(mds)
-
-    nonzero_ids_unknowns = [findall(!iszero, md[1:n_unknowns]) for md in mds]
-    nonzero_ids_params = [findall(!iszero, md[n_unknowns+1:end]) for md in mds]
-
-    evaluated_mons = zeros(CC, n_mds, n_sols, n_instances)
-    for i in 1:n_instances
-        params = view(parameters, :, i)
-        sols = view(solutions, :, :, i)
-        for (j, md) in enumerate(mds)
-            params_part = params[nonzero_ids_params[j]]
-            md_params_part = md[n_unknowns+1:end][nonzero_ids_params[j]]
-            params_eval = prod(params_part.^md_params_part)
-            sols_part = view(sols, nonzero_ids_unknowns[j], :)
-            md_sols_part = md[1:n_unknowns][nonzero_ids_unknowns[j]]
-            evaluated_mons[j, :, i] = vec(prod(sols_part.^md_sols_part, dims=1)).*params_eval
-        end
-    end
-    return evaluated_mons
-end
-
 function interpolate(
     A::AbstractMatrix{CC},
     mons::MonomialVector;
@@ -254,7 +222,7 @@ function interpolate(
 end
 
 function interpolate_vanishing_polynomials(
-    samples::VarietySamples,
+    samples::Samples,
     vars::Vector{Variable},
     mons::MonomialVector;
     tol::Float64=1e-5