From 7475ece7285e5b5d9864e572101b460db642a5cf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 28 Nov 2023 16:55:37 +0100 Subject: [PATCH 1/3] Add support for orthonormal basis --- src/dependence.jl | 38 ++++++++++++++++++++++++++------------ test/null.jl | 17 ++++++++++------- 2 files changed, 36 insertions(+), 19 deletions(-) diff --git a/src/dependence.jl b/src/dependence.jl index 77f1941..001fb8e 100644 --- a/src/dependence.jl +++ b/src/dependence.jl @@ -219,6 +219,21 @@ function Base.convert( return BasisDependence{StaircaseDependence}(d, d.basis) end +_full_basis(basis::MB.AbstractPolynomialBasis) = basis, eachindex(basis) + +function _full_basis(basis::MB.AbstractMonomialBasis) + full = MB.maxdegree_basis( + typeof(basis), + MP.variables(basis), + MP.maxdegree(basis.monomials), + ) + index_map = Int[ + something(_index(basis, full[i]), 0) + for i in eachindex(full) + ] + return full, index_map +end + """ function BasisDependence{StaircaseDependence}( is_dependent::Function, @@ -236,29 +251,28 @@ Foundations of Computational Mathematics 8 (2008): 607-647. """ function BasisDependence{StaircaseDependence}( r, - basis::MB.MonomialBasis{M}, -) where {M} - if isempty(basis.monomials) + basis::MB.AbstractPolynomialBasis, +) + if isempty(basis) return BasisDependence(basis, StaircaseDependence[]) end - function dependence(mono) - i = _index(basis, mono) + vars = MP.variables(basis) + full_basis, index_map = _full_basis(basis) + function dependence(full_index) + i = index_map[full_index] return if isnothing(i) TRIVIAL else is_dependent!(r, i) ? DEPENDENT : INDEPENDENT end end - vars = MP.variables(basis) - full_basis = - MB.maxdegree_basis(typeof(basis), vars, MP.maxdegree(basis.monomials)) d = StaircaseDependence[] # This sieve of [LLR08, Algorithm 1] is a performance improvement but not only. # It also ensures that the standard monomials have the "staircase structure". function is_corner_multiple(mono, indices, dependence) for i in eachindex(dependence) if is_dependent(dependence[i]) && - MP.divides(full_basis.monomials[indices[i]], mono) + MP.divides(full_basis.monomials[indices[i]], mono) # TODO tol return true end end @@ -266,10 +280,10 @@ function BasisDependence{StaircaseDependence}( end keep = Int[] # Compute standard monomials and corners - for (i, mono) in enumerate(full_basis.monomials) - if !is_corner_multiple(mono, keep, d) + for i in eachindex(full_basis) + if !is_corner_multiple(full_basis[i], keep, d) push!(keep, i) - push!(d, StaircaseDependence(true, dependence(mono))) + push!(d, StaircaseDependence(true, dependence(i))) end end full_basis = typeof(full_basis)(full_basis.monomials[keep]) diff --git a/test/null.jl b/test/null.jl index 41318f6..77b2223 100644 --- a/test/null.jl +++ b/test/null.jl @@ -18,13 +18,16 @@ function test_partial_commutative_fix(x, y) 0 0 1 0 # x * y 1 0 0 0 # x^2 ] - basis = MB.MonomialBasis(MP.monomials((x, y), 0:2)) - null = MM.MacaulayNullspace(matrix, basis, 1e-8) - D = MM.StaircaseDependence - solver = MM.StaircaseSolver{Float64}(max_partial_iterations = 1) - shift_solver = MM.ShiftNullspace{D}(solver) - sol = MM.solve(null, shift_solver) - return testelements(sol, [[1, 1], [-1, 1]], 1e-8) + monos = MP.monomials((x, y), 0:2) + for basis in [MB.MonomialBasis(monos), MB.OrthonormalCoefficientsBasis(monos)] + null = MM.MacaulayNullspace(matrix, basis, 1e-8) + D = MM.StaircaseDependence + solver = MM.StaircaseSolver{Float64}(max_partial_iterations = 1) + shift_solver = MM.ShiftNullspace{D}(solver) + sol = MM.solve(null, shift_solver) + testelements(sol, [[1, 1], [-1, 1]], 1e-8) + end + return end function test_dependent_border(x, y) From fc9bc3c143cd7b6aaf536630e250edfa9149d897 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 30 Nov 2023 11:53:27 +0100 Subject: [PATCH 2/3] Add support for orthogonal basis --- src/dependence.jl | 90 ++++++++++++++++++++++++++++++----------------- test/null.jl | 2 +- 2 files changed, 58 insertions(+), 34 deletions(-) diff --git a/src/dependence.jl b/src/dependence.jl index 001fb8e..db4ea56 100644 --- a/src/dependence.jl +++ b/src/dependence.jl @@ -260,58 +260,82 @@ function BasisDependence{StaircaseDependence}( full_basis, index_map = _full_basis(basis) function dependence(full_index) i = index_map[full_index] - return if isnothing(i) + return if iszero(i) TRIVIAL else is_dependent!(r, i) ? DEPENDENT : INDEPENDENT end end - d = StaircaseDependence[] # This sieve of [LLR08, Algorithm 1] is a performance improvement but not only. # It also ensures that the standard monomials have the "staircase structure". - function is_corner_multiple(mono, indices, dependence) - for i in eachindex(dependence) - if is_dependent(dependence[i]) && - MP.divides(full_basis.monomials[indices[i]], mono) # TODO tol + corners = Int[] + function is_corner_multiple(mono) + for corner in corners + if MP.divides(full_basis[corner], mono) # TODO tol return true end end return false end - keep = Int[] # Compute standard monomials and corners - for i in eachindex(full_basis) - if !is_corner_multiple(full_basis[i], keep, d) - push!(keep, i) - push!(d, StaircaseDependence(true, dependence(i))) + keep = falses(length(full_basis)) + d = Vector{StaircaseDependence}(undef, length(full_basis)) + for (i, mono) in enumerate(full_basis) + if !is_corner_multiple(mono) + keep[i] = true + dep = dependence(i) + d[i] = StaircaseDependence(true, dep) + if is_dependent(dep) + push!(corners, i) + end end end - full_basis = typeof(full_basis)(full_basis.monomials[keep]) - new_basis = MB.MonomialBasis( - eltype(basis.monomials)[ - full_basis.monomials[i] * shift for - i in eachindex(d) if !is_dependent(d[i]) for shift in vars - ], - ) - full_basis, I1, I2 = MB.merge_bases(full_basis, new_basis) - deps = Vector{StaircaseDependence}(undef, length(full_basis.monomials)) - for (i, mono) in enumerate(full_basis.monomials) - if iszero(I1[i]) - @assert !iszero(I2[i]) - if is_corner_multiple(mono, 1:(i-1), view(deps, 1:(i-1))) - std = false - else - # If it was not seen before, it means it is outside the basis - # so it is trivial standard - @assert isnothing(_index(basis, mono)) - std = true + new_basis = eltype(MB.generators(full_basis))[] + new_deps = StaircaseDependence[] + for (i, mono) in enumerate(full_basis) + if keep[i] && is_standard(d[i]) + for shift in vars + shifted = mono * shift + j = _index(full_basis, shifted) + if isnothing(j) + push!(new_basis, shifted) + push!(new_deps, StaircaseDependence( + is_corner_multiple(shifted), + TRIVIAL, + )) + else + if keep[j] + continue + end + keep[j] = true + d[j] = StaircaseDependence(false, dependence(j)) + end end - deps[i] = StaircaseDependence(std, dependence(mono)) + end + end + I = findall(keep) + bd = BasisDependence(full_basis[I], d[I]) + if isempty(new_basis) + return bd + else + return MB.merge_bases( + bd, + BasisDependence(MB.MonomialBasis(new_basis), new_deps), + ) + end +end + +function MB.merge_bases(b1::BasisDependence, b2::BasisDependence) + basis, I1, I2 = MB.merge_bases(b1.basis, b2.basis) + deps = Vector{StaircaseDependence}(undef, length(basis)) + for i in eachindex(basis) + if iszero(I1[i]) + deps[i] = b2.dependence[I2[i]] else - deps[i] = d[I1[i]] + deps[i] = b1.dependence[I1[i]] end end - return BasisDependence(full_basis, deps) + return BasisDependence(basis, deps) end function _exponents(monos, i) diff --git a/test/null.jl b/test/null.jl index 77b2223..57c5528 100644 --- a/test/null.jl +++ b/test/null.jl @@ -19,7 +19,7 @@ function test_partial_commutative_fix(x, y) 1 0 0 0 # x^2 ] monos = MP.monomials((x, y), 0:2) - for basis in [MB.MonomialBasis(monos), MB.OrthonormalCoefficientsBasis(monos)] + for basis in [MB.MonomialBasis(monos)] #, MB.OrthonormalCoefficientsBasis(monos)] null = MM.MacaulayNullspace(matrix, basis, 1e-8) D = MM.StaircaseDependence solver = MM.StaircaseSolver{Float64}(max_partial_iterations = 1) From 9d629d8ebe5fb7cfcce1f74ba6e9d81a5a1d25ea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 30 Nov 2023 14:51:43 +0100 Subject: [PATCH 3/3] WIP --- test/null.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/null.jl b/test/null.jl index 57c5528..77b2223 100644 --- a/test/null.jl +++ b/test/null.jl @@ -19,7 +19,7 @@ function test_partial_commutative_fix(x, y) 1 0 0 0 # x^2 ] monos = MP.monomials((x, y), 0:2) - for basis in [MB.MonomialBasis(monos)] #, MB.OrthonormalCoefficientsBasis(monos)] + for basis in [MB.MonomialBasis(monos), MB.OrthonormalCoefficientsBasis(monos)] null = MM.MacaulayNullspace(matrix, basis, 1e-8) D = MM.StaircaseDependence solver = MM.StaircaseSolver{Float64}(max_partial_iterations = 1)