Merge pull request #231 from JuliaSymbolics/s/rpf

WIP rational polynomial form

Author: shashi

Project.toml

@@ -1,10 +1,11 @@
 name = "SymbolicUtils"
 uuid = "d1185830-fcd6-423d-90d6-eec64667417b"
 authors = ["Shashi Gowda"]
-version = "0.13.4"
+version = "0.14.0"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
+Bijections = "e2ed5e7c-b2de-5872-ae92-c73ca462fb04"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9"

page/_layout/sidebar.html

@@ -28,6 +28,7 @@ <h1 class="sidebar-title"><a href="/">SymbolicUtils<span style="opacity:0.4">.jl
     </style>
     <nav class="sidebar-nav" style="opacity: 0.9">
       <a class="sidebar-nav-item {{ispage index.html}}active{{end}}" href="/">Getting Started</a>
+      <a class="sidebar-nav-item {{ispage representation.html}}active{{end}}" href="/representation/">Term representation</a>
       <a class="sidebar-nav-item {{ispage rewrite/*}}active{{end}}" href="/rewrite/">Term Rewriting</a>
       <a class="sidebar-nav-item {{ispage interface/*}}active{{end}}" href="/interface/">Interfacing</a>
       <a class="sidebar-nav-item {{ispage codegen/*}}active{{end}}" href="/codegen/">Code generation</a>

page/index.md

