"""
Example of sensitivity analyses on the wing weight model
=========================================================
"""

# %%
#
# This example is a brief overview of the use of the most usual sensitivity analysis techniques and how to call them:
#
# - PCC: Partial Correlation Coefficients
# - PRCC: Partial Rank Correlation Coefficients
# - SRC: Standard Regression Coefficients
# - SRRC: Standard Rank Regression Coefficients
# - Pearson coefficients
# - Spearman coefficients
# - Taylor expansion importance factors
# - Sobol' indices
# - Rank-based estimation of Sobol' indices
# - HSIC : Hilbert-Schmidt Independence Criterion
#
# We present the methods on the :ref:`WingWeight function<use-case-wingweight>` and use the same notations.

# %%
# Definition of the model
# -----------------------
#
# We load the model from the usecases module.
#
#
import openturns as ot
import openturns.viewer as otv
from openturns.usecases.wingweight_function import WingWeightModel

m = WingWeightModel()

# %%
# Cross cuts of the function
# --------------------------
#
# Let's have a look on 2D cross cuts of the wing weight function.
# For each 2D cross cut, the other variables are fixed to the input distribution mean values.
# This graph allows one to have a first idea of the variations of the function in pair of dimensions.
# The colors of each contour plot are comparable.

lowerBound = m.distribution.getRange().getLowerBound()
upperBound = m.distribution.getRange().getUpperBound()

nX = ot.ResourceMap.GetAsUnsignedInteger("Evaluation-DefaultPointNumber")
description = m.distribution.getDescription()
description.add("")
m.model.setDescription(description)
m.model.setName("wing weight model")
grid = m.model.drawCrossCuts(
    m.distribution.getMean(),
    lowerBound,
    upperBound,
    [nX] * m.model.getInputDimension(),
    False,
    True,
    176.0,
    363.0,
)
grid.setTitle("")
# Get View object to manipulate the underlying figure
# Here we decide the colormap and the number of levels used for all contours
view = otv.View(grid, contour_kw={"cmap": "hsv", "levels": 55})

axes = view.getAxes()
fig = view.getFigure()
fig.set_size_inches(12, 12)  # reduce the size

# Setup a large colorbar
fig.colorbar(
    view.getSubviews()[1][0].getContourSets()[0], ax=axes[:-2, -1], fraction=0.3
)
# Hide unwanted axes labels
for i in range(len(axes)):
    for j in range(i + 1):
        if i < len(axes) - 1:
            axes[i][j].xaxis.set_ticklabels([])
        if j > 0:
            axes[i][j].yaxis.set_ticklabels([])
fig.subplots_adjust(top=0.99, bottom=0.05, left=0.06, right=0.99)

# %%
# We can see that the variables :math:`t_c, N_z, A, W_{dg}` seem to be influent on the wing weight whereas :math:`\Lambda, \ell, q, W_p, W_{fw}` have less influence on the function.

# %%
# Data generation
# ---------------
#
# We create the input and output data for the estimation of the different sensitivity coefficients and we get the input variables description:

inputNames = m.distribution.getDescription()

size = 500
inputDesign = m.distribution.getSample(size)
outputDesign = m.model(inputDesign)

# %%
# Let's estimate the PCC, PRCC, SRC, SRRC, Pearson and Spearman coefficients, display and analyze them.
# We create a :class:`~openturns.CorrelationAnalysis` model.

corr_analysis = ot.CorrelationAnalysis(inputDesign, outputDesign)

# %%
# PCC coefficients
# ----------------
# We compute here PCC coefficients using the :class:`~openturns.CorrelationAnalysis`.

# %%
pcc_indices = corr_analysis.computePCC()
print(pcc_indices)

# %%
#

# %%
graph = ot.SobolIndicesAlgorithm.DrawCorrelationCoefficients(
    pcc_indices, inputNames, "PCC coefficients - Wing weight"
)
view = otv.View(graph)

