From b46b905e8e8c0f57508642360add67c2929e74d7 Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Mon, 26 Dec 2022 14:47:31 +0100 Subject: [PATCH 01/23] First step to implement complex variables functionality --- src/DynamicPolynomials.jl | 6 ++-- src/comp.jl | 10 +++++-- src/monovec.jl | 19 +++++++++++++ src/subs.jl | 42 ++++++++++++++++++++------- src/var.jl | 60 +++++++++++++++++++++++++++++++++++++-- 5 files changed, 120 insertions(+), 17 deletions(-) diff --git a/src/DynamicPolynomials.jl b/src/DynamicPolynomials.jl index 69d8f27..72662f5 100644 --- a/src/DynamicPolynomials.jl +++ b/src/DynamicPolynomials.jl @@ -40,8 +40,10 @@ MP.polynomialtype(::Union{PolyType{C}, Type{<:PolyType{C}}}, ::Type{T}) where {C _vars(p::AbstractArray{<:PolyType}) = mergevars(_vars.(p))[1] MP.variables(p::Union{PolyType, MonomialVector, AbstractArray{<:PolyType}}) = _vars(p) # tuple(_vars(p)) MP.nvariables(p::Union{PolyType, MonomialVector, AbstractArray{<:PolyType}}) = length(_vars(p)) -MP.similarvariable(p::Union{PolyType{C}, Type{<:PolyType{C}}}, ::Type{Val{V}}) where {C, V} = PolyVar{C}(string(V)) -MP.similarvariable(p::Union{PolyType{C}, Type{<:PolyType{C}}}, V::Symbol) where {C} = PolyVar{C}(string(V)) +MP.similarvariable(p::Type{<:PolyType{C}}, ::Type{Val{V}}) where {C, V} = PolyVar{C}(string(V)) +MP.similarvariable(p::Type{<:PolyType{C}}, V::Symbol) where {C} = PolyVar{C}(string(V)) +MP.similarvariable(p::PolyVar{C}, ::Type{Val{V}}) where {C,V} = PolyVar{C}(string(V), p.kind) +MP.similarvariable(p::PolyVar{C}, V::Symbol) where {C} = similarvariable(p, Val{V}) include("promote.jl") diff --git a/src/comp.jl b/src/comp.jl index ae03217..fc17113 100644 --- a/src/comp.jl +++ b/src/comp.jl @@ -22,10 +22,16 @@ end # Comparison of PolyVar function (==)(x::PolyVar{C}, y::PolyVar{C}) where C - x.id == y.id + x.id == y.id && x.kind == y.kind end -Base.isless(x::PolyVar{C}, y::PolyVar{C}) where C = isless(y.id, x.id) +function Base.isless(x::PolyVar{C}, y::PolyVar{C}) where {C} + if x.id == y.id + return isless(y.kind, x.kind) + else + return isless(y.id, x.id) + end +end # Comparison of Monomial diff --git a/src/monovec.jl b/src/monovec.jl index dc1546c..953e8f7 100644 --- a/src/monovec.jl +++ b/src/monovec.jl @@ -70,6 +70,25 @@ end MultivariatePolynomials.extdegree(x::MonomialVector) = isempty(x) ? (0, 0) : extrema(sum.(x.Z)) MultivariatePolynomials.mindegree(x::MonomialVector) = isempty(x) ? 0 : minimum(sum.(x.Z)) MultivariatePolynomials.maxdegree(x::MonomialVector) = isempty(x) ? 0 : maximum(sum.(x.Z)) +# 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 MP.$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 _vars(m::Union{Monomial, MonomialVector}) = m.vars diff --git a/src/subs.jl b/src/subs.jl index be74738..d1df033 100644 --- a/src/subs.jl +++ b/src/subs.jl @@ -7,17 +7,39 @@ function Base.setindex!(sv::SafeValues, v, i::Int) sv.values[i] = v end -function fillmap!(vals, vars::Vector{PolyVar{false}}, s::MP.Substitution) - for j in eachindex(vars) - if vars[j] == s.first - vals[j] = s.second +function fillmap!(vals, vars::Vector{PolyVar{C}}, s::MP.Substitution) 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. + if !iscomplex(s.first) + for j in eachindex(vars) + if vars[j] == s.first + vals[j] = s.second + end + end + else + for j in eachindex(vars) + if vars[j].id == s.first.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 == vars[j].kind + vals[j] = real(s.second) # just to make sure the type is correct + end + # else: # assignment to parts of a complex variable, but the full appears in the polynomial + # We don't support this for the moment (would require more complex logic, since two substitutions may give the + # full value, or one part is left over - but can vals hold such a leftover?). + end end - end -end -function fillmap!(vals, vars::Vector{PolyVar{true}}, 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 0f66f8f..ec04d70 100644 --- a/src/var.jl +++ b/src/var.jl @@ -1,4 +1,4 @@ -export PolyVar, @polyvar, @ncpolyvar +export PolyVar, @polyvar, @ncpolyvar, @polycvar export polyvecvar @@ -41,16 +41,31 @@ macro ncpolyvar(args...) :($(foldl((x,y) -> :($x; $y), exprs, init=:())); $(Expr(:tuple, esc.(vars)...))) end +macro polycvar(args...) + vars, exprs = buildpolyvars(PolyVarComplex{true}, args) + return :($(foldl((x, y) -> :($x; $y), exprs, init=:())); $(Expr(:tuple, esc.(vars)...))) +end + +@enum ComplexKind cpNone cpFull cpConj cpReal cpImag + struct PolyVar{C} <: AbstractVariable id::Int name::String + kind::ComplexKind - function PolyVar{C}(name::AbstractString) where {C} + function PolyVar{C}(name::AbstractString, kind::ComplexKind=cpNone) where {C} # gensym returns something like Symbol("##42") # we first remove "##" and then parse it into an Int id = parse(Int, string(gensym())[3:end]) - new(id, convert(String, name)) + new(id, convert(String, name), kind) end + function PolyVar{C}(id::Int, name::String, kind::ComplexKind) where {C} + new(id, name, kind) + end +end + +struct PolyVarComplex{C} + PolyVarComplex{C}(name::AbstractString) where {C} = PolyVar{C}(name, cpFull) end Base.hash(x::PolyVar, u::UInt) = hash(x.id, u) @@ -71,6 +86,45 @@ _vars(v::PolyVar) = [v] iscomm(::Type{PolyVar{C}}) where {C} = C +MP.iscomplex(x::PolyVar{C}) where {C} = x.kind == cpFull || x.kind == cpConj +MP.isrealpart(x::PolyVar{C}) where {C} = x.kind == cpReal +MP.isimagpart(x::PolyVar{C}) where {C} = x.kind == cpImag +MP.isconj(x::PolyVar{C}) where {C} = x.kind == cpConj +MP.ordvar(x::PolyVar{C}) where {C} = x.kind == cpNone || x.kind == cpFull ? x : PolyVar{C}(x.id, x.name, cpFull) + +function Base.conj(x::PolyVar{C}) where {C} + if x.kind == cpFull + return PolyVar{C}(x.id, x.name, cpConj) + elseif x.kind == cpConj + return PolyVar{C}(x.id, x.name, cpFull) + else + return x + end +end +# for efficiency reasons +function Base.conj(x::Monomial{true}) + cv = conj.(x.vars) + perm = sortperm(cv, rev=true) + return Monomial{true}(cv[perm], x.z[perm]) +end + +function Base.real(x::PolyVar{C}) where {C} + if x.kind == cpFull || x.kind == cpConj + return PolyVar{C}(x.id, x.name, cpReal) + else + return x + end +end +function Base.imag(x::PolyVar{C}) where {C} + if x.kind == cpFull + return PolyVar{C}(x.id, x.name, cpImag) + elseif x.kind == cpConj + return Term{C,Int}(-1, PolyVar{C}(x.id, x.name, cpImag)) + else + return MA.Zero() + end +end + function mergevars_to!(vars::Vector{PV}, varsvec::Vector{Vector{PV}}) where {PV<:PolyVar} empty!(vars) n = length(varsvec) From a5cb35ddfc784dae8b76fd50dd2ce7ed8143640b Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Mon, 26 Dec 2022 14:53:44 +0100 Subject: [PATCH 02/23] Move monomial function to correct file --- src/mono.jl | 7 +++++++ src/var.jl | 6 ------ 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/mono.jl b/src/mono.jl index 699ba1b..7815cd9 100644 --- a/src/mono.jl +++ b/src/mono.jl @@ -80,3 +80,10 @@ MP.monomial(m::Monomial) = m # end # i > length(m1.z) #end + +# for efficiency reasons +function Base.conj(x::Monomial{true}) + cv = conj.(x.vars) + perm = sortperm(cv, rev=true) + return Monomial{true}(cv[perm], x.z[perm]) +end \ No newline at end of file diff --git a/src/var.jl b/src/var.jl index ec04d70..1c67dff 100644 --- a/src/var.jl +++ b/src/var.jl @@ -101,12 +101,6 @@ function Base.conj(x::PolyVar{C}) where {C} return x end end -# for efficiency reasons -function Base.conj(x::Monomial{true}) - cv = conj.(x.vars) - perm = sortperm(cv, rev=true) - return Monomial{true}(cv[perm], x.z[perm]) -end function Base.real(x::PolyVar{C}) where {C} if x.kind == cpFull || x.kind == cpConj From 91583e2aad255e9d3542afb969dd089110d84c12 Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Fri, 21 Apr 2023 10:32:06 +0200 Subject: [PATCH 03/23] Rename ordvar to ordinary_variable --- src/var.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/var.jl b/src/var.jl index 1c67dff..60374c2 100644 --- a/src/var.jl +++ b/src/var.jl @@ -90,7 +90,7 @@ MP.iscomplex(x::PolyVar{C}) where {C} = x.kind == cpFull || x.kind == cpConj MP.isrealpart(x::PolyVar{C}) where {C} = x.kind == cpReal MP.isimagpart(x::PolyVar{C}) where {C} = x.kind == cpImag MP.isconj(x::PolyVar{C}) where {C} = x.kind == cpConj -MP.ordvar(x::PolyVar{C}) where {C} = x.kind == cpNone || x.kind == cpFull ? x : PolyVar{C}(x.id, x.name, cpFull) +MP.ordinary_variable(x::PolyVar{C}) where {C} = x.kind == cpNone || x.kind == cpFull ? x : PolyVar{C}(x.id, x.name, cpFull) function Base.conj(x::PolyVar{C}) where {C} if x.kind == cpFull From b3d8d736df52b53a0f2fc682bdb8ab8a5046fa4e Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Tue, 9 May 2023 18:21:04 +0200 Subject: [PATCH 04/23] Specialized complex routines for MonomialVector --- src/monovec.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/monovec.jl b/src/monovec.jl index 953e8f7..2e97da8 100644 --- a/src/monovec.jl +++ b/src/monovec.jl @@ -67,9 +67,9 @@ function Base.iterate(x::MonomialVector, state::Int) state < length(x) ? (x[state+1], state+1) : nothing end -MultivariatePolynomials.extdegree(x::MonomialVector) = isempty(x) ? (0, 0) : extrema(sum.(x.Z)) -MultivariatePolynomials.mindegree(x::MonomialVector) = isempty(x) ? 0 : minimum(sum.(x.Z)) -MultivariatePolynomials.maxdegree(x::MonomialVector) = isempty(x) ? 0 : maximum(sum.(x.Z)) +MP.extdegree(x::MonomialVector) = isempty(x) ? (0, 0) : extrema(sum.(x.Z)) +MP.mindegree(x::MonomialVector) = isempty(x) ? 0 : minimum(sum.(x.Z)) +MP.maxdegree(x::MonomialVector) = isempty(x) ? 0 : maximum(sum.(x.Z)) # Complex-valued degrees for monomial vectors for (fun, call, def, ret) in [ (:extdegree_complex, :extrema, (0, 0), :((min(v1, v2), max(v1, v2)))), @@ -89,6 +89,9 @@ for (fun, call, def, ret) in [ end end) end +# faster complex-related functions +MP.iscomplex(x::MonomialVector) = any(iscomplex, x.vars) +MP.conj(x::MonomialVector) = MonomialVector(conj.(x.vars), x.Z) # TODO: do we want to alias Z or copy Z? _vars(m::Union{Monomial, MonomialVector}) = m.vars From 0cd4cfdb891d883093d378d8ff0a2f94613aeeed Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Tue, 9 May 2023 18:22:50 +0200 Subject: [PATCH 05/23] fix --- src/monovec.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/monovec.jl b/src/monovec.jl index 2e97da8..4df106c 100644 --- a/src/monovec.jl +++ b/src/monovec.jl @@ -91,7 +91,7 @@ for (fun, call, def, ret) in [ end # faster complex-related functions MP.iscomplex(x::MonomialVector) = any(iscomplex, x.vars) -MP.conj(x::MonomialVector) = MonomialVector(conj.(x.vars), x.Z) # TODO: do we want to alias Z or copy Z? +Base.conj(x::MonomialVector) = MonomialVector(conj.(x.vars), x.Z) # TODO: do we want to alias Z or copy Z? _vars(m::Union{Monomial, MonomialVector}) = m.vars From 098ee294c94b9e0130991de8530cb097caeb73c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 May 2023 13:41:17 +0200 Subject: [PATCH 06/23] coefficienttype -> coefficient_type --- .github/workflows/ci.yml | 2 +- src/mult.jl | 2 +- src/subs.jl | 2 +- test/poly.jl | 6 +++--- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 141f953..c85f1c4 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -52,6 +52,6 @@ jobs: with: depwarn: error - uses: julia-actions/julia-processcoverage@v1 - - uses: codecov/codecov-action@v1 + - uses: codecov/codecov-action@v3 with: file: lcov.info diff --git a/src/mult.jl b/src/mult.jl index 3f97e0d..5718066 100644 --- a/src/mult.jl +++ b/src/mult.jl @@ -105,7 +105,7 @@ function Base.:(*)(p::Polynomial{true, S}, q::Polynomial{true, T}) where {S, T} if iszero(p) || iszero(q) zero(PT) else - polynomialclean(_mul(MP.coefficienttype(PT), p, q)...) + polynomialclean(_mul(MP.coefficient_type(PT), p, q)...) end end function MA.operate_to!(p::Polynomial{false, T}, ::typeof(*), q1::MP.AbstractPolynomialLike, q2::MP.AbstractPolynomialLike) where T diff --git a/src/subs.jl b/src/subs.jl index be74738..cc890c8 100644 --- a/src/subs.jl +++ b/src/subs.jl @@ -93,7 +93,7 @@ function _subs(::MP.Eval, p::Polynomial{C, T}, vals::AbstractVector{S}) where {C end end function _subs(::MP.Subs, p::Polynomial{C, T}, vals::AbstractVector{S}) where {C, T, S} - Tout = MA.promote_operation(*, T, MP.coefficienttype(S)) + Tout = MA.promote_operation(*, T, MP.coefficient_type(S)) q = zero_with_variables(Polynomial{C, Tout}, mergevars_of(PolyVar{C}, vals)[1]) for i in 1:length(p.a) MA.operate!(+, q, p.a[i] * monoeval(p.x.Z[i], vals)) diff --git a/test/poly.jl b/test/poly.jl index adbeea4..974b816 100644 --- a/test/poly.jl +++ b/test/poly.jl @@ -1,9 +1,9 @@ @testset "Term and Polynomial tests" begin @testset "Term" begin @polyvar x -# @test coefficienttype(1x) == Int -# @test coefficienttype(1.0x^2) == Float64 -# @test coefficienttype(Term{true, Int}) == Int +# @test coefficient_type(1x) == Int +# @test coefficient_type(1.0x^2) == Float64 +# @test coefficient_type(Term{true, Int}) == Int @test zero_term(Term{false, Int}).α == 0 @test one(Term{true, Int}).α == 1 @polyvar x From 43a0a619a89999c91e47910ceec0e2b8ec14e4d9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 10 May 2023 16:20:25 +0200 Subject: [PATCH 07/23] Remove incorrect method --- src/DynamicPolynomials.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/src/DynamicPolynomials.jl b/src/DynamicPolynomials.jl index 9652863..4f95cd9 100644 --- a/src/DynamicPolynomials.jl +++ b/src/DynamicPolynomials.jl @@ -35,7 +35,6 @@ MP.term_type(::Union{PolyType{C}, Type{<:PolyType{C}}}, ::Type{T}) where {C, T} MP.term_type(::Type{Polynomial{C}}) where {C} = Term{C} MP.polynomial_type(::Type{Term{C}}) where {C} = Polynomial{C} MP.polynomial_type(::Type{Term{C, T}}) where {T, C} = Polynomial{C, T} -MP.polynomial_type(::Type{T}, ::Type{<:DMonomialLike{C}}) where {T, C} = Polynomial{C, T} MP.polynomial_type(::Union{PolyType{C}, Type{<:PolyType{C}}}, ::Type{T}) where {C, T} = Polynomial{C, T} _vars(p::AbstractArray{<:PolyType}) = mergevars(_vars.(p))[1] MP.variables(p::Union{PolyType, MonomialVector, AbstractArray{<:PolyType}}) = _vars(p) # tuple(_vars(p)) From 214d169110af4ed2b600bea41dec8a29aae0629b Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Sun, 28 May 2023 21:31:22 +0200 Subject: [PATCH 08/23] Add a (sometimes working) implementation for incomplete substitution --- src/subs.jl | 45 ++++++++++++++++++++++++++++++++++++++------- 1 file changed, 38 insertions(+), 7 deletions(-) diff --git a/src/subs.jl b/src/subs.jl index 918bcb6..dab002e 100644 --- a/src/subs.jl +++ b/src/subs.jl @@ -16,7 +16,15 @@ function fillmap!( # 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. - if !iscomplex(s.first) + # 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 @@ -25,7 +33,8 @@ function fillmap!( end else for j in eachindex(vars) - if vars[j].id == s.first.id + 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 @@ -37,12 +46,34 @@ function fillmap!( else vals[j] = imag(value) end - elseif s.first.kind == vars[j].kind - vals[j] = real(s.second) # just to make sure the type is correct + 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 - # else: # assignment to parts of a complex variable, but the full appears in the polynomial - # We don't support this for the moment (would require more complex logic, since two substitutions may give the - # full value, or one part is left over - but can vals hold such a leftover?). end end end From c5c1c99e330466beb59706a175b78a2192aa3b71 Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Wed, 31 May 2023 18:16:18 +0200 Subject: [PATCH 09/23] Fix old syntax --- src/mono.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/mono.jl b/src/mono.jl index 4fd826e..2cd2d1b 100644 --- a/src/mono.jl +++ b/src/mono.jl @@ -102,8 +102,8 @@ MP.monomial(m::Monomial) = m #end # for efficiency reasons -function Base.conj(x::Monomial{true}) +function Base.conj(x::Monomial{V,M}) where {V<:Commutative,M} cv = conj.(x.vars) perm = sortperm(cv, rev = true) - return Monomial{true}(cv[perm], x.z[perm]) + return Monomial{V,M}(cv[perm], x.z[perm]) end From 31a8d4ce98581359fe60f079f0e7d9b5d5490536 Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Fri, 7 Jul 2023 10:47:33 +0200 Subject: [PATCH 10/23] Format --- test/poly.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/test/poly.jl b/test/poly.jl index 9cd4f28..a5553aa 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] From a19c0b03a904639894d31380519ff7d3ac47ec3b Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Thu, 17 Aug 2023 18:00:41 +0200 Subject: [PATCH 11/23] Remove TODO --- src/monomial_vector.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/monomial_vector.jl b/src/monomial_vector.jl index 98a07fb..562d284 100644 --- a/src/monomial_vector.jl +++ b/src/monomial_vector.jl @@ -113,7 +113,7 @@ for (fun, call, def, ret) in [ end # faster complex-related functions MultivariatePolynomials.iscomplex(x::MonomialVector) = any(iscomplex, x.vars) -Base.conj(x::MonomialVector) = MonomialVector(conj.(x.vars), x.Z) # TODO: do we want to alias Z or copy Z? +Base.conj(x::MonomialVector) = MonomialVector(conj.(x.vars), x.Z) MP.variables(m::Union{Monomial,MonomialVector}) = m.vars From 6b055bd2e981a16662871c3999cfc021fbc4c5de Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Thu, 24 Aug 2023 17:15:49 +0200 Subject: [PATCH 12/23] Add complex kind parameter to VT("name") constructor --- src/var.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/var.jl b/src/var.jl index 2eeba67..317ff42 100644 --- a/src/var.jl +++ b/src/var.jl @@ -150,9 +150,13 @@ struct Variable{V,M} <: AbstractVariable 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}, From a61aaa2a98d86b50f65cefed9e38abce3c7ff854 Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Fri, 1 Sep 2023 18:03:53 +0200 Subject: [PATCH 13/23] Remove iscomplex --- src/DynamicPolynomials.jl | 2 +- src/monomial_vector.jl | 8 ++++++-- src/var.jl | 2 +- 3 files changed, 8 insertions(+), 4 deletions(-) diff --git a/src/DynamicPolynomials.jl b/src/DynamicPolynomials.jl index f698120..453e192 100644 --- a/src/DynamicPolynomials.jl +++ b/src/DynamicPolynomials.jl @@ -72,7 +72,7 @@ function MP.similar_variable( return MP.similar_variable(P, S) end function MP.similar_variable(p::PolyType{V,M}, s::Symbol) where {V,M} - return Variable(string(s), V, M, iscomplex(p) ? cpFull : cpNone) + 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, diff --git a/src/monomial_vector.jl b/src/monomial_vector.jl index a810990..309e0d1 100644 --- a/src/monomial_vector.jl +++ b/src/monomial_vector.jl @@ -112,7 +112,7 @@ for (fun, call, def, ret) in [ end) end # faster complex-related functions -MultivariatePolynomials.iscomplex(x::MonomialVector) = any(iscomplex, x.vars) +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 @@ -144,7 +144,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/var.jl b/src/var.jl index 317ff42..1395c19 100644 --- a/src/var.jl +++ b/src/var.jl @@ -195,7 +195,7 @@ MP.variables(v::Variable) = [v] iscomm(::Type{Variable{C}}) where {C} = C -MP.iscomplex(x::Variable{C}) where {C} = x.kind == cpFull || x.kind == cpConj +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 From f87fb6226ea173851183a1bb62081712d79eac2c Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Thu, 30 Nov 2023 11:06:37 +0100 Subject: [PATCH 14/23] Remove dev dependency of MP, 0.5.3 has landed --- .github/workflows/ci.yml | 5 ----- Project.toml | 2 +- 2 files changed, 1 insertion(+), 6 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 09ffa12..4b564d9 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -42,11 +42,6 @@ jobs: ${{ runner.os }}-test-${{ env.cache-name }}- ${{ runner.os }}-test- ${{ runner.os }}- - - name: MP#complex - shell: julia --project=@. {0} - run: | - using Pkg - Pkg.add(PackageSpec(url="https://github.com/projekter/MultivariatePolynomials.jl.git", rev="master")) - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 with: diff --git a/Project.toml b/Project.toml index 622dac7..e780dbc 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" From da7fd62516e68974a2d2a97e26648af0d045d2f6 Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Wed, 21 Feb 2024 09:33:58 +0100 Subject: [PATCH 15/23] Add ComplexKind docstring --- src/var.jl | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/var.jl b/src/var.jl index 8ab1b8b..31012fb 100644 --- a/src/var.jl +++ b/src/var.jl @@ -143,6 +143,23 @@ 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: +- `cpNone`: the variable is real-valued, [`conj`](@ref) and [`real`](@ref) are identities, [`imag`](@ref) is zero. +- `cpFull`: the variable is a declared complex-valued variable with distinct conjugate, real, and imaginary part. + [`ordinary_variable`](@ref) is an identity. +- `cpConj`: 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. +- `cpReal`: 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. +- `cpImag`: 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 cpNone cpFull cpConj cpReal cpImag struct Variable{V,M} <: AbstractVariable From f8d3214d30c86f786311c515ce537bbd415ffec4 Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Wed, 21 Feb 2024 11:36:14 +0100 Subject: [PATCH 16/23] Rename ComplexKind instances --- src/DynamicPolynomials.jl | 4 +-- src/subs.jl | 28 +++++++++---------- src/var.jl | 58 +++++++++++++++++++-------------------- 3 files changed, 45 insertions(+), 45 deletions(-) diff --git a/src/DynamicPolynomials.jl b/src/DynamicPolynomials.jl index 43a03e5..4f0d690 100644 --- a/src/DynamicPolynomials.jl +++ b/src/DynamicPolynomials.jl @@ -73,10 +73,10 @@ function MP.similar_variable( return MP.similar_variable(P, S) end function MP.similar_variable(p::PolyType{V,M}, s::Symbol) where {V,M} - return Variable(string(s), V, M, isreal(p) ? cpNone : cpFull) + 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, cpNone) # we cannot infer this from the type, + return Variable(string(s), V, M, REAL) # we cannot infer this from the type, end include("promote.jl") diff --git a/src/subs.jl b/src/subs.jl index 5c9de9d..6ee73bf 100644 --- a/src/subs.jl +++ b/src/subs.jl @@ -24,7 +24,7 @@ function fillmap!( # - 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 + if s.first.kind == REAL for j in eachindex(vars) if vars[j] == s.first vals[j] = s.second @@ -35,41 +35,41 @@ function fillmap!( 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 + if s.first.kind == COMPLEX || s.first.kind == CONJ + value = s.first.kind == CONJ ? conj(s.second) : s.second + if vars[j].kind == COMPLEX vals[j] = value - elseif vars[j].kind == cpConj + elseif vars[j].kind == CONJ vals[j] = conj(value) - elseif vars[j].kind == cpReal + elseif vars[j].kind == REAL_PART vals[j] = real(value) else vals[j] = imag(value) end - elseif s.first.kind == cpReal + elseif s.first.kind == REAL_PART 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 + if vars[j].kind == COMPLEX vals[j] = value + im * imag(vals[j]) - elseif vars[j].kind == cpConj + elseif vars[j].kind == CONJ vals[j] = value - im * imag(vals[j]) - elseif vars[j].kind == cpReal + elseif vars[j].kind == REAL_PART vals[j] = value end # else we know the real part but use the imaginary part; do nothing else - @assert(s.first.kind == cpImag) + @assert(s.first.kind == IMAG_PART) 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 + if vars[j].kind == COMPLEX vals[j] = real(vals[j]) + im * value - elseif vars[j].kind == cpConj + elseif vars[j].kind == CONJ vals[j] = real(vals[j]) - im * value - elseif vars[j].kind == cpImag + elseif vars[j].kind == IMAG_PART vals[j] = value end # else we know the imaginary part but use the real part; do nothing diff --git a/src/var.jl b/src/var.jl index 31012fb..9245f63 100644 --- a/src/var.jl +++ b/src/var.jl @@ -73,7 +73,7 @@ function _extract_kw_args(args, variable_order, complex_kind) 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 + complex_kind = arg.args[2] == :true ? COMPLEX : REAL else error("Unrecognized keyword argument `$(arg.args[1])`") end @@ -87,7 +87,7 @@ 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, complex_kind = - _extract_kw_args(args, :($(Commutative{CreationOrder})), cpNone) + _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 = :())); @@ -96,7 +96,7 @@ end macro ncpolyvar(args...) pos_args, variable_order, monomial_order, complex_kind = - _extract_kw_args(args, :($(NonCommutative{CreationOrder})), cpNone) + _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 = :())); @@ -105,7 +105,7 @@ end macro polycvar(args...) pos_args, variable_order, monomial_order, complex_kind = - _extract_kw_args(args, :($(Commutative{CreationOrder})), cpFull) + _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 = :())); @@ -147,20 +147,20 @@ end ComplexKind The type used to identify whether a variable is complex-valued. It has the following instances: -- `cpNone`: the variable is real-valued, [`conj`](@ref) and [`real`](@ref) are identities, [`imag`](@ref) is zero. -- `cpFull`: the variable is a declared complex-valued variable with distinct conjugate, real, and imaginary part. +- `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. -- `cpConj`: the variable is the conjugate of a declared complex-valued variable. [`ordinary_variable`](@ref) is equal to +- `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. -- `cpReal`: the variable is the real part of a complex-valued variable. It is real-valued itself: [`conj`](@ref) and +- `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. -- `cpImag`: the variable is the imaginary part of a declared complex-valued variable (or the negative imaginary part of the +- `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 cpNone cpFull cpConj cpReal cpImag +@enum ComplexKind REAL COMPLEX CONJ REAL_PART IMAG_PART struct Variable{V,M} <: AbstractVariable name::String @@ -169,7 +169,7 @@ struct Variable{V,M} <: AbstractVariable function Variable{V,M}( name::AbstractString, - kind::ComplexKind = cpNone, + kind::ComplexKind = REAL, ) where {V<:AbstractVariableOrdering,M<:MP.AbstractMonomialOrdering} return new{V,M}(convert(String, name), kind, instantiate(V)) end @@ -178,7 +178,7 @@ struct Variable{V,M} <: AbstractVariable name::AbstractString, ::Type{V}, ::Type{M}, - kind::ComplexKind = cpNone, + kind::ComplexKind = REAL, ) where {V<:AbstractVariableOrdering,M<:MP.AbstractMonomialOrdering} return new{V,M}(convert(String, name), kind, instantiate(V)) end @@ -187,8 +187,8 @@ struct Variable{V,M} <: AbstractVariable from::Variable{V,M}, kind::ComplexKind, ) where {V<:AbstractVariableOrdering,M<:MP.AbstractMonomialOrdering} - (from.kind == cpNone && kind == cpNone) || - (from.kind != cpNone && kind != cpNone) || + (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 @@ -214,36 +214,36 @@ 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 +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 == cpNone || x.kind == cpFull ? x : Variable(x, cpFull) + return x.kind == REAL || x.kind == COMPLEX ? x : Variable(x, COMPLEX) end function Base.conj(x::Variable) - if x.kind == cpFull - return Variable(x, cpConj) - elseif x.kind == cpConj - return Variable(x, cpFull) + 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 == cpFull || x.kind == cpConj - return Variable(x, cpReal) + 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 == cpFull - return Variable(x, cpImag) - elseif x.kind == cpConj - return _Term{V,M,Int}(-1, Monomial(Variable(x, cpImag))) + 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 From 7e075850425c3de9193c73ebd2f0159964d5cd57 Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Sat, 4 May 2024 10:50:42 +0200 Subject: [PATCH 17/23] typo --- src/DynamicPolynomials.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/DynamicPolynomials.jl b/src/DynamicPolynomials.jl index 4f0d690..35c56a0 100644 --- a/src/DynamicPolynomials.jl +++ b/src/DynamicPolynomials.jl @@ -76,7 +76,7 @@ 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, + return Variable(string(s), V, M, REAL) # we cannot infer this from the type end include("promote.jl") From 8cf9c1215bf5e4c7641314769df217223737f449 Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Sat, 4 May 2024 10:50:55 +0200 Subject: [PATCH 18/23] change substitution rules --- src/subs.jl | 36 +++++++++++++----------------------- 1 file changed, 13 insertions(+), 23 deletions(-) diff --git a/src/subs.jl b/src/subs.jl index 6ee73bf..1792097 100644 --- a/src/subs.jl +++ b/src/subs.jl @@ -13,17 +13,15 @@ function fillmap!( s::MP.Substitution, ) 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. + # 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 @@ -51,28 +49,20 @@ function fillmap!( "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 == COMPLEX - vals[j] = value + im * imag(vals[j]) - elseif vars[j].kind == CONJ - vals[j] = value - im * imag(vals[j]) - elseif vars[j].kind == REAL_PART + if vars[j].kind == REAL_PART vals[j] = value end - # else we know the real part but use the imaginary part; do nothing + # don't do a partial substitution else @assert(s.first.kind == IMAG_PART) 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 == COMPLEX - vals[j] = real(vals[j]) + im * value - elseif vars[j].kind == CONJ - vals[j] = real(vals[j]) - im * value - elseif vars[j].kind == IMAG_PART + if vars[j].kind == IMAG_PART vals[j] = value end - # else we know the imaginary part but use the real part; do nothing + # don't do a partial substitution end end end From 6c230f65dd97b7af6d2f4477c8ad946e3d527e99 Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Sat, 4 May 2024 10:52:17 +0200 Subject: [PATCH 19/23] Rename polycvar macro --- src/var.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/var.jl b/src/var.jl index 9245f63..3e4fc66 100644 --- a/src/var.jl +++ b/src/var.jl @@ -1,4 +1,4 @@ -export Variable, @polyvar, @ncpolyvar, @polycvar +export Variable, @polyvar, @ncpolyvar, @complex_polyvar function polyarrayvar( variable_order, @@ -103,7 +103,7 @@ macro ncpolyvar(args...) $(Expr(:tuple, esc.(vars)...))) end -macro polycvar(args...) +macro complex_polyvar(args...) pos_args, variable_order, monomial_order, complex_kind = _extract_kw_args(args, :($(Commutative{CreationOrder})), COMPLEX) vars, exprs = From a6e1d801ebd4e5156d23653da4d3af6623cefb01 Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Wed, 8 May 2024 18:36:51 +0200 Subject: [PATCH 20/23] Error for partial substitutions --- src/subs.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/subs.jl b/src/subs.jl index 1792097..a7a5e62 100644 --- a/src/subs.jl +++ b/src/subs.jl @@ -51,8 +51,9 @@ function fillmap!( value = real(s.second) # just to make sure the type is correct if vars[j].kind == REAL_PART vals[j] = value + elseif vars[j].kind != IMAG_PART + error("Found complex variable with substitution of real part - not implemented") end - # don't do a partial substitution else @assert(s.first.kind == IMAG_PART) isreal(s.second) || error( @@ -61,8 +62,9 @@ function fillmap!( value = real(s.second) # just to make sure the type is correct if vars[j].kind == IMAG_PART vals[j] = value + elseif vars[j].kind != REAL_PART + error("Found complex variable with substitution of imaginary part - not implemented") end - # don't do a partial substitution end end end From f9337ffeea0a93df76334d9779efdb01365edcfc Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Wed, 8 May 2024 18:39:31 +0200 Subject: [PATCH 21/23] Format --- src/subs.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/subs.jl b/src/subs.jl index a7a5e62..6c0c9ec 100644 --- a/src/subs.jl +++ b/src/subs.jl @@ -52,7 +52,9 @@ function fillmap!( if vars[j].kind == REAL_PART vals[j] = value elseif vars[j].kind != IMAG_PART - error("Found complex variable with substitution of real part - not implemented") + error( + "Found complex variable with substitution of real part - not implemented", + ) end else @assert(s.first.kind == IMAG_PART) @@ -63,7 +65,9 @@ function fillmap!( if vars[j].kind == IMAG_PART vals[j] = value elseif vars[j].kind != REAL_PART - error("Found complex variable with substitution of imaginary part - not implemented") + error( + "Found complex variable with substitution of imaginary part - not implemented", + ) end end end From e0ad44dd6409f9276ca701ba91b476bfcb61a9c8 Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Wed, 8 May 2024 20:48:21 +0200 Subject: [PATCH 22/23] Complex substitution test cases --- test/poly.jl | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/test/poly.jl b/test/poly.jl index a5553aa..50c736a 100644 --- a/test/poly.jl +++ b/test/poly.jl @@ -108,5 +108,26 @@ ) @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 + @test pc(conj(z) => 2 - 3im) == 798 + 9im + @test real(pc)(z => 2 + 3im) == 798 + err = ArgumentError( + "Variable `zᵢ` was not assigned a value. Use `subs` to substitute only a subset of the variables.", + ) + @test_throws err real(pc)(real(z) => 2) + @test subs(real(pc), real(z) => 2) == + 12 - 224imag(z) - 6imag(z)^2 + 56imag(z)^3 + @test real(pc)([real(z), imag(z)] => [2, 3]) == 798 + err = ErrorException( + "Found complex variable with substitution of real part - not implemented", + ) + @test_throws err subs(pc, real(z) => 2) + err = ErrorException( + "Found complex variable with substitution of imaginary part - not implemented", + ) + @test_throws err subs(pc, imag(z) => 3) end end From 0845707acb299d79e970d4cd2e621dbbc04e80b4 Mon Sep 17 00:00:00 2001 From: Benjamin Desef Date: Thu, 9 May 2024 13:05:18 +0200 Subject: [PATCH 23/23] More strict substitution rules: only allow ordinary variable --- src/subs.jl | 47 ++++++++++++----------------------------------- test/poly.jl | 20 +++++--------------- 2 files changed, 17 insertions(+), 50 deletions(-) diff --git a/src/subs.jl b/src/subs.jl index 6c0c9ec..f856522 100644 --- a/src/subs.jl +++ b/src/subs.jl @@ -30,45 +30,22 @@ function fillmap!( 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 s.first.kind == COMPLEX || s.first.kind == CONJ - value = s.first.kind == CONJ ? conj(s.second) : s.second - if vars[j].kind == COMPLEX - vals[j] = value - elseif vars[j].kind == CONJ - vals[j] = conj(value) - elseif vars[j].kind == REAL_PART - vals[j] = real(value) - else - vals[j] = imag(value) - end - elseif s.first.kind == REAL_PART - 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 == REAL_PART - vals[j] = value - elseif vars[j].kind != IMAG_PART - error( - "Found complex variable with substitution of real part - not implemented", - ) - end + 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 - @assert(s.first.kind == IMAG_PART) - 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 == IMAG_PART - vals[j] = value - elseif vars[j].kind != REAL_PART - error( - "Found complex variable with substitution of imaginary part - not implemented", - ) - end + vals[j] = imag(s.second) end end end diff --git a/test/poly.jl b/test/poly.jl index 50c736a..a02f05d 100644 --- a/test/poly.jl +++ b/test/poly.jl @@ -112,22 +112,12 @@ @complex_polyvar z pc = z^3 + 2real(z) - 7imag(z^4) @test pc(z => 2 + 3im) == 798 + 9im - @test pc(conj(z) => 2 - 3im) == 798 + 9im - @test real(pc)(z => 2 + 3im) == 798 err = ArgumentError( - "Variable `zᵢ` was not assigned a value. Use `subs` to substitute only a subset of the variables.", - ) - @test_throws err real(pc)(real(z) => 2) - @test subs(real(pc), real(z) => 2) == - 12 - 224imag(z) - 6imag(z)^2 + 56imag(z)^3 - @test real(pc)([real(z), imag(z)] => [2, 3]) == 798 - err = ErrorException( - "Found complex variable with substitution of real part - not implemented", + "Substitution with complex variables requires the ordinary_variable in the substitution specification" ) - @test_throws err subs(pc, real(z) => 2) - err = ErrorException( - "Found complex variable with substitution of imaginary part - not implemented", - ) - @test_throws err subs(pc, imag(z) => 3) + @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