diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 00000000..453925c3 --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1 @@ +style = "sciml" \ No newline at end of file diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml new file mode 100644 index 00000000..2a3517a0 --- /dev/null +++ b/.github/workflows/FormatCheck.yml @@ -0,0 +1,42 @@ +name: format-check + +on: + push: + branches: + - 'master' + - 'release-' + tags: '*' + pull_request: + +jobs: + build: + runs-on: ${{ matrix.os }} + strategy: + matrix: + julia-version: [1] + julia-arch: [x86] + os: [ubuntu-latest] + steps: + - uses: julia-actions/setup-julia@latest + with: + version: ${{ matrix.julia-version }} + + - uses: actions/checkout@v1 + - name: Install JuliaFormatter and format + # This will use the latest version by default but you can set the version like so: + # + # julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter", version="0.13.0"))' + run: | + julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))' + julia -e 'using JuliaFormatter; format(".", verbose=true)' + - name: Format check + run: | + julia -e ' + out = Cmd(`git diff --name-only`) |> read |> String + if out == "" + exit(0) + else + @error "Some files have not been formatted !!!" + write(stdout, out) + exit(1) + end' diff --git a/docs/make.jl b/docs/make.jl index e09b99de..c6e93932 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -2,18 +2,14 @@ using Documenter, Integrals include("pages.jl") -makedocs( - sitename="Integrals.jl", - authors="Chris Rackauckas", - modules=[Integrals,Integrals.SciMLBase], - clean=true,doctest=false, - format = Documenter.HTML(analytics = "UA-90474609-3", - assets = ["assets/favicon.ico"], - canonical="https://integrals.sciml.ai/stable/"), - pages=pages -) +makedocs(sitename = "Integrals.jl", + authors = "Chris Rackauckas", + modules = [Integrals, Integrals.SciMLBase], + clean = true, doctest = false, + format = Documenter.HTML(analytics = "UA-90474609-3", + assets = ["assets/favicon.ico"], + canonical = "https://integrals.sciml.ai/stable/"), + pages = pages) -deploydocs( - repo = "github.com/SciML/Integrals.jl.git"; - push_preview = true -) +deploydocs(repo = "github.com/SciML/Integrals.jl.git"; + push_preview = true) diff --git a/docs/pages.jl b/docs/pages.jl index f7d52a7d..d9db8209 100644 --- a/docs/pages.jl +++ b/docs/pages.jl @@ -1,15 +1,9 @@ pages = [ "Home" => "index.md", - "Tutorials" => Any[ - "tutorials/numerical_integrals.md", - "tutorials/differentiating_integrals.md" - ], - "Basics" => Any[ - "basics/IntegralProblem.md", - "basics/solve.md", - "basics/FAQ.md" - ], - "Solvers" => Any[ - "solvers/IntegralSolvers.md" - ] -] \ No newline at end of file + "Tutorials" => Any["tutorials/numerical_integrals.md", + "tutorials/differentiating_integrals.md"], + "Basics" => Any["basics/IntegralProblem.md", + "basics/solve.md", + "basics/FAQ.md"], + "Solvers" => Any["solvers/IntegralSolvers.md"], +] diff --git a/lib/IntegralsCuba/src/IntegralsCuba.jl b/lib/IntegralsCuba/src/IntegralsCuba.jl index 48e2f851..6cdb54d1 100644 --- a/lib/IntegralsCuba/src/IntegralsCuba.jl +++ b/lib/IntegralsCuba/src/IntegralsCuba.jl @@ -9,11 +9,12 @@ struct CubaSUAVE <: AbstractCubaAlgorithm end struct CubaDivonne <: AbstractCubaAlgorithm end struct CubaCuhre <: AbstractCubaAlgorithm end -function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgorithm, sensealg, - lb, ub, p, args...; - reltol=1e-8, abstol=1e-8, - maxiters=alg isa CubaSUAVE ? 1000000 : typemax(Int), - kwargs...) +function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgorithm, + sensealg, + lb, ub, p, args...; + reltol = 1e-8, abstol = 1e-8, + maxiters = alg isa CubaSUAVE ? 1000000 : typemax(Int), + kwargs...) prob = transformation_if_inf(prob) #intercept for infinite transformation p = p if lb isa Number && prob.batch == 0 @@ -36,7 +37,8 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori end else f = function (x, dx) - dx .= prob.f(scale_x!(_x, ub, lb, x), p) .* prod((y) -> y[1] - y[2], zip(ub, lb)) + dx .= prob.f(scale_x!(_x, ub, lb, x), p) .* + prod((y) -> y[1] - y[2], zip(ub, lb)) end end else @@ -50,11 +52,13 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori else if prob.f([lb ub], p) isa Vector f = function (x, dx) - dx .= prob.f(scale_x(ub, lb, x), p)' .* prod((y) -> y[1] - y[2], zip(ub, lb)) + dx .= prob.f(scale_x(ub, lb, x), p)' .* + prod((y) -> y[1] - y[2], zip(ub, lb)) end else f = function (x, dx) - dx .= prob.f(scale_x(ub, lb, x), p) .* prod((y) -> y[1] - y[2], zip(ub, lb)) + dx .= prob.f(scale_x(ub, lb, x), p) .* + prod((y) -> y[1] - y[2], zip(ub, lb)) end end end @@ -67,11 +71,13 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori else if prob.f([lb ub], p) isa Vector f = function (x, dx) - dx .= prob.f(scale_x(ub, lb, x), p)' .* prod((y) -> y[1] - y[2], zip(ub, lb)) + dx .= prob.f(scale_x(ub, lb, x), p)' .* + prod((y) -> y[1] - y[2], zip(ub, lb)) end else f = function (x, dx) - dx .= prob.f(scale_x(ub, lb, x), p) .* prod((y) -> y[1] - y[2], zip(ub, lb)) + dx .= prob.f(scale_x(ub, lb, x), p) .* + prod((y) -> y[1] - y[2], zip(ub, lb)) end end end @@ -83,21 +89,21 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori nvec = prob.batch == 0 ? 1 : prob.batch if alg isa CubaVegas - out = Cuba.vegas(f, ndim, prob.nout; rtol=reltol, - atol=abstol, nvec=nvec, - maxevals=maxiters, kwargs...) + out = Cuba.vegas(f, ndim, prob.nout; rtol = reltol, + atol = abstol, nvec = nvec, + maxevals = maxiters, kwargs...) elseif alg isa CubaSUAVE - out = Cuba.suave(f, ndim, prob.nout; rtol=reltol, - atol=abstol, nvec=nvec, - maxevals=maxiters, kwargs...) + out = Cuba.suave(f, ndim, prob.nout; rtol = reltol, + atol = abstol, nvec = nvec, + maxevals = maxiters, kwargs...) elseif alg isa CubaDivonne - out = Cuba.divonne(f, ndim, prob.nout; rtol=reltol, - atol=abstol, nvec=nvec, - maxevals=maxiters, kwargs...) + out = Cuba.divonne(f, ndim, prob.nout; rtol = reltol, + atol = abstol, nvec = nvec, + maxevals = maxiters, kwargs...) elseif alg isa CubaCuhre - out = Cuba.cuhre(f, ndim, prob.nout; rtol=reltol, - atol=abstol, nvec=nvec, - maxevals=maxiters, kwargs...) + out = Cuba.cuhre(f, ndim, prob.nout; rtol = reltol, + atol = abstol, nvec = nvec, + maxevals = maxiters, kwargs...) end if isinplace(prob) || prob.batch != 0 @@ -111,9 +117,9 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori end SciMLBase.build_solution(prob, alg, val, out.error, - chi=out.probability, retcode=:Success) + chi = out.probability, retcode = :Success) end export CubaVegas, CubaSUAVE, CubaDivonne, CubaCuhre -end \ No newline at end of file +end diff --git a/lib/IntegralsCuba/test/runtests.jl b/lib/IntegralsCuba/test/runtests.jl index e69de29b..8b137891 100644 --- a/lib/IntegralsCuba/test/runtests.jl +++ b/lib/IntegralsCuba/test/runtests.jl @@ -0,0 +1 @@ + diff --git a/lib/IntegralsCubature/src/IntegralsCubature.jl b/lib/IntegralsCubature/src/IntegralsCubature.jl index f496a685..f2108956 100644 --- a/lib/IntegralsCubature/src/IntegralsCubature.jl +++ b/lib/IntegralsCubature/src/IntegralsCubature.jl @@ -9,11 +9,11 @@ struct CubatureJLh <: AbstractCubatureJLAlgorithm end struct CubatureJLp <: AbstractCubatureJLAlgorithm end function Integrals.__solvebp_call(prob::IntegralProblem, - alg::AbstractCubatureJLAlgorithm, - sensealg, lb, ub, p, args...; - reltol=1e-8, abstol=1e-8, - maxiters=typemax(Int), - kwargs...) + alg::AbstractCubatureJLAlgorithm, + sensealg, lb, ub, p, args...; + reltol = 1e-8, abstol = 1e-8, + maxiters = typemax(Int), + kwargs...) prob = transformation_if_inf(prob) #intercept for infinite transformation nout = prob.nout if nout == 1 @@ -27,23 +27,23 @@ function Integrals.__solvebp_call(prob::IntegralProblem, if lb isa Number if alg isa CubatureJLh _val, err = Cubature.hquadrature(f, lb, ub; - reltol=reltol, abstol=abstol, - maxevals=maxiters) + reltol = reltol, abstol = abstol, + maxevals = maxiters) else _val, err = Cubature.pquadrature(f, lb, ub; - reltol=reltol, abstol=abstol, - maxevals=maxiters) + reltol = reltol, abstol = abstol, + maxevals = maxiters) end val = prob.f(lb, p) isa Number ? _val : [_val] else if alg isa CubatureJLh _val, err = Cubature.hcubature(f, lb, ub; - reltol=reltol, abstol=abstol, - maxevals=maxiters) + reltol = reltol, abstol = abstol, + maxevals = maxiters) else _val, err = Cubature.pcubature(f, lb, ub; - reltol=reltol, abstol=abstol, - maxevals=maxiters) + reltol = reltol, abstol = abstol, + maxevals = maxiters) end if isinplace(prob) || !isa(prob.f(lb, p), Number) @@ -75,22 +75,22 @@ function Integrals.__solvebp_call(prob::IntegralProblem, if lb isa Number if alg isa CubatureJLh _val, err = Cubature.hquadrature_v(f, lb, ub; - reltol=reltol, abstol=abstol, - maxevals=maxiters) + reltol = reltol, abstol = abstol, + maxevals = maxiters) else _val, err = Cubature.pquadrature_v(f, lb, ub; - reltol=reltol, abstol=abstol, - maxevals=maxiters) + reltol = reltol, abstol = abstol, + maxevals = maxiters) end else if alg isa CubatureJLh _val, err = Cubature.hcubature_v(f, lb, ub; - reltol=reltol, abstol=abstol, - maxevals=maxiters) + reltol = reltol, abstol = abstol, + maxevals = maxiters) else _val, err = Cubature.pcubature_v(f, lb, ub; - reltol=reltol, abstol=abstol, - maxevals=maxiters) + reltol = reltol, abstol = abstol, + maxevals = maxiters) end end val = _val isa Number ? [_val] : _val @@ -105,22 +105,22 @@ function Integrals.__solvebp_call(prob::IntegralProblem, if lb isa Number if alg isa CubatureJLh val, err = Cubature.hquadrature(nout, f, lb, ub; - reltol=reltol, abstol=abstol, - maxevals=maxiters) + reltol = reltol, abstol = abstol, + maxevals = maxiters) else val, err = Cubature.pquadrature(nout, f, lb, ub; - reltol=reltol, abstol=abstol, - maxevals=maxiters) + reltol = reltol, abstol = abstol, + maxevals = maxiters) end else if alg isa CubatureJLh val, err = Cubature.hcubature(nout, f, lb, ub; - reltol=reltol, abstol=abstol, - maxevals=maxiters) + reltol = reltol, abstol = abstol, + maxevals = maxiters) else val, err = Cubature.pcubature(nout, f, lb, ub; - reltol=reltol, abstol=abstol, - maxevals=maxiters) + reltol = reltol, abstol = abstol, + maxevals = maxiters) end end else @@ -137,29 +137,29 @@ function Integrals.__solvebp_call(prob::IntegralProblem, if lb isa Number if alg isa CubatureJLh val, err = Cubature.hquadrature_v(nout, f, lb, ub; - reltol=reltol, abstol=abstol, - maxevals=maxiters) + reltol = reltol, abstol = abstol, + maxevals = maxiters) else val, err = Cubature.pquadrature_v(nout, f, lb, ub; - reltol=reltol, abstol=abstol, - maxevals=maxiters) + reltol = reltol, abstol = abstol, + maxevals = maxiters) end else if alg isa CubatureJLh val, err = Cubature.hcubature_v(nout, f, lb, ub; - reltol=reltol, abstol=abstol, - maxevals=maxiters) + reltol = reltol, abstol = abstol, + maxevals = maxiters) else val, err = Cubature.pcubature_v(nout, f, lb, ub; - reltol=reltol, abstol=abstol, - maxevals=maxiters) + reltol = reltol, abstol = abstol, + maxevals = maxiters) end end end end - SciMLBase.build_solution(prob, alg, val, err, retcode=:Success) + SciMLBase.build_solution(prob, alg, val, err, retcode = :Success) end export CubatureJLh, CubatureJLp -end \ No newline at end of file +end diff --git a/lib/IntegralsCubature/test/runtests.jl b/lib/IntegralsCubature/test/runtests.jl index e69de29b..8b137891 100644 --- a/lib/IntegralsCubature/test/runtests.jl +++ b/lib/IntegralsCubature/test/runtests.jl @@ -0,0 +1 @@ + diff --git a/src/Integrals.jl b/src/Integrals.jl index 9bf95712..5157412a 100644 --- a/src/Integrals.jl +++ b/src/Integrals.jl @@ -1,8 +1,8 @@ module Integrals -using Reexport, MonteCarloIntegration, QuadGK, HCubature +using Reexport, MonteCarloIntegration, QuadGK, HCubature @reexport using SciMLBase -using Zygote, ReverseDiff, ForwardDiff , LinearAlgebra +using Zygote, ReverseDiff, ForwardDiff, LinearAlgebra import ChainRulesCore import ChainRulesCore: NoTangent @@ -15,7 +15,7 @@ struct VEGAS <: SciMLBase.AbstractIntegralAlgorithm ncalls::Int end -VEGAS(;nbins = 100,ncalls = 1000) = VEGAS(nbins,ncalls) +VEGAS(; nbins = 100, ncalls = 1000) = VEGAS(nbins, ncalls) abstract type QuadSensitivityAlg end struct ReCallVJP{V} @@ -28,239 +28,244 @@ struct ReverseDiffVJP compile::Bool end -function scale_x!(_x,ub,lb,x) +function scale_x!(_x, ub, lb, x) _x .= (ub .- lb) .* x .+ lb _x end -function scale_x(ub,lb,x) +function scale_x(ub, lb, x) (ub .- lb) .* x .+ lb end function v_inf(t) - return t ./ (1 .- t.^2) + return t ./ (1 .- t .^ 2) end - -function v_semiinf(t , a , upto_inf) - if upto_inf == true - return a .+ (t ./ (1 .- t)) - else - return a .+ (t ./ (1 .+ t)) - end +function v_semiinf(t, a, upto_inf) + if upto_inf == true + return a .+ (t ./ (1 .- t)) + else + return a .+ (t ./ (1 .+ t)) + end end -function transfor_inf_number(t , p , f , lb , ub) - if lb == -Inf && ub == Inf - j = (1 .+ t.^2 )/(1 .- t.^2).^2 - return f(v_inf(t) , p)*(j) - elseif lb != -Inf && ub == Inf - a = lb - j = 1 ./ ((1 .- t).^2) - return f(v_semiinf(t , a , 1) , p)*(j) - elseif lb == -Inf && ub != Inf - a = ub - j = 1 ./ ((1 .+ t ).^2) - return f(v_semiinf(t , a , 0) , p)*(j) - end +function transfor_inf_number(t, p, f, lb, ub) + if lb == -Inf && ub == Inf + j = (1 .+ t .^ 2) / (1 .- t .^ 2) .^ 2 + return f(v_inf(t), p) * (j) + elseif lb != -Inf && ub == Inf + a = lb + j = 1 ./ ((1 .- t) .^ 2) + return f(v_semiinf(t, a, 1), p) * (j) + elseif lb == -Inf && ub != Inf + a = ub + j = 1 ./ ((1 .+ t) .^ 2) + return f(v_semiinf(t, a, 0), p) * (j) + end end -function transform_inf(t , p , f , lb , ub) - if lb isa Number && ub isa Number - return transfor_inf_number(t , p , f , lb , ub) - end - - lbb = lb .== -Inf - ubb = ub .== Inf - _none = .!lbb .& .!ubb - _inf = lbb .& ubb - semiup = .!lbb .& ubb - semilw = lbb .& .!ubb - - function v(t) - return t.*_none + v_inf(t).*_inf + v_semiinf(t , lb , 1).*semiup + v_semiinf(t , ub , 0).*semilw - end +function transform_inf(t, p, f, lb, ub) + if lb isa Number && ub isa Number + return transfor_inf_number(t, p, f, lb, ub) + end + + lbb = lb .== -Inf + ubb = ub .== Inf + _none = .!lbb .& .!ubb + _inf = lbb .& ubb + semiup = .!lbb .& ubb + semilw = lbb .& .!ubb + + function v(t) + return t .* _none + v_inf(t) .* _inf + v_semiinf(t, lb, 1) .* semiup + + v_semiinf(t, ub, 0) .* semilw + end jac = Zygote.@ignore ForwardDiff.jacobian(x -> v(x), t) - j = det(jac) - f(v(t) , p)*(j) + j = det(jac) + f(v(t), p) * (j) end function transformation_if_inf(prob) - if (prob.lb isa Number && prob.ub isa Number && (prob.ub == Inf || prob.lb == -Inf)) - if prob.lb == -Inf && prob.ub == Inf - g = prob.f - h(t , p) = transform_inf(t , p , g , prob.lb , prob.ub) - lb = -1.00 - ub = 1.00 - _prob = IntegralProblem(h , lb ,ub, prob.p , nout = prob.nout , batch = prob.batch , prob.kwargs...) - return _prob - elseif prob.lb != -Inf && prob.ub == Inf - g = prob.f - h_(t , p) = transform_inf(t , p , g , prob.lb , prob.ub) - lb = 0.00 - ub = 1.00 - _prob = IntegralProblem(h_ , lb ,ub, prob.p , nout = prob.nout , batch = prob.batch , prob.kwargs...) - return _prob - elseif prob.lb == -Inf && prob.ub != Inf - g = prob.f - _h(t , p) = transform_inf(t , p , g , prob.lb , prob.ub) - lb = -1.00 - ub = 0.00 - _prob = IntegralProblem(_h , lb ,ub, prob.p , nout = prob.nout , batch = prob.batch , prob.kwargs...) - return _prob - end - end - if -Inf in prob.lb || Inf in prob.ub - lbb = prob.lb .== -Inf - ubb = prob.ub .== Inf - _none = .!lbb .& .!ubb - _inf = lbb .& ubb - _semiup = .!lbb .& ubb - _semilw = lbb .& .!ubb - - _lb = 0.00.*_semiup + -1.00.*_inf + -1.00.*_semilw + _none.*prob.lb - _ub = 1.00.*_semiup + 1.00.*_inf + 0.00.*_semilw + _none.*prob.ub - g = prob.f - hs(t , p) = transform_inf(t , p , g , prob.lb , prob.ub) - _prob = IntegralProblem(hs, _lb ,_ub, prob.p , nout = prob.nout , batch = prob.batch , prob.kwargs...) - return _prob - end - return prob + if (prob.lb isa Number && prob.ub isa Number && (prob.ub == Inf || prob.lb == -Inf)) + if prob.lb == -Inf && prob.ub == Inf + g = prob.f + h(t, p) = transform_inf(t, p, g, prob.lb, prob.ub) + lb = -1.00 + ub = 1.00 + _prob = IntegralProblem(h, lb, ub, prob.p, nout = prob.nout, batch = prob.batch, + prob.kwargs...) + return _prob + elseif prob.lb != -Inf && prob.ub == Inf + g = prob.f + h_(t, p) = transform_inf(t, p, g, prob.lb, prob.ub) + lb = 0.00 + ub = 1.00 + _prob = IntegralProblem(h_, lb, ub, prob.p, nout = prob.nout, + batch = prob.batch, prob.kwargs...) + return _prob + elseif prob.lb == -Inf && prob.ub != Inf + g = prob.f + _h(t, p) = transform_inf(t, p, g, prob.lb, prob.ub) + lb = -1.00 + ub = 0.00 + _prob = IntegralProblem(_h, lb, ub, prob.p, nout = prob.nout, + batch = prob.batch, prob.kwargs...) + return _prob + end + end + if -Inf in prob.lb || Inf in prob.ub + lbb = prob.lb .== -Inf + ubb = prob.ub .== Inf + _none = .!lbb .& .!ubb + _inf = lbb .& ubb + _semiup = .!lbb .& ubb + _semilw = lbb .& .!ubb + + _lb = 0.00 .* _semiup + -1.00 .* _inf + -1.00 .* _semilw + _none .* prob.lb + _ub = 1.00 .* _semiup + 1.00 .* _inf + 0.00 .* _semilw + _none .* prob.ub + g = prob.f + hs(t, p) = transform_inf(t, p, g, prob.lb, prob.ub) + _prob = IntegralProblem(hs, _lb, _ub, prob.p, nout = prob.nout, batch = prob.batch, + prob.kwargs...) + return _prob + end + return prob end -function SciMLBase.solve(prob::IntegralProblem,::Nothing,sensealg,lb,ub,p,args...; - reltol = 1e-8, abstol = 1e-8, kwargs...) +function SciMLBase.solve(prob::IntegralProblem, ::Nothing, sensealg, lb, ub, p, args...; + reltol = 1e-8, abstol = 1e-8, kwargs...) if lb isa Number - __solve(prob,QuadGKJL();reltol=reltol,abstol=abstol,kwargs...) + __solve(prob, QuadGKJL(); reltol = reltol, abstol = abstol, kwargs...) elseif length(lb) > 8 && reltol < 1e-4 || abstol < 1e-4 - __solve(prob,VEGAS();reltol=reltol,abstol=abstol,kwargs...) + __solve(prob, VEGAS(); reltol = reltol, abstol = abstol, kwargs...) else - __solve(prob,HCubatureJL();reltol=reltol,abstol=abstol,kwargs...) + __solve(prob, HCubatureJL(); reltol = reltol, abstol = abstol, kwargs...) end end function SciMLBase.solve(prob::IntegralProblem, - alg::SciMLBase.AbstractIntegralAlgorithm, - args...; sensealg = ReCallVJP(ZygoteVJP()), kwargs...) - prob = transformation_if_inf(prob) - __solvebp(prob,alg,sensealg,prob.lb,prob.ub,prob.p,args...;kwargs...) + alg::SciMLBase.AbstractIntegralAlgorithm, + args...; sensealg = ReCallVJP(ZygoteVJP()), kwargs...) + prob = transformation_if_inf(prob) + __solvebp(prob, alg, sensealg, prob.lb, prob.ub, prob.p, args...; kwargs...) end # Give a layer to intercept with AD -__solvebp(args...;kwargs...) = __solvebp_call(args...;kwargs...) +__solvebp(args...; kwargs...) = __solvebp_call(args...; kwargs...) -function __solvebp_call(prob::IntegralProblem,::QuadGKJL,sensealg,lb,ub,p,args...; - reltol = 1e-8, abstol = 1e-8, - maxiters = typemax(Int), - kwargs...) +function __solvebp_call(prob::IntegralProblem, ::QuadGKJL, sensealg, lb, ub, p, args...; + reltol = 1e-8, abstol = 1e-8, + maxiters = typemax(Int), + kwargs...) if isinplace(prob) || lb isa AbstractArray || ub isa AbstractArray error("QuadGKJL only accepts one-dimensional quadrature problems.") end @assert prob.batch == 0 @assert prob.nout == 1 p = p - f = x -> prob.f(x,p) - val,err = quadgk(f, lb, ub, - rtol=reltol, atol=abstol, - kwargs...) - SciMLBase.build_solution(prob,QuadGKJL(),val,err,retcode = :Success) + f = x -> prob.f(x, p) + val, err = quadgk(f, lb, ub, + rtol = reltol, atol = abstol, + kwargs...) + SciMLBase.build_solution(prob, QuadGKJL(), val, err, retcode = :Success) end -function __solvebp_call(prob::IntegralProblem,::HCubatureJL,sensealg,lb,ub,p,args...; - reltol = 1e-8, abstol = 1e-8, - maxiters = typemax(Int), - kwargs...) +function __solvebp_call(prob::IntegralProblem, ::HCubatureJL, sensealg, lb, ub, p, args...; + reltol = 1e-8, abstol = 1e-8, + maxiters = typemax(Int), + kwargs...) p = p if isinplace(prob) dx = zeros(prob.nout) - f = (x) -> (prob.f(dx,x,p); dx) + f = (x) -> (prob.f(dx, x, p); dx) else - f = (x) -> prob.f(x,p) + f = (x) -> prob.f(x, p) end @assert prob.batch == 0 if lb isa Number - val,err = hquadrature(f, lb, ub; - rtol=reltol, atol=abstol, - maxevals=maxiters, kwargs...) + val, err = hquadrature(f, lb, ub; + rtol = reltol, atol = abstol, + maxevals = maxiters, kwargs...) else - val,err = hcubature(f, lb, ub; - rtol=reltol, atol=abstol, - maxevals=maxiters, kwargs...) + val, err = hcubature(f, lb, ub; + rtol = reltol, atol = abstol, + maxevals = maxiters, kwargs...) end - SciMLBase.build_solution(prob,HCubatureJL(),val,err,retcode = :Success) + SciMLBase.build_solution(prob, HCubatureJL(), val, err, retcode = :Success) end -function __solvebp_call(prob::IntegralProblem,alg::VEGAS,sensealg,lb,ub,p,args...; - reltol = 1e-8, abstol = 1e-8, - maxiters = typemax(Int), - kwargs...) +function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, lb, ub, p, args...; + reltol = 1e-8, abstol = 1e-8, + maxiters = typemax(Int), + kwargs...) p = p @assert prob.nout == 1 if prob.batch == 0 if isinplace(prob) - dx = zeros(prob.nout) - f = (x) -> (prob.f(dx,x,p); dx) + dx = zeros(prob.nout) + f = (x) -> (prob.f(dx, x, p); dx) else - f = (x) -> prob.f(x,p) + f = (x) -> prob.f(x, p) end else if isinplace(prob) - dx = zeros(prob.batch) - f = (x) -> (prob.f(dx,x',p); dx) + dx = zeros(prob.batch) + f = (x) -> (prob.f(dx, x', p); dx) else - f = (x) -> prob.f(x',p) + f = (x) -> prob.f(x', p) end end - val,err,chi = vegas(f, lb, ub, rtol=reltol, atol=abstol, - maxiter = maxiters, nbins = alg.nbins, - ncalls = alg.ncalls, batch=prob.batch != 0, kwargs...) - SciMLBase.build_solution(prob,alg,val,err,chi=chi,retcode = :Success) + val, err, chi = vegas(f, lb, ub, rtol = reltol, atol = abstol, + maxiter = maxiters, nbins = alg.nbins, + ncalls = alg.ncalls, batch = prob.batch != 0, kwargs...) + SciMLBase.build_solution(prob, alg, val, err, chi = chi, retcode = :Success) end -function ChainRulesCore.rrule(::typeof(__solvebp),prob,alg,sensealg,lb,ub,p,args...;kwargs...) - out = __solvebp_call(prob,alg,sensealg,lb,ub,p,args...;kwargs...) +function ChainRulesCore.rrule(::typeof(__solvebp), prob, alg, sensealg, lb, ub, p, args...; + kwargs...) + out = __solvebp_call(prob, alg, sensealg, lb, ub, p, args...; kwargs...) function quadrature_adjoint(Δ) - y = typeof(Δ) <: Array{<:Number,0} ? Δ[1] : Δ + y = typeof(Δ) <: Array{<:Number, 0} ? Δ[1] : Δ if isinplace(prob) dx = zeros(prob.nout) - _f = (x) -> prob.f(dx,x,p) + _f = (x) -> prob.f(dx, x, p) if sensealg.vjp isa ZygoteVJP - dfdp = function (dx,x,p) - _,back = Zygote.pullback(p) do p - _dx = Zygote.Buffer(x, prob.nout, size(x,2)) - prob.f(_dx,x,p) + dfdp = function (dx, x, p) + _, back = Zygote.pullback(p) do p + _dx = Zygote.Buffer(x, prob.nout, size(x, 2)) + prob.f(_dx, x, p) copy(_dx) end - z = zeros(size(x,2)) - for idx in 1:size(x,2) + z = zeros(size(x, 2)) + for idx in 1:size(x, 2) z[1] = 1 - dx[:,idx] = back(z)[1] - z[idx]=0 + dx[:, idx] = back(z)[1] + z[idx] = 0 end end elseif sensealg.vjp isa ReverseDiffVJP error("TODO") end else - _f = (x) -> prob.f(x,p) + _f = (x) -> prob.f(x, p) if sensealg.vjp isa ZygoteVJP if prob.batch > 0 - dfdp = function (x,p) - _,back = Zygote.pullback(p->prob.f(x,p),p) + dfdp = function (x, p) + _, back = Zygote.pullback(p -> prob.f(x, p), p) - out = zeros(length(p),size(x,2)) - z = zeros(size(x,2)) - for idx in 1:size(x,2) + out = zeros(length(p), size(x, 2)) + z = zeros(size(x, 2)) + for idx in 1:size(x, 2) z[idx] = 1 - out[:,idx] = back(z)[1] - z[idx]=0 + out[:, idx] = back(z)[1] + z[idx] = 0 end out end else - dfdp = function (x,p) - _,back = Zygote.pullback(p->prob.f(x,p),p) + dfdp = function (x, p) + _, back = Zygote.pullback(p -> prob.f(x, p), p) back(y)[1] end end @@ -270,81 +275,90 @@ function ChainRulesCore.rrule(::typeof(__solvebp),prob,alg,sensealg,lb,ub,p,args end end - dp_prob = IntegralProblem(dfdp,lb,ub,p;nout=length(p),batch = prob.batch,kwargs...) + dp_prob = IntegralProblem(dfdp, lb, ub, p; nout = length(p), batch = prob.batch, + kwargs...) if p isa Number - dp = __solvebp_call(dp_prob,alg,sensealg,lb,ub,p,args...;kwargs...)[1] + dp = __solvebp_call(dp_prob, alg, sensealg, lb, ub, p, args...; kwargs...)[1] else - dp = __solvebp_call(dp_prob,alg,sensealg,lb,ub,p,args...;kwargs...).u + dp = __solvebp_call(dp_prob, alg, sensealg, lb, ub, p, args...; kwargs...).u end if lb isa Number dlb = -_f(lb) dub = _f(ub) - return (NoTangent(),NoTangent(),NoTangent(),NoTangent(),dlb,dub,dp,ntuple(x->NoTangent(),length(args))...) + return (NoTangent(), NoTangent(), NoTangent(), NoTangent(), dlb, dub, dp, + ntuple(x -> NoTangent(), length(args))...) else - return (NoTangent(),NoTangent(),NoTangent(),NoTangent(),NoTangent(),NoTangent(),dp,ntuple(x->NoTangent(),length(args))...) + return (NoTangent(), NoTangent(), NoTangent(), NoTangent(), NoTangent(), + NoTangent(), dp, ntuple(x -> NoTangent(), length(args))...) end end - out,quadrature_adjoint + out, quadrature_adjoint end -ZygoteRules.@adjoint function ZygoteRules.literal_getproperty( - sol::SciMLBase.IntegralSolution,::Val{:u}) - sol.u,Δ->(SciMLBase.build_solution(sol.prob,sol.alg,Δ,sol.resid),) +ZygoteRules.@adjoint function ZygoteRules.literal_getproperty(sol::SciMLBase.IntegralSolution, + ::Val{:u}) + sol.u, Δ -> (SciMLBase.build_solution(sol.prob, sol.alg, Δ, sol.resid),) end - ### Forward-Mode AD Intercepts # Direct AD on solvers with QuadGK and HCubature -function __solvebp(prob,alg::QuadGKJL,sensealg,lb,ub,p::AbstractArray{<:ForwardDiff.Dual{T,V,P},N},args...;kwargs...) where {T,V,P,N} - __solvebp_call(prob,alg,sensealg,lb,ub,p,args...;kwargs...) +function __solvebp(prob, alg::QuadGKJL, sensealg, lb, ub, + p::AbstractArray{<:ForwardDiff.Dual{T, V, P}, N}, args...; + kwargs...) where {T, V, P, N} + __solvebp_call(prob, alg, sensealg, lb, ub, p, args...; kwargs...) end -function __solvebp(prob,alg::HCubatureJL,sensealg,lb,ub,p::AbstractArray{<:ForwardDiff.Dual{T,V,P},N},args...;kwargs...) where {T,V,P,N} - __solvebp_call(prob,alg,sensealg,lb,ub,p,args...;kwargs...) +function __solvebp(prob, alg::HCubatureJL, sensealg, lb, ub, + p::AbstractArray{<:ForwardDiff.Dual{T, V, P}, N}, args...; + kwargs...) where {T, V, P, N} + __solvebp_call(prob, alg, sensealg, lb, ub, p, args...; kwargs...) end # Manually split for the pushforward -function __solvebp(prob,alg,sensealg,lb,ub,p::AbstractArray{<:ForwardDiff.Dual{T,V,P},N},args...;kwargs...) where {T,V,P,N} - primal = __solvebp_call(prob,alg,sensealg,lb,ub,ForwardDiff.value.(p),args...;kwargs...) +function __solvebp(prob, alg, sensealg, lb, ub, + p::AbstractArray{<:ForwardDiff.Dual{T, V, P}, N}, args...; + kwargs...) where {T, V, P, N} + primal = __solvebp_call(prob, alg, sensealg, lb, ub, ForwardDiff.value.(p), args...; + kwargs...) - nout = prob.nout*P + nout = prob.nout * P if isinplace(prob) - 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) if prob.batch > 0 - dx = similar(dualp, prob.nout, size(x,2)) + dx = similar(dualp, prob.nout, size(x, 2)) else dx = similar(dualp, prob.nout) end - prob.f(dx,x,dualp) + prob.f(dx, x, dualp) - ys = reinterpret(ForwardDiff.Dual{T,V,P}, dx) - idx=0 + ys = reinterpret(ForwardDiff.Dual{T, V, P}, dx) + idx = 0 for y in ys for p in ForwardDiff.partials(y) - out[idx+=1] = p + out[idx += 1] = p end end return out end else - dfdp = function (x,p) - dualp = reinterpret(ForwardDiff.Dual{T,V,P}, p) - ys = prob.f(x,dualp) + dfdp = function (x, p) + dualp = reinterpret(ForwardDiff.Dual{T, V, P}, p) + ys = prob.f(x, dualp) if prob.batch > 0 - out = similar(p, V, nout, size(x,2)) + out = similar(p, V, nout, size(x, 2)) else out = similar(p, V, nout) end - idx=0 + idx = 0 for y in ys for p in ForwardDiff.partials(y) - out[idx+=1] = p + out[idx += 1] = p end end @@ -353,17 +367,18 @@ function __solvebp(prob,alg,sensealg,lb,ub,p::AbstractArray{<:ForwardDiff.Dual{T end rawp = copy(reinterpret(V, p)) - dp_prob = IntegralProblem(dfdp,lb,ub,rawp;nout=nout,batch = prob.batch,kwargs...) - dual = __solvebp_call(dp_prob,alg,sensealg,lb,ub,rawp,args...;kwargs...) + dp_prob = IntegralProblem(dfdp, lb, ub, rawp; nout = nout, batch = prob.batch, + kwargs...) + dual = __solvebp_call(dp_prob, alg, sensealg, lb, ub, rawp, args...; kwargs...) res = similar(p, prob.nout) partials = reinterpret(typeof(first(res).partials), dual.u) for idx in eachindex(res) - res[idx] = ForwardDiff.Dual{T,V,P}(primal.u[idx], partials[idx]) + res[idx] = ForwardDiff.Dual{T, V, P}(primal.u[idx], partials[idx]) end if primal.u isa Number res = first(res) end - SciMLBase.build_solution(prob,alg,res,primal.resid) + SciMLBase.build_solution(prob, alg, res, primal.resid) end export QuadGKJL, HCubatureJL, VEGAS diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl index a83e3c82..cdef9fcd 100644 --- a/test/derivative_tests.jl +++ b/test/derivative_tests.jl @@ -1,245 +1,243 @@ -using Integrals, Zygote, FiniteDiff, ForwardDiff, DiffEqSensitivity -using IntegralsCuba, IntegralsCubature -using Test - -### One Dimensional - -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]) -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) - -#@test dlb1 ≈ dlb3 -#@test dub1 ≈ dub3 -@test dp1 ≈ dp3 - -### N-dimensional - -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]) -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) - -function testf(lb,ub,p) - prob = IntegralProblem(f,lb,ub,p,nout=2) - sin(sum(solve(prob,CubaCuhre(),reltol=1e-6,abstol=1e-6))) -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] - -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=1) - 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] # 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_broken dp2 ≈ dp3 #passes - -### 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)) -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 - -### Batch multi dim -f(x,p) = x[1,:]*p[1].+p[2]*p[3] - -lb =[1.0,1.0] -ub = [3.0,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=1) - 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) - -@test dp1 ≈ dp3 -@test dp2 ≈ dp3 - -### Batch multi dim nout -f(x,p) = x*p[1].+p[2]*p[3] - -lb =[1.0,1.0] -ub = [3.0,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)) -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 - - -## iip Batch mulit dim -function g(dx, x,p) - dx.= sum(x*p[1].+p[2]*p[3], dims=1) -end - -lb =[1.0,1.0] -ub = [3.0,3.0] -p = [2.0, 3.0, 4.0] -prob = IntegralProblem(g,lb,ub,p) - -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] -end - -testf3(lb,ub,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 - -## iip Batch mulit dim nout -function g(dx, x,p) - dx.= x*p[1].+p[2]*p[3] -end - -lb =[1.0,1.0] -ub = [3.0,3.0] -p = [2.0, 3.0, 4.0] -prob = IntegralProblem(g,lb,ub,p) - -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)) -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 +using Integrals, Zygote, FiniteDiff, ForwardDiff, DiffEqSensitivity +using IntegralsCuba, IntegralsCubature +using Test + +### One Dimensional + +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]) +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) + +#@test dlb1 ≈ dlb3 +#@test dub1 ≈ dub3 +@test dp1 ≈ dp3 + +### N-dimensional + +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]) +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) + +function testf(lb, ub, p) + prob = IntegralProblem(f, lb, ub, p, nout = 2) + sin(sum(solve(prob, CubaCuhre(), reltol = 1e-6, abstol = 1e-6))) +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] + +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 = 1) + 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] # 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_broken dp2 ≈ dp3 #passes + +### 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)) +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 + +### Batch multi dim +f(x, p) = x[1, :] * p[1] .+ p[2] * p[3] + +lb = [1.0, 1.0] +ub = [3.0, 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 = 1) + 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) + +@test dp1 ≈ dp3 +@test dp2 ≈ dp3 + +### Batch multi dim nout +f(x, p) = x * p[1] .+ p[2] * p[3] + +lb = [1.0, 1.0] +ub = [3.0, 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)) +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 + +## iip Batch mulit dim +function g(dx, x, p) + dx .= sum(x * p[1] .+ p[2] * p[3], dims = 1) +end + +lb = [1.0, 1.0] +ub = [3.0, 3.0] +p = [2.0, 3.0, 4.0] +prob = IntegralProblem(g, lb, ub, p) + +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] +end + +testf3(lb, ub, 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 + +## iip Batch mulit dim nout +function g(dx, x, p) + dx .= x * p[1] .+ p[2] * p[3] +end + +lb = [1.0, 1.0] +ub = [3.0, 3.0] +p = [2.0, 3.0, 4.0] +prob = IntegralProblem(g, lb, ub, p) + +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)) +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 diff --git a/test/inf_integral_tests.jl b/test/inf_integral_tests.jl index aa2f9f2a..60561ec7 100644 --- a/test/inf_integral_tests.jl +++ b/test/inf_integral_tests.jl @@ -1,33 +1,32 @@ -using Integrals , Distributions , Test +using Integrals, Distributions, Test -μ = [0.00 , 0.00] -Σ = [0.4 0.0 ; 0.00 0.4] -d = MvNormal(μ , Σ) -m(x , p) = pdf(d, x) -prob = IntegralProblem(m, [-Inf , -Inf] , [Inf , Inf]) -sol = solve(prob,HCubatureJL(),reltol=1e-3,abstol=1e-3) -@test (1.00 - sol.u)^2 < 1e-6 +μ = [0.00, 0.00] +Σ = [0.4 0.0; 0.00 0.4] +d = MvNormal(μ, Σ) +m(x, p) = pdf(d, x) +prob = IntegralProblem(m, [-Inf, -Inf], [Inf, Inf]) +sol = solve(prob, HCubatureJL(), reltol = 1e-3, abstol = 1e-3) +@test (1.00 - sol.u)^2 < 1e-6 -μ = [0.00 , 0.00] -Σ = [0.4 0.0 ; 0.00 0.4] -d = MvNormal(μ , Σ) -m(x , p) = pdf(d, x) -prob = IntegralProblem(m, [-Inf , 0.00] , [Inf , Inf]) -sol = solve(prob,HCubatureJL(),reltol=1e-3,abstol=1e-3) -@test (0.500 - sol.u)^2 < 1e-6 +μ = [0.00, 0.00] +Σ = [0.4 0.0; 0.00 0.4] +d = MvNormal(μ, Σ) +m(x, p) = pdf(d, x) +prob = IntegralProblem(m, [-Inf, 0.00], [Inf, Inf]) +sol = solve(prob, HCubatureJL(), reltol = 1e-3, abstol = 1e-3) +@test (0.500 - sol.u)^2 < 1e-6 -f(x , p) = pdf(Normal(0.00 , 1.00) , x) -prob = IntegralProblem(f, -Inf , Inf) -sol = solve(prob,HCubatureJL(),reltol=1e-3,abstol=1e-3) -@test (1.00 - sol.u)^2 < 1e-6 +f(x, p) = pdf(Normal(0.00, 1.00), x) +prob = IntegralProblem(f, -Inf, Inf) +sol = solve(prob, HCubatureJL(), reltol = 1e-3, abstol = 1e-3) +@test (1.00 - sol.u)^2 < 1e-6 +f(x, p) = pdf(Normal(0.00, 1.00), x) +prob = IntegralProblem(f, -Inf, 0.00) +sol = solve(prob, HCubatureJL(), reltol = 1e-3, abstol = 1e-3) +@test (0.50 - sol.u)^2 < 1e-6 -f(x , p) = pdf(Normal(0.00 , 1.00) , x) -prob = IntegralProblem(f, -Inf , 0.00) -sol = solve(prob,HCubatureJL(),reltol=1e-3,abstol=1e-3) -@test (0.50 - sol.u)^2 < 1e-6 - -f(x , p) = (1 / (x^2 + 1)) -prob = IntegralProblem(f, 0 , Inf) -sol = solve(prob,HCubatureJL(),reltol=1e-3,abstol=1e-3) -@test (pi/2 - sol.u)^2 < 1e-6 +f(x, p) = (1 / (x^2 + 1)) +prob = IntegralProblem(f, 0, Inf) +sol = solve(prob, HCubatureJL(), reltol = 1e-3, abstol = 1e-3) +@test (pi / 2 - sol.u)^2 < 1e-6 diff --git a/test/interface_tests.jl b/test/interface_tests.jl index f35d5a91..9138bc8d 100644 --- a/test/interface_tests.jl +++ b/test/interface_tests.jl @@ -1,310 +1,324 @@ -using Integrals -using IntegralsCuba, IntegralsCubature -using Test - -max_dim_test = 2 -max_nout_test = 2 -reltol=1e-3 -abstol=1e-3 - -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=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 )) - -integrands = [ - (x,p) -> 1.0, - (x,p) -> x isa Number ? cos(x) : prod(cos.(x)) - ] -iip_integrands = [ (dx,x,p)-> (dx .= f(x,p)) for f ∈ integrands] - -integrands_v = [ - (x,p,nout) -> collect(1.0:nout) - (x,p,nout) -> integrands[2](x,p)*collect(1.0:nout) - ] -iip_integrands_v = [ (dx,x,p,nout)-> (dx .= f(x,p,nout)) for f ∈ integrands_v] - -exact_sol = [ - (ndim, nout, lb, ub) -> prod(ub-lb), - (ndim, nout, lb, ub) -> prod(sin.(ub)-sin.(lb)) - ] - -exact_sol_v = [ - (ndim, nout, lb, ub) -> prod(ub-lb) * collect(1.0:nout), - (ndim, nout, lb, ub) -> exact_sol[2](ndim,nout,lb,ub) * collect(1:nout) - ] - -batch_f(f) = (pts,p) -> begin - fevals = zeros(size(pts,2)) - for i = 1:size(pts, 2) - x = pts[:,i] - fevals[i] = f(x,p) - end - fevals -end - -batch_iip_f(f) = (fevals,pts,p) -> begin - for i = 1:size(pts, 2) - x = pts[:,i] - fevals[i] = f(x,p) - end - nothing -end - -batch_f_v(f, nout) = (pts,p) -> begin - fevals = zeros(nout, size(pts,2)) - for i = 1:size(pts, 2) - x = pts[:,i] - fevals[:,i] = f(x,p,nout) - end - fevals -end - -batch_iip_f_v(f,nout) = (fevals,pts,p) -> begin - for i = 1:size(pts, 2) - x = pts[:,i] - fevals[:,i] = f(x,p, nout) - end - nothing -end - -@testset "Standard Single Dimension Integrands" begin - lb,ub = (1.0,3.0) - nout = 1 - dim = 1 - for alg in algs - if alg_req[alg].min_dim > 1 - continue - 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) - @test sol.u ≈ exact_sol[i](dim,nout,lb,ub) rtol = 1e-2 - end - end -end - -@testset "Standard Integrands" begin - nout = 1 - for alg in algs - req = alg_req[alg] - for i in 1:length(integrands) - for dim = 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 - continue - end - @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob,alg,reltol=reltol,abstol=abstol) - @test sol.u ≈ exact_sol[i](dim,nout,lb,ub) rtol = 1e-2 - end - end - end -end - -@testset "In-place Standard Integrands" begin - nout = 1 - for alg in algs - req = alg_req[alg] - for i in 1:length(iip_integrands) - for dim = 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 - 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) - else - sol = solve(prob,alg,reltol=reltol,abstol=abstol) - end - @test sol.u ≈ [exact_sol[i](dim,nout,lb,ub)] rtol = 1e-2 - end - end - end -end - -@testset "Batched Single Dimension Integrands" begin - (lb,ub) = (1.0,3.0) - (dim, nout) = (1,1) - for alg in algs - req = alg_req[alg] - for i in 1:length(integrands) - prob = IntegralProblem(batch_f(integrands[i]),lb,ub,batch=10) - if req.min_dim > 1 || !req.allows_batch - continue - end - @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob,alg,reltol=reltol,abstol=abstol) - @test sol.u ≈ [exact_sol[i](dim,nout,lb,ub)] rtol = 1e-2 - end - end -end - -@testset "Batched Standard Integrands" begin - nout = 1 - for alg in algs - req = alg_req[alg] - for i in 1:length(integrands) - for dim = 1:max_dim_test - (lb,ub) = (ones(dim),3ones(dim)) - prob = IntegralProblem(batch_f(integrands[i]),lb,ub,batch=10) - if dim > req.max_dim || dim < req.min_dim || !req.allows_batch - continue - end - @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob,alg,reltol=reltol,abstol=abstol) - @test sol.u ≈ [exact_sol[i](dim,nout,lb,ub)] rtol = 1e-2 - end - end - end -end - -@testset "In-Place Batched Standard Integrands" begin - nout = 1 - for alg in algs - req = req = alg_req[alg] - for i in 1:length(iip_integrands) - for dim = 1:max_dim_test - (lb,ub) = (ones(dim),3ones(dim)) - prob = IntegralProblem(batch_iip_f(integrands[i]),lb,ub,batch=10) - if dim > req.max_dim || dim < req.min_dim || !req.allows_batch || !req.allows_iip - continue - end - @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob,alg,reltol=reltol,abstol=abstol) - @test sol.u ≈ [exact_sol[i](dim,nout,lb,ub)] rtol = 1e-2 - end - end - end -end - -######## Vector Valued Integrands -@testset "Standard Single Dimension Vector Integrands" begin - lb,ub = (1.0,3.0) - dim = 1 - for alg in algs - req = alg_req[alg] - for i in 1:length(integrands_v) - for nout = 1:max_nout_test - prob = IntegralProblem((x,p) -> integrands_v[i](x,p,nout),lb,ub, nout = nout) - if req.min_dim > 1 || req.nout < nout - continue - end - @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob,alg,reltol=reltol,abstol=abstol) - @test sol.u ≈ exact_sol_v[i](dim,nout,lb,ub) rtol = 1e-2 - end - end - end -end - -@testset "Standard Vector Integrands" begin - for alg in algs - req = alg_req[alg] - for i in 1:length(integrands_v) - for dim = 1:max_dim_test - lb, ub = (ones(dim), 3ones(dim)) - for nout = 1:max_nout_test - if dim > req.max_dim || dim < req.min_dim || req.nout < nout || alg isa QuadGKJL #QuadGKJL requires 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) - @test sol.u ≈ exact_sol_v[i](dim,nout,lb,ub) rtol = 1e-2 - end - end - end - end -end - -@testset "In-place Standard Vector Integrands" begin - for alg in algs - req = alg_req[alg] - for i in 1:length(iip_integrands_v) - for dim = 1:max_dim_test - lb, ub = (ones(dim), 3ones(dim)) - for nout = 1:max_nout_test - 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 - 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) - else - sol = solve(prob,alg,reltol=reltol,abstol=abstol) - end - @test sol.u ≈ exact_sol_v[i](dim,nout,lb,ub) rtol = 1e-2 - end - end - end - end -end - -@testset "Batched Single Dimension Vector Integrands" begin - (lb,ub) = (1.0,3.0) - (dim, nout) = (1,2) - for alg in algs - req = alg_req[alg] - for i in 1:length(integrands_v) - prob = IntegralProblem(batch_f_v(integrands_v[i],nout),lb,ub,batch=10,nout = nout) - if req.min_dim > 1 || !req.allows_batch || req.nout < nout - continue - end - @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob,alg,reltol=reltol,abstol=abstol) - @test sol.u ≈ exact_sol_v[i](dim,nout,lb,ub) rtol = 1e-2 - end - end -end - -@testset "Batched Standard Vector Integrands" begin - nout = 2 - for alg in algs - req = alg_req[alg] - for i in 1:length(integrands_v) - for dim = 1:max_dim_test - (lb,ub) = (ones(dim),3ones(dim)) - prob = IntegralProblem(batch_f_v(integrands_v[i],nout),lb,ub,batch=10,nout = nout) - if dim > req.max_dim || dim < req.min_dim || !req.allows_batch || req.nout < nout - continue - end - @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob,alg,reltol=reltol,abstol=abstol) - @test sol.u ≈ exact_sol_v[i](dim,nout,lb,ub) rtol = 1e-2 - end - end - end -end - -@testset "In-Place Batched Standard Vector Integrands" begin - nout = 2 - for alg in algs - req = req = alg_req[alg] - for i in 1:length(iip_integrands_v) - for dim = 1:max_dim_test - (lb,ub) = (ones(dim),3ones(dim)) - prob = IntegralProblem(batch_iip_f_v(integrands_v[i],nout),lb,ub,batch=10,nout = nout) - if dim > req.max_dim || dim < req.min_dim || !req.allows_batch || !req.allows_iip || req.nout < nout - continue - end - @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout" - sol = solve(prob,alg,reltol=reltol,abstol=abstol) - @test sol.u ≈ exact_sol_v[i](dim,nout,lb,ub) rtol = 1e-2 - end - end - end -end +using Integrals +using IntegralsCuba, IntegralsCubature +using Test + +max_dim_test = 2 +max_nout_test = 2 +reltol = 1e-3 +abstol = 1e-3 + +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 = 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)) + +integrands = [ + (x, p) -> 1.0, + (x, p) -> x isa Number ? cos(x) : prod(cos.(x)), +] +iip_integrands = [(dx, x, p) -> (dx .= f(x, p)) for f in integrands] + +integrands_v = [(x, p, nout) -> collect(1.0:nout) + (x, p, nout) -> integrands[2](x, p) * collect(1.0:nout)] +iip_integrands_v = [(dx, x, p, nout) -> (dx .= f(x, p, nout)) for f in integrands_v] + +exact_sol = [ + (ndim, nout, lb, ub) -> prod(ub - lb), + (ndim, nout, lb, ub) -> prod(sin.(ub) - sin.(lb)), +] + +exact_sol_v = [ + (ndim, nout, lb, ub) -> prod(ub - lb) * collect(1.0:nout), + (ndim, nout, lb, ub) -> exact_sol[2](ndim, nout, lb, ub) * collect(1:nout), +] + +batch_f(f) = (pts, p) -> begin + fevals = zeros(size(pts, 2)) + for i in 1:size(pts, 2) + x = pts[:, i] + fevals[i] = f(x, p) + end + fevals +end + +batch_iip_f(f) = (fevals, pts, p) -> begin + for i in 1:size(pts, 2) + x = pts[:, i] + fevals[i] = f(x, p) + end + nothing +end + +batch_f_v(f, nout) = (pts, p) -> begin + fevals = zeros(nout, size(pts, 2)) + for i in 1:size(pts, 2) + x = pts[:, i] + fevals[:, i] = f(x, p, nout) + end + fevals +end + +batch_iip_f_v(f, nout) = (fevals, pts, p) -> begin + for i in 1:size(pts, 2) + x = pts[:, i] + fevals[:, i] = f(x, p, nout) + end + nothing +end + +@testset "Standard Single Dimension Integrands" begin + lb, ub = (1.0, 3.0) + nout = 1 + dim = 1 + for alg in algs + if alg_req[alg].min_dim > 1 + continue + 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) + @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2 + end + end +end + +@testset "Standard Integrands" begin + nout = 1 + for alg in algs + req = alg_req[alg] + for i in 1:length(integrands) + 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 + continue + end + @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout" + sol = solve(prob, alg, reltol = reltol, abstol = abstol) + @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2 + end + end + end +end + +@testset "In-place Standard Integrands" begin + nout = 1 + for alg in algs + req = alg_req[alg] + for i in 1:length(iip_integrands) + 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 + 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) + else + sol = solve(prob, alg, reltol = reltol, abstol = abstol) + end + @test sol.u≈[exact_sol[i](dim, nout, lb, ub)] rtol=1e-2 + end + end + end +end + +@testset "Batched Single Dimension Integrands" begin + (lb, ub) = (1.0, 3.0) + (dim, nout) = (1, 1) + for alg in algs + req = alg_req[alg] + for i in 1:length(integrands) + prob = IntegralProblem(batch_f(integrands[i]), lb, ub, batch = 10) + if req.min_dim > 1 || !req.allows_batch + continue + end + @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout" + sol = solve(prob, alg, reltol = reltol, abstol = abstol) + @test sol.u≈[exact_sol[i](dim, nout, lb, ub)] rtol=1e-2 + end + end +end + +@testset "Batched Standard Integrands" begin + nout = 1 + for alg in algs + req = alg_req[alg] + for i in 1:length(integrands) + for dim in 1:max_dim_test + (lb, ub) = (ones(dim), 3ones(dim)) + prob = IntegralProblem(batch_f(integrands[i]), lb, ub, batch = 10) + if dim > req.max_dim || dim < req.min_dim || !req.allows_batch + continue + end + @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout" + sol = solve(prob, alg, reltol = reltol, abstol = abstol) + @test sol.u≈[exact_sol[i](dim, nout, lb, ub)] rtol=1e-2 + end + end + end +end + +@testset "In-Place Batched Standard Integrands" begin + nout = 1 + for alg in algs + req = req = alg_req[alg] + for i in 1:length(iip_integrands) + for dim in 1:max_dim_test + (lb, ub) = (ones(dim), 3ones(dim)) + prob = IntegralProblem(batch_iip_f(integrands[i]), lb, ub, batch = 10) + if dim > req.max_dim || dim < req.min_dim || !req.allows_batch || + !req.allows_iip + continue + end + @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout" + sol = solve(prob, alg, reltol = reltol, abstol = abstol) + @test sol.u≈[exact_sol[i](dim, nout, lb, ub)] rtol=1e-2 + end + end + end +end + +######## Vector Valued Integrands +@testset "Standard Single Dimension Vector Integrands" begin + lb, ub = (1.0, 3.0) + dim = 1 + for alg in algs + req = alg_req[alg] + for i in 1:length(integrands_v) + for nout in 1:max_nout_test + prob = IntegralProblem((x, p) -> integrands_v[i](x, p, nout), lb, ub, + nout = nout) + if req.min_dim > 1 || req.nout < nout + continue + end + @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout" + sol = solve(prob, alg, reltol = reltol, abstol = abstol) + @test sol.u≈exact_sol_v[i](dim, nout, lb, ub) rtol=1e-2 + end + end + end +end + +@testset "Standard Vector Integrands" begin for alg in algs + req = alg_req[alg] + for i in 1:length(integrands_v) + for dim in 1:max_dim_test + 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 #QuadGKJL requires 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) + @test sol.u≈exact_sol_v[i](dim, nout, lb, ub) rtol=1e-2 + end + end + end +end end + +@testset "In-place Standard Vector Integrands" begin for alg in algs + req = alg_req[alg] + for i in 1:length(iip_integrands_v) + for dim in 1:max_dim_test + lb, ub = (ones(dim), 3ones(dim)) + for nout in 1:max_nout_test + 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 + 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) + else + sol = solve(prob, alg, reltol = reltol, abstol = abstol) + end + @test sol.u≈exact_sol_v[i](dim, nout, lb, ub) rtol=1e-2 + end + end + end +end end + +@testset "Batched Single Dimension Vector Integrands" begin + (lb, ub) = (1.0, 3.0) + (dim, nout) = (1, 2) + for alg in algs + req = alg_req[alg] + for i in 1:length(integrands_v) + prob = IntegralProblem(batch_f_v(integrands_v[i], nout), lb, ub, batch = 10, + nout = nout) + if req.min_dim > 1 || !req.allows_batch || req.nout < nout + continue + end + @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout" + sol = solve(prob, alg, reltol = reltol, abstol = abstol) + @test sol.u≈exact_sol_v[i](dim, nout, lb, ub) rtol=1e-2 + end + end +end + +@testset "Batched Standard Vector Integrands" begin + nout = 2 + for alg in algs + req = alg_req[alg] + for i in 1:length(integrands_v) + for dim in 1:max_dim_test + (lb, ub) = (ones(dim), 3ones(dim)) + prob = IntegralProblem(batch_f_v(integrands_v[i], nout), lb, ub, batch = 10, + nout = nout) + if dim > req.max_dim || dim < req.min_dim || !req.allows_batch || + req.nout < nout + continue + end + @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout" + sol = solve(prob, alg, reltol = reltol, abstol = abstol) + @test sol.u≈exact_sol_v[i](dim, nout, lb, ub) rtol=1e-2 + end + end + end +end + +@testset "In-Place Batched Standard Vector Integrands" begin + nout = 2 + for alg in algs + req = req = alg_req[alg] + for i in 1:length(iip_integrands_v) + for dim in 1:max_dim_test + (lb, ub) = (ones(dim), 3ones(dim)) + prob = IntegralProblem(batch_iip_f_v(integrands_v[i], nout), lb, ub, + batch = 10, nout = nout) + if dim > req.max_dim || dim < req.min_dim || !req.allows_batch || + !req.allows_iip || req.nout < nout + continue + end + @info "Alg = $alg, Integrand = $i, Dimension = $dim, Output Dimension = $nout" + sol = solve(prob, alg, reltol = reltol, abstol = abstol) + @test sol.u≈exact_sol_v[i](dim, nout, lb, ub) rtol=1e-2 + end + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index e634f3a8..77d97b8f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,7 +4,7 @@ using Test function dev_subpkg(subpkg) subpkg_path = joinpath(dirname(@__DIR__), "lib", subpkg) - Pkg.develop(PackageSpec(path=subpkg_path)) + Pkg.develop(PackageSpec(path = subpkg_path)) end dev_subpkg("IntegralsCuba")