# Global sensitivity analysis - following time series

**URL:** https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770
**Category:** Modelling & Simulations
**Created:** [March 9, 2020, 8:52pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770 "2020-03-09T20:52:30Z")
**Posts on this page:** 16
**Page:** 1

<div class="post-metadata">

### Author: ![alewolf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alewolf/32/9000_2.png) [@alewolf](https://discourse.julialang.org/u/alewolf)
#### Post date: [March 9, 2020, 8:52pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770/1 "2020-03-09T20:52:30Z")

</div>

I tested the example in [https://docs.juliadiffeq.org/stable/analysis/global\_sensitivity/#GSA-examples-1](https://docs.juliadiffeq.org/stable/analysis/global_sensitivity/#GSA-examples-1).  
As expected it all went well.  
This example, calculates sensitivity with regard to the mean over time on u1 and the maximum over time on u2.  
I am interesting on looking how sensitivities change over time - to do this I made a simple change:

```julia
using DiffEqSensitivity, Statistics, OrdinaryDiffEq #load packages
function f(du,u,p,td
  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)
  prob1 = remake(prob;p=p)
  sol = solve(prob1,Tsit5();saveat=t)
  #[mean(sol[1,:]), maximum(sol[2,:])] ---> just use sol as function result
end

```

I was expecting this would allow to follow the time course of sensitivities.  
it does not work and returns the error at the bottom.  
If I change the function to include sensitivities only fro u1 or u2 (below) it works,

```julia
f1 = function (p)
  prob1 = remake(prob;p=p)
  sol = solve(prob1,Tsit5();saveat=t)
  #[mean(sol[1,:]), maximum(sol[2,:])] ---> just use sol as function result
  sol[1,:]# or sol[2,:] as output will work but not together (not just sol)
end

```

I don’t mind running f1 once for sol[1,:] and then sol[2,:], but then they will be visiting different sections of the parameter space (2 independent runs) so I was hoping to run all sensitivities at once. I remember this was doable and documented in old versions of the package - is this a bug or some intended change.

(old code - where senstivities over time where returned for all states/parameters.

```julia
m = DiffEqSensitivity.morris_sensitivity(Problem,solver,t,p_ranges,p_steps;len_trajectory=150,total_num_trajectory=100,num_trajectory=15, relative_scale=rel_scale)

```

Thanks,

A

This is the error of running code as is:

DimensionMismatch(“tried to assign 2×200 array to 2×200 destination”)

Stacktrace:  
[1] throw\_setindex\_mismatch(::ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{,Tuple{}}},DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}},DiffEqBase.DEStats}, ::Tuple{Int64,Int64}) at .\indices.jl:169  
[2] setindex\_shape\_check at .\indices.jl:226 [inlined]  
[3] macro expansion at .\multidimensional.jl:722 [inlined]  
[4] \_unsafe\_setindex!(::IndexLinear, ::Array{Float64,2}, ::ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{,Tuple{}}},DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}},DiffEqBase.DEStats}, ::Base.Slice{Base.OneTo{Int64}}, ::UnitRange{Int64}) at .\multidimensional.jl:717  
[5] \_setindex! at .\multidimensional.jl:712 [inlined]  
[6] setindex! at .\abstractarray.jl:1074 [inlined]  
[7] \_typed\_hcat(::Type{Float64}, ::Array{ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{,Tuple{}}},DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}},DiffEqBase.DEStats},1}) at .\abstractarray.jl:1343  
[8] reduce(::typeof(hcat), ::Array{ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{,Tuple{}}},DiffEqBase.StandardODEProblem},Tsit5,OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Tsit5Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.Tsit5ConstantCache{Float64,Float64}}},DiffEqBase.DEStats},1}) at .\abstractarray.jl:1378  
[9] #gsa#103(::Bool, ::typeof(gsa), ::getfield(Main, Symbol(“##35#36”)), ::Morris, ::Array{Array{Int64,1},1}) at C:\Users\awolf-yadlin.juliapro\JuliaPro\_v1.2.0-2\packages\DiffEqSensitivity\9ybQN\src\global\_sensitivity\morris\_sensitivity.jl:87  
[10] gsa(::Function, ::Morris, ::Array{Array{Int64,1},1}) at C:\Users\awolf-yadlin.juliapro\JuliaPro\_v1.2.0-2\packages\DiffEqSensitivity\9ybQN\src\global\_sensitivity\morris\_sensitivity.jl:69  
[11] top-level scope at In[74]:1

