From 75bf7a1f7d9f899462d448a1473f2597ace1ea15 Mon Sep 17 00:00:00 2001 From: lxvm Date: Fri, 3 Nov 2023 18:34:55 -0400 Subject: [PATCH 1/9] initial commit --- test/derivative_tests.jl | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 87b326f7..6470c7d5 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -2,6 +2,29 @@ using Integrals, Zygote, FiniteDiff, ForwardDiff#, SciMLSensitivity using Cuba, Cubature using Test +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, + allows_iip = false), + 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, + allows_iip = true), + CubatureJLh() => (nout = Inf, allows_batch = true, min_dim = 1, + max_dim = Inf, allows_iip = true), + 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, + allows_iip = true), + CubaSUAVE() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, + allows_iip = true), + 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, + allows_iip = true)) + + ### One Dimensional f(x, p) = sum(sin.(p[1] * x)) From 6c0aa44febcec16455497d878fd6667ba9c704ed Mon Sep 17 00:00:00 2001 From: lxvm Date: Sat, 30 Dec 2023 22:51:49 -0800 Subject: [PATCH 2/9] overhaul tests --- test/derivative_tests.jl | 400 ++++++++++++++++++--------------------- 1 file changed, 186 insertions(+), 214 deletions(-) diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 6470c7d5..4538ddb6 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -2,265 +2,237 @@ using Integrals, Zygote, FiniteDiff, ForwardDiff#, SciMLSensitivity using Cuba, Cubature using Test -algs = [QuadGKJL(), HCubatureJL(), CubatureJLh(), CubatureJLp(), #VEGAS(), #CubaVegas(), - CubaSUAVE(), CubaDivonne(), CubaCuhre()] +max_dim_test = 2 +max_nout_test = 2 -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, - max_dim = Inf, allows_iip = true), - VEGAS() => (nout = 1, allows_batch = true, min_dim = 2, max_dim = Inf, - allows_iip = true), - CubatureJLh() => (nout = Inf, allows_batch = true, min_dim = 1, - max_dim = Inf, allows_iip = true), - 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, - allows_iip = true), - CubaSUAVE() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, +reltol = 1e-3 +abstol = 1e-3 + +alg_req = Dict( + QuadGKJL() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = 1, allows_iip = true), - CubaDivonne() => (nout = Inf, allows_batch = true, min_dim = 2, + HCubatureJL() => (nout = Inf, allows_batch = false, min_dim = 1, max_dim = Inf, allows_iip = true), - CubaCuhre() => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf, - allows_iip = true)) - + # VEGAS() => (nout = 1, allows_batch = true, min_dim = 2, max_dim = Inf, + # allows_iip = true), + # CubatureJLh() => (nout = Inf, allows_batch = true, min_dim = 1, + # max_dim = Inf, allows_iip = true), + # 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, + # allows_iip = true), + # CubaSUAVE() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, + # allows_iip = true), + # 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, + # allows_iip = true), +) +# helper function / test runner +scalarize_solution = ( + sol -> sin(sum(sol)), + sol -> sin(sol[1]), +) + +do_tests = function (; f, scalarize, lb, ub, p, alg, abstol, reltol) + testf = function (lb, ub, p) + prob = IntegralProblem(f, (lb, ub), p) + scalarize(solve(prob, alg; reltol, abstol)) + end + testf(lb, ub, p) + + dlb1, dub1, dp1 = Zygote.gradient(testf, lb, ub, p) + + f_lb = lb -> testf(lb, ub, p) + f_ub = ub -> testf(lb, ub, p) + + dlb = lb isa AbstractArray ? :gradient : :derivative + dub = ub isa AbstractArray ? :gradient : :derivative + + dlb2 = getproperty(FiniteDiff, Symbol(:finite_difference_, dlb))(f_lb, lb) + dub2 = getproperty(FiniteDiff, Symbol(:finite_difference_, dub))(f_ub, ub) + + if lb isa Number + @test dlb1 ≈ dlb2 atol=abstol rtol=reltol + @test dub1 ≈ dub2 atol=abstol rtol=reltol + else # TODO: implement multivariate limit derivatives in ZygoteExt + @test_broken dlb1 ≈ dlb2 atol=abstol rtol=reltol + @test_broken dub1 ≈ dub2 atol=abstol rtol=reltol + end + + # TODO: implement limit derivatives in ForwardDiffExt + @test_broken dlb2 ≈ getproperty(ForwardDiff, dlb)(dfdlb, lb) atol=abstol rtol=reltol + @test_broken dub2 ≈ getproperty(ForwardDiff, dub)(dfdub, ub) atol=abstol rtol=reltol + + f_p = p -> testf(lb, ub, p) + + dp = p isa AbstractArray ? :gradient : :derivative + + dp2 = getproperty(FiniteDiff, Symbol(:finite_difference_, dp))(f_p, p) + dp3 = getproperty(ForwardDiff, dp)(f_p, p) + + @test dp1 ≈ dp2 atol=abstol rtol=reltol + @test dp1 ≈ dp3 atol=abstol rtol=reltol + + return +end ### One Dimensional +f_1d_scalar = (x, p) -> sum(q -> sin(q*x), p) +for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution) + req.nout > 1 || continue + req.min_dim > 0 || continue -f(x, p) = sum(sin.(p[1] * x)) -lb = 1.0 -ub = 3.0 -p = 2.0 -prob = IntegralProblem(f, lb, ub, p) -sol = solve(prob, QuadGKJL(), reltol = 1e-3, abstol = 1e-3) - -function testf(lb, ub, p) - prob = IntegralProblem(f, lb, ub, p) - sin(solve(prob, QuadGKJL(), reltol = 1e-3, abstol = 1e-3)[1]) + @info "One-dimensional, scalar, oop derivative test" alg scalarize=i + do_tests(; f=f_1d_scalar, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) end -dlb1, dub1, dp1 = Zygote.gradient(testf, lb, ub, p) -dlb2 = FiniteDiff.finite_difference_derivative(lb -> testf(lb, ub, p), lb) -dub2 = FiniteDiff.finite_difference_derivative(ub -> testf(lb, ub, p), ub) -dp2 = FiniteDiff.finite_difference_derivative(p -> testf(lb, ub, p), p) -# dlb3 = ForwardDiff.derivative(lb->testf(lb,ub,p),lb) -# dub3 = ForwardDiff.derivative(ub->testf(lb,ub,p),ub) -dp3 = ForwardDiff.derivative(p -> testf(lb, ub, p), p) +## One-dimensional nout +f_1d_nout = (x, p) -> map(q -> q*x, p) +for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test + req.nout > 1 || continue + req.min_dim > 0 || continue -#@test dlb1 ≈ dlb3 -#@test dub1 ≈ dub3 -@test dp1 ≈ dp3 + @info "One-dimensional, multivariate, oop derivative test" alg scalarize=i nout + do_tests(; f=f_1d_nout, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) +end ### N-dimensional +f_nd_scalar = (x, p) -> prod(y -> f_1d_scalar(y, p), x) +for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test + req.nout > 1 || continue + req.min_dim <= dim <= req.max_dim || continue -f(x, p) = sum(sin.(x .* p)) -lb = ones(2) -ub = 3ones(2) -p = [1.5, 2.0] -prob = IntegralProblem(f, lb, ub, p) -sol = solve(prob, CubaCuhre(), reltol = 1e-3, abstol = 1e-3) - -function testf(lb, ub, p) - prob = IntegralProblem(f, lb, ub, p) - sin(solve(prob, CubaCuhre(), reltol = 1e-6, abstol = 1e-6)[1]) -end - -function testf2(lb, ub, p) - prob = IntegralProblem(f, lb, ub, p) - sin(solve(prob, HCubatureJL(), reltol = 1e-6, abstol = 1e-6)[1]) + @info "Multi-dimensional, scalar, oop derivative test" alg scalarize=i, dim + do_tests(; f=f_nd_scalar, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) end -dlb1, dub1, dp1 = Zygote.gradient(testf, lb, ub, p) -dlb2 = FiniteDiff.finite_difference_gradient(lb -> testf(lb, ub, p), lb) -dub2 = FiniteDiff.finite_difference_gradient(ub -> testf(lb, ub, p), ub) -dp2 = FiniteDiff.finite_difference_gradient(p -> testf(lb, ub, p), p) - -@test_broken dlb1 ≈ dlb2 -@test_broken dub1 ≈ dub2 -@test dp1 ≈ dp2 - -# dlb3 = ForwardDiff.gradient(lb->testf(lb,ub,p),lb) -# dub3 = ForwardDiff.gradient(ub->testf(lb,ub,p),ub) -dp3 = ForwardDiff.gradient(p -> testf2(lb, ub, p), p) - -#@test dlb1 ≈ dlb3 -#@test dub1 ≈ dub3 -@test dp1 ≈ dp3 - -# dlb4 = ForwardDiff.gradient(lb->testf(lb,ub,p),lb) -# dub4 = ForwardDiff.gradient(ub->testf(lb,ub,p),ub) -dp4 = ForwardDiff.gradient(p -> testf(lb, ub, p), p) - -#@test dlb1 ≈ dlb4 -#@test dub1 ≈ dub4 -@test dp1 ≈ dp4 - -### N-dimensional N-out - -f(x, p) = sin.(x .* p) -lb = ones(2) -ub = 3ones(2) -p = [1.5, 2.0] -prob = IntegralProblem(f, lb, ub, p, nout = 2) -sol = solve(prob, CubaCuhre(), reltol = 1e-3, abstol = 1e-3) +### N-dimensional nout +f_nd_nout = (x, p) -> mapreduce(y -> f_1d_nout(y, p), +, x) +for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test + req.nout > 1 || continue + req.min_dim <= dim <= req.max_dim || continue -function testf(lb, ub, p) - prob = IntegralProblem(f, lb, ub, p, nout = 2) - sin(sum(solve(prob, CubaCuhre(), reltol = 1e-6, abstol = 1e-6))) + @info "Multi-dimensional, multivariate, oop derivative test" alg scalarize=i dim nout + do_tests(; f=f_nd_nout, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) end -function testf2(lb, ub, p) - prob = IntegralProblem(f, lb, ub, p, nout = 2) - sin(sum(solve(prob, HCubatureJL(), reltol = 1e-6, abstol = 1e-6))) -end - -dlb1, dub1, dp1 = Zygote.gradient(testf, lb, ub, p) -dlb2 = FiniteDiff.finite_difference_gradient(lb -> testf(lb, ub, p), lb) -dub2 = FiniteDiff.finite_difference_gradient(ub -> testf(lb, ub, p), ub) -dp2 = FiniteDiff.finite_difference_gradient(p -> testf(lb, ub, p), p) - -@test_broken dlb1 ≈ dlb2 -@test_broken dub1 ≈ dub2 -@test dp1 ≈ dp2 - -# dlb3 = ForwardDiff.gradient(lb->testf(lb,ub,p),lb) -# dub3 = ForwardDiff.gradient(ub->testf(lb,ub,p),ub) -dp3 = ForwardDiff.gradient(p -> testf2(lb, ub, p), p) - -#@test dlb1 ≈ dlb3 -#@test dub1 ≈ dub3 -@test dp1 ≈ dp3 - -# dlb4 = ForwardDiff.gradient(lb->testf(lb,ub,p),lb) -# dub4 = ForwardDiff.gradient(ub->testf(lb,ub,p),ub) -dp4 = ForwardDiff.gradient(p -> testf(lb, ub, p), p) - -#@test dlb1 ≈ dlb4 -#@test dub1 ≈ dub4 -@test dp1 ≈ dp4 - -### Batch Single dim -f(x, p) = x * p[1] .+ p[2] * p[3] # scalar integrand - -lb = 1.0 -ub = 3.0 -p = [2.0, 3.0, 4.0] -prob = IntegralProblem(f, lb, ub, p) +### One Dimensional +f_1d_scalar_iip = (y, x, p) -> y .= f_1d_scalar(x, p) +for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution) + req.nout > 1 || continue + req.min_dim > 0 || continue -function testf3(lb, ub, p; f = f) - prob = IntegralProblem(f, lb, ub, p, batch = 10, nout = 1) - solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3)[1] + @info "One-dimensional, scalar, iip derivative test" alg scalarize=i + do_tests(; f=IntegralFunction(f_1d_scalar_iip, zeros(1)), scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) end -dp1 = ForwardDiff.gradient(p -> testf3(lb, ub, p), p) -dp2 = Zygote.gradient(p -> testf3(lb, ub, p), p)[1] # TODO fix: LoadError: DimensionMismatch("variable with size(x) == (1, 15) cannot have a gradient with size(dx) == (15,)") -dp3 = FiniteDiff.finite_difference_gradient(p -> testf3(lb, ub, p), p) - -@test dp1 ≈ dp3 #passes -@test dp2 ≈ dp3 #passes +## One-dimensional nout +f_1d_nout_iip = (y, x, p) -> y .= f_1d_nout(x, p) +for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test + req.nout > 1 || continue + req.min_dim > 0 || continue -### Batch single dim, nout -f(x, p) = (x' * p[1] .+ p[2] * p[3]) .* [1; 2] - -lb = 1.0 -ub = 3.0 -p = [2.0, 3.0, 4.0] -prob = IntegralProblem(f, lb, ub, p) - -function testf3(lb, ub, p; f = f) - prob = IntegralProblem(f, lb, ub, p, batch = 10, nout = 2) - sum(solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3)) + @info "One-dimensional, multivariate, iip derivative test" alg scalarize=i nout + do_tests(; f=IntegralFunction(f_1d_nout, zeros(nout)), scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) 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) - -@test dp1 ≈ dp3 #passes -@test dp2 ≈ dp3 #passes +### N-dimensional +f_nd_scalar_iip = (y, x, p) -> y .= f_nd_scalar(x, p) +for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test + req.nout > 1 || continue + req.min_dim <= dim <= req.max_dim || continue -### Batch multi dim -f(x, p) = x[1, :] * p[1] .+ p[2] * p[3] + @info "Multi-dimensional, scalar, iip derivative test" alg scalarize=i, dim + do_tests(; f=IntegralFunction(f_nd_scalar_iip, zeros(1)), scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) +end -lb = [1.0, 1.0] -ub = [3.0, 3.0] -p = [2.0, 3.0, 4.0] -prob = IntegralProblem(f, lb, ub, p) +### N-dimensional nout iip +f_nd_nout_iip = (y, x, p) -> y .= f_nd_nout(x, p) +for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test + req.nout > 1 || continue + req.min_dim <= dim <= req.max_dim || continue -function testf3(lb, ub, p; f = f) - prob = IntegralProblem(f, lb, ub, p, batch = 10, nout = 1) - solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3)[1] + @info "Multi-dimensional, multivariate, iip derivative test" alg scalarize=i dim nout + do_tests(; f=IntegralFunction(f_nd_nout, zeros(nout)), scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) 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) +### Batch, One Dimensional +bf_1d_scalar = (x, p) -> map(y -> f_1d_scalar(y, p), x) +for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution) + req.nout > 1 || continue + req.min_dim > 0 || continue -@test dp1 ≈ dp3 -@test dp2 ≈ dp3 - -### Batch multi dim nout -f(x, p) = x * p[1] .+ p[2] * p[3] + @info "Batched, one-dimensional, scalar, oop derivative test" alg scalarize=i + do_tests(; f=BatchIntegralFunction(bf_1d_scalar), scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) +end -lb = [1.0, 1.0] -ub = [3.0, 3.0] -p = [2.0, 3.0, 4.0] -prob = IntegralProblem(f, lb, ub, p) +## Batch, One-dimensional nout +bf_1d_nout = (x, p) -> stack(y -> f_1d_nout(y, p), x) +for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test + req.nout > 1 || continue + req.min_dim > 0 || continue -function testf3(lb, ub, p; f = f) - prob = IntegralProblem(f, lb, ub, p, batch = 10, nout = 2) - sum(solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3)) + @info "Batched, one-dimensional, multivariate, oop derivative test" alg scalarize=i nout + do_tests(; f=BatchIntegralFunction(bf_1d_nout), scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) 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) - -@test dp1 ≈ dp3 -@test dp2 ≈ dp3 +### Batch, N-dimensional +bf_nd_scalar = (x, p) -> map(y -> f_nd_scalar(y, p), eachslice(x; dims=ndims(x))) +for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test + req.nout > 1 || continue + req.min_dim <= dim <= req.max_dim || continue -## iip Batch multi dim -function g(dx, x, p) - dx .= sum(x * p[1] .+ p[2] * p[3], dims = 1) + @info "Batched, multi-dimensional, scalar, oop derivative test" alg scalarize=i, dim + do_tests(; f=BatchIntegralFunction(bf_nd_scalar), scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) end -lb = [1.0, 1.0] -ub = [3.0, 3.0] -p = [2.0, 3.0, 4.0] -prob = IntegralProblem(g, lb, ub, p) +### Batch, N-dimensional nout +bf_nd_nout = (x, p) -> stack(y -> f_nd_nout(y, p), eachslice(x; dims=ndims(x))) +for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test + req.nout > 1 || continue + req.min_dim <= dim <= req.max_dim || continue -function testf3(lb, ub, p; f = g) - prob = IntegralProblem(f, lb, ub, p, batch = 10, nout = 1) - solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3)[1] + @info "Batch, multi-dimensional, multivariate, oop derivative test" alg scalarize=i dim nout + do_tests(; f=BatchIntegralFunction(bf_nd_nout), scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) end -testf3(lb, ub, p) +### Batch, one-dimensional +bf_1d_scalar_iip = (y, x, p) -> y .= bf_1d_scalar(x, p) +for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution) + req.nout > 1 || continue + req.min_dim > 0 || continue -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) + @info "Batched, one-dimensional, scalar, iip derivative test" alg scalarize=i + do_tests(; f=BatchIntegralFunction(bf_1d_scalar_iip, zeros(1, 0)), scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) +end -@test dp1 ≈ dp3 -@test dp2 ≈ dp3 +## Batch, one-dimensional nout +bf_1d_nout_iip = (y, x, p) -> y .= bf_1d_nout(x, p) +for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test + req.nout > 1 || continue + req.min_dim > 0 || continue -## iip Batch mulit dim nout -function g(dx, x, p) - dx .= x * p[1] .+ p[2] * p[3] + @info "Batched, one-dimensional, multivariate, iip derivative test" alg scalarize=i nout + do_tests(; f=BatchIntegralFunction(bf_1d_nout, zeros(nout, 0)), scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) end -lb = [1.0, 1.0] -ub = [3.0, 3.0] -p = [2.0, 3.0, 4.0] -prob = IntegralProblem(g, lb, ub, p) +### Batch, N-dimensional +bf_nd_scalar_iip = (y, x, p) -> y .= bf_nd_scalar(x, p) +for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test + req.nout > 1 || continue + req.min_dim <= dim <= req.max_dim || continue -function testf3(lb, ub, p; f = g) - prob = IntegralProblem(f, lb, ub, p, batch = 10, nout = 2) - sum(solve(prob, CubatureJLh(); reltol = 1e-3, abstol = 1e-3)) + @info "Batched, multi-dimensional, scalar, iip derivative test" alg scalarize=i, dim + do_tests(; f=BatchIntegralFunction(bf_nd_scalar_iip, zeros(1, 0)), scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) 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) +### Batch, N-dimensional nout iip +bf_nd_nout_iip = (y, x, p) -> y .= bf_nd_nout(x, p) +for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test + req.nout > 1 || continue + req.min_dim <= dim <= req.max_dim || continue -@test dp1 ≈ dp3 -@test dp2 ≈ dp3 + @info "Batched, multi-dimensional, multivariate, iip derivative test" alg scalarize=i dim nout + do_tests(; f=BatchIntegralFunction(bf_nd_nout, zeros(nout, 0)), scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) +end From 93120a30fc4395edaef4cf5aa8b414f39f844273 Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 1 Jan 2024 10:16:13 -0800 Subject: [PATCH 3/9] remove unused test dep --- Project.toml | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 868f2451..ebe9c3b2 100644 --- a/Project.toml +++ b/Project.toml @@ -48,7 +48,6 @@ QuadGK = "2.9" Reexport = "1.0" SafeTestsets = "0.1" SciMLBase = "2.6" -SciMLSensitivity = "7.41" StaticArrays = "1" Test = "1" Zygote = "0.6.60" @@ -66,10 +65,9 @@ FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" -SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" [targets] -test = ["Aqua", "Arblib", "SciMLSensitivity", "StaticArrays", "FiniteDiff", "Pkg", "SafeTestsets", "Test", "Distributions", "ForwardDiff", "Zygote", "ChainRulesCore", "FastGaussQuadrature", "Cuba", "Cubature"] +test = ["Aqua", "Arblib", "StaticArrays", "FiniteDiff", "Pkg", "SafeTestsets", "Test", "Distributions", "ForwardDiff", "Zygote", "ChainRulesCore", "FastGaussQuadrature", "Cuba", "Cubature"] From 33ae84cabc4aae93909959219945412f34eba037 Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 1 Jan 2024 10:20:42 -0800 Subject: [PATCH 4/9] fix forwarddiff bugs --- ext/IntegralsForwardDiffExt.jl | 26 ++++--- test/derivative_tests.jl | 120 ++++++++++++++++++++------------- 2 files changed, 89 insertions(+), 57 deletions(-) diff --git a/ext/IntegralsForwardDiffExt.jl b/ext/IntegralsForwardDiffExt.jl index 5aec07c3..300b6063 100644 --- a/ext/IntegralsForwardDiffExt.jl +++ b/ext/IntegralsForwardDiffExt.jl @@ -3,7 +3,8 @@ using Integrals isdefined(Base, :get_extension) ? (using ForwardDiff) : (using ..ForwardDiff) ### Forward-Mode AD Intercepts -# Direct AD on solvers with QuadGK and HCubature +#= Direct AD on solvers with QuadGK and HCubature +# incompatible with iip since types must change function Integrals.__solvebp(cache, alg::QuadGKJL, sensealg, domain, p::AbstractArray{<:ForwardDiff.Dual{T, V, P}, N}; kwargs...) where {T, V, P, N} @@ -15,11 +16,14 @@ function Integrals.__solvebp(cache, alg::HCubatureJL, sensealg, domain, kwargs...) where {T, V, P, N} Integrals.__solvebp_call(cache, alg, sensealg, domain, p; kwargs...) end +=# + +# TODO: add the pushforward for derivative w.r.t lb, and ub (and then combinations?) # Manually split for the pushforward function Integrals.__solvebp(cache, alg, sensealg, domain, - p::AbstractArray{<:ForwardDiff.Dual{T, V, P}, N}; - kwargs...) where {T, V, P, N} + p::Union{D,AbstractArray{<:D}}; + kwargs...) where {T, V, P, D<:ForwardDiff.Dual{T, V, P}} # we need the output type to avoid perturbation confusion while unwrapping nested duals # We compute a vector-valued integral of the primal and dual simultaneously @@ -32,11 +36,11 @@ function Integrals.__solvebp(cache, alg, sensealg, domain, len, size(cache.f.integrand_prototype)...) - dfdp_ = function (out, x, p) - dualp = reinterpret(ForwardDiff.Dual{T, V, P}, p) + dfdp_ = function (out, x, _p) + dualp = reinterpret(ForwardDiff.Dual{T, V, P}, _p) dout = reinterpret(reshape, DT, out) - cache.f(dout, x, dualp) - return out + cache.f(dout, x, p isa D ? only(dualp) : reshape(dualp, size(p))) + return end dfdp = if cache.f isa BatchIntegralFunction BatchIntegralFunction{true}(dfdp_, dual_prototype) @@ -55,9 +59,9 @@ function Integrals.__solvebp(cache, alg, sensealg, domain, DT = y isa AbstractArray ? eltype(y) : typeof(y) elt = unwrap_dualvaltype(DT) - dfdp_ = function (x, p) - dualp = reinterpret(ForwardDiff.Dual{T, V, P}, p) - ys = cache.f(x, dualp) + dfdp_ = function (x, _p) + dualp = reinterpret(ForwardDiff.Dual{T, V, P}, _p) + ys = cache.f(x, p isa D ? only(dualp) : reshape(dualp, size(p))) ys_ = ys isa AbstractArray ? ys : [ys] # we need to reshape in order for batching to be consistent return reinterpret(reshape, elt, ys_) @@ -70,7 +74,7 @@ function Integrals.__solvebp(cache, alg, sensealg, domain, end ForwardDiff.can_dual(elt) || ForwardDiff.throw_cannot_dual(elt) - rawp = copy(reinterpret(V, p)) + rawp = p isa D ? reinterpret(V, [p]) : copy(reinterpret(V, vec(p))) prob = Integrals.build_problem(cache) dp_prob = remake(prob, f = dfdp, p = rawp) diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 4538ddb6..068b7488 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -15,18 +15,18 @@ alg_req = Dict( max_dim = Inf, allows_iip = true), # VEGAS() => (nout = 1, allows_batch = true, min_dim = 2, max_dim = Inf, # allows_iip = true), - # CubatureJLh() => (nout = Inf, allows_batch = true, min_dim = 1, - # max_dim = Inf, allows_iip = true), - # CubatureJLp() => (nout = Inf, allows_batch = true, min_dim = 1, - # max_dim = Inf, allows_iip = true), + CubatureJLh() => (nout = Inf, allows_batch = true, min_dim = 1, + max_dim = Inf, allows_iip = true), + 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, # allows_iip = true), - # CubaSUAVE() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, - # allows_iip = true), - # 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, - # allows_iip = true), + CubaSUAVE() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, + allows_iip = true), + 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, + allows_iip = true), ) # helper function / test runner scalarize_solution = ( @@ -41,7 +41,7 @@ do_tests = function (; f, scalarize, lb, ub, p, alg, abstol, reltol) end testf(lb, ub, p) - dlb1, dub1, dp1 = Zygote.gradient(testf, lb, ub, p) + # dlb1, dub1, dp1 = Zygote.gradient(testf, lb, ub, p) f_lb = lb -> testf(lb, ub, p) f_ub = ub -> testf(lb, ub, p) @@ -52,13 +52,13 @@ do_tests = function (; f, scalarize, lb, ub, p, alg, abstol, reltol) dlb2 = getproperty(FiniteDiff, Symbol(:finite_difference_, dlb))(f_lb, lb) dub2 = getproperty(FiniteDiff, Symbol(:finite_difference_, dub))(f_ub, ub) - if lb isa Number - @test dlb1 ≈ dlb2 atol=abstol rtol=reltol - @test dub1 ≈ dub2 atol=abstol rtol=reltol - else # TODO: implement multivariate limit derivatives in ZygoteExt - @test_broken dlb1 ≈ dlb2 atol=abstol rtol=reltol - @test_broken dub1 ≈ dub2 atol=abstol rtol=reltol - end + # if lb isa Number + # @test dlb1 ≈ dlb2 atol=abstol rtol=reltol + # @test dub1 ≈ dub2 atol=abstol rtol=reltol + # else # TODO: implement multivariate limit derivatives in ZygoteExt + # @test_broken dlb1 ≈ dlb2 atol=abstol rtol=reltol + # @test_broken dub1 ≈ dub2 atol=abstol rtol=reltol + # end # TODO: implement limit derivatives in ForwardDiffExt @test_broken dlb2 ≈ getproperty(ForwardDiff, dlb)(dfdlb, lb) atol=abstol rtol=reltol @@ -71,14 +71,49 @@ do_tests = function (; f, scalarize, lb, ub, p, alg, abstol, reltol) dp2 = getproperty(FiniteDiff, Symbol(:finite_difference_, dp))(f_p, p) dp3 = getproperty(ForwardDiff, dp)(f_p, p) - @test dp1 ≈ dp2 atol=abstol rtol=reltol - @test dp1 ≈ dp3 atol=abstol rtol=reltol + # @test dp1 ≈ dp2 atol=abstol rtol=reltol + @test dp2 ≈ dp3 atol=abstol rtol=reltol return end -### One Dimensional f_1d_scalar = (x, p) -> sum(q -> sin(q*x), p) +f_1d_nout = (x, p) -> map(q -> q*x, p) +f_nd_scalar = (x, p) -> prod(y -> f_1d_scalar(y, p), x) +f_nd_nout = (x, p) -> mapreduce(y -> f_1d_nout(y, p), +, x) + +f_1d_scalar_iip = (y, x, p) -> y .= f_1d_scalar(x, p) +f_1d_nout_iip = (y, x, p) -> y .= f_1d_nout(x, p) +f_nd_scalar_iip = (y, x, p) -> y .= f_nd_scalar(x, p) +f_nd_nout_iip = (y, x, p) -> y .= f_nd_nout(x, p) + +bf_helper = (f, x, p) -> begin + elt = typeof(zero(eltype(x))*zero(eltype(p))) # output type of above functions + if p isa AbstractArray + # p and f_*_nout are of size nout + # this is like a call to stack that should also work for empty arrays + out = similar(p, elt, size(p)..., size(x, ndims(x))) + for (v,y) in zip(eachslice(out; dims=ndims(out)), eachslice(x; dims=ndims(x))) + v .= f(x isa AbstractVector ? only(y) : y, p) + end + out + else + elt[f(x isa AbstractVector ? only(y) : y, p) for y in eachslice(x; dims=ndims(x))] + end +end + +bf_1d_scalar = (x, p) -> bf_helper(f_1d_scalar, x, p) +bf_1d_nout = (x, p) -> bf_helper(f_1d_nout, x, p) +bf_nd_scalar = (x, p) -> bf_helper(f_nd_scalar, x, p) +bf_nd_nout = (x, p) -> bf_helper(f_nd_nout, x, p) + +bf_1d_nout_iip = (y, x, p) -> y .= bf_1d_nout(x, p) +bf_1d_scalar_iip = (y, x, p) -> y .= bf_1d_scalar(x, p) +bf_nd_scalar_iip = (y, x, p) -> y .= bf_nd_scalar(x, p) +bf_nd_nout_iip = (y, x, p) -> y .= bf_nd_nout(x, p) + + +### One Dimensional for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution) req.nout > 1 || continue req.min_dim > 0 || continue @@ -88,7 +123,6 @@ for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution end ## One-dimensional nout -f_1d_nout = (x, p) -> map(q -> q*x, p) for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test req.nout > 1 || continue req.min_dim > 0 || continue @@ -98,17 +132,15 @@ for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution end ### N-dimensional -f_nd_scalar = (x, p) -> prod(y -> f_1d_scalar(y, p), x) for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue - @info "Multi-dimensional, scalar, oop derivative test" alg scalarize=i, dim + @info "Multi-dimensional, scalar, oop derivative test" alg scalarize=i dim do_tests(; f=f_nd_scalar, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) end ### N-dimensional nout -f_nd_nout = (x, p) -> mapreduce(y -> f_1d_nout(y, p), +, x) for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue @@ -118,7 +150,6 @@ for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution end ### One Dimensional -f_1d_scalar_iip = (y, x, p) -> y .= f_1d_scalar(x, p) for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution) req.nout > 1 || continue req.min_dim > 0 || continue @@ -128,38 +159,35 @@ for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution end ## One-dimensional nout -f_1d_nout_iip = (y, x, p) -> y .= f_1d_nout(x, p) for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test req.nout > 1 || continue req.min_dim > 0 || continue @info "One-dimensional, multivariate, iip derivative test" alg scalarize=i nout - do_tests(; f=IntegralFunction(f_1d_nout, zeros(nout)), scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) + do_tests(; f=IntegralFunction(f_1d_nout_iip, zeros(nout)), scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) end ### N-dimensional -f_nd_scalar_iip = (y, x, p) -> y .= f_nd_scalar(x, p) for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue - @info "Multi-dimensional, scalar, iip derivative test" alg scalarize=i, dim + @info "Multi-dimensional, scalar, iip derivative test" alg scalarize=i dim do_tests(; f=IntegralFunction(f_nd_scalar_iip, zeros(1)), scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) end ### N-dimensional nout iip -f_nd_nout_iip = (y, x, p) -> y .= f_nd_nout(x, p) for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue @info "Multi-dimensional, multivariate, iip derivative test" alg scalarize=i dim nout - do_tests(; f=IntegralFunction(f_nd_nout, zeros(nout)), scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) + do_tests(; f=IntegralFunction(f_nd_nout_iip, zeros(nout)), scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) end ### Batch, One Dimensional -bf_1d_scalar = (x, p) -> map(y -> f_1d_scalar(y, p), x) for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution) + req.allows_batch || continue req.nout > 1 || continue req.min_dim > 0 || continue @@ -168,8 +196,8 @@ for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution end ## Batch, One-dimensional nout -bf_1d_nout = (x, p) -> stack(y -> f_1d_nout(y, p), x) for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test + req.allows_batch || continue req.nout > 1 || continue req.min_dim > 0 || continue @@ -178,18 +206,18 @@ for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution end ### Batch, N-dimensional -bf_nd_scalar = (x, p) -> map(y -> f_nd_scalar(y, p), eachslice(x; dims=ndims(x))) for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test + req.allows_batch || continue req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue - @info "Batched, multi-dimensional, scalar, oop derivative test" alg scalarize=i, dim + @info "Batched, multi-dimensional, scalar, oop derivative test" alg scalarize=i dim do_tests(; f=BatchIntegralFunction(bf_nd_scalar), scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) end ### Batch, N-dimensional nout -bf_nd_nout = (x, p) -> stack(y -> f_nd_nout(y, p), eachslice(x; dims=ndims(x))) for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test + req.allows_batch || continue req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue @@ -198,41 +226,41 @@ for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution end ### Batch, one-dimensional -bf_1d_scalar_iip = (y, x, p) -> y .= bf_1d_scalar(x, p) for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution) + req.allows_batch || continue req.nout > 1 || continue req.min_dim > 0 || continue @info "Batched, one-dimensional, scalar, iip derivative test" alg scalarize=i - do_tests(; f=BatchIntegralFunction(bf_1d_scalar_iip, zeros(1, 0)), scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) + do_tests(; f=BatchIntegralFunction(bf_1d_scalar_iip, zeros(0)), scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) end ## Batch, one-dimensional nout -bf_1d_nout_iip = (y, x, p) -> y .= bf_1d_nout(x, p) for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test + req.allows_batch || continue req.nout > 1 || continue req.min_dim > 0 || continue @info "Batched, one-dimensional, multivariate, iip derivative test" alg scalarize=i nout - do_tests(; f=BatchIntegralFunction(bf_1d_nout, zeros(nout, 0)), scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) + do_tests(; f=BatchIntegralFunction(bf_1d_nout_iip, zeros(nout, 0)), scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) end ### Batch, N-dimensional -bf_nd_scalar_iip = (y, x, p) -> y .= bf_nd_scalar(x, p) for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test + req.allows_batch || continue req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue - @info "Batched, multi-dimensional, scalar, iip derivative test" alg scalarize=i, dim - do_tests(; f=BatchIntegralFunction(bf_nd_scalar_iip, zeros(1, 0)), scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) + @info "Batched, multi-dimensional, scalar, iip derivative test" alg scalarize=i dim + do_tests(; f=BatchIntegralFunction(bf_nd_scalar_iip, zeros(0)), scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) end ### Batch, N-dimensional nout iip -bf_nd_nout_iip = (y, x, p) -> y .= bf_nd_nout(x, p) for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test + req.allows_batch || continue req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue @info "Batched, multi-dimensional, multivariate, iip derivative test" alg scalarize=i dim nout - do_tests(; f=BatchIntegralFunction(bf_nd_nout, zeros(nout, 0)), scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) + do_tests(; f=BatchIntegralFunction(bf_nd_nout_iip, zeros(nout, 0)), scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) end From c6f2b7d07e0212210800e87ee2bf7be2044af300 Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 1 Jan 2024 10:33:59 -0800 Subject: [PATCH 5/9] fix dimension test --- test/derivative_tests.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 068b7488..7c1d7b6b 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -116,7 +116,7 @@ bf_nd_nout_iip = (y, x, p) -> y .= bf_nd_nout(x, p) ### One Dimensional for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution) req.nout > 1 || continue - req.min_dim > 0 || continue + req.min_dim == 1 || continue @info "One-dimensional, scalar, oop derivative test" alg scalarize=i do_tests(; f=f_1d_scalar, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) @@ -125,7 +125,7 @@ end ## One-dimensional nout for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test req.nout > 1 || continue - req.min_dim > 0 || continue + req.min_dim == 1 || continue @info "One-dimensional, multivariate, oop derivative test" alg scalarize=i nout do_tests(; f=f_1d_nout, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) @@ -152,7 +152,7 @@ end ### One Dimensional for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution) req.nout > 1 || continue - req.min_dim > 0 || continue + req.min_dim == 1 || continue @info "One-dimensional, scalar, iip derivative test" alg scalarize=i do_tests(; f=IntegralFunction(f_1d_scalar_iip, zeros(1)), scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) @@ -161,7 +161,7 @@ end ## One-dimensional nout for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test req.nout > 1 || continue - req.min_dim > 0 || continue + req.min_dim == 1 || continue @info "One-dimensional, multivariate, iip derivative test" alg scalarize=i nout do_tests(; f=IntegralFunction(f_1d_nout_iip, zeros(nout)), scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) @@ -189,7 +189,7 @@ end for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution) req.allows_batch || continue req.nout > 1 || continue - req.min_dim > 0 || continue + req.min_dim == 1 || continue @info "Batched, one-dimensional, scalar, oop derivative test" alg scalarize=i do_tests(; f=BatchIntegralFunction(bf_1d_scalar), scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) @@ -199,7 +199,7 @@ end for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test req.allows_batch || continue req.nout > 1 || continue - req.min_dim > 0 || continue + req.min_dim == 1 || continue @info "Batched, one-dimensional, multivariate, oop derivative test" alg scalarize=i nout do_tests(; f=BatchIntegralFunction(bf_1d_nout), scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) @@ -229,7 +229,7 @@ end for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution) req.allows_batch || continue req.nout > 1 || continue - req.min_dim > 0 || continue + req.min_dim == 1 || continue @info "Batched, one-dimensional, scalar, iip derivative test" alg scalarize=i do_tests(; f=BatchIntegralFunction(bf_1d_scalar_iip, zeros(0)), scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) @@ -239,7 +239,7 @@ end for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test req.allows_batch || continue req.nout > 1 || continue - req.min_dim > 0 || continue + req.min_dim == 1 || continue @info "Batched, one-dimensional, multivariate, iip derivative test" alg scalarize=i nout do_tests(; f=BatchIntegralFunction(bf_1d_nout_iip, zeros(nout, 0)), scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) From d026912e3fe1a269141e08c30db8ef4f33e91e2f Mon Sep 17 00:00:00 2001 From: lxvm Date: Mon, 1 Jan 2024 10:43:14 -0800 Subject: [PATCH 6/9] don't test derivatives for MC algorithms --- test/derivative_tests.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 7c1d7b6b..37d19880 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -21,12 +21,12 @@ alg_req = Dict( max_dim = Inf, allows_iip = true), # 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, - allows_iip = true), - 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, - allows_iip = true), + # CubaSUAVE() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = Inf, + # allows_iip = true), + # 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, + # allows_iip = true), ) # helper function / test runner scalarize_solution = ( From d517d67942e0460a9359cccf64ab560b49751f9b Mon Sep 17 00:00:00 2001 From: lxvm Date: Wed, 3 Jan 2024 01:33:50 -0800 Subject: [PATCH 7/9] rrule bugfixes --- ext/IntegralsZygoteExt.jl | 25 ++++++++++++++++--------- 1 file changed, 16 insertions(+), 9 deletions(-) diff --git a/ext/IntegralsZygoteExt.jl b/ext/IntegralsZygoteExt.jl index d5dcc034..90d6f882 100644 --- a/ext/IntegralsZygoteExt.jl +++ b/ext/IntegralsZygoteExt.jl @@ -1,4 +1,5 @@ module IntegralsZygoteExt +using LinearAlgebra: dot using Integrals if isdefined(Base, :get_extension) using Zygote @@ -11,6 +12,7 @@ else end ChainRulesCore.@non_differentiable Integrals.checkkwargs(kwargs...) ChainRulesCore.@non_differentiable Integrals.isinplace(f, args...) # fixes #99 +ChainRulesCore.@non_differentiable Integrals.init_cacheval(alg, prob) function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, sensealg, domain, p; @@ -24,23 +26,25 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal # https://juliadiff.org/ChainRulesCore.jl/dev/design/many_tangents.html#manytypes if isinplace(cache) # zygote doesn't support mutation, so we build an oop pullback - dx = similar(cache.f.integrand_prototype) - _f = x -> cache.f(dx, x, p) if sensealg.vjp isa Integrals.ZygoteVJP if cache.f isa BatchIntegralFunction + dx = similar(cache.f.integrand_prototype, size(cache.f.integrand_prototype)[begin:end-1]..., 1) + _f = x -> (cache.f(dx, x, p); dx) # TODO: let the user pass a batched jacobian so we can return a BatchIntegralFunction dfdp_ = function (x, p) x_ = x isa AbstractArray ? reshape(x, size(x)..., 1) : [x] z, back = Zygote.pullback(p) do p - _dx = Zygote.Buffer(dx, size(dx)[begin:(end - 1)]..., 1) + _dx = Zygote.Buffer(dx) cache.f(_dx, x_, p) copy(_dx) end return back(z .= (Δ isa AbstractArray ? reshape(Δ, size(Δ)..., 1) : - [Δ]))[1] + Δ))[1] end dfdp = IntegralFunction{false}(dfdp_, nothing) else + dx = similar(cache.f.integrand_prototype) + _f = x -> (cache.f(dx, x, p); dx) dfdp_ = function (x, p) _, back = Zygote.pullback(p) do p _dx = Zygote.Buffer(dx) @@ -62,8 +66,7 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal dfdp_ = function (x, p) x_ = x isa AbstractArray ? reshape(x, size(x)..., 1) : [x] z, back = Zygote.pullback(p -> cache.f(x_, p), p) - return back(z .= (Δ isa AbstractArray ? reshape(Δ, size(Δ)..., 1) : - [Δ]))[1] + return back(Δ isa AbstractArray ? reshape(Δ, size(Δ)..., 1) : [Δ])[1] end dfdp = IntegralFunction{false}(dfdp_, nothing) else @@ -98,13 +101,15 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal lb, ub = domain if lb isa Number - dlb = cache.f isa BatchIntegralFunction ? -_f([lb]) : -_f(lb) - dub = cache.f isa BatchIntegralFunction ? _f([ub]) : _f(ub) + # TODO replace evaluation at endpoint (which anyone can do without Integrals.jl) + # with integration of dfdx uing the same quadrature + dlb = cache.f isa BatchIntegralFunction ? -batch_unwrap(_f([lb])) : -_f(lb) + dub = cache.f isa BatchIntegralFunction ? batch_unwrap(_f([ub])) : _f(ub) return (NoTangent(), NoTangent(), NoTangent(), NoTangent(), - Tangent{typeof(domain)}(dlb, dub), + Tangent{typeof(domain)}(dot(dlb, Δ), dot(dub, Δ)), dp) else # we need to compute 2*length(lb) integrals on the faces of the hypercube, as we @@ -123,6 +128,8 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal out, quadrature_adjoint end +batch_unwrap(x::AbstractArray) = dropdims(x; dims=ndims(x)) + Zygote.@adjoint function Zygote.literal_getproperty(sol::SciMLBase.IntegralSolution, ::Val{:u}) sol.u, Δ -> (SciMLBase.build_solution(sol.prob, sol.alg, Δ, sol.resid),) From b024f0d2bed488567fe369d48d4cf2428e2ea1e4 Mon Sep 17 00:00:00 2001 From: lxvm Date: Wed, 3 Jan 2024 01:34:08 -0800 Subject: [PATCH 8/9] rrule tests --- test/derivative_tests.jl | 249 +++++++++++++++++++++++---------------- 1 file changed, 146 insertions(+), 103 deletions(-) diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index 37d19880..f04840e3 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -1,5 +1,6 @@ using Integrals, Zygote, FiniteDiff, ForwardDiff#, SciMLSensitivity using Cuba, Cubature +using FastGaussQuadrature using Test max_dim_test = 2 @@ -9,6 +10,8 @@ reltol = 1e-3 abstol = 1e-3 alg_req = Dict( + QuadratureRule(gausslegendre, n=50) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false, + allows_iip = false), QuadGKJL() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = 1, allows_iip = true), HCubatureJL() => (nout = Inf, allows_batch = false, min_dim = 1, @@ -28,12 +31,67 @@ alg_req = Dict( # CubaCuhre() => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf, # allows_iip = true), ) -# helper function / test runner + +# integrands should have same shape as parameters, independent of dimensionality +integrands =( + (x, p) -> map(q -> prod(y -> sin(y*q), x), p), +) + +# function to turn the output into a scalar / test different tangent types scalarize_solution = ( sol -> sin(sum(sol)), sol -> sin(sol[1]), ) +# we will be able to use broadcasting for this after https://github.com/FluxML/Zygote.jl/pull/1488 +function buffer_copyto!(y, x) + for (j, i) in zip(eachindex(y), eachindex(x)) + y[j] = x[i] + end + return y +end +function f_helper!(f, y, x, p) + buffer_copyto!(y, f(x, p)) + return +end + +# the Zygote implementation is inconsistent about 0-d so we hijack it +struct Scalar{T<:Real} <: Real + x::T +end +Base.iterate(a::Scalar) = (a.x, nothing) +Base.iterate(::Scalar, _) = nothing +Base.IteratorSize(::Type{Scalar{T}}) where {T} = Base.HasShape{0}() +Base.eltype(::Type{Scalar{T}}) where {T} = T +Base.length(a::Scalar) = 1 +Base.size(::Scalar) = () +Base.:+(a::Scalar, b::Scalar) = Scalar(a.x+b.x) +Base.:*(a::Number, b::Scalar) = a*b.x +Base.:*(a::Scalar, b::Number) = a.x*b +Base.:*(a::Scalar, b::Scalar) = Scalar(a.x*b.x) +Base.zero(a::Scalar) = Scalar(zero(a.x)) +Base.map(f, a::Scalar) = map(f, a.x) +(::Type{T})(a::Scalar) where {T<:Real} = T(a.x) +struct ScalarAxes end # the implementation doesn't preserve singleton axes +Base.axes(::Scalar) = ScalarAxes() +Base.iterate(::ScalarAxes) = nothing +Base.reshape(A::AbstractArray, ::ScalarAxes) = Scalar(only(A)) + +# here we assume f evaluated at scalar inputs gives a scalar output +# p will be able to be a number after https://github.com/FluxML/Zygote.jl/pull/1489 +# p will be able to be a 0-array after https://github.com/FluxML/Zygote.jl/pull/1491 +# p can't be either without both prs +function batch_helper(f, x, p) + t = f(zero(eltype(x)), zero(eltype(eltype(p)))) + typeof(t).([f(y, q) for q in p, y in eachslice(x; dims=ndims(x))]) +end + +function batch_helper!(f, y, x, p) + buffer_copyto!(y, batch_helper(f, x, p)) + return +end + +# helper function / test runner do_tests = function (; f, scalarize, lb, ub, p, alg, abstol, reltol) testf = function (lb, ub, p) prob = IntegralProblem(f, (lb, ub), p) @@ -41,7 +99,7 @@ do_tests = function (; f, scalarize, lb, ub, p, alg, abstol, reltol) end testf(lb, ub, p) - # dlb1, dub1, dp1 = Zygote.gradient(testf, lb, ub, p) + dlb1, dub1, dp1 = Zygote.gradient(testf, lb, ub, p isa Number && f isa BatchIntegralFunction ? Scalar(p) : p) f_lb = lb -> testf(lb, ub, p) f_ub = ub -> testf(lb, ub, p) @@ -52,13 +110,13 @@ do_tests = function (; f, scalarize, lb, ub, p, alg, abstol, reltol) dlb2 = getproperty(FiniteDiff, Symbol(:finite_difference_, dlb))(f_lb, lb) dub2 = getproperty(FiniteDiff, Symbol(:finite_difference_, dub))(f_ub, ub) - # if lb isa Number - # @test dlb1 ≈ dlb2 atol=abstol rtol=reltol - # @test dub1 ≈ dub2 atol=abstol rtol=reltol - # else # TODO: implement multivariate limit derivatives in ZygoteExt - # @test_broken dlb1 ≈ dlb2 atol=abstol rtol=reltol - # @test_broken dub1 ≈ dub2 atol=abstol rtol=reltol - # end + if lb isa Number + @test dlb1 ≈ dlb2 atol=abstol rtol=reltol + @test dub1 ≈ dub2 atol=abstol rtol=reltol + else # TODO: implement multivariate limit derivatives in ZygoteExt + @test_broken dlb1 ≈ dlb2 atol=abstol rtol=reltol + @test_broken dub1 ≈ dub2 atol=abstol rtol=reltol + end # TODO: implement limit derivatives in ForwardDiffExt @test_broken dlb2 ≈ getproperty(ForwardDiff, dlb)(dfdlb, lb) atol=abstol rtol=reltol @@ -71,196 +129,181 @@ do_tests = function (; f, scalarize, lb, ub, p, alg, abstol, reltol) dp2 = getproperty(FiniteDiff, Symbol(:finite_difference_, dp))(f_p, p) dp3 = getproperty(ForwardDiff, dp)(f_p, p) - # @test dp1 ≈ dp2 atol=abstol rtol=reltol + @test dp1 ≈ dp2 atol=abstol rtol=reltol @test dp2 ≈ dp3 atol=abstol rtol=reltol return end -f_1d_scalar = (x, p) -> sum(q -> sin(q*x), p) -f_1d_nout = (x, p) -> map(q -> q*x, p) -f_nd_scalar = (x, p) -> prod(y -> f_1d_scalar(y, p), x) -f_nd_nout = (x, p) -> mapreduce(y -> f_1d_nout(y, p), +, x) - -f_1d_scalar_iip = (y, x, p) -> y .= f_1d_scalar(x, p) -f_1d_nout_iip = (y, x, p) -> y .= f_1d_nout(x, p) -f_nd_scalar_iip = (y, x, p) -> y .= f_nd_scalar(x, p) -f_nd_nout_iip = (y, x, p) -> y .= f_nd_nout(x, p) - -bf_helper = (f, x, p) -> begin - elt = typeof(zero(eltype(x))*zero(eltype(p))) # output type of above functions - if p isa AbstractArray - # p and f_*_nout are of size nout - # this is like a call to stack that should also work for empty arrays - out = similar(p, elt, size(p)..., size(x, ndims(x))) - for (v,y) in zip(eachslice(out; dims=ndims(out)), eachslice(x; dims=ndims(x))) - v .= f(x isa AbstractVector ? only(y) : y, p) - end - out - else - elt[f(x isa AbstractVector ? only(y) : y, p) for y in eachslice(x; dims=ndims(x))] - end -end - -bf_1d_scalar = (x, p) -> bf_helper(f_1d_scalar, x, p) -bf_1d_nout = (x, p) -> bf_helper(f_1d_nout, x, p) -bf_nd_scalar = (x, p) -> bf_helper(f_nd_scalar, x, p) -bf_nd_nout = (x, p) -> bf_helper(f_nd_nout, x, p) - -bf_1d_nout_iip = (y, x, p) -> y .= bf_1d_nout(x, p) -bf_1d_scalar_iip = (y, x, p) -> y .= bf_1d_scalar(x, p) -bf_nd_scalar_iip = (y, x, p) -> y .= bf_nd_scalar(x, p) -bf_nd_nout_iip = (y, x, p) -> y .= bf_nd_nout(x, p) - - +#= ### One Dimensional -for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution) +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution) req.nout > 1 || continue - req.min_dim == 1 || continue + req.min_dim <= 1 || continue - @info "One-dimensional, scalar, oop derivative test" alg scalarize=i - do_tests(; f=f_1d_scalar, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) + @info "One-dimensional, scalar, oop derivative test" alg integrand=j scalarize=i + do_tests(; f, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) end ## One-dimensional nout -for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test req.nout > 1 || continue - req.min_dim == 1 || continue + req.min_dim <= 1 || continue - @info "One-dimensional, multivariate, oop derivative test" alg scalarize=i nout - do_tests(; f=f_1d_nout, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) + @info "One-dimensional, multivariate, oop derivative test" alg integrand=j scalarize=i nout + do_tests(; f, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) end ### N-dimensional -for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue - @info "Multi-dimensional, scalar, oop derivative test" alg scalarize=i dim - do_tests(; f=f_nd_scalar, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) + @info "Multi-dimensional, scalar, oop derivative test" alg integrand=j scalarize=i dim + do_tests(; f, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) end ### N-dimensional nout -for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue - @info "Multi-dimensional, multivariate, oop derivative test" alg scalarize=i dim nout - do_tests(; f=f_nd_nout, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) + @info "Multi-dimensional, multivariate, oop derivative test" alg integrand=j scalarize=i dim nout + do_tests(; f, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) end ### One Dimensional -for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution) +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution) + req.allows_iip || continue req.nout > 1 || continue - req.min_dim == 1 || continue + req.min_dim <= 1 || continue - @info "One-dimensional, scalar, iip derivative test" alg scalarize=i - do_tests(; f=IntegralFunction(f_1d_scalar_iip, zeros(1)), scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) + @info "One-dimensional, scalar, iip derivative test" alg integrand=j scalarize=i + fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(1)) + do_tests(; f=fiip, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) end ## One-dimensional nout -for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test + req.allows_iip || continue req.nout > 1 || continue - req.min_dim == 1 || continue + req.min_dim <= 1 || continue - @info "One-dimensional, multivariate, iip derivative test" alg scalarize=i nout - do_tests(; f=IntegralFunction(f_1d_nout_iip, zeros(nout)), scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) + @info "One-dimensional, multivariate, iip derivative test" alg integrand=j scalarize=i nout + fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(nout)) + do_tests(; f=fiip, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) end ### N-dimensional -for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test + req.allows_iip || continue req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue - @info "Multi-dimensional, scalar, iip derivative test" alg scalarize=i dim - do_tests(; f=IntegralFunction(f_nd_scalar_iip, zeros(1)), scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) + @info "Multi-dimensional, scalar, iip derivative test" alg integrand=j scalarize=i dim + fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(1)) + do_tests(; f=fiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) end ### N-dimensional nout iip -for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test + req.allows_iip || continue req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue - @info "Multi-dimensional, multivariate, iip derivative test" alg scalarize=i dim nout - do_tests(; f=IntegralFunction(f_nd_nout_iip, zeros(nout)), scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) + @info "Multi-dimensional, multivariate, iip derivative test" alg integrand=j scalarize=i dim nout + fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(nout)) + do_tests(; f=fiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) end - +=# ### Batch, One Dimensional -for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution) +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution) req.allows_batch || continue req.nout > 1 || continue - req.min_dim == 1 || continue + req.min_dim <= 1 || continue - @info "Batched, one-dimensional, scalar, oop derivative test" alg scalarize=i - do_tests(; f=BatchIntegralFunction(bf_1d_scalar), scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) + @info "Batched, one-dimensional, scalar, oop derivative test" alg integrand=j scalarize=i + bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p)) + do_tests(; f=bf, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) end ## Batch, One-dimensional nout -for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test req.allows_batch || continue req.nout > 1 || continue - req.min_dim == 1 || continue + req.min_dim <= 1 || continue - @info "Batched, one-dimensional, multivariate, oop derivative test" alg scalarize=i nout - do_tests(; f=BatchIntegralFunction(bf_1d_nout), scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) + @info "Batched, one-dimensional, multivariate, oop derivative test" alg integrand=j scalarize=i nout + bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p)) + do_tests(; f=bf, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) end ### Batch, N-dimensional -for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test req.allows_batch || continue req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue - @info "Batched, multi-dimensional, scalar, oop derivative test" alg scalarize=i dim - do_tests(; f=BatchIntegralFunction(bf_nd_scalar), scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) + @info "Batched, multi-dimensional, scalar, oop derivative test" alg integrand=j scalarize=i dim + bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p)) + do_tests(; f=bf, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) end ### Batch, N-dimensional nout -for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test req.allows_batch || continue req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue - @info "Batch, multi-dimensional, multivariate, oop derivative test" alg scalarize=i dim nout - do_tests(; f=BatchIntegralFunction(bf_nd_nout), scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) + @info "Batch, multi-dimensional, multivariate, oop derivative test" alg integrand=j scalarize=i dim nout + bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p)) + do_tests(; f=bf, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) end ### Batch, one-dimensional -for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution) +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution) req.allows_batch || continue + req.allows_iip || continue req.nout > 1 || continue - req.min_dim == 1 || continue + req.min_dim <= 1 || continue - @info "Batched, one-dimensional, scalar, iip derivative test" alg scalarize=i - do_tests(; f=BatchIntegralFunction(bf_1d_scalar_iip, zeros(0)), scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) + @info "Batched, one-dimensional, scalar, iip derivative test" alg integrand=j scalarize=i + bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(0)) + do_tests(; f=bfiip, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol) end ## Batch, one-dimensional nout -for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test req.allows_batch || continue + req.allows_iip || continue req.nout > 1 || continue - req.min_dim == 1 || continue + req.min_dim <= 1 || continue - @info "Batched, one-dimensional, multivariate, iip derivative test" alg scalarize=i nout - do_tests(; f=BatchIntegralFunction(bf_1d_nout_iip, zeros(nout, 0)), scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) + @info "Batched, one-dimensional, multivariate, iip derivative test" alg integrand=j scalarize=i nout + bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(nout, 0)) + do_tests(; f=bfiip, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol) end ### Batch, N-dimensional -for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test req.allows_batch || continue + req.allows_iip || continue req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue - @info "Batched, multi-dimensional, scalar, iip derivative test" alg scalarize=i dim - do_tests(; f=BatchIntegralFunction(bf_nd_scalar_iip, zeros(0)), scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) + @info "Batched, multi-dimensional, scalar, iip derivative test" alg integrand=j scalarize=i dim + bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(0)) + do_tests(; f=bfiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol) end ### Batch, N-dimensional nout iip -for (alg, req) in pairs(alg_req), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test +for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test req.allows_batch || continue + req.allows_iip || continue req.nout > 1 || continue req.min_dim <= dim <= req.max_dim || continue - @info "Batched, multi-dimensional, multivariate, iip derivative test" alg scalarize=i dim nout - do_tests(; f=BatchIntegralFunction(bf_nd_nout_iip, zeros(nout, 0)), scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) + @info "Batched, multi-dimensional, multivariate, iip derivative test" alg integrand=j scalarize=i dim nout + bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(nout, 0)) + do_tests(; f=bfiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) end From f5648fb908109aa32751aabb18f4693c98edb3d3 Mon Sep 17 00:00:00 2001 From: lxvm Date: Wed, 3 Jan 2024 01:41:17 -0800 Subject: [PATCH 9/9] restore tests --- test/derivative_tests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index f04840e3..061549d8 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -135,7 +135,7 @@ do_tests = function (; f, scalarize, lb, ub, p, alg, abstol, reltol) return end -#= + ### One Dimensional for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution) req.nout > 1 || continue @@ -215,7 +215,7 @@ for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(nout)) do_tests(; f=fiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol) end -=# + ### Batch, One Dimensional for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution) req.allows_batch || continue