# %%
# PRCC coefficients
# -----------------
# We compute here PRCC coefficients using the :class:`~openturns.CorrelationAnalysis`.

# %%
prcc_indices = corr_analysis.computePRCC()
print(prcc_indices)

# %%
graph = ot.SobolIndicesAlgorithm.DrawCorrelationCoefficients(
    prcc_indices, inputNames, "PRCC coefficients - Wing weight"
)
view = otv.View(graph)

# %%
# SRC coefficients
# -------------------
# We compute here SRC coefficients using the :class:`~openturns.CorrelationAnalysis`.

# %%
src_indices = corr_analysis.computeSRC()
print(src_indices)

# %%
graph = ot.SobolIndicesAlgorithm.DrawCorrelationCoefficients(
    src_indices, inputNames, "SRC coefficients - Wing weight"
)
view = otv.View(graph)

# %%
# Normalized squared SRC coefficients (coefficients are made to sum to 1) :

# %%
squared_src_indices = corr_analysis.computeSquaredSRC(True)
print(squared_src_indices)

# %%
# And their associated graph:

# %%
graph = ot.SobolIndicesAlgorithm.DrawCorrelationCoefficients(
    squared_src_indices, inputNames, "Squared SRC coefficients - Wing weight"
)
view = otv.View(graph)

# %%
#

# %%
# SRRC coefficients
# --------------------
# We compute here SRRC coefficients using the :class:`~openturns.CorrelationAnalysis`.

# %%
srrc_indices = corr_analysis.computeSRRC()
print(srrc_indices)

# %%
graph = ot.SobolIndicesAlgorithm.DrawCorrelationCoefficients(
    srrc_indices, inputNames, "SRRC coefficients - Wing weight"
)
view = otv.View(graph)

# %%
# Pearson coefficients
# -----------------------
# We compute here the Pearson :math:`\rho` coefficients using the :class:`~openturns.CorrelationAnalysis`.

# %%
pearson_correlation = corr_analysis.computeLinearCorrelation()
print(pearson_correlation)

# %%
title_pearson_graph = "Pearson correlation coefficients - Wing weight"
graph = ot.SobolIndicesAlgorithm.DrawCorrelationCoefficients(
    pearson_correlation, inputNames, title_pearson_graph
)
view = otv.View(graph)

# %%
# Spearman coefficients
# -----------------------
# We compute here the Spearman :math:`\rho_s` coefficients using the :class:`~openturns.CorrelationAnalysis`.

# %%
spearman_correlation = corr_analysis.computeSpearmanCorrelation()
print(spearman_correlation)

# %%
title_spearman_graph = "Spearman correlation coefficients - Wing weight"
graph = ot.SobolIndicesAlgorithm.DrawCorrelationCoefficients(
    spearman_correlation, inputNames, title_spearman_graph
)
view = otv.View(graph)

# %%
#
# The different computed correlation estimators show that the variables :math:`S_w, A, N_z, t_c` seem to be the most correlated with the wing weight in absolute value.
# Pearson and Spearman coefficients do not reveal any linear nor monotonic correlation as no coefficients are equal to +/- 1.
# Coefficients about :math:`t_c` are negative revealing a negative correlation with the wing weight, that is consistent with the model expression.


# %%
# Taylor expansion importance factors
# -----------------------------------
# We compute here the Taylor expansion importance factors using :class:`~openturns.TaylorExpansionMoments`.

# %%

# %%
# We create a distribution-based RandomVector.
X = ot.RandomVector(m.distribution)

# %%
# We create a composite RandomVector Y from X and m.model.
Y = ot.CompositeRandomVector(m.model, X)

# %%
# We create a Taylor expansion method to approximate moments.
taylor = ot.TaylorExpansionMoments(Y)

# %%
# We get the importance factors.
print(taylor.getImportanceFactors())

# %%
# We draw the importance factors
graph = taylor.drawImportanceFactors()
graph.setTitle("Taylor expansion imporfance factors - Wing weight")
view = otv.View(graph)