---

<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 11, 2020, 7:29pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770/2 "2020-03-11T19:29:58Z")

</div>

> [@alewolf](#):
>
> If I change the function to include sensitivities only fro u1 or u2 (below) it works,

It wants a vector. It’s a vector → vector function which then computes the global sensitivity of each output w.r.t.the inputs. You could use that to reshape to a timeseries on your own though.

---

<div class="post-metadata">

### Author: ![alewolf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alewolf/32/9000_2.png) [@alewolf](https://discourse.julialang.org/u/alewolf)
#### Post date: [March 12, 2020, 3:04am UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770/3 "2020-03-12T03:04:30Z")

</div>

On the same page-  
I tried to run the sobol parallel code you have in the documentation. it did not work. There mightbe a bug there.

```julia
N = 10000
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(N,lb,ub,sampler)

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

f2 = 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) # Case1
   # sol = solve(ensemble_prob,Tsit5(),EnsembleThreads();saveat=t, trajectories=500 ) # Case2
  # Now sol[i] is the solution for the ith set of parameters
  out = zeros(size(p,1),2)
  for i in 1:size(p,1)
    out[i,1] = mean(sol[i][1,:])
    out[i,2] = maximum(sol[i][2,:])
  end
  out
end

sobol_result = gsa(f2,Sobol(),A,B,batch=true)

```

returns in case1 and case2 (with and w/o trajectories respectively)

Is there something I am doing wrong, I copied the example verbatim - except for naming the fucntion f2.

Thanks,

Ale

Case1:

UndefKeywordError: keyword argument trajectories not assigned

Stacktrace:  
[1] (::DiffEqBase.var"#kw##\_\_solve")(::NamedTuple{(:saveat,),Tuple{Array{Float64,1}}}, ::typeof(DiffEqBase.\_\_solve), ::EnsembleProblem{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},var"#prob\_func#9"{Array{Float64,2}},DiffEqBase.var"#333#339",DiffEqBase.var"#335#341",Array{Any,1}}, ::Tsit5, ::EnsembleThreads) at .\none:0  
[2] #solve#445(::Base.Iterators.Pairs{Symbol,Array{Float64,1},Tuple{Symbol},NamedTuple{(:saveat,),Tuple{Array{Float64,1}}}}, ::typeof(solve), ::EnsembleProblem{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},var"#prob\_func#9"{Array{Float64,2}},DiffEqBase.var"#333#339",DiffEqBase.var"#335#341",Array{Any,1}}, ::Tsit5, ::Vararg{Any,N} where N) at C:\Users\awolf-yadlin.juliapro\JuliaPro\_v1.3.1-1\packages\DiffEqBase\mDFok\src\solve.jl:77  
[3] (::DiffEqBase.var"#kw##solve")(::NamedTuple{(:saveat,),Tuple{Array{Float64,1}}}, ::typeof(solve), ::EnsembleProblem{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},var"#prob\_func#9"{Array{Float64,2}},DiffEqBase.var"#333#339",DiffEqBase.var"#335#341",Array{Any,1}}, ::Tsit5, ::EnsembleThreads) at .\none:0  
[4] (::var"#7#8")(::Array{Float64,2}) at .\In[9]:4  
[5] #gsa#125(::Bool, ::Symbol, ::typeof(gsa), ::var"#7#8", ::Sobol, ::Array{Float64,2}, ::Array{Float64,2}) at C:\Users\awolf-yadlin.juliapro\JuliaPro\_v1.3.1-1\packages\DiffEqSensitivity\qOHEV\src\global\_sensitivity\sobol\_sensitivity.jl:33  
[6] (::DiffEqSensitivity.var"#kw##gsa")(::NamedTuple{(:batch,),Tuple{Bool}}, ::typeof(gsa), ::Function, ::Sobol, ::Array{Float64,2}, ::Array{Float64,2}) at .\none:0  
[7] top-level scope at In[11]:1

