gsa tutorial updated

docs/src/tutorial/gsa.md

@@ -1,25 +1,33 @@
 # Global Sensitivity Analysis
 
-This tutorial describes methods of **Global Sensitivity Analysis (GSA)** by using **HetaSimulator** together with **GlobalSensitivity** package. The content is based on [GlobalSensitivity tutorial](https://docs.sciml.ai/GlobalSensitivity/stable/tutorials/juliacon21/).
+This tutorial describes methods of **Global Sensitivity Analysis (GSA)** by using **HetaSimulator** together with **GlobalSensitivity** package. The content is based on [GlobalSensitivity tutorial](https://docs.sciml.ai/GlobalSensitivity/stable/tutorials/juliacon21/). Current version of the tutorial is based on **HetaSimulator v0.6.2** and **GlobalSensitivity v2.7.0**
 
 In order to run GSA we need to install **GlobalSensitivity** package in addition to **HetaSimulator** and **Plots**.
-We will also use **QuasiMonteCarlo** package for sampling and **CairoMakie** for plots.
+We will also use **QuasiMonteCarlo** package for sampling.
 
 ```julia
 julia> ]
-(gsa) pkg> add GlobalSensitivity QuasiMonteCarlo CairoMakie
+(gsa) pkg> add GlobalSensitivity QuasiMonteCarlo
 ```
 
 ## Working example
 
 All GSA methods are illustrated by one-compartment [Ethanol PK Model](https://github.com/insysbio/alco). 
-Download the content of */alco/src* and *alco/data-mumenthaler-2000/scenarios.csv* to your working directory.
+Download the following files to your working directory:
+```
+wget https://raw.githubusercontent.com/insysbio/alco/main/src/index.heta
+wget https://raw.githubusercontent.com/insysbio/alco/main/src/mod1.heta
+wget https://raw.githubusercontent.com/insysbio/alco/main/src/mod2.heta
+wget https://raw.githubusercontent.com/insysbio/alco/main/src/mod3.heta
+wget https://raw.githubusercontent.com/insysbio/alco/main/src/qsp-units.heta
+wget https://raw.githubusercontent.com/insysbio/alco/main/data-mumenthaler-2000/scenarios.csv
+```
 
-Load the platform into the Julia environment. You should provide the path to the modeling platform as the first argument to `load_platform`.
+Load the platform into the Julia environment. You should provide the path to the `index.heta` file as the first argument to `load_platform`.
 We will use the same working directory where the `index.heta` file is located.
 
 ```julia
-using HetaSimulator, GlobalSensitivity, CairoMakie, Plots
+using HetaSimulator, GlobalSensitivity, Plots
 
 p = load_platform(".")
 
@@ -32,33 +40,53 @@ According to Model1 scenario we will study the sensitivity of *BrAC* observable
 ```julia
 scn1 = scenarios(p)[:scn1]
 s = sim(scn1)
-Plots.plot(s, vars=[:BrAC])
+plot(s, vars=[:BrAC])
 ```
 
-![plot-default](./gsa-files/brac-default.png)
+![plot-default](./gsa-fig01.png)
 
-To run GSA we need to define a function which takes in parameter values and outputs an observable. The observable can be any function of timeseries (e.g. mean, maximum, minimum, AUC, value at a certain timepoint). We can also define vector output to study sensitivity of different observables.
-Here we are interested in maximum alcohol concentration in the compartment.
+To run GSA we need to define a function which takes in parameter values and outputs an observable. The observable can be any function of timeseries (e.g. mean, maximum, minimum, AUC, value at a certain timepoint). We can also define vector output to study sensitivity of different observables or an observed variable at different timepoints.
+
+Let's have a look at some examples of such functions. 
 
 ```julia
 pnames = [:k_a, :Vd, :Vmax, :Km]
+nump = length(pnames)
 
+# Concentration value at t=8 hours
+function C8hours_func(params)
+  sol = sim(scn1; parameters= pnames .=> params)
+  return sol(8.0,:BrAC)
+end
+
+# Concentration values at two timepoints t=4,8 hours
+function C8hours_func(params)
+  sol = sim(scn1; parameters= pnames .=> params)
+  return sol.([4.0, 8.0],:BrAC)
+end
+```
+In this tutorial we will observe maximum alcohol concentration in the compartment.
+```julia
 function Cmax_func(params)
   sol = sim(scn1; parameters= pnames .=> params)
   brac = sol.(times(sol),:BrAC)
   return maximum(brac)
 end
 ```
 
-Now we define parameters bounds based mean/sd values reported in the original paper [Ethanol Pharmacokinetics in White Women](https://doi.org/10.1111/j.1530-0277.2000.tb02103.x). (<span style="color:red">**TODO**</span>: Note about distributions)
+
+Now we define parameters bounds based mean/sd values reported in the original paper [Ethanol Pharmacokinetics in White Women](https://doi.org/10.1111/j.1530-0277.2000.tb02103.x). Parameters bound imply Uniform distribution. Usage of other distributions (Normal, LogNormal, etc) for parameters sampling will be demonstrated later in the tutorial.
 
 ```julia
 bounds = [[0.01, 0.5], [0.01, 1.], [0.01, 0.5], [0.01, 5.0]]
 ```
 
+
 ## GSA Methods
 
-### Regression/Correlation coefficients
+### Regression/Correlation coefficients method
+
+[RegressionGSA method details](https://docs.sciml.ai/GlobalSensitivity/stable/methods/regression/)
 
 The least computationally demanding method to estimate global sensitivities is to compute correlation/regression coefficients. This method is based on the assumption, that the output can be approximated by a linear model of parameters (e.g. `BrACmax = a1*k_a + a2∗Vd + a3∗Vmax + a4∗Km`). This assumption is reasonable if **the output is monotonic in each of the input parameters**, otherwise the reported coefficients can be missleading. 
 
@@ -70,45 +98,54 @@ Usually **Standard Regression Coefficients (SRC)** and **Partial Correlation Coe
 nsamples = 1000
 reg_sens = GlobalSensitivity.gsa(Cmax_func, RegressionGSA(true), bounds, samples = nsamples)
 
-fig = Figure(resolution = (600, 400))
-ax, hm = CairoMakie.heatmap(fig[1, 1],
-    reg_sens.partial_correlation,
-    axis = (xticksvisible = false, yticksvisible = false, yticks = ([1,2,3,4], String.(pnames)),
-        xticklabelsvisible = false, title = "Partial correlation"))
-Colorbar(fig[1, 2], hm)
-ax, hm = CairoMakie.heatmap(fig[2, 1],
-    reg_sens.standard_regression,
-    axis = (xticksvisible = false, yticksvisible = false, yticks = ([1,2,3,4], String.(pnames)),
-        xticklabelsvisible = false, title = "Standard regression"))
-Colorbar(fig[2, 2], hm)
-fig
+p = plot(layout=(2,1), size=(400,600))
+heatmap!(p[1],
+    reg_sens.partial_correlation, fc =cgrad([:blue, :orange]),
+    yticks=false, xticks = (1:nump, String.(pnames)),
+    title = "Partial correlation")
+
+heatmap!(p[2],
+    reg_sens.standard_regression,fc =cgrad([:blue, :orange]),
+    yticks=false, xticks = (1:nump, String.(pnames)),
+    title = "Standard correlation")
+p
 ```
-![plot-regression](./gsa-files/regression.png)
+![plot-regression](./gsa-fig02.png)
+
+**(Note!)** Regression-based GSA can also be applied to the results of multiple simulations provided to `gsa()` function as `Matrices`. See [GlobalSensitivity docs](https://docs.sciml.ai/GlobalSensitivity/stable/methods/regression/). This allows us generate parameters in advance from pre-defined Distributions.
 
-**(Note!)** Regression-based GSA can also be used with the results of multiple simulations provided to `gsa()` function as DataFrames. See [GlobalSensitivity docs](https://docs.sciml.ai/GlobalSensitivity/stable/methods/regression/).
+```julia
+# Matrix of input parameters
+X = zeros(length(pnames), nsamples)
+X[1,:] .= rand(LogNormal(-3.08, 0.77), nsamples) # k_a
+X[2,:] .= rand(Uniform(0.01, 1.), nsamples) # Vd
+X[3,:] .= rand(Uniform(0.01, 0.5), nsamples) # Vmax
+X[4,:] .= rand(LogNormal(0.46, 0.49), nsamples) # Km
+
+# output max concentration values
+Y = [Cmax_func(X[:, i]) for i in 1:size(X,2)]'
+
+reg_sens = GlobalSensitivity.gsa(X, Y, RegressionGSA(true))
+```
 
 ### Sobol (and eFAST) method
 
+[Sobol method details](https://docs.sciml.ai/GlobalSensitivity/stable/methods/sobol/)
+
 Sobol method is the GSA method of variance decomposition based on mathematical result by I.M. Sobol. Variance of the output is decomposed into contributions from each parameter and their interactions. Two indices are used to measure these contributions. The contribution of each parameter is represented by **first-order index** and the contribution of the parameter together with all its interactions with other parameters is represented by **total-order index**. 
 
-The requirement for the number of samples is the same as for Correlatin/Regression coefficients method: The number of samples required by the method should be larger than 2^P (where P is the number of input parameters).
+The requirement for the number of samples is the same as for Correlatin/Regression coefficients method: **The number of samples required by the method should be larger than 2^P (where P is the number of input parameters).**
 
 ```julia
 nsamples = 1000
 sobol_sens1 = GlobalSensitivity.gsa(Cmax_func, Sobol(), bounds, samples = nsamples)
 
-fig = Figure(resolution = (600, 400))
-barplot(fig[1, 1],
-    [1, 2, 3, 4],
-    sobol_sens1.S1,
-    color = :green,
-    axis = (xticks = (1:4, String.(pnames)),
-        title = "Cmax (Sobol)", ylabel = "First order"))
-barplot(fig[2, 1], [1, 2, 3, 4], sobol_sens1.ST, color = :green,
-    axis = (xticks = (1:4, String.(pnames)), ylabel = "Total order"))
-fig
+p = plot(layout=(2,1), size=(400,600), xticks = (1:nump, String.(pnames)))
+bar!(p[1], 1:nump, sobol_sens1.S1, title = "Sobol First order indeces", label=false)
+bar!(p[2], 1:nump, sobol_sens1.ST, title = "Sobol Total order indeces", label=false)
+p
 ```
-![plot-sobol1](./gsa-files/sobol1.png)
+![plot-sobol1](./gsa-fig03.png)
 
 **(Note 1 !)** It is recommended to use low-discrepancy sequence instead of a pure Monte-Carlo to effectively sample the search space . One can choose `SobolSample()`, `LatinHypercubeSample()`, `HaltonSample()` or other quasi-Monte Carlo samplers from `QuasiMonteCarlo` package. 
 
@@ -123,15 +160,35 @@ A,B = QuasiMonteCarlo.generate_design_matrices(nsamples,lb,ub,sampler)
 sobol_sens2 = GlobalSensitivity.gsa(Cmax_func, Sobol(), A, B)
 ```
 
-**(Note 2 !)** To monitor the convergence of Sobol indices it is usefull to set different sample sizes and visualize how the indices' values stabilize with the increase of sample size. 
+**(Note 2 !)** To ensure the convergence of Sobol indices and choose reasonable sample size it is useful to set different sample sizes and monitor (plot) how the indices' values stabilize with the increase of sample size. 
+
+**(Note 3 !)** As of GlobalSensitivity@v2.7.0 `Sobol()` implementation supports  Uniform input ranges and not parameters' distributions. If you prefer to input parameters' distributions you can choose `eFAST()`, which is based on `Sobol()` method. Please refer to the [eFAST docs](https://docs.sciml.ai/GlobalSensitivity/stable/methods/efast/) for detailed example.
+
+### Morris method (Elementary effects method)
 
-<span style="color:red">**TODO**</span>
+[Morris method details](https://docs.sciml.ai/GlobalSensitivity/stable/methods/morris/)
 
+Morris method estimates parameters' contribution to the output be computing individual "elementary effects" (EE), which are local sensitivities at different points in the parameters space. The method runs multiple trajectories through the grid of points in parameter space, computes EEs and outputs mean and std of EE per parameter. The sensitivities are calculated per parameter, so the method doesn't take into account interactions between parameters. Morris method doesn't have a strong mathematical foundation (as Sobol method). Due to its simplicity, it is commonly used as a pre-scan method with all the input parameters of interest. Parameters identified as noninfluential can be then fixed to run computationally demanding but more reliable variance-based methods (`Sobol()`, `eFAST()`)
+
+The requirement for the number of samples is the same as for Correlatin/Regression coefficients method: **The number of samples required by the method should be larger than 2^P (where P is the number of input parameters).**
+
+The number of estimated points in the parameter space is controlled by `total_num_trajectory` and `num_trajectory` arguments. See [Morris method docs](https://docs.sciml.ai/GlobalSensitivity/stable/methods/morris/) for details. 
 ```julia
+morris_sens = GlobalSensitivity.gsa(Cmax_func, Morris(total_num_trajectory=5000,num_trajectory=600), bounds)
 
+bar(1:nump, morris_sens.means_star[1,:], title = "Morris method", label=false)
 ```
+![plot-morris](./gsa-fig04.png)
+
+**(Note 1 !)** Morris method also outputs std of EE per parameter: `morris_sens.variances`. High std reported by Morris method can be a marker of parameters interaction or high nonlinearity between the parameters and the output. 
+
+**(Note 2 !)** As in `Sobol()` method it useful to set different number of trajectories and monitor how the mean values stabilize with the increase of sample size. 
 
-**(Note 3 !)** It is natural to speed-up GSA (e.g. Sobol method) by running simulations in parallel. To do so we need to rewrite `Cmax_func` and use `batch` keyword argument in the `gsa()` function.
+**(Note 3 !)** As of GlobalSensitivity@v2.7.0 `Morris()` implementation doesn't support `generate_design_matrices` interface (see Sobol method notes), so only default sampling scheme (Monte-Carlo) can be used.
+
+## Parallel setup
+
+It is natural to speed-up GSA methods (e.g. Sobol method) by running simulations in parallel. To do so we need to rewrite `Cmax_func` and use `batch` keyword argument in the `gsa()` function.
 
 ```julia
 function Cmax_batch_func(params_batch)
@@ -147,4 +204,4 @@ function Cmax_batch_func(params_batch)
   end
 
   sobol_sens = GlobalSensitivity.gsa(Cmax_batch_func, Sobol(), A, B, batch=true)
-```
+```