Complex-valued variables

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) ? cpNone : cpFull)
+end
+function MP.similar_variable(::Type{<:PolyType{V,M}}, s::Symbol) where {V,M}
+    return Variable(string(s), V, M, cpNone) # 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

@@ -91,6 +91,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
 
@@ -121,7 +143,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,73 @@ 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 may also assign a complex value to its conjugate, or just the real or imaginary parts
+    # Any combination of z, conj(z), real(z), imag(z), imag(conj(z)) can occur in either the polynomial or the substitution,
+    # and we must handle all of them correctly.
+    # Note: This may or may not work... Issues can arise if the substitutions contain the real and imaginary (or only one of
+    # those) of a variable separately whenever vals is not of the correct type:
+    # - Unless subs() originally had a polynomial-valued rhs, vals will be scalars, monomials, or terms. So when we try to
+    #   assign a polynomial to its value (which is necessary, as the one-step substitution of only the real or only the
+    #   imaginary part is incomplete), this is an impossible conversion.
+    # - The coefficients in vals might not be complex-valued; but to assign only a part of the variable, we necessarily need to
+    #   introduce an explicit imaginary coefficient to the value.
+    # Currently, we don't do anything to catch these errors.
+    if s.first.kind == cpNone
+        for j in eachindex(vars)
+            if vars[j] == s.first
+                vals[j] = s.second
+                C == Commutative && break
+            end
+        end
+    else
+        for j in eachindex(vars)
+            if vars[j].variable_order.order.id ==
+               s.first.variable_order.order.id
+                if s.first.kind == cpFull || s.first.kind == cpConj
+                    value = s.first.kind == cpConj ? conj(s.second) : s.second
+                    if vars[j].kind == cpFull
+                        vals[j] = value
+                    elseif vars[j].kind == cpConj
+                        vals[j] = conj(value)
+                    elseif vars[j].kind == cpReal
+                        vals[j] = real(value)
+                    else
+                        vals[j] = imag(value)
+                    end
+                elseif s.first.kind == cpReal
+                    isreal(s.second) || error(
+                        "Cannot assign a complex value to the real part of an expression",
+                    )
+                    value = real(s.second) # just to make sure the type is correct
+                    if vars[j].kind == cpFull
+                        vals[j] = value + im * imag(vals[j])
+                    elseif vars[j].kind == cpConj
+                        vals[j] = value - im * imag(vals[j])
+                    elseif vars[j].kind == cpReal
+                        vals[j] = value
+                    end
+                    # else we know the real part but use the imaginary part; do nothing
+                else
+                    @assert(s.first.kind == cpImag)
+                    isreal(s.second) || error(
+                        "Cannot assign a complex value to the imaginary part of an expression",
+                    )
+                    value = real(s.second) # just to make sure the type is correct
+                    if vars[j].kind == cpFull
+                        vals[j] = real(vals[j]) + im * value
+                    elseif vars[j].kind == cpConj
+                        vals[j] = real(vals[j]) - im * value
+                    elseif vars[j].kind == cpImag
+                        vals[j] = value
+                    end
+                    # else we know the imaginary part but use the real part; do nothing
+                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, @polycvar
 
-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 ? cpFull : cpNone
             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})), cpNone)
+    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})), cpNone)
+    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 polycvar(args...)
+    pos_args, variable_order, monomial_order, complex_kind =
+        _extract_kw_args(args, :($(Commutative{CreationOrder})), cpFull)
+    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,37 @@ function instantiate(::Type{NonCommutative{O}}) where {O}
     return NonCommutative(instantiate(O))
 end
 
+@enum ComplexKind cpNone cpFull cpConj cpReal cpImag
+
 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 = cpNone,
+    ) 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 = cpNone,
+    ) where {V<:AbstractVariableOrdering,M<:MP.AbstractMonomialOrdering}
+        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}
-        return new{V,M}(convert(String, name), instantiate(V))
+        (from.kind == cpNone && kind == cpNone) ||
+            (from.kind != cpNone && kind != cpNone) ||
+            error("Cannot change the complex type of an existing variable")
+        return new{V,M}(from.name, kind, from.variable_order)
     end
 end
 
@@ -150,6 +197,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 != cpFull && x.kind != cpConj
+MP.isrealpart(x::Variable{C}) where {C} = x.kind == cpReal
+MP.isimagpart(x::Variable{C}) where {C} = x.kind == cpImag
+MP.isconj(x::Variable{C}) where {C} = x.kind == cpConj
+function MP.ordinary_variable(x::Variable)
+    return x.kind == cpNone || x.kind == cpFull ? x : Variable(x, cpFull)
+end
+
+function Base.conj(x::Variable)
+    if x.kind == cpFull
+        return Variable(x, cpConj)
+    elseif x.kind == cpConj
+        return Variable(x, cpFull)
+    else
+        return x
+    end
+end
+
+function Base.real(x::Variable)
+    if x.kind == cpFull || x.kind == cpConj
+        return Variable(x, cpReal)
+    else
+        return x
+    end
+end
+function Base.imag(x::Variable{V,M}) where {V,M}
+    if x.kind == cpFull
+        return Variable(x, cpImag)
+    elseif x.kind == cpConj
+        return _Term{V,M,Int}(-1, Monomial(Variable(x, cpImag)))
+    else
+        return MA.Zero()
+    end
+end
+
 function mergevars_to!(
     vars::Vector{PV},
     varsvec::Vector{Vector{PV}},