diff --git a/src/border.jl b/src/border.jl index 66a4ee8..1e7c33b 100644 --- a/src/border.jl +++ b/src/border.jl @@ -37,13 +37,26 @@ function BorderBasis{AnyDependence}(b::BorderBasis{<:StaircaseDependence}) end function BorderBasis{StaircaseDependence}(b::BorderBasis{<:AnyDependence}) - d = StaircaseDependence(_ -> true, b.dependent) + basis, I, _ = + MB.merge_bases(b.dependence.independent, b.dependence.dependent) + d = StaircaseDependence(basis) do i + return iszero(I[i]) + end dependent = convert(AnyDependence, d).dependent rows = _indices(b.dependence.independent, d.standard) cols = _indices(b.dependence.dependent, dependent) return BorderBasis(d, b.matrix[rows, cols]) end +function solve( + b::BorderBasis{<:AnyDependence,T}, + solver::SemialgebraicSets.AbstractMultiplicationMatricesSolver = MultivariateMoments.SemialgebraicSets.ReorderedSchurMultiplicationMatricesSolver{ + T, + }(), +) where {T} + return solve(BorderBasis{StaircaseDependence}(b), solver) +end + function solve( b::BorderBasis{<:StaircaseDependence,T}, solver::SemialgebraicSets.AbstractMultiplicationMatricesSolver = MultivariateMoments.SemialgebraicSets.ReorderedSchurMultiplicationMatricesSolver{ @@ -52,8 +65,9 @@ function solve( ) where {T} d = b.dependence dependent = convert(AnyDependence, d).dependent - vars = MP.variables(d.standard) - m = length(d.standard) + vars = MP.variables(d) + standard, _, I = MB.merge_bases(d.trivial_standard, d.standard) + m = length(standard) # If a monomial `border` is not in `dependent` and it is not a corner # then it can be divided by another monomial `x^α in dependent`. # So there exists a monomial `x^β` such that `border = x^(α + β)`. @@ -73,18 +87,24 @@ function solve( mult = Matrix{T}[zeros(T, m, m) for _ in eachindex(vars)] completed_border = Dict{eltype(dependent.monomials),Vector{T}}() function known_border_coefficients(border) - return !isnothing(_index(d.standard, border)) || + return !isnothing(_index(standard, border)) || !isnothing(_index(dependent, border)) || haskey(completed_border, border) end function border_coefficients(border) - k = _index(d.standard, border) + k = _index(standard, border) if !isnothing(k) return SparseArrays.sparsevec([k], [one(T)], m) end k = _index(dependent, border) if !isnothing(k) - return b.matrix[:, k] + v = zeros(T, m) + for i in eachindex(I) + if !iszero(I[i]) + v[i] = b.matrix[I[i], k] + end + end + return v end if haskey(completed_border, border) return completed_border[border] @@ -119,7 +139,7 @@ function solve( # monomial order so that we know that if `try_add_to_border` # fails, it will fail again if we run this for loop again with the same `o`. # That allows to limit the number of iteration of the outer loop by `length(vars)` - for std in d.standard.monomials + for std in standard.monomials for (k, shift) in enumerate(vars) if k in view(order, 1:o) continue @@ -131,12 +151,12 @@ function solve( if k in view(order, 1:o) continue end - if all(d.standard.monomials) do std + if all(standard.monomials) do std return known_border_coefficients(shift * std) end o += 1 order[o] = k - for (col, std) in enumerate(d.standard.monomials) + for (col, std) in enumerate(standard.monomials) mult[k][:, col] = border_coefficients(shift * std) end end @@ -229,9 +249,12 @@ function solve(b::BorderBasis{E}, solver::AlgebraicBorderSolver{D}) where {D,E} end """ - struct AlgebraicFallbackBorderSolver{M<:Union{Nothing,SS.AbstractMultiplicationMatricesSolver},S<:Union{Nothing,SS.AbstractAlgebraicSolver}} + struct AlgebraicFallbackBorderSolver{ + M<:Union{Nothing,SS.AbstractMultiplicationMatricesSolver}, + S<:AlgebraicBorderSolver, + } multiplication_matrices_solver::M - algebraic_solver::A + algebraic_solver::S end Solve with `multiplication_matrices_solver` and if it fails, falls back to diff --git a/src/dependence.jl b/src/dependence.jl index 5a08402..9dcda78 100644 --- a/src/dependence.jl +++ b/src/dependence.jl @@ -26,15 +26,20 @@ end function Base.show(io::IO, d::AnyDependence) println(io, "AnyDependence") - println(io, "with independent basis:") - println(io, d.independent) - println(io, "and dependent basis:") - println(io, d.dependent) + if !isempty(d.independent.monomials) + println(io, "with independent basis:") + println(io, d.independent) + end + if !isempty(d.dependent.monomials) + println(io, "and dependent basis:") + println(io, d.dependent) + end return end """ struct StaircaseDependence{B<:MB.AbstractPolynomialBasis} <: AbstractDependence + trivial_standard::B standard::B corners::B dependent_border::B @@ -48,6 +53,7 @@ by multiplying an independent with a variable. If the number of variables is 2 or 3, it can be visualized with Plots.jl. """ struct StaircaseDependence{B<:MB.AbstractPolynomialBasis} <: AbstractDependence + trivial_standard::B standard::B corners::B dependent_border::B @@ -63,14 +69,26 @@ end function Base.show(io::IO, d::StaircaseDependence) println(io, "StaircaseDependence") - println(io, "with standard basis:") - println(io, d.standard) - println(io, "and corners basis:") - println(io, d.corners) - println(io, "and dependent border basis:") - println(io, d.dependent_border) - println(io, "and independent border basis:") - print(io, d.independent_border) + if !isempty(d.trivial_standard.monomials) + println(io, "with trivial standard basis:") + println(io, d.trivial_standard) + end + if !isempty(d.standard.monomials) + println(io, "with standard basis:") + println(io, d.standard) + end + if !isempty(d.corners.monomials) + println(io, "and corners basis:") + println(io, d.corners) + end + if !isempty(d.dependent_border.monomials) + println(io, "and dependent border basis:") + println(io, d.dependent_border) + end + if !isempty(d.independent_border.monomials) + println(io, "and independent border basis:") + print(io, d.independent_border) + end return end @@ -79,6 +97,26 @@ function Base.convert(::Type{AnyDependence}, d::StaircaseDependence) return AnyDependence(d.standard, dependent) end +function AnyDependence( + is_dependent::Function, + basis::MB.MonomialBasis{M}, +) where {M} + independent = M[] + dependent = M[] + for i in eachindex(basis.monomials) + mono = basis.monomials[i] + if is_dependent(i) + push!(dependent, mono) + else + push!(independent, mono) + end + end + return AnyDependence( + MB.MonomialBasis(independent), + MB.MonomialBasis(dependent), + ) +end + """ function StaircaseDependence( is_dependent::Function, @@ -98,6 +136,7 @@ function StaircaseDependence( is_dependent::Function, basis::MB.MonomialBasis{M}, ) where {M} + trivial_standard = M[] standard = M[] corners = M[] vars = MP.variables(basis) @@ -106,7 +145,9 @@ function StaircaseDependence( # It also ensures that the standard monomials have the "staircase structure". if !any(Base.Fix2(MP.divides, mono), corners) i = _index(basis, mono) - if isnothing(i) || is_dependent(i) + if isnothing(i) + push!(trivial_standard, mono) + elseif !is_dependent(i) push!(standard, mono) else push!(corners, mono) @@ -121,7 +162,7 @@ function StaircaseDependence( if isnothing(_monomial_index(standard, border)) && isnothing(_monomial_index(corners, border)) i = _index(basis, mono) - if isnothing(i) || is_dependent(i) + if isnothing(i) || !is_dependent(i) push!(independent_border, border) else push!(dependent_border, border) @@ -130,6 +171,7 @@ function StaircaseDependence( end end return StaircaseDependence( + MB.MonomialBasis(trivial_standard), MB.MonomialBasis(standard), MB.MonomialBasis(corners), MB.MonomialBasis(dependent_border), @@ -158,36 +200,37 @@ function _deg(deg, b, args...) end end -__combine(::Function, ::Nothing) = nothing -__combine(::Function, a) = a -__combine(f::Function, ::Nothing, args...) = __combine(f, args...) -function __combine(f::Function, a, args...) - d = __combine(f, args...) +__reduce(::Function, ::Nothing) = nothing +__reduce(::Function, a) = a +__reduce(f::Function, ::Nothing, args...) = __reduce(f, args...) +function __reduce(f::Function, a, args...) + d = __reduce(f, args...) if isnothing(d) return a else f(a, d) end end -function _combine(f, args...) - d = __combine(f, args...) +function _reduce(f, args...) + d = __reduce(f, args...) if isnothing(d) error("Cannot compute `$(f)degree` as all bases are empty") end return d end -function _combine_deg(combine, deg, d::AnyDependence, args...) - return _combine( +function _map_reduce(combine, deg, d::AnyDependence, args...) + return _reduce( combine, _deg(deg, d.independent, args...), _deg(deg, d.dependent, args...), ) end -function _combine_deg(combine, deg, d::StaircaseDependence, args...) - return _combine( +function _map_reduce(combine, deg, d::StaircaseDependence, args...) + return _reduce( combine, + _deg(deg, d.trivial_standard, args...), _deg(deg, d.standard, args...), _deg(deg, d.corners, args...), _deg(deg, d.dependent_border, args...), @@ -195,12 +238,20 @@ function _combine_deg(combine, deg, d::StaircaseDependence, args...) ) end +function _reduce_variables(v, w) + return MP.variables(prod(v) * prod(w)) +end + +function MP.variables(d::AbstractDependence) + return _map_reduce(_reduce_variables, MP.variables, d) +end + function MP.mindegree(d::AbstractDependence, args...) - return _combine_deg(min, MP.mindegree, d, args...) + return _map_reduce(min, MP.mindegree, d, args...) end function MP.maxdegree(d::AbstractDependence, args...) - return _combine_deg(max, MP.maxdegree, d, args...) + return _map_reduce(max, MP.maxdegree, d, args...) end function _ticks(d::AbstractDependence, v) @@ -216,17 +267,21 @@ RecipesBase.@recipe function f(m::AnyDependence) if length(t) >= 3 zticks --> t[3] end - RecipesBase.@series begin - seriestype --> :scatter - markershape --> :circle - label := "Independent" - _split_exponents(m.independent) + if !isempty(m.independent.monomials) + RecipesBase.@series begin + seriestype --> :scatter + markershape --> :circle + label := "Independent" + _split_exponents(m.independent) + end end - RecipesBase.@series begin - seriestype --> :scatter - markershape --> :rect - label := "Dependent" - _split_exponents(m.dependent) + if !isempty(m.dependent.monomials) + RecipesBase.@series begin + seriestype --> :scatter + markershape --> :rect + label := "Dependent" + _split_exponents(m.dependent) + end end end @@ -239,28 +294,44 @@ RecipesBase.@recipe function f(m::StaircaseDependence) if length(t) >= 3 zticks --> t[3] end - RecipesBase.@series begin - seriestype --> :scatter - markershape --> :circle - label := "Standard" - _split_exponents(m.standard) + if !isempty(m.trivial_standard.monomials) + RecipesBase.@series begin + seriestype --> :scatter + markershape --> :circle + label := "Trivial standard" + _split_exponents(m.trivial_standard) + end + end + if !isempty(m.standard.monomials) + RecipesBase.@series begin + seriestype --> :scatter + markershape --> :circle + label := "Standard" + _split_exponents(m.standard) + end end - RecipesBase.@series begin - seriestype --> :scatter - markershape --> :rect - label := "Corners" - _split_exponents(m.corners) + if !isempty(m.corners.monomials) + RecipesBase.@series begin + seriestype --> :scatter + markershape --> :rect + label := "Corners" + _split_exponents(m.corners) + end end - RecipesBase.@series begin - seriestype --> :scatter - markershape --> :rect - label := "Dependent border" - _split_exponents(m.dependent_border) + if !isempty(m.dependent_border.monomials) + RecipesBase.@series begin + seriestype --> :scatter + markershape --> :rect + label := "Dependent border" + _split_exponents(m.dependent_border) + end end - RecipesBase.@series begin - seriestype --> :scatter - markershape --> :circle - label := "Independent border" - _split_exponents(m.independent_border) + if !isempty(m.independent_border.monomials) + RecipesBase.@series begin + seriestype --> :scatter + markershape --> :circle + label := "Independent border" + _split_exponents(m.independent_border) + end end end diff --git a/src/echelon.jl b/src/echelon.jl index 1c3cc1d..5dd02db 100644 --- a/src/echelon.jl +++ b/src/echelon.jl @@ -26,7 +26,9 @@ function build_system(U::AbstractMatrix, basis::MB.MonomialBasis) end """ - struct Echelon + struct Echelon{D<:AbstractDependence,S<:Union{Nothing,SemialgebraicSets.AbstractAlgebraicSolver}} <: + MacaulayNullspaceSolver + solver::S end Given a [`MacaulayNullspace`](@ref), computes its echelon form @@ -78,15 +80,23 @@ Foundations of Computational Mathematics 8 (2008): 607-647. *Numerical polynomial algebra.* Society for Industrial and Applied Mathematics, 2004. """ -struct Echelon{S<:Union{Nothing,SemialgebraicSets.AbstractAlgebraicSolver}} <: - MacaulayNullspaceSolver +struct Echelon{ + D<:AbstractDependence, + S<:Union{Nothing,SemialgebraicSets.AbstractAlgebraicSolver}, +} <: MacaulayNullspaceSolver + fallback::Bool solver::S end -Echelon() = Echelon(nothing) - -function border_basis_solver(solver::Echelon) - return AlgebraicBorderSolver{AnyDependence}(solver.solver) +function Echelon{D}(; + fallback::Bool = false, + solver::Union{Nothing,SemialgebraicSets.AbstractAlgebraicSolver} = nothing, +) where {D} + return Echelon{D,typeof(solver)}(fallback, solver) end +Echelon(; kws...) = Echelon{AnyDependence}(; kws...) + +# TODO remove +Echelon(solver) = Echelon(; solver) # TODO remove in next breaking release function compute_support!( @@ -101,7 +111,10 @@ end import RowEchelon -function border_basis_and_solver(null::MacaulayNullspace, solver::Echelon) +function border_basis_and_solver( + null::MacaulayNullspace, + solver::Echelon{D}, +) where {D} # Ideally, we should determine the pivots with a greedy sieve algorithm [LLR08, Algorithm 1] # so that we have more chance that low order monomials are in β and then more chance # so that the pivots form an order ideal. We just use `rref` which does not implement the sieve @@ -118,9 +131,14 @@ function border_basis_and_solver(null::MacaulayNullspace, solver::Echelon) RowEchelon.rref!(Z, null.accuracy / sqrt(size(Z, 2))) #r, vals = solve_system(U', μ.x) # TODO determine what is better between rank_check and sqrt(rank_check) here - solver = AlgebraicBorderSolver{AnyDependence}( + algebraic_solver = AlgebraicBorderSolver{D}( algorithm = SS.Buchberger(√null.accuracy), solver = solver.solver, ) + if solver.fallback + solver = AlgebraicFallbackBorderSolver(; algebraic_solver) + else + solver = algebraic_solver + end return build_system(Z, null.basis), solver end diff --git a/src/flat.jl b/src/flat.jl index 249a2e7..544e7ef 100644 --- a/src/flat.jl +++ b/src/flat.jl @@ -6,6 +6,10 @@ end SemialgebraicSets.is_zero_dimensional(::ZeroDimensionalVariety) = true Base.length(v::ZeroDimensionalVariety) = length(v.elements) Base.iterate(v::ZeroDimensionalVariety, args...) = iterate(v.elements, args...) +function Base.show(io::IO, V::ZeroDimensionalVariety) + println(io, "ZeroDimensionalVariety with elements:") + return show(io, V.elements) +end # Decomposition of the pencil of matrices function decompose( diff --git a/src/shift.jl b/src/shift.jl index 574a905..eacdefb 100644 --- a/src/shift.jl +++ b/src/shift.jl @@ -11,7 +11,7 @@ RankDependence(matrix, check) = RankDependence(matrix, check, Int[], 0) function is_dependent!(r::RankDependence, row) if r.old_rank == size(r.matrix, 2) - return false + return true end rows = vcat(r.independent_rows, row) new_rank = LinearAlgebra.rank(r.matrix[rows, :], r.check) @@ -26,7 +26,7 @@ function is_dependent!(r::RankDependence, row) if independent r.independent_rows = rows end - return independent + return !independent end function AnyDependence(null::MacaulayNullspace, rank_check::RankCheck) @@ -39,23 +39,12 @@ function StaircaseDependence(null::MacaulayNullspace, rank_check::RankCheck) return StaircaseDependence(Base.Fix1(is_dependent!, r), null.basis) end -function _indices_or_ignore(in::MB.MonomialBasis, from::MB.MonomialBasis) - indices = Int[] - for mono in from.monomials - i = _index(in, mono) - if !isnothing(i) - push!(indices, i) - end - end - return indices -end - function _indices(in::MB.MonomialBasis, from::MB.MonomialBasis) return Int[_index(in, mono) for mono in from.monomials] end function BorderBasis(d::AnyDependence, null::MacaulayNullspace) - indep_rows = _indices_or_ignore(null.basis, d.independent) + indep_rows = _indices(null.basis, d.independent) dep_rows = _indices(null.basis, d.dependent) @assert length(indep_rows) == size(null.matrix, 2) return BorderBasis( @@ -64,12 +53,17 @@ function BorderBasis(d::AnyDependence, null::MacaulayNullspace) ) end +function column_compression(A, rows) + return LinearAlgebra.svd(A[rows, :]).U +end + function BorderBasis(d::StaircaseDependence, null::MacaulayNullspace) - indep_rows = _indices_or_ignore(null.basis, d.standard) + indep_rows = _indices(null.basis, d.standard) dependent = convert(AnyDependence, d).dependent dep_rows = _indices(null.basis, dependent) - if length(indep_rows) < size(null.matrix, 2) - error("Column compression not supported yet") + U = null.matrix + if length(indep_rows) < size(U, 2) + U = column_compression(U, indep_rows) end return BorderBasis( d, diff --git a/test/commutativetests.jl b/test/commutativetests.jl index a5d9ec0..a8d868b 100644 --- a/test/commutativetests.jl +++ b/test/commutativetests.jl @@ -6,3 +6,4 @@ include("extract.jl") include("atomic.jl") include("hermitian_poly.jl") include("dependence.jl") +include("null.jl") diff --git a/test/dependence.jl b/test/dependence.jl index 2b16fd2..b40107b 100644 --- a/test/dependence.jl +++ b/test/dependence.jl @@ -12,8 +12,10 @@ function test_degree_error(x, y, z) M = typeof(x * y * z) m = b(M[]) a = MM.AnyDependence(m, m) + @test sprint(show, a) == "AnyDependence\n" @test isempty(a) - s = MM.StaircaseDependence(m, m, m, m) + s = MM.StaircaseDependence(m, m, m, m, m) + @test sprint(show, s) == "StaircaseDependence\n" @test isempty(s) for (f, name) in [(MP.mindegree, "min"), (MP.maxdegree, "max")] err = ErrorException( @@ -27,6 +29,7 @@ function test_degree_error(x, y, z) end function _test_recipe(dep, ticks, args, names, indep) + @test sprint(show, dep) isa String d = Dict{Symbol,Any}() r = RB.apply_recipe(d, dep) @test length(d) == 1 + length(ticks) @@ -70,11 +73,23 @@ function test_recipe(x, y, z) [true, false], ) return _test_recipe( - MM.StaircaseDependence(b(a .* z^0), b(c .* z^0), b(d), b(e)), + MM.StaircaseDependence( + b(a .* z^0), + b([x^0 * y^0 * z]), + b(c .* z^0), + b(d), + b(e), + ), [0:1, 0:3, 0:4], - [(A..., [0]), (C..., [0, 0]), D, E], - ["Standard", "Corners", "Dependent border", "Independent border"], - [true, false, false, true], + [(A..., [0]), ([0], [0], [1]), (C..., [0, 0]), D, E], + [ + "Trivial standard", + "Standard", + "Corners", + "Dependent border", + "Independent border", + ], + [true, true, false, false, true], ) end diff --git a/test/null.jl b/test/null.jl new file mode 100644 index 0000000..84aeca2 --- /dev/null +++ b/test/null.jl @@ -0,0 +1,62 @@ +module TestNull + +using Test +import MultivariatePolynomials as MP +import MultivariateBases as MB +import MultivariateMoments as MM + +b(x) = MB.MonomialBasis(x) + +function test_univariate_infinity(x) + Z = Float64[ + 1 0 + 2 0 + 0 1 + ] + for (null, sol) in [ + (MM.MacaulayNullspace(Z, b([1, x, x^2]), 1e-8), 2), + (MM.MacaulayNullspace(Z[1:2, 1:1], b([1, x]), 1e-8), 2), + (MM.MacaulayNullspace(zeros(1, 0), b([x]), 1e-8), 0), + ] + @testset "$D" for D in [MM.AnyDependence, MM.StaircaseDependence] + @testset "$name" for (solver, name) in [ + (MM.ShiftNullspace{D}(), "shift"), + (MM.Echelon{D}(fallback = false), "echelon no fallback"), + (MM.Echelon{D}(fallback = true), "echelon fallback"), + ] + if name == "echelon no fallback" && size(null.matrix) == (1, 0) + # FIXME The system is `x = 0` hence it thing that it is homogeneous and it finds no sol + # The homogeneous check should be better... It should probably be explicitly stated upfront + continue + end + @test collect(MM.solve(null, solver)) == [[sol]] + end + end + end +end + +function runtests(args...) + for name in names(@__MODULE__; all = true) + if startswith("$name", "test_") + @testset "$(name)" begin + getfield(@__MODULE__, name)(args...) + end + end + end +end + +end + +using Test + +import DynamicPolynomials +@testset "DynamicPolynomials" begin + DynamicPolynomials.@polyvar x + TestNull.runtests(x) +end + +import TypedPolynomials +@testset "TypedPolynomials" begin + TypedPolynomials.@polyvar x + TestNull.runtests(x) +end