benchmark_hsic.py

import openturns as ot
from openturns.usecases.wingweight_function import WingWeightModel
import time
import pandas as pd

ot.Log.Show(ot.Log.NONE)
m = WingWeightModel()

inputNames = m.distributionX.getDescription()

stat = "u-stat" # or v-stat
if stat == "u-stat":
    estimatorType = ot.HSICUStat()
    csv_file = "results_ustat_branche.csv"
else:
    estimatorType = ot.HSICVStat()
    csv_file = "results_vstat_branche.csv"

results = {}

# Generating once the models
covarianceModelCollection = []
for i in range(m.dim):
    inputCovariance = ot.SquaredExponential(1)
    covarianceModelCollection.append(inputCovariance)

outputCovariance = ot.SquaredExponential(1)
covarianceModelCollection.append(outputCovariance)

for sizeHSIC in [100, 200, 500, 1000, 2000, 2500]:

    size_results = []
    ot.RandomGenerator.SetSeed(0)
    inputDesignHSIC = m.distributionX.getSample(sizeHSIC)
    outputDesignHSIC = m.model(inputDesignHSIC)

    for i in range(m.dim):
        Xi = inputDesignHSIC.getMarginal(i)
        covarianceModelCollection[i].setScale(Xi.computeStandardDeviation())

    # We define a covariance kernel associated to the output variable.
    covarianceModelCollection[m.dim].setScale(outputDesignHSIC.computeStandardDeviation())

    # We now build the HSIC estimator:
    print("Size = ", sizeHSIC)
    tic = time.time()
    globHSIC = ot.HSICEstimatorGlobalSensitivity(
        covarianceModelCollection, inputDesignHSIC, outputDesignHSIC, estimatorType
    )
    toc = time.time()
    print("Instanciation time = ", toc - tic)
    size_results.append(toc - tic)
    # We get the R2-HSIC indices:
    tic = time.time()
    R2HSICIndices = globHSIC.getR2HSICIndices()
    toc = time.time()
    size_results.append(toc - tic)
    print("R2-HSIC Indices: ", R2HSICIndices)
    print("Time for R2 estimate = ", toc - tic)

    # and the HSIC indices:
    tic = time.time()
    HSICIndices = globHSIC.getHSICIndices()
    toc = time.time()
    size_results.append(toc - tic)

    print("HSIC Indices: ", HSICIndices)
    print("Elapsed time : ", toc - tic)

    # We have an asymptotic estimate of the value for this estimator.
    tic = time.time()
    pvas = globHSIC.getPValuesAsymptotic()
    toc = time.time()
    size_results.append(toc - tic)
    print("p-value (asymptotic): ", pvas)
    print("Elapsed time : ", toc - tic)

    # The p-value by permutation.
    if sizeHSIC <= 2500:
        tic = time.time()
        pvperm = globHSIC.getPValuesPermutation()
        toc = time.time()
        print("p-value (permutation): ", pvperm)
        print("Elapsed time : ", toc - tic)
        size_results.append(toc - tic)
    else:
        size_results.append(1e6)
    print("\n   ")
    results[sizeHSIC] = size_results

results = pd.DataFrame.from_dict(results)
results.index = ["init", "r2-hsic", "hsic", "pval-asy", "pval-perm"]
results = results.T
results.to_csv(csv_file)
print(results.to_markdown())
import openturns as ot
from openturns.usecases.wingweight_function import WingWeightModel
import time
import pandas as pd

ot.Log.Show(ot.Log.NONE)
m = WingWeightModel()

inputNames = m.distributionX.getDescription()

stat = "u-stat" # or v-stat
if stat == "u-stat":
    estimatorType = ot.HSICUStat()
    csv_file = "results_ustat_branche.csv"
else:
    estimatorType = ot.HSICVStat()
    csv_file = "results_vstat_branche.csv"

results = {}

# Generating once the models
covarianceModelCollection = []
for i in range(m.dim):
    inputCovariance = ot.SquaredExponential(1)
    covarianceModelCollection.append(inputCovariance)

outputCovariance = ot.SquaredExponential(1)
covarianceModelCollection.append(outputCovariance)

for sizeHSIC in [100, 200, 500, 1000, 2000, 2500]:

    size_results = []
    ot.RandomGenerator.SetSeed(0)
    inputDesignHSIC = m.distributionX.getSample(sizeHSIC)
    outputDesignHSIC = m.model(inputDesignHSIC)

    for i in range(m.dim):
        Xi = inputDesignHSIC.getMarginal(i)
        covarianceModelCollection[i].setScale(Xi.computeStandardDeviation())

    # We define a covariance kernel associated to the output variable.
    covarianceModelCollection[m.dim].setScale(outputDesignHSIC.computeStandardDeviation())

    # We now build the HSIC estimator:
    print("Size = ", sizeHSIC)
    tic = time.time()
    globHSIC = ot.HSICEstimatorGlobalSensitivity(
        covarianceModelCollection, inputDesignHSIC, outputDesignHSIC, estimatorType
    )
    toc = time.time()
    print("Instanciation time = ", toc - tic)
    size_results.append(toc - tic)
    # We get the R2-HSIC indices:
    tic = time.time()
    R2HSICIndices = globHSIC.getR2HSICIndices()
    toc = time.time()
    size_results.append(toc - tic)
    print("R2-HSIC Indices: ", R2HSICIndices)
    print("Time for R2 estimate = ", toc - tic)

    # and the HSIC indices:
    tic = time.time()
    HSICIndices = globHSIC.getHSICIndices()
    toc = time.time()
    size_results.append(toc - tic)

    print("HSIC Indices: ", HSICIndices)
    print("Elapsed time : ", toc - tic)

    # We have an asymptotic estimate of the value for this estimator.
    tic = time.time()
    pvas = globHSIC.getPValuesAsymptotic()
    toc = time.time()
    size_results.append(toc - tic)
    print("p-value (asymptotic): ", pvas)
    print("Elapsed time : ", toc - tic)

    # The p-value by permutation.
    if sizeHSIC <= 2500:
        tic = time.time()
        pvperm = globHSIC.getPValuesPermutation()
        toc = time.time()
        print("p-value (permutation): ", pvperm)
        print("Elapsed time : ", toc - tic)
        size_results.append(toc - tic)
    else:
        size_results.append(1e6)
    print("\n   ")
    results[sizeHSIC] = size_results

results = pd.DataFrame.from_dict(results)
results.index = ["init", "r2-hsic", "hsic", "pval-asy", "pval-perm"]
results = results.T
results.to_csv(csv_file)
print(results.to_markdown())