# Parallelize GSA on continuous outputs

**URL:** <https://discourse.julialang.org/t/parallelize-gsa-on-continuous-outputs/95964>\
**Category:** Modelling & Simulations\
**Tags:** diffeq\
**Created:** [March 12, 2023, 4:15pm UTC](https://discourse.julialang.org/t/parallelize-gsa-on-continuous-outputs/95964 "2023-03-12T16:15:15Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![xztt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xztt/32/33843_2.png) [@xztt](https://discourse.julialang.org/u/xztt)\
**Post date:** [March 12, 2023, 4:15pm UTC](https://discourse.julialang.org/t/parallelize-gsa-on-continuous-outputs/95964/1 "2023-03-12T16:15:15Z")

</div>

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

```julia
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

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 12, 2023, 5:07pm UTC](https://discourse.julialang.org/t/parallelize-gsa-on-continuous-outputs/95964/2 "2023-03-12T17:07:54Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![xztt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xztt/32/33843_2.png) [@xztt](https://discourse.julialang.org/u/xztt)\
**Post date:** [March 12, 2023, 5:41pm UTC](https://discourse.julialang.org/t/parallelize-gsa-on-continuous-outputs/95964/3 "2023-03-12T17:41:24Z")

</div>

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.

> [@Time Series sensitivity analysis for multiple outputs](https://discourse.julialang.org/t/time-series-sensitivity-analysis-for-multiple-outputs/80631):
>
> Hi, function circ\_local\_time(x) \_prob = remake(prob,u0=x[1:3],p=x[4:end]) sol = solve(\_prob,Tsit5(), reltol=1e-6, abstol=1e-6,saveat=savetime) sol[LV.p] end St = ForwardDiff.jacobian(circ\_local\_time,[u0;p]) x = LinRange(0,1,100) plot(x,St[:,12]) When I run the above code I can see how the sensitivity of the parameter changes during a cycle. However, when I add another output eg function circ\_local\_time(x) \_prob = remake(prob,u0=x[1:3],p=x[4:end]) sol = solve(\_prob,Tsit5()…

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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 12, 2023, 5:47pm UTC](https://discourse.julialang.org/t/parallelize-gsa-on-continuous-outputs/95964/4 "2023-03-12T17:47:46Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![xztt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xztt/32/33843_2.png) [@xztt](https://discourse.julialang.org/u/xztt)\
**Post date:** [March 13, 2023, 9:33am UTC](https://discourse.julialang.org/t/parallelize-gsa-on-continuous-outputs/95964/5 "2023-03-13T09:33:08Z")

</div>

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

```julia
  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

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 13, 2023, 10:46am UTC](https://discourse.julialang.org/t/parallelize-gsa-on-continuous-outputs/95964/6 "2023-03-13T10:46:30Z")

</div>

What’s the issue?

---

<div class="post-metadata">

**Author:** ![xztt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xztt/32/33843_2.png) [@xztt](https://discourse.julialang.org/u/xztt)\
**Post date:** [March 13, 2023, 11:05am UTC](https://discourse.julialang.org/t/parallelize-gsa-on-continuous-outputs/95964/7 "2023-03-13T11:05:12Z")

</div>

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.

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![xztt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xztt/32/33843_2.png) [@xztt](https://discourse.julialang.org/u/xztt)\
**Post date:** [March 13, 2023, 2:48pm UTC](https://discourse.julialang.org/t/parallelize-gsa-on-continuous-outputs/95964/8 "2023-03-13T14:48:19Z")

</div>

I have been playing around with a couple of things,

If I write

```julia
  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

```julia
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

```julia
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?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 14, 2023, 1:22am UTC](https://discourse.julialang.org/t/parallelize-gsa-on-continuous-outputs/95964/9 "2023-03-14T01:22:47Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Vaibhavdixit02](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vaibhavdixit02/32/2916_2.png) [@Vaibhavdixit02](https://discourse.julialang.org/u/Vaibhavdixit02)\
**Post date:** [March 14, 2023, 7:42pm UTC](https://discourse.julialang.org/t/parallelize-gsa-on-continuous-outputs/95964/10 "2023-03-14T19:42:36Z")

</div>

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

```julia
 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.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [March 14, 2023, 8:06pm UTC](https://discourse.julialang.org/t/parallelize-gsa-on-continuous-outputs/95964/11 "2023-03-14T20:06:33Z")

</div>

> [@Vaibhavdixit02](#):
>
> 3 dimensional output is not expected to work.

Is there a specific reason we cannot just implicitly reshape?

---

<div class="post-metadata">

**Author:** ![Vaibhavdixit02](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/vaibhavdixit02/32/2916_2.png) [@Vaibhavdixit02](https://discourse.julialang.org/u/Vaibhavdixit02)\
**Post date:** [March 15, 2023, 5:04am UTC](https://discourse.julialang.org/t/parallelize-gsa-on-continuous-outputs/95964/12 "2023-03-15T05:04:17Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![TS-CUBED](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ts-cubed/32/25292_2.png) [@TS-CUBED](https://discourse.julialang.org/u/TS-CUBED)\
**Post date:** [March 16, 2023, 5:29pm UTC](https://discourse.julialang.org/t/parallelize-gsa-on-continuous-outputs/95964/13 "2023-03-16T17:29:51Z")

</div>

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.