diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
index 9c79359..959ad88 100644
--- a/.JuliaFormatter.toml
+++ b/.JuliaFormatter.toml
@@ -1,2 +1,3 @@
 style = "sciml"
-format_markdown = true
\ No newline at end of file
+format_markdown = true
+format_docstrings = true
\ No newline at end of file
diff --git a/docs/pages.jl b/docs/pages.jl
index d8c8a7d..594e50e 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -6,5 +6,5 @@ pages = ["index.md",
         "basics/SampledIntegralProblem.md",
         "basics/solve.md",
         "basics/FAQ.md"],
-    "Solvers" => Any["solvers/IntegralSolvers.md"],
+    "Solvers" => Any["solvers/IntegralSolvers.md"]
 ]
diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md
index baaff81..4cd33aa 100644
--- a/docs/src/basics/FAQ.md
+++ b/docs/src/basics/FAQ.md
@@ -40,8 +40,9 @@ For badly-behaved integrands, such as (nearly) singular and highly oscillatory
 functions, most algorithms will fail to converge and either throw an error or
 silently return the incorrect result. In some cases Integrals.jl can provide an
 error code when things go wrong, but otherwise you can always check if the error
-estimate for the integral is less than the requested tolerance, e.g. `sol.resid
-< max(abstol, reltol*norm(sol.u))`. Sometimes using a larger tolerance or higher
+estimate for the integral is less than the requested tolerance, e.g.
+`sol.resid < max(abstol, reltol*norm(sol.u))`.
+Sometimes using a larger tolerance or higher
 precision arithmetic may help.
 
 ## How can I integrate arbitrarily-spaced data?
@@ -60,4 +61,4 @@ Otherwise, feel free to open an issue or pull request.
 
 ## Can I take derivatives with respect to the limits of integration?
 
-Currently this is not implemented.
\ No newline at end of file
+Currently this is not implemented.
diff --git a/docs/src/basics/IntegralFunction.md b/docs/src/basics/IntegralFunction.md
index eee1709..58a67db 100644
--- a/docs/src/basics/IntegralFunction.md
+++ b/docs/src/basics/IntegralFunction.md
@@ -3,4 +3,4 @@
 ```@docs
 SciMLBase.IntegralFunction
 SciMLBase.BatchIntegralFunction
-```
\ No newline at end of file
+```
diff --git a/docs/src/basics/SampledIntegralProblem.md b/docs/src/basics/SampledIntegralProblem.md
index fd2aa84..f856c68 100644
--- a/docs/src/basics/SampledIntegralProblem.md
+++ b/docs/src/basics/SampledIntegralProblem.md
@@ -12,7 +12,7 @@ Say, by some means we have generated a dataset `x` and `y`:
 ```@example 1
 using Integrals # hide
 f = x -> x^2
-x = range(0, 1, length=20)
+x = range(0, 1, length = 20)
 y = f.(x)
 ```
 
@@ -63,9 +63,9 @@ using Integrals # hide
 f1 = x -> x^2
 f2 = x -> x^3
 f3 = x -> x^4
-x = range(0, 1, length=20)
+x = range(0, 1, length = 20)
 y = [f1.(x) f2.(x) f3.(x)]
-problem = SampledIntegralProblem(y, x; dim=1)
+problem = SampledIntegralProblem(y, x; dim = 1)
 method = SimpsonsRule()
 solve(problem, method)
 ```
diff --git a/docs/src/tutorials/differentiating_integrals.md b/docs/src/tutorials/differentiating_integrals.md
index 763f8a0..e2c4b5e 100644
--- a/docs/src/tutorials/differentiating_integrals.md
+++ b/docs/src/tutorials/differentiating_integrals.md
@@ -16,7 +16,7 @@ function testf(p)
     sin(solve(prob, CubaCuhre(), reltol = 1e-6, abstol = 1e-6)[1])
 end
 testf(p)
-dp1 = Zygote.gradient(testf,p)
+dp1 = Zygote.gradient(testf, p)
 dp2 = FiniteDiff.finite_difference_gradient(testf, p)
 dp3 = ForwardDiff.gradient(testf, p)
 dp1[1] ≈ dp2 ≈ dp3
diff --git a/ext/IntegralsArblibExt.jl b/ext/IntegralsArblibExt.jl
index 37811bd..e0c408a 100644
--- a/ext/IntegralsArblibExt.jl
+++ b/ext/IntegralsArblibExt.jl
@@ -3,9 +3,9 @@ module IntegralsArblibExt
 using Arblib
 using Integrals
 
-function Integrals.__solvebp_call(prob::IntegralProblem, alg::ArblibJL, sensealg, domain, p;
-    reltol = 1e-8, abstol = 1e-8, maxiters = nothing)
-
+function Integrals.__solvebp_call(
+        prob::IntegralProblem, alg::ArblibJL, sensealg, domain, p;
+        reltol = 1e-8, abstol = 1e-8, maxiters = nothing)
     lb_, ub_ = domain
     lb, ub = map(first, domain)
     if !isone(length(lb_)) || !isone(length(ub_))
@@ -17,16 +17,18 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::ArblibJL, sensealg
         res = Acb(0)
         y_ = similar(prob.f.integrand_prototype, typeof(res))
         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)
+        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)
+        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
 end
 
diff --git a/ext/IntegralsCubaExt.jl b/ext/IntegralsCubaExt.jl
index e391c65..f5c1697 100644
--- a/ext/IntegralsCubaExt.jl
+++ b/ext/IntegralsCubaExt.jl
@@ -2,8 +2,8 @@ module IntegralsCubaExt
 
 using Integrals, Cuba
 import Integrals: transformation_if_inf,
-    scale_x, scale_x!, CubaVegas, AbstractCubaAlgorithm,
-    CubaSUAVE, CubaDivonne, CubaCuhre
+                  scale_x, scale_x!, CubaVegas, AbstractCubaAlgorithm,
+                  CubaSUAVE, CubaDivonne, CubaCuhre
 
 function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgorithm,
         sensealg,
@@ -22,7 +22,8 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgori
         nvec = min(maxiters, prob.f.max_batch)
         # nvec == 1 in Cuba will change vectors to matrices, so we won't support it when
         # batching
-        nvec > 1 || throw(ArgumentError("BatchIntegralFunction must take multiple batch points"))
+        nvec > 1 ||
+            throw(ArgumentError("BatchIntegralFunction must take multiple batch points"))
 
         if mid isa Real
             _x = zeros(typeof(mid), nvec)
diff --git a/ext/IntegralsFastGaussQuadratureExt.jl b/ext/IntegralsFastGaussQuadratureExt.jl
index f3868bf..1dc5c39 100644
--- a/ext/IntegralsFastGaussQuadratureExt.jl
+++ b/ext/IntegralsFastGaussQuadratureExt.jl
@@ -33,7 +33,7 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::Integrals.GaussLeg
         sensealg, domain, p;
         reltol = nothing, abstol = nothing,
         maxiters = nothing) where {C}
