From 121bd58ba0e7e7c90535980dbaab3db9d61d38e4 Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Mon, 13 May 2024 22:04:36 +0200 Subject: [PATCH] Complex-valued variables (#121) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * 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 --- Project.toml | 2 +- src/DynamicPolynomials.jl | 10 +-- src/comp.jl | 9 ++- src/mono.jl | 7 ++ src/monomial_vector.jl | 28 +++++++- src/subs.jl | 53 ++++++++++----- src/var.jl | 133 +++++++++++++++++++++++++++++++++----- test/poly.jl | 21 +++++- 8 files changed, 220 insertions(+), 43 deletions(-) diff --git a/Project.toml b/Project.toml index 40cb98a..b9c52d1 100644 --- a/Project.toml +++ b/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" diff --git a/src/DynamicPolynomials.jl b/src/DynamicPolynomials.jl index 5df9ab0..3cc41dd 100644 --- a/src/DynamicPolynomials.jl +++ b/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") diff --git a/src/comp.jl b/src/comp.jl index d9c0fbd..62e7edd 100644 --- a/src/comp.jl +++ b/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 diff --git a/src/mono.jl b/src/mono.jl index b3d97c6..02ef0c9 100644 --- a/src/mono.jl +++ b/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 diff --git a/src/monomial_vector.jl b/src/monomial_vector.jl index 991f4b1..b9d3e51 100644 --- a/src/monomial_vector.jl +++ b/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 diff --git a/src/subs.jl b/src/subs.jl index 216ca00..f856522 100644 --- a/src/subs.jl +++ b/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) diff --git a/src/var.jl b/src/var.jl index b7eb350..3e4fc66 100644 --- a/src/var.jl +++ b/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,6 +42,7 @@ function buildpolyvar(var, variable_order, monomial_order) $(esc(varname)) = polyarrayvar( $(variable_order), $(monomial_order), + $(complex_kind), $prefix, $(esc.(var.args[2:end])...), ) @@ -35,11 +50,12 @@ function buildpolyvar(var, variable_order, monomial_order) 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 @@ -47,7 +63,7 @@ function buildpolyvars(args, variable_order, monomial_order) 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,6 +72,8 @@ 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 @@ -63,22 +81,33 @@ function _extract_kw_args(args, variable_order) 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}}, diff --git a/test/poly.jl b/test/poly.jl index 80c55a3..cb0d7b2 100644 --- a/test/poly.jl +++ b/test/poly.jl @@ -55,10 +55,16 @@ p.x[2] == x @inferred DynamicPolynomials.Polynomial(i -> float(i), [x, x * x]) - @inferred DynamicPolynomials.Polynomial(i -> float(i), MonomialVector([x * x, x])) + @inferred DynamicPolynomials.Polynomial( + i -> float(i), + MonomialVector([x * x, x]), + ) for p in ( DynamicPolynomials.Polynomial(i -> float(i), [x, x * x]), - DynamicPolynomials.Polynomial(i -> float(i), MonomialVector([x * x, x])), + DynamicPolynomials.Polynomial( + i -> float(i), + MonomialVector([x * x, x]), + ), ) @test typeof(p) == PTF @test p.a == [1.0, 2.0] @@ -130,5 +136,16 @@ ) @test_throws err t(y => 1) @test_throws err p(y => 1) + + @complex_polyvar z + pc = z^3 + 2real(z) - 7imag(z^4) + @test pc(z => 2 + 3im) == 798 + 9im + err = ArgumentError( + "Substitution with complex variables requires the ordinary_variable in the substitution specification" + ) + @test_throws err pc(conj(z) => 2 - 3im) + @test_throws err pc(real(z) => 2) + @test_throws err pc(imag(z) => 3) + @test real(pc)(z => 2 + 3im) == 798 end end