diff --git a/Project.toml b/Project.toml
index 34abb67..552f89b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MonteCarloIntegration = "4886b29c-78c9-11e9-0a6e-41e1f4161f7b"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 
@@ -46,8 +47,8 @@ HCubature = "1.5.2"
 LinearAlgebra = "1.10"
 MCIntegration = "0.4.2"
 MonteCarloIntegration = "0.2"
-Pkg = "1.10"
 QuadGK = "2.9"
+Random = "1.10"
 Reexport = "1.0"
 SafeTestsets = "0.1"
 SciMLBase = "2.6"
@@ -67,11 +68,10 @@ FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 MCIntegration = "ea1e2de9-7db7-4b42-91ee-0cd1bf6df167"
-Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "Arblib", "StaticArrays", "FiniteDiff", "Pkg", "SafeTestsets", "Test", "Distributions", "ForwardDiff", "Zygote", "ChainRulesCore", "FastGaussQuadrature", "Cuba", "Cubature", "MCIntegration"]
+test = ["Aqua", "Arblib", "StaticArrays", "FiniteDiff", "SafeTestsets", "Test", "Distributions", "ForwardDiff", "Zygote", "ChainRulesCore", "FastGaussQuadrature", "Cuba", "Cubature", "MCIntegration"]
diff --git a/docs/src/basics/solve.md b/docs/src/basics/solve.md
index d513a1e..c572b00 100644
--- a/docs/src/basics/solve.md
+++ b/docs/src/basics/solve.md
@@ -3,7 +3,3 @@
 ```@docs
 solve(prob::IntegralProblem, alg::SciMLBase.AbstractIntegralAlgorithm)
 ```
