Complex-valued variables (#121)

* First step to implement complex variables functionality

* Move monomial function to correct file

* Rename ordvar to ordinary_variable

* Specialized complex routines for MonomialVector

* fix

* coefficienttype -> coefficient_type

* Remove incorrect method

* Add a (sometimes working) implementation for incomplete substitution

* Fix old syntax

* Format

* Remove TODO

* Add complex kind parameter to VT("name") constructor

* Remove iscomplex

* Remove dev dependency of MP, 0.5.3 has landed

* Add ComplexKind docstring

* Rename ComplexKind instances

* typo

* change substitution rules

* Rename polycvar macro

* Error for partial substitutions

* Format

* Complex substitution test cases

* More strict substitution rules: only allow ordinary variable

---------

Co-authored-by: Benoît Legat <benoit.legat@gmail.com>

Project.toml

@@ -13,7 +13,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-MultivariatePolynomials = "0.5"
+MultivariatePolynomials = "0.5.3"
 MutableArithmetics = "1"
 Reexport = "1"
 julia = "1"

src/DynamicPolynomials.jl

@@ -72,11 +72,11 @@ function MP.similar_variable(
 ) where {V,M,S}
     return MP.similar_variable(P, S)
 end
-function MP.similar_variable(
-    ::Union{PolyType{V,M},Type{<:PolyType{V,M}}},
-    s::Symbol,
-) where {V,M}
-    return Variable(string(s), V, M)
+function MP.similar_variable(p::PolyType{V,M}, s::Symbol) where {V,M}
+    return Variable(string(s), V, M, isreal(p) ? REAL : COMPLEX)
+end
+function MP.similar_variable(::Type{<:PolyType{V,M}}, s::Symbol) where {V,M}
+    return Variable(string(s), V, M, REAL) # we cannot infer this from the type
 end
 
 include("promote.jl")

src/comp.jl

@@ -27,14 +27,19 @@ function (==)(
     x::Variable{<:AnyCommutative{CreationOrder}},
     y::Variable{<:AnyCommutative{CreationOrder}},
 )
-    return x.variable_order.order.id == y.variable_order.order.id
+    return x.variable_order.order.id == y.variable_order.order.id &&
+           x.kind == y.kind
 end
 
 function Base.isless(
     x::Variable{<:AnyCommutative{CreationOrder}},
     y::Variable{<:AnyCommutative{CreationOrder}},
 )
-    return isless(y.variable_order.order.id, x.variable_order.order.id)
+    if x.variable_order.order.id == y.variable_order.order.id
+        return isless(y.kind, x.kind)
+    else
+        return isless(y.variable_order.order.id, x.variable_order.order.id)
+    end
 end
 
 # Comparison of Monomial

src/mono.jl

@@ -113,3 +113,10 @@ function __add_variables!(
     mono.z[map] = tmp
     return
 end
+
+# for efficiency reasons
+function Base.conj(x::Monomial{V,M}) where {V<:Commutative,M}
+    cv = conj.(x.vars)
+    perm = sortperm(cv, rev = true)
+    return Monomial{V,M}(cv[perm], x.z[perm])
+end

src/monomial_vector.jl

@@ -88,6 +88,28 @@ end
 function MultivariatePolynomials.maxdegree(x::MonomialVector)
     return isempty(x) ? 0 : maximum(sum.(x.Z))
 end
+# Complex-valued degrees for monomial vectors
+for (fun, call, def, ret) in [
+    (:extdegree_complex, :extrema, (0, 0), :((min(v1, v2), max(v1, v2)))),
+    (:mindegree_complex, :minimum, 0, :(min(v1, v2))),
+    (:maxdegree_complex, :maximum, 0, :(max(v1, v2))),
+]
+    eval(quote
+        function MultivariatePolynomials.$fun(x::MonomialVector)
+            isempty(x) && return $def
+            vars = variables(x)
+            @assert(!any(isrealpart, vars) && !any(isimagpart, vars))
+            grouping = isconj.(vars)
+            v1 = $call(sum(z[grouping]) for z in x.Z)
+            grouping = map(!, grouping)
+            v2 = $call(sum(z[grouping]) for z in x.Z)
+            return $ret
+        end
+    end)
+end
+# faster complex-related functions
+Base.isreal(x::MonomialVector) = all(isreal, x.vars)
+Base.conj(x::MonomialVector) = MonomialVector(conj.(x.vars), x.Z)
 
 MP.variables(m::Union{Monomial,MonomialVector}) = m.vars
 
@@ -118,7 +140,11 @@ end
 # TODO replace by MP function
 function _error_for_negative_degree(deg)
     if deg < 0
-        throw(ArgumentError("The degree should be a nonnegative number but the provided degree `$deg` is negative."))
+        throw(
+            ArgumentError(
+                "The degree should be a nonnegative number but the provided degree `$deg` is negative.",
+            ),
+        )
     end
 end
 

src/subs.jl

@@ -9,23 +9,46 @@ end
 
 function fillmap!(
     vals,
-    vars::Vector{<:Variable{<:NonCommutative}},
+    vars::Vector{<:Variable{C}},
     s::MP.Substitution,
-)
-    for j in eachindex(vars)
-        if vars[j] == s.first
-            vals[j] = s.second
+) where {C}
+    # We may assign a complex or real variable to its value (ordinary substitution).
+    # We follow the following rules:
+    # - If a single substitution rule determines the value of a real variable, just substitute it.
+    # - If a single substitution rule determines the value of a complex variable or its conjugate, substitute the appropriate
+    #   value whereever something related to this variable is found (i.e., the complex variable, the conjugate variable, or
+    #   its real or imaginary part)
+    # - If a single substitution rule determines the value of the real or imaginary part of a complex variable alone, then only
+    #   replace the real or imaginary parts if they occur explicitly. Don't do a partial substitution, i.e., `z` with the rule
+    #   `zᵣ => 1` is left alone and not changed into `1 + im*zᵢ`. Even if both the real and imaginary parts are substituted as
+    #   two individual rules (which we don't know of in this method), `z` will not be replaced.
+    if s.first.kind == REAL
+        for j in eachindex(vars)
+            if vars[j] == s.first
+                vals[j] = s.second
+                C == Commutative && break
+            end
+        end
+    else
+        s.first.kind == COMPLEX || throw(
+            ArgumentError(
+                "Substitution with complex variables requires the ordinary_variable in the substitution specification",
+            ),
+        )
+        for j in eachindex(vars)
+            if vars[j].variable_order.order.id ==
+               s.first.variable_order.order.id
+                if vars[j].kind == COMPLEX
+                    vals[j] = s.second
+                elseif vars[j].kind == CONJ
+                    vals[j] = conj(s.second)
+                elseif vars[j].kind == REAL_PART
+                    vals[j] = real(s.second)
+                else
+                    vals[j] = imag(s.second)
+                end
+            end
         end
-    end
-end
-function fillmap!(
-    vals,
-    vars::Vector{<:Variable{<:Commutative}},
-    s::MP.Substitution,
-)
-    j = findfirst(isequal(s.first), vars)
-    if j !== nothing
-        vals[j] = s.second
     end
 end
 function fillmap!(vals, vars, s::MP.AbstractMultiSubstitution)

src/var.jl

@@ -1,20 +1,34 @@
-export Variable, @polyvar, @ncpolyvar
+export Variable, @polyvar, @ncpolyvar, @complex_polyvar
 
-function polyarrayvar(variable_order, monomial_order, prefix, indices...)
+function polyarrayvar(
+    variable_order,
+    monomial_order,
+    complex_kind,
+    prefix,
+    indices...,
+)
     return map(
         i -> Variable(
             "$(prefix)[$(join(i, ","))]",
             variable_order,
             monomial_order,
+            complex_kind,
         ),
         Iterators.product(indices...),
     )
 end
 
-function buildpolyvar(var, variable_order, monomial_order)
+function buildpolyvar(var, variable_order, monomial_order, complex_kind)
     if isa(var, Symbol)
         var,
-        :($(esc(var)) = $Variable($"$var", $variable_order, $monomial_order))
+        :(
+            $(esc(var)) = $Variable(
+                $"$var",
+                $variable_order,
+                $monomial_order,
+                $complex_kind,
+            )
+        )
     else
         isa(var, Expr) || error("Expected $var to be a variable name")
         Base.Meta.isexpr(var, :ref) ||
@@ -28,26 +42,28 @@ function buildpolyvar(var, variable_order, monomial_order)
             $(esc(varname)) = polyarrayvar(
                 $(variable_order),
                 $(monomial_order),
+                $(complex_kind),
                 $prefix,
                 $(esc.(var.args[2:end])...),
             )
         )
     end
 end
 
