Merge pull request #197 from agerlach/divonne_options

Add missing Divonne options v2

ext/IntegralsCubaExt.jl

@@ -1,8 +1,9 @@
 module IntegralsCubaExt
 
 using Integrals, Cuba
-import Integrals: transformation_if_inf, scale_x, scale_x!, CubaVegas, AbstractCubaAlgorithm,
-        CubaSUAVE, CubaDivonne, CubaCuhre
+import Integrals: transformation_if_inf,
+    scale_x, scale_x!, CubaVegas, AbstractCubaAlgorithm,
+    CubaSUAVE, CubaDivonne, CubaCuhre
 
 function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgorithm,
         sensealg,
@@ -116,4 +117,4 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori
         chi = out.probability, retcode = ReturnCode.Success)
 end
 
-end
+end

ext/IntegralsCubatureExt.jl

@@ -2,23 +2,23 @@ module IntegralsCubatureExt
 
 using Integrals, Cubature
 
-import Integrals: transformation_if_inf, scale_x, scale_x!, CubatureJLh, CubatureJLp,
-        AbstractCubatureJLAlgorithm
+import Integrals: transformation_if_inf,
+    scale_x, scale_x!, CubatureJLh, CubatureJLp,
+    AbstractCubatureJLAlgorithm
 import Cubature: INDIVIDUAL, PAIRED, L1, L2, LINF
 
 Integrals.CubatureJLh(; error_norm = Cubature.INDIVIDUAL) = CubatureJLh(error_norm)
 Integrals.CubatureJLp(; error_norm = Cubature.INDIVIDUAL) = CubatureJLp(error_norm)
 
 function Integrals.__solvebp_call(prob::IntegralProblem,
-    alg::AbstractCubatureJLAlgorithm,
-    sensealg, domain, p;
-    reltol = 1e-8, abstol = 1e-8,
-    maxiters = typemax(Int))
-
+        alg::AbstractCubatureJLAlgorithm,
+        sensealg, domain, p;
+        reltol = 1e-8, abstol = 1e-8,
+        maxiters = typemax(Int))
     lb, ub = domain
     mid = (lb + ub) / 2
 
-       # we get to pick fdim or not based on the IntegralFunction and its output dimensions
+    # we get to pick fdim or not based on the IntegralFunction and its output dimensions
     y = if prob.f isa BatchIntegralFunction
         isinplace(prob.f) ? prob.f.integrand_prototype :
         mid isa Number ? prob.f(eltype(mid)[], p) :
