Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Complex-valued variables #213

Merged
merged 41 commits into from
Oct 6, 2023
Merged
Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
41 commits
Select commit Hold shift + click to select a range
a6e9b84
First draft for complex-valued interface
projekter Jun 11, 2022
a84d217
Adjust imaginary part
projekter Jun 12, 2022
26a471f
Merge remote-tracking branch 'upstream/master'
projekter Dec 26, 2022
73a0bc6
Slight updates
projekter Dec 26, 2022
98905be
Add exports
projekter Dec 26, 2022
9ab15da
Rename ordvar to ordinary_variable
projekter Apr 21, 2023
9679a99
Change assert to error
projekter May 9, 2023
2606f9c
Code style
projekter May 9, 2023
5df604c
Remove real_coefficient_type_polynomial_like
projekter May 9, 2023
9280e2d
Change comparisons to comparison functions
projekter May 9, 2023
8e0424f
Explain _show in comment
projekter May 9, 2023
b897034
Shortcuts for conj
projekter May 9, 2023
064f0e5
Add explicit isreal/iscomplex for monomial vectors
projekter May 9, 2023
f4606e7
Fix
projekter May 9, 2023
7a4c107
Improve real/imag of APL
projekter May 9, 2023
5889093
Replace zip by powers
projekter May 9, 2023
ef5ad13
Fix
projekter May 10, 2023
8e463a8
Merge branch 'master' of https://github.com/projekter/MultivariatePol…
projekter May 10, 2023
0687b19
Explain choice of function signature
projekter May 19, 2023
cc76af5
Merge branch 'master' of https://github.com/JuliaAlgebra/Multivariate…
projekter May 19, 2023
76fa39a
Format
projekter May 19, 2023
718a92e
Documentation
projekter May 22, 2023
7ba9cc6
Merge branch 'master' of https://github.com/JuliaAlgebra/Multivariate…
projekter May 28, 2023
8c53052
Fix iscomplex detection
projekter May 28, 2023
e7d57d6
Docstring update, new function spelling
projekter May 28, 2023
d5784cc
Test for complex case
projekter May 28, 2023
a72912f
Remove missing references to Base functions
projekter May 28, 2023
1ddcd0c
Format test
projekter May 28, 2023
cf96129
Merge branch 'master' of https://github.com/JuliaAlgebra/Multivariate…
projekter Jul 7, 2023
7223c0c
Format
projekter Jul 7, 2023
cfa73d0
Compatibility in tests
projekter Jul 10, 2023
12ad263
Fix adjoint issues with new interface
projekter Jul 10, 2023
60f1a58
Revert APL adjoint implementation (mutable coefficients would be alia…
projekter Jul 22, 2023
25b8e02
Add copy_if_mutable safeguards
projekter Aug 17, 2023
a09ca39
Remove complex branch CI
projekter Aug 17, 2023
ebc3c6e
Update src/complex.jl
projekter Aug 24, 2023
7624bd7
Update src/complex.jl
projekter Aug 24, 2023
e479b43
Update src/complex.jl
projekter Aug 24, 2023
a8ac55e
Avoid splatting
projekter Aug 24, 2023
c19afba
Remove iscomplex
projekter Sep 1, 2023
3190908
Improve documentation
projekter Sep 12, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions src/MultivariatePolynomials.jl
Original file line number Diff line number Diff line change
Expand Up @@ -80,6 +80,7 @@ include("hash.jl")
include("promote.jl")
include("conversion.jl")

include("complex.jl")
include("operators.jl")
include("comparison.jl")

Expand Down
285 changes: 285 additions & 0 deletions src/complex.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,285 @@
export iscomplex, isrealpart, isimagpart, isconj, ordinary_variable, degree_complex, halfdegree, mindegree_complex,
minhalfdegree, maxdegree_complex, maxhalfdegree, extdegree_complex, exthalfdegree

"""
iscomplex(x::AbstractVariable)

Return whether a given variable was declared as a complex-valued variable (also their conjugates are complex, but their real
and imaginary parts are not).
By default, all variables are real-valued.
"""
iscomplex(::AbstractVariable) = false
projekter marked this conversation as resolved.
Show resolved Hide resolved
Base.isreal(v::AbstractPolynomialLike) = !iscomplex(v)
"""
isrealpart(x::AbstractVariable)

Return whether the given variable is the real part of a complex-valued variable.

See also [`iscomplex`](@ref iscomplex), [`isimagpart`](@ref isimagpart), [`isconj`](@ref isconj).
"""
isrealpart(::AbstractVariable) = false
projekter marked this conversation as resolved.
Show resolved Hide resolved
"""
isimagpart(x::AbstractVariable)

Return whether the given variable is the imaginary part of a complex-valued variable.

See also [`iscomplex`](@ref iscomplex), [`isrealpart`](@ref isrealpart), [`isconj`](@ref isconj).
"""
isimagpart(::AbstractVariable) = false
"""
isconj(x::AbstractVariable)

Return whether the given variable is obtained by conjugating a user-defined complex-valued variable.

See also [`iscomplex`](@ref iscomplex), [`isrealpart`](@ref isrealpart), [`isimagpart`](@ref isimagpart).
"""
isconj(::AbstractVariable) = false
"""
ordinary_variable(x::Union{AbstractVariable, AbstractVector{<:AbstractVariable}})

Given some (complex-valued) variable that was transformed by conjugation, taking its real part, or taking its
imaginary part, return the original variable as it was defined by the user.

See also [`conj`](@ref conj), [`real`](@ref), [`imag`](@ref).
"""
ordinary_variable(x::AbstractVariable) = x

"""
conj(x::AbstractVariable)

Return the complex conjugate of a given variable if it was declared as a complex variable; else return the
variable unchanged.

conj(x::AbstractMonomial)
conj(x::AbstractVector{<:AbstractMonomial})
conj(x::AbstractTerm)
conj(x::AbstractPolynomial)

Return the complex conjugate of `x` by applying conjugation to all coefficients and variables.

See also [`iscomplex`](@ref iscomplex), [`isconj`](@ref isconj).
"""
Base.conj(x::AbstractVariable) = x

"""
real(x::AbstractVariable)

Return the real part of a given variable if it was declared as a complex variable; else return the variable
unchanged.

real(x::AbstractMonomial)
real(x::AbstractVector{<:AbstractMonomial})
real(x::AbstractTerm)
real(x::AbstractPolynomial)

Return the real part of `x` by applying `real` to all coefficients and variables; for this purpose, every complex-valued
variable is decomposed into its real- and imaginary parts.

See also [`iscomplex`](@ref iscomplex), [`isrealpart`](@ref isrealpart), [`imag`](@ref).
"""
Base.real(x::AbstractVariable) = x

"""
imag(x::AbstractVariable)

Return the imaginary part of a given variable if it was declared as a complex variable; else return zero.

imag(x::AbstractMonomial)
imag(x::AbstractVector{<:AbstractMonomial})
imag(x::AbstractTerm)
imag(x::AbstractPolynomial)

Return the imaginary part of `x` by applying `imag` to all coefficients and variables; for this purpose, every complex-valued
variable is decomposed into its real- and imaginary parts.

See also [`iscomplex`](@ref iscomplex), [`isimagpart`](@ref isimagpart), [`real`](@ref).
"""
Base.imag(::AbstractVariable) = MA.Zero()

# extend to higher-level elements. We make all those type-stable (but we need convert, as the construction method may
# deliver simpler types than the inputs if they were deliberately casted, e.g., term to monomial)
function iscomplex(p::AbstractPolynomialLike)
for v in variables(p)
if !iszero(maxdegree(p, v)) && iscomplex(v)
return true
end
end
return false
end
Base.conj(x::M) where {M<:AbstractMonomial} =
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would add a shortcut here in case isreal(x) is true

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same for the three methods below

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added the shortcuts; my initial rationale for not doing so was that the checks themselves are not just a query of a flag, but instead an iteration over potentially all the involved variables, so this is not a super cheap operation.
There's also the question of the function contract: Now it may return either the same object (for a no-op) or a new operation. Since MP is abstract, we don't know whether the underlying implementations are immutable/have addresses.
I added an isreal/iscomplex check for monomial vectors as well to have this explicitly and updated this in DynamicPolynomials, where it can be done much more efficiently.

convert(M, reduce(*, conj(var)^exp for (var, exp) in zip(variables(x), exponents(x)); init=constantmonomial(x)))
projekter marked this conversation as resolved.
Show resolved Hide resolved
Base.conj(x::V) where {V<:AbstractVector{<:AbstractMonomial}} = monovec(conj.(x))
Base.conj(x::T) where {T<:AbstractTerm} = convert(T, conj(coefficient(x)) * conj(monomial(x)))
Base.conj(x::P) where {P<:AbstractPolynomial} = nterms(x) == 0 ? x : convert(P, sum(conj(t) for t in x))

real_coefficient_type_polynomial_like(::AbstractPolynomialLike{R}) where {R<:Real} = R
projekter marked this conversation as resolved.
Show resolved Hide resolved
real_coefficient_type_polynomial_like(::AbstractPolynomialLike{Complex{R}}) where {R<:Real} = R
# Real and imaginary parts are harder to realize. The real part of a monomial can easily be a polynomial.
for fun in [:real, :imag]
eval(quote
function Base.$fun(x::Union{AbstractMonomial,AbstractTerm})
projekter marked this conversation as resolved.
Show resolved Hide resolved
# We replace every complex variable by its decomposition into real and imaginary part
substs = [var => real(var) + 1im * imag(var) for var in variables(x) if iscomplex(var)]
# subs throws an error if it doesn't receive at least one substitution
full_version = length(substs) > 0 ? subs(x, substs...) : polynomial(x)
projekter marked this conversation as resolved.
Show resolved Hide resolved
# Now everything that is imaginary can be recognized by its coefficient
return convert(
polynomialtype(full_version, real_coefficient_type_polynomial_like(full_version)),
mapcoefficients!($fun, full_version)
)
end
Base.$fun(x::AbstractVector{<:AbstractMonomial}) = map(Base.$fun, x)
Base.$fun(x::AbstractPolynomial{T}) where {T} = nterms(x) == 0 ?
projekter marked this conversation as resolved.
Show resolved Hide resolved
projekter marked this conversation as resolved.
Show resolved Hide resolved
zero(polynomialtype(x, real_coefficient_type_polynomial_like(x))) : sum($fun(t) for t in x)
end)
end

# Also give complex-valued degree definitions. We choose not to overwrite degree, as this will lead to issues in monovecs
# and their sorting. So now there are two ways to calculate degrees: strictly by considering all variables independently,
# and also by looking at their complex structure.
"""
degree_complex(t::AbstractTermLike)

Return the _total complex degree_ of the monomial of the term `t`, i.e., the maximum of the total degree of the declared
variables in `t` and the total degree of the conjugate variables in `t`.
To be well-defined, the monomial must not contain real parts or imaginary parts of variables.

degree_complex(t::AbstractTermLike, v::AbstractVariable)

Returns the exponent of the variable `v` or its conjugate in the monomial of the term `t`, whatever is larger.

See also [`isconj`](@ref isconj).
"""
function degree_complex(t::AbstractTermLike)
vars = variables(t)
@assert(!any(isrealpart, vars) && !any(isimagpart, vars))
grouping = isconj.(vars)
exps = exponents(t)
return max(sum(exps[grouping]), sum(exps[map(!, grouping)]))
end
degree_complex(t::AbstractTermLike, var::AbstractVariable) = degree_complex(monomial(t), var)
degree_complex(v::AbstractVariable, var::AbstractVariable) = (v == var ? 1 : 0)
function degree_complex(m::AbstractMonomial, v::AbstractVariable)
deg = 0
deg_c = 0
c_v = conj(v)
for (var, exp) in powers(m)
@assert(!isrealpart(var) && !isimagpart(var))
if var == v
deg += exp
elseif var == c_v
deg_c += exp
end
end
return max(deg, deg_c)
end

"""
halfdegree(t::AbstractTermLike)

Return the equivalent of `ceil(degree(t)/2)`` for real-valued terms or `degree_complex(t)` for terms with only complex
variables; however, respect any mixing between complex and real-valued variables.
"""
function halfdegree(t::AbstractTermLike)
realdeg = 0
cpdeg = 0
conjdeg = 0
for (var, exp) in powers(t)
if iscomplex(var)
if isconj(var)
conjdeg += exp
else
@assert(!isrealpart(var) && !isimagpart(var))
cpdeg += exp
end
else
realdeg += exp
end
end
return ((realdeg + 1) >> 1) + max(cpdeg, conjdeg)
end

"""
mindegree_complex(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}})

Return the minimal total complex degree of the monomials of `p`, i.e., `minimum(degree_complex, terms(p))`.

mindegree_complex(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, v::AbstractVariable)

Return the minimal complex degree of the monomials of `p` in the variable `v`, i.e., `minimum(degree_complex.(terms(p), v))`.
"""
mindegree_complex(X::AbstractVector{<:AbstractTermLike}, args...) =
isempty(X) ? 0 : minimum(t -> degree_complex(t, args...), X, init=0)
mindegree_complex(p::AbstractPolynomialLike, args...) = mindegree_complex(terms(p), args...)

"""
minhalfdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}})

Return the minmal half degree of the monomials of `p`, i.e., `minimum(halfdegree, terms(p))`
"""
minhalfdegree(X::AbstractVector{<:AbstractTermLike}, args...) = isempty(X) ? 0 : minimum(halfdegree, X, init=0)
minhalfdegree(p::AbstractPolynomialLike) = minhalfdegree(terms(p))

"""
maxdegree_complex(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}})

Return the maximal total complex degree of the monomials of `p`, i.e., `maximum(degree_complex, terms(p))`.

maxdegree_complex(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, v::AbstractVariable)

Return the maximal complex degree of the monomials of `p` in the variable `v`, i.e., `maximum(degree_complex.(terms(p), v))`.
"""
maxdegree_complex(X::AbstractVector{<:AbstractTermLike}, args...) =
mapreduce(t -> degree_complex(t, args...), max, X, init=0)
maxdegree_complex(p::AbstractPolynomialLike, args...) = maxdegree_complex(terms(p), args...)

"""
maxhalfdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}})

Return the maximal half degree of the monomials of `p`, i.e., `maximum(halfdegree, terms(p))`
"""
maxhalfdegree(X::AbstractVector{<:AbstractTermLike}, args...) = isempty(X) ? 0 : maximum(halfdegree, X, init=0)
maxhalfdegree(p::AbstractPolynomialLike) = maxhalfdegree(terms(p))

"""
extdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}})

Returns the extremal total complex degrees of the monomials of `p`, i.e., `(mindegree_complex(p), maxdegree_complex(p))`.

extdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}}, v::AbstractVariable)

Returns the extremal complex degrees of the monomials of `p` in the variable `v`, i.e.,
`(mindegree_complex(p, v), maxdegree_complex(p, v))`.
"""
extdegree_complex(p::Union{AbstractPolynomialLike,AbstractVector{<:AbstractTermLike}}, args...) =
(mindegree_complex(p, args...), maxdegree_complex(p, args...))

"""
exthalfdegree(p::Union{AbstractPolynomialLike, AbstractVector{<:AbstractTermLike}})

Return the extremal half degree of the monomials of `p`, i.e., `(minhalfdegree(p), maxhalfdegree(p))`
"""
exthalfdegree(p::Union{AbstractPolynomialLike,AbstractVector{<:AbstractTermLike}}) = (minhalfdegree(p), maxhalfdegree(p))

function ordinary_variable(x::AbstractVector{<:AbstractVariable})
# let's assume the number of elements in x is small, else a conversion to a dict (probably better OrderedDict) would
# be better
results = similar(x, 0)
sizehint!(results, length(x))
j = 0
for el in x
ov = ordinary_variable(el)
found = false
for i in 1:j
if results[i] == ov
found = true
break
end
end
if !found
push!(results, ov)
j += 1
end
end
return results
end
8 changes: 5 additions & 3 deletions src/operators.jl
Original file line number Diff line number Diff line change
Expand Up @@ -322,11 +322,13 @@ function _promote_terms(t1::AbstractTerm{LinearAlgebra.UniformScaling{T}}, t2::A
return promote(t1, t2)
end

LinearAlgebra.adjoint(v::AbstractVariable) = v
LinearAlgebra.adjoint(m::AbstractMonomial) = m
LinearAlgebra.adjoint(t::AbstractTerm) = _term(LinearAlgebra.adjoint(coefficient(t)), monomial(t))
LinearAlgebra.adjoint(v::AbstractVariable) = conj(v)
LinearAlgebra.adjoint(m::AbstractMonomial) = conj(m)
LinearAlgebra.adjoint(t::AbstractTerm) = _term(LinearAlgebra.adjoint(coefficient(t)), LinearAlgebra.adjoint(monomial(t)))
LinearAlgebra.adjoint(p::AbstractPolynomialLike) = polynomial(map(LinearAlgebra.adjoint, terms(p)))
LinearAlgebra.adjoint(r::RationalPoly) = adjoint(numerator(r)) / adjoint(denominator(r))
LinearAlgebra.hermitian_type(::Type{T}) where {T<:AbstractPolynomialLike} = T
LinearAlgebra.hermitian(v::AbstractPolynomialLike, ::Symbol) = (@assert(!iscomplex(v)); v)
projekter marked this conversation as resolved.
Show resolved Hide resolved

LinearAlgebra.transpose(v::AbstractVariable) = v
LinearAlgebra.transpose(m::AbstractMonomial) = m
Expand Down
23 changes: 20 additions & 3 deletions src/show.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,14 +22,31 @@ Base.show(io::IO, p::TypesWithShow) = show(io, MIME"text/plain"(), p)
# VARIABLES
function _show(io::IO, mime::MIME, var::AbstractVariable)
base, indices = name_base_indices(var)
if isempty(indices)
print(io, base)
if isconj(var)
for c in base
print(io, c, '\u0305')
end
else
print(io, base)
end
if !isempty(indices)
print_subscript(io, mime, indices)
end
isrealpart(var) && print(io, "ᵣ")
isimagpart(var) && print(io, "ᵢ")
end

function _show(io::IO, mime::MIME"text/print", var::AbstractVariable)
if isconj(var)
for c in name(var)
print(io, c, '\u0305')
projekter marked this conversation as resolved.
Show resolved Hide resolved
end
else
print(io, name(var))
end
isrealpart(var) && print(io, "ᵣ")
projekter marked this conversation as resolved.
Show resolved Hide resolved
isimagpart(var) && print(io, "ᵢ")
end
_show(io::IO, mime::MIME"text/print", var::AbstractVariable) = print(io, name(var))

function print_subscript(io::IO, ::MIME"text/latex", index)
print(io, "_{", join(index, ","), "}")
Expand Down