-    if !all(isone∘length, domain)
+    if !all(isone ∘ length, domain)
         error("GaussLegendre only accepts one-dimensional quadrature problems.")
     end
     @assert prob.f isa IntegralFunction
diff --git a/ext/IntegralsForwardDiffExt.jl b/ext/IntegralsForwardDiffExt.jl
index 300b606..49abe0b 100644
--- a/ext/IntegralsForwardDiffExt.jl
+++ b/ext/IntegralsForwardDiffExt.jl
@@ -22,8 +22,8 @@ end
 
 # Manually split for the pushforward
 function Integrals.__solvebp(cache, alg, sensealg, domain,
-        p::Union{D,AbstractArray{<:D}};
-        kwargs...) where {T, V, P, D<:ForwardDiff.Dual{T, V, P}}
+        p::Union{D, AbstractArray{<:D}};
+        kwargs...) where {T, V, P, D <: ForwardDiff.Dual{T, V, P}}
 
     # we need the output type to avoid perturbation confusion while unwrapping nested duals
     # We compute a vector-valued integral of the primal and dual simultaneously
diff --git a/ext/IntegralsMCIntegrationExt.jl b/ext/IntegralsMCIntegrationExt.jl
index fb9a68c..e4d6779 100644
--- a/ext/IntegralsMCIntegrationExt.jl
+++ b/ext/IntegralsMCIntegrationExt.jl
@@ -3,10 +3,10 @@ module IntegralsMCIntegrationExt
 using MCIntegration, Integrals
 
 function Integrals.__solvebp_call(prob::IntegralProblem, alg::VEGASMC, sensealg, domain, p;
-    reltol = nothing, abstol = nothing, maxiters = 10)
+        reltol = nothing, abstol = nothing, maxiters = 10)
     lb, ub = domain
     mid = vec(collect((lb + ub) / 2))
-    var = Continuous(vec([tuple(a,b) for (a,b) in zip(lb, ub)]))
+    var = Continuous(vec([tuple(a, b) for (a, b) in zip(lb, ub)]))
 
     if prob.f isa BatchIntegralFunction
         error("VEGASMC doesn't support batching. See https://github.com/numericalEFT/MCIntegration.jl/issues/29")
@@ -16,7 +16,7 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::VEGASMC, sensealg,
             f_ = (x, f, c) -> begin
                 n = 0
                 for v in x
-                    mid[n+=1] = first(v)
+                    mid[n += 1] = first(v)
                 end
                 prob.f(f0, mid, p)
                 f .= vec(f0)
@@ -26,21 +26,22 @@ function Integrals.__solvebp_call(prob::IntegralProblem, alg::VEGASMC, sensealg,
             f_ = (x, c) -> begin
                 n = 0
                 for v in x
-                    mid[n+=1] = first(v)
+                    mid[n += 1] = first(v)
                 end
                 fx = prob.f(mid, p)
                 fx isa AbstractArray ? vec(fx) : fx
             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),
-            solver=:vegasmc, niter=maxiters, verbose=-2,  print=-2, alg.kws...)
+        res = integrate(f_; var, dof, inplace = isinplace(prob), type = eltype(f0),
+            solver = :vegasmc, niter = maxiters, verbose = -2, print = -2, alg.kws...)
         out, err, chi = if f0 isa Number
             map(only, (res.mean, res.stdev, res.chi2))
         else
             map(a -> reshape(a, size(f0)), (res.mean, res.stdev, res.chi2))
         end
-        SciMLBase.build_solution(prob, alg, out, err, chi=chi, retcode = ReturnCode.Success)
+        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 90d6f88..1b38956 100644
--- a/ext/IntegralsZygoteExt.jl
+++ b/ext/IntegralsZygoteExt.jl
@@ -28,7 +28,8 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal
             # zygote doesn't support mutation, so we build an oop pullback
             if sensealg.vjp isa Integrals.ZygoteVJP
                 if cache.f isa BatchIntegralFunction
-                    dx = similar(cache.f.integrand_prototype, size(cache.f.integrand_prototype)[begin:end-1]..., 1)
+                    dx = similar(cache.f.integrand_prototype,
+                        size(cache.f.integrand_prototype)[begin:(end - 1)]..., 1)
                     _f = x -> (cache.f(dx, x, p); dx)
                     # TODO: let the user pass a batched jacobian so we can return a BatchIntegralFunction
                     dfdp_ = function (x, p)
@@ -104,7 +105,7 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal
             # TODO replace evaluation at endpoint (which anyone can do without Integrals.jl)
             # with integration of dfdx uing the same quadrature
             dlb = cache.f isa BatchIntegralFunction ? -batch_unwrap(_f([lb])) : -_f(lb)
-            dub = cache.f isa BatchIntegralFunction ?  batch_unwrap(_f([ub])) : _f(ub)
+            dub = cache.f isa BatchIntegralFunction ? batch_unwrap(_f([ub])) : _f(ub)
             return (NoTangent(),
                 NoTangent(),
                 NoTangent(),
@@ -128,7 +129,7 @@ function ChainRulesCore.rrule(::typeof(Integrals.__solvebp), cache, alg, senseal
     out, quadrature_adjoint
 end
 
-batch_unwrap(x::AbstractArray) = dropdims(x; dims=ndims(x))
+batch_unwrap(x::AbstractArray) = dropdims(x; dims = ndims(x))
 
 Zygote.@adjoint function Zygote.literal_getproperty(sol::SciMLBase.IntegralSolution,
         ::Val{:u})
diff --git a/src/Integrals.jl b/src/Integrals.jl
index c503827..91d4772 100644
--- a/src/Integrals.jl
+++ b/src/Integrals.jl
@@ -66,7 +66,6 @@ end
 # 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
@@ -106,10 +105,10 @@ function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, domain, p
                 BatchIntegrand((y, x) -> prob.f(y, x, p), similar(bu))
             else
                 fsize = size(bu)[begin:(end - 1)]
-                BatchIntegrand{Array{eltype(bu),ndims(bu)-1}}() do y, x
+                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)))
+                    map!(collect, y, eachslice(y_; dims = ndims(bu)))
                     return nothing
                 end
             end
@@ -121,8 +120,8 @@ function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, domain, p
             f = if u isa AbstractVector
                 BatchIntegrand((y, x) -> y .= prob.f(x, p), u)
             else
-                BatchIntegrand{Array{eltype(u),ndims(u)-1}}() do y, x
-                    map!(collect, y, eachslice(prob.f(x, p); dims=ndims(u)))
+                BatchIntegrand{Array{eltype(u), ndims(u) - 1}}() do y, x
+                    map!(collect, y, eachslice(prob.f(x, p); dims = ndims(u)))
                     return nothing
                 end
             end
