diff --git a/src/border.jl b/src/border.jl index 0704b32..a9e515d 100644 --- a/src/border.jl +++ b/src/border.jl @@ -71,18 +71,21 @@ struct StaircaseSolver{ max_iterations::Int rank_check::R solver::M + print_level::Int end function StaircaseSolver{T}(; max_partial_iterations::Int = 0, max_iterations::Int = -1, rank_check::RankCheck = LeadingRelativeRankTol(Base.rtoldefault(T)), solver = SS.ReorderedSchurMultiplicationMatricesSolver{T}(), + print_level::Int = 1, ) where {T} return StaircaseSolver{T,typeof(rank_check),typeof(solver)}( max_partial_iterations, max_iterations, rank_check, solver, + print_level, ) end @@ -225,7 +228,11 @@ function solve( if solver.max_iterations == 0 Uperp = nothing else - Uperp = commutation_fix(mult, solver.solver.ε) + Uperp = commutation_fix( + mult, + √solver.solver.ε; + print_level = solver.print_level, + ) end com_fix = if isnothing(Uperp) nothing @@ -267,7 +274,7 @@ function solve( end end -function commutation_fix(matrices, ε) +function commutation_fix(matrices, ε; print_level) if isempty(first(matrices)) return end @@ -301,6 +308,11 @@ function commutation_fix(matrices, ε) if isnothing(leading_F) return else + if print_level >= 1 + @info( + "Found multiplication matrices with commutation error of `$leading_ratio` which is larger than the tolerance of `$ε`. Adding this to the equations and continuing." + ) + end return leading_F.U[:, 2:end] end end diff --git a/src/shift.jl b/src/shift.jl index 012e063..a77c0d8 100644 --- a/src/shift.jl +++ b/src/shift.jl @@ -106,31 +106,47 @@ IFAC Proceedings Volumes 45.16 (2012): 1203-1208. struct ShiftNullspace{D,S,C<:RankCheck} <: MacaulayNullspaceSolver solver::S check::C + print_level::Int end # Because the matrix is orthogonal, we know the SVD of the whole matrix is # `ones(...)` so an `AbsoluteRankTol` would be fine here. # However, since we also know that the first row (which correspond to the # constant monomial) should be a standard monomial, `LeadingRelativeRankTol` # ensures that we will take it. -function ShiftNullspace{D}(check::RankCheck) where {D} - return ShiftNullspace{D,Nothing,typeof(check)}(nothing, check) +function ShiftNullspace{D}(check::RankCheck; print_level = 1) where {D} + return ShiftNullspace{D,Nothing,typeof(check)}(nothing, check, print_level) end -function ShiftNullspace{D}(solver::StaircaseSolver) where {D} +function ShiftNullspace{D}(solver::StaircaseSolver; print_level = 1) where {D} check = solver.rank_check - return ShiftNullspace{D,typeof(solver),typeof(check)}(solver, check) + return ShiftNullspace{D,typeof(solver),typeof(check)}( + solver, + check, + print_level, + ) end -function ShiftNullspace{D}() where {D} - return ShiftNullspace{D}(LeadingRelativeRankTol(Base.rtoldefault(Float64))) +function ShiftNullspace{D}(; kws...) where {D} + return ShiftNullspace{D}( + LeadingRelativeRankTol(Base.rtoldefault(Float64)); + kws..., + ) end ShiftNullspace(args...) = ShiftNullspace{StaircaseDependence}(args...) function _solver( + null::MacaulayNullspace, shift::ShiftNullspace{StaircaseDependence,Nothing}, ::Type{T}, ) where {T} - return StaircaseSolver{T}(rank_check = shift.check) + return StaircaseSolver{T}(; + rank_check = shift.check, + solver = SS.ReorderedSchurMultiplicationMatricesSolver( + convert(T, null.accuracy), + ), + print_level = shift.print_level, + ) end -function _solver(shift::ShiftNullspace, ::Type) + +function _solver(::MacaulayNullspace, shift::ShiftNullspace, ::Type) return shift.solver end @@ -138,5 +154,5 @@ function border_basis_and_solver( null::MacaulayNullspace{T}, shift::ShiftNullspace{D}, ) where {T,D} - return BorderBasis{D}(null, shift.check), _solver(shift, T) + return BorderBasis{D}(null, shift.check), _solver(null, shift, T) end diff --git a/test/extract.jl b/test/extract.jl index b189f5d..53eda31 100644 --- a/test/extract.jl +++ b/test/extract.jl @@ -362,6 +362,29 @@ function large_norm(rank_check) @test nothing === atomic_measure(ν, rank_check) end +function no_com(solver, T) + Mod.@polyvar x[1:2] + Q = T[ + 1.0000000000000002 2.9999999777389235 1.408602019734018 8.999999911981313 4.225805987406138 3.043010098670093 + 2.9999999777389235 8.999999911981313 4.225805987406138 26.999999824988187 12.677417999642149 9.129030026074997 + 1.408602019734018 4.225805987406138 3.043010098670093 12.677417999642149 9.129030026074997 9.580642414414385 + 8.999999911981313 26.999999824988187 12.677417999642149 80.99999955990646 38.03225428611012 27.387090350285558 + 4.225805987406138 12.677417999642149 9.129030026074997 38.03225428611012 27.387090350285558 28.74192618075039 + 3.043010098670093 9.129030026074997 9.580642414414385 27.387090350285558 28.74192618075039 35.73117167739159 + ] + ν = moment_matrix(Q, monomials(x, 0:2)) + η = AtomicMeasure( + x, + [ + WeightedDiracMeasure(T[1, 3], T(0.863799337)), + WeightedDiracMeasure(T[4, 3], T(0.136200673)), + ], + ) + rank_check = LeadingRelativeRankTol(1e-7) + atoms = atomic_measure(ν, rank_check, solver) + @test atoms ≈ η rtol = 1e-4 +end + _short(x::String) = x[1:min(10, length(x))] _short(x) = _short(string(x)) @@ -443,6 +466,12 @@ function test_extract() jcg14_6_1(6e-3) jcg14_6_1(8e-4, false) large_norm(1e-2) + @testset "No comm fix $(_short(solver)) $T" for solver in + [ShiftNullspace()], + T in [Float64] + #, BigFloat] + no_com(solver, T) + end return end