Case2 (adding trajectories)

TaskFailedException:  
BoundsError: attempt to access 4×60000 Array{Float64,2} at index [5, Base.Slice(Base.OneTo(60000))]  
Stacktrace:  
[1] throw\_boundserror(::Array{Float64,2}, ::Tuple{Int64,Base.Slice{Base.OneTo{Int64}}}) at .\abstractarray.jl:538  
[2] checkbounds at .\abstractarray.jl:503 [inlined]  
[3] \_getindex at .\multidimensional.jl:669 [inlined]  
[4] getindex at .\abstractarray.jl:981 [inlined]  
[5] prob\_func at .\In[12]:2 [inlined]  
[6] macro expansion at C:\Users\awolf-yadlin.juliapro\JuliaPro\_v1.3.1-1\packages\DiffEqBase\mDFok\src\ensemble\basic\_ensemble\_solve.jl:155 [inlined]  
[7] (::DiffEqBase.var"#364#threadsfor\_fun#368"{EnsembleProblem{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},var"#prob\_func#12"{Array{Float64,2}},DiffEqBase.var"#333#339",DiffEqBase.var"#335#341",Array{Any,1}},Tsit5,UnitRange{Int64},Tuple{Pair{Symbol,Array{Float64,1}}},Array{Any,1},Base.OneTo{Int64}})(::Bool) at .\threadingconstructs.jl:61  
[8] (::DiffEqBase.var"#364#threadsfor\_fun#368"{EnsembleProblem{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},var"#prob\_func#12"{Array{Float64,2}},DiffEqBase.var"#333#339",DiffEqBase.var"#335#341",Array{Any,1}},Tsit5,UnitRange{Int64},Tuple{Pair{Symbol,Array{Float64,1}}},Array{Any,1},Base.OneTo{Int64}})() at .\threadingconstructs.jl:28

Stacktrace:  
[1] wait(::Task) at .\task.jl:251  
[2] macro expansion at .\threadingconstructs.jl:69 [inlined]  
[3] solve\_batch(::EnsembleProblem{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},var"#prob\_func#12"{Array{Float64,2}},DiffEqBase.var"#333#339",DiffEqBase.var"#335#341",Array{Any,1}}, ::Tsit5, ::EnsembleThreads, ::UnitRange{Int64}, ::Int64, ::Pair{Symbol,Array{Float64,1}}) at C:\Users\awolf-yadlin.juliapro\JuliaPro\_v1.3.1-1\packages\DiffEqBase\mDFok\src\ensemble\basic\_ensemble\_solve.jl:152  
[4] macro expansion at C:\Users\awolf-yadlin.juliapro\JuliaPro\_v1.3.1-1\packages\DiffEqBase\mDFok\src\ensemble\basic\_ensemble\_solve.jl:91 [inlined]  
[5] macro expansion at .\util.jl:212 [inlined]  
[6] #\_\_solve#356(::Int64, ::Int64, ::Int64, ::Base.Iterators.Pairs{Symbol,Array{Float64,1},Tuple{Symbol},NamedTuple{(:saveat,),Tuple{Array{Float64,1}}}}, ::typeof(DiffEqBase.\_\_solve), ::EnsembleProblem{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},var"#prob\_func#12"{Array{Float64,2}},DiffEqBase.var"#333#339",DiffEqBase.var"#335#341",Array{Any,1}}, ::Tsit5, ::EnsembleThreads) at C:\Users\awolf-yadlin.juliapro\JuliaPro\_v1.3.1-1\packages\DiffEqBase\mDFok\src\ensemble\basic\_ensemble\_solve.jl:85  
[7] (::DiffEqBase.var"#kw##\_\_solve")(::NamedTuple{(:saveat, :trajectories),Tuple{Array{Float64,1},Int64}}, ::typeof(DiffEqBase.\_\_solve), ::EnsembleProblem{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},var"#prob\_func#12"{Array{Float64,2}},DiffEqBase.var"#333#339",DiffEqBase.var"#335#341",Array{Any,1}}, ::Tsit5, ::EnsembleThreads) at .\none:0  
[8] #solve#445 at C:\Users\awolf-yadlin.juliapro\JuliaPro\_v1.3.1-1\packages\DiffEqBase\mDFok\src\solve.jl:77 [inlined]  
[9] (::DiffEqBase.var"#kw##solve")(::NamedTuple{(:saveat, :trajectories),Tuple{Array{Float64,1},Int64}}, ::typeof(solve), ::EnsembleProblem{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,Array{Float64,1},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},var"#prob\_func#12"{Array{Float64,2}},DiffEqBase.var"#333#339",DiffEqBase.var"#335#341",Array{Any,1}}, ::Tsit5, ::EnsembleThreads) at .\none:0  
[10] (::var"#10#11")(::Array{Float64,2}) at .\In[12]:4  
[11] #gsa#125(::Bool, ::Symbol, ::typeof(gsa), ::var"#10#11", ::Sobol, ::Array{Float64,2}, ::Array{Float64,2}) at C:\Users\awolf-yadlin.juliapro\JuliaPro\_v1.3.1-1\packages\DiffEqSensitivity\qOHEV\src\global\_sensitivity\sobol\_sensitivity.jl:33  
[12] (::DiffEqSensitivity.var"#kw##gsa")(::NamedTuple{(:batch,),Tuple{Bool}}, ::typeof(gsa), ::Function, ::Sobol, ::Array{Float64,2}, ::Array{Float64,2}) at .\none:0  
[13] top-level scope at In[14]:1