-
-Additionally, the extra keyword arguments are splatted to the library calls, so
-see the documentation of the integrator library for all the extra details.
-These extra keyword arguments are not guaranteed to act uniformly.
diff --git a/ext/IntegralsArblibExt.jl b/ext/IntegralsArblibExt.jl
index e0c408a..0eac48a 100644
--- a/ext/IntegralsArblibExt.jl
+++ b/ext/IntegralsArblibExt.jl
@@ -15,21 +15,19 @@ function Integrals.__solvebp_call(
 
     if isinplace(prob)
         res = Acb(0)
-        y_ = similar(prob.f.integrand_prototype, typeof(res))
+        @assert res isa eltype(prob.f.integrand_prototype) "Arblib require inplace prototypes to store Acb elements"
+        y_ = similar(prob.f.integrand_prototype)
         f_ = (y, x; kws...) -> (prob.f(y_, x, p; kws...); Arblib.set!(y, only(y_)))
         val = Arblib.integrate!(f_, res, lb, ub, atol = abstol, rtol = reltol,
             check_analytic = alg.check_analytic, take_prec = alg.take_prec,
             warn_on_no_convergence = alg.warn_on_no_convergence, opts = alg.opts)
-        SciMLBase.build_solution(
-            prob, alg, val, get_radius(val), retcode = ReturnCode.Success)
     else
         f_ = (x; kws...) -> only(prob.f(x, p; kws...))
         val = Arblib.integrate(f_, lb, ub, atol = abstol, rtol = reltol,
             check_analytic = alg.check_analytic, take_prec = alg.take_prec,
             warn_on_no_convergence = alg.warn_on_no_convergence, opts = alg.opts)
-        SciMLBase.build_solution(
-            prob, alg, val, get_radius(val), retcode = ReturnCode.Success)
     end
+    SciMLBase.build_solution(prob, alg, val, get_radius(val), retcode = ReturnCode.Success)
 end
 
 function get_radius(ball)
diff --git a/ext/IntegralsCubaExt.jl b/ext/IntegralsCubaExt.jl
index f5c1697..82bad2a 100644
--- a/ext/IntegralsCubaExt.jl
+++ b/ext/IntegralsCubaExt.jl
@@ -18,8 +18,12 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori
         throw(ArgumentError("Cuba.jl only supports real-valued integrands"))
     # we could support other types by multiplying by the jacobian determinant at the end
 
-    if prob.f isa BatchIntegralFunction
-        nvec = min(maxiters, prob.f.max_batch)
+    f = prob.f
+    prototype = Integrals.get_prototype(prob)
+    if f isa BatchIntegralFunction
+        fsize = size(prototype)[begin:(end - 1)]
+        ncomp = prod(fsize)
+        nvec = min(maxiters, f.max_batch)
         # nvec == 1 in Cuba will change vectors to matrices, so we won't support it when
         # batching
         nvec > 1 ||
@@ -33,24 +37,21 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori
             scale = x -> scale_x!(view(_x, :, 1:size(x, 2)), ub, lb, x)
         end
 
-        if isinplace(prob)
-            fsize = size(prob.f.integrand_prototype)[begin:(end - 1)]
-            y = similar(prob.f.integrand_prototype, fsize..., nvec)
-            ax = map(_ -> (:), fsize)
-            f = function (x, dx)
-                dy = @view(y[ax..., begin:(begin + size(dx, 2) - 1)])
-                prob.f(dy, scale(x), p)
-                dx .= reshape(dy, :, size(dx, 2)) .* vol
+        if isinplace(f)
+            ax = ntuple(_ -> (:), length(fsize))
+            _f = let y_ = similar(prototype, fsize..., nvec)
+                function (u, _y)
+                    y = @view(y_[ax..., begin:(begin + size(_y, 2) - 1)])
+                    f(y, scale(u), p)
+                    _y .= reshape(y, size(_y)) .* vol
+                end
             end
         else
-            y = mid isa Number ? prob.f(typeof(mid)[], p) :
-                prob.f(Matrix{typeof(mid)}(undef, length(mid), 0), p)
-            fsize = size(y)[begin:(end - 1)]
-            f = (x, dx) -> dx .= reshape(prob.f(scale(x), p), :, size(dx, 2)) .* vol
+            _f = (u, y) -> y .= reshape(f(scale(u), p), size(y)) .* vol
         end
-        ncomp = prod(fsize)
     else
         nvec = 1
+        ncomp = length(prototype)
 
         if mid isa Real
             scale = x -> scale_x(ub, lb, only(x))
@@ -59,31 +60,33 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori
             scale = x -> scale_x!(_x, ub, lb, x)
         end
 
-        if isinplace(prob)
-            y = similar(prob.f.integrand_prototype)
-            f = (x, dx) -> dx .= vec(prob.f(y, scale(x), p)) .* vol
+        if isinplace(f)
+            _f = let y = similar(prototype)
+                (u, _y) -> begin
+                    f(y, scale(u), p)
+                    _y .= vec(y) .* vol
+                end
+            end
         else
-            y = prob.f(mid, p)
-            f = (x, dx) -> dx .= Iterators.flatten(prob.f(scale(x), p)) .* vol
+            _f = (u, y) -> y .= Iterators.flatten(f(scale(u), p)) .* vol
         end
-        ncomp = length(y)
     end
 
-    if alg isa CubaVegas
-        out = Cuba.vegas(f, ndim, ncomp; rtol = reltol,
+    out = if alg isa CubaVegas
+        Cuba.vegas(_f, ndim, ncomp; rtol = reltol,
             atol = abstol, nvec = nvec,
             maxevals = maxiters,
             flags = alg.flags, seed = alg.seed, minevals = alg.minevals,
             nstart = alg.nstart, nincrease = alg.nincrease,
             gridno = alg.gridno)
     elseif alg isa CubaSUAVE
-        out = Cuba.suave(f, ndim, ncomp; rtol = reltol,
+        Cuba.suave(_f, ndim, ncomp; rtol = reltol,
             atol = abstol, nvec = nvec,
             maxevals = maxiters,
             flags = alg.flags, seed = alg.seed, minevals = alg.minevals,
             nnew = alg.nnew, nmin = alg.nmin, flatness = alg.flatness)
     elseif alg isa CubaDivonne
-        out = Cuba.divonne(f, ndim, ncomp; rtol = reltol,
+        Cuba.divonne(_f, ndim, ncomp; rtol = reltol,
             atol = abstol, nvec = nvec,
             maxevals = maxiters,
             flags = alg.flags, seed = alg.seed, minevals = alg.minevals,
@@ -91,26 +94,26 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori
             maxpass = alg.maxpass, border = alg.border,
             maxchisq = alg.maxchisq, mindeviation = alg.mindeviation)
     elseif alg isa CubaCuhre
-        out = Cuba.cuhre(f, ndim, ncomp; rtol = reltol,
+        Cuba.cuhre(_f, ndim, ncomp; rtol = reltol,
             atol = abstol, nvec = nvec,
             maxevals = maxiters,
             flags = alg.flags, minevals = alg.minevals, key = alg.key)
     end
 
     # out.integral is a Vector{Float64}, but we want to return it to the shape of the integrand
-    if prob.f isa BatchIntegralFunction
-        if y isa AbstractVector
-            val = out.integral[1]
+    val = if f isa BatchIntegralFunction
+        if prototype isa AbstractVector
+            out.integral[1]
         else
-            val = reshape(out.integral, fsize)
+            reshape(out.integral, fsize)
         end
     else
-        if y isa Real
-            val = out.integral[1]
-        elseif y isa AbstractVector
-            val = out.integral
+        if prototype isa Real
+            out.integral[1]
+        elseif prototype isa AbstractVector
+            out.integral
         else
-            val = reshape(out.integral, size(y))
+            reshape(out.integral, size(prototype))
         end
     end
 
diff --git a/ext/IntegralsCubatureExt.jl b/ext/IntegralsCubatureExt.jl
index 9033291..57327cc 100644
--- a/ext/IntegralsCubatureExt.jl
+++ b/ext/IntegralsCubatureExt.jl
@@ -13,143 +13,141 @@ function Integrals.__solvebp_call(prob::IntegralProblem,
     mid = (lb + ub) / 2
 
     # we get to pick fdim or not based on the IntegralFunction and its output dimensions
-    y = if prob.f isa BatchIntegralFunction
-        isinplace(prob.f) ? prob.f.integrand_prototype :
-        mid isa Number ? prob.f(eltype(mid)[], p) :
-        prob.f(Matrix{eltype(mid)}(undef, length(mid), 0), p)
-    else
-        # we evaluate the oop function to decide whether the output should be vectorized
-        isinplace(prob.f) ? prob.f.integrand_prototype : prob.f(mid, p)
-    end
+    f = prob.f
+    prototype = Integrals.get_prototype(prob)
 
-    @assert eltype(y)<:Real "Cubature.jl is only compatible with real-valued integrands"
+    @assert eltype(prototype)<:Real "Cubature.jl is only compatible with real-valued integrands"
 
-    if prob.f isa BatchIntegralFunction
-        if y isa AbstractVector # this branch could be omitted since the following one should work similarly
-            if isinplace(prob)
+    if f isa BatchIntegralFunction
+        if prototype isa AbstractVector # this branch could be omitted since the following one should work similarly
+            if isinplace(f)
                 # dx is a Vector, but we provide the integrand a vector of the same type as
                 # y, which needs to be resized since the number of batch points changes.
-                dy = similar(y)
-                f = (x, dx) -> begin
-                    resize!(dy, length(dx))
-                    prob.f(dy, x, p)
-                    dx .= dy
+                _f = let y = similar(prototype)
+                    (u, v) -> begin
+                        resize!(y, length(v))
+                        f(y, u, p)
+                        v .= y
+                    end
                 end
             else
-                f = (x, dx) -> (dx .= prob.f(x, p))
+                _f = (u, v) -> (v .= f(u, p))
             end
             if mid isa Number
                 if alg isa CubatureJLh
-                    val, err = Cubature.hquadrature_v(f, lb, ub;
+                    val, err = Cubature.hquadrature_v(_f, lb, ub;
                         reltol = reltol, abstol = abstol,
                         maxevals = maxiters)
                 else
-                    val, err = Cubature.pquadrature_v(f, lb, ub;
+                    val, err = Cubature.pquadrature_v(_f, lb, ub;
                         reltol = reltol, abstol = abstol,
                         maxevals = maxiters)
                 end
             else
                 if alg isa CubatureJLh
-                    val, err = Cubature.hcubature_v(f, lb, ub;
+                    val, err = Cubature.hcubature_v(_f, lb, ub;
                         reltol = reltol, abstol = abstol,
                         maxevals = maxiters)
                 else
-                    val, err = Cubature.pcubature_v(f, lb, ub;
+                    val, err = Cubature.pcubature_v(_f, lb, ub;
                         reltol = reltol, abstol = abstol,
                         maxevals = maxiters)
                 end
             end
-        elseif y isa AbstractArray
-            bfsize = size(y)[begin:(end - 1)]
-            bfdim = prod(bfsize)
-            if isinplace(prob)
+        elseif prototype isa AbstractArray
+            fsize = size(prototype)[begin:(end - 1)]
+            fdim = prod(fsize)
+            if isinplace(f)
                 # dx is a Matrix, but to provide a buffer of the same type as y, we make
                 # would like to make views of a larger buffer, but CubatureJL doesn't set
                 # a hard limit for max_batch, so we allocate a new buffer with the needed size
-                f = (x, dx) -> begin
-                    dy = similar(y, bfsize..., size(dx, 2))
-                    prob.f(dy, x, p)
-                    dx .= reshape(dy, bfdim, size(dx, 2))
+                _f = let fsize = fsize
+                    (u, v) -> begin
+                        y = similar(prototype, fsize..., size(v, 2))
+                        f(y, u, p)
+                        v .= reshape(y, fdim, size(v, 2))
+                    end
                 end
             else
-                f = (x, dx) -> (dx .= reshape(prob.f(x, p), bfdim, size(dx, 2)))
+                _f = (u, v) -> (v .= reshape(f(u, p), fdim, size(v, 2)))
             end
             if mid isa Number
                 if alg isa CubatureJLh
-                    val_, err = Cubature.hquadrature_v(bfdim, f, lb, ub;
+                    val_, err = Cubature.hquadrature_v(fdim, _f, lb, ub;
                         reltol = reltol, abstol = abstol,
                         maxevals = maxiters, error_norm = alg.error_norm)
                 else
-                    val_, err = Cubature.pquadrature_v(bfdim, f, lb, ub;
+                    val_, err = Cubature.pquadrature_v(fdim, _f, lb, ub;
                         reltol = reltol, abstol = abstol,
                         maxevals = maxiters, error_norm = alg.error_norm)
                 end
             else
                 if alg isa CubatureJLh
-                    val_, err = Cubature.hcubature_v(bfdim, f, lb, ub;
+                    val_, err = Cubature.hcubature_v(fdim, _f, lb, ub;
                         reltol = reltol, abstol = abstol,
                         maxevals = maxiters, error_norm = alg.error_norm)
                 else
-                    val_, err = Cubature.pcubature_v(bfdim, f, lb, ub;
+                    val_, err = Cubature.pcubature_v(fdim, _f, lb, ub;
                         reltol = reltol, abstol = abstol,
                         maxevals = maxiters, error_norm = alg.error_norm)
                 end
             end
-            val = reshape(val_, bfsize...)
+            val = reshape(val_, fsize...)
         else
             error("BatchIntegralFunction integrands must be arrays for Cubature.jl")
         end
     else
-        if y isa Real
+        if prototype isa Real
             # no inplace in this case, since the integrand_prototype would be mutable
-            f = x -> prob.f(x, p)
+            _f = u -> f(u, p)
             if lb isa Number
                 if alg isa CubatureJLh
-                    val, err = Cubature.hquadrature(f, lb, ub;
+                    val, err = Cubature.hquadrature(_f, lb, ub;
                         reltol = reltol, abstol = abstol,
                         maxevals = maxiters)
                 else
-                    val, err = Cubature.pquadrature(f, lb, ub;
+                    val, err = Cubature.pquadrature(_f, lb, ub;
                         reltol = reltol, abstol = abstol,
                         maxevals = maxiters)
                 end
             else
                 if alg isa CubatureJLh
-                    val, err = Cubature.hcubature(f, lb, ub;
+                    val, err = Cubature.hcubature(_f, lb, ub;
                         reltol = reltol, abstol = abstol,
                         maxevals = maxiters)
                 else
-                    val, err = Cubature.pcubature(f, lb, ub;
+                    val, err = Cubature.pcubature(_f, lb, ub;
                         reltol = reltol, abstol = abstol,
                         maxevals = maxiters)
                 end
             end
-        elseif y isa AbstractArray
-            fsize = size(y)
-            fdim = length(y)
+        elseif prototype isa AbstractArray
+            fsize = size(prototype)
+            fdim = length(prototype)
             if isinplace(prob)
-                dy = similar(y)
-                f = (x, v) -> (prob.f(dy, x, p); v .= vec(dy))
+                _f = let y = similar(prototype)
+                    (u, v) -> (f(y, u, p); v .= vec(y))
+                end
             else
-                f = (x, v) -> (v .= vec(prob.f(x, p)))
+                _f = (u, v) -> (v .= vec(f(u, p)))
             end
             if mid isa Number
                 if alg isa CubatureJLh
-                    val_, err = Cubature.hquadrature(fdim, f, lb, ub;
+                    val_, err = Cubature.hquadrature(fdim, _f, lb, ub;
                         reltol = reltol, abstol = abstol,
                         maxevals = maxiters, error_norm = alg.error_norm)
                 else
-                    val_, err = Cubature.pquadrature(fdim, f, lb, ub;
+                    val_, err = Cubature.pquadrature(fdim, _f, lb, ub;
                         reltol = reltol, abstol = abstol,
                         maxevals = maxiters, error_norm = alg.error_norm)
                 end
             else
                 if alg isa CubatureJLh
-                    val_, err = Cubature.hcubature(fdim, f, lb, ub;
+                    val_, err = Cubature.hcubature(fdim, _f, lb, ub;
                         reltol = reltol, abstol = abstol,
                         maxevals = maxiters, error_norm = alg.error_norm)
                 else
-                    val_, err = Cubature.pcubature(fdim, f, lb, ub;
+                    val_, err = Cubature.pcubature(fdim, _f, lb, ub;
                         reltol = reltol, abstol = abstol,
                         maxevals = maxiters, error_norm = alg.error_norm)
                 end
diff --git a/ext/IntegralsFastGaussQuadratureExt.jl b/ext/IntegralsFastGaussQuadratureExt.jl
index 1dc5c39..704145d 100644
--- a/ext/IntegralsFastGaussQuadratureExt.jl
+++ b/ext/IntegralsFastGaussQuadratureExt.jl
@@ -15,7 +15,13 @@ Integrals.gausslegendre(n) = FastGaussQuadrature.gausslegendre(n)
 function gauss_legendre(f::F, p, lb, ub, nodes, weights) where {F}
     scale = (ub - lb) / 2
     shift = (lb + ub) / 2
-    I = mapreduce((w, x) -> w * f(scale * x + shift, p), +, weights, nodes)
+    # TODO reuse Integrals.evalrule instead
+    x0, xs = Iterators.peel(nodes)
+    w0, ws = Iterators.peel(weights)
+    I = w0 * f(scale * x0 + shift, p)
+    for (w, x) in zip(ws, xs)
+        I += w * f(scale * x + shift, p)
+    end
     return scale * I
 end
 function composite_gauss_legendre(f::F, p, lb, ub, nodes, weights, subintervals) where {F}
diff --git a/ext/IntegralsForwardDiffExt.jl b/ext/IntegralsForwardDiffExt.jl
index ee840cd..f09309c 100644
--- a/ext/IntegralsForwardDiffExt.jl
+++ b/ext/IntegralsForwardDiffExt.jl
@@ -3,6 +3,13 @@ using Integrals
 isdefined(Base, :get_extension) ? (using ForwardDiff) : (using ..ForwardDiff)
 ### Forward-Mode AD Intercepts
 
+function Integrals._evaluate!(f, y, u,
+        p::Union{D, AbstractArray{<:D}}) where {T, V, P, D <: ForwardDiff.Dual{T, V, P}}
+    dy = similar(y, replace_dualvaltype(eltype(p), eltype(y)))
+    f(dy, u, p)
+    return dy
+end
+
 # Default to direct AD on solvers
 function Integrals.__solvebp(cache, alg, sensealg, domain,
         p::Union{D, AbstractArray{<:D}};
@@ -19,9 +26,16 @@ function Integrals.__solvebp(cache, alg, sensealg, domain,
             IntegralFunction{true}(cache.f.f, dprototype)
         end
         prob = Integrals.build_problem(cache)
-        dprob = remake(prob, f = df)
-        dcache = init(
-            dprob, alg; sensealg = sensealg, do_inf_transformation = Val(false), kwargs...)
+        dcache = Integrals.IntegralCache(cache.iip,
+            df,
+            domain,
+            p,
+            cache.prob_kwargs,
+            alg,
+            sensealg,
+            cache.kwargs,
+            cache.cacheval
+        )
         Integrals.__solvebp_call(dcache, alg, sensealg, domain, p; kwargs...)
     else
         Integrals.__solvebp_call(cache, alg, sensealg, domain, p; kwargs...)
@@ -90,13 +104,11 @@ function Integrals.__solvebp(
 
     prob = Integrals.build_problem(cache)
     dp_prob = remake(prob, f = dfdp, p = rawp)
-    # the infinity transformation was already applied to f so we don't apply it to dfdp
     dp_cache = init(dp_prob,
         alg;
         sensealg = sensealg,
-        do_inf_transformation = Val(false),
         cache.kwargs...)
-    dual = Integrals.__solvebp_call(dp_cache, alg, sensealg, domain, rawp; kwargs...)
+    dual = solve!(dp_cache)
 
     res = reinterpret(reshape, DT, dual.u)
     # unwrap the dual when the primal would return a scalar
diff --git a/ext/IntegralsMCIntegrationExt.jl b/ext/IntegralsMCIntegrationExt.jl
index e4d6779..4255bac 100644
--- a/ext/IntegralsMCIntegrationExt.jl
+++ b/ext/IntegralsMCIntegrationExt.jl
@@ -2,46 +2,54 @@ module IntegralsMCIntegrationExt
 
 using MCIntegration, Integrals
 
+_oftype(::Number, x) = only(x)
+_oftype(y, x) = oftype(y, x)
+
 function Integrals.__solvebp_call(prob::IntegralProblem, alg::VEGASMC, sensealg, domain, p;
         reltol = nothing, abstol = nothing, maxiters = 10)
     lb, ub = domain
-    mid = vec(collect((lb + ub) / 2))
+    mid = (lb + ub) / 2
+    tmp = vec(collect(mid))
     var = Continuous(vec([tuple(a, b) for (a, b) in zip(lb, ub)]))
 
-    if prob.f isa BatchIntegralFunction
+    f = prob.f
+    if f isa BatchIntegralFunction
         error("VEGASMC doesn't support batching. See https://github.com/numericalEFT/MCIntegration.jl/issues/29")
     else
+        prototype = Integrals.get_prototype(prob)
         if isinplace(prob)
-            f0 = similar(prob.f.integrand_prototype)
-            f_ = (x, f, c) -> begin
-                n = 0
-                for v in x
-                    mid[n += 1] = first(v)
+            _f = let y = similar(prototype)
+                (u, _y, c) -> begin
+                    n = 0
+                    for v in u
+                        tmp[n += 1] = first(v)
+                    end
+                    f(y, _oftype(mid, tmp), p)
+                    _y .= vec(y)
                 end
-                prob.f(f0, mid, p)
-                f .= vec(f0)
             end
         else
-            f0 = prob.f(mid, p)
-            f_ = (x, c) -> begin
+            _f = (u, c) -> begin
                 n = 0
-                for v in x
-                    mid[n += 1] = first(v)
+                for v in u
+                    tmp[n += 1] = first(v)
                 end
-                fx = prob.f(mid, p)
-                fx isa AbstractArray ? vec(fx) : fx
+                y = f(_oftype(mid, tmp), p)
+                y isa AbstractArray ? vec(y) : y
             end
         end
-        dof = ones(Int, length(f0)) # each composite Continuous var gets 1 dof
-        res = integrate(f_; var, dof, inplace = isinplace(prob), type = eltype(f0),
+        dof = ones(Int, length(prototype)) # each composite Continuous var gets 1 dof
+        res = integrate(_f; var, dof, inplace = isinplace(prob), type = eltype(prototype),
             solver = :vegasmc, niter = maxiters, verbose = -2, print = -2, alg.kws...)
-        out, err, chi = if f0 isa Number
+        # the package itself is not type-stable
+        out::typeof(prototype), err::typeof(prototype), chi2 = if prototype isa Number
             map(only, (res.mean, res.stdev, res.chi2))
         else
-            map(a -> reshape(a, size(f0)), (res.mean, res.stdev, res.chi2))
+            map(a -> reshape(a, size(prototype)), (res.mean, res.stdev, res.chi2))
         end
-        SciMLBase.build_solution(
-            prob, alg, out, err, chi = chi, retcode = ReturnCode.Success)
+        chi::typeof(prototype) = map(sqrt, chi2)
+        SciMLBase.build_solution(prob, alg, out, err, chi = chi,
+            retcode = ReturnCode.Success)
     end
 end
 
diff --git a/ext/IntegralsZygoteExt.jl b/ext/IntegralsZygoteExt.jl
index 1b38956..bb5d457 100644
--- a/ext/IntegralsZygoteExt.jl
+++ b/ext/IntegralsZygoteExt.jl
@@ -14,6 +14,43 @@ 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.__solve), cache::Integrals.IntegralCache,
+        alg::Integrals.ChangeOfVariables, sensealg, udomain, p;
+        kwargs...)
+    _cache, vdomain = Integrals._change_variables(cache, alg, sensealg, udomain, p)
+    sol, back = Zygote.pullback((args...) -> Integrals.__solve(args...; kwargs...),
+        _cache, alg.alg, sensealg, vdomain, p)
+    function change_of_variables_pullback(Δ)
+        return (NoTangent(), back(Δ)...)
+    end
+    prob = Integrals.build_problem(cache)
+    _sol = SciMLBase.build_solution(
+        prob, alg.alg, sol.u, sol.resid, chi = sol.chi, retcode = sol.retcode, stats = sol.stats)
+    return _sol, change_of_variables_pullback
+end
+
+# we will need to implement the following adjoints when we compute ∂f/∂u
+function ChainRulesCore.rrule(::typeof(Integrals.substitute_v), args...)
+    function substitute_v_pullback(_)
+        return NoTangent(), ntuple(_ -> NoTangent(), length(args))...
+    end
+    return Integrals.substitute_v(args...), substitute_v_pullback
+end
+function ChainRulesCore.rrule(::typeof(Integrals.substitute_bv), args...)
+    function substitute_bv_pullback(_)
+        return NoTangent(), ntuple(_ -> NoTangent(), length(args))...
+    end
+    return Integrals.substitute_bv(args...), substitute_bv_pullback
+end
+function ChainRulesCore.rrule(::typeof(Integrals._evaluate!), f, y, u, p)
+    out, back = Zygote.pullback(y, u, p) do y, u, p
+        b = Zygote.Buffer(y)
+        f(b, u, p)
+        return copy(b)
+    end
+    out, Δ -> (NoTangent(), NoTangent(), back(Δ)...)
+end
+
 function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, sensealg, domain,
         p;
         kwargs...)
@@ -89,16 +126,10 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal
         dp_cache = init(dp_prob,
             alg;
             sensealg = sensealg,
-            do_inf_transformation = Val(false),
             cache.kwargs...)
 
         project_p = ProjectTo(p)
-        dp = project_p(Integrals.__solvebp_call(dp_cache,
-            alg,
-            sensealg,
-            domain,
-            p;
-            kwargs...).u)
+        dp = project_p(solve!(dp_cache).u)
 
         lb, ub = domain
         if lb isa Number
diff --git a/src/Integrals.jl b/src/Integrals.jl
index 91d4772..5edee6d 100644
--- a/src/Integrals.jl
+++ b/src/Integrals.jl
@@ -7,7 +7,9 @@ end
 using Reexport, MonteCarloIntegration, QuadGK, HCubature
 @reexport using SciMLBase
 using LinearAlgebra
+using Random
 
+include("algorithms_meta.jl")
 include("common.jl")
 include("algorithms.jl")
 include("algorithms_sampled.jl")
@@ -63,158 +65,215 @@ function checkkwargs(kwargs...)
     return nothing
 end
 
+# Give a layer for meta algorithms above AD
+__solve(args...; kwargs...) = __solvebp(args...; kwargs...)
 # Give a layer to intercept with AD
 __solvebp(args...; kwargs...) = __solvebp_call(args...; kwargs...)
 
-function quadgk_prob_types(f, lb::T, ub::T, p, nrm) where {T}
-    DT = float(T)   # we need to be careful to infer the same result as `evalrule`
-    RT = Base.promote_op(*, DT, Base.promote_op(f, DT, typeof(p)))    # kernel
-    NT = Base.promote_op(nrm, RT)
-    return DT, RT, NT
+function init_cacheval(alg::ChangeOfVariables, prob::IntegralProblem)
+    f, domain = alg.fu2gv(prob.f, prob.domain)
+    cache_alg = init_cacheval(alg.alg, remake(prob, f = f, domain = domain))
+    return (alg = cache_alg,)
 end
-function init_cacheval(alg::QuadGKJL, prob::IntegralProblem)
-    lb, ub = map(first, prob.domain)
-    DT, RT, NT = quadgk_prob_types(prob.f, lb, ub, prob.p, alg.norm)
-    return (isconcretetype(RT) ? QuadGK.alloc_segbuf(DT, RT, NT) : nothing)
+
+function __solve(cache::IntegralCache, alg::ChangeOfVariables, sensealg, udomain, p;
+        kwargs...)
+    _cache, vdomain = _change_variables(cache, alg, sensealg, udomain, p)
+    sol = __solve(_cache, alg.alg, sensealg, vdomain, p; kwargs...)
+    prob = build_problem(cache)
+    return SciMLBase.build_solution(
+        prob, alg.alg, sol.u, sol.resid, chi = sol.chi, retcode = sol.retcode, stats = sol.stats)
 end
-function refresh_cacheval(cacheval, alg::QuadGKJL, prob)
+
+function _change_variables(cache, alg, sensealg, udomain, p)
+    cacheval = cache.cacheval.alg
+    g, vdomain = alg.fu2gv(cache.f, udomain)
+    _cache = IntegralCache(Val(isinplace(g)),
+        g,
+        vdomain,
+        p,
+        cache.prob_kwargs,
+        alg.alg,
+        sensealg,
+        cache.kwargs,
+        cacheval)
+    return _cache, vdomain
+end
+
+function get_prototype(prob::IntegralProblem)
+    f = prob.f
+    f.integrand_prototype !== nothing && return f.integrand_prototype
+    isinplace(f) && throw(ArgumentError("in-place integral functions require a prototype"))
+    lb, ub = prob.domain
+    mid = (lb + ub) / 2
+    p = prob.p
+    if f isa BatchIntegralFunction
+        mid isa Number ? f(eltype(mid)[], p) :
+        f(Matrix{eltype(mid)}(undef, length(mid), 0), p)
+    else
+        f(mid, p)
+    end
+end
+
+function init_cacheval(alg::QuadGKJL, prob::IntegralProblem)
+    alg.buffer === nothing && return
     lb, ub = map(first, prob.domain)
-    DT, RT, NT = quadgk_prob_types(prob.f, lb, ub, prob.p, alg.norm)
-    isconcretetype(RT) || return nothing
-    T = QuadGK.Segment{DT, RT, NT}
-    return (cacheval isa Vector{T} ? cacheval : QuadGK.alloc_segbuf(DT, RT, NT))
+    mid = (lb + ub) / 2
+    prototype = get_prototype(prob) * mid
+    TX = typeof(mid)
+    TI = typeof(prototype)
+    TE = typeof(alg.norm(prototype))
+    QuadGK.alloc_segbuf(TX, TI, TE)
 end
 
 function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, domain, p;
         reltol = 1e-8, abstol = 1e-8,
         maxiters = typemax(Int))
     prob = build_problem(cache)
-    lb_, ub_ = domain
     lb, ub = map(first, domain)
-    if !isone(length(lb_)) || !isone(length(ub_))
+    if any(!isone ∘ length, domain)
         error("QuadGKJL only accepts one-dimensional quadrature problems.")
     end
 
-    mid = ((lb + ub) / 2)
-    if prob.f isa BatchIntegralFunction
-        if isinplace(prob)
+    f = cache.f
+    mid = (lb + ub) / 2
+    if f isa BatchIntegralFunction
+        if isinplace(f)
             # quadgk only works with vector buffers. If the buffer is an array, we have to
             # turn it into a vector of arrays
-            bu = prob.f.integrand_prototype
-            f = if bu isa AbstractVector
-                BatchIntegrand((y, x) -> prob.f(y, x, p), similar(bu))
+            prototype = f.integrand_prototype
+            _f = if prototype isa AbstractVector
+                BatchIntegrand((y, u) -> f(y, u, p), similar(prototype))
             else
-                fsize = size(bu)[begin:(end - 1)]
-                BatchIntegrand{Array{eltype(bu), ndims(bu) - 1}}() do y, x
-                    y_ = similar(bu, fsize..., length(y))
-                    prob.f(y_, x, p)
-                    map!(collect, y, eachslice(y_; dims = ndims(bu)))
-                    return nothing
+                fsize = size(prototype)[begin:(end - 1)]
+                BatchIntegrand{Array{eltype(prototype), length(fsize)}}() do v, u
+                    let y = similar(v, eltype(eltype(v)), fsize..., length(v))
+                        f(y, u, p)
+                        map!(collect, v, eachslice(y; dims = ndims(y)))
+                    end
+                    return
                 end
             end
-            val, err = quadgk(f, lb, ub, segbuf = cache.cacheval, maxevals = maxiters,
+            val, err = quadgk(_f, lb, ub, segbuf = cache.cacheval, maxevals = maxiters,
                 rtol = reltol, atol = abstol, order = alg.order, norm = alg.norm)
-            SciMLBase.build_solution(prob, alg, val, err, retcode = ReturnCode.Success)
         else
-            u = prob.f(typeof(mid)[], p)
-            f = if u isa AbstractVector
-                BatchIntegrand((y, x) -> y .= prob.f(x, p), u)
+            prototype = f(typeof(mid)[], p)
+            _f = if prototype isa AbstractVector
+                BatchIntegrand((y, u) -> y .= f(u, p), prototype)
             else
-                BatchIntegrand{Array{eltype(u), ndims(u) - 1}}() do y, x
-                    map!(collect, y, eachslice(prob.f(x, p); dims = ndims(u)))
-                    return nothing
+                BatchIntegrand{Array{eltype(prototype), ndims(prototype) - 1}}() do v, u
+                    y = f(u, p)
+                    map!(collect, v, eachslice(y; dims = ndims(y)))
+                    return
                 end
             end
-            val, err = quadgk(f, lb, ub, segbuf = cache.cacheval, maxevals = maxiters,
+            val, err = quadgk(_f, lb, ub, segbuf = cache.cacheval, maxevals = maxiters,
                 rtol = reltol, atol = abstol, order = alg.order, norm = alg.norm)
-            SciMLBase.build_solution(prob, alg, val, err, retcode = ReturnCode.Success)
         end
     else
-        if isinplace(prob)
-            result = prob.f.integrand_prototype * mid   # result may have different units than prototype
-            f = (y, x) -> prob.f(y, x, p)
-            val, err = quadgk!(
-                f, result, lb, ub, segbuf = cache.cacheval, maxevals = maxiters,
+        if isinplace(f)
+            result = f.integrand_prototype * mid   # result may have different units than prototype
+            _f = (y, u) -> f(y, u, p)
+            val, err = quadgk!(_f, result, lb, ub, segbuf = cache.cacheval,
+                maxevals = maxiters,
                 rtol = reltol, atol = abstol, order = alg.order, norm = alg.norm)
-            SciMLBase.build_solution(prob, alg, val, err, retcode = ReturnCode.Success)
         else
-            f = x -> prob.f(x, p)
-            val, err = quadgk(f, lb, ub, segbuf = cache.cacheval, maxevals = maxiters,
+            _f = u -> f(u, p)
+            val, err = quadgk(_f, lb, ub, segbuf = cache.cacheval, maxevals = maxiters,
                 rtol = reltol, atol = abstol, order = alg.order, norm = alg.norm)
-            SciMLBase.build_solution(prob, alg, val, err, retcode = ReturnCode.Success)
         end
     end
+    SciMLBase.build_solution(prob, alg, val, err, retcode = ReturnCode.Success)
 end
 
-function __solvebp_call(prob::IntegralProblem, alg::HCubatureJL, sensealg, domain, p;
+function init_cacheval(alg::HCubatureJL, prob::IntegralProblem)
+    alg.buffer === nothing && return
+    prototype = get_prototype(prob)
+    lb, ub = map(x -> x isa Number ? tuple(x) : x, prob.domain)
+    HCubature.hcubature_buffer(x -> prototype, lb, ub; norm = alg.norm)
+end
+function __solvebp_call(cache::IntegralCache, alg::HCubatureJL, sensealg, domain, p;
         reltol = 1e-8, abstol = 1e-8,
         maxiters = typemax(Int))
+    prob = build_problem(cache)
     lb, ub = domain
+    f = cache.f
 
-    @assert prob.f isa IntegralFunction
-    if isinplace(prob)
+    @assert f isa IntegralFunction
+    if isinplace(f)
         # allocate a new output array at each evaluation since HCubature.jl doesn't support
         # inplace ops
-        f = x -> (dx = similar(prob.f.integrand_prototype); prob.f(dx, x, prob.p); dx)
+        prototype = f.integrand_prototype
+        _f = let f = f.f
+            u -> (y = similar(prototype); f(y, u, p); y)
+        end
     else
-        f = x -> prob.f(x, prob.p)
+        _f = let f = f.f
+            u -> f(u, p)
+        end
     end
 
-    if lb isa Number
-        val, err = hquadrature(f, lb, ub;
-            rtol = reltol, atol = abstol,
+    val, err = if lb isa Number
+        hquadrature(_f, lb, ub;
+            rtol = reltol, atol = abstol, buffer = cache.cacheval,
             maxevals = maxiters, norm = alg.norm, initdiv = alg.initdiv)
     else
-        val, err = hcubature(f, lb, ub;
-            rtol = reltol, atol = abstol,
-            maxevals = maxiters, norm = alg.norm, initdiv = alg.initdiv)
+        ret = get_prototype(prob) * (prod(ub - lb) / 2) # this calculation for type stability with vector endpoints
+        hcubature(_f, lb, ub;
+            rtol = reltol, atol = abstol, buffer = cache.cacheval,
+            maxevals = maxiters, norm = alg.norm,
+            initdiv = alg.initdiv)::Tuple{typeof(ret), typeof(alg.norm(ret))}
     end
     SciMLBase.build_solution(prob, alg, val, err, retcode = ReturnCode.Success)
 end
 
 function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, domain, p;
-        reltol = 1e-8, abstol = 1e-8,
+        reltol = 1e-4, abstol = 1e-4,
         maxiters = 1000)
     lb, ub = domain
     mid = (lb + ub) / 2
-    if prob.f isa BatchIntegralFunction
+    f = prob.f
+    alg.seed !== nothing && Random.seed!(alg.seed)
+    if f isa BatchIntegralFunction
+        # MonteCarloIntegration v0.0.x passes points as rows of a matrix
+        # MonteCarloIntegration v0.1 passes batches as a vector of views of
+        # a matrix with points as columns of a matrix
+        # see https://github.com/ranjanan/MonteCarloIntegration.jl/issues/16
+        # This is an ugly hack that is compatible with both
+        wrangle = x -> begin
+            xx = eltype(x) <: SubArray ? parent(first(x)) : x'
+            mid isa Number ? vec(xx) : xx
+        end
         if isinplace(prob)
             y = similar(prob.f.integrand_prototype,
                 size(prob.f.integrand_prototype)[begin:(end - 1)]...,
                 prob.f.max_batch)
-            # MonteCarloIntegration v0.0.x passes points as rows of a matrix
-            # MonteCarloIntegration v0.1 passes batches as a vector of views of
-            # a matrix with points as columns of a matrix
-            # see https://github.com/ranjanan/MonteCarloIntegration.jl/issues/16
-            # This is an ugly hack that is compatible with both
-            f = x -> (prob.f(y, eltype(x) <: SubArray ? parent(first(x)) : x', p); vec(y))
+            _f = x -> (f(y, wrangle(x), p); vec(y))
         else
-            y = prob.f(
-                mid isa Number ? typeof(mid)[] :
-                Matrix{eltype(mid)}(undef, length(mid), 0),
-                p)
-            f = x -> prob.f(eltype(x) <: SubArray ? parent(first(x)) : x', p)
+            y = mid isa Number ? f(typeof(mid)[], p) :
+                f(Matrix{eltype(mid)}(undef, length(mid), 0), p)
+            _f = x -> vec(f(wrangle(x), p))
         end
     else
         if isinplace(prob)
             @assert prob.f.integrand_prototype isa
                     AbstractArray{<:Real}&&length(prob.f.integrand_prototype) == 1 "VEGAS only supports Float64-valued integrands"
             y = similar(prob.f.integrand_prototype)
-            f = x -> (prob.f(y, x, p); only(y))
+            _f = x -> (prob.f(y, mid isa Number ? only(x) : x, p); only(y))
         else
             y = prob.f(mid, p)
-            f = x -> only(prob.f(x, prob.p))
+            _f = x -> only(prob.f(mid isa Number ? only(x) : x, prob.p))
         end
     end
 
-    if prob.f isa BatchIntegralFunction
+    if f isa BatchIntegralFunction
         @assert prod(size(y)[begin:(end - 1)]) == 1&&eltype(y) <: Real "VEGAS only supports Float64-valued scalar integrands"
     else
         @assert length(y) == 1&&eltype(y) <: Real "VEGAS only supports Float64-valued scalar integrands"
     end
 
-    ncalls = prob.f isa BatchIntegralFunction ? prob.f.max_batch : alg.ncalls
-    out = vegas(f, lb, ub, rtol = reltol, atol = abstol,
+    ncalls = alg.ncalls
+    out = vegas(_f, lb, ub, rtol = reltol, atol = abstol,
         maxiter = maxiters, nbins = alg.nbins, debug = alg.debug,
         ncalls = ncalls, batch = prob.f isa BatchIntegralFunction)
     val, err, chi = out isa Tuple ? out :
diff --git a/src/algorithms.jl b/src/algorithms.jl
index 4222b99..4841f0a 100644
--- a/src/algorithms.jl
+++ b/src/algorithms.jl
@@ -1,10 +1,12 @@
 """
-    QuadGKJL(; order = 7, norm=norm)
+    QuadGKJL(; order = 7, norm = norm, buffer = nothing)
 
 One-dimensional Gauss-Kronrod integration from QuadGK.jl.
 This method also takes the optional arguments `order` and `norm`.
 Which are the order of the integration rule
-and the norm for calculating the error, respectively
+and the norm for calculating the error, respectively.
+Lastly, the `buffer` keyword, if set, will allocate a buffer to reuse
+for multiple integrals, which may require evaluating the integrand.
 
 ## References
 
@@ -18,20 +20,23 @@ pages={1133--1145},
 year={1997}
 }
 """
-struct QuadGKJL{F} <: SciMLBase.AbstractIntegralAlgorithm where {F}
+struct QuadGKJL{F, B} <: SciMLBase.AbstractIntegralAlgorithm
     order::Int
     norm::F
+    buffer::B
 end
-QuadGKJL(; order = 7, norm = norm) = QuadGKJL(order, norm)
+QuadGKJL(; order = 7, norm = norm, buffer = nothing) = QuadGKJL(order, norm, buffer)
 
 """
-    HCubatureJL(; norm=norm, initdiv=1)
+    HCubatureJL(; norm=norm, initdiv=1, buffer = nothing)
 
 Multidimensional "h-adaptive" integration from HCubature.jl.
 This method also takes the optional arguments `initdiv` and `norm`.
 Which are the initial number of segments
 each dimension of the integration domain is divided into,
 and the norm for calculating the error, respectively.
+Lastly, the `buffer` keyword, if set, will allocate a buffer to reuse
+for multiple integrals, which may require evaluating the integrand.
 
 ## References
 
@@ -46,14 +51,17 @@ year={1980},
 publisher={Elsevier}
 }
 """
-struct HCubatureJL{F} <: SciMLBase.AbstractIntegralAlgorithm where {F}
+struct HCubatureJL{F, B} <: SciMLBase.AbstractIntegralAlgorithm
     initdiv::Int
     norm::F
+    buffer::B
+end
+function HCubatureJL(; initdiv = 1, norm = norm, buffer = nothing)
+    HCubatureJL(initdiv, norm, buffer)
 end
-HCubatureJL(; initdiv = 1, norm = norm) = HCubatureJL(initdiv, norm)
 
 """
-    VEGAS(; nbins = 100, ncalls = 1000, debug=false)
+    VEGAS(; nbins = 100, ncalls = 1000, debug=false, seed = nothing)
 
 Multidimensional adaptive Monte Carlo integration from MonteCarloIntegration.jl.
 Importance sampling is used to reduce variance.
@@ -80,12 +88,15 @@ year={1978},
 publisher={Elsevier}
 }
 """
-struct VEGAS <: SciMLBase.AbstractIntegralAlgorithm
+struct VEGAS{S} <: SciMLBase.AbstractIntegralAlgorithm
     nbins::Int
     ncalls::Int
     debug::Bool
+    seed::S
+end
+function VEGAS(; nbins = 100, ncalls = 1000, debug = false, seed = nothing)
+    VEGAS(nbins, ncalls, debug, seed)
 end
-VEGAS(; nbins = 100, ncalls = 1000, debug = false) = VEGAS(nbins, ncalls, debug)
 
 """
     GaussLegendre{C, N, W}
diff --git a/src/algorithms_meta.jl b/src/algorithms_meta.jl
new file mode 100644
index 0000000..8e9fc30
--- /dev/null
+++ b/src/algorithms_meta.jl
@@ -0,0 +1,15 @@
+abstract type AbstractIntegralMetaAlgorithm <: SciMLBase.AbstractIntegralAlgorithm end
+
+"""
+    ChangeOfVariables(fu2gv, alg)
+
+Apply a change of variables from `∫ f(u,p) du` to an equivalent integral `∫ g(v,p) dv` using
+a helper function `fu2gv(f, u_domain) -> (g, v_domain)` where `f` and `g` should be
+integral functions. Acts as a wrapper to algorithm `alg`
+"""
+# internal algorithm
+struct ChangeOfVariables{T, A <: SciMLBase.AbstractIntegralAlgorithm} <:
+       AbstractIntegralMetaAlgorithm
+    fu2gv::T
+    alg::A
+end
diff --git a/src/common.jl b/src/common.jl
index cd31327..4c82fea 100644
--- a/src/common.jl
+++ b/src/common.jl
@@ -19,15 +19,17 @@ function SciMLBase.init(prob::IntegralProblem{iip},
         sensealg = ReCallVJP(ZygoteVJP()),
         do_inf_transformation = nothing, kwargs...) where {iip}
     checkkwargs(kwargs...)
-    prob = transformation_if_inf(prob, do_inf_transformation)
-    cacheval = init_cacheval(alg, prob)
+    do_inf_transformation === nothing ||
+        @warn "do_inf_transformation is deprecated. All integral problems are transformed"
+    _alg = alg isa ChangeOfVariables ? alg : ChangeOfVariables(transformation_if_inf, alg)
+    cacheval = init_cacheval(_alg, prob)
 
     IntegralCache{iip,
         typeof(prob.f),
         typeof(prob.domain),
         typeof(prob.p),
         typeof(prob.kwargs),
-        typeof(alg),
+        typeof(_alg),
         typeof(sensealg),
         typeof(kwargs),
         typeof(cacheval)}(Val(iip),
@@ -35,7 +37,7 @@ function SciMLBase.init(prob::IntegralProblem{iip},
         prob.domain,
         prob.p,
         prob.kwargs,
-        alg,
+        _alg,
         sensealg,
         kwargs,
         cacheval)
@@ -102,7 +104,7 @@ function SciMLBase.solve(prob::IntegralProblem,
 end
 
 function SciMLBase.solve!(cache::IntegralCache)
-    __solvebp(cache, cache.alg, cache.sensealg, cache.domain, cache.p;
+    __solve(cache, cache.alg, cache.sensealg, cache.domain, cache.p;
         cache.kwargs...)
 end
 
diff --git a/src/infinity_handling.jl b/src/infinity_handling.jl
index 25cfa69..b822362 100644
--- a/src/infinity_handling.jl
+++ b/src/infinity_handling.jl
@@ -1,7 +1,94 @@
-_oftype(x, y) = oftype(x, y)
-_oftype(::SubArray, y) = y
+# generic change of variables code
 
-function substitute_bounds(lb, ub)
+function substitute_u(u2v, lb, ub)
+    bounds = map(u2v, lb, ub)
+    lb_sub = lb isa Number ? first(bounds) : map(first, bounds)
+    ub_sub = ub isa Number ? last(bounds) : map(last, bounds)
+    return (lb_sub, ub_sub)
+end
+
+# without batching point container type should match the inputs
+substitute_v(v2ujac, v, lb::Number, ub::Number) = v2ujac(only(v), lb, ub)
+function substitute_v(v2ujac, v, lb::AbstractVector, ub::AbstractVector)
+    xjac = map((l, u, v) -> v2ujac(v, l, u), lb, ub, v) # ordering may influence container type
+    x = map(first, xjac)
+    jac = prod(last, xjac)
+    return x, jac
+end
+
+function substitute_f(f::IntegralFunction{false}, v2ujac, lb, ub)
+    _f = f.f
+    IntegralFunction{false}(f.integrand_prototype) do v, p
+        u, jac = substitute_v(v2ujac, v, lb, ub)
+        return _f(u, p) * jac
+    end
+end
+function substitute_f(f::IntegralFunction{true}, v2ujac, lb, ub)
+    _f = f.f
+    prototype = similar(f.integrand_prototype)
+    vol = prod((ub - lb) / 2) # just to get the type of the jacobian determinant
+    IntegralFunction{true}(prototype * vol) do y, v, p
+        u, jac = substitute_v(v2ujac, v, lb, ub)
+        _y = _evaluate!(f, prototype, u, p)
+        y .= _y .* jac
+        return
+    end
+end
+
+# with batching the point container type is assumed to be mutable
+function substitute_bv(v2ujac, v::AbstractArray, lb::Number, ub::Number)
+    x = similar(v, typeof(one(eltype(v)) * (first(lb) + first(ub))))
+    jac = similar(x)
+    for i in axes(v, 1)
+        x[i], jac[i] = v2ujac(v[i], lb, ub)
+    end
+    return x, jac
+end
+function substitute_bv(v2ujac, v::AbstractArray, lb::AbstractVector, ub::AbstractVector)
+    x = similar(v, typeof(one(eltype(v)) * (first(lb) + first(ub))))
+    jac = similar(v, typeof(zero(eltype(v)) * prod(lb)), size(v)[end])
+    idx = CartesianIndices(axes(v)[begin:(end - 1)])
+    for i in axes(v)[end]
+        _jac = one(eltype(jac))
+        for (ii, l, u) in zip(idx, lb, ub)
+            x[ii, i], j = v2ujac(v[ii, i], l, u)
+            _jac *= j
+        end
+        jac[i] = _jac
+    end
+    return x, jac
+end
+
+function substitute_f(f::BatchIntegralFunction{false}, v2ujac, lb, ub)
+    _f = f.f
+    BatchIntegralFunction{false}(f.integrand_prototype, max_batch = f.max_batch) do v, p
+        u, jac = substitute_bv(v2ujac, v, lb, ub)
+        y = _f(u, p)
+        return y .* reshape(jac, ntuple(d -> d == ndims(y) ? length(jac) : 1, ndims(y)))
+    end
+end
+function substitute_f(f::BatchIntegralFunction{true}, v2ujac, lb, ub)
+    _f = f.f
+    prototype = similar(f.integrand_prototype)
+    vol = prod((ub - lb) / 2) # just to get the type of the jacobian determinant
+    BatchIntegralFunction{true}(prototype * vol, max_batch = f.max_batch) do y, v, p
+        u, jac = substitute_bv(v2ujac, v, lb, ub)
+        _prototype = similar(prototype, size(y))
+        _y = _evaluate!(_f, _prototype, u, p)
+        y .= _y .* reshape(jac, ntuple(d -> d == ndims(y) ? length(jac) : 1, ndims(y)))
+        return
+    end
+end
+
+# we need this function for autodiff compatibility where internal buffers need special types
+function _evaluate!(f, y, u, p)
+    f(y, u, p)
+    return y
+end
+
+# specific changes of variables
+
+function u2t(lb, ub)
     mid = (lb + ub) / 2 # floating-point promotion
     if isinf(lb) && isinf(ub)
         lb_sub = flipsign(one(mid), lb)
@@ -19,7 +106,7 @@ function substitute_bounds(lb, ub)
     return lb_sub, ub_sub # unitless
 end
 
-function substitute_t(t::Number, lb::Number, ub::Number)
+function t2ujac(t::Number, lb::Number, ub::Number)
     u = oneunit(eltype(lb))
     # apply correct units
     if isinf(lb) && isinf(ub)
@@ -32,114 +119,17 @@ function substitute_t(t::Number, lb::Number, ub::Number)
         den = inv(1 - flipsign(t, ub))
         lb + t * den * u, den^2 * u
     else
-        den = (ub - lb) * oftype(t, 0.5)
+        den = (ub - lb) * oftype(float(one(u)), 0.5)
         lb + (1 + t) * den, den
     end
 end
-function substitute_t(t::AbstractVector, lb::AbstractVector, ub::AbstractVector)
-    x = similar(t, typeof(one(eltype(t)) * (first(lb) + first(ub))))
-    jac = one(eltype(t))
-    for i in eachindex(lb)
-        x[i], dj = substitute_t(t[i], lb[i], ub[i])
-        jac *= dj
-    end
-    return _oftype(t, x), jac
-end
-
-function substitute_f(f::F, t, p, lb, ub) where {F}
-    x, jac = substitute_t(t, lb, ub)
-    return f(x, p) * jac
-end
-function substitute_f(f::F, dt, t, p, lb, ub) where {F}
-    x, jac = substitute_t(t, lb, ub)
-    f(dt, x, p)
-    dt .*= jac
-    return
-end
-
-function substitute_t(t::AbstractVector, lb::Number, ub::Number)
-    x = similar(t, typeof(one(eltype(t)) * (lb + ub)))
-    jac = similar(x)
-    for (i, ti) in enumerate(t)
-        x[i], jac[i] = substitute_t(ti, lb, ub)
-    end
-    return x, jac
-end
-function substitute_t(t::AbstractArray, lb::AbstractVector, ub::AbstractVector)
-    x = similar(t, typeof(one(eltype(t)) * (first(lb) + first(ub))))
-    jac = similar(x, size(t, ndims(t)))
-    for (i, it) in enumerate(axes(t)[end])
-        x[axes(x)[begin:(end - 1)]..., i], jac[i] = substitute_t(
-            t[axes(t)[begin:(end - 1)]..., it], lb, ub)
-    end
-    return x, jac
-end
-
-function substitute_batchf(f::F, t, p, lb, ub) where {F}
-    x, jac = substitute_t(t, lb, ub)
-    r = f(x, p)
-    return r .* reshape(jac, ntuple(d -> d == ndims(r) ? length(jac) : 1, ndims(r)))
-end
-function substitute_batchf(f::F, dt, t, p, lb, ub) where {F}
-    x, jac = substitute_t(t, lb, ub)
-    f(dt, x, p)
-    for (i, j) in zip(axes(dt)[end], jac)
-        for idt in CartesianIndices(axes(dt)[begin:(end - 1)])
-            dt[idt, i] *= j
-        end
-    end
-    return
-end
 
-function transformation_if_inf(prob, ::Val{true})
-    lb, ub = promote(prob.domain...)
-    f = prob.f
-    bounds = map(substitute_bounds, lb, ub)
-    lb_sub = lb isa Number ? first(bounds) : map(first, bounds)
-    ub_sub = ub isa Number ? last(bounds) : map(last, bounds)
-    f_sub = if isinplace(prob)
-        if f isa BatchIntegralFunction
-            BatchIntegralFunction{true}(
-                let f = f.f
-                    (dt, t, p) -> substitute_batchf(f, dt, t, p, lb, ub)
-                end,
-                f.integrand_prototype,
-                max_batch = f.max_batch)
-        else
-            IntegralFunction{true}(let f = f.f
-                    (dt, t, p) -> substitute_f(f, dt, t, p, lb, ub)
-                end,
-                f.integrand_prototype)
-        end
-    else
-        if f isa BatchIntegralFunction
-            BatchIntegralFunction{false}(let f = f.f
-                    (t, p) -> substitute_batchf(f, t, p, lb, ub)
-                end,
-                f.integrand_prototype)
-        else
-            IntegralFunction{false}(let f = f.f
-                    (t, p) -> substitute_f(f, t, p, lb, ub)
-                end,
-                f.integrand_prototype)
-        end
-    end
-    return remake(prob, f = f_sub, domain = (lb_sub, ub_sub))
+function transformation_if_inf(f, domain)
+    lb, ub = promote(domain...)
+    tdomain = substitute_u(u2t, lb, ub)
+    g = substitute_f(f, t2ujac, lb, ub)
+    return g, tdomain
 end
 
-function transformation_if_inf(prob, ::Nothing)
-    lb, ub = prob.domain
-    if any(isinf, lb) || any(isinf, ub)
-        return transformation_if_inf(prob, Val(true))
-    else
-        return transformation_if_inf(prob, Val(false))
-    end
-end
-
-function transformation_if_inf(prob, ::Val{false})
-    return prob
-end
-
-function transformation_if_inf(prob, do_inf_transformation = nothing)
-    transformation_if_inf(prob, do_inf_transformation)
-end
+# to implement more transformations, define functions u2v and v2ujac and then a wrapper
+# similar to transformation_if_inf to pass to ChangeOfVariables
diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl
index ba22e8f..9f0cd35 100644
--- a/test/derivative_tests.jl
+++ b/test/derivative_tests.jl
@@ -18,8 +18,8 @@ alg_req = Dict(
     QuadGKJL() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = 1,
         allows_iip = true),
     HCubatureJL() => (nout = Inf, allows_batch = false, min_dim = 1,
-        max_dim = Inf, allows_iip = true), CubatureJLh() => (
-        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))
diff --git a/test/gaussian_quadrature_tests.jl b/test/gaussian_quadrature_tests.jl
index d3a2387..ad8eb41 100644
--- a/test/gaussian_quadrature_tests.jl
+++ b/test/gaussian_quadrature_tests.jl
@@ -87,7 +87,7 @@ alg = GaussLegendre(subintervals = 1)
 alg = GaussLegendre(subintervals = 17)
 @test sol.u ≈ sqrt(π) / 2
 
-# Make sure broadcasting correctly handles the argument p 
+# Make sure broadcasting correctly handles the argument p
 f = (x, p) -> 1 + x + x^p[1] - cos(x * p[2]) + exp(x) * p[3]
 p = [0.3, 1.3, -0.5]
 prob = IntegralProblem(f, 2.0, 6.3, p)
diff --git a/test/inf_integral_tests.jl b/test/inf_integral_tests.jl
index c0763e1..ee41083 100644
--- a/test/inf_integral_tests.jl
+++ b/test/inf_integral_tests.jl
@@ -7,12 +7,12 @@ abstol = 1e-3
 alg_req = Dict(
     QuadratureRule(gausslegendre, n = 50) => (
         nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false,
-        allows_iip = false, allows_inf = true), QuadGKJL() => (
-        nout = Inf, allows_batch = true, min_dim = 1, max_dim = 1,
+        allows_iip = false, allows_inf = true),
+    QuadGKJL() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = 1,
         allows_iip = true, allows_inf = true),
     HCubatureJL() => (nout = Inf, allows_batch = false, min_dim = 1,
-        max_dim = Inf, allows_iip = true, allows_inf = true), CubatureJLh() => (
-        nout = Inf, allows_batch = true, min_dim = 1,
+        max_dim = Inf, allows_iip = true, allows_inf = true),
+    CubatureJLh() => (nout = Inf, allows_batch = true, min_dim = 1,
         max_dim = Inf, allows_iip = true, allows_inf = true)
 )
 # GaussLegendre(n=50) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false,
@@ -148,11 +148,11 @@ end
 
 do_tests = function (; f, domain, alg, abstol, reltol, solution)
     prob = IntegralProblem(f, domain)
-    sol = solve(prob, alg; reltol, abstol, do_inf_transformation = Val(true))
+    sol = solve(prob, alg; reltol, abstol)
     @test abs(only(sol) - solution) < max(abstol, reltol * abs(solution))
-    cache = @test_nowarn @inferred init(prob, alg; do_inf_transformation = Val(true))
+    cache = @test_nowarn @inferred init(prob, alg)
     @test_nowarn @inferred solve!(cache)
-    @test_nowarn @inferred solve(prob, alg; do_inf_transformation = Val(true))
+    @test_nowarn @inferred solve(prob, alg)
 end
 
 # IntegralFunction{false}
@@ -201,3 +201,24 @@ for (alg, req) in pairs(alg_req), (j, (; f, domain, solution)) in enumerate(prob
     bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(0))
     do_tests(; f = bfiip, domain, solution, alg, abstol, reltol)
 end
+
+@testset "Caching interface" begin
+    # two distinct semi-infinite transformations should still work as expected
+    f = (x, p) -> pdf(Normal(0.00, 1.00), x)
+    domain = (0.0, -Inf)
+    solution = -0.5
+    prob = IntegralProblem(f, domain)
+    alg = QuadGKJL()
+    cache = init(prob, alg; abstol, reltol)
+    sol = solve!(cache)
+    @test abs(only(sol.u) - solution) < max(abstol, reltol * abs(solution))
+    @test sol.prob == IntegralProblem(f, domain)
+    @test sol.alg == alg
+    domain = (-Inf, 0.0)
+    solution = 0.5
+    cache.domain = domain
+    sol = solve!(cache)
+    @test abs(only(sol.u) - solution) < max(abstol, reltol * abs(solution))
+    @test sol.prob == IntegralProblem(f, domain)
+    @test sol.alg == alg
+end
diff --git a/test/interface_tests.jl b/test/interface_tests.jl
index 272aed5..920681e 100644
--- a/test/interface_tests.jl
+++ b/test/interface_tests.jl
@@ -17,8 +17,8 @@ alg_req = Dict(
         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),
+    VEGAS(seed = 109) => (nout = 1, allows_batch = true, min_dim = 2, max_dim = Inf,
+        allows_iip = true),
     VEGASMC(seed = 42) => (nout = Inf, allows_batch = false, min_dim = 1, max_dim = Inf,
         allows_iip = true),
     CubatureJLh() => (nout = Inf, allows_batch = true, min_dim = 1,
@@ -34,7 +34,7 @@ alg_req = Dict(
     CubaCuhre() => (nout = Inf, allows_batch = true, min_dim = 2, max_dim = Inf,
         allows_iip = true),
     ArblibJL() => (
-        nout = 1, allows_batch = false, min_dim = 1, max_dim = 1, allows_iip = true)
+        nout = 1, allows_batch = false, min_dim = 1, max_dim = 1, allows_iip = false)
 )
 
 integrands = [
@@ -44,7 +44,9 @@ integrands = [
 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)]
+                let f = integrands[2]
+                    (x, p, nout) -> f(x, p) * collect(1.0:nout)
+                end]
 iip_integrands_v = [(dx, x, p, nout) -> (dx .= f(x, p, nout)) for f in integrands_v]
 
 exact_sol = [
@@ -103,7 +105,7 @@ end
         for i in 1:length(integrands)
             prob = IntegralProblem(integrands[i], lb, ub)
             @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-            sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+            sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol)
             @test sol.alg == alg
             @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2
         end
@@ -121,7 +123,7 @@ end
                     continue
                 end
                 @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol)
                 @test sol.alg == alg
                 @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2
             end
@@ -140,7 +142,7 @@ end
                     continue
                 end
                 @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol)
                 @test sol.alg == alg
                 if sol.u isa Number
                     @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2
@@ -162,7 +164,7 @@ end
                 continue
             end
             @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-            sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+            sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol)
             @test sol.alg == alg
             @test sol.u[1]≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2
         end
@@ -180,7 +182,7 @@ end
                     continue
                 end
                 @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol)
                 @test sol.alg == alg
                 if sol.u isa Number
                     @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2
@@ -204,7 +206,7 @@ end
                     continue
                 end
                 @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol)
                 @test sol.alg == alg
                 if sol.u isa Number
                     @test sol.u≈exact_sol[i](dim, nout, lb, ub) rtol=1e-2
@@ -223,13 +225,15 @@ end
     for (alg, req) in pairs(alg_req)
         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,
+                prob = IntegralProblem(let f = integrands_v[i], nout = nout
+                        (x, p) -> f(x, p, nout)
+                    end, lb, ub,
                     nout = nout)
                 if req.min_dim > 1 || req.nout < nout
                     continue
                 end
                 @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol)
                 @test sol.alg == alg
                 if nout == 1
                     @test sol.u[1]≈exact_sol_v[i](dim, nout, lb, ub)[1] rtol=1e-2
@@ -250,10 +254,12 @@ end
                     if dim > req.max_dim || dim < req.min_dim || req.nout < nout
                         continue
                     end
-                    prob = IntegralProblem((x, p) -> integrands_v[i](x, p, nout), lb, ub,
+                    prob = IntegralProblem(let f = integrands_v[i], nout = nout
+                            (x, p) -> f(x, p, nout)
+                        end, lb, ub,
                         nout = nout)
                     @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                    sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                    sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol)
                     @test sol.alg == alg
                     if nout == 1
                         @test sol.u[1]≈exact_sol_v[i](dim, nout, lb, ub)[1] rtol=1e-2
@@ -272,15 +278,16 @@ end
             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),
+                    prob = IntegralProblem(let f = iip_integrands_v[i]
+                            (dx, x, p) -> f(dx, x, p, nout)
+                        end,
                         lb, ub, nout = nout)
                     if dim > req.max_dim || dim < req.min_dim || req.nout < nout ||
                        !req.allows_iip
                         continue
                     end
                     @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                    sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                    sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol)
                     @test sol.alg == alg
                     if nout == 1
                         @test sol.u[1]≈exact_sol_v[i](dim, nout, lb, ub)[1] rtol=1e-2
@@ -304,7 +311,7 @@ end
                 continue
             end
             @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-            sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+            sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol)
             @test sol.alg == alg
             @test sol.u≈exact_sol_v[i](dim, nout, lb, ub) rtol=1e-2
         end
@@ -325,7 +332,7 @@ end
                     continue
                 end
                 @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol)
                 @test sol.alg == alg
                 @test sol.u≈exact_sol_v[i](dim, nout, lb, ub) rtol=1e-2
             end
@@ -346,7 +353,7 @@ end
                     continue
                 end
                 @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout"
-                sol = solve(prob, alg, reltol = reltol, abstol = abstol)
+                sol = @inferred solve(prob, alg, reltol = reltol, abstol = abstol)
                 @test sol.alg == alg
                 @test sol.u≈exact_sol_v[i](dim, nout, lb, ub) rtol=1e-2
             end
diff --git a/test/nested_ad_tests.jl b/test/nested_ad_tests.jl
index d4506aa..65840ea 100644
--- a/test/nested_ad_tests.jl
+++ b/test/nested_ad_tests.jl
@@ -8,8 +8,8 @@ function my_integration(p)
     return solve(my_problem, CubatureJLh(), reltol = 1e-3, abstol = 1e-3)  # Errors
 end
 my_solution = my_integration(my_parameters)
-@test ForwardDiff.jacobian(my_integration, my_parameters) == [0.0 8.0]
-@test ForwardDiff.jacobian(x -> ForwardDiff.jacobian(my_integration, x), my_parameters) ==
+@test ForwardDiff.jacobian(my_integration, my_parameters) ≈ [0.0 8.0]
+@test ForwardDiff.jacobian(x -> ForwardDiff.jacobian(my_integration, x), my_parameters) ≈
       [0.0 0.0
        0.0 4.0]
 
diff --git a/test/qa.jl b/test/qa.jl
index ccf1183..9992475 100644
--- a/test/qa.jl
+++ b/test/qa.jl
@@ -1,4 +1,5 @@
 using Integrals, Aqua
+using Test
 @testset "Aqua" begin
     Aqua.find_persistent_tasks_deps(Integrals)
     Aqua.test_ambiguities(Integrals, recursive = false)
diff --git a/test/runtests.jl b/test/runtests.jl
index 090b8fd..bfdbfe9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,4 +1,3 @@
-using Pkg
 using SafeTestsets
 using Test