# Parallelize GSA on continuous outputs

Hi there,

I am just wondering how I could adapt the GSA example in the docs such that it performs Sobol analysis on the continuous outputs instead of the discrete outputs. sol[1,:] & sol[2,:]

I have tried multiple different things however nothing appears to be going in my favor.

I am sure it’s a simple adaptation so I apologise and thank you in advance

Thanks,
Harry

``````t = collect(range(0, stop=10, length=200))

f1 = function (p)
prob_func(prob,i,repeat) = remake(prob;p=p[:,i])
ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
# 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] = maximum(sol[i][2,:])
end
out
end
``````

What do you mean by continuous outputs here? This looks like just discrete

Yes the example is discrete, however I am looking for SA on the time series outputs of sol[1,:] and sol[2,:]. I am just wondering how I could adapt the above example to incorporate this.

I am looking for something similar to what is discussed in the above link. The only difference being I use GSA and it’s run in parallel.

You mean, just change out to match the time series size and don’t take the mean? That should just work, did you try it?

Yes I tried changing it however I feel I am getting my sizes mixed up somewhere

``````  f1 = function (p)
prob_func(prob,i,repeat) = remake(prob;p=p[:,i])
ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
# Now sol[i] is the solution for the ith set of parameters
out = zeros(2,size(p,2),200)
for i in 1:size(p,2)
out[1,i,:] = sol[i][1,:]
out[2,i,:] = sol[i][2,:]
end
out
end
``````

Thanks
Harry

What’s the issue?

So the below returns just a 4 element vector containing NANs, I have included the corresponsing serial run which retrurns a 400x4 vector i.e giving me a sensitivity at each time point.

``````using DifferentialEquations, Plots, LinearAlgebra, Distributions, Random, GlobalSensitivity, QuasiMonteCarlo
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))

f1 = function (p) # Para
prob_func(prob,i,repeat) = remake(prob;p=p[:,i])
ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
# Now sol[i] is the solution for the ith set of parameters
out = zeros(2,size(p,2),200)
for i in 1:size(p,2)
out[1,i,:] = sol[i][1,:]
out[2,i,:] = sol[i][2,:]
end
out
end
## Serial function
f1 = function (p) #Serial
prob1 = remake(prob;p=p)
sol = solve(prob1,Tsit5();saveat=t)
[sol[1,:]; sol[2,:]]
end

samples = 50
lb = [1.0, 1.0, 1.0, 1.0]
ub = [5.0, 5.0, 5.0, 5.0]
sampler = SobolSample()
A,B = QuasiMonteCarlo.generate_design_matrices(samples,lb,ub,sampler)

sobol_result = gsa(f1,Sobol(),A,B, batch = true)
``````

I have been playing around with a couple of things,

If I write

``````  f1_batch = function (p) # Para
prob_func(prob,i,repeat) = remake(prob;p=p[:,i])
ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
# Now sol[i] is the solution for the ith set of parameters
# create the output array:
out = [[sol[i][1,:]; sol[i][2,:]] for i in 1:size(p,2)]
end
``````

This one generates an error with

``````ERROR: MethodError: no method matching *(::Vector{Float64}, ::Vector{Float64})

Some of the types have been truncated in the stacktrace for improved reading. To emit complete information
in the stack trace, evaluate `TruncatedStacktraces.VERBOSE[] = true` and re-run the code.

Closest candidates are:
*(::Any, ::Any, ::Any, ::Any...)
@ Base operators.jl:578
*(::StridedMatrix{T}, ::StridedVector{S}) where {T<:Union{Float32, Float64, ComplexF64, ComplexF32}, S<:Real}
@ LinearAlgebra /opt/julia-1.9.0-rc1/share/julia/stdlib/v1.9/LinearAlgebra/src/matmul.jl:49
@ LinearAlgebra /opt/julia-1.9.0-rc1/share/julia/stdlib/v1.9/LinearAlgebra/src/lq.jl:269
...

Stacktrace:
[3] getindex
[4] copy
[5] materialize
[6] (::GlobalSensitivity.var"#36#57"{Vector{Vector{Vector{Float64}}}, Vector{Vector{Float64}}, Vector{Vector{Float64}}})(k::Int64)
@ GlobalSensitivity ./none:0
[7] iterate
@ ./generator.jl:47 [inlined]
[8] collect
@ ./array.jl:782 [inlined]
[9] gsa_sobol_all_y_analysis(method::Sobol, all_y::Vector{Vector{Float64}}, d::Int64, n::Int64, Ei_estimator::Symbol, y_size::Nothing, #unused#::Val{false})
@ GlobalSensitivity ~/.julia/packages/GlobalSensitivity/4fbWQ/src/sobol_sensitivity.jl:179
[10] gsa(f::var"#142#144", method::Sobol, A::Matrix{Float64}, B::Matrix{Float64}; batch::Bool, Ei_estimator::Symbol, distributed::Val{false}, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ GlobalSensitivity ~/.julia/packages/GlobalSensitivity/4fbWQ/src/sobol_sensitivity.jl:139
[11] top-level scope
@ ./timing.jl:273 [inlined]
[12] top-level scope
``````

Then trying to compare the parallel to the serial function

``````f1_batch(A)[50] = f1(A)[50]
``````

Gives me a result which is consistent for all i not just 50.

So the error when I run `sobol_result_batch = gsa(f1_batch,Sobol(),A,B, batch=true)` implies as if this is an error with with broadcasting in the batching?

@Vaibhavdixit02 do we have a test with a 3 dimensional output?

3 dimensional output is not expected to work. You can achieve what you are trying @xztt by modifying the `f1` a bit such that output remains a matrix

`````` f1 = function (p) # Para
prob_func(prob,i,repeat) = remake(prob;p=p[:,i])
ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
# Now sol[i] is the solution for the ith set of parameters
out = zeros(400,size(p,2))
for i in 1:size(p,2)
out[1:200,i] = Array(sol[i][1,:]')
out[201:400,i] = Array(sol[i][2,:]')
end
out
end
``````

now the indices will be a 400 * 4 matrix with each column representing the values for both the states, but it can easily be reshaped into the num_states*num_times output for each parameter.

1 Like

Is there a specific reason we cannot just implicitly reshape?

We could, what I meant there was currently there is no explicit handling of such output.

Thanks, this works nicely. We misunderstood how the batch output should be structured (and arrived at the conclusion that it might not be possible to have 3d outputs).

That’s why I tried to get the out[i] to contain the arrays, but messed that up. Your Array-Fu is much better than mine

This saves us a lot of work setting up plan B of running the simulations outside the GSA (in parallel) and then doing the GSA in serial on the solution array.

Fingers crossed it’ll work on the real case as well.

Disclaimer: I’m one of @xztt’s thesis supervisors.

2 Likes