-function buildpolyvars(args, variable_order, monomial_order)
+function buildpolyvars(args, variable_order, monomial_order, complex_kind)
     vars = Symbol[]
     exprs = []
     for arg in args
-        var, expr = buildpolyvar(arg, variable_order, monomial_order)
+        var, expr =
+            buildpolyvar(arg, variable_order, monomial_order, complex_kind)
         push!(vars, var)
         push!(exprs, expr)
     end
     return vars, exprs
 end
 
 # Inspired from `JuMP.Containers._extract_kw_args`
-function _extract_kw_args(args, variable_order)
+function _extract_kw_args(args, variable_order, complex_kind)
     positional_args = Any[]
     monomial_order = :($(MP.Graded{MP.LexOrder}))
     for arg in args
@@ -56,29 +72,42 @@ function _extract_kw_args(args, variable_order)
                 variable_order = esc(arg.args[2])
             elseif arg.args[1] == :monomial_order
                 monomial_order = esc(arg.args[2])
+            elseif arg.args[1] == :complex
+                complex_kind = arg.args[2] == :true ? COMPLEX : REAL
             else
                 error("Unrecognized keyword argument `$(arg.args[1])`")
             end
         else
             push!(positional_args, arg)
         end
     end
-    return positional_args, variable_order, monomial_order
+    return positional_args, variable_order, monomial_order, complex_kind
 end
 
 # Variable vector x returned garanteed to be sorted so that if p is built with x then vars(p) == x
 macro polyvar(args...)
