Real and complex coupled types as parametric types

jarl · July 30, 2018, 1:35pm

I want a function that takes two parameters, one (potentially) complex and one real parameter. They should have the same precision / underlying format. I tried the following but get an error I do not understand:

julia> function f(z::T, myreal::real(T)) where{T<:Number}
           return z+myreal;
       end
ERROR: MethodError: no method matching real(::TypeVar)
Closest candidates are:
  real(::Complex) at complex.jl:54
  real(::Real) at complex.jl:67
  real(::Type{T<:Real}) where T<:Real at complex.jl:98
  ...

This is confusing to me, since I am allowed to use real(T) in a different way:

function f2(z::T, y=one(real(T))) where{T<:Number}
    return z+y;
end

The following code gives the behaviour I want, except that the operations are carried out in Complex even if z is Real:

f3(z::T,myreal::T) where {T<:Real} = f3(complex(z),myreal)
function f3(z::Complex{T}, myreal::T) where {T<:Real}
    return z+myreal;
end

julia version 0.6.4

platawiec · July 30, 2018, 1:57pm

Try:

julia> function f(z::Union{T, Complex{T}}, myreal::T) where{T<:Real}
           return typeof(z);
       end
f (generic function with 1 method)

julia> f(1.0, 0.0)
Float64

julia> f(1.0im, 0.0)
Complex{Float64}

stevengj · July 30, 2018, 2:20pm

This thread was marked “performance,” but declaring argument types has no impact on performance in Julia. The compiler specializes the function for the caller’s argument types even if no types were declared. The type declarations are only useful as a “filter” to say which functions are applicable to which types.

If you declare argument types at all, you should make them as generic as possible, e.g. just use Number, to make your function as flexible as possible (e.g. to work with dual numbers for automatic differentiation).

jarl · July 30, 2018, 3:02pm

Thanks. @platawiec’s solution does partially solves my problem. I fear my attempt to construct a mwe made it unclear where the performance aspect lies (or maybe I misunderstood something).

A bit more elaborate example: I want the first argument to determine which arithmetic my algorithm should run, and default to Complex128. I also have some kwargs.

So, a and b should be the same:

a=myfun(3.0)
b=myfun(Complex128,3.0)

This is how I tried to achieve this (drastically simplified from my application, dont pay so much attention to the contents of the function):

myfun(x::Float64;kwargs...)=myfun(Complex128,x;kwargs...)
function myfun(::Type{T_arithmetic},x::Float64;
               x1::Vector{<:Number}=ones(T_arithmetic,2)
               ) where T_arithmetic<:Number
    z=zeros(T_arithmetic,2)
    z[1]=x; 
    z+=x1; # no inplace 
    return z;
end

The call myfun(Complex128,3.0) is type stable according to (@code_warntype), but myfun(3.0) is not. If I replace <:Number with T_arithmetic in the kwarg it becomes type-stable. However, I want it to be real vector. If I replace <:Number with Vector{real(T_arithmetic)} I get the same error as above.

This type instability can cause a performance penalty, right?

Type stability check:

julia> @code_warntype(myfun(3.0))
Variables:
  #self#::#myfun
  x::Float64

