Merge pull request #271 from blegat/multivariate_polynomials

Allow any implementation of the MultivariatePolynomials API

Author: shashi

Project.toml

@@ -1,18 +1,19 @@
 name = "SymbolicUtils"
 uuid = "d1185830-fcd6-423d-90d6-eec64667417b"
 authors = ["Shashi Gowda"]
-version = "0.12.1"
+version = "0.13.0"
 
 [deps]
-AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 NaNMath = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -21,14 +22,15 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
-AbstractAlgebra = "0.9, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18"
 AbstractTrees = "0.3"
 ChainRulesCore = "0.9, 0.10"
 Combinatorics = "1.0"
 ConstructionBase = "1.1"
 DataStructures = "0.18"
+DynamicPolynomials = "0.3"
 IfElse = "0.1"
 LabelledArrays = "1.5"
+MultivariatePolynomials = "0.3"
 NaNMath = "0.3"
 Setfield = "0.7"
 SpecialFunctions = "0.10, 1.0"

benchmark/goldstein_price.jl

@@ -0,0 +1,19 @@
+using SymbolicUtils, Test
+using SymbolicUtils: Term, symtype
+using BenchmarkTools
+
+function goldstein_price()
+    @syms x1 x2
+    f1 = x1+x2+1
+    f2 = 19-14*x1+3*x1^2-14*x2+6*x1*x2+3*x2^2
+    f3 = 2*x1-3*x2
+    f4 = 18-32*x1+12*x1^2+48*x2-36*x1*x2+27*x2^2
+
+    # f(x) is the Goldstein-Price function
+    f = (1+f1^2*f2)*(30+f3^2*f4)
+    @btime expand($f)
+    f2 = f^2
+    @btime expand($f2)
+    f3 = f^3
+    @btime expand($f3)
+end

src/SymbolicUtils.jl

@@ -36,8 +36,9 @@ include("rule.jl")
 include("matchers.jl")
 
 # Convert to an efficient multi-variate polynomial representation
-import AbstractAlgebra.Generic: MPoly, PolynomialRing, ZZ, exponent_vector
-using AbstractAlgebra: ismonomial, symbols
+import MultivariatePolynomials
+const MP = MultivariatePolynomials
+import DynamicPolynomials
 export expand
 include("abstractalgebra.jl")
 

src/abstractalgebra.jl

@@ -9,7 +9,7 @@ with a symbol, store the symbol => term mapping in `dict`.
 function labels! end
 
 # Turn a Term into a multivariate polynomial
-function labels!(dicts, t::Sym)
+function labels!(dicts, t::Sym,  variable_type::Type)
     sym2term, term2sym = dicts
     if !haskey(term2sym, t)
         sym2term[t] = t
@@ -18,14 +18,14 @@ function labels!(dicts, t::Sym)
     return t
 end
 
-function labels!(dicts, t)
-    if t isa Integer
+function labels!(dicts, t, variable_type::Type)
+    if t isa Number
         return t
     elseif istree(t) && (operation(t) == (*) || operation(t) == (+) || operation(t) == (-))
         tt = arguments(t)
-        return similarterm(t, operation(t), map(x->labels!(dicts, x), tt), symtype(t))
+        return similarterm(t, operation(t), map(x->labels!(dicts, x, variable_type), tt), symtype(t))
     elseif istree(t) && operation(t) == (^) && length(arguments(t)) > 1 && isnonnegint(arguments(t)[2])
-        return similarterm(t, operation(t), map(x->labels!(dicts, x), arguments(t)), symtype(t))
+        return similarterm(t, operation(t), map(x->labels!(dicts, x, variable_type), arguments(t)), symtype(t))
     else
         sym2term, term2sym = dicts
         if haskey(term2sym, t)
@@ -36,7 +36,7 @@ function labels!(dicts, t)
             sym = Sym{symtype(t)}(gensym(nameof(operation(t))))
             dicts2 = _dicts(dicts[2])
             sym2term[sym] = similarterm(t, operation(t),
-                                        map(x->to_mpoly(x, dicts)[1], arguments(t)),
+                                        map(x->to_mpoly(x, variable_type, dicts)[1], arguments(t)),
                                         symtype(t))
         else
             sym = Sym{symtype(t)}(gensym("literal"))