# %%
#
# The Taylor expansion importance factors is consistent with the previous estimators as :math:`S_w, A, N_z, t_c` seem to be the most influent variables.
# To analyze the relevance of the previous indices, a Sobol' analysis is now carried out.

# %%
# Sobol' indices
# --------------
# We compute the Sobol' indices from both sampling approach and  Polynomial Chaos Expansion.

# %%
sizeSobol = 1000
sie = ot.SobolIndicesExperiment(m.distribution, sizeSobol)
inputDesignSobol = sie.generate()
inputNames = m.distribution.getDescription()
inputDesignSobol.setDescription(inputNames)
inputDesignSobol.getSize()

# %%
# We see that 12000 function evaluations are required to estimate the first order and total Sobol' indices.

# %%
# Then, we evaluate the outputs corresponding to this design of experiments.

# %%
outputDesignSobol = m.model(inputDesignSobol)

# %%
# We estimate the Sobol' indices with the :class:`~openturns.SaltelliSensitivityAlgorithm`.

# %%
sensitivityAnalysis = ot.SaltelliSensitivityAlgorithm(
    inputDesignSobol, outputDesignSobol, sizeSobol
)

# %%
# The `getFirstOrderIndices` and `getTotalOrderIndices` methods respectively return estimates of all first order and total Sobol' indices.

# %%
print("First order indices:", sensitivityAnalysis.getFirstOrderIndices())

# %%
print("Total order indices:", sensitivityAnalysis.getTotalOrderIndices())


# %%
# The `draw` method produces the following graph. The vertical bars represent the 95% confidence intervals of the estimates.

# %%
graph = sensitivityAnalysis.draw()
graph.setTitle("Sobol indices with Saltelli - wing weight")
view = otv.View(graph)

# %%
# We see that several Sobol' indices are negative, that is inconsistent with the theory. Therefore, a larger number of samples is required to get consistent indices
sizeSobol = 10000
sie = ot.SobolIndicesExperiment(m.distribution, sizeSobol)
inputDesignSobol = sie.generate()
inputNames = m.distribution.getDescription()
inputDesignSobol.setDescription(inputNames)
inputDesignSobol.getSize()
outputDesignSobol = m.model(inputDesignSobol)

sensitivityAnalysis = ot.SaltelliSensitivityAlgorithm(
    inputDesignSobol, outputDesignSobol, sizeSobol
)

sensitivityAnalysis.getFirstOrderIndices()
sensitivityAnalysis.getTotalOrderIndices()

graph = sensitivityAnalysis.draw()
graph.setTitle("Sobol indices with Saltelli - wing weight")
view = otv.View(graph)

# %%
# It improves the accuracy of the estimation but, for very low indices, Saltelli scheme is not accurate since several confidence intervals provide negative lower bounds.

# %%
# Now, we estimate the Sobol' indices using Polynomial Chaos Expansion.
# We create a Functional Chaos Expansion.
sizePCE = 800
inputDesignPCE = m.distribution.getSample(sizePCE)
outputDesignPCE = m.model(inputDesignPCE)

algo = ot.FunctionalChaosAlgorithm(inputDesignPCE, outputDesignPCE, m.distribution)

algo.run()
result = algo.getResult()

# %%
# Then, we exploit the surrogate model to compute the Sobol' indices.
sensitivityAnalysis = ot.FunctionalChaosSobolIndices(result)
sensitivityAnalysis

# %%
firstOrder = [sensitivityAnalysis.getSobolIndex(i) for i in range(m.dim)]
totalOrder = [sensitivityAnalysis.getSobolTotalIndex(i) for i in range(m.dim)]
graph = ot.SobolIndicesAlgorithm.DrawSobolIndices(inputNames, firstOrder, totalOrder)
graph.setTitle("Sobol indices by Polynomial Chaos Expansion - wing weight")
view = otv.View(graph)


