mbaudin47/gist:dfe71c411bf551c47fa0bb4188b7c5ef

## gistfile1.txt
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Create a FunctionalChaosResult with extra-features.
Compute Sobol' first order (closed) and total indices of a group of variables.
Michaël Baudin, 2023
"""
import openturns as ot


class SuperFunctionalChaosResult:
    def __init__(self, polynomialChaosResult, verboseThreshold=0.01):
        """
        Create a FunctionalChaosResult with extra-features.

        Parameters
        ----------
        polynomialChaosResult : ot.PolynomialChaosResult
            The polynomial chaos result.
        verboseThreshold : float, strictly positive
            The threshold value used to print coefficients.

        Returns
        -------
        None.

        """
        self.polynomialChaosResult = polynomialChaosResult
        # Get enumerate function
        basis = polynomialChaosResult.getOrthogonalBasis()
        self.enumerateFunction = basis.getEnumerateFunction()
        self.coefficients = polynomialChaosResult.getCoefficients()
        self.nbCoeffs = self.coefficients.getSize()
        self.indices = polynomialChaosResult.getIndices()
        self.variance = self._computeVariance()
        self.verboseThreshold = verboseThreshold
        return None

    def _computeVariance(self):
        """
        Returns the variance of the chaos.

        This is a Python implementation of
        `ot.FunctionalChaosRandomVector.getCovariance()`.

        Parameters
        ----------

        Returns
        -------
        variance : float
            The variance of the polynomial chaos.
        """
        #
        variance = 0.0
        for i in range(self.nbCoeffs):
            # Loop over the coefficients
            coefi = self.coefficients[i][0]
            j = self.indices[i]
            multiIndex = self.enumerateFunction(j)
            degreei = sum(multiIndex)
            if degreei > 0:
                variance += coefi**2
        return variance

    def getSobolGroupedIndex(self, group, verbose=False):
        """
        Returns the first order indices of a group of variables.

        This is a Python implementation of
        `ot.FunctionalChaosSobolIndices.getSobolGroupedIndex()`.

        Parameters
        ----------
        group : list(int)
            The indices of the variables in the group.
        verbose : bool
            Set to True to print intermediate messages

        Returns
        -------
        S : float
            The first order Sobol" index of the group.
        """
        #
        S = 0.0
        for i in range(self.nbCoeffs):
            # Loop over the coefficients
            j = self.indices[i]
            coefi = self.coefficients[i][0]
            multiIndex = self.enumerateFunction(j)
            status = self._first_indicator(multiIndex, group)
            fractionOfVariance = coefi**2 / self.variance
            if status and verbose and fractionOfVariance > self.verboseThreshold:
                print(
                    "(getSobolGroupedIndex) %s, FoVar=%.4f"
                    % (str(multiIndex), fractionOfVariance)
                )
            if status:
                S += coefi**2
        #
        S = S / self.variance
        return S

    def getSobolGroupedTotalIndex(self, group, verbose=False):
        """
        Returns the total order indices of a group of variables.

        This is a Python implementation of
        `ot.FunctionalChaosSobolIndices.getSobolGroupedTotalIndex()`.

        Parameters
        ----------
        group : list(int)
            The indices of the variables in the group.
        verbose : bool
            Set to True to print intermediate messages

        Returns
        -------
        S : float
            The total order Sobol" index of the group.
        """
        #
        S = 0.0
        for i in range(self.nbCoeffs):
            # Loop over the coefficients
            j = self.indices[i]
            coefi = self.coefficients[i][0]
            multiIndex = self.enumerateFunction(j)
            status = self._total_indicator(multiIndex, group)
            fractionOfVariance = coefi**2 / self.variance
            if status and verbose and fractionOfVariance > self.verboseThreshold:
                print(
                    "(getSobolGroupedTotalIndex) %s, FoVar=%.4f"
                    % (str(multiIndex), fractionOfVariance)
                )
            if status:
                S += coefi**2
        #
        S = S / self.variance
        return S

    def _first_indicator(self, multiIndex, group):
        """
        Return True if the multiindex is OK for first order Sobol" index of the group.

        Returns True if the given multiindice must be taken into
        account for the computation fo the first order Sobol" indices
        of the group.
        The rule is the following:
        if all variables not in the group have a zero degree and
        if at least one variable in the group have a nonzero degree,
        then the multiindice must be taken into account.

        Parameters
        ----------
        multiIndex : list(int)
            The multiindex.
        group : list(int)
            The group

        Returns
        -------
        status : bool
            Set to True if the multiindex is to be included to compute
            the first order Sobol index of the group.

        Example
        -------
        The variables are X0, X1, X2, X3.

        The group is (X1, X2).

        indicator([0,0,0,0],[1,2]) : False (exception)
        indicator([0,1,1,0],[1,2]) : True (because X1 and X2 have a nonzero degree
                                           and X0 and X3 have a zero degree)
        indicator([0,1,2,0],[1,2]) : True
        indicator([0,2,1,0],[1,2]) : True
        indicator([0,1,0,0],[1,2]) : False (because X2 has a zero degree)
        indicator([1,1,1,0],[1,2]) : False (because X0 has a nonzero degree)
        indicator([1,0,0,0],[1,2]) : False
        indicator([1,0,1,0],[1,2]) : False
        indicator([0,0,0,1],[1,2]) : False
        """
        d = len(multiIndex)
        degreei = sum(multiIndex)
        # Skip the constant term
        if degreei == 0:
            status = False
        else:
            # List of variables with nonzero degree in the multindice
            nonzeroIndiceList = []
            for i in range(d):
                if multiIndex[i] > 0:
                    nonzeroIndiceList.append(i)
            # Check that each variable in the nonzero degree list is in the group
            status = True
            for varIndex in nonzeroIndiceList:
                if varIndex not in group:
                    status = False
                    break
        return status

    def _total_indicator(self, multiIndex, group):
        """
        Return True if the multiindex is OK for total order Sobol" index of the group.

        Returns True if the given multiindice must be taken into
        account for the computation fo the total order Sobol" indices
        of the group.
        The rule is the following:
        if any variable within the group has a nonzero degree,
        then the multiindice must be taken into account.

        Parameters
        ----------
        multiIndex : list(int)
            The multiindex.
        group : list(int)
            The group

        Returns
        -------
        status : bool
            Set to True if the multiindex is to be included to compute
            the first order Sobol index of the group.


        Example
        -------
        The variables are X0, X1, X2, X3.

        The group is (X1, X2).

        indicator([1,1,0,0],[1,2]) : True (because X1 has a nonzero degree)
        indicator([1,0,0,0],[1,2]) : False (because X1 and X2 have a zero degree)
        indicator([1,0,1,0],[1,2]) : True (because X2 has a nonzero degree)
        indicator([0,0,0,1],[1,2]) : False (because X1 and X2 have a zero degree)
        """
        degreei = sum(multiIndex)
        # Skip the constant term
        if degreei == 0:
            status = False
        else:
            # Check that any variable in in the group has a nonzero degree
            status = False
            for varIndex in group:
                if multiIndex[varIndex] > 0:
                    status = True
                    break
        return status

    def draw_Sobol_indices(self):
        """
        Draw Sobol indices.

        Returns
        -------
        graph : ot.Graph
            The Sobol" indices.

        """
        metamodel = self.polynomialChaosResult.getMetaModel()
        input_names = metamodel.getInputDescription()
        dimension = metamodel.getInputDimension()
        pce_SA = ot.FunctionalChaosSobolIndices(self.polynomialChaosResult)
        first_order = [pce_SA.getSobolIndex(i) for i in range(dimension)]
        total_order = [pce_SA.getSobolTotalIndex(i) for i in range(dimension)]
        graph = ot.SobolIndicesAlgorithm.DrawSobolIndices(
            input_names, first_order, total_order
        )
        return graph

    def printCoefficientsTable(self):
        """
        Print the coefficients of the polynomial chaos.

        Only coefficients with a part of variance
        greater than verboseThreshold are printed.
        The part of variance of a coefficient a[k] is :

            PoV[k] = a[k] ** 2 / variance(PCE)

        where

        variance(PCE) = a[1] ** 2 + ... + a[number_of_coeffs] ** 2

        that is, the sum of squared coefficients, except the first a[0].

        Parameters
        ----------
        polynomialChaosResult : ot.PolynomialChaosResult
            The polynomial chaos result.
        """
        print(self)
        return

    def __str__(self):
        """
        Print the coefficients of the polynomial chaos.

        Only coefficients with a part of variance
        greater than verboseThreshold are printed.
        The part of variance of a coefficient a[k] is :

            PoV[k] = a[k] ** 2 / variance(PCE)

        where

        variance(PCE) = a[1] ** 2 + ... + a[number_of_coeffs] ** 2

        that is, the sum of squared coefficients, except the first a[0].

        Parameters
        ----------
        polynomialChaosResult : ot.PolynomialChaosResult
            The polynomial chaos result.
        """
        coefficients = self.polynomialChaosResult.getCoefficients()
        basis = self.polynomialChaosResult.getOrthogonalBasis()
        enumerateFunction = basis.getEnumerateFunction()
        indices = self.polynomialChaosResult.getIndices()
        nbCoeffs = indices.getSize()
        s = ""
        s += "Number of coefficients : %d\n" % (nbCoeffs)
        s += "# Indice, Multi-indice, Degree : Value \n"
        for k in range(nbCoeffs):
            absoluteIndex = indices[k]
            multiindex = enumerateFunction(absoluteIndex)
            c = coefficients[k][0]
            fractionOfVariance = c**2 / self.variance
            if fractionOfVariance > self.verboseThreshold:
                s += "#%d (%d), %s : %.4f \n" % (k, absoluteIndex, multiindex, c)
        return s

    def _repr_html_(self):
        """Get HTML representation."""
        coefficients = self.polynomialChaosResult.getCoefficients()
        basis = self.polynomialChaosResult.getOrthogonalBasis()
        enumerateFunction = basis.getEnumerateFunction()
        indices = self.polynomialChaosResult.getIndices()
        nbCoeffs = indices.getSize()
        # Table of attributes
        inputDimension = self.polynomialChaosResult.getMetaModel().getInputDimension()
        outputDimension = self.polynomialChaosResult.getMetaModel().getOutputDimension()
        basisSize = self.polynomialChaosResult.getIndices().getSize()
        relativeErrors = self.polynomialChaosResult.getRelativeErrors()
        residuals = self.polynomialChaosResult.getResiduals()
        # Header
        html = ""
        html += "<ul>\n"
        html += "  <li>Input dimension : %d</li>\n" % (inputDimension)
        html += "  <li>Output dimension : %d</li>\n" % (outputDimension)
        html += "  <li>Basis size : %d</li>\n" % (basisSize)
        html += "  <li>Relative errors : %s</li>\n" % (relativeErrors)
        html += "  <li>Residuals : %s</li>\n" % (residuals)
        html += "</ul>\n"
        # Table of coefficients
        html += "<table>\n"
        # Header
        html += "<tr>\n"
        html += "  <th>Index</th>\n"
        html += "  <th>Global index</th>\n"
        html += "  <th>Multi-index</th>\n"
        html += "  <th>Coefficient</th>\n"
        html += "</tr>\n"
        # Content
        for k in range(nbCoeffs):
            absoluteIndex = indices[k]
            multiindex = enumerateFunction(absoluteIndex)
            html += "<tr>\n"
            html += "  <td>%d</td>\n" % (k)
            html += "  <td>%d</td>\n" % (absoluteIndex)
            html += "  <td>%s</td>\n" % (multiindex)
            html += "  <td>%s</td>\n" % (coefficients[k])
            html += "</tr>\n"
        html += "</table>\n"

        return html
	#!/usr/bin/env python3
	# -- coding: utf-8 --
	"""
	Create a FunctionalChaosResult with extra-features.
	Compute Sobol' first order (closed) and total indices of a group of variables.
	Michaël Baudin, 2023
	"""
	import openturns as ot


	class SuperFunctionalChaosResult:
	def __init__(self, polynomialChaosResult, verboseThreshold=0.01):
	"""
	Create a FunctionalChaosResult with extra-features.

	Parameters
	----------
	polynomialChaosResult : ot.PolynomialChaosResult
	The polynomial chaos result.
	verboseThreshold : float, strictly positive
	The threshold value used to print coefficients.

	Returns
	-------
	None.

	"""
	self.polynomialChaosResult = polynomialChaosResult
	# Get enumerate function
	basis = polynomialChaosResult.getOrthogonalBasis()
	self.enumerateFunction = basis.getEnumerateFunction()
	self.coefficients = polynomialChaosResult.getCoefficients()
	self.nbCoeffs = self.coefficients.getSize()
	self.indices = polynomialChaosResult.getIndices()
	self.variance = self._computeVariance()
	self.verboseThreshold = verboseThreshold
	return None

	def _computeVariance(self):
	"""
	Returns the variance of the chaos.

	This is a Python implementation of
	`ot.FunctionalChaosRandomVector.getCovariance()`.

	Parameters
	----------

	Returns
	-------
	variance : float
	The variance of the polynomial chaos.
	"""
	#
	variance = 0.0
	for i in range(self.nbCoeffs):
	# Loop over the coefficients
	coefi = self.coefficients[i][0]
	j = self.indices[i]
	multiIndex = self.enumerateFunction(j)
	degreei = sum(multiIndex)
	if degreei > 0:
	variance += coefi**2
	return variance

	def getSobolGroupedIndex(self, group, verbose=False):
	"""
	Returns the first order indices of a group of variables.

	This is a Python implementation of
	`ot.FunctionalChaosSobolIndices.getSobolGroupedIndex()`.

	Parameters
	----------
	group : list(int)
	The indices of the variables in the group.
	verbose : bool
	Set to True to print intermediate messages

	Returns
	-------
	S : float
	The first order Sobol" index of the group.
	"""
	#
	S = 0.0
	for i in range(self.nbCoeffs):
	# Loop over the coefficients
	j = self.indices[i]
	coefi = self.coefficients[i][0]
	multiIndex = self.enumerateFunction(j)
	status = self._first_indicator(multiIndex, group)
	fractionOfVariance = coefi**2 / self.variance
	if status and verbose and fractionOfVariance > self.verboseThreshold:
	print(
	"(getSobolGroupedIndex) %s, FoVar=%.4f"
	% (str(multiIndex), fractionOfVariance)
	)
	if status:
	S += coefi**2
	#
	S = S / self.variance
	return S

	def getSobolGroupedTotalIndex(self, group, verbose=False):
	"""
	Returns the total order indices of a group of variables.

	This is a Python implementation of
	`ot.FunctionalChaosSobolIndices.getSobolGroupedTotalIndex()`.

	Parameters
	----------
	group : list(int)
	The indices of the variables in the group.
	verbose : bool
	Set to True to print intermediate messages

	Returns
	-------
	S : float
	The total order Sobol" index of the group.
	"""
	#
	S = 0.0
	for i in range(self.nbCoeffs):
	# Loop over the coefficients
	j = self.indices[i]
	coefi = self.coefficients[i][0]
	multiIndex = self.enumerateFunction(j)
	status = self._total_indicator(multiIndex, group)
	fractionOfVariance = coefi**2 / self.variance
	if status and verbose and fractionOfVariance > self.verboseThreshold:
	print(
	"(getSobolGroupedTotalIndex) %s, FoVar=%.4f"
	% (str(multiIndex), fractionOfVariance)
	)
	if status:
	S += coefi**2
	#
	S = S / self.variance
	return S

	def _first_indicator(self, multiIndex, group):
	"""
	Return True if the multiindex is OK for first order Sobol" index of the group.

	Returns True if the given multiindice must be taken into
	account for the computation fo the first order Sobol" indices
	of the group.
	The rule is the following:
	if all variables not in the group have a zero degree and
	if at least one variable in the group have a nonzero degree,
	then the multiindice must be taken into account.

	Parameters
	----------
	multiIndex : list(int)
	The multiindex.
	group : list(int)
	The group

	Returns
	-------
	status : bool
	Set to True if the multiindex is to be included to compute
	the first order Sobol index of the group.

	Example
	-------
	The variables are X0, X1, X2, X3.

	The group is (X1, X2).

	indicator([0,0,0,0],[1,2]) : False (exception)
	indicator([0,1,1,0],[1,2]) : True (because X1 and X2 have a nonzero degree
	and X0 and X3 have a zero degree)
	indicator([0,1,2,0],[1,2]) : True
	indicator([0,2,1,0],[1,2]) : True
	indicator([0,1,0,0],[1,2]) : False (because X2 has a zero degree)
	indicator([1,1,1,0],[1,2]) : False (because X0 has a nonzero degree)
	indicator([1,0,0,0],[1,2]) : False
	indicator([1,0,1,0],[1,2]) : False
	indicator([0,0,0,1],[1,2]) : False
	"""
	d = len(multiIndex)
	degreei = sum(multiIndex)
	# Skip the constant term
	if degreei == 0:
	status = False
	else:
	# List of variables with nonzero degree in the multindice
	nonzeroIndiceList = []
	for i in range(d):
	if multiIndex[i] > 0:
	nonzeroIndiceList.append(i)
	# Check that each variable in the nonzero degree list is in the group
	status = True
	for varIndex in nonzeroIndiceList:
	if varIndex not in group:
	status = False
	break
	return status

	def _total_indicator(self, multiIndex, group):
	"""
	Return True if the multiindex is OK for total order Sobol" index of the group.

	Returns True if the given multiindice must be taken into
	account for the computation fo the total order Sobol" indices
	of the group.
	The rule is the following:
	if any variable within the group has a nonzero degree,
	then the multiindice must be taken into account.

	Parameters
	----------
	multiIndex : list(int)
	The multiindex.
	group : list(int)
	The group

	Returns
	-------
	status : bool
	Set to True if the multiindex is to be included to compute
	the first order Sobol index of the group.


	Example
	-------
	The variables are X0, X1, X2, X3.

	The group is (X1, X2).

	indicator([1,1,0,0],[1,2]) : True (because X1 has a nonzero degree)
	indicator([1,0,0,0],[1,2]) : False (because X1 and X2 have a zero degree)
	indicator([1,0,1,0],[1,2]) : True (because X2 has a nonzero degree)
	indicator([0,0,0,1],[1,2]) : False (because X1 and X2 have a zero degree)
	"""
	degreei = sum(multiIndex)
	# Skip the constant term
	if degreei == 0:
	status = False
	else:
	# Check that any variable in in the group has a nonzero degree
	status = False
	for varIndex in group:
	if multiIndex[varIndex] > 0:
	status = True
	break
	return status

	def draw_Sobol_indices(self):
	"""
	Draw Sobol indices.

	Returns
	-------
	graph : ot.Graph
	The Sobol" indices.

	"""
	metamodel = self.polynomialChaosResult.getMetaModel()
	input_names = metamodel.getInputDescription()
	dimension = metamodel.getInputDimension()
	pce_SA = ot.FunctionalChaosSobolIndices(self.polynomialChaosResult)
	first_order = [pce_SA.getSobolIndex(i) for i in range(dimension)]
	total_order = [pce_SA.getSobolTotalIndex(i) for i in range(dimension)]
	graph = ot.SobolIndicesAlgorithm.DrawSobolIndices(
	input_names, first_order, total_order
	)
	return graph

	def printCoefficientsTable(self):
	"""
	Print the coefficients of the polynomial chaos.

	Only coefficients with a part of variance
	greater than verboseThreshold are printed.
	The part of variance of a coefficient a[k] is :

	PoV[k] = a[k] ** 2 / variance(PCE)

	where

	variance(PCE) = a[1] 2 + ... + a[number_of_coeffs] 2

	that is, the sum of squared coefficients, except the first a[0].

	Parameters
	----------
	polynomialChaosResult : ot.PolynomialChaosResult
	The polynomial chaos result.
	"""
	print(self)
	return

	def __str__(self):
	"""
	Print the coefficients of the polynomial chaos.

	Only coefficients with a part of variance
	greater than verboseThreshold are printed.
	The part of variance of a coefficient a[k] is :

	PoV[k] = a[k] ** 2 / variance(PCE)

	where

	variance(PCE) = a[1] 2 + ... + a[number_of_coeffs] 2

	that is, the sum of squared coefficients, except the first a[0].

	Parameters
	----------
	polynomialChaosResult : ot.PolynomialChaosResult
	The polynomial chaos result.
	"""
	coefficients = self.polynomialChaosResult.getCoefficients()
	basis = self.polynomialChaosResult.getOrthogonalBasis()
	enumerateFunction = basis.getEnumerateFunction()
	indices = self.polynomialChaosResult.getIndices()
	nbCoeffs = indices.getSize()
	s = ""
	s += "Number of coefficients : %d\n" % (nbCoeffs)
	s += "# Indice, Multi-indice, Degree : Value \n"
	for k in range(nbCoeffs):
	absoluteIndex = indices[k]
	multiindex = enumerateFunction(absoluteIndex)
	c = coefficients[k][0]
	fractionOfVariance = c**2 / self.variance
	if fractionOfVariance > self.verboseThreshold:
	s += "#%d (%d), %s : %.4f \n" % (k, absoluteIndex, multiindex, c)
	return s

	def _repr_html_(self):
	"""Get HTML representation."""
	coefficients = self.polynomialChaosResult.getCoefficients()
	basis = self.polynomialChaosResult.getOrthogonalBasis()
	enumerateFunction = basis.getEnumerateFunction()
	indices = self.polynomialChaosResult.getIndices()
	nbCoeffs = indices.getSize()
	# Table of attributes
	inputDimension = self.polynomialChaosResult.getMetaModel().getInputDimension()
	outputDimension = self.polynomialChaosResult.getMetaModel().getOutputDimension()
	basisSize = self.polynomialChaosResult.getIndices().getSize()
	relativeErrors = self.polynomialChaosResult.getRelativeErrors()
	residuals = self.polynomialChaosResult.getResiduals()
	# Header
	html = ""
	html += "<ul>\n"
	html += " <li>Input dimension : %d</li>\n" % (inputDimension)
	html += " <li>Output dimension : %d</li>\n" % (outputDimension)
	html += " <li>Basis size : %d</li>\n" % (basisSize)
	html += " <li>Relative errors : %s</li>\n" % (relativeErrors)
	html += " <li>Residuals : %s</li>\n" % (residuals)
	html += "</ul>\n"
	# Table of coefficients
	html += "<table>\n"
	# Header
	html += "<tr>\n"
	html += " <th>Index</th>\n"
	html += " <th>Global index</th>\n"
	html += " <th>Multi-index</th>\n"
	html += " <th>Coefficient</th>\n"
	html += "</tr>\n"
	# Content
	for k in range(nbCoeffs):
	absoluteIndex = indices[k]
	multiindex = enumerateFunction(absoluteIndex)
	html += "<tr>\n"
	html += " <td>%d</td>\n" % (k)
	html += " <td>%d</td>\n" % (absoluteIndex)
	html += " <td>%s</td>\n" % (multiindex)
	html += " <td>%s</td>\n" % (coefficients[k])
	html += "</tr>\n"
	html += "</table>\n"

	return html