@@ -49,7 +49,7 @@ function labels!(dicts, t)
     end
 end
 
-ismpoly(x) = x isa MPoly || x isa Integer
+ismpoly(x) = x isa MP.AbstractPolynomialLike || x isa Number
 isnonnegint(x) = x isa Integer && x >= 0
 
 _dicts(t2s=OrderedDict{Any, Sym}()) = (OrderedDict{Sym, Any}(), t2s)
@@ -71,70 +71,82 @@ let
     MPOLY_MAKER = Fixpoint(Postwalk(PassThrough(RestartedChain(mpoly_rules)), similarterm=simterm))
 
     global to_mpoly
-    function to_mpoly(t, dicts=_dicts())
+    function to_mpoly(t, variable_type::Type=DynamicPolynomials.PolyVar{true}, dicts=_dicts())
         # term2sym is only used to assign the same
         # symbol for the same term -- in other words,
         # it does common subexpression elimination
         t = MPOLY_CLEANUP(t)
         sym2term, term2sym = dicts
-        labeled = labels!((sym2term, term2sym), t)
+        labeled = labels!((sym2term, term2sym), t, variable_type)
 
         if isempty(sym2term)
             return MPOLY_MAKER(labeled), Dict{Sym,Any}()
         end
 
         ks = sort(collect(keys(sym2term)), lt=<ₑ)
-        R, vars = PolynomialRing(ZZ, String.(nameof.(ks)))
+        vars = MP.similarvariable.(variable_type, nameof.(ks))
 
-        replace_with_poly = Dict{Sym,MPoly}(zip(ks, vars))
+        replace_with_poly = Dict{Sym,eltype(vars)}(zip(ks, vars))
         t_poly = substitute(labeled, replace_with_poly, fold=false)
         MPOLY_MAKER(t_poly), sym2term
     end
 end
 
 function to_term(reference, x, dict)
-    syms = Dict(zip(nameof.(keys(dict)), keys(dict)))
+    syms = Dict(zip(string.(nameof.(keys(dict))), keys(dict)))
     dict = copy(dict)
     for (k, v) in dict
         dict[k] = _to_term(reference, v, dict, syms)
     end
-    _to_term(reference, x, dict, syms)
+    return _to_term(reference, x, dict, syms)
+    #return substitute(t, dict, fold=false)
 end
 
-function _to_term(reference, x::MPoly, dict, syms)
-
-    function mul_coeffs(exps, ring)
-        l = length(syms)
-        ss = symbols(ring)
-        monics = [e == 1 ? syms[ss[i]] : syms[ss[i]]^e for (i, e) in enumerate(exps) if !iszero(e)]
-        if length(monics) == 1
-            return monics[1]
-        elseif length(monics) == 0
-            return 1
-        else
-            return similarterm(reference, *, monics, symtype(reference))
+_to_term(reference, x::Number, dict, syms) = x
+_to_term(reference, var::MP.AbstractVariable, dict, syms) = substitute(syms[MP.name(var)], dict, fold=false)
+function _to_term(reference, mono::MP.AbstractMonomialLike, dict, syms)
+    monics = [
+        begin
+            t = _to_term(reference, var, dict, syms)
+            exp == 1 ? t : t^exp
         end
+        for (var, exp) in MP.powers(mono) if !iszero(exp)
+    ]
+    if length(monics) == 1
+        return monics[1]
+    elseif isempty(monics)
+        return 1
+    else
+        return similarterm(reference, *, monics, symtype(reference))
     end
+end
 
