From 746ce581e76f7e2e235f4b82a4f0fe704a1a2ac9 Mon Sep 17 00:00:00 2001
From: Pablo Zubieta <pabloferz@yahoo.com.mx>
Date: Tue, 1 Aug 2017 17:35:46 -0500
Subject: [PATCH 1/3] Improve choosing dividing-axis

---
 src/genz-malik.jl | 36 +++++++++++++++++++++++-------------
 1 file changed, 23 insertions(+), 13 deletions(-)

diff --git a/src/genz-malik.jl b/src/genz-malik.jl
index f2ca0c5..6e75efc 100644
--- a/src/genz-malik.jl
+++ b/src/genz-malik.jl
@@ -11,7 +11,7 @@ with k components equal to λ and other components equal to zero.
 """
 function combos(k::Integer, λ::T, ::Type{Val{n}}) where {n, T<:Number}
     combos = Combinatorics.combinations(1:n, k)
-    p = Array{SVector{n,T}}(length(combos))
+    p = Vector{SVector{n,T}}(length(combos))
     v = MVector{n,T}()
     for (i,c) in enumerate(combos)
         v[:] = 0
@@ -31,7 +31,7 @@ with k components equal to ±λ and other components equal to zero
 function signcombos(k::Integer, λ::T, ::Type{Val{n}}) where {n, T<:Number}
     combos = Combinatorics.combinations(1:n, k)
     twoᵏ = 1 << k
-    p = Array{SVector{n,T}}(length(combos) * twoᵏ)
+    p = Vector{SVector{n,T}}(length(combos) * twoᵏ)
     v = MVector{n,T}()
     for (i,c) in enumerate(combos)
         j = (i-1)*twoᵏ + 1
@@ -114,17 +114,17 @@ error `E` (via the given `norm`), and the suggested coordinate `k` ∈ `1:n`
 to subdivide next.
 """
 function (g::GenzMalik{n,T})(f, a::SVector{n}, b::SVector{n}, norm=vecnorm) where {n,T}
-    c = (a+b)*T(0.5)
-    Δ = (b-a)*T(0.5)
-    V = abs(prod(Δ))
+    c = T(0.5).*(a.+b)
+    Δ = T(0.5).*abs.(b.-a)
+    V = prod(Δ)
 
     f₁ = f(c)
 
     f₂ = zero(f₁)
     f₃ = zero(f₁)
-    kdivide = 1
-    maxdivdiff = zero(norm(f₁))
     twelvef₁ = 12f₁
+    maxdivdiff = zero(norm(f₁))
+    divdiff = MVector{n,typeof(maxdivdiff)}()
     for i = 1:n
         p₂ = Δ .* g.p[1][i]
         f₂ᵢ = f(c + p₂) + f(c - p₂)
@@ -134,11 +134,7 @@ function (g::GenzMalik{n,T})(f, a::SVector{n}, b::SVector{n}, norm=vecnorm) wher
         f₃ += f₃ᵢ
         # fourth divided difference: f₃ᵢ-2f₁ - 7*(f₂ᵢ-2f₁),
         # where 7 = (λ₃/λ₂)^2 [see van Dooren and de Ridder]
-        divdiff = norm(f₃ᵢ + twelvef₁ - 7*f₂ᵢ)
-        if divdiff > maxdivdiff
-            kdivide = i
-            maxdivdiff = divdiff
-        end
+        divdiff[i] = norm(f₃ᵢ + twelvef₁ - 7*f₂ᵢ)
     end
 
     f₄ = zero(f₁)
@@ -153,5 +149,19 @@ function (g::GenzMalik{n,T})(f, a::SVector{n}, b::SVector{n}, norm=vecnorm) wher
 
     I = V * (g.w[1]*f₁ + g.w[2]*f₂ + g.w[3]*f₃ + g.w[4]*f₄ + g.w[5]*f₅)
     I′ = V * (g.w′[1]*f₁ + g.w′[2]*f₂ + g.w′[3]*f₃ + g.w′[4]*f₄)
-    return I, norm(I - I′), kdivide
+    E = norm(I - I′)
+
+    # choose axis
+    kdivide = 1
+    δf = 0.001 * E / V
+    for i = 1:n
+        if (δ = divdiff[i] - maxdivdiff) > δf
+            kdivide = i
+            maxdivdiff = divdiff[i]
+        elseif abs(δ) < δf && Δ[i] > Δ[kdivide]
+            kdivide = i
+        end
+    end
+
+    return I, E, kdivide
 end

From 0cca5bbafb43071f5d3b3e8434b22f814ee2f884 Mon Sep 17 00:00:00 2001
From: Pablo Zubieta <pabloferz@yahoo.com.mx>
Date: Tue, 1 Aug 2017 18:41:39 -0500
Subject: [PATCH 2/3] Add a test for issue #4

---
 test/runtests.jl | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/test/runtests.jl b/test/runtests.jl
index 2e3a9d0..a4e67d4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -44,3 +44,10 @@ end
       end
       @test HCubature.countevals(HCubature.Trivial()) == 1
 end
+
+@testset "axischoosing" begin
+    #Issue #4
+    let f = x -> 1.0 + (x[1] * x[3] * sin(x[2]))^2
+        @test hcubature(f, (0.0,0.0,-0.2), (0.2,2π,0.2), rtol=1e-6)[1] ≈ (22502//140625)*π rtol=1e-6
+    end
+end

From fa0c58e8a102a0d32dc7ad43d2fd3d31551211ee Mon Sep 17 00:00:00 2001
From: Pablo Zubieta <pabloferz@yahoo.com.mx>
Date: Thu, 16 Nov 2017 16:52:51 -0600
Subject: [PATCH 3/3] Improve choosing dividing-axis

---
 src/genz-malik.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/genz-malik.jl b/src/genz-malik.jl
index 6e75efc..96e0f37 100644
--- a/src/genz-malik.jl
+++ b/src/genz-malik.jl
@@ -153,7 +153,7 @@ function (g::GenzMalik{n,T})(f, a::SVector{n}, b::SVector{n}, norm=vecnorm) wher
 
     # choose axis
     kdivide = 1
-    δf = 0.001 * E / V
+    δf = E / (10^n * V)
     for i = 1:n
         if (δ = divdiff[i] - maxdivdiff) > δf
             kdivide = i