-    pos_args, variable_order, monomial_order =
-        _extract_kw_args(args, :($(Commutative{CreationOrder})))
-    vars, exprs = buildpolyvars(pos_args, variable_order, monomial_order)
+    pos_args, variable_order, monomial_order, complex_kind =
+        _extract_kw_args(args, :($(Commutative{CreationOrder})), REAL)
+    vars, exprs =
+        buildpolyvars(pos_args, variable_order, monomial_order, complex_kind)
     return :($(foldl((x, y) -> :($x; $y), exprs, init = :()));
     $(Expr(:tuple, esc.(vars)...)))
 end
 
 macro ncpolyvar(args...)
-    pos_args, variable_order, monomial_order =
-        _extract_kw_args(args, :($(NonCommutative{CreationOrder})))
-    vars, exprs = buildpolyvars(pos_args, variable_order, monomial_order)
+    pos_args, variable_order, monomial_order, complex_kind =
+        _extract_kw_args(args, :($(NonCommutative{CreationOrder})), REAL)
+    vars, exprs =
+        buildpolyvars(pos_args, variable_order, monomial_order, complex_kind)
+    return :($(foldl((x, y) -> :($x; $y), exprs, init = :()));
+    $(Expr(:tuple, esc.(vars)...)))
+end
+
+macro complex_polyvar(args...)
+    pos_args, variable_order, monomial_order, complex_kind =
+        _extract_kw_args(args, :($(Commutative{CreationOrder})), COMPLEX)
+    vars, exprs =
+        buildpolyvars(pos_args, variable_order, monomial_order, complex_kind)
     return :($(foldl((x, y) -> :($x; $y), exprs, init = :()));
     $(Expr(:tuple, esc.(vars)...)))
 end
@@ -114,19 +143,54 @@ function instantiate(::Type{NonCommutative{O}}) where {O}
     return NonCommutative(instantiate(O))
 end
 