-    monoms = [mul_coeffs(exponent_vector(x, i), x.parent) for i in 1:x.length]
-    if length(monoms) == 0
-        return 0
-    elseif length(monoms) == 1
-        t = !isone(x.coeffs[1]) ?  monoms[1] * Int(x.coeffs[1]) : monoms[1]
+function _to_term(reference, term::MP.AbstractTermLike, dict, syms)
+    coef = MP.coefficient(term)
+    mono = _to_term(reference, MP.monomial(term), dict, syms)
+    if isone(coef)
+        return mono
     else
-        t = similarterm(reference,
-                        +,
-                        map((x,y)->isone(y) ? x : Int(y)*x,
-                            monoms, x.coeffs[1:length(monoms)]),
-                        symtype(reference))
+        return MP.coefficient(term) * mono
     end
+end
 
-    substitute(t, dict, fold=false)
+function _to_term(reference, x::MP.AbstractPolynomialLike, dict, syms)
+    if MP.nterms(x) == 0
+        return 0
+    elseif MP.nterms(x) == 1
+        return _to_term(reference, first(MP.terms(x)), dict, syms)
+    else
+        terms = map(MP.terms(x)) do term
+            _to_term(reference, term, dict, syms)
+        end
+        return similarterm(reference, +, terms, symtype(reference))
+    end
 end
 
 function _to_term(reference, x, dict, vars)
     if istree(x)
-        t=similarterm(x, operation(x), _to_term.((reference,), arguments(x), (dict,), (vars,)), symtype(x))
+        t = similarterm(x, operation(x), _to_term.((reference,), arguments(x), (dict,), (vars,)), symtype(x))
     else
         if haskey(dict, x)
             return dict[x]
@@ -144,18 +156,21 @@ function _to_term(reference, x, dict, vars)
     end
 end
 
-<ₑ(a::MPoly, b::MPoly) = false
+<ₑ(a::MP.AbstractPolynomialLike, b::MP.AbstractPolynomialLike) = false
 
 """
-    expand(expr)
+    expand(expr, variable_type::Type=DynamicPolynomials.PolyVar{true})
 
-Expand expressions by distributing multiplication over addition.
+Expand expressions by distributing multiplication over addition, e.g.,
+`a*(b+c)` becomes `ab+ac`.
 
-`a*(b+c)` becomes `ab+ac`. `expand` uses [AbstractAlgebra.jl](https://nemocas.github.io/AbstractAlgebra.jl/latest/) to construct
-dense Multi-variate polynomial to do this very fast.
+`expand` uses replace symbols and non-algebraic expressions by variables of type
+`variable_type` to compute the distribution using a specialized sparse
+multivariate polynomials implementation.
+`variable_type` can be any subtype of `MultivariatePolynomials.AbstractVariable`.
 """
-function expand(x)
-    to_term(x, to_mpoly(x)...)
+function expand(expr, variable_type::Type=DynamicPolynomials.PolyVar{true})
+    to_term(expr, to_mpoly(expr, variable_type)...)
 end
 
 Base.@deprecate polynormalize(x) expand(x)

src/types.jl

@@ -494,7 +494,7 @@ end
 function show_mul(io, args)
     length(args) == 1 && return print_arg(io, *, args[1])
 
-    paren_scalar = args[1] isa Complex || args[1] isa Rational
+    paren_scalar = args[1] isa Complex || args[1] isa Rational || (args[1] isa Number && !isfinite(args[1]))
     minus = args[1] isa Number && args[1] == -1
     unit = args[1] isa Number && args[1] == 1
     nostar = !paren_scalar && args[1] isa Number && !(args[2] isa Number)
@@ -867,7 +867,10 @@ istree(a::Pow) = true
 
 operation(a::Pow) = ^
 
-arguments(a::Pow) = [a.base, a.exp]
+# Use `Union` to avoid promoting the base and exponent to the same type.
+# For instance, if `a.base` is a multivariate polynomial and  `a.exp` is a number,
+# we don't want to promote `a.exp` to a multivariate polynomial.
+arguments(a::Pow) = Union{typeof(a.base), typeof(a.exp)}[a.base, a.exp]
 
 Base.hash(p::Pow, u::UInt) = hash(p.exp, hash(p.base, u))