I’m trying to solve a global optimization problem concerning the parameters of the “Two Compartment Model”. Given an initial parameter search space that’s split into smaller boxes of a predefined maximum width, I’m interested in those boxes that contain zero in their partial derivatives of the error function, in my case a simple sum of least squares. The hull of those boxes – which should be significantly smaller than the initial search space – should then contain the optimal parameters for the problem with given measurement data and a known exact solution for y_2 only.
I already implemented this problem in C++, but now I wanted to try it in Julia, because the search for suitable boxes (the step that should be executed on the GPU) depends on automatic differentiation, which is not that easy to do in CUDA when working with intervals, at least to my knowledge.
The combination of IntervalArithmetic.jl, ForwardDiff.jl and CUDA.jl looked very promising, and theoretically it should have been as easy as writing code that works on the CPU, then wrap my data in a CuArray and Julia handles the rest. But I guess I was mistaken, because I can’t get it to work. This is the Code:
#!/usr/bin/env julia
using ForwardDiff
using IntervalArithmetic
using CUDA
using CRlibm
setrounding(Interval, :accurate)
SPLIT_MAX = 0.05
TOL = 1e-10
# measurement data for y2
y2m = [0.0532, 0.0478, 0.0410, 0.0328,
0.0323, 0.0148, 0.0216, 0.0127,
0.0099, 0.0081, 0.0065, 0.0043,
0.0013, 0.0015, 0.0060, 0.0126]
# y2 of the ODE
function y2(t, p1, p2, p3)
D = sqrt((p1-p2+p3)^2 + 4*p2*p3)
a = p3 / D
ps = p1 + p2 + p3
l1 = 0.5 * (ps-D)
l2 = 0.5 * (ps+D)
return a * (exp(-l1*t)-exp(-l2*t))
end
# same as above, but with a box as argument
Ob(p) = sum((y2(t, p[1], p[2], p[3]) - y2m[t])^2 for t in 1:16)
# typecast from box to vector so that gradient() works
Vector(B::IntervalBox) = [X for X in B]
dO = p -> IntervalBox(ForwardDiff.gradient(Ob, Vector(p)))
# monotonicity test
monotest(p) = contains_zero(dO(p))
# bisect initial search space until boxes are no wider than SPLIT_MAX
function split(domain)
to_split = [domain]
splits = []
while !isempty(to_split)
current_box = pop!(to_split)
if diam(current_box) < SPLIT_MAX
push!(splits, current_box)
else
append!(to_split, bisect(current_box))
end
end
return splits
end
# search space
P1 = 0.01..2.0
P2 = 0.05..3.0
P3 = 0.05..3.0
P = IntervalBox(P1, P2, P3)
println("Search space:\n\t", P)
splits = split(P)
println("Searching in ", length(splits), " boxes ",
"with a maximum width of ", SPLIT_MAX, " ...")
# find boxes that contain possible optimum
# XXX on the CPU, this works ...
# mask = monotest.(splits)
# ... but this does not work on the GPU
d_splits = CuArray(splits)
mask = monotest.(d_splits)
# only keep those boxes
interesting = splits[mask]
# get the error for every suitable box
errors = Ob.(interesting)
# get boxes with the smallest acceptable error
best_idx = findall(x->abs(x.lo*mid(x))<=TOL, errors)
best_boxes = [(errors[i], interesting[i]) for i in best_idx]
println("Best boxes:")
for best in best_boxes
println(best[1], best[2])
end
println("Hull:\n\t", reduce(∪, [b[2] for b in best_boxes]))
… which results in the following error:
ERROR: LoadError: CuArray only supports element types that are stored inline
Stacktrace:
[1] CuArray{Any, 1, CUDA.Mem.DeviceBuffer}(#unused#::UndefInitializer, dims::Tuple{Int64})
@ CUDA ~/.julia/packages/CUDA/0IDh2/src/array.jl:36
[2] CuArray
@ ~/.julia/packages/CUDA/0IDh2/src/array.jl:287 [inlined]
[3] CuArray
@ ~/.julia/packages/CUDA/0IDh2/src/array.jl:292 [inlined]
[4] CuArray(A::Base.Vector{Any})
@ CUDA ~/.julia/packages/CUDA/0IDh2/src/array.jl:301
[5] top-level scope
@ ~/twocomp.jl:79
in expression starting at ~/twocomp.jl:79
It appears that CUDA doesn’t like the list of IntervalBoxes. I already tried splitting my one list of 3d boxes into three lists of 1d boxes, changing the monotest.(d_splits) into broadcast(monotest, d_s1, d_s2, d_s3) (where d_sn are my lists of boxes), but still no success. Instead of the CuArray error I get errors about type instability.
What am I missing here?
My Julia version is 1.6.5 (that’s what Fedora 34 Server Edition provides atm) with CUDA v3.6.2, IntervalArithmetic v0.20.2 and ForwardDiff v0.10.24.