Body:
  begin 
      return $(Expr(:invoke, MethodInstance for #myfun#1(::Array{Any,1}, ::Function, ::Float64), :(Main.#myfun#1), :($(Expr(:foreigncall, :(:jl_alloc_array_1d), Array{Any,1}, svec(Any, Int64), Array{Any,1}, 0, 0, 0))), :(#self#), :(x)))
  end::Any

julia> @code_warntype(myfun(Complex128,3.0))
Variables:
  #self# <optimized out>
  #temp# <optimized out>
  x::Float64
  z::Array{Complex{Float64},1}

Body:
  begin 
      SSAValue(1) = $(Expr(:invoke, MethodInstance for fill!(::Array{Complex{Float64},1}, ::Complex{Float64}), :(Base.fill!), :($(Expr(:foreigncall, :(:jl_alloc_array_1d), Array{Complex{Float64},1}, svec(Any, Int64), Array{Complex{Float64},1}, 0, 2, 0))), :($(Expr(:new, Complex{Float64}, :((Base.sitofp)(Float64, 1)::Float64), :((Base.sitofp)(Float64, 0)::Float64))))))
      $(Expr(:inbounds, false))
      # meta: location /home/jarl/.julia/v0.6/NonlinearEigenproblems/src/mynewton.jl #myfun#2 7
      z::Array{Complex{Float64},1} = $(Expr(:invoke, MethodInstance for fill!(::Array{Complex{Float64},1}, ::Complex{Float64}), :(Base.fill!), :($(Expr(:foreigncall, :(:jl_alloc_array_1d), Array{Complex{Float64},1}, svec(Any, Int64), Array{Complex{Float64},1}, 0, 2, 0))), :($(Expr(:new, Complex{Float64}, :((Base.sitofp)(Float64, 0)::Float64), :((Base.sitofp)(Float64, 0)::Float64)))))) # line 8:
      (Base.arrayset)(z::Array{Complex{Float64},1}, $(Expr(:new, Complex{Float64}, :(x), :((Base.sitofp)(Float64, 0)::Float64))), 1)::Array{Complex{Float64},1} # line 9:
      z::Array{Complex{Float64},1} = $(Expr(:invoke, MethodInstance for +(::Array{Complex{Float64},1}, ::Array{Complex{Float64},1}), :(Main.+), :(z), SSAValue(1)))
      # meta: pop location
      $(Expr(:inbounds, :pop))
      return z::Array{Complex{Float64},1}
  end::Array{Complex{Float64},1}

dpsanders · July 30, 2018, 3:33pm

Something like this? You need a conversion in there somewhere.

julia> myfun(x, y) = myfun(Complex128(x), y)
myfun (generic function with 1 method)

julia> myfun(x::Complex128, y) = x + y
myfun (generic function with 2 methods)

julia> myfun(1, 2)
3.0 + 0.0im

jarl · July 30, 2018, 4:06pm

Thanks. From my perspective, this is like my second example in the original post. When I run f(1,2) I would like to compute 3 with integer arithmetic. (In particular I don’t want 3.0+0.0im and also not 3+0im.)

DNF · July 30, 2018, 5:27pm

What about something like this:

function myfun(::Type{S}, x::T; x1::Vector{T}=zeros(T, 2)) where {S<:Union{T, Complex{T}}} where {T<:Real} 
    z = zeros(S, 2)
    z[1] = x
    z += x1  # no inplace 
    return z
end
myfun(x::Number; kwargs...) = myfun(Complex128, x; kwargs...)

julia> myfun(3.0)
Complex{Float64}
Array{Complex{Float64},1}
2-element Array{Complex{Float64},1}:
 3.0+0.0im
 0.0+0.0im

julia> myfun(Float64, 3.0)
Float64
Array{Float64,1}
2-element Array{Float64,1}:
 3.0
 0.0

?

dpsanders · July 30, 2018, 7:05pm

Then I didn’t understand what you meant by “wanting to use complex128 by default”.
What behaviour are you looking for? Please give example output.

jarl · July 30, 2018, 7:29pm

The response of @DNF solves the problem completely. I did not know one could use where twice. Awasome. Thanks everyone for all great quick answers.

jarl · July 30, 2018, 7:30pm

Thanks. I think the confusion arose because I had two minimal examples with inconsistent notation.

jarl · July 31, 2018, 8:53am

When I adapted @DNF’s solution to my application, I encountered another type error, which I would like to provide here for completeness. I managed to get around it in a quite unnatural way. It would be nice with some insight why this happens, or if there are more natural work-arounds.

If I remove the x::T in the parameter list (not kwarg-list) the type T seems not visible in the kwargs:

This works:

  function myfun1(::Type{S}, x::T; x1::Vector{T}=zeros(T, 2)) where {S<:Union{T, Complex{T}}} where {T<:Real} 
    z = zeros(S, 2)
    z += x1  
    return z
end

This throws an error:

julia> function myfun2(::Type{S}; x1::Vector{T}=zeros(T, 2)) where {S<:Union{T, Complex{T}}} where {T<:Real} 
    z = zeros(S, 2)
    z += x1  
    return z
end
ERROR: UndefVarError: T not defined

The work-around consists of adding an optional dummy parameter of type T:

julia> function myfun3(::Type{S}, dummy::T=zero(T); x1::Vector{T}=zeros(T, 2)) where {S<:Union{T, Complex{T}}} where {T<:Real} 
    z = zeros(S, 2)
    z += x1  
    return z
end

The same happens in julia 0.7.0-beta2.0.

DNF · August 4, 2018, 8:15am

I, too, have observed this behaviour, and would like to know what is going on. I also had to tweak my own solution a bit to make it work. Several times I got the T not defined error, for reasons I could not understand.

What’s going on here?

Azamat · August 4, 2018, 10:32pm

You cannot dispatch on keyword arguments in Julia.