---

<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, 2020, 12:15pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770/4 "2020-03-12T12:15:46Z")

</div>

If you batch you have to change your function.

---

<div class="post-metadata">

### Author: ![alewolf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alewolf/32/9000_2.png) [@alewolf](https://discourse.julialang.org/u/alewolf)
#### Post date: [March 12, 2020, 4:23pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770/5 "2020-03-12T16:23:24Z")

</div>

```julia

f2 = 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) # Case1
   # sol = solve(ensemble_prob,Tsit5(),EnsembleThreads();saveat=t, trajectories=500 ) # Case2
  # Now sol[i] is the solution for the ith set of parameters
  out = zeros(size(p,1),2)
  for i in 1:size(p,1)
    out[i,1] = mean(sol[i][1,:])
    out[i,2] = maximum(sol[i][2,:])
  end
  out
end

```

Isn’t this the function you are using for batch in the example?  
This is the one I am using, I have tried to modifiy it, but I can’t get this to work - I did include the trajectories argument, in the ensemble defiinition otherwise it doesn’t even attempt to solve

---

<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, 2020, 2:02am UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770/6 "2020-03-13T02:02:05Z")

</div>

oh that should be `sol = solve(ensemble_prob,Tsit5(),EnsembleThreads();saveat=t,trajectories=size(p,1))`.

---

<div class="post-metadata">

### Author: ![alewolf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alewolf/32/9000_2.png) [@alewolf](https://discourse.julialang.org/u/alewolf)
#### Post date: [March 13, 2020, 4:58am UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770/7 "2020-03-13T04:58:46Z")

</div>

> [@ChrisRackauckas](#):
>
> solve(ensemble\_prob,Tsit5(),EnsembleThreads();saveat=t,trajectories=size(p,1))

Sorry to keep at it I tried the fix you suggested -  
Still does not work for me, I am showing full code and error below.

Still get bounds error:

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

f2 = 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,1))
  # Now sol[i] is the solution for the ith set of parameters
  out = zeros(size(p,1),2)
  for i in 1:size(p,1)
    out[i,1] = mean(sol[i][1,:])
    out[i,2] = maximum(sol[i][2,:])
  end
  out
end

sobol_result = gsa(f2,Sobol(),A,B,batch=true)