@@ -134,7 +133,8 @@ function __solvebp_call(cache::IntegralCache, alg::QuadGKJL, sensealg, domain, p
         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,
+            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
@@ -189,8 +189,9 @@ function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, domain, p;
             # 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))
         else
-            y = prob.f(mid isa Number ? typeof(mid)[] :
-                       Matrix{eltype(mid)}(undef, length(mid), 0),
+            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)
         end
@@ -221,7 +222,8 @@ function __solvebp_call(prob::IntegralProblem, alg::VEGAS, sensealg, domain, p;
     SciMLBase.build_solution(prob, alg, val, err, chi = chi, retcode = ReturnCode.Success)
 end
 
-export QuadGKJL, HCubatureJL, VEGAS, VEGASMC, GaussLegendre, QuadratureRule, TrapezoidalRule, SimpsonsRule
+export QuadGKJL, HCubatureJL, VEGAS, VEGASMC, GaussLegendre, QuadratureRule,
+       TrapezoidalRule, SimpsonsRule
 export CubaVegas, CubaSUAVE, CubaDivonne, CubaCuhre
 export CubatureJLh, CubatureJLp
 export ArblibJL
diff --git a/src/algorithms.jl b/src/algorithms.jl
index 9f1bef7..4222b99 100644
--- a/src/algorithms.jl
+++ b/src/algorithms.jl
@@ -127,9 +127,8 @@ function GaussLegendre(; n = 250, subintervals = 1, nodes = nothing, weights = n
     return GaussLegendre(nodes, weights, subintervals)
 end
 
-
 """
-  QuadratureRule(q; n=250)
+QuadratureRule(q; n=250)
 
 Algorithm to construct and evaluate a quadrature rule `q` of `n` points computed from the
 inputs as `x, w = q(n)`. It assumes the nodes and weights are for the standard interval
diff --git a/src/algorithms_extension.jl b/src/algorithms_extension.jl
index b79eaa8..995f712 100644
--- a/src/algorithms_extension.jl
+++ b/src/algorithms_extension.jl
@@ -123,13 +123,15 @@ end
 
 function CubaVegas(; flags = 0, seed = 0, minevals = 0, nstart = 1000, nincrease = 500,
         gridno = 0)
-    isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) && error("CubaVegas requires `using Cuba`")
+    isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) &&
+        error("CubaVegas requires `using Cuba`")
     return CubaVegas(flags, seed, minevals, nstart, nincrease, gridno)
 end
 
 function CubaSUAVE(; flags = 0, seed = 0, minevals = 0, nnew = 1000, nmin = 2,
         flatness = 25.0)
-    isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) && error("CubaSUAVE requires `using Cuba`")
+    isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) &&
+        error("CubaSUAVE requires `using Cuba`")
     return CubaSUAVE(flags, seed, minevals, nnew, nmin, flatness)
 end
 
@@ -138,13 +140,15 @@ function CubaDivonne(; flags = 0, seed = 0, minevals = 0,
         maxchisq = 10.0, mindeviation = 0.25,
         xgiven = zeros(Cdouble, 0, 0),
         nextra = 0, peakfinder = C_NULL)
-    isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) && error("CubaDivonne requires `using Cuba`")
+    isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) &&
+        error("CubaDivonne requires `using Cuba`")
     return CubaDivonne(flags, seed, minevals, key1, key2, key3, maxpass, border, maxchisq,
         mindeviation, xgiven, nextra, peakfinder)
 end
 
 function CubaCuhre(; flags = 0, minevals = 0, key = 0)
-    isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) && error("CubaCuhre requires `using Cuba`")
+    isnothing(Base.get_extension(@__MODULE__, :IntegralsCubaExt)) &&
+        error("CubaCuhre requires `using Cuba`")
     return CubaCuhre(flags, minevals, key)
 end
 
@@ -174,8 +178,9 @@ publisher={Elsevier}
 struct CubatureJLh <: AbstractCubatureJLAlgorithm
     error_norm::Int32
 end
-function CubatureJLh(; error_norm=0)
-    isnothing(Base.get_extension(@__MODULE__, :IntegralsCubatureExt)) && error("CubatureJLh requires `using Cubature`")
+function CubatureJLh(; error_norm = 0)
+    isnothing(Base.get_extension(@__MODULE__, :IntegralsCubatureExt)) &&
+        error("CubatureJLh requires `using Cubature`")
     return CubatureJLh(error_norm)
 end
 
@@ -193,13 +198,12 @@ Defaults to `Cubature.INDIVIDUAL`, other options are
 struct CubatureJLp <: AbstractCubatureJLAlgorithm
     error_norm::Int32
 end
-function CubatureJLp(; error_norm=0)
-    isnothing(Base.get_extension(@__MODULE__, :IntegralsCubatureExt)) && error("CubatureJLp requires `using Cubature`")
+function CubatureJLp(; error_norm = 0)
+    isnothing(Base.get_extension(@__MODULE__, :IntegralsCubatureExt)) &&
+        error("CubatureJLp requires `using Cubature`")
     return CubatureJLp(error_norm)
 end
 
-
-
 """
     ArblibJL(; check_analytic=false, take_prec=false, warn_on_no_convergence=false, opts=C_NULL)
 
@@ -221,8 +225,10 @@ struct ArblibJL{O} <: AbstractIntegralExtensionAlgorithm
     warn_on_no_convergence::Bool
     opts::O
 end
-function ArblibJL(; check_analytic=false, take_prec=false, warn_on_no_convergence=false, opts=C_NULL)
-    isnothing(Base.get_extension(@__MODULE__, :IntegralsArblibExt)) && error("ArblibJL requires `using Arblib`")
+function ArblibJL(; check_analytic = false, take_prec = false,
+        warn_on_no_convergence = false, opts = C_NULL)
+    isnothing(Base.get_extension(@__MODULE__, :IntegralsArblibExt)) &&
+        error("ArblibJL requires `using Arblib`")
     return ArblibJL(check_analytic, take_prec, warn_on_no_convergence, opts)
 end
 
@@ -232,14 +238,15 @@ end
 Markov-chain based Vegas algorithm from MCIntegration.jl
 
 Refer to
-[`MCIntegration.integrate`](https://numericaleft.github.io/MCIntegration.jl/dev/lib/montecarlo/#MCIntegration.integrate-Tuple{Function})
+[`MCIntegration.integrate`](https://numericaleft.github.io/MCIntegration.jl/dev/lib/montecarlo/#MCIntegration.integrate-Tuple%7BFunction%7D)
 for documentation on the keywords, which are passed directly to the solver with a set of
 defaults that works for conforming integrands.
 """
-struct VEGASMC{K<:NamedTuple} <: AbstractIntegralExtensionAlgorithm
+struct VEGASMC{K <: NamedTuple} <: AbstractIntegralExtensionAlgorithm
     kws::K
 end
 function VEGASMC(; kws...)
-    isnothing(Base.get_extension(@__MODULE__, :IntegralsMCIntegrationExt)) && error("VEGASMC requires `using MCIntegration`")
+    isnothing(Base.get_extension(@__MODULE__, :IntegralsMCIntegrationExt)) &&
+        error("VEGASMC requires `using MCIntegration`")
     return VEGASMC(NamedTuple(kws))
 end
diff --git a/src/algorithms_sampled.jl b/src/algorithms_sampled.jl
index 6cd5dd5..9e71237 100644
--- a/src/algorithms_sampled.jl
+++ b/src/algorithms_sampled.jl
@@ -5,7 +5,6 @@ abstract type AbstractSampledIntegralAlgorithm <: SciMLBase.AbstractIntegralAlgo
 
 Struct for evaluating an integral via the trapezoidal rule.
 
-
 Example with sampled data:
 
 ```
@@ -27,7 +26,6 @@ Struct for evaluating an integral via the Simpson's composite 1/3-3/8
 rule over `AbstractRange`s (evenly spaced points) and
 Simpson's composite 1/3 rule for non-equidistant grids.
 
-
 Example with equidistant data:
 
 ```
diff --git a/src/infinity_handling.jl b/src/infinity_handling.jl
index ef03f69..0cf5c54 100644
--- a/src/infinity_handling.jl
+++ b/src/infinity_handling.jl
@@ -37,7 +37,7 @@ function substitute_t(t::Number, lb::Number, ub::Number)
     end
 end
 function substitute_t(t::AbstractVector, lb::AbstractVector, ub::AbstractVector)
-    x = similar(t, typeof(one(eltype(t))*(first(lb)+first(ub))))
+    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])
@@ -58,7 +58,7 @@ function substitute_f(f::F, dt, t, p, lb, ub) where {F}
 end
 
 function substitute_t(t::AbstractVector, lb::Number, ub::Number)
-    x = similar(t, typeof(one(eltype(t))*(lb+ub)))
+    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)
@@ -66,10 +66,11 @@ function substitute_t(t::AbstractVector, lb::Number, ub::Number)
     return x, jac
 end
 function substitute_t(t::AbstractArray, lb::AbstractVector, ub::AbstractVector)
-    x = similar(t, typeof(one(eltype(t))*(first(lb)+first(ub))))
+    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)
+        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
@@ -77,13 +78,13 @@ 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)))
+    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])
+        for idt in CartesianIndices(axes(dt)[begin:(end - 1)])
             dt[idt, i] *= j
         end
     end
@@ -95,22 +96,31 @@ function transformation_if_inf(prob, ::Val{true})
     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)
+    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,
+            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,
+            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,
+            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,
+            IntegralFunction{false}(let f = f.f
+                    (t, p) -> substitute_f(f, t, p, lb, ub)
+                end,
                 f.integrand_prototype)
         end
     end
diff --git a/src/simpsons.jl b/src/simpsons.jl
index 7c45b0f..06b5fb7 100644
--- a/src/simpsons.jl
+++ b/src/simpsons.jl
@@ -49,11 +49,13 @@ end
     isodd(j) &&
         return (x[begin + j + 1] - x[begin + j - 1])^3 / (x[begin + j] - x[begin + j - 1]) /
                (x[begin + j + 1] - x[begin + j]) / 6
-    iseven(j) && 
+    iseven(j) &&
         return (x[begin + j] - x[begin + j - 2]) / 6 *
-           (2 - (x[begin + j - 1] - x[begin + j - 2]) / (x[begin + j] - x[begin + j - 1])) +
-           (x[begin + j + 2] - x[begin + j]) / 6 *
-           (2 - (x[begin + j + 2] - x[begin + j + 1]) / (x[begin + j + 1] - x[begin + j]))
+               (2 -
+                (x[begin + j - 1] - x[begin + j - 2]) / (x[begin + j] - x[begin + j - 1])) +
+               (x[begin + j + 2] - x[begin + j]) / 6 *
+               (2 -
+                (x[begin + j + 2] - x[begin + j + 1]) / (x[begin + j + 1] - x[begin + j]))
 end
 
 function find_weights(x::AbstractVector, ::SimpsonsRule)
diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl
index ffd9e81..ba22e8f 100644
--- a/test/derivative_tests.jl
+++ b/test/derivative_tests.jl
@@ -10,39 +10,39 @@ reltol = 1e-3
 abstol = 1e-3
 
 alg_req = Dict(
-    QuadratureRule(gausslegendre, n=50) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false,
+    QuadratureRule(gausslegendre, n = 50) => (
+        nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false,
         allows_iip = false),
-    GaussLegendre(n=50) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false,
+    GaussLegendre(n = 50) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false,
         allows_iip = false),
     QuadGKJL() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = 1,
         allows_iip = true),
     HCubatureJL() => (nout = Inf, allows_batch = false, min_dim = 1,
-        max_dim = Inf, allows_iip = true),
-    # VEGAS() => (nout = 1, allows_batch = true, min_dim = 2, max_dim = Inf,
-        # allows_iip = true),
-    CubatureJLh() => (nout = Inf, allows_batch = true, min_dim = 1,
+        max_dim = Inf, allows_iip = true), 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),
-)
+        max_dim = Inf, allows_iip = true))
+# VEGAS() => (nout = 1, allows_batch = true, min_dim = 2, 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 should have same shape as parameters, independent of dimensionality
-integrands =(
-    (x, p) -> map(q -> prod(y -> sin(y*q), x), p),
+integrands = (
+    (x, p) -> map(q -> prod(y -> sin(y * q), x), p),
 )
 
 # function to turn the output into a scalar / test different tangent types
 scalarize_solution = (
     sol -> sin(sum(sol)),
-    sol -> sin(sol[1]),
+    sol -> sin(sol[1])
 )
 
 # we will be able to use broadcasting for this after https://github.com/FluxML/Zygote.jl/pull/1488
@@ -58,7 +58,7 @@ function f_helper!(f, y, x, p)
 end
 
 # the Zygote implementation is inconsistent about 0-d so we hijack it
-struct Scalar{T<:Real} <: Real
+struct Scalar{T <: Real} <: Real
     x::T
 end
 Base.iterate(a::Scalar) = (a.x, nothing)
@@ -67,13 +67,13 @@ Base.IteratorSize(::Type{Scalar{T}}) where {T} = Base.HasShape{0}()
 Base.eltype(::Type{Scalar{T}}) where {T} = T
 Base.length(a::Scalar) = 1
 Base.size(::Scalar) = ()
-Base.:+(a::Scalar, b::Scalar) = Scalar(a.x+b.x)
-Base.:*(a::Number, b::Scalar) = a*b.x
-Base.:*(a::Scalar, b::Number) = a.x*b
-Base.:*(a::Scalar, b::Scalar) = Scalar(a.x*b.x)
+Base.:+(a::Scalar, b::Scalar) = Scalar(a.x + b.x)
+Base.:*(a::Number, b::Scalar) = a * b.x
+Base.:*(a::Scalar, b::Number) = a.x * b
+Base.:*(a::Scalar, b::Scalar) = Scalar(a.x * b.x)
 Base.zero(a::Scalar) = Scalar(zero(a.x))
 Base.map(f, a::Scalar) = map(f, a.x)
-(::Type{T})(a::Scalar) where {T<:Real} = T(a.x)
+(::Type{T})(a::Scalar) where {T <: Real} = T(a.x)
 struct ScalarAxes end # the implementation doesn't preserve singleton axes
 Base.axes(::Scalar) = ScalarAxes()
 Base.iterate(::ScalarAxes) = nothing
@@ -85,7 +85,7 @@ Base.reshape(A::AbstractArray, ::ScalarAxes) = Scalar(only(A))
 # p can't be either without both prs
 function batch_helper(f, x, p)
     t = f(zero(eltype(x)), zero(eltype(eltype(p))))
-    typeof(t).([f(y, q) for q in p, y in eachslice(x; dims=ndims(x))])
+    typeof(t).([f(y, q) for q in p, y in eachslice(x; dims = ndims(x))])
 end
 
 function batch_helper!(f, y, x, p)
@@ -101,7 +101,8 @@ do_tests = function (; f, scalarize, lb, ub, p, alg, abstol, reltol)
     end
     testf(lb, ub, p)
 
-    dlb1, dub1, dp1 = Zygote.gradient(testf, lb, ub, p isa Number && f isa BatchIntegralFunction ? Scalar(p) : p)
+    dlb1, dub1, dp1 = Zygote.gradient(
+        testf, lb, ub, p isa Number && f isa BatchIntegralFunction ? Scalar(p) : p)
 
     f_lb = lb -> testf(lb, ub, p)
     f_ub = ub -> testf(lb, ub, p)
@@ -113,16 +114,16 @@ do_tests = function (; f, scalarize, lb, ub, p, alg, abstol, reltol)
     dub2 = getproperty(FiniteDiff, Symbol(:finite_difference_, dub))(f_ub, ub)
 
     if lb isa Number
-        @test dlb1 ≈ dlb2 atol=abstol rtol=reltol
-        @test dub1 ≈ dub2 atol=abstol rtol=reltol
+        @test dlb1≈dlb2 atol=abstol rtol=reltol
+        @test dub1≈dub2 atol=abstol rtol=reltol
     else # TODO: implement multivariate limit derivatives in ZygoteExt
-        @test_broken dlb1 ≈ dlb2 atol=abstol rtol=reltol
-        @test_broken dub1 ≈ dub2 atol=abstol rtol=reltol
+        @test_broken dlb1≈dlb2 atol=abstol rtol=reltol
+        @test_broken dub1≈dub2 atol=abstol rtol=reltol
     end
 
     # TODO: implement limit derivatives in ForwardDiffExt
-    @test_broken dlb2 ≈ getproperty(ForwardDiff, dlb)(dfdlb, lb) atol=abstol rtol=reltol
-    @test_broken dub2 ≈ getproperty(ForwardDiff, dub)(dfdub, ub) atol=abstol rtol=reltol
+    @test_broken dlb2≈getproperty(ForwardDiff, dlb)(dfdlb, lb) atol=abstol rtol=reltol
+    @test_broken dub2≈getproperty(ForwardDiff, dub)(dfdub, ub) atol=abstol rtol=reltol
 
     f_p = p -> testf(lb, ub, p)
 
@@ -131,15 +132,16 @@ do_tests = function (; f, scalarize, lb, ub, p, alg, abstol, reltol)
     dp2 = getproperty(FiniteDiff, Symbol(:finite_difference_, dp))(f_p, p)
     dp3 = getproperty(ForwardDiff, dp)(f_p, p)
 
-    @test dp1 ≈ dp2 atol=abstol rtol=reltol
-    @test dp2 ≈ dp3 atol=abstol rtol=reltol
+    @test dp1≈dp2 atol=abstol rtol=reltol
+    @test dp2≈dp3 atol=abstol rtol=reltol
 
     return
 end
 
-
 ### One Dimensional
-for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution)
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands),
+    (i, scalarize) in enumerate(scalarize_solution)
+
     req.nout > 1 || continue
     req.min_dim <= 1 || continue
 
@@ -148,16 +150,21 @@ for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize
 end
 
 ## One-dimensional nout
-for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands),
+    (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test
+
     req.nout > 1 || continue
     req.min_dim <= 1 || continue
 
     @info "One-dimensional, multivariate, oop derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i nout
-    do_tests(; f, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol)
+    do_tests(;
+        f, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol)
 end
 
 ### N-dimensional
-for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands),
+    (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test
+
     req.nout > 1 || continue
     req.min_dim <= dim <= req.max_dim || continue
 
@@ -166,104 +173,134 @@ for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize
 end
 
 ### N-dimensional nout
-for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands),
+    (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test,
+    nout in 1:max_nout_test
+
     req.nout > 1 || continue
     req.min_dim <= dim <= req.max_dim || continue
 
     @info "Multi-dimensional, multivariate, oop derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i dim nout
-    do_tests(; f, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol)
+    do_tests(; f, scalarize, lb = ones(dim), ub = 3ones(dim),
+        p = [2.0i for i in 1:nout], alg, abstol, reltol)
 end
 
 ### One Dimensional
-for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution)
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands),
+    (i, scalarize) in enumerate(scalarize_solution)
+
     req.allows_iip || continue
     req.nout > 1 || continue
     req.min_dim <= 1 || continue
 
     @info "One-dimensional, scalar, iip derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i
     fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(1))
-    do_tests(; f=fiip, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol)
+    do_tests(; f = fiip, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol)
 end
 
 ## One-dimensional nout
-for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands),
+    (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test
+
     req.allows_iip || continue
     req.nout > 1 || continue
     req.min_dim <= 1 || continue
 
     @info "One-dimensional, multivariate, iip derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i nout
     fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(nout))
-    do_tests(; f=fiip, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol)
+    do_tests(; f = fiip, scalarize, lb = 1.0, ub = 3.0,
+        p = [2.0i for i in 1:nout], alg, abstol, reltol)
 end
 
 ### N-dimensional
-for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands),
+    (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test
+
     req.allows_iip || continue
     req.nout > 1 || continue
     req.min_dim <= dim <= req.max_dim || continue
 
     @info "Multi-dimensional, scalar, iip derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i dim
     fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(1))
-    do_tests(; f=fiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol)
+    do_tests(;
+        f = fiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol)
 end
 
 ### N-dimensional nout iip
-for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands),
+    (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test,
+    nout in 1:max_nout_test
+
     req.allows_iip || continue
     req.nout > 1 || continue
     req.min_dim <= dim <= req.max_dim || continue
 
     @info "Multi-dimensional, multivariate, iip derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i dim nout
     fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(nout))
-    do_tests(; f=fiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol)
+    do_tests(; f = fiip, scalarize, lb = ones(dim), ub = 3ones(dim),
+        p = [2.0i for i in 1:nout], alg, abstol, reltol)
 end
 
 ### Batch, One Dimensional
-for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution)
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands),
+    (i, scalarize) in enumerate(scalarize_solution)
+
     req.allows_batch || continue
     req.nout > 1 || continue
     req.min_dim <= 1 || continue
 
     @info "Batched, one-dimensional, scalar, oop derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i
     bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p))
-    do_tests(; f=bf, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol)
+    do_tests(; f = bf, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol)
 end
 
 ## Batch, One-dimensional nout
-for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands),
+    (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test
+
     req.allows_batch || continue
     req.nout > 1 || continue
     req.min_dim <= 1 || continue
 
     @info "Batched, one-dimensional, multivariate, oop derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i nout
     bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p))
-    do_tests(; f=bf, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol)
+    do_tests(; f = bf, scalarize, lb = 1.0, ub = 3.0,
+        p = [2.0i for i in 1:nout], alg, abstol, reltol)
 end
 
 ### Batch, N-dimensional
-for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands),
+    (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test
+
     req.allows_batch || continue
     req.nout > 1 || continue
     req.min_dim <= dim <= req.max_dim || continue
 
     @info "Batched, multi-dimensional, scalar, oop derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i dim
     bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p))
-    do_tests(; f=bf, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol)
+    do_tests(;
+        f = bf, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol)
 end
 
 ### Batch, N-dimensional nout
-for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands),
+    (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test,
+    nout in 1:max_nout_test
+
     req.allows_batch || continue
     req.nout > 1 || continue
     req.min_dim <= dim <= req.max_dim || continue
 
     @info "Batch, multi-dimensional, multivariate, oop derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i dim nout
     bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p))
-    do_tests(; f=bf, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol)
+    do_tests(; f = bf, scalarize, lb = ones(dim), ub = 3ones(dim),
+        p = [2.0i for i in 1:nout], alg, abstol, reltol)
 end
 
 ### Batch, one-dimensional
-for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution)
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands),
+    (i, scalarize) in enumerate(scalarize_solution)
+
     req.allows_batch || continue
     req.allows_iip || continue
     req.nout > 1 || continue
@@ -271,11 +308,13 @@ for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize
 
     @info "Batched, one-dimensional, scalar, iip derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i
     bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(0))
-    do_tests(; f=bfiip, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol)
+    do_tests(; f = bfiip, scalarize, lb = 1.0, ub = 3.0, p = 2.0, alg, abstol, reltol)
 end
 
 ## Batch, one-dimensional nout
-for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands),
+    (i, scalarize) in enumerate(scalarize_solution), nout in 1:max_nout_test
+
     req.allows_batch || continue
     req.allows_iip || continue
     req.nout > 1 || continue
@@ -283,11 +322,14 @@ for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize
 
     @info "Batched, one-dimensional, multivariate, iip derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i nout
     bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(nout, 0))
-    do_tests(; f=bfiip, scalarize, lb = 1.0, ub = 3.0, p = [2.0i for i in 1:nout], alg, abstol, reltol)
+    do_tests(; f = bfiip, scalarize, lb = 1.0, ub = 3.0,
+        p = [2.0i for i in 1:nout], alg, abstol, reltol)
 end
 
 ### Batch, N-dimensional
-for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands),
+    (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test
+
     req.allows_batch || continue
     req.allows_iip || continue
     req.nout > 1 || continue
@@ -295,11 +337,15 @@ for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize
 
     @info "Batched, multi-dimensional, scalar, iip derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i dim
     bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(0))
-    do_tests(; f=bfiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = 2.0, alg, abstol, reltol)
+    do_tests(; f = bfiip, scalarize, lb = ones(dim),
+        ub = 3ones(dim), p = 2.0, alg, abstol, reltol)
 end
 
 ### Batch, N-dimensional nout iip
-for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test, nout in 1:max_nout_test
+for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands),
+    (i, scalarize) in enumerate(scalarize_solution), dim in 1:max_dim_test,
+    nout in 1:max_nout_test
+
     req.allows_batch || continue
     req.allows_iip || continue
     req.nout > 1 || continue
@@ -307,5 +353,6 @@ for (alg, req) in pairs(alg_req), (j, f) in enumerate(integrands), (i, scalarize
 
     @info "Batched, multi-dimensional, multivariate, iip derivative test" alg=nameof(typeof(alg)) integrand=j scalarize=i dim nout
     bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(nout, 0))
-    do_tests(; f=bfiip, scalarize, lb = ones(dim), ub = 3ones(dim), p = [2.0i for i in 1:nout], alg, abstol, reltol)
+    do_tests(; f = bfiip, scalarize, lb = ones(dim), ub = 3ones(dim),
+        p = [2.0i for i in 1:nout], alg, abstol, reltol)
 end
diff --git a/test/inf_integral_tests.jl b/test/inf_integral_tests.jl
index c93ea51..9164fd4 100644
--- a/test/inf_integral_tests.jl
+++ b/test/inf_integral_tests.jl
@@ -5,116 +5,116 @@ abstol = 1e-3
 
 # not all quadratures are compatible with infinities if they evaluate the endpoints
 alg_req = Dict(
-    QuadratureRule(gausslegendre, n=50) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false,
-        allows_iip = false, allows_inf=true),
-    # 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 = true, allows_inf=true),
+    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 = true, allows_inf = true),
     HCubatureJL() => (nout = Inf, allows_batch = false, min_dim = 1,
-        max_dim = Inf, allows_iip = true, allows_inf=true),
-    # VEGAS() => (nout = 1, allows_batch = true, min_dim = 2, max_dim = Inf,
-        # allows_iip = true),
-    CubatureJLh() => (nout = Inf, allows_batch = true, min_dim = 1,
-        max_dim = Inf, allows_iip = true, allows_inf=true),
-    # CubatureJLp() => (nout = Inf, allows_batch = true, min_dim = 1,
-    #     max_dim = Inf, allows_iip = true, allows_inf=false),
-    # 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),
+        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,
+#     allows_iip = false, allows_inf=true),
+# VEGAS() => (nout = 1, allows_batch = true, min_dim = 2, max_dim = Inf,
+# allows_iip = true),
+# CubatureJLp() => (nout = Inf, allows_batch = true, min_dim = 1,
+#     max_dim = Inf, allows_iip = true, allows_inf=false),
+# 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),
 problems = (
     (; # 1. multi-variate infinite limits: Gaussian
-        f=(x, p) -> pdf(MvNormal([0.00, 0.00], [0.4 0.0; 0.00 0.4]), x),
-        domain=(@SVector[-Inf, -Inf], @SVector[Inf, Inf]),
-        solution=1.00,
+        f = (x, p) -> pdf(MvNormal([0.00, 0.00], [0.4 0.0; 0.00 0.4]), x),
+        domain = (@SVector[-Inf, -Inf], @SVector[Inf, Inf]),
+        solution = 1.00
     ),
     (; # 2. multi-variate flipped infinite limits: Gaussian
-        f=(x, p) -> pdf(MvNormal([0.00, 0.00], [0.4 0.0; 0.00 0.4]), x),
-        domain=(@SVector[Inf, Inf], @SVector[-Inf, -Inf]),
-        solution=1.00,
+        f = (x, p) -> pdf(MvNormal([0.00, 0.00], [0.4 0.0; 0.00 0.4]), x),
+        domain = (@SVector[Inf, Inf], @SVector[-Inf, -Inf]),
+        solution = 1.00
     ),
     (; # 3. multi-variate mixed infinite/semi-infinite upper limit: Gaussian
-        f=(x, p) -> pdf(MvNormal([0.00, 0.00], [0.4 0.0; 0.00 0.4]), x),
-        domain=(@SVector[-Inf, 0], @SVector[Inf, Inf]),
-        solution=0.5,
+        f = (x, p) -> pdf(MvNormal([0.00, 0.00], [0.4 0.0; 0.00 0.4]), x),
+        domain = (@SVector[-Inf, 0], @SVector[Inf, Inf]),
+        solution = 0.5
     ),
     (; # 4. multi-variate mixed infinite/semi-infinite lower limit: Gaussian
-        f=(x, p) -> pdf(MvNormal([0.00, 0.00], [0.4 0.0; 0.00 0.4]), x),
-        domain=(@SVector[-Inf, -Inf], @SVector[Inf, 0]),
-        solution=0.5,
+        f = (x, p) -> pdf(MvNormal([0.00, 0.00], [0.4 0.0; 0.00 0.4]), x),
+        domain = (@SVector[-Inf, -Inf], @SVector[Inf, 0]),
+        solution = 0.5
     ),
     (; # 5. single-variable infinite limit: Gaussian
-        f=(x, p) -> pdf(Normal(0.00, 1.00), x),
-        domain=(-Inf, Inf),
-        solution=1.0,
+        f = (x, p) -> pdf(Normal(0.00, 1.00), x),
+        domain = (-Inf, Inf),
+        solution = 1.0
     ),
     (; # 6. single-variable flipped infinite limit: Gaussian
-        f=(x, p) -> pdf(Normal(0.00, 1.00), x),
-        domain=(Inf, -Inf),
-        solution=-1.0,
+        f = (x, p) -> pdf(Normal(0.00, 1.00), x),
+        domain = (Inf, -Inf),
+        solution = -1.0
     ),
     (; # 7. single-variable semi-infinite upper limit: Gaussian
-        f=(x, p) -> pdf(Normal(0.00, 1.00), x),
-        domain=(0.00, Inf),
-        solution=0.5,
+        f = (x, p) -> pdf(Normal(0.00, 1.00), x),
+        domain = (0.00, Inf),
+        solution = 0.5
     ),
     (; # 8. single-variable flipped, semi-infinite upper limit: Gaussian
-        f=(x, p) -> pdf(Normal(0.00, 1.00), x),
-        domain=(0.00, -Inf),
-        solution=-0.5,
+        f = (x, p) -> pdf(Normal(0.00, 1.00), x),
+        domain = (0.00, -Inf),
+        solution = -0.5
     ),
     (; # 9. single-variable semi-infinite lower limit: Gaussian
-        f=(x, p) -> pdf(Normal(0.00, 1.00), x),
-        domain=(-Inf, 0.00),
-        solution=0.5,
+        f = (x, p) -> pdf(Normal(0.00, 1.00), x),
+        domain = (-Inf, 0.00),
+        solution = 0.5
     ),
     (; # 10. single-variable flipped, semi-infinite lower limit: Gaussian
-        f=(x, p) -> pdf(Normal(0.00, 1.00), x),
-        domain=(Inf, 0.00),
-        solution=-0.5,
+        f = (x, p) -> pdf(Normal(0.00, 1.00), x),
+        domain = (Inf, 0.00),
+        solution = -0.5
     ),
     (; # 11. single-variable infinite limit: Lorentzian
-        f=(x, p) -> 1 / (x^2 + 1),
-        domain=(-Inf, Inf),
-        solution=pi/1,
+        f = (x, p) -> 1 / (x^2 + 1),
+        domain = (-Inf, Inf),
+        solution = pi / 1
     ),
     (; # 12. single-variable shifted, semi-infinite lower limit: Lorentzian
-        f=(x, p) -> 1 / ((x-2)^2 + 1),
-        domain=(-Inf, 2),
-        solution=pi/2,
+        f = (x, p) -> 1 / ((x - 2)^2 + 1),
+        domain = (-Inf, 2),
+        solution = pi / 2
     ),
     (; # 13. single-variable shifted, semi-infinite upper limit: Lorentzian
-        f=(x, p) -> 1 / ((x-2)^2 + 1),
-        domain=(2, Inf),
-        solution=pi/2,
+        f = (x, p) -> 1 / ((x - 2)^2 + 1),
+        domain = (2, Inf),
+        solution = pi / 2
     ),
     (; # 14. single-variable flipped, shifted, semi-infinite lower limit: Lorentzian
-        f=(x, p) -> 1 / ((x-2)^2 + 1),
-        domain=(Inf, 2),
-        solution=-pi/2,
+        f = (x, p) -> 1 / ((x - 2)^2 + 1),
+        domain = (Inf, 2),
+        solution = -pi / 2
     ),
     (; # 15. single-variable flipped, shifted, semi-infinite upper limit: Lorentzian
-        f=(x, p) -> 1 / ((x-2)^2 + 1),
-        domain=(2, -Inf),
-        solution=-pi/2,
+        f = (x, p) -> 1 / ((x - 2)^2 + 1),
+        domain = (2, -Inf),
+        solution = -pi / 2
     ),
     (; # 16. single-variable finite limits: constant
-        f=(x, p) -> 1.0,
-        domain=(1, 3),
-        solution=2,
+        f = (x, p) -> 1.0,
+        domain = (1, 3),
+        solution = 2
     ),
     (; # 17. single-variable flipped, finite limits: constant
-        f=(x, p) -> 1.0,
-        domain=(3, 1),
-        solution=-2,
-    ),
+        f = (x, p) -> 1.0,
+        domain = (3, 1),
+        solution = -2
+    )
 )
 
 function f_helper!(f, y, x, p)
@@ -123,7 +123,7 @@ function f_helper!(f, y, x, p)
 end
 
 function batch_helper(f, x, p)
-    map(i -> f(x[axes(x)[begin:end-1]..., i], p), axes(x)[end])
+    map(i -> f(x[axes(x)[begin:(end - 1)]..., i], p), axes(x)[end])
 end
 
 function batch_helper!(f, y, x, p)
@@ -134,7 +134,7 @@ 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))
-    @test abs(only(sol) - solution) < max(abstol, reltol*abs(solution))
+    @test abs(only(sol) - solution) < max(abstol, reltol * abs(solution))
     cache = @test_nowarn @inferred init(prob, alg; do_inf_transformation = Val(true))
     @test_nowarn @inferred solve!(cache)
     @test_nowarn @inferred solve(prob, alg; do_inf_transformation = Val(true))
@@ -159,7 +159,7 @@ for (alg, req) in pairs(alg_req), (j, (; f, domain, solution)) in enumerate(prob
 
     @info "iip infinity test" alg=nameof(typeof(alg)) problem=j
     fiip = IntegralFunction((y, x, p) -> f_helper!(f, y, x, p), zeros(size(solution)))
-    do_tests(; f=fiip, domain, solution, alg, abstol, reltol)
+    do_tests(; f = fiip, domain, solution, alg, abstol, reltol)
 end
 
 # BatchIntegralFunction{false}
@@ -171,7 +171,7 @@ for (alg, req) in pairs(alg_req), (j, (; f, domain, solution)) in enumerate(prob
 
     @info "Batched, oop infinity test" alg=nameof(typeof(alg)) problem=j
     bf = BatchIntegralFunction((x, p) -> batch_helper(f, x, p))
-    do_tests(; f=bf, domain, solution, alg, abstol, reltol)
+    do_tests(; f = bf, domain, solution, alg, abstol, reltol)
 end
 
 # BatchIntegralFunction{true}
@@ -184,5 +184,5 @@ for (alg, req) in pairs(alg_req), (j, (; f, domain, solution)) in enumerate(prob
 
     @info "Batched, iip infinity test" alg=nameof(typeof(alg)) problem=j
     bfiip = BatchIntegralFunction((y, x, p) -> batch_helper!(f, y, x, p), zeros(0))
-    do_tests(; f=bfiip, domain, solution, alg, abstol, reltol)
+    do_tests(; f = bfiip, domain, solution, alg, abstol, reltol)
 end
diff --git a/test/interface_tests.jl b/test/interface_tests.jl
index 4eea513..272aed5 100644
--- a/test/interface_tests.jl
+++ b/test/interface_tests.jl
@@ -10,7 +10,8 @@ abstol = 1e-3
 alg_req = Dict(
     QuadGKJL() => (nout = Inf, allows_batch = true, min_dim = 1, max_dim = 1,
         allows_iip = true),
-    QuadratureRule(gausslegendre, n=50) => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false,
+    QuadratureRule(gausslegendre, n = 50) => (
+        nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false,
         allows_iip = false),
     GaussLegendre() => (nout = Inf, min_dim = 1, max_dim = 1, allows_batch = false,
         allows_iip = false),
@@ -18,7 +19,7 @@ alg_req = Dict(
         max_dim = Inf, allows_iip = true),
     # VEGAS() => (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,
+    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,
         max_dim = Inf, allows_iip = true),
@@ -32,27 +33,28 @@ alg_req = Dict(
         max_dim = Inf, allows_iip = true),
     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),
+    ArblibJL() => (
+        nout = 1, allows_batch = false, min_dim = 1, max_dim = 1, allows_iip = true)
 )
 
 integrands = [
     (x, p) -> 1.0,
-    (x, p) -> x isa Number ? cos(x) : prod(cos.(x)),
+    (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)]
+                (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)),
+    (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),
+    (ndim, nout, lb, ub) -> exact_sol[2](ndim, nout, lb, ub) * collect(1:nout)
 ]
 
 batch_f(f) = (pts, p) -> begin
@@ -67,7 +69,7 @@ end
 batch_iip_f(f) = (fevals, pts, p) -> begin
     ax = axes(pts)
     for i in ax[end]
-        x = pts[ax[begin:(end-1)]..., i]
+        x = pts[ax[begin:(end - 1)]..., i]
         fevals[i] = f(x, p)
     end
     nothing
@@ -270,10 +272,11 @@ 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(
+                        (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 ||
-                        !req.allows_iip
+                       !req.allows_iip
                         continue
                     end
                     @info "Alg = $(nameof(typeof(alg))), Integrand = $i, Dimension = $dim, Output Dimension = $nout"
@@ -354,7 +357,8 @@ end
 @testset "Allowed keyword test" begin
     f(u, p) = sum(sin.(u))
     prob = IntegralProblem(f, ones(3), 3ones(3))
-    @test_throws Integrals.CommonKwargError((:relztol => 1e-3, :abstol => 1e-3)) solve(prob,
+    @test_throws Integrals.CommonKwargError((:relztol => 1e-3, :abstol => 1e-3)) solve(
+        prob,
         HCubatureJL();
         relztol = 1e-3,
         abstol = 1e-3)
diff --git a/test/nested_ad_tests.jl b/test/nested_ad_tests.jl
index 079d8fa..d4506aa 100644
--- a/test/nested_ad_tests.jl
+++ b/test/nested_ad_tests.jl
@@ -11,7 +11,7 @@ my_solution = my_integration(my_parameters)
 @test ForwardDiff.jacobian(my_integration, my_parameters) == [0.0 8.0]
 @test ForwardDiff.jacobian(x -> ForwardDiff.jacobian(my_integration, x), my_parameters) ==
       [0.0 0.0
-    0.0 4.0]
+       0.0 4.0]
 
 ff(x, p) = sum(sin.(x .* p))
 lb = ones(2)
diff --git a/test/runtests.jl b/test/runtests.jl
index 8d15e5b..090b8fd 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,7 +2,6 @@ using Pkg
 using SafeTestsets
 using Test
 
-
 @time @safetestset "Quality Assurance" include("qa.jl")
 @time @safetestset "Interface Tests" include("interface_tests.jl")
 @time @safetestset "Derivative Tests" include("derivative_tests.jl")
diff --git a/test/sampled_tests.jl b/test/sampled_tests.jl
index 19a6787..9c7398b 100644
--- a/test/sampled_tests.jl
+++ b/test/sampled_tests.jl
@@ -8,7 +8,7 @@ using Integrals, Test
     grid2 = rand(npoints) .* (ub - lb) .+ lb
     grid2 = [lb; sort(grid2); ub]
 
-    grid3 = rand(npoints+1).*(ub-lb) .+ lb # also test odd number of points
+    grid3 = rand(npoints + 1) .* (ub - lb) .+ lb # also test odd number of points
     grid3 = [lb; sort(grid3); ub]
 
     exact_sols = [1 / 6 * (ub^6 - lb^6), sin(ub) - sin(lb)]