@@ -176,12 +176,12 @@ function Integrals.__solvebp_call(prob::IntegralProblem,
         if prob.batch == 0
             if isinplace(prob)
                 dx = zeros(eltype(lb), prob.nout)
-@@ -181,6 +334,7 @@ function Integrals.__solvebp_call(prob::IntegralProblem,
+    @@ -181,6 +334,7 @@ function Integrals.__solvebp_call(prob::IntegralProblem,
             end
         end
     end
     =#
     SciMLBase.build_solution(prob, alg, val, err, retcode = ReturnCode.Success)
 end
 
-end
+end

ext/IntegralsForwardDiffExt.jl

@@ -85,7 +85,7 @@ function Integrals.__solvebp(cache, alg, sensealg, domain,
     res = reinterpret(reshape, DT, dual.u)
     # unwrap the dual when the primal would return a scalar
     out = if (cache.f isa BatchIntegralFunction && y isa AbstractVector) ||
-        !(y isa AbstractArray)
+             !(y isa AbstractArray)
         only(res)
     else
         res

src/Integrals.jl

@@ -169,7 +169,8 @@ function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, domain, p;
     out = vegas(f, lb, ub, rtol = reltol, atol = abstol,
         maxiter = maxiters, nbins = alg.nbins, debug = alg.debug,
         ncalls = ncalls, batch = prob.f isa BatchIntegralFunction)
-    val, err, chi = out isa Tuple ? out : (out.integral_estimate, out.standard_deviation, out.chi_squared_average)
+    val, err, chi = out isa Tuple ? out :
+                    (out.integral_estimate, out.standard_deviation, out.chi_squared_average)
     SciMLBase.build_solution(prob, alg, val, err, chi = chi, retcode = ReturnCode.Success)
 end
 

src/algorithms.jl

@@ -246,8 +246,8 @@ year={1981},
 publisher={ACM New York, NY, USA}
 }
 """
-struct CubaDivonne{R1, R2, R3} <:
-       AbstractCubaAlgorithm where {R1 <: Real, R2 <: Real, R3 <: Real}
+struct CubaDivonne{R1, R2, R3, R4} <:
+       AbstractCubaAlgorithm where {R1 <: Real, R2 <: Real, R3 <: Real, R4 <: Real}
     flags::Int
     seed::Int
     minevals::Int
@@ -258,6 +258,9 @@ struct CubaDivonne{R1, R2, R3} <:
     border::R1
     maxchisq::R2
     mindeviation::R3
+    xgiven::Matrix{R4}
+    nextra::Int
+    peakfinder::Ptr{Cvoid}
 end
 """
     CubaCuhre()
@@ -293,9 +296,11 @@ function CubaSUAVE(; flags = 0, seed = 0, minevals = 0, nnew = 1000, nmin = 2,
 end
 function CubaDivonne(; flags = 0, seed = 0, minevals = 0,
         key1 = 47, key2 = 1, key3 = 1, maxpass = 5, border = 0.0,
-        maxchisq = 10.0, mindeviation = 0.25)
+        maxchisq = 10.0, mindeviation = 0.25,
+        xgiven = zeros(Cdouble, 0, 0),
+        nextra = 0, peakfinder = C_NULL)
     CubaDivonne(flags, seed, minevals, key1, key2, key3, maxpass, border, maxchisq,
-        mindeviation)
+        mindeviation, xgiven, nextra, peakfinder)
 end
 CubaCuhre(; flags = 0, minevals = 0, key = 0) = CubaCuhre(flags, minevals, key)
 
@@ -325,7 +330,6 @@ struct CubatureJLh <: AbstractCubatureJLAlgorithm
     error_norm::Int32
 end
 
-
 """
     CubatureJLp()
 

test/interface_tests.jl

@@ -7,26 +7,26 @@ max_nout_test = 2
 reltol = 1e-3
 abstol = 1e-3
 
-algs = [QuadGKJL(), HCubatureJL(), CubatureJLh(), CubatureJLp(), VEGAS(), #CubaVegas(),
-    CubaSUAVE(), CubaDivonne(), CubaCuhre()]
+algs = [QuadGKJL, HCubatureJL, CubatureJLh, CubatureJLp, VEGAS, #CubaVegas,
+    CubaSUAVE, CubaDivonne, CubaCuhre]
 
-alg_req = Dict(QuadGKJL() => (nout = 1, allows_batch = false, min_dim = 1, max_dim = 1,
+alg_req = Dict(QuadGKJL => (nout = 1, allows_batch = false, min_dim = 1, max_dim = 1,
         allows_iip = false),
-    HCubatureJL() => (nout = Inf, allows_batch = false, min_dim = 1,
+    HCubatureJL => (nout = Inf, allows_batch = false, min_dim = 1,
         max_dim = Inf, allows_iip = true),
-    VEGAS() => (nout = 1, allows_batch = true, min_dim = 2, max_dim = Inf,
+    VEGAS => (nout = 1, allows_batch = true, min_dim = 2, max_dim = Inf,
         allows_iip = true),
-    CubatureJLh() => (nout = Inf, allows_batch = true, min_dim = 1,
+    CubatureJLh => (nout = Inf, allows_batch = true, min_dim = 1,
         max_dim = Inf, allows_iip = true),
-    CubatureJLp() => (nout = Inf, allows_batch = true, min_dim = 1,
+    CubatureJLp => (nout = Inf, allows_batch = true, min_dim = 1,
         max_dim = Inf, allows_iip = true),
-    CubaVegas() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf,
+    CubaVegas => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf,
         allows_iip = true),
-    CubaSUAVE() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf,
+    CubaSUAVE => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf,
         allows_iip = true),
-    CubaDivonne() => (nout = Inf, allows_batch = true, min_dim = 2,
+    CubaDivonne => (nout = Inf, allows_batch = true, min_dim = 2,
         max_dim = Inf, allows_iip = true),
-    CubaCuhre() => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf,
+    CubaCuhre => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf,
         allows_iip = true))
 
 integrands = [
@@ -96,7 +96,7 @@ end
         for i in 1:length(integrands)
             prob = IntegralProblem(integrands[i], lb, ub)
             @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-            sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+            sol = solve(prob, alg(), reltol = reltol, abstol = abstol)
             @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2
         end
     end
@@ -110,11 +110,11 @@ end
             for dim in 1:max_dim_test
                 lb, ub = (ones(dim), 3ones(dim))
                 prob = IntegralProblem(integrands[i], lb, ub)
-                if dim > req.max_dim || dim < req.min_dim || alg isa QuadGKJL  #QuadGKJL requires numbers, not single element arrays
+                if dim > req.max_dim || dim < req.min_dim || alg() isa QuadGKJL  #QuadGKJL requires numbers, not single element arrays
                     continue
                 end
                 @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                sol = solve(prob, alg(), reltol = reltol, abstol = abstol)
                 @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2
             end
         end
@@ -129,14 +129,14 @@ end
             for dim in 1:max_dim_test
                 lb, ub = (ones(dim), 3ones(dim))
                 prob = IntegralProblem(iip_integrands[i], lb, ub)
-                if dim > req.max_dim || dim < req.min_dim || alg isa QuadGKJL  #QuadGKJL requires numbers, not single element arrays
+                if dim > req.max_dim || dim < req.min_dim || alg() isa QuadGKJL  #QuadGKJL requires numbers, not single element arrays
                     continue
                 end
                 @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                if alg isa HCubatureJL && dim == 1 # HCubature library requires finer tol to pass test. When requiring array outputs for iip integrands
-                    sol = solve(prob, alg, reltol = 1e-5, abstol = 1e-5)
+                if alg() isa HCubatureJL && dim == 1 # HCubature library requires finer tol to pass test. When requiring array outputs for iip integrands
+                    sol = solve(prob, alg(), reltol = 1e-5, abstol = 1e-5)
                 else
-                    sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                    sol = solve(prob, alg(), reltol = reltol, abstol = abstol)
                 end
                 if sol.u isa Number
                     @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2
@@ -159,7 +159,7 @@ end
                 continue
             end
             @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-            sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+            sol = solve(prob, alg(), reltol = reltol, abstol = abstol)
             @test sol.u[1]≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2
         end
     end
@@ -177,7 +177,7 @@ end
                     continue
                 end
                 @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                sol = solve(prob, alg(), reltol = reltol, abstol = abstol)
                 if sol.u isa Number
                     @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2
                 else
@@ -201,7 +201,7 @@ end
                     continue
                 end
                 @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                sol = solve(prob, alg(), reltol = reltol, abstol = abstol)
                 if sol.u isa Number
                     @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2
                 else
@@ -226,7 +226,7 @@ end
                     continue
                 end
                 @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                sol = solve(prob, alg(), reltol = reltol, abstol = abstol)
                 if nout == 1
                     @test sol.u[1]≈exact_sol_v[i](dim, nout, lb, ub)[1] rtol=1e-2
                 else
@@ -245,14 +245,14 @@ end
                 lb, ub = (ones(dim), 3ones(dim))
                 for nout in 1:max_nout_test
                     if dim > req.max_dim || dim < req.min_dim || req.nout < nout ||
-                       alg isa QuadGKJL || alg isa VEGAS
+                       alg() isa QuadGKJL || alg() isa VEGAS
                         #QuadGKJL and VEGAS require numbers, not single element arrays
                         continue
                     end
                     prob = IntegralProblem((x, p) -> integrands_v[i](x, p, nout), lb, ub,
                         nout = nout)
                     @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                    sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                    sol = solve(prob, alg(), reltol = reltol, abstol = abstol)
                     if nout == 1
                         @test sol.u[1]≈exact_sol_v[i](dim, nout, lb, ub)[1] rtol=1e-2
                     else
@@ -274,14 +274,14 @@ end
                     prob = IntegralProblem((dx, x, p) -> iip_integrands_v[i](dx, x, p, nout),
                         lb, ub, nout = nout)
                     if dim > req.max_dim || dim < req.min_dim || req.nout < nout ||
-                       alg isa QuadGKJL  #QuadGKJL requires numbers, not single element arrays
+                       alg() isa QuadGKJL  #QuadGKJL requires numbers, not single element arrays
                         continue
                     end
                     @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout"
                     if alg isa HCubatureJL && dim == 1 # HCubature library requires finer tol to pass test. When requiring array outputs for iip integrands
-                        sol = solve(prob, alg, reltol = 1e-5, abstol = 1e-5)
+                        sol = solve(prob, alg(), reltol = 1e-5, abstol = 1e-5)
                     else
-                        sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                        sol = solve(prob, alg(), reltol = reltol, abstol = abstol)
                     end
                     if nout == 1
                         @test sol.u[1]≈exact_sol_v[i](dim, nout, lb, ub)[1] rtol=1e-2
@@ -306,7 +306,7 @@ end
                 continue
             end
             @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-            sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+            sol = solve(prob, alg(), reltol = reltol, abstol = abstol)
             @test sol.u≈exact_sol_v[i](dim, nout, lb, ub) rtol=1e-2
         end
     end
@@ -327,7 +327,7 @@ end
                     continue
                 end
                 @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                sol = solve(prob, alg(), reltol = reltol, abstol = abstol)
                 @test sol.u≈exact_sol_v[i](dim, nout, lb, ub) rtol=1e-2
             end
         end
@@ -348,7 +348,7 @@ end
                     continue
                 end
                 @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                sol = solve(prob, alg(), reltol = reltol, abstol = abstol)
                 @test sol.u≈exact_sol_v[i](dim, nout, lb, ub) rtol=1e-2
             end
         end
@@ -375,7 +375,7 @@ end
         end
         for i in 1:length(integrands)
             prob = IntegralProblem(integrands[i], lb, ub, p)
-            cache = init(prob, alg, reltol = reltol, abstol = abstol)
+            cache = init(prob, alg(), reltol = reltol, abstol = abstol)
             @test solve!(cache).u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2
             cache.lb = lb = 0.5
             @test solve!(cache).u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2

test/nested_ad_tests.jl

@@ -5,42 +5,43 @@ my_function(x, p) = x^2 + p[1]^3 * x + p[2]^2
 function my_integration(p)
     my_problem = IntegralProblem(my_function, -1.0, 1.0, p)
     # return solve(my_problem, HCubatureJL(), reltol=1e-3, abstol=1e-3)  # Works
-    return solve(my_problem, CubatureJLh(), reltol=1e-3, abstol=1e-3)  # Errors
+    return solve(my_problem, CubatureJLh(), reltol = 1e-3, abstol = 1e-3)  # Errors
 end
 my_solution = my_integration(my_parameters)
 @test ForwardDiff.jacobian(my_integration, my_parameters) == [0.0 8.0]
-@test ForwardDiff.jacobian(x->ForwardDiff.jacobian(my_integration, x), my_parameters) == [0.0 0.0
-                                                                                          0.0 4.0]
+@test ForwardDiff.jacobian(x -> ForwardDiff.jacobian(my_integration, x), my_parameters) ==
+      [0.0 0.0
+    0.0 4.0]
 
-ff(x,p) = sum(sin.(x .* p))
+ff(x, p) = sum(sin.(x .* p))
 lb = ones(2)
 ub = 3ones(2)
-p = [1.5,2.0]
+p = [1.5, 2.0]
 
 function testf(p)
-    prob = IntegralProblem(ff,lb,ub,p)
-    sin(solve(prob,CubaCuhre(),reltol=1e-6,abstol=1e-6)[1])
+    prob = IntegralProblem(ff, lb, ub, p)
+    sin(solve(prob, CubaCuhre(), reltol = 1e-6, abstol = 1e-6)[1])
 end
 
-hp1 = FiniteDiff.finite_difference_hessian(testf,p)
-hp2 = ForwardDiff.hessian(testf,p)
-@test hp1 ≈ hp2 atol=1e-4
+hp1 = FiniteDiff.finite_difference_hessian(testf, p)
+hp2 = ForwardDiff.hessian(testf, p)
+@test hp1≈hp2 atol=1e-4
 
-ff2(x,p) = x*p[1].+p[2]*p[3]
-lb =1.0
+ff2(x, p) = x * p[1] .+ p[2] * p[3]
+lb = 1.0
 ub = 3.0
 p = [2.0, 3.0, 4.0]
 _ff3 = BatchIntegralFunction(ff2)
-prob = IntegralProblem(_ff3,lb,ub,p)
+prob = IntegralProblem(_ff3, lb, ub, p)
 
-function testf3(lb,ub,p; f=_ff3)
-    prob = IntegralProblem(_ff3,lb,ub,p)
-    solve(prob, CubatureJLh(); reltol=1e-3,abstol=1e-3)[1]
+function testf3(lb, ub, p; f = _ff3)
+    prob = IntegralProblem(_ff3, lb, ub, p)
+    solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3)[1]
 end
 
-dp1 = ForwardDiff.gradient(p->testf3(lb,ub,p),p)
-dp2 = Zygote.gradient(p->testf3(lb,ub,p),p)[1]
-dp3 = FiniteDiff.finite_difference_gradient(p->testf3(lb,ub,p),p)
+dp1 = ForwardDiff.gradient(p -> testf3(lb, ub, p), p)
+dp2 = Zygote.gradient(p -> testf3(lb, ub, p), p)[1]
+dp3 = FiniteDiff.finite_difference_gradient(p -> testf3(lb, ub, p), p)
 
 @test dp1 ≈ dp3
 @test dp2 ≈ dp3

test/runtests.jl

@@ -8,4 +8,4 @@ using Test
 @time @safetestset "Gaussian Quadrature Tests" include("gaussian_quadrature_tests.jl")
 @time @safetestset "Sampled Integration Tests" include("sampled_tests.jl")
 @time @safetestset "QuadratureFunction Tests" include("quadrule_tests.jl")
-@time @safetestset "Nested AD Tests" include("nested_ad_tests.jl")
+@time @safetestset "Nested AD Tests" include("nested_ad_tests.jl")