```

BoundsError: attempt to access 4×2 Array{Float64,2} at index [Base.Slice(Base.OneTo(4)), 1:20000]

Stacktrace:  
[1] throw\_boundserror(::Array{Float64,2}, ::Tuple{Base.Slice{Base.OneTo{Int64}},UnitRange{Int64}}) at .\abstractarray.jl:538  
[2] checkbounds at .\abstractarray.jl:503 [inlined]  
[3] \_getindex at .\multidimensional.jl:669 [inlined]  
[4] getindex(::Array{Float64,2}, ::Function, ::UnitRange{Int64}) at .\abstractarray.jl:981  
[5] #gsa#125(::Bool, ::Symbol, ::typeof(gsa), ::var"#7#8", ::Sobol, ::Array{Float64,2}, ::Array{Float64,2}) at C:\Users\awolf-yadlin.juliapro\JuliaPro\_v1.3.1-1\packages\DiffEqSensitivity\qOHEV\src\global\_sensitivity\sobol\_sensitivity.jl:66  
[6] (::DiffEqSensitivity.var"#kw##gsa")(::NamedTuple{(:batch,),Tuple{Bool}}, ::typeof(gsa), ::Function, ::Sobol, ::Array{Float64,2}, ::Array{Float64,2}) at .\none:0  
[7] top-level scope at In[14]:1

---

<div class="post-metadata">

### Author: ![alewolf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alewolf/32/9000_2.png) [@alewolf](https://discourse.julialang.org/u/alewolf)
#### Post date: [March 13, 2020, 9:57pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770/8 "2020-03-13T21:57:08Z")

</div>

Should the code have a double loop?  
This does not crash, but I am not 100% sure I am looking at the right thing anymore. Results are veru different from serial analysis

```julia
function f(du,u,p,t)
# @inbounds begin
  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
# nothing
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))

f2 = function (p)
  prob_func(prob,i,repeat) = remake(prob;p=p[i,:])
    println(prob,p)
  ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
  sol = solve(ensemble_prob,Tsit5(),EnsembleThreads();saveat=t,trajectories=size(p,1))
  # Now sol[i] is the solution for the ith set of parameters
  out = zeros(size(p,2),size(p,1),2)
  println(size(p))
  for j in 1:size(p,2)
    for i in 1:size(p,1)
      if i==1
          println(mean(sol[i][1,:]))
          println(sol[i][1,:])    
          println(size(sol[i][1,:]), typeof(sol[i][1,:]))
      end
      out[j,i,1] = mean(sol[i][1,:])
      out[j,i,2] = maximum(sol[i][2,:])
    end
  end
    println(out)
    println(size(out), typeof(out))
  out
end

sobol_result = gsa(f2,Sobol(),A,B,batch=true)

```

---

<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, 2020, 8:29am UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770/9 "2020-03-14T08:29:29Z")

</div>

oh I think I had the indexing flipped:

```julia
using OrdinaryDiffEq, DiffEqSensitivity
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))

f2 = 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

using QuasiMonteCarlo, Statistics

N = 10000
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(N,lb,ub,sampler)

sobol_result = gsa(f2,Sobol(),A,B,batch=true)

```

---

<div class="post-metadata">

### Author: ![alewolf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alewolf/32/9000_2.png) [@alewolf](https://discourse.julialang.org/u/alewolf)
#### Post date: [March 14, 2020, 3:03pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770/10 "2020-03-14T15:03:12Z")

</div>

Thanks!

I’ll try this as soon as I have 5 min today.

I still get quite confused with indexing for parallel or ensemble computation. So I appreciate you taking the time to help 🙂

Works perfect!

A

---

<div class="post-metadata">

### Author: ![alewolf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alewolf/32/9000_2.png) [@alewolf](https://discourse.julialang.org/u/alewolf)
#### Post date: [March 18, 2020, 5:07pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770/11 "2020-03-18T17:07:35Z")

</div>

Last question on this subject -  
I am able to tun Sobol on my real model (44 equations, 65 params)  
As long as I keep N below 10,000. This seems like a small number for such a big set of params. When I tried 50000, my RAM collapsed (usage 100%). I am using a system, provided at work, with 64GB RAM. Is there a way other than increasing RAM (not very feasible now) to avoid overwhelming the RAM? ie decrease memory use when running GSA? I’d appreciate any ideas!

Thanks,

A

---

<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 18, 2020, 5:16pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770/12 "2020-03-18T17:16:37Z")

</div>

> [@alewolf](#):
>
> As long as I keep N below 10,000. This seems like a small number for such a big set of params. When I tried 50000, my RAM collapsed (usage 100%). I am using a system, provided at work, with 64GB RAM

How many timepoints saved, i.e. the size of the output space? 50,000 \* output space would be the size of the allocations.

---

<div class="post-metadata">

### Author: ![alewolf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alewolf/32/9000_2.png) [@alewolf](https://discourse.julialang.org/u/alewolf)
#### Post date: [March 18, 2020, 5:31pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770/13 "2020-03-18T17:31:19Z")

</div>

Here is the code:  
It is similar to the one you were helping me on the small example.  
I am doing sensitivity against the max of one variable and the mean of other, I am saving at 21 times to do that.

```julia

