Merge pull request #129 from ArnoStrouwen/inf

start generalizing inf

Author: ChrisRackauckas

src/Integrals.jl

@@ -56,7 +56,6 @@ function checkkwargs(kwargs...)
     end
     return nothing
 end
-
 """
 ```julia
 solve(prob::IntegralProblem, alg::SciMLBase.AbstractIntegralAlgorithm; kwargs...)

src/algorithms.jl

@@ -89,15 +89,15 @@ VEGAS(; nbins = 100, ncalls = 1000, debug = false) = VEGAS(nbins, ncalls, debug)
 Struct for evaluating an integral via (composite) Gauss-Legendre quadrature.
 The field `C` will be `true` if `subintervals > 1`, and `false` otherwise.
 
-The fields `nodes::N` and `weights::W` are defined by 
-`nodes, weights = gausslegendre(n)` for a given number of nodes `n`. 
+The fields `nodes::N` and `weights::W` are defined by
+`nodes, weights = gausslegendre(n)` for a given number of nodes `n`.
 
-The field `subintervals::Int64 = 1` (with default value `1`) defines the 
-number of intervals to partition the original interval of integration 
+The field `subintervals::Int64 = 1` (with default value `1`) defines the
+number of intervals to partition the original interval of integration
 `[a, b]` into, splitting it into `[xⱼ, xⱼ₊₁]` for `j = 1,…,subintervals`,
-where `xⱼ = a + (j-1)h` and `h = (b-a)/subintervals`. Gauss-Legendre 
-quadrature is then applied on each subinterval. For example, if 
-`[a, b] = [-1, 1]` and `subintervals = 2`, then Gauss-Legendre 
+where `xⱼ = a + (j-1)h` and `h = (b-a)/subintervals`. Gauss-Legendre
+quadrature is then applied on each subinterval. For example, if
+`[a, b] = [-1, 1]` and `subintervals = 2`, then Gauss-Legendre
 quadrature will be applied separately on `[-1, 0]` and `[0, 1]`,
 summing the two results.
 """

src/infinity_handling.jl

@@ -49,21 +49,46 @@ function substitute_f_vector(t, p, f, lb, ub)
     end
     f(x, p) * prod(jac_diag)
 end
+function substitute_f_vector_iip(dt, t, p, f, lb, ub)
+    x = similar(t)
+    jac_diag = similar(t)
+    for i in eachindex(lb)
+        if isinf(lb[i]) && isinf(ub[i])
+            x[i] = t[i] / (1 - t[i]^2)
+            jac_diag[i] = (1 + t[i]^2) / (1 - t[i]^2)^2
+        elseif isinf(lb[i])
+            x[i] = ub[i] + (t[i] / (1 + t[i]))
+            jac_diag[i] = 1 / ((1 + t[i])^2)
+        elseif isinf(ub[i])
+            x[i] = lb[i] + (t[i] / (1 - t[i]))
+            jac_diag[i] = 1 / ((1 - t[i])^2)
+        else
+            x[i] = t[i]
+            jac_diag[i] = one(lb[i])
+        end
+    end
+    f(dt, x, p)
+    dt .*= prod(jac_diag)
+end
+
 function transformation_if_inf(prob, ::Val{true})
     lb = prob.lb
     ub = prob.ub
     f = prob.f
     if lb isa Number
         lb_sub, ub_sub = substitute_bounds(lb, ub)
         f_sub = (t, p) -> substitute_f_scalar(t, p, f, lb, ub)
-        return remake(prob, f = f_sub, lb = lb_sub, ub = ub_sub)
     else
         bounds = substitute_bounds.(lb, ub)
         lb_sub = first.(bounds)
         ub_sub = last.(bounds)
-        f_sub = (t, p) -> substitute_f_vector(t, p, f, lb, ub)
-        return remake(prob, f = f_sub, lb = lb_sub, ub = ub_sub)
+        if isinplace(prob)
+            f_sub = (dt, t, p) -> substitute_f_vector_iip(dt, t, p, f, lb, ub)
+        else
+            f_sub = (t, p) -> substitute_f_vector(t, p, f, lb, ub)
+        end
     end
+    return remake(prob, f = f_sub, lb = lb_sub, ub = ub_sub)
 end
 
 function transformation_if_inf(prob, ::Nothing)

test/inf_integral_tests.jl

@@ -18,6 +18,14 @@ sol = solve(prob, HCubatureJL(), reltol = 1e-3, abstol = 1e-3)
 @test (0.500 - sol.u)^2 < 1e-6
 @test_nowarn @inferred Integrals.transformation_if_inf(prob, Val(true))
 
+μ = [0.00, 0.00]
+Σ = [0.4 0.0; 0.00 0.4]
+d = MvNormal(μ, Σ)
+m_iip(dx, x, p) = dx .= pdf(d, x)
+prob = IntegralProblem(m_iip, [-Inf, 0.00], [Inf, Inf])
+sol = solve(prob, HCubatureJL(), reltol = 1e-3, abstol = 1e-3)
+@test (0.500 - sol.u[1])^2 < 1e-6
+
 f(x, p) = pdf(Normal(0.00, 1.00), x)
 prob = IntegralProblem(f, -Inf, Inf)
 sol = solve(prob, HCubatureJL(), reltol = 1e-3, abstol = 1e-3)