JuliaMath / HCubature.jl

pure-Julia multidimensional h-adaptive integration
Other
148 stars 23 forks source link

Type instability when using a vector for the limits. #35

Closed caetanotc closed 3 years ago

caetanotc commented 3 years ago

Hi, I noticed that I was getting a type instability when using hcubature() inside a function I created. I did some tests and that just happens when I use a vector for a and b, the same didn't happen with a tuple.

Here is what I get from @code_warntype:

julia> using HCubature

julia> f(x,y) = x + y
f (generic function with 1 method)

julia> @code_warntype hcubature(x->f(x[1],x[2]),[0,0],[1,1])
Variables
  #self#::Core.Compiler.Const(HCubature.hcubature, false)
  f::Core.Compiler.Const(var"#9#10"(), false)
  a::Array{Int64,1}
  b::Array{Int64,1}
  norm::typeof(LinearAlgebra.norm)
  rtol::Int64
  atol::Int64
  maxevals::Int64
  initdiv::Int64

Body::Tuple{Any,Any}
1 ─ %1 = HCubature.norm::Core.Compiler.Const(LinearAlgebra.norm, false)
│        (norm = %1)
│        (rtol = 0)
│        (atol = 0)
│        (maxevals = HCubature.typemax(HCubature.Int))
│        (initdiv = 1)
│   %7 = HCubature.:(var"#hcubature#3")(norm, rtol::Core.Compiler.Const(0, false), atol::Core.Compiler.Const(0, false), maxevals::Core.Compiler.Const(9223372036854775807, false), initdiv::Core.Compiler.Const(1, false), #self#, f, a, b)::Tuple{Any,Any}
└──      return %7

julia> @code_warntype hcubature(x->f(x[1],x[2]),(0,0),(1,1))
Variables
  #self#::Core.Compiler.Const(HCubature.hcubature, false)
  f::Core.Compiler.Const(var"#11#12"(), false)
  a::Tuple{Int64,Int64}
  b::Tuple{Int64,Int64}
  norm::typeof(LinearAlgebra.norm)
  rtol::Int64
  atol::Int64
  maxevals::Int64
  initdiv::Int64

Body::Tuple{Float64,Float64}
1 ─ %1 = HCubature.norm::Core.Compiler.Const(LinearAlgebra.norm, false)
│        (norm = %1)
│        (rtol = 0)
│        (atol = 0)
│        (maxevals = HCubature.typemax(HCubature.Int))
│        (initdiv = 1)
│   %7 = HCubature.:(var"#hcubature#3")(norm, rtol::Core.Compiler.Const(0, false), atol::Core.Compiler.Const(0, false), maxevals::Core.Compiler.Const(9223372036854775807, false), initdiv::Core.Compiler.Const(1, false), #self#, f, a, b)::Tuple{Float64,Float64}
└──      return %7

In the documentation, it says that both options can be used. Any ideas on why this is happening?

stevengj commented 3 years ago

In the documentation, it says that both options can be used. Any ideas on why this is happening?

Both options can be used, but a tuple allows Julia's compiler to infer the length, whereas a vector does not. Internally, HCubature uses StaticArrays.jl for speed, but that means that a separate version of hcubature is compiled for different dimensionalities of the integral, and using a non-inferable length means that dynamic dispatch is required.

There's nothing wrong with dynamic dispatch — since multidimensional integration is typically expensive, the overhead of dynamically determining the dimensionality and dispatching to the appropriate compiled version is probably low.