Merge pull request #5 from aurelio-amerio/new_method

Added HCubature as a new integration method

Project.toml

@@ -1,16 +1,18 @@
 name = "MultiQuad"
 uuid = "8fca5bbe-a785-4bb7-9793-c0c8ddca5d4b"
 authors = ["Aurelio Amerio <aure.amerio@gmail.com>"]
-version = "1.0.1"
+version = "1.1.0"
 
 [deps]
 Cuba = "8a292aeb-7a57-582c-b821-06e4c11590b1"
+HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 Cuba = "2.0.0"
+HCubature = "1.4.0"
 QuadGK = "2.3.1"
 Unitful = "0.18, 1.0"
 julia = "1.0"

src/MultiQuad.jl

@@ -1,13 +1,19 @@
 module MultiQuad
 
-using QuadGK, Cuba, Unitful
+using QuadGK, Cuba, HCubature, Unitful
 
 export quad, dblquad, tplquad
 
 @doc raw"""
     quad(arg::Function, x1, x2; method = :quadgk, kwargs...)
 
 Performs the integral ``\int_{x1}^{x2}f(x)dx``
+
+Available integration methods:
+- `:suave`
+- `:vegas`
+- `:quadgk`
+
 See [QuadGK](https://github.com/JuliaMath/QuadGK.jl) and [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments.
 
 # Examples
@@ -63,7 +69,16 @@ end
     dblquad(arg::Function, x1, x2, y1::Function, y2::Function; method = :cuhre, kwargs...)
 
 Performs the integral ``\int_{x1}^{x2}\int_{y1(x)}^{y2(x)}f(y,x)dydx``
-See [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments.
+
+Available integration methods:
+- `:cuhre`
+- `:divonne`
+- `:suave`
+- `:vegas`
+- `:hcubature`
+
+See [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments fro the `:cuhre`, `:divonne`, `:suave` and `:vegas` methods.
+See [HCubature](https://github.com/stevengj/HCubature.jl) for all the available keywords for the `:hcubature` method.
 
 # Examples
 ```jldoctest
@@ -98,6 +113,8 @@ function dblquad(
         integrate = suave
     elseif method == :vegas
         integrate = vegas
+    elseif method == :hcubature
+        integrate = hcubature
     else
         ex = ErrorException("Integration method $method is not supported!")
         throw(ex)
@@ -110,11 +127,21 @@ function dblquad(
 
     arg2(a, b) = ustrip(units, (x2 - x1) * arg1(a, (x2 - x1) * b + x1))::Float64
 
-    function integrand(x, f)
-        f[1] = arg2(x[1], x[2])
-    end
+    if method == :hcubature
+        function integrand(arr)
+            return arg2(arr[1], arr[2])
+        end
+
+        min_arr = [0, 0]
+        max_arr = [1, 1]
+        result, err = integrate(integrand, min_arr, max_arr; kwargs...)
+    else
+        function integrand2(x, f)
+            f[1] = arg2(x[1], x[2])
+        end
 
-    result, err = integrate(integrand, 2, 1; kwargs...)
+        result, err = integrate(integrand2, 2, 1; kwargs...)
+    end
 
     if units == Unitful.NoUnits
         return result[1], err[1]
@@ -127,7 +154,16 @@ end
     tplquad(arg::Function, x1, x2, y1::Function, y2::Function, z1::Function, z2::Function; method = :cuhre, kwargs...)
 
 Performs the integral ``\int_{x1}^{x2}\int_{y1(x)}^{y2(x)}\int_{z1(x,y)}^{z2(x,y)}f(z,y,x)dzdydx``
-See [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments.
+
+Available integration methods:
+- `:cuhre`
+- `:divonne`
+- `:suave`
+- `:vegas`
+- `:hcubature`
+
+See [Cuba.jl](https://giordano.github.io/Cuba.jl/stable/) for all the available keyword arguments fro the `:cuhre`, `:divonne`, `:suave` and `:vegas` methods.
+See [HCubature](https://github.com/stevengj/HCubature.jl) for all the available keywords for the `:hcubature` method.
 
 # Examples
 ```jldoctest
@@ -164,6 +200,8 @@ function tplquad(
         integrate = suave
     elseif method == :vegas
         integrate = vegas
+    elseif method == :hcubature
+        integrate = hcubature
     else
         ex = ErrorException("Integration method $method is not supported!")
         throw(ex)
@@ -181,11 +219,21 @@ function tplquad(
     arg2(a, b, c) =
         ustrip(units, (x2 - x1) * arg1(a, b, (x2 - x1) * c + x1))::Float64
 
-    function integrand(x, f)
-        f[1] = arg2(x[1], x[2], x[3])
-    end
+    if method == :hcubature
+        function integrand(arr)
+            return arg2(arr[1], arr[2], arr[3])
+        end
+
+        min_arr = [0, 0, 0]
+        max_arr = [1, 1, 1]
+        result, err = integrate(integrand, min_arr, max_arr; kwargs...)
+    else
+        function integrand2(x, f)
+            f[1] = arg2(x[1], x[2], x[3])
+        end
 
-    result, err = integrate(integrand, 3, 1; kwargs...)
+        result, err = integrate(integrand2, 3, 1; kwargs...)
+    end
 
     if units == Unitful.NoUnits
         return result[1], err[1]

test/runtests.jl

@@ -58,8 +58,26 @@ end
 
     func1(y, x) = y * x^3
     res1 = 241664 / 5
-    int1, err1 = dblquad(func1, xmin, xmax, ymin, ymax, rtol = rtol)
+    int1, err1 = dblquad(
+        func1,
+        xmin,
+        xmax,
+        ymin,
+        ymax,
+        rtol = rtol,
+        method = :cuhre,
+    )
+    int1b, err1b = dblquad(
+        func1,
+        xmin,
+        xmax,
+        ymin,
+        ymax,
+        rtol = rtol,
+        method = :hcubature,
+    )
     @test isapprox(int1, res1, atol = err1)
+    @test isapprox(int1b, res1, atol = err1b)
     @test typeof(dblquad(
         func1,
         xmin,
@@ -95,8 +113,18 @@ end
     func2(y, x) = y * x^3
 
     int2, err2 = dblquad(func2, xmin2, xmax2, ymin2, ymax2, rtol = rtol)
+    int2b, err2b = dblquad(
+        func2,
+        xmin2,
+        xmax2,
+        ymin2,
+        ymax2,
+        rtol = rtol,
+        method = :hcubature,
+    )
 
     @test isapprox(int2, res1 * u"m^10", atol = err2)
+    @test isapprox(int2b, res1 * u"m^10", atol = err2b)
     @test typeof(dblquad(
         func2,
         xmin2,
@@ -149,7 +177,19 @@ end
     func1(z, y, x) = y * x^3 * sin(z)
     res = 241664 / 5 * (cos(3) - cos(4))
     int1, err1 = tplquad(func1, xmin, xmax, ymin, ymax, zmin, zmax, rtol = rtol)
+    int1b, err1b = tplquad(
+        func1,
+        xmin,
+        xmax,
+        ymin,
+        ymax,
+        zmin,
+        zmax,
+        rtol = rtol,
+        method = :hcubature,
+    )
     @test isapprox(int1, res, atol = err1)
+    @test isapprox(int1b, res, atol = err1b)
     @test typeof(tplquad(
         func1,
         xmin,
@@ -202,7 +242,19 @@ end
         zmax2,
         rtol = rtol,
     )
+    int2b, err2b = tplquad(
+        func2,
+        xmin2,
+        xmax2,
+        ymin2,
+        ymax2,
+        zmin2,
+        zmax2,
+        rtol = rtol,
+        method = :hcubature,
+    )
     @test isapprox(int2, res * u"m^10", atol = err2)
+    @test isapprox(int2b, res * u"m^10", atol = err2b)
     @test typeof(tplquad(
         func2,
         xmin2,