/
IntegralsCuba.jl
125 lines (114 loc) · 4.3 KB
/
IntegralsCuba.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
module IntegralsCuba
using Integrals, Cuba
import Integrals: transformation_if_inf, scale_x, scale_x!
abstract type AbstractCubaAlgorithm <: SciMLBase.AbstractIntegralAlgorithm end
struct CubaVegas <: AbstractCubaAlgorithm end
struct CubaSUAVE <: AbstractCubaAlgorithm end
struct CubaDivonne <: AbstractCubaAlgorithm end
struct CubaCuhre <: AbstractCubaAlgorithm end
function Integrals.__solvebp_call(prob::IntegralProblem, alg::AbstractCubaAlgorithm,
sensealg,
lb, ub, p, args...;
reltol = 1e-8, abstol = 1e-8,
maxiters = alg isa CubaSUAVE ? 1000000 : typemax(Int),
kwargs...)
prob = transformation_if_inf(prob) #intercept for infinite transformation
p = p
if lb isa Number && prob.batch == 0
_x = Float64[lb]
elseif lb isa Number
_x = zeros(length(lb), prob.batch)
elseif prob.batch == 0
_x = zeros(length(lb))
else
_x = zeros(length(lb), prob.batch)
end
ub = ub
lb = lb
if prob.batch == 0
if isinplace(prob)
f = function (x, dx)
prob.f(dx, scale_x!(_x, ub, lb, x), p)
dx .*= prod((y) -> y[1] - y[2], zip(ub, lb))
end
else
f = function (x, dx)
dx .= prob.f(scale_x!(_x, ub, lb, x), p) .*
prod((y) -> y[1] - y[2], zip(ub, lb))
end
end
else
if lb isa Number
if isinplace(prob)
f = function (x, dx)
#todo check scale_x!
prob.f(dx', scale_x!(view(_x, 1:length(x)), ub, lb, x), p)
dx .*= prod((y) -> y[1] - y[2], zip(ub, lb))
end
else
if prob.f([lb ub], p) isa Vector
f = function (x, dx)
dx .= prob.f(scale_x(ub, lb, x), p)' .*
prod((y) -> y[1] - y[2], zip(ub, lb))
end
else
f = function (x, dx)
dx .= prob.f(scale_x(ub, lb, x), p) .*
prod((y) -> y[1] - y[2], zip(ub, lb))
end
end
end
else
if isinplace(prob)
f = function (x, dx)
prob.f(dx, scale_x(ub, lb, x), p)
dx .*= prod((y) -> y[1] - y[2], zip(ub, lb))
end
else
if prob.f([lb ub], p) isa Vector
f = function (x, dx)
dx .= prob.f(scale_x(ub, lb, x), p)' .*
prod((y) -> y[1] - y[2], zip(ub, lb))
end
else
f = function (x, dx)
dx .= prob.f(scale_x(ub, lb, x), p) .*
prod((y) -> y[1] - y[2], zip(ub, lb))
end
end
end
end
end
ndim = length(lb)
nvec = prob.batch == 0 ? 1 : prob.batch
if alg isa CubaVegas
out = Cuba.vegas(f, ndim, prob.nout; rtol = reltol,
atol = abstol, nvec = nvec,
maxevals = maxiters, kwargs...)
elseif alg isa CubaSUAVE
out = Cuba.suave(f, ndim, prob.nout; rtol = reltol,
atol = abstol, nvec = nvec,
maxevals = maxiters, kwargs...)
elseif alg isa CubaDivonne
out = Cuba.divonne(f, ndim, prob.nout; rtol = reltol,
atol = abstol, nvec = nvec,
maxevals = maxiters, kwargs...)
elseif alg isa CubaCuhre
out = Cuba.cuhre(f, ndim, prob.nout; rtol = reltol,
atol = abstol, nvec = nvec,
maxevals = maxiters, kwargs...)
end
if isinplace(prob) || prob.batch != 0
val = out.integral
else
if prob.nout == 1 && prob.f(lb, p) isa Number
val = out.integral[1]
else
val = out.integral
end
end
SciMLBase.build_solution(prob, alg, val, out.error,
chi = out.probability, retcode = :Success)
end
export CubaVegas, CubaSUAVE, CubaDivonne, CubaCuhre
end