# Hcubature - New to Julia

**URL:** <https://discourse.julialang.org/t/hcubature-new-to-julia/77155>\
**Category:** New to Julia\
**Tags:** integral\
**Created:** [February 27, 2022, 3:22pm UTC](https://discourse.julialang.org/t/hcubature-new-to-julia/77155 "2022-02-27T15:22:14Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![maria96](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maria96/32/205557_2.png) [@maria96](https://discourse.julialang.org/u/maria96)\
**Post date:** [February 27, 2022, 3:22pm UTC](https://discourse.julialang.org/t/hcubature-new-to-julia/77155/1 "2022-02-27T15:22:14Z")

</div>

Hi all! I’m very new to Julia so this may be quite simple I don’t know.  
I have two functions written in Julia that I want to integrate and both should give a value very close to 1. They both produce the same type of variables but one integral works and the other doesn’t. Can someone help me?

**The one that works**

```julia
using Cubature
function joedensity(u,v,alpha)
    w = 1-u
    z = 1-v
    (w^alpha+z^alpha-(w*z)^alpha)^(1/alpha-2)*(w*z)^(alpha-1)*(alpha-1+w^alpha+z^alpha-(w*z)^alpha)
end

u = [0.1, 0.5, 0.3, 0.7, 0.9]
v = [0.01, 0.2, 0.8, 0.3, 0.99]
alpha = 3
joedensity.(u,v,alpha)
#= Results
5-element Vector{Float64}:
 2.407462338226518
 1.0414427450537544
 0.2669954330753848
 0.5695056921157302
 0.1997673764704445
=#
hcubature(x -> joedensity.(x[1],x[2],alpha),[0,0],[1,1],abstol=0)
#= Results
(1.0000000000808127, 9.994730188502498e-9)
=#

```

**The one that doesn’t**

_this function calls other functions that I defined, but it is working exactly as I want_

```julia
using Cubature
function cop2019(u,v,delta,thetaW)
    y = zeros(2,length(u))
    y[1,:] .= quanty.(u,delta)
    y[2,:] .= quanty.(v,delta)
    map(x -> jointdist(x,delta,thetaW), eachcol(y)) ./ (marg.(y[1,:],delta) .* marg.(y[2,:],delta))
end 

u = [0.1, 0.5, 0.3, 0.7, 0.9]
v = [0.01, 0.2, 0.8, 0.3, 0.99]
delta = 0.2
thetaW = 0.4;
cop2019(u,v,delta,thetaW)
#= Results
5-element Vector{Float64}:
 2.044469250947409
 0.9381167412987401
 0.7335531873911258
 0.840984518753334
 2.3684851388756365
=#
hcubature(x -> cop2019(x[1],x[2],delta,thetaW),[0,0],[1,1],abstol=0)
#= Part of the error
MethodError: Cannot `convert` an object of type Vector{Float64} to an object of type Float64
Closest candidates are:
  convert(::Type{T}, ::Base.TwicePrecision) where T<:Number at ~/julia-1.7.2/share/julia/base/twiceprecision.jl:262
  convert(::Type{T}, ::AbstractChar) where T<:Number at ~/julia-1.7.2/share/julia/base/char.jl:185
  convert(::Type{T}, ::CartesianIndex{1}) where T<:Number at ~/julia-1.7.2/share/julia/base/multidimensional.jl:136
  ...
=#

```

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [February 27, 2022, 4:00pm UTC](https://discourse.julialang.org/t/hcubature-new-to-julia/77155/2 "2022-02-27T16:00:20Z")

</div>

It expects your function to return a scalar, but it is returning a vector.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [February 27, 2022, 4:50pm UTC](https://discourse.julialang.org/t/hcubature-new-to-julia/77155/3 "2022-02-27T16:50:42Z")

</div>

Yeah, you apparently have the impression that you need to work on vector input and return vectors as well. But that’s not how this works. Your function should take scalar input and return scalar output.

You often see languages like matlab an python requiring vectorized functions for quadrature, for performance reasons. But that’s not optimal, adaptive quadrature/cubature functions can do better if they work on individual samples, and can also avoid array allocations.

---

<div class="post-metadata">

**Author:** ![Elrod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elrod/32/22461_2.png) [@Elrod](https://discourse.julialang.org/u/Elrod)\
**Post date:** [February 27, 2022, 6:01pm UTC](https://discourse.julialang.org/t/hcubature-new-to-julia/77155/4 "2022-02-27T18:01:31Z")

</div>

It may be worth noting that [Cubature.jl](https://github.com/JuliaMath/Cubature.jl#parallelizing-the-integrand-evaluation), as opposed to HCubature.jl, does support using functions that evaluate batches of points in vectors.

---

<div class="post-metadata">

**Author:** ![maria96](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maria96/32/205557_2.png) [@maria96](https://discourse.julialang.org/u/maria96)\
**Post date:** [February 27, 2022, 10:35pm UTC](https://discourse.julialang.org/t/hcubature-new-to-julia/77155/5 "2022-02-27T22:35:05Z")

</div>

Hey, thanks! I changed the cop2019 function to now return return a scalar. The double integral now runs although it gives a DomainError that I’m not sure why it’s happening.

I’m rewriting all my R code in Julia in the hopes it’s a bit faster (so far I haven’t seen that happen, which is a shame) and there it works with the vectors u and v I defined. I will keep trying

---

<div class="post-metadata">

**Author:** ![maria96](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maria96/32/205557_2.png) [@maria96](https://discourse.julialang.org/u/maria96)\
**Post date:** [February 27, 2022, 10:36pm UTC](https://discourse.julialang.org/t/hcubature-new-to-julia/77155/6 "2022-02-27T22:36:11Z")

</div>

Hi, thanks! I think I’m using Cubature.jl. Should I try with HCubature.jl instead?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 27, 2022, 10:41pm UTC](https://discourse.julialang.org/t/hcubature-new-to-julia/77155/7 "2022-02-27T22:41:48Z")

</div>

> [@maria96](#):
>
> Hi, thanks! I think I’m using Cubature.jl. Should I try with HCubature.jl instead?

HCubature is pure-Julia code and is a bit more flexible in some ways. e.g. it can work more types of return values, such as `SVector`-valued functions. Cubature, on the other hand, supports batched integrand evaluation, which can be useful for parallelization, and also has `pcubature` routines. (For hcubature, they use essentially the same algorithm.)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 27, 2022, 10:45pm UTC](https://discourse.julialang.org/t/hcubature-new-to-julia/77155/8 "2022-02-27T22:45:48Z")

</div>

> [@maria96](#):
>
> The double integral now runs although it gives a DomainError that I’m not sure why it’s happening.

You have a `x^(1/alpha-2)` exponentiation, i.e. a fractional power (whose value is only known at runtime). That will give a `DomainError` in Julia for negative real `x`, unless you switch to complex numbers.

For a fixed exponent, like `alpha == 3`, it will be much faster to hardcode the exponent. (And you should [avoid globals](https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-global-variables) anyway.) And `x^(1/3-2)` can be computed efficiently by `cbrt(x)/x^2` and it will work even for negative real `x`.