t = collect(range(0, stop=1, length=21))
condition1(u,t,integrator) = u[38] < 1
affect1!(integrator) = terminate!(integrator)
cb1 = DiscreteCallback(condition1, affect1!)

f2 = function (p)
  prob_func(PKPD_Prob,i,repeat) = remake(PKPD_Prob;p=p[:,i])
  ensemble_prob = EnsembleProblem(PKPD_Prob,prob_func=prob_func)
  sol = DifferentialEquations.solve(ensemble_prob,Rodas5(),EnsembleThreads(), callback=cb1;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][33,:])
    out[2,i] = maximum(sol[i][44,:])
  end
  out
end

#p is my original parameters value

p_ranges=fill(Float64[],65)
[p_ranges[i] = [p[i]/2,p[i]*2] for i in 1:16] 
p_ranges[17] = [p[17],p[17]*4]   
p_ranges[18] = [p[18]/1,p[18]*1]     
p_ranges[19] = [p[19]/5,p[19]*5] 
p_ranges[20] = [p[20],p[20]*3]   
p_ranges[21] = [p[21]/5,p[21]*5]   
[p_ranges[i] = [p[i]/2,p[i]*2] for i in 22:36] 
[p_ranges[i] = [p[i]/3,p[i]*3] for i in 37:40] 
[p_ranges[i] = [p[i]/3,p[i]*3] for i in 41:42] 
[p_ranges[i] = [p[i]/3,p[i]*3] for i in 43:48] 
[p_ranges[i] = [p[i]/2,p[i]*5] for i in [49,51,52,54]] 
[p_ranges[i] = [p[i]/2,p[i]*5] for i in [50,53]] 
[p_ranges[i] = [p[i]/2,p[i]*2] for i in 55:56] 
[p_ranges[i] = [p[i]/3,p[i]*3] for i in 57:60] 
[p_ranges[i] = [p[i]/5,p[i]*5] for i in 61:65] 

N = 25000
lb = [x[1] for x in p_ranges]
ub = [x[2] for x in p_ranges]
sampler = SobolSample()
A,B = QuasiMonteCarlo.generate_design_matrices(N,lb,ub,sampler)

println("GSA PARAMS LOADED")

sobol_result = gsa(f2,Sobol(),A,B,batch=true)

```

---

<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 18, 2020, 5:51pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770/14 "2020-03-18T17:51:41Z")

</div>

Open an issue, but I don’t plan to get to it soon. You can do the math I suggested above and see if like it’s a multiple of 4 or something. If it’s a factor higher than that, there’s something to investigate.

---

<div class="post-metadata">

### Author: ![alewolf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alewolf/32/9000_2.png) [@alewolf](https://discourse.julialang.org/u/alewolf)
#### Post date: [March 18, 2020, 5:55pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770/15 "2020-03-18T17:55:36Z")

</div>

Will try, doing the math as above and will open issue.  
You want the issue opened in DiffEqSensitivity? or DiffEqEnsemble?

---

<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 18, 2020, 5:56pm UTC](https://discourse.julialang.org/t/global-sensitivity-analysis-following-time-series/35770/16 "2020-03-18T17:56:16Z")

</div>

DiffEqSensitivity