Parallelize GSA on continuous outputs

Specific Domains Modelling & Simulations
1

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 :slight_smile:

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)
  sol = solve(ensemble_prob,Tsit5(),EnsembleThreads();saveat=t,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] = maximum(sol[i][2,:])
  end
  out
end
2

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

3

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.

4

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?

5

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)
    sol = solve(ensemble_prob,Tsit5(),EnsembleThreads();saveat=t,trajectories=size(p,2))
    # 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

6

What’s the issue?

7

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)
    sol = solve(ensemble_prob,Tsit5(),EnsembleThreads();saveat=t,trajectories=size(p,2))
    # 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)
8

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)
    sol = solve(ensemble_prob,Tsit5(),EnsembleThreads();saveat=t,trajectories=size(p,2))
    # 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
  *(::StridedVecOrMat, ::Adjoint{<:Any, <:LinearAlgebra.LQPackedQ})
   @ LinearAlgebra /opt/julia-1.9.0-rc1/share/julia/stdlib/v1.9/LinearAlgebra/src/lq.jl:269
  ...

Stacktrace:
  [1] _broadcast_getindex_evalf
    @ ./broadcast.jl:683 [inlined]
  [2] _broadcast_getindex
    @ ./broadcast.jl:656 [inlined]
  [3] getindex
    @ ./broadcast.jl:610 [inlined]
  [4] copy
    @ ./broadcast.jl:912 [inlined]
  [5] materialize
    @ ./broadcast.jl:873 [inlined]
  [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
    @ ~/Downloads/MWE_GSA_TS.jl:0

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?

9

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

10

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)
         sol = solve(ensemble_prob,Tsit5(),EnsembleThreads();saveat=t,trajectories=size(p,2))
         # 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
11

Is there a specific reason we cannot just implicitly reshape?

12

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

13

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 :wink:

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