# %%
# Furthermore, first order Sobol' indices can also been estimated in a data-driven way using a rank-based sensitivity algorithm.
# In such a way, the estimation of sensitivity indices does not involve any surrogate model.
sizeRankSobol = 800
inputDesignRankSobol = m.distribution.getSample(sizeRankSobol)
outputDesignankSobol = m.model(inputDesignRankSobol)
myRankSobol = ot.RankSobolSensitivityAlgorithm(
    inputDesignRankSobol, outputDesignankSobol
)
indicesrankSobol = myRankSobol.getFirstOrderIndices()
print("First order indices:", indicesrankSobol)
graph = myRankSobol.draw()
graph.setTitle("Sobol indices by rank-based estimation - wing weight")
view = otv.View(graph)

# %%
#
# The Sobol' indices confirm the previous analyses, in terms of ranking of the most influent variables.
# We also see that five variables have a quasi null total Sobol' indices, that indicates almost no influence on the wing weight.
# There is no discrepancy between first order and total Sobol' indices, that indicates no or very low interaction between the variables in the variance of the output.
# As the most important variables act only through decoupled first degree contributions, the hypothesis of a linear dependence between the input variables and the weight is legitimate.
# This explains why both squared SRC and Taylor give the exact same results even if the first one is based on a :math:`\mathcal{L}^2` linear approximation
# and the second one is based on a linear expansion around the mean value of the input variables.


# %%
# HSIC indices
# ------------

# %%
# We then estimate the HSIC indices using a data-driven approach.
sizeHSIC = 250
inputDesignHSIC = m.distribution.getSample(sizeHSIC)
outputDesignHSIC = m.model(inputDesignHSIC)

covarianceModelCollection = []

# %%
for i in range(m.dim):
    Xi = inputDesignHSIC.getMarginal(i)
    inputCovariance = ot.SquaredExponential(1)
    inputCovariance.setScale(Xi.computeStandardDeviation())
    covarianceModelCollection.append(inputCovariance)

# %%
# We define a covariance kernel associated to the output variable.
outputCovariance = ot.SquaredExponential(1)
outputCovariance.setScale(outputDesignHSIC.computeStandardDeviation())
covarianceModelCollection.append(outputCovariance)

# %%
# In this paragraph, we perform the analysis on the raw data: that is
# the global HSIC estimator.
estimatorType = ot.HSICUStat()

# %%
# We now build the HSIC estimator:
globHSIC = ot.HSICEstimatorGlobalSensitivity(
    covarianceModelCollection, inputDesignHSIC, outputDesignHSIC, estimatorType
)

# %%
# We get the R2-HSIC indices:
R2HSICIndices = globHSIC.getR2HSICIndices()
print("\n Global HSIC analysis")
print("R2-HSIC Indices: ", R2HSICIndices)

# %%
# and the HSIC indices:
HSICIndices = globHSIC.getHSICIndices()
print("HSIC Indices: ", HSICIndices)

# %%
# The p-value by permutation.
pvperm = globHSIC.getPValuesPermutation()
print("p-value (permutation): ", pvperm)

# %%
# We have an asymptotic estimate of the value for this estimator.
pvas = globHSIC.getPValuesAsymptotic()
print("p-value (asymptotic): ", pvas)

# %%
# We vizualise the results.
graph1 = globHSIC.drawHSICIndices()
view1 = otv.View(graph1)

graph2 = globHSIC.drawPValuesAsymptotic()
view2 = otv.View(graph2)

graph3 = globHSIC.drawR2HSICIndices()
view3 = otv.View(graph3)

graph4 = globHSIC.drawPValuesPermutation()
view4 = otv.View(graph4)

# %%
# The HSIC indices go in the same way as the other estimators in terms the most influent variables.
# The variables :math:`W_{fw}, q, l, W_p` seem to be independent to the output as the corresponding p-values are high.
# We can also see that the asymptotic p-values and p-values estimated by permutation are quite similar.

# %%
otv.View.ShowAll()