Large refactoring with StarAlgebras

.github/workflows/ci.yml

@@ -41,7 +41,7 @@ jobs:
         run: |
           using Pkg
           Pkg.add([
-              PackageSpec(name="DynamicPolynomials", rev="master"),
+              PackageSpec(name="DynamicPolynomials", rev="release-0.7"),
               PackageSpec(name="TypedPolynomials", rev="master"),
               PackageSpec(name="StarAlgebras", rev="bl/term"),
           ])

AGENTS.md

@@ -0,0 +1,57 @@
+
+This repository as well as DynamicPolynomials have accumulated a lot of technical depth due to bad design decisions and it is time to get rid of it.
+The bad design decisions were
+1. Try to support `*(::MP.AbstractPolynomialLike, ::Any)` this creates a lot of invalidation and promotion rules are a mess. We should only support `*(::MP.AbstractPolynomialLike{T}, ::T)` and have `*(::MP.AbstractPolynomialLike, ::Number)` that convert the number to `T`. Same for other operators like `+`, `-` etc..
+2. Implement `+(p::MP.AbstractPolynomialLike, q::MP.AbstractPolynomialLike)` where `p` and `q` have different variables. Not simply by first promoting them to the same set of variables and then use a simple implementation of the sum but instead to have some complicated implementation dealing with monomials over different variables. We should just directly promote with `StarAlgebras.promote_bases` now and only implement `+` over polynomials of the same algebra
+3. Having different implementation of polynomials, `MP.Polynomial` using a list of terms and `DynamicPolynomials.Polynomial` having a separate list of coefficients and list of monomials. We should just use StarAlgebras.AlgebraElement now.
+
+So the plan is:
+1. Move `MP.Term` to `StarAlgebras.Term`, so `MP.AbstractTermLike` won't be an abstract type anymore but a union of `AbstractMonomialLike` and `SA.Term`
+2. Remove `MP.Polynomial`, and replace it by the implementation of polynomials done in MultivariateBases using StarAlgebras. Doing so, we can just merge MultivariateBases into MultivariatePolynomials, `MP.AbstractPolynomialLike` will then be a union of that StarAlgebras.AlgebraElement and `AbstractTermLike`.
+
+Terminology:
+
+MOI: MathOptInterface
+MA: MutableArithmetics
+DP: DynamicPolynomials
+TP: TypedPolynomials
+MP: MultivariatePolynomials
+SA: StarAlgebras
+MB: MultivariateBases
+
+All packages should be in ~/.julia/dev
+
+MA has a large test suite that is used by JuMP, MOI, MP and therefore also DP and TP. It is in MP/test. It's nice but it takes a lot of time to run. First makes sure everything is done before running it for the final touches.
+This refactor is supposed to get rid of a lot of code and technical depth, don't start writing a lot of code before first asking me.
+
+The promotion rules are a mess. Because we had to support promotion with Any, this messed up badly with Julia's internals.
+In Julia, normally you just need to implement promote_rule(::A, ::B), not promote_rule(::B, ::A). But because I implemented promotion with ::Any, I had to do both and it created such a mess!!
+It was also a rabbithole where I then needed also APL{<:Any} and things like that, quite complicated...
+I want things to be easier now. Before, I belived if x is a variable, promote_type(typeof(x), typeof("x")) would be something like Term{String}. We don't want that anymore, let it just be Any[].
+You see the vibe ? Let's start simple and then be really careful for test we used to have whether we really want this test to still pass or we just want to simplify.
+
+About backward compat, we are going to make a breaking release so I prefer being more breaking and having a simpler code than the opposite.
+
+A lot of code in MP are like sum(t::Vector{<:AbstractTerm}) = dot(coefficient.(t), monomial.(t)) and then dot(c::AbstractVector, m::AbstractVector{<:AbstractMonomialLike}) = sum(c .* m).
+These create stackoverflow in case DP or TP forget to implement one of the two! This was because, as explained above, TP used the default polynomials that were a vector of terms and DP
+was a vector of coefficients separated to a vector of monomials. Now it all gets simpler because it will be a separate vector of coefficients and monomials since we'll just be using SA.AlgebraElement!!
+So we can simplify. Also, because SA.AlgebraElement <: MA.AbstractMutable, we already have a fallback for Base.sum so we might just be able to remove these!
+These were written before MutableArithmetics existed. I am maintaining MA, SA, JuMP, MOI, MP, DP, TP and I don't want to maintain duplicates. All these packages define mutable objects.
+MA allows me to have already a lot of code in common of all of them. Then, JuMP, SA and MOI will have different implementations of basically "sum of terms" but that's fine.
+What I don't want is MP, TP and DP to also have their own version, they should just use SA.AlgebraElement!
+
+It has already been started by another AI agent on the branches:
+DP: bl/sa_term
+MP: bl/sa_poly
+SA: bl/term
+
+Just continue on these branches, feel free to simplify them if you see anything better, it was done a while ago by an older AI, you are smarter ;)
+
+In the printing of stack-traces, you can see that because we just have generic types in SA that we parametrize, things are getting very long and it's getting very difficult to debug.
+This gives the user extra flexibility to try new these, and in these cases, it will be nice to have precise stacktraces, but we also want the common cases (like what you get with DP.@polyvar x y; 2 * x + y)
+to have a very small type like DP.Polynomial{Int}
+We can solve it easily by having a "const Polynomial = ..." in DynamicPolynomials, don't hesitate to do this early on, it will help you be more context-efficient.
+
+Note that MB currently already defines a polynomial using SA.AlgebraElement. So we can just basically steal the code he has.
+He also has code that handles the interaction between MP.APL and SA.AlgebraElement, this will go away since everything will be an AlgebraElement now !
+This means that in MB, a lot of code will go away since it will move to MP, what will be left is just defining new bases essentially.

