JuliaSymbolics / Symbolics.jl

Symbolic programming for the next generation of numerical software
https://docs.sciml.ai/Symbolics/stable/
Other
1.36k stars 152 forks source link

Is it possible to get rid of `Complex{Num}`? #884

Open lairez opened 1 year ago

lairez commented 1 year ago

The insistance of Symbolics.jl to separate the real and the imaginary parts of complex numbers is causing a host of issues when we expect complex numbers to be mere numbers. (Being a bit provocative, you could declare Complex <: Real, that should not cause too much problems...)

Two examples.

First, a substitution that does not work as expected.

julia> @variables X
1-element Vector{Num}:
 X

julia> substitute([X], Dict([X=>1]))
1-element Vector{Num}:
 1.0

julia> substitute([im*X], Dict([X=>1]))
1-element Vector{Complex{Num}}:
 X*im

Second, a example that shows that separating real and imaginary parts may not be the good way for computations, so Symbolics.jl should not impose it.

julia> @variables X::Complex
1-element Vector{Complex{Num}}:
 X

julia> 1/(1-X^10)
(1 + 16((real(X)^2 - (imag(X)^2))^2 - 4(imag(X)^2)*(real(X)^2))*(real(X)^2 - (imag(X)^2))*(imag(X)^2)*(real(X)^2) - (((real(X)^2 - (imag(X)^2))^2 - 4(imag(X)^2)*(real(X)^2))^2 - 16((real(X)^2 - (imag(X)^2))^2)*(imag(X)^2)*(real(X)^2))*(real(X)^2 - (imag(X)^2))) / ((-2(((real(X)^2 - (imag(X)^2))^2 - 4(imag(X)^2)*(real(X)^2))^2 - 16((real(X)^2 - (imag(X)^2))^2)*(imag(X)^2)*(real(X)^2))*imag(X)*real(X) - 8((real(X)^2 - (imag(X)^2))^2 - 4(imag(X)^2)*(real(X)^2))*((real(X)^2 - (imag(X)^2))^2)*imag(X)*real(X))^2 + (1 + 16((real(X)^2 - (imag(X)^2))^2 - 4(imag(X)^2)*(real(X)^2))*(real(X)^2 - (imag(X)^2))*(imag(X)^2)*(real(X)^2) - (((real(X)^2 - (imag(X)^2))^2 - 4(imag(X)^2)*(real(X)^2))^2 - 16((real(X)^2 - (imag(X)^2))^2)*(imag(X)^2)*(real(X)^2))*(real(X)^2 - (imag(X)^2)))^2) + (im*(2(((real(X)^2 - (imag(X)^2))^2 - 4(imag(X)^2)*(real(X)^2))^2 - 16((real(X)^2 - (imag(X)^2))^2)*(imag(X)^2)*(real(X)^2))*imag(X)*real(X) + 8((real(X)^2 - (imag(X)^2))^2 - 4(imag(X)^2)*(real(X)^2))*((real(X)^2 - (imag(X)^2))^2)*imag(X)*real(X))) / ((-2(((real(X)^2 - (imag(X)^2))^2 - 4(imag(X)^2)*(real(X)^2))^2 - 16((real(X)^2 - (imag(X)^2))^2)*(imag(X)^2)*(real(X)^2))*imag(X)*real(X) - 8((real(X)^2 - (imag(X)^2))^2 - 4(imag(X)^2)*(real(X)^2))*((real(X)^2 - (imag(X)^2))^2)*imag(X)*real(X))^2 + (1 + 16((real(X)^2 - (imag(X)^2))^2 - 4(imag(X)^2)*(real(X)^2))*(real(X)^2 - (imag(X)^2))*(imag(X)^2)*(real(X)^2) - (((real(X)^2 - (imag(X)^2))^2 - 4(imag(X)^2)*(real(X)^2))^2 - 16((real(X)^2 - (imag(X)^2))^2)*(imag(X)^2)*(real(X)^2))*(real(X)^2 - (imag(X)^2)))^2)

Is there any way to avoid this behavior and have complex numbers behave the same hay as real numbers?

aravindh-krishnamoorthy commented 1 year ago

Personally, I think Complex{Num} is a cool idea and in fact, if I had to develop my own symbolic library, I'd want to offload as much as possible to the "real libraries."

The great thing about Complex{Num} is that it is Symbolics' Num wrapped in Base's Complex. I feel this is an elegant way to implement complex numbers. But, I'm not aware of sufficient design details of Symbolics.jl to suggest how to continue with complex numbers.

To say the truth, I also couldn't work out how to use substitute with Complex{Num} but also realised that substitute is mostly offloaded to SymbolicUtils.jl. Anyway, at least the following works:

julia> @variables z::Complex
1-element Vector{Complex{Num}}:
 z

julia> Symbolics.substitute(z, Dict(z => 1+im))
1 + im

But the following is broken:

