# Global sensitivity analysis - following time series

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

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

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

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]

Stacktrace:
[2] macro expansion at .\threadingconstructs.jl:69 [inlined]
[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]
[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
[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

If you batch you have to change your function.

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

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

1 Like

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

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

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

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

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

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

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