@@ -20,6 +20,7 @@ where appropriate -->
 - Fast expressions
 - A [rule-based rewriting language](/rewrite/#rule-based_rewriting).
 - A [combinator library](/rewrite/#composing_rewriters) for making rewriters.
+- [Efficient representation](/representation/) of numeric expressions
 - Type promotion:
   - Symbols (`Sym`s) carry type information. ([read more](#creating_symbolic_expressions))
   - Compound expressions composed of `Sym`s propagate type information. ([read more](#expression_interface))
@@ -131,7 +132,7 @@ g(2//5, g(1, β))
 
 ## Expression interface
 
-Symbolic expressions are of type `Term{T}`, `Add{T}`, `Mul{T}` or `Pow{T}` and denote some function call where one or more arguments are themselves such expressions or `Sym`s.
+Symbolic expressions are of type `Term{T}`, `Add{T}`, `Mul{T}`, `Pow{T}` or `Div{T}` and denote some function call where one or more arguments are themselves such expressions or `Sym`s. See more about the representation [here](/representation/).
 
 All the expression types support the following:
 
@@ -155,8 +156,10 @@ The rules with which the canonical form of `Symbolic{<:Number}` terms are constr
 
  -  `0 + x`, `1 * x` and `x^1` always gives `x`
  -  `0 * x` always gives `0` and `x ^ 0` gives `1`
- -  `-x`, `1/x` and `x\1` get transformed into `(-1)*x`, `x^(-1)` and `x^(-1)`.
+ -  `-x` becomes `(-1)*x`
  -  commutativity and associativity over `+` and `*` are assumed. Re-ordering of terms will be done under a [total order](https://github.com/JuliaSymbolics/SymbolicUtils.jl/blob/master/src/ordering.jl)
+ - `p/q * x` or `x * p/q` results in `(p*x)/q`
+ - `p/q * x/y` results in `(p*x)/(q*y)`
  -  `x + ... + x` will be fused into `n*x` with type `Mul`
  -  `x * ... * x` will be fused into `x^n` with type `Pow`
  -  sum of `Add`'s are fused

page/representation.md

@@ -0,0 +1,86 @@
+## Term representation and simplification
+
+### Preliminary representation
+
+Expressions are stored as a tree where the nodes store their children in some efficient data structure. A node `n` is a non-leaf node `istree(n)` is true. For such nodes, `gethead` and `getargs` must be defined. These return the head and args of the term represented by the subtree `n`.
+
+A generic term is represented with the `Term` type. It simply holds the function `f` being called and a vector of arguments in the order they were passed to `f`. The operators `+` (and `-`), `*`, `/` and `^` create terms with special storage.
+
+Linear combinations such as $\alpha_1  x_1 + \alpha_2 x_2 +...+ \alpha_n x_n$
+are represented by `Add(Dict(x_1 => \alpha_1, x_2 => \alpha_2, ..., x_n => \alpha_n))`. Now $x_n$ may themselves be other types mentioned in this section, but may not be an `Add`. When an `Add` is added to an `Add`, we merge their dictionaries and adding up matching coefficients.
+
+Similarly, $x_1^{m_1}x_2^{m_2}\ellipsisx_{m_n}$ is represented by
+`Mul(Dict(x_1 => m_1, x_2 => m_2,..., x_n => m_n))`. $x_i$ may not themselves be `Mul`, multiplying a Mul with another Mul returns a flattened Mul.
+
+$p / q$ is represented by `Div(p, q)`. The result of `*` on `Div` is maintainted as a `Div`. For example, `Div(p_1, q_1) * Div(p_2, q_2)` results in `Div(p_1 * p_2, q_1 * q_2)` and so on. The effect is, in `Div(p, q)`, `p` or `q` or, if they are Mul, any of their multiplicands is not a Div. So `Mul`s must always be nested inside a `Div` and can never show up immediately wrapping it.
+
+
+### Polynomial representation
+
+Packages like DynamicPolynomials.jl provide even more efficient representations for polynomials. They also have efficient operations such as multi-variate polynomial GCD. However, DynamicPolynomials can only represent flat polynomials. For example, `(x-3)*(x+5)` can only be represented as `(x^2) + 15 - 8x`. Moreover, DynamicPolynomials does not have ways to represent generic Terms such as `sin(x-y)` in the tree.
+
+To reconcile these differences while being able to use the efficient representation of DynamicPolynomials we have the `PolyForm` type. This type holds a polynomial and the mappings necessary to present them as a SymbolicUtils expression.  The mappings constructed for the conversion are 1) a bijection from DynamicPolynomials PolyVar type to a Symbolics `Sym`, and 2) a mapping from `Sym`s to non-polynomial terms that the `Sym`s stand-in for. These terms may themselves contain PolyForm if there are polynomials inside them. The mappings are transiently global, that is, when all references to the mappings go out of scope, they are released and re-created.
+
+```julia
+julia> @syms x y
+(x, y)
+
+julia> PolyForm((x-3)*(y-5))
+x*y + 15 - 5x - 3y
+```
+
+Further, generic terms such as `sin` or `cos` terms are replaced with a symbol when representing the polynomial.
+
+```julia
+julia> p = PolyForm((sin(x) + cos(x))^2)
+(cos(x)^2) + 2cos(x)*sin(x) + (sin(x)^2)
+
+julia> p.p # this is the actual DynamicPolynomial stored
+cos_3658410937268741549² + 2cos_3658410937268741549sin_10720964503106793468 + sin_10720964503106793468²
+```
+
+By default, polynomials inside generic terms are not also converted to PolyForm. For example,
+
+```julia
+julia> PolyForm(sin((x-3)*(y-5)))
+sin((x - 3)*(y - 5))
+```
+But you can pass in the `recurse=true` keyword argument to get this going.
+
+```julia
+julia> PolyForm(sin((x-3)*(y-5)), recurse=true)
+sin(x*y + 15 - 5x - 3y)
+```
+
+Polynomials are constructed by first turning symbols and generic terms into Polynomial variables and then applying the `+`, `*`, `^` operations on these variables. You can control which of these operations are actually constructed with the `Fs` keyword argument. It is a Union type of the functions to apply. For example, let's say you want to turn terms into polynomials by only applying the `+` and `^` operations, and want to preserve `*` operations as-is, you could pass in `Fs=Union{typeof(+), typeof(^)}`
+
+```julia
+julia> PolyForm((x+y)^2*(x-y), Fs=Union{typeof(+), typeof(^)}, recurse=true)
+((x - (y))*((x^2) + 2x*y + (y^2)))
+```
+
+Note that in this case `recurse=true` was necessary as otherwise the polynomialization would not descend into the `*` operation as it is now considered a generic operation.
+
+### Simplifying fractions
+
+`simplify_fractions(expr)` recurses through `expr` finding `Div`s and simplifying them using polynomial divison.
+
+First the factors of the numerators and the denominators are converted into PolyForm objects, then numerators and denominators are divided by their respective pairwise GCDs. The conversion of the numerator and denominator into PolyForm is set up so that `simplify_fractions` does not result in increase in the expression size due to polynomial expansion. Specifically, the factors are individually converted into PolyForm objects, and any powers of polynomial is not expanded, but the divison process repeatedly divides them as many times as the power.
+
+
+```julia
+julia> simplify_fractions((x*y+5x+3y+15)/((x+3)*(x-4)))
+(5.0 + y) / (x - 4)
+
+julia> simplify_fractions((x*y+5x+3y+15)^2/((x+3)*(x-4)))
+((5.0 + y)*(15 + 5x + x*y + 3y)) / (x - 4)
+
+julia> simplify_fractions(3/(x+3) + x/(x+3))
+1
+
+julia> simplify_fractions(sin(x)/cos(x) + cos(x)/sin(x))
+(cos(x)^2 + sin(x)^2) / (cos(x)*sin(x))
+
+julia> simplify(ans)
+1 / (cos(x)*sin(x))
+```

src/SymbolicUtils.jl

@@ -45,6 +45,7 @@ const MP = MultivariatePolynomials
 import DynamicPolynomials
 export expand
 include("abstractalgebra.jl")
+include("polyform.jl")
 
 # Term ordering
 include("ordering.jl")

src/abstractalgebra.jl

@@ -1,181 +0,0 @@
-# Polynomial Normal Form
-
-"""
-    labels!(dict, t)
-
-Find all terms that are not + and * and replace them
-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,  variable_type::Type)
-    sym2term, term2sym = dicts
-    if !haskey(term2sym, t)
-        sym2term[t] = t
-        term2sym[t] = t
-    end
-    return t
-end
-
-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, 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, variable_type), arguments(t)), symtype(t))
-    else
-        sym2term, term2sym = dicts
-        if haskey(term2sym, t)
-            return term2sym[t]
-        end
-        if istree(t)
-            tt = arguments(t)
-            sym = Sym{symtype(t)}(gensym(nameof(operation(t))))
-            dicts2 = _dicts(dicts[2])
-            sym2term[sym] = similarterm(t, operation(t),
-                                        map(x->to_mpoly(x, variable_type, dicts)[1], arguments(t)),
-                                        symtype(t))
-        else
-            sym = Sym{symtype(t)}(gensym("literal"))
-            sym2term[sym] = t
-        end
-
-        term2sym[t] = sym
-
-        return sym
-    end
-end
-
-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)
-
-let
-    mpoly_preprocess = [@rule(identity(~x) => ~x)
-                        @rule(zero(~x) => 0)
-                        @rule(one(~x) => 1)]
-
-    simterm(x, f, args;metadata=nothing) = similarterm(x,f,args, symtype(x); metadata=metadata)
-    mpoly_rules = [@rule(~x::ismpoly - ~y::ismpoly => ~x + -1 * (~y))
-                   @rule(-(~x) => -1 * ~x)
-                   @acrule(~x::ismpoly + ~y::ismpoly => ~x + ~y)
-                   @rule(+(~x) => ~x)
-                   @acrule(~x::ismpoly * ~y::ismpoly => ~x * ~y)
-                   @rule(*(~x) => ~x)
-                   @rule((~x::ismpoly)^(~a::isnonnegint) => (~x)^(~a))]
-    global const MPOLY_CLEANUP = Fixpoint(Postwalk(PassThrough(RestartedChain(mpoly_preprocess)), similarterm=simterm))
-    MPOLY_MAKER = Fixpoint(Postwalk(PassThrough(RestartedChain(mpoly_rules)), similarterm=simterm))
-
-    global to_mpoly
-    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, variable_type)
-
-        if isempty(sym2term)
-            return MPOLY_MAKER(labeled), Dict{Sym,Any}()
-        end
-
-        ks = sort(collect(keys(sym2term)), lt=<ₑ)
-        vars = MP.similarvariable.(variable_type, nameof.(ks))
-
-        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(string.(nameof.(keys(dict))), keys(dict)))
-    dict = copy(dict)
-    for (k, v) in dict
-        dict[k] = _to_term(reference, v, dict, syms)
-    end
-    return _to_term(reference, x, dict, syms)
-    #return substitute(t, dict, fold=false)
-end
-
-_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
-
-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
-        return MP.coefficient(term) * mono
-    end
-end
-
-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))
-    else
-        if haskey(dict, x)
-            return dict[x]
-        else
-            return x
-        end
-    end
-end
-
-<ₑ(a::MP.AbstractPolynomialLike, b::MP.AbstractPolynomialLike) = false
-
-"""
-    expand(expr, variable_type::Type=DynamicPolynomials.PolyVar{true})
-
-Expand expressions by distributing multiplication over addition, e.g.,
-`a*(b+c)` becomes `ab+ac`.
-
-`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(expr, variable_type::Type=DynamicPolynomials.PolyVar{true})
-    to_term(expr, to_mpoly(expr, variable_type)...)
-end
-
-## Hack to fix https://github.com/JuliaAlgebra/MultivariatePolynomials.jl/issues/169
-
-Base.promote_rule(::Type{S}, ::Type{T}) where {S<:Symbolic, T<:MP.AbstractPolynomialLike}= Any
-Base.promote_rule(::Type{T}, ::Type{S}) where {S<:Symbolic, T<:MP.AbstractPolynomialLike}= Any
-
-Base.@deprecate polynormalize(x) expand(x)

src/methods.jl

@@ -80,6 +80,7 @@ end
 @number_methods(Add, term(f, a), term(f, a, b), skipbasics)
 @number_methods(Mul, term(f, a), term(f, a, b), skipbasics)
 @number_methods(Pow, term(f, a), term(f, a, b), skipbasics)
+@number_methods(Div, term(f, a), term(f, a, b), skipbasics)
 
 for f in diadic
     @eval promote_symtype(::$(typeof(f)),