julia> Symbolics.substitute(2*z, Dict(z => 1+im))
ERROR: type ComplexTerm has no field metadata
Stacktrace:
 [1] getproperty
   @ ./Base.jl:38 [inlined]
 [2] metadata(s::Symbolics.ComplexTerm{Real})
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/H684H/src/types.jl:599
 [3] substitute(expr::Symbolics.ComplexTerm{Real}, dict::Dict{Symbolics.ComplexTerm{Real}, Complex{Int64}}; fold::Bool)
   @ SymbolicUtils ~/.julia/packages/SymbolicUtils/H684H/src/substitute.jl:34
 [4] substitute
   @ ~/.julia/packages/SymbolicUtils/H684H/src/substitute.jl:16 [inlined]
 [5] #16#17
   @ ~/.julia/dev/Symbolics/src/num.jl:87 [inlined]
 [6] (::Symbolics.var"#16#19"{Symbolics.var"#16#17#20"{Dict{Symbolics.ComplexTerm{Real}, Complex{Int64}}}})(expr::Symbolics.ComplexTerm{Real})
   @ Symbolics ~/.julia/dev/Symbolics/src/num.jl:87
 [7] substitute(expr::Complex{Num}, s::Dict{Complex{Num}, Complex{Int64}}; kw::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ Symbolics ~/.julia/dev/Symbolics/src/Symbolics.jl:156
 [8] substitute(expr::Complex{Num}, s::Dict{Complex{Num}, Complex{Int64}})
   @ Symbolics ~/.julia/dev/Symbolics/src/Symbolics.jl:156
 [9] top-level scope
   @ REPL[6]:1

See also #883

amilsted commented 1 year ago

Personally, I think Complex{Num} is a cool idea and in fact, if I had to develop my own symbolic library, I'd want to offload as much as possible to the "real libraries."

I think it's problematic to separate real and imaginary parts by default - there is specialized code for handling numerical operations on complex numbers natively. Where this exists, I want to use it! It would be nice if splitting into real(z) + im*imag(z) were optional.

ChrisRackauckas commented 1 year ago

@shashi comments?

xtalax commented 1 year ago

Piggybacking this thread to say that the ComplexTerm bug would be fixed by #789 and this being unmerged is the last remaining blocker that I know of to solving complex PDEs like the Schrödinger equation in MOL

m0Cey commented 1 year ago

Hello everyone! I'm new to GitHub, so I didn't quite understand if the pull request (is this a right name?) from the last comment already have effect on the package, but this piece of code:

@variables a::Complex{Num}
a_c = conj(a)
substitute(a_c, Dict([a => 1im]))

gives the same error to OP

ERROR: type ComplexTerm has no field metadata

So, can you clarify for me please if it's still WIP fix or am I doing something wrong? P. S. The output generated from Symbolics v5.5.1.

shashi commented 1 year ago

I merged that PR which fixes this "no field metadata" issue.

I think Complex{Num} is here to stay. If there was an AbstractComplex type in Base we would definitely be doing something like what we do for Reals. It's probably an issue that first needs to go into the host language and libraries using Complex already.

corbett5 commented 1 year ago

After updating to v5.5.3 such that I have the fix from #789, performing substitute on a complex expression no longer fails, but it does not work as expected:

using Symbolics;
@variables z::Complex{Num};
substitute(2 * z, Dict(z => 1))

The output is

2real(z) + 2im*imag(z)
Sollovin commented 7 months ago

I think I also met the problem when I'm trying to obtain the coefficient matrix of a series of complex linear equations.

For example,

using Symbolics
@variables x1 x2 x3
@variables xc1::Complex xc2::Complex xc3::Complex

X = vcat(x1, x2, x3)
Xc = vcat(xc1, xc2, xc3)
A = [[1 2 3]; [3.2 2.5 3.3]; [1.1 1.4 2.9]]
Ac = [[1+1im 2 3]; [3.2 2.5 3.3]; [1.1 1.4 2.9]]
Y = vcat(5, 6, 7)

Symbolics.A_b(A*X - Y .~ 0, X, true) # Works.

Symbolics.A_b(Ac*Xc - Y .~ 0, Xc, true) # Do not work.

Because Ac*Xc - Y .~ 0 gives out

3-element Vector{Vector{Equation}}:
 [-5 - imag(xc1) + real(xc1) + 3.0real(xc3) + 2.0real(xc2) ~ 0, imag(xc1) + 3.0imag(xc3) + real(xc1) + 2.0imag(xc2) ~ 0]
 [-6 + 3.2real(xc1) + 3.3real(xc3) + 2.5real(xc2) ~ 0, 3.2imag(xc1) + 3.3imag(xc3) + 2.5imag(xc2) ~ 0]
 [-7 + 1.1real(xc1) + 2.9real(xc3) + 1.4real(xc2) ~ 0, 1.1imag(xc1) + 2.9imag(xc3) + 1.4imag(xc2) ~ 0]