Global sensitivity analysis - following time series

Specific Domains Modelling & Simulations
1

I tested the example in 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:

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,

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.

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

2

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.

1 Like
3

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.

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

4

If you batch you have to change your function.

5 

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

6

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

1 Like
7

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:

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

8

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

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

oh I think I had the indexing flipped:

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)
1 Like
10

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

Works perfect!

A

11

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

12

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

13

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.


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

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.

1 Like
15

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

16

DiffEqSensitivity