Project.toml

@@ -5,7 +5,6 @@ version = "0.5.19"
 repo = "https://github.com/JuliaAlgebra/MultivariatePolynomials.jl"
 
 [deps]
-DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
@@ -18,7 +17,6 @@ MultivariatePolynomialsChainRulesCoreExt = "ChainRulesCore"
 
 [compat]
 ChainRulesCore = "1"
-DataStructures = "0.19"
 MutableArithmetics = "0.3, 1"
 StarAlgebras = "0.3"
 julia = "1.10"

docs/src/types.md

@@ -95,14 +95,14 @@ leading_term
 leading_coefficient
 leading_monomial
 deg_num_leading_terms
-remove_leading_term
+SA.remove_leading_term
 remove_monomials
 filter_terms
 OfDegree
 monic
 map_coefficients
-map_coefficients!
-map_coefficients_to!
+SA.map_coefficients!
+SA.map_coefficients_to!
 conj(::_APL)
 real(::_APL)
 imag(::_APL)

src/MultivariatePolynomials.jl

@@ -2,66 +2,11 @@ module MultivariatePolynomials
 
 import LinearAlgebra
 
-import DataStructures
-
 import MutableArithmetics as MA
 
 import StarAlgebras as SA
 
-"""
-    AbstractPolynomialLike{T}
-
-Abstract type for a value that can act like a polynomial. For instance, an
-`AbstractTerm{T}` is an `AbstractPolynomialLike{T}` since it can act as a
-polynomial of only one term.
-"""
-abstract type AbstractPolynomialLike{T} <: MA.AbstractMutable end
-
-"""
-    AbstractMonomialLike
-
-Abstract type for a value that can act like a monomial. For instance, an `AbstractVariable` is an `AbstractMonomialLike` since it can act as a monomial of one variable with degree `1`.
-"""
-abstract type AbstractMonomialLike <: AbstractPolynomialLike{Int} end
-
-"""
-    AbstractVariable <: AbstractMonomialLike
-
-Abstract type for a variable.
-"""
-abstract type AbstractVariable <: AbstractMonomialLike end
-
-"""
-    AbstractMonomial <: AbstractMonomialLike
-
-Abstract type for a monomial, i.e. a product of variables elevated to a nonnegative integer power.
-"""
-abstract type AbstractMonomial <: AbstractMonomialLike end
-
-"""
-    AbstractTermLike{T}
-
-Union type for values that can act like a term. This includes `AbstractMonomialLike`
-(which acts as a term with coefficient `1`) and `SA.Term{T}`.
-"""
-const AbstractTermLike{T} = Union{AbstractMonomialLike,SA.Term{T}}
-
-"""
-    AbstractTerm{T}
-
-Type alias for a term of coefficient type `T`, i.e. the product between a
-value of type `T` and a monomial. This is [`StarAlgebras.Term{T}`](@ref).
-"""
-const AbstractTerm{T} = SA.Term{T}
-
-"""
-    AbstractPolynomial{T} <: AbstractPolynomialLike{T}
-
-Abstract type for a polynomial of coefficient type `T`, i.e. a sum of `AbstractTerm{T}`s.
-"""
-abstract type AbstractPolynomial{T} <: AbstractPolynomialLike{T} end
-
-const _APL{T} = Union{AbstractPolynomialLike{T},SA.Term{T}}
+include("types.jl")
 
 include("zip.jl")
 include("lazy_iterators.jl")
@@ -91,9 +36,16 @@ include("division.jl")
 include("gcd.jl")
 include("det.jl")
 
-include("default_term.jl")
-include("sequences.jl")
-include("default_polynomial.jl")
+
+# Polynomial bases and algebra representation shared with MultivariateBases.
+include("mb_interface.jl")
+include("mb_variables.jl")
+include("mb_polynomial.jl")
+include("mb_bases.jl")
+include("mb_mstructures.jl")
+include("mb_monomial_basis.jl")
+include("mb_algebra.jl")
+include("mb_arithmetic.jl")
 
 include("deprecate.jl")
 

src/comparison.jl

@@ -51,6 +51,31 @@ end
 left_constant_eq(α, p::_APL; comp = (==)) = right_term_eq(p, α; comp)
 right_constant_eq(p::_APL, α; comp = (==)) = right_term_eq(p, α; comp)
 
