jump-dev / JuMP.jl

Modeling language for Mathematical Optimization (linear, mixed-integer, conic, semidefinite, nonlinear)
http://jump.dev/JuMP.jl/
Other
2.17k stars 390 forks source link

Add a method for Complex(a,b) #3733

Closed araujoms closed 2 months ago

araujoms commented 2 months ago

Could you please add a method Complex(a,b) = a + im*b for real JuMP variables a and b?

The reason is that some time ago I generalized the function __svec_to_smat_complex! from https://github.com/jump-dev/Hypatia.jl/blob/master/src/Cones/arrayutilities.jl , to work with JuMP variables as well by replacing Complex(a,b) with a + im*b. But as it turns out this slightly degrades performance for numerical types, so I'd like to undo my hack.

odow commented 2 months ago

I think defining Complex(a, b) to return something other than ::Complex is legal, but bad form in Julia, so I would prefer not to do this.

Do you mean complex(r, i)?

help?> complex
search: complex Complex ComplexF64 ComplexF32 ComplexF16 ComplexPlane ComplexVariable precompile __precompile__ Highs_compilationDate

  complex(r, [i])

  Convert real numbers or arrays to complex. i defaults to zero.

  Examples
  ≡≡≡≡≡≡≡≡

  julia> complex(7)
  7 + 0im

  julia> complex([1, 2, 3])
  3-element Vector{Complex{Int64}}:
   1 + 0im
   2 + 0im
   3 + 0im

  ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

  complex(T::Type)

  Return an appropriate type which can represent a value of type T as a complex number. Equivalent to typeof(complex(zero(T))).

  Examples
  ≡≡≡≡≡≡≡≡

  julia> complex(Complex{Int})
  Complex{Int64}

  julia> complex(Int)
  Complex{Int64}

julia> complex(1, 2)
1 + 2im

Is this really an improvement on a + im * b? Do you have a benchmark on a realistic instance? I don't really see how this could make a difference.

araujoms commented 2 months ago

I really do mean Complex(a,b). No, I don't have a benchmark, but this is a function used deep inside almost every cone in Hypatia, so every microsecond matters.

@chriscoey I believe you're the one who wrote the original function with Complex(a,b) instead of a + im*b, care to comment?

odow commented 2 months ago

Why does code in Hypatia need to work for JuMP variables?

In your own code, can't you just do:

my_complex(a, b) = Complex(a, b)
my_complex(a::AbstractJuMPScalar, b::AbstractJuMPScalar) = a + im * b

Complex(a,b) instead of a + im*b

a + im*b does

tmp = (im * b)::ComplexF64
(a + tmp)::ComplexF64

so it introduces an additional temporary.

araujoms commented 2 months ago

Because several of the examples use JuMP.

And no, I'm not defining a my_complex function in order to add a method to JuMP. That's hideous and makes the code harder to read. Moreover, Complex is a really basic function, I'm certain that such a method will also be useful for other people.

chriscoey commented 2 months ago

I did Complex(a, b) in my original code for the reason Oscar mentioned (avoid a temp), and personally I think it's just as readable. Hypatia doesn't use JuMP internally, but some of the Hypatia examples do in order to help build the MOI models.

odow commented 2 months ago

I think we can accept a PR for:

Base.complex(x::AbstractJuMPScalar, y::AbstractJuMPScalar) = x + im * y
Base.complex(x, y::AbstractJuMPScalar) = x + im * y
Base.complex(x::AbstractJuMPScalar, y) = x + im * y

but I'm not going to support Base.Complex. That's a constructor and we can't return an instance of the type.

araujoms commented 2 months ago

Great, that works for me, I'll do a PR. Shouldn't it be Union{GenericVariableRef{<:Real}, GenericAffExpr{<:Real}, GenericQuadExpr{<:Real}} instead of AbstractJuMPScalar, though?

odow commented 2 months ago

I guess so. I assumed complex would just do the equivalent, but it throws a MethodError, so I guess we should match.

julia> complex(1+2im, 3)
ERROR: MethodError: no method matching complex(::Complex{Int64}, ::Int64)

Closest candidates are:
  complex(::Real, ::Real)
   @ Base complex.jl:175
  complex(::Complex)
   @ Base complex.jl:173

Stacktrace:
 [1] top-level scope
   @ REPL[150]:1

julia> complex(1+2im, 0)
ERROR: MethodError: no method matching complex(::Complex{Int64}, ::Int64)

Closest candidates are:
  complex(::Real, ::Real)
   @ Base complex.jl:175
  complex(::Complex)
   @ Base complex.jl:173

Stacktrace:
 [1] top-level scope
   @ REPL[151]:1

julia> complex(1+2im, 0im)
ERROR: MethodError: no method matching complex(::Complex{Int64}, ::Complex{Int64})

Closest candidates are:
  complex(::Complex)
   @ Base complex.jl:173

Stacktrace:
 [1] top-level scope
   @ REPL[152]:1