+"""
+    ComplexKind
+
+The type used to identify whether a variable is complex-valued. It has the following instances:
+- `REAL`: the variable is real-valued, [`conj`](@ref) and [`real`](@ref) are identities, [`imag`](@ref) is zero.
+- `COMPLEX`: the variable is a declared complex-valued variable with distinct conjugate, real, and imaginary part.
+  [`ordinary_variable`](@ref) is an identity.
+- `CONJ`: the variable is the conjugate of a declared complex-valued variable. [`ordinary_variable`](@ref) is equal to
+  [`conj`](@ref). It will print with an over-bar.
+- `REAL_PART`: the variable is the real part of a complex-valued variable. It is real-valued itself: [`conj`](@ref) and
+  [`real`](@ref) are identities, [`imag`](@ref) is zero. [`ordinary_variable`](@ref) will give back the original complex-valued
+  _declared_ variable from which the real part was extracted. It will print with an `ᵣ` subscript.
+- `IMAG_PART`: the variable is the imaginary part of a declared complex-valued variable (or the negative imaginary part of the
+  conjugate of a declared complex-valued variable). It is real-valued itself: [`conj`](@ref) and [`real`](@ref) are identities,
+  [`imag`](@ref) is zero. [`ordinary_variable`](@ref) will give back the original complex-valued _declared_ variable from which
+  the imaginary part was extracted. It will print with an `ᵢ` subscript.
+"""
+@enum ComplexKind REAL COMPLEX CONJ REAL_PART IMAG_PART
+
 struct Variable{V,M} <: AbstractVariable
     name::String
+    kind::ComplexKind
     variable_order::V
 
-    function Variable{V,M}(name::AbstractString) where {V<:AbstractVariableOrdering,M<:MP.AbstractMonomialOrdering}
-        return new{V,M}(convert(String, name), instantiate(V))
+    function Variable{V,M}(
+        name::AbstractString,
+        kind::ComplexKind = REAL,
+    ) where {V<:AbstractVariableOrdering,M<:MP.AbstractMonomialOrdering}
+        return new{V,M}(convert(String, name), kind, instantiate(V))
     end
+
     function Variable(
         name::AbstractString,
         ::Type{V},
         ::Type{M},
+        kind::ComplexKind = REAL,
     ) where {V<:AbstractVariableOrdering,M<:MP.AbstractMonomialOrdering}
-        return new{V,M}(convert(String, name), instantiate(V))
+        return new{V,M}(convert(String, name), kind, instantiate(V))
+    end
+
+    function Variable(
+        from::Variable{V,M},
+        kind::ComplexKind,
+    ) where {V<:AbstractVariableOrdering,M<:MP.AbstractMonomialOrdering}
+        (from.kind == REAL && kind == REAL) ||
+            (from.kind != REAL && kind != REAL) ||
+            error("Cannot change the complex type of an existing variable")
+        return new{V,M}(from.name, kind, from.variable_order)
     end
 end
 
@@ -150,6 +214,41 @@ MP.ordering(::Type{Variable{V,M}}) where {V,M} = M
 
 iscomm(::Type{Variable{C}}) where {C} = C
 
+Base.isreal(x::Variable{C}) where {C} = x.kind != COMPLEX && x.kind != CONJ
+MP.isrealpart(x::Variable{C}) where {C} = x.kind == REAL_PART
+MP.isimagpart(x::Variable{C}) where {C} = x.kind == IMAG_PART
+MP.isconj(x::Variable{C}) where {C} = x.kind == CONJ
+function MP.ordinary_variable(x::Variable)
+    return x.kind == REAL || x.kind == COMPLEX ? x : Variable(x, COMPLEX)
+end
+
+function Base.conj(x::Variable)
+    if x.kind == COMPLEX
+        return Variable(x, CONJ)
+    elseif x.kind == CONJ
+        return Variable(x, COMPLEX)
+    else
+        return x
+    end
+end
+
+function Base.real(x::Variable)
+    if x.kind == COMPLEX || x.kind == CONJ
+        return Variable(x, REAL_PART)
+    else
+        return x
+    end
+end
+function Base.imag(x::Variable{V,M}) where {V,M}
+    if x.kind == COMPLEX
+        return Variable(x, IMAG_PART)
+    elseif x.kind == CONJ
+        return _Term{V,M,Int}(-1, Monomial(Variable(x, IMAG_PART)))
+    else
+        return MA.Zero()
+    end
+end
+
 function mergevars_to!(
     vars::Vector{PV},
     varsvec::Vector{Vector{PV}},