Mutual Information Sensitivity (#176)

* Add KS Rank sensitivity method

* Delete settings.json

* Update ks_rank_sensitivity.jl

add raw string annotation

* Update tests

KSRank is stochastic and slight variations in the outcome measures are possible. New tests check whether parameters are correctly identified as sensitive.

* Update tests and docs

Make tests less sensitive to random variability in the sensitivity criterion. Also indicate the random variability influence in the documentation.

* Formatting

Use JuliaFormatter

* Update docs

Add more detailed explanation of f(Y) function signature, and fix typos in math.

* Update ks_rank_sensitivity.jl

Add description of returned KSRankResult struct

* Create docs page for KS Rank method

* Update ks_rank_sensitivity.jl

Add @doc signature

* Update ks_rank_sensitivity.jl

Change math $ into ``

* Update ks_rank_sensitivity.jl

Use correct link

* Rename ks_rank to rsa and fix docs

Renamed KS_Rank to RSA in line with how the method is known in literature.

* Change name of docs page

* Bump lower bound of StatsBase to 0.33.7

fix compatibility issues with SortingAlgorithms.jl

* Make tests faster and modify default

Modify default argument of n_dummy_parameters to 10

* format

* Add mutual information method

Add mutual information sensitivity analysis method from Lüdtke, N., Panzeri, S., Brown, M., Broomhead, D. S., Knowles, J., Montemurro, M. A., & Kell, D. B. (2007). Information-theoretic sensitivity analysis: a general method for credit assignment in complex networks. Journal of The Royal Society Interface, 5(19), 223–235. https://doi.org/10.1098/RSIF.2007.1079

* Update mutual_information_method.jl

* Update dependency to ComplexityMeasures

As `ComplexityMeasures` is being maintained, choose this dependency instead of `InformationMeasures` which has not been updated in 4 years.

* Update src/mutual_information_sensitivity.jl

Co-authored-by: George Datseris <[email protected]>

* Update mutual_information_sensitivity.jl

* Update mutual_information_sensitivity.jl

* Update mutual_information_sensitivity.jl

Update docstring to include additional detail about the discretization entropy.

* Bump version to 2.7

Also apply formatter and update readme

* Simplify and remove higher orders

Second order and total order give some weird results. For now, I've implemented a very basic version of the first order sensitivities using mutual information. The outcomes also correspond to the sobol first order outcomes.

* Update mutual_information_method.jl

* Update mutual_information_method.jl

* Apply formatter

* Update Project.toml

* Update mutual_information_sensitivity.jl

---------

Co-authored-by: George Datseris <[email protected]>

Project.toml

@@ -1,10 +1,11 @@
 name = "GlobalSensitivity"
 uuid = "af5da776-676b-467e-8baf-acd8249e4f0f"
 authors = ["Vaibhavdixit02 <[email protected]>"]
-version = "2.6.2"
+version = "2.7.0"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+ComplexityMeasures = "ab4b797d-85ee-42ba-b621-05d793b346a2"
 Copulas = "ae264745-0b69-425e-9d9d-cf662c5eec93"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
@@ -23,6 +24,7 @@ Trapz = "592b5752-818d-11e9-1e9a-2b8ca4a44cd1"
 [compat]
 Aqua = "0.8"
 Combinatorics = "1"
+ComplexityMeasures = "3.6"
 Copulas = "0.1.22"
 Distributions = "0.25.87"
 FFTW = "1.3"

README.md

@@ -10,7 +10,7 @@
 [![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle)
 [![DOI](https://joss.theoj.org/papers/10.21105/joss.04561/status.svg)](https://doi.org/10.21105/joss.04561)
 
-GlobalSensitivity.jl package contains implementation of some the most popular GSA methods. Currently it supports Delta Moment-Independent, DGSM, EASI, eFAST, Morris, Fractional Factorial, RBD-FAST, Sobol and Regression based sensitivity methods.
+GlobalSensitivity.jl package contains implementation of some the most popular GSA methods. Currently it supports Delta Moment-Independent, DGSM, EASI, eFAST, Morris, Mutual Information, Fractional Factorial, RBD-FAST, RSA, Sobol and Regression based sensitivity methods.
 
 ## Tutorials and Documentation
 

docs/pages.jl

@@ -12,5 +12,6 @@ pages = [
         "methods/easi.md",
         "methods/fractional.md",
         "methods/rbdfast.md",
-        "methods/rsa.md"]
+        "methods/rsa.md",
+        "methods/mutualinformation.md"]
 ]

docs/src/methods/mutualinformation.md

@@ -0,0 +1,5 @@
+# Mutual Information Method
+
+```@docs
+MutualInformation
+```

src/GlobalSensitivity.jl

@@ -5,6 +5,7 @@ using QuasiMonteCarlo, ForwardDiff, KernelDensity, Trapz
 using Parameters: @unpack
 using FFTW, Distributions, StatsBase
 using Copulas, Combinatorics, ThreadsX
+using ComplexityMeasures: entropy, ValueHistogram, StateSpaceSet
 
 abstract type GSAMethod end
 
@@ -19,6 +20,7 @@ include("rbd-fast_sensitivity.jl")
 include("fractional_factorial_sensitivity.jl")
 include("shapley_sensitivity.jl")
 include("rsa_sensitivity.jl")
+include("mutual_information_sensitivity.jl")
 
 """
     gsa(f, method::GSAMethod, param_range; samples, batch=false)
@@ -61,7 +63,7 @@ function gsa(f, method::GSAMethod, param_range; samples, batch = false) end
 export gsa
 
 export Sobol, Morris, RegressionGSA, DGSM, eFAST, DeltaMoment, EASI, FractionalFactorial,
-       RBDFAST, Shapley, RSA
+       RBDFAST, Shapley, RSA, MutualInformation
 # Added for shapley_sensitivity
 
 end # module

src/mutual_information_sensitivity.jl

@@ -0,0 +1,160 @@
+@doc raw"""
+
+MutualInformation(; n_bootstraps = 1000, conf_level = 0.95)
+
+- `n_bootstraps`: Number of bootstraps to be used for estimation of null distribution. Default is `1000`.
+- `conf_level`: Confidence level for the minimum bound estimation. Default is `0.95`.
+
+
+## Method Details
+
+The sensitivity analysis based on mutual information is an alternative approach to sensitivity analysis based on information theoretic measures. In this method,
+the output uncertainty is quantified by the entropy of the output distribution, instead of taking a variance-based approach. The Shannon entropy of the output is 
+given by:
+
+```math
+H(Y) = -\sum_y p(y) \log p(y)
+```
+Where ``p(y)`` is the probability density function of the output ``Y``. By fixing an input ``X_i``, the conditional entropy of the output ``Y`` is given by:
+
+```math
+H(Y|X_i) = -\sum_{x} p(x) H(Y|X_i = x)
+```
+
+The mutual information between the input ``X_i`` and the output ``Y`` is then given by:
+
+```math
+I(X_i;Y) = H(Y) - H(Y|X_i) = H(X) + H(Y) - H(X,Y)
+```
+
+Where ``H(X,Y)`` is the joint entropy of the input and output. The mutual information can be used to calculate the sensitivity indices of the input parameters. 
+
+### Sensitivity Indices
+The sensitivity indices are calculated as the mutual information between the input ``X_i`` and the output ``Y`` where the ``\alpha`` quantile of
+null distribution of the output is subtracted from the mutual information:
+
+```math
+S_{1,i} = I(X_i;Y) - Q(I(X_i; Y_\text{null}), \alpha)
+```
+
+Using mutual information for sensitivity analysis is introduced in Lüdtke et al. (2007)[^1] and also present in Datseris & Parlitz (2022)[^2]. 
+
+## API
+
+    gsa(f, method::MutualInformation, p_range; samples, batch = false)
+
+Returns a `MutualInformationResult` object containing the resulting sensitivity indices for the parameters and the corresponding confidence intervals.
+The `MutualInformationResult` object contains the following fields:
+- `S`: Sensitivity indices.
+- `mutual_information`: Computed mutual information values
+- `bounds`: Computed upper bounds of the null distribution of mutual information.
+
+### Example
+
+```julia
+using GlobalSensitivity
+
+function ishi_batch(X)
+    A= 7
+    B= 0.1
+    @. sin(X[1,:]) + A*sin(X[2,:])^2+ B*X[3,:]^4 *sin(X[1,:])
+end
+function ishi(X)
+    A= 7
+    B= 0.1
+    sin(X[1]) + A*sin(X[2])^2+ B*X[3]^4 *sin(X[1])
+end
+
+lb = -ones(4)*π
+ub = ones(4)*π
+
+res1 = gsa(ishi,MutualInformation(),[[lb[i],ub[i]] for i in 1:4],samples=15000)
+res2 = gsa(ishi_batch,MutualInformation(),[[lb[i],ub[i]] for i in 1:4],samples=15000,batch=true)
+```
+
+### References
+[^1]: Lüdtke, N., Panzeri, S., Brown, M., Broomhead, D. S., Knowles, J., Montemurro, M. A., & Kell, D. B. (2007). Information-theoretic sensitivity analysis: a general method for credit assignment in complex networks. Journal of The Royal Society Interface, 5(19), 223–235.
+[^2]: Datseris, G., & Parlitz, U. (2022). Nonlinear Dynamics, Ch. 7, pg. 105-119.
+"""
+struct MutualInformation <: GSAMethod
+    n_bootstraps::Int
+    conf_level::Real
+end
+
+function MutualInformation(; n_bootstraps = 1000, conf_level = 0.95)
+    MutualInformation(n_bootstraps, conf_level)
+end
+
+struct MutualInformationResult{T}
+    S::T
+    mutual_information::T
+    bounds::T
+end
+
+function _compute_bounds(Xi, Y, conf_level, n_bootstraps)
+
+    # perform permutations of Y and calculate mutual information
+    mi_values = zeros(n_bootstraps)
+    est = ValueHistogram(Int(round(sqrt(length(Y)))))
+    Y_perm = copy(Y)
+    entropy_Xi = entropy(est, Xi)
+    entropy_Y = entropy(est, Y)
+    for i in 1:n_bootstraps
+        shuffle!(Y_perm)
+        mi_values[i] = entropy_Xi + entropy_Y -
+                       entropy(est, StateSpaceSet(Xi, Y_perm))
+    end
+
+    return quantile(mi_values, conf_level)
+end
+
+function _compute_mi(X::AbstractArray, Y::AbstractVector, method::MutualInformation)
+    # K is the number of variables, samples is the number of simulations
+    K = size(X, 1)
+
+    if method.n_bootstraps > size(X, 2)
+        throw(ArgumentError("Number of bootstraps must be less than or equal to the number of samples"))
+    end
+
+    est = ValueHistogram(Int(round(sqrt(size(Y, 1)))))
+    entropy_Y = entropy(est, Y)
+
+    sensitivities = zeros(K)
+    bounds = zeros(K)
+
+    # calculate mutual information
+    @inbounds for i in 1:K
+        Xi = @view X[i, :]
+        sensitivities[i] = entropy(est, Xi) + entropy_Y -
+                           entropy(est, StateSpaceSet(Xi, Y))
+        bounds[i] = _compute_bounds(Xi, Y, method.conf_level, method.n_bootstraps)
+    end
+
+    sensitivities, bounds
+end
+
+function gsa(f, method::MutualInformation, p_range; samples, batch = false)
+    lb = [i[1] for i in p_range]
+    ub = [i[2] for i in p_range]
+
+    X = QuasiMonteCarlo.sample(samples, lb, ub, QuasiMonteCarlo.SobolSample())
+
+    if batch
+        Y = f(X)
+        multioutput = Y isa AbstractMatrix
+    else
+        Y = [f(X[:, j]) for j in axes(X, 2)]
+        multioutput = !(eltype(Y) <: Number)
+        if eltype(Y) <: RecursiveArrayTools.AbstractVectorOfArray
+            y_size = size(Y[1])
+            Y = vec.(Y)
+            desol = true
+        end
+    end
+
+    mutual_information_values, bounds = _compute_mi(X, Y, method)
+
+    mi_sensitivity = max.(0.0, mutual_information_values .- bounds)
+
+    return MutualInformationResult(mi_sensitivity, mutual_information_values, bounds)
+end

test/mutual_information_method.jl

@@ -0,0 +1,64 @@
+using GlobalSensitivity, Test, QuasiMonteCarlo
+
+function ishi_batch(X)
+    A = 7
+    B = 0.1
+    @. sin(X[1, :]) + A * sin(X[2, :])^2 + B * X[3, :]^4 * sin(X[1, :])
+end
+function ishi(X)
+    A = 7
+    B = 0.1
+    sin(X[1]) + A * sin(X[2])^2 + B * X[3]^4 * sin(X[1])
+end
+
+function linear_int(X)
+    A = 7
+    B = 5.0
+    A * X[1] + B * X[2] * X[3]
+end
+
+function linear_batch(X)
+    A = 7
+    B = 0.1
+    @. A * X[1, :] + B * X[2, :]
+end
+function linear(X)
+    A = 7
+    B = 0.1
+    A * X[1] + B * X[2]
+end
+
+lb = -ones(4) * π
+ub = ones(4) * π
+
+res1 = gsa(
+    ishi, MutualInformation(), [[lb[i], ub[i]] for i in 1:4], samples = 10_000)
+
+res2 = gsa(
+    ishi_batch, MutualInformation(), [[lb[i], ub[i]] for i in 1:4], samples = 10_000, batch = true)
+
+res_sobol = gsa(
+    ishi, Sobol(order = [0, 1, 2]), [[lb[i], ub[i]] for i in 1:4], samples = 10_000)
+print(res1.mutual_information)
+
+@test res1.mutual_information≈[
+    0.8149166301300568, 1.1100302491434046, 0.6926030791690287, 0.5325449915119265] atol=1e-3
+@test res2.mutual_information≈[
+    0.8149166301300568, 1.1100302491434046, 0.6926030791690287, 0.5325449915119265] atol=1e-3
+
+@test sortperm(res1.S) == [4, 3, 1, 2]
+@test sortperm(res2.S) == [4, 3, 1, 2]
+
+@test sortperm(res1.S) == sortperm(abs.(res_sobol.S1))
+@test sortperm(res2.S) == sortperm(abs.(res_sobol.S1))
+
+res1 = gsa(
+    linear, MutualInformation(), [[lb[i], ub[i]] for i in 1:4], samples = 10_000)
+res2 = gsa(
+    linear_batch, MutualInformation(), [[lb[i], ub[i]] for i in 1:4], batch = true,
+    samples = 10_000)
+
+@test res1.mutual_information≈[
+    5.413269699769483, 0.5971993581613084, 0.6037656606447346, 0.6470271152889264] atol=1e-3
+@test res2.mutual_information≈[
+    5.413269699769483, 0.5971993581613084, 0.6037656606447346, 0.6470271152889264] atol=1e-3

test/runtests.jl

@@ -15,5 +15,6 @@ const GROUP = get(ENV, "GROUP", "All")
         @time @safetestset "Easi Method" include("easi_method.jl")
         @time @safetestset "Shapley Method" include("shapley_method.jl")
         @time @safetestset "RSA Method" include("rsa_method.jl")
+        @time @safetestset "Mutual Information Method" include("mutual_information_method.jl")
     end
 end