When performing a GSA on the Serial function runs quicker than the parallel counterpart. I have included the same test with an ODE version to highlight how the multithreaded version should execute more efficiently.
using DifferentialEquations
# Define DDE model
begin
function bc_model(du, u, h, p, t)
p0, q0, v0, d0, p1, q1, v1, d1, d2, beta0, beta1, tau = p
hist3 = h(p, t - tau)[3]
du[1] = (v0 / (1 + beta0 * (hist3^2))) * (p0 - q0) * u[1] - d0 * u[1]
du[2] = (v0 / (1 + beta0 * (hist3^2))) * (1 - p0 + q0) * u[1] +
(v1 / (1 + beta1 * (hist3^2))) * (p1 - q1) * u[2] - d1 * u[2]
du[3] = (v1 / (1 + beta1 * (hist3^2))) * (1 - p1 + q1) * u[2] - d2 * u[3]
end
h(p, t) = ones(3)
tau = 1
lags = [tau]
p0 = 0.2;
q0 = 0.3;
v0 = 1;
d0 = 5;
p1 = 0.2;
q1 = 0.3;
v1 = 1;
d1 = 1;
d2 = 1;
beta0 = 1;
beta1 = 1;
p = [p0, q0, v0, d0, p1, q1, v1, d1, d2, beta0, beta1, tau]
tspan = (0.0, 10.0)
u0 = [1.0, 1.0, 1.0]
prob = DDEProblem(bc_model, u0, h, tspan, p; constant_lags = lags)
alg = MethodOfSteps(Tsit5())
@time sol = solve(prob, alg)
end
using GlobalSensitivity, QuasiMonteCarlo, Statistics
samples = 30000
lb = 0.9*prob.p
ub = 1.1*prob.p
sampler = SobolSample()
A,B = QuasiMonteCarlo.generate_design_matrices(samples,lb,ub,sampler)
# Serial
f1_s = function (p)
prob1 = remake(prob;p=p)
sol = solve(prob1,MethodOfSteps(Tsit5()))
[mean(sol[1,:]), mean(sol[2,:])]
end
@time res = gsa(f1_s,Sobol(),A,B)
# Roughly 40s
## Para
f1_p = function (p)
prob_func(prob,i,repeat) = remake(prob;p=p[:,i], constant_lags = lags)
ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
sol = solve(ensemble_prob,MethodOfSteps(Tsit5()),EnsembleThreads();trajectories=size(p,2))
# Now sol[i] is the solution for the ith set of parameters
out = zeros(2,size(p,2))
for i in 1:size(p,2)
out[1,i] = mean(sol[i][1,:])
out[2,i] = mean(sol[i][2,:])
end
out
end
@time res = gsa(f1_p,Sobol(),A,B,batch = true)
# Roughly 170s
## ODE EXamples ##
function f(du,u,p,t)
du[1] = p[1]*u[1] - p[2]*u[1]*u[2] #prey
du[2] = -p[3]*u[2] + p[4]*u[1]*u[2] #predator
end
u0 = [1.0;1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob = ODEProblem(f,u0,tspan,p)
t = collect(range(0, stop=10, length=200))
using GlobalSensitivity, QuasiMonteCarlo, Statistics
samples = 50000
lb = 0.9*prob.p
ub = 1.1*prob.p
sampler = SobolSample()
A,B = QuasiMonteCarlo.generate_design_matrices(samples,lb,ub,sampler)
# Serial
f1_s = function (p)
prob1 = remake(prob;p=p)
sol = solve(prob1,Tsit5(), saveat = t)
[mean(sol[1,:]), mean(sol[2,:])]
end
@time res = gsa(f1_s,Sobol(),A,B)
# Roughly 35s
## Para
f1_p = function (p)
prob_func(prob,i,repeat) = remake(prob;p=p[:,i])
ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
sol = solve(ensemble_prob,Tsit5(), saveat = t,EnsembleThreads();trajectories=size(p,2))
# Now sol[i] is the solution for the ith set of parameters
out = zeros(2,size(p,2))
for i in 1:size(p,2)
out[1,i] = mean(sol[i][1,:])
out[2,i] = mean(sol[i][2,:])
end
out
end
@time res = gsa(f1_p,Sobol(),A,B,batch = true)
# Roughly 18.5s - Quicker as expected with multiple threads
Describe the bug 🐞
When performing a GSA on the Serial function runs quicker than the parallel counterpart. I have included the same test with an ODE version to highlight how the multithreaded version should execute more efficiently.
Expected behavior
The multithreaded version should run more efficiently than the serial version
Minimal Reproducible Example 👇
Environment (please complete the following information):
using Pkg; Pkg.status()using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)versioninfo()Additional context
This is on the example from the DDE website