rewrite SampledSystem struct

MultivariatePolynomialSystems · Nov 17, 2023 · 4cb7428 · 4cb7428
1 parent de4fe6f
commit 4cb7428
Show file tree

Hide file tree

Showing 7 changed files with 253 additions and 206 deletions.
diff --git a/.gitignore b/.gitignore
@@ -1,4 +1,5 @@
 .DS_Store
 /Manifest.toml
 /docs/Manifest.toml
-/docs/build/
+/docs/build/
+.vscode
diff --git a/README.md b/README.md
@@ -3,7 +3,7 @@
 [![License](https://img.shields.io/badge/License-MIT-yellow?style=flat&label=License&color=%23FFA500%09)](https://github.com/MultivariatePolynomialSystems/DecomposingPolynomialSystems.jl/blob/main/LICENSE)
 [![Dev](https://img.shields.io/badge/docs-dev-blue.svg)](https://multivariatepolynomialsystems.github.io/DecomposingPolynomialSystems.jl/dev)
 
-DecomposingPolynomialSystems.jl is a Julia package that computes the symmetries that fix the parameters (specifically, the group of deck transformations) of a parametric polynomial system with finitely many solutions for generic parameters with a view towards decomposing the given polynomial system.
+DecomposingPolynomialSystems.jl is a Julia package that computes the symmetries that fix the parameters (specifically, the group of deck transformations) of a parametric polynomial system with finitely many solutions with a view towards decomposing the given polynomial system.
 
 ## Installation
 

diff --git a/src/DecomposingPolynomialSystems.jl b/src/DecomposingPolynomialSystems.jl
@@ -15,8 +15,8 @@ include("monomials.jl")
 include("expression_map.jl")
 include("sampled_system.jl")
 include("scalings.jl")
-include("interpolation.jl")
-include("deck_transformations.jl")
+# include("interpolation.jl")
+# include("deck_transformations.jl")
 # include("invariants.jl")
 # include("implicitization.jl")
 # include("decompose.jl")

diff --git a/src/deck_transformations.jl b/src/deck_transformations.jl
@@ -96,16 +96,17 @@ end
 function _deck_vandermonde_matrix(
     deck_permutation::Vector{Int},
     function_id::Int,
-    solutions::AbstractArray{T, 3},
-    eval_num_mons::AbstractArray{T, 3},
-    eval_denom_mons::AbstractArray{T, 3}
-) where {T<:Complex}
-
+    solutions::AbstractArray{ComplexF64, 3},
+    eval_num_mons::AbstractArray{ComplexF64, 3},
+    eval_denom_mons::AbstractArray{ComplexF64, 3}
+)
     _, n_sols, n_instances = size(solutions)
     n_num_mons = size(eval_num_mons, 1)
     n_denom_mons = size(eval_denom_mons, 1)
 
-    A = zeros(T, n_instances*n_sols, n_num_mons+n_denom_mons)
+    # TODO: make n_rows close to n_cols by randomly picking solution ids
+    # TODO: sol_ids = rand(1:nsols, Int(ceil((n_num_mons+n_denom_mons)/n_instances)))
+    A = zeros(ComplexF64, n_instances*n_sols, n_num_mons+n_denom_mons)
     @assert size(A, 1) >= size(A, 2)
 
     for i in 1:n_sols
@@ -208,91 +209,76 @@ end
 
 function symmetries_fixing_parameters_graded!(
     F::SampledSystem,
-    scalings::ScalingGroup,
-    mon_classes::Dict{Vector{Int}, MonomialVector{T}};
+    scalings::ScalingGroup;
+    degree_bound::Integer=1,
     tols::Tolerances=Tolerances(),
     logging::Bool=false
-) where {T<:Integer}
+)
 
-    max_n_mons = max(length.(collect(values(mon_classes)))...)  # size of the largest class
     n_unknowns, n_sols, _ = size(F.samples.solutions)  # TODO: what if n_sols is huge?
-    n_instances = Int(ceil(2/n_sols*max_n_mons))
-
     C = F.deck_permutations
     symmetries = _init_symmetries(length(C), unknowns(F))
 
-    sample_system!(F, n_instances)
-
-    for (num_deg, num_mons) in mon_classes
-        eval_num_mons = nothing
-        for i in 1:n_unknowns
-            denom_deg = _denom_deg(num_deg, scalings.grading, i)  # i-th variable
-            denom_mons = get(mon_classes, denom_deg, nothing)
-            if !isnothing(denom_mons)
-                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)
-                    end
-                    eval_denom_mons = HC.evaluate(denom_mons, F.samples)
-                    for (j, symmetry) in enumerate(symmetries)
-                        if ismissing(symmetry[i])
-                            symmetry[i] = _interpolate_deck_function(
-                                C[j],
-                                i,
-                                F.samples.solutions,
-                                eval_num_mons,
-                                eval_denom_mons,
-                                num_mons,
-                                denom_mons,
-                                tols;
-                                logging=logging
-                            )
-                            if logging && !ismissing(symmetry[i])
-                                printstyled(
-                                    "Good representative for the ",
-                                    to_ordinal(j),
-                                    " symmetry, variable ",
-                                    unknowns(F)[i],
-                                    ":\n",
-                                    color=:red
+    for d in 1:degree_bound
+        mons = MonomialVector{Int8}(scalings.vars; degree=d)
+        mon_classes = to_classes(mons, scalings.grading)
+
+        max_n_mons = max(length.(collect(values(mon_classes)))...)  # size of the largest class
+        n_instances = Int(ceil(2/n_sols*max_n_mons))
+        sample_system!(F, n_instances)
+
+        for (num_deg, num_mons) in mon_classes
+            eval_num_mons = nothing
+            for i in 1:n_unknowns
+                denom_deg = _denom_deg(num_deg, scalings.grading, i)  # i-th variable
+                denom_mons = get(mon_classes, denom_deg, nothing)
+                if !isnothing(denom_mons)
+                    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)
+                        end
+                        eval_denom_mons = HC.evaluate(denom_mons, F.samples)
+                        for (j, symmetry) in enumerate(symmetries)
+                            if ismissing(symmetry[i])
+                                symmetry[i] = _interpolate_deck_function(
+                                    C[j],
+                                    i,
+                                    F.samples.solutions,
+                                    eval_num_mons,
+                                    eval_denom_mons,
+                                    num_mons,
+                                    denom_mons,
+                                    tols;
+                                    logging=logging
                                 )
-                                println(symmetry[i])
+                                if logging && !ismissing(symmetry[i])
+                                    printstyled(
+                                        "Good representative for the ",
+                                        to_ordinal(j),
+                                        " symmetry, variable ",
+                                        unknowns(F)[i],
+                                        ":\n",
+                                        color=:red
+                                    )
+                                    println(symmetry[i])
+                                end
                             end
                         end
                     end
                 end
             end
-        end
 
-        if _all_interpolated(symmetries)
-            logging && printstyled("--- All symmetries are interpolated ---\n", color=:blue)
-            return DeckTransformationGroup(symmetries, F)
+            if _all_interpolated(symmetries)
+                logging && printstyled("--- All symmetries are interpolated ---\n", color=:blue)
+                return DeckTransformationGroup(symmetries, F)
+            end
         end
     end
 
     return DeckTransformationGroup(symmetries, F)
 end
 
-function symmetries_fixing_parameters_graded!(
-    F::SampledSystem,
-    scalings::ScalingGroup;
-    degree_bound::Integer=1,
-    tols::Tolerances=Tolerances(),
-    logging::Bool=false
-)
-
-    mons = MonomialVector{Int8}(scalings.vars; degree=degree_bound)
-    mon_classes = to_classes(mons, scalings.grading)
-    return symmetries_fixing_parameters_graded!(
-        F,
-        scalings,
-        mon_classes;
-        tols=tols,
-        logging=logging
-    )
-end
-
 function symmetries_fixing_parameters_dense!(
     F::SampledSystem; 
     degree_bound::Integer=1,

diff --git a/src/interpolation.jl b/src/interpolation.jl
@@ -20,17 +20,18 @@ function rational_function(
     n_num_mons, n_denom_mons = length(num_mons.mds), length(denom_mons.mds)
     @assert length(coeffs) == n_num_mons + n_denom_mons
     num, denom = coeffs[1:n_num_mons], coeffs[n_num_mons+1:end]
-    if sum((!iszero).(denom)) == 1
-        id = findfirst(!iszero, denom)
-        num /= denom[id]
-        sparsify!(num, tol; digits=1)
-        denom[id] = 1
-    end
-    num, denom = simplify_numbers(num), simplify_numbers(denom)
     if norm(num) < tol
         @warn "Numerator close to zero"
     end
     @assert norm(denom) > tol
+    nonzero_ids = findall(!iszero, denom)
+    nonzero_denom = denom[nonzero_ids]
+    if all(x->x==nonzero_denom[1], nonzero_denom)
+        num /= nonzero_denom[1]
+        sparsify!(num, tol; digits=1)
+        denom[nonzero_ids] = ones(length(nonzero_ids))
+    end
+    num, denom = simplify_numbers(num), simplify_numbers(denom)
     if logging
         println("numerator = ", sum(to_expressions(num_mons).*num))
         println("denominator = ", sum(to_expressions(denom_mons).*denom))