+for comp in (:(==), :isequal)
+    @eval begin
+        function Base.$comp(
+            p::Union{AbstractPolynomial{T},AbstractTerm{T}},
+            a::Union{T,Number},
+        ) where {T}
+            return right_constant_eq(p, a; comp = $comp)
+        end
+        function Base.$comp(
+            a::Union{T,Number},
+            p::Union{AbstractPolynomial{T},AbstractTerm{T}},
+        ) where {T}
+            return left_constant_eq(a, p; comp = $comp)
+        end
+        Base.$comp(m::AbstractMonomialLike, a::Number) =
+            right_constant_eq(m, a; comp = $comp)
+        Base.$comp(a::Number, m::AbstractMonomialLike) =
+            left_constant_eq(a, m; comp = $comp)
+        Base.$comp(p::AbstractPolynomial, m::AbstractMonomialLike) =
+            right_term_eq(p, term(m); comp = $comp)
+        Base.$comp(m::AbstractMonomialLike, p::AbstractPolynomial) =
+            right_term_eq(p, term(m); comp = $comp)
+    end
+end
+
 function Base.:(==)(mono::AbstractMonomial, v::AbstractVariable)
     return isone(degree(mono)) && variable(mono) == v
 end
@@ -317,8 +342,9 @@ function ordering end
 ordering(::Type{<:AbstractMonomial}) = Graded{LexOrder}
 ordering(::Type{P}) where {P} = ordering(monomial_type(P))
 ordering(p::AbstractPolynomialLike) = ordering(typeof(p))
-# Useful for instance to ask ordering given the list
-# of variables
+# For type-level: derive ordering from Vector's element type
+ordering(::Type{<:AbstractVector{T}}) where {T} = ordering(T)
+# Useful for instance to ask ordering given the list of variables
 ordering(::AbstractVector{T}) where {T} = ordering(T)
 ordering(t::Tuple) = ordering(first(t))
 

src/complex.jl

@@ -203,7 +203,7 @@ for fun in [:real, :imag]
                         full_version,
                         real(coefficient_type(full_version)),
                     ),
-                    map_coefficients!($fun, full_version),
+                    SA.map_coefficients!($fun, full_version),
                 )
             end
             function Base.$fun(x::AbstractVector{<:AbstractMonomial})

src/conversion.jl

@@ -1,16 +1,12 @@
-function convert_constant end
-Base.convert(::Type{P}, α) where {P<:AbstractPolynomialLike} = convert_constant(P, α)
-function Base.convert(::Type{SA.Term{T,M}}, α) where {T,M}
-    return convert_constant(SA.Term{T,M}, α)
-end
-function convert_constant(::Type{SA.Term{T,M}}, α) where {T,M}
-    return term(convert(T, α), constant_monomial(SA.Term{T,M}))
-end
-function convert_constant(::Type{PT}, α) where {PT<:AbstractPolynomial}
-    return convert(PT, convert(term_type(PT), α))
+function Base.convert(
+    ::Type{P},
+    m::AbstractMonomialLike,
+) where {P<:AbstractPolynomial}
+    return convert(P, term(m))
 end
-function Base.convert(::Type{P}, p::_APL) where {T,P<:AbstractPolynomial{T}}
-    return error("`convert` not implemented for $P")
+
+function Base.convert(::Type{P}, c::Number) where {P<:AbstractPolynomial}
+    return convert(P, constant_term(c, P))
 end
 
 function Base.convert(
@@ -59,21 +55,7 @@ function Base.convert(
         throw(InexactError(:convert, T, p))
     end
 end
-# Disambiguation: SA.Term{T,M} from AbstractPolynomial (more specific than the α catch-all)
-function Base.convert(
-    ::Type{SA.Term{T,M}},
-    p::AbstractPolynomial,
-) where {T,M}
-    if iszero(nterms(p))
-        convert(SA.Term{T,M}, zero_term(p))
-    elseif isone(nterms(p))
-        convert(SA.Term{T,M}, leading_term(p))
-    else
-        throw(InexactError(:convert, SA.Term{T,M}, p))
-    end
-end
-
-MA.scaling(p::AbstractPolynomialLike{T}) where {T} = convert(T, p)
+MA.scaling(p::AbstractPolynomialLike) = convert(coefficient_type(p), p)
 # Conversion polynomial -> constant
 # We don't define a method for `Base.convert` to reduce invalidations;
 # see https://github.com/JuliaAlgebra/MultivariatePolynomials.jl/pull/172
@@ -89,11 +71,6 @@ function convert_to_constant(::Type{S}, p::_APL) where {S}
     return s
 end
 Base.convert(::Type{T}, p::_APL) where {T<:Number} = convert_to_constant(T, p)
-function convert_to_constant(p::_APL{S}) where {S}
-    return convert_to_constant(S, p)
+function convert_to_constant(p::_APL)
+    return convert_to_constant(coefficient_type(p), p)
 end
-
-# Also covers, e.g., `convert(_APL, ::P)` where `P<:_APL`
-Base.convert(::Type{PT}, p::PT) where {PT<:_APL} = p
-# Disambiguation: identity conversion for AbstractPolynomialLike
-Base.convert(::Type{PT}, p::PT) where {PT<:AbstractPolynomialLike} = p