Estimate the expected value of test functions using quasi-Monte Carlo.

sobol_convergence_with_or_not_zero.py

"""
Estimate the expected value of the test functions using quasi-Monte Carlo.

See the effect of:
- scrambling the Sobol' sequence,
- including or excluding the first Sobol' point (zero).

See:

- The first zero point should be moved back into the low discrepancy sequences
  https://github.com/openturns/openturns/issues/2686

- ENH: add stats.qmc module with quasi Monte Carlo functionality
  https://github.com/scipy/scipy/pull/10844

- Integration convergence using Sobol' sequence: removing the first point
  Pamphile Roy
  https://gist.github.com/tupui/fb6e219b1dd2316b7498ebce231bfff5

References
----------
- Owen, A. B. (2020, August). On dropping the first Sobol' point.
  In International conference on Monte Carlo and quasi-Monte Carlo methods
  in scientific computing (pp. 71-86). Cham: Springer International Publishing.
- Kucherenko and Daniel Albrecht and Andrea Saltelli.
  Exploring multi-dimensional spaces: a Comparison of Latin Hypercube and
  Quasi Monte Carlo Sampling Techniques. arXiv 1505.02350, 2015.
"""

# %%
import openturns as ot
from openturns.usecases import ishigami_function
import openturns.viewer as otv
import time
import numpy as np
from scipy.stats import qmc


# %%
def computeAbsoluteError(
    test_problem,
    sample_size_min=1,
    sample_size_factor=2.0,
    iteration_maximum=100,
    maximum_elapsed_time=2.0,
    verbose=False,
):
    """
    Compute the convergence of absolute errors of E[Y] using MC and QMC sampling.

    This function estimates the absolute error of the expected value of a model
    output using Monte Carlo (MC) and several Quasi-Monte Carlo (QMC) sampling
    methods. It performs iterative experiments with increasing sample sizes until
    a time limit or a maximum iteration count is reached.

    Parameters
    ----------
    test_problem : dict
        Dictionary describing the test problem with the following keys:
            - "model" : callable
                Function mapping input samples to model outputs.
            - "distribution" : openturns.Distribution
                Input distribution of the model.
            - "expectation" : float
                Exact or reference expected value of the model output.
    sample_size_min : int, optional
        Initial number of samples. Must be positive.
    sample_size_factor : float, optional
        Multiplicative factor used to increase the sample size at each iteration.
        Must be greater than 1.0.
    iteration_maximum : int, optional
        Maximum number of iterations. Must be positive.
    maximum_elapsed_time : float, optional
        Maximum allowed total computation time in seconds.
    verbose : bool, optional
        If True, display progress information at each iteration.

    Returns
    -------
    tuple of lists
        A tuple containing six lists:
        - sample_size_list : list of int
            Sample sizes used at each iteration.
        - absolute_error_list_MC : list of float
            Absolute errors computed with standard Monte Carlo sampling.
        - absolute_error_list_sobol_with_zero : list of float
            Absolute errors computed with Sobol' sequence including zero.
        - absolute_error_list_sobol_no_zero : list of float
            Absolute errors computed with Sobol' sequence excluding zero.
        - absolute_error_list_scrambled_sobol_with_zero : list of float
            Absolute errors computed with scrambled Sobol' sequence including zero.
        - absolute_error_list_scrambled_sobol_no_zero : list of float
            Absolute errors computed with scrambled Sobol' sequence excluding zero.

    Raises
    ------
    ValueError
        If any input argument is invalid (for example, non-positive sample size
        or iteration count).

    Examples
    --------
    >>> problem = {
    ...     "model": my_model,
    ...     "distribution": my_distribution,
    ...     "expectation": exact_expectation
    ... }
    >>> results = computeAbsoluteError(problem, verbose=True)
    >>> sizes, err_mc, err_sob0, err_sob1, err_scr0, err_scr1 = results
    >>> len(sizes)
    10
    """

    def computeAbsoluteErrorForProblem(sample_size, get_experiment, test_problem):
        model = test_problem["model"]
        distribution = test_problem["distribution"]
        expectation = test_problem["expectation"]
        inputSample, weight = get_experiment(sample_size, distribution)
        outputSample = model(inputSample)
        estimated_expectation = weight.dot(outputSample.asPoint())
        absolute_error = abs(estimated_expectation - expectation)
        return absolute_error

    #
    sample_size_list = []
    absolute_error_list_MC = []
    absolute_error_list_sobol_with_zero = []
    absolute_error_list_sobol_no_zero = []
    absolute_error_list_scrambled_sobol_with_zero = []
    absolute_error_list_scrambled_sobol_no_zero = []
    sample_size = sample_size_min
    for i in range(iteration_maximum):
        sample_size_list.append(sample_size)
        t1 = time.time()
        #
        absolute_error_MC = computeAbsoluteErrorForProblem(
            sample_size, get_montecarlo_experiment, test_problem
        )
        absolute_error_list_MC.append(absolute_error_MC)
        #
        absolute_error_sobol_with_zero = computeAbsoluteErrorForProblem(
            sample_size, get_sobol_with_zero, test_problem
        )
        absolute_error_list_sobol_with_zero.append(absolute_error_sobol_with_zero)
        #
        absolute_error_sobol_no_zero = computeAbsoluteErrorForProblem(
            sample_size, get_sobol_no_zero, test_problem
        )
        absolute_error_list_sobol_no_zero.append(absolute_error_sobol_no_zero)
        #
        absolute_error_scrambled_sobol_with_zero = computeAbsoluteErrorForProblem(
            sample_size, get_scrambled_sobol_with_zero, test_problem
        )
        absolute_error_list_scrambled_sobol_with_zero.append(absolute_error_scrambled_sobol_with_zero)
        #
        absolute_error_scrambled_sobol_no_zero = computeAbsoluteErrorForProblem(
            sample_size, get_sobol_no_zero, test_problem
        )
        absolute_error_list_scrambled_sobol_no_zero.append(absolute_error_scrambled_sobol_no_zero)
        #
        t2 = time.time()
        elapsed_time = t2 - t1
        if verbose:
            print(
                f"n={sample_size}, "
                f"elapsed={elapsed_time:.2f} (s), "
                f"Err(MC)={absolute_error_MC:.4e}, "
                f"Err(Sobol' with 0)={absolute_error_sobol_with_zero:.4e}, "
                f"Err(Sobol' no 0)={absolute_error_sobol_no_zero:.4e}, "
                f"Err(scrambled Sobol' with 0)={absolute_error_scrambled_sobol_with_zero:.4e}, "
                f"Err(scrambled Sobol' no 0)={absolute_error_scrambled_sobol_no_zero:.4e}"
            )
        if elapsed_time > maximum_elapsed_time:
            break
        # Update the sample size
        sample_size = max(1 + sample_size, int(sample_size_factor * sample_size))
    return (
        sample_size_list,
        absolute_error_list_MC,
        absolute_error_list_sobol_with_zero,
        absolute_error_list_sobol_no_zero,
        absolute_error_list_scrambled_sobol_with_zero,
        absolute_error_list_scrambled_sobol_no_zero,
    )


# %%
def polynomialDataFitting(maximum_degree, x_train, y_train):
    """
    Computes the polynomial curve fitting with given total degree

    This is for learning purposes only: please consider a serious metamodel
    for real applications, e.g. polynomial chaos or kriging.

    Parameters
    ----------
    maximum_degree : int
        The maximum polynomial degree
    x_train : ot.Sample()
        The input training sample
    y_train : ot.Sample()
        The output training sample

    Returns
    -------
    responseSurface : ot.Function
        The polynomial
    basis : ot.Basis()
        The polynomial basis
    myLeastSquares : ot.LinearLeastSquares()
        The linear least squares
    """
    polynomialCollection = [
        "x^%d" % (degree) for degree in range(1, maximum_degree + 1)
    ]
    basis = ot.SymbolicFunction(["x"], polynomialCollection)
    designMatrix = basis(x_train)
    myLeastSquares = ot.LinearLeastSquares(designMatrix, y_train)
    myLeastSquares.run()
    responseSurface = myLeastSquares.getMetaModel()

    return (responseSurface, basis, myLeastSquares)


# %%
def linearSample(xmin, xmax, npoints):
    """
    Returns a sample from linear grid

    Parameters
    ----------
    xmin : float
        The minimum value
    xmax : float
        The maximum value
    npoints : int
        The number of values

    Returns
    -------
    vertices : ot.Sample(npoints, 1)
        The list of values
    """
    step = (xmax - xmin) / (npoints - 1)
    rg = ot.RegularGrid(xmin, step, npoints)
    vertices = rg.getVertices()
    return vertices


# %%
def computeLeastSquaresFit(sample_size_list, absolute_error, maximum_degree=1):
    """
    Compute a least-squares polynomial fit between sample size and absolute error.

    This function filters out non-positive errors, applies a logarithmic
    transformation to both sample size and absolute error, and fits a polynomial
    model of a specified degree using least squares regression. It then computes
    predicted values and the slope of the fitted model.

    Parameters
    ----------
    sample_size_list : list or array-like of float
        List of sample sizes. Must contain strictly positive values.
    absolute_error : list or array-like of float
        Corresponding list of absolute errors. Non-positive values are ignored.
    maximum_degree : int, optional
        Maximum polynomial degree to fit. Must be a non-negative integer.
        Default is 1.

    Returns
    -------
    sampleSizeListForLS : ot.Sample
        Sample of x-values (sample sizes) used for least-squares prediction.
    ypredicted_test : ot.Sample
        Predicted y-values (absolute errors) obtained from the fitted model.
    slope : float
        Estimated slope coefficient of the fitted polynomial in log-log space.

    Raises
    ------
    ValueError
        If all absolute errors are non-positive.
    TypeError
        If inputs are not list-like or contain non-numeric elements.

    Examples
    --------
    >>> sample_size_list = [10, 100, 1000]
    >>> absolute_error = [0.1, 0.03, 0.01]
    >>> sampleSizeListForLS, ypredicted_test, slope = computeLeastSquaresFit(
    ...     sample_size_list, absolute_error, maximum_degree=1
    ... )
    >>> slope
    -0.5
    """
    size = len(absolute_error)
    sample_size_list_filtered = []
    absolute_error_list_filtered = []
    for i in range(size):
        if absolute_error[i] > 0.0:
            sample_size_list_filtered.append(sample_size_list[i])
            absolute_error_list_filtered.append(absolute_error[i])
    # Construit les échantillons
    sample_size_list_filtered = ot.Sample.BuildFromPoint(sample_size_list_filtered)
    absolute_error_list_filtered = ot.Sample.BuildFromPoint(
        absolute_error_list_filtered
    )

    # Ajuste un modèle polynomial de degré maximum_degree

    responseSurface, basis, myLeastSquares = polynomialDataFitting(
        maximum_degree,
        np.log10(sample_size_list_filtered),
        np.log10(absolute_error_list_filtered),
    )
    # Graphics
    coef = myLeastSquares.getLinear()
    number_of_points_fit = 20
    sampleSizeListForLS = linearSample(
        sample_size_list_filtered[0, 0],
        sample_size_list_filtered[-1, 0],
        number_of_points_fit,
    )
    x_test = np.log10(sampleSizeListForLS)
    ypredicted_test = ot.Sample.BuildFromPoint(
        [10 ** v[0] for v in responseSurface(basis(x_test))]
    )
    slope = coef[0, 0]
    return sampleSizeListForLS, ypredicted_test, slope


# %%
def plotAbsoluteErrorCloudAndSlope(
    sample_size_list,
    absolute_error_list_MC,
    absolute_error_list_sobol_with_zero,
    absolute_error_list_sobol_no_zero,
    absolute_error_list_scrambled_sobol_with_zero,
    absolute_error_list_scrambled_sobol_no_zero,
):
    """
    Plot absolute error convergence and least-squares slopes on a log-log graph.

    This function generates a log-log plot showing the convergence of absolute
    errors obtained from Monte Carlo (MC) and various Quasi-Monte Carlo (QMC)
    methods. For each method, it plots the scatter of absolute errors versus
    sample size along with a least-squares fit curve estimating the slope of
    convergence.

    Parameters
    ----------
    sample_size_list : list of int
        Sample sizes used at each iteration.
    absolute_error_list_MC : list of float
        Absolute errors obtained with standard Monte Carlo sampling.
    absolute_error_list_sobol_with_zero : list of float
        Absolute errors obtained with Sobol' sequence including zero.
    absolute_error_list_sobol_no_zero : list of float
        Absolute errors obtained with Sobol' sequence excluding zero.
    absolute_error_list_scrambled_sobol_with_zero : list of float
        Absolute errors obtained with scrambled Sobol' sequence including zero.
    absolute_error_list_scrambled_sobol_no_zero : list of float
        Absolute errors obtained with scrambled Sobol' sequence excluding zero.

    Returns
    -------
    ot.Graph
        OpenTURNS Graph object representing the log-log plot of absolute error
        convergence with least-squares slope estimates for each sampling method.

    Raises
    ------
    ValueError
        If the input lists have inconsistent lengths or contain invalid values.

    Examples
    --------
    >>> graph = plotAbsoluteErrorCloudAndSlope(
    ...     sample_size_list,
    ...     absolute_error_list_MC,
    ...     absolute_error_list_sobol_with_zero,
    ...     absolute_error_list_sobol_no_zero,
    ...     absolute_error_list_scrambled_sobol_with_zero,
    ...     absolute_error_list_scrambled_sobol_no_zero,
    ... )
    >>> view = ot.View(graph)
    >>> view.show()
    """

    def plotCloudAndSlope(
        sample_size_list, absolute_error_list, pointStyle, label, color, curveStyle
    ):
        graph = ot.Graph("", "", "", True)
        sampleSizeListForLS_full, ypredicted_test_full, slope_full = (
            computeLeastSquaresFit(sample_size_list, absolute_error_list)
        )
        cloud = ot.Cloud(sample_size_list, absolute_error_list)
        cloud.setPointStyle(pointStyle)
        cloud.setLegend(f"{label} ({slope_full:.2f})")
        cloud.setColor(color)
        graph.add(cloud)

        # Linear fit to Full
        curve = ot.Curve(sampleSizeListForLS_full, ypredicted_test_full)
        curve.setLineStyle(curveStyle)
        curve.setColor(color)
        graph.add(curve)

        return graph

    palette = ot.DrawableImplementation.BuildDefaultPalette(5)
    graph = ot.Graph("", "", "", True)
    # MC
    pointStyle = "+"
    label = "MC"
    color = palette[0]
    curveStyle = "dashed"
    curve = plotCloudAndSlope(
        sample_size_list, absolute_error_list_MC, pointStyle, label, color, curveStyle
    )
    graph.add(curve)
    # QMC, Sobol' with zero
    pointStyle = "o"
    label = "Sobol' with 0"
    color = palette[1]
    curveStyle = "dashed"
    curve = plotCloudAndSlope(
        sample_size_list,
        absolute_error_list_sobol_with_zero,
        pointStyle,
        label,
        color,
        curveStyle,
    )
    graph.add(curve)
    # QMC, Sobol' no zero
    pointStyle = "v"
    label = "Sobol' no 0"
    color = palette[2]
    curveStyle = "dashed"
    curve = plotCloudAndSlope(
        sample_size_list,
        absolute_error_list_sobol_no_zero,
        pointStyle,
        label,
        color,
        curveStyle,
    )
    graph.add(curve)
    # QMC, scrambled Sobol' with zero
    pointStyle = "^"
    label = "Scr. Sobol' with 0"
    color = palette[3]
    curveStyle = "dashed"
    curve = plotCloudAndSlope(
        sample_size_list,
        absolute_error_list_scrambled_sobol_with_zero,
        pointStyle,
        label,
        color,
        curveStyle,
    )
    graph.add(curve)
    # QMC, scrambled Sobol' no zero
    pointStyle = "x"
    label = "Scr. Sobol' no 0"
    color = palette[4]
    curveStyle = "dashed"
    curve = plotCloudAndSlope(
        sample_size_list,
        absolute_error_list_scrambled_sobol_no_zero,
        pointStyle,
        label,
        color,
        curveStyle,
    )
    graph.add(curve)
    #
    graph.setXTitle("Sample size")
    graph.setYTitle("Absolute error")
    graph.setLogScale(ot.GraphImplementation.LOGXY)
    # Set bounds
    plot_sample_size_min = int(min(sample_size_list) / 2)
    plot_sample_size_max = int(max(sample_size_list) * 2)
    all_absolute_errors = (
        absolute_error_list_MC
        + absolute_error_list_sobol_with_zero
        + absolute_error_list_sobol_no_zero
        + absolute_error_list_scrambled_sobol_with_zero
        + absolute_error_list_scrambled_sobol_no_zero
    )
    plot_error_max = max(all_absolute_errors) * 10.0
    non_zero_errors = [v for v in all_absolute_errors if v != 0.0]
    plot_error_min = min(non_zero_errors) / 10.0
    graph.setBoundingBox(
        ot.Interval(
            [plot_sample_size_min, plot_error_min],
            [plot_sample_size_max, plot_error_max],
        )
    )

    graph.setLegendPosition("upper left")
    graph.setLegendCorner((1.0, 1.0))
    return graph


# %%
def computeAndPlotCoefficientConvergence(
    test_problem,
    sample_size_min=1,
    sample_size_factor=2.0,
    iteration_maximum=100,
    maximum_elapsed_time=1.0,
    verbose=False,
):
    """
    Compute and plot the convergence of a Fourier coefficient using MC and QMC data.

    This function computes the convergence of a Fourier coefficient expectation
    for a given test problem using both Monte Carlo (MC) and Quasi-Monte Carlo
    (QMC) sampling methods. It then generates a log-log plot displaying the
    absolute error convergence and the least-squares slope for each method.

    Parameters
    ----------
    test_problem : dict
        Dictionary describing the test problem with the following keys:
            - "model" : callable
                Function mapping input samples to model outputs.
            - "distribution" : openturns.Distribution
                Input distribution of the model.
            - "expectation" : float
                Exact or reference expected value of the model output.
    sample_size_min : int, optional
        Initial sample size. Must be positive.
    sample_size_factor : float, optional
        Factor used to increase the sample size at each iteration.
    iteration_maximum : int, optional
        Maximum number of iterations. Must be positive.
    maximum_elapsed_time : float, optional
        Maximum allowed computation time in seconds.
    verbose : bool, optional
        If True, print progress information at each iteration.

    Returns
    -------
    ot.Graph
        OpenTURNS Graph object displaying the convergence of absolute errors
        and least-squares slope estimates for MC and QMC methods.

    Raises
    ------
    ValueError
        If any parameter is invalid (for example, non-positive sample size or
        iteration count).

    Examples
    --------
    >>> problem = {
    ...     "model": my_model,
    ...     "distribution": my_distribution,
    ...     "expectation": exact_expectation
    ... }
    >>> graph = computeAndPlotCoefficientConvergence(problem, verbose=True)
    >>> ot.View(graph).show()
    """
    (
        sample_size_list,
        absolute_error_list_MC,
        absolute_error_list_sobol_with_zero,
        absolute_error_list_sobol_no_zero,
        absolute_error_list_scrambled_sobol_with_zero,
        absolute_error_list_scrambled_sobol_no_zero,
    ) = computeAbsoluteError(
        test_problem,
        sample_size_min=sample_size_min,
        sample_size_factor=sample_size_factor,
        iteration_maximum=iteration_maximum,
        maximum_elapsed_time=maximum_elapsed_time,
        verbose=verbose,
    )
    graph = plotAbsoluteErrorCloudAndSlope(
        sample_size_list,
        absolute_error_list_MC,
        absolute_error_list_sobol_with_zero,
        absolute_error_list_sobol_no_zero,
        absolute_error_list_scrambled_sobol_with_zero,
        absolute_error_list_scrambled_sobol_no_zero
    )
    return graph


# %%
def get_montecarlo_experiment(sample_size, distribution):
    """
    Generate a Monte Carlo experiment and return samples with weights.

    This function builds a Monte Carlo experiment based on the provided
    distribution and generates a sample of the given size along with the
    corresponding weights.

    Parameters
    ----------
    sample_size : int
        Number of samples to generate. Must be positive.
    distribution : openturns.Distribution
        Input distribution used for sampling.

    Returns
    -------
    tuple
        A tuple (sample, weight) where:
        - sample : openturns.Sample
            Generated sample of size `sample_size`.
        - weight : openturns.Point
            Corresponding weights for the generated sample.
    """
    experiment = ot.MonteCarloExperiment(distribution, sample_size)
    sample, weight = experiment.generateWithWeights()
    return sample, weight


# %%
def get_sobol_with_zero(sample_size, distribution):
    """
    Generate a Sobol' sequence experiment including the zero point.

    This function constructs a low-discrepancy experiment using a Sobol' sequence
    starting with the zero point. It returns the generated sample and the
    associated weights.

    Parameters
    ----------
    sample_size : int
        Number of samples to generate. Must be positive.
    distribution : openturns.Distribution
        Input distribution used for sampling.

    Returns
    -------
    tuple
        A tuple (sample, weight) where:
        - sample : openturns.Sample
            Generated Sobol' sample including the zero point.
        - weight : openturns.Point
            Corresponding weights for the generated sample.
    """
    ot.ResourceMap.SetAsUnsignedInteger("SobolSequence-InitialSeed", 0)
    experiment = ot.LowDiscrepancyExperiment(
        ot.SobolSequence(), distribution, sample_size, True
    )
    sample, weight = experiment.generateWithWeights()
    return sample, weight


# %%
def get_sobol_no_zero(sample_size, distribution):
    """
    Generate a Sobol' sequence experiment excluding the zero point.

    This function constructs a low-discrepancy experiment using a Sobol' sequence
    that skips the initial zero point. It returns the generated sample and the
    associated weights.

    Parameters
    ----------
    sample_size : int
        Number of samples to generate. Must be positive.
    distribution : openturns.Distribution
        Input distribution used for sampling.

    Returns
    -------
    tuple
        A tuple (sample, weight) where:
        - sample : openturns.Sample
            Generated Sobol' sample excluding the zero point.
        - weight : openturns.Point
            Corresponding weights for the generated sample.
    """
    ot.ResourceMap.SetAsUnsignedInteger("SobolSequence-InitialSeed", 1)
    experiment = ot.LowDiscrepancyExperiment(
        ot.SobolSequence(), distribution, sample_size, True
    )
    sample, weight = experiment.generateWithWeights()
    return sample, weight


# %%
# Sobol' scrambled with zero using Scipy
def get_scrambled_sobol_with_zero(sample_size, distribution):
    """
    Generate a scrambled Sobol' sequence including the zero point using SciPy.

    This function creates a scrambled Sobol' low-discrepancy sequence with the
    initial point included. The sequence is generated with SciPy's `qmc.Sobol`
    engine, then transformed to follow the target distribution using
    OpenTURNS' `DistributionTransformation`.

    Parameters
    ----------
    sample_size : int
        Number of samples to generate. Must be positive.
    distribution : openturns.Distribution
        Target distribution used to transform the unit Sobol' sample.

    Returns
    -------
    tuple
        A tuple (sample, weight) where:
        - sample : openturns.Sample
            Generated scrambled Sobol' sample transformed to the target
            distribution.
        - weight : openturns.Point
            Corresponding uniform weights of size `sample_size`.
    """
    dimension = distribution.getDimension()
    transformation = ot.DistributionTransformation(
        ot.IndependentCopula(dimension), distribution
    )
    engine = qmc.Sobol(d=dimension, scramble=True)
    unit_sample = engine.random(sample_size)
    sample = transformation(unit_sample)
    weight = (1.0 / sample_size) * ot.Point(sample_size, 1.0)
    return sample, weight



# %%
# Sobol' scrambled without zero using Scipy
def get_scrambled_sobol_no_zero(sample_size, distribution):
    """
    Generate a scrambled Sobol' sequence excluding the zero point using SciPy.

    This function creates a scrambled Sobol' low-discrepancy sequence with the
    initial point removed. The sequence is generated using SciPy's `qmc.Sobol`
    engine and transformed to follow the target distribution via OpenTURNS'
    `DistributionTransformation`.

    Parameters
    ----------
    sample_size : int
        Number of samples to generate. Must be positive.
    distribution : openturns.Distribution
        Target distribution used to transform the unit Sobol' sample.

    Returns
    -------
    tuple
        A tuple (sample, weight) where:
        - sample : openturns.Sample
            Generated scrambled Sobol' sample excluding the zero point and
            transformed to the target distribution.
        - weight : openturns.Point
            Corresponding uniform weights of size `sample_size`.
    """
    dimension = distribution.getDimension()
    transformation = ot.DistributionTransformation(
        ot.IndependentCopula(dimension), distribution
    )
    engine = qmc.Sobol(d=dimension, scramble=True)
    unit_sample = engine.random(sample_size + 1)
    unit_sample = unit_sample[1:]  # Remove the first point
    sample = transformation(unit_sample)
    weight = (1.0 / sample_size) * ot.Point(sample_size, 1.0)
    return sample, weight

# %%
# Check sampling
sample_size = 2 ** 3
distribution = ot.JointDistribution([ot.Uniform(0.0, 1.0)] * 5)
print("Monte-Carlo")
sample, weight = get_montecarlo_experiment(sample_size, distribution)
print(weight)
print(sample)
print("Sobol' with 0")
sample, weight = get_sobol_with_zero(sample_size, distribution)
print(weight)
print(sample)
print("Sobol' no 0")
sample, weight = get_sobol_no_zero(sample_size, distribution)
print(weight)
print(sample)
print("Scrambled Sobol' with 0")
sample, weight = get_scrambled_sobol_with_zero(sample_size, distribution)
print(weight)
print(sample)
print("Scrambled Sobol' no 0")
sample, weight = get_scrambled_sobol_no_zero(sample_size, distribution)
print(weight)
print(sample)


# %%
# See the convergence of MC, QMC0 and QMC1 on Ishigami
im = ishigami_function.IshigamiModel()
test_problem_ishigami = {
    "expectation": im.expectation,
    "distribution": im.inputDistribution,
    "model": im.model,
    "name": "Ishigami",
}


# %%
# Sum of 5 centered exps.
# See eq. (1) in (Owen, 2020).
def sumOfExponentials(x):
    s = np.sum(np.exp(x) - np.exp(1) + 1.0)
    return [s]


dimension = 5
model = ot.PythonFunction(dimension, 1, sumOfExponentials)
expectation = 0.0
inputDistribution = ot.JointDistribution([ot.Uniform(0.0, 1.0)] * dimension)
test_problem_sum_exp = {
    "expectation": expectation,
    "distribution": inputDistribution,
    "model": model,
    "name": "Sum of 5 exp",
}


# %%
# Sum of 5 units.
# See eq. (4) in (Owen, 2020).
def sumOfUnits(x):
    s = np.sum(x) ** 2
    return [s]


dimension = 5
model = ot.PythonFunction(dimension, 1, sumOfUnits)
expectation = dimension / 3.0 + dimension * (dimension - 1) / 4
inputDistribution = ot.JointDistribution([ot.Uniform(0.0, 1.0)] * dimension)
test_problem_sum_unit = {
    "expectation": expectation,
    "distribution": inputDistribution,
    "model": model,
    "name": "Sum of 5 units",
}


# %%
# Product of 5 exponentials
# See eq. (4) in (Owen, 2020).
def productOfExponentials(x):
    s = np.prod(np.exp(x) - np.exp(1) + 1.0)
    return [s]


dimension = 5
model = ot.PythonFunction(dimension, 1, productOfExponentials)
expectation = 0.0
inputDistribution = ot.JointDistribution([ot.Uniform(0.0, 1.0)] * dimension)
test_problem_product_exp = {
    "expectation": expectation,
    "distribution": inputDistribution,
    "model": model,
    "name": "Product of 5 exp",
}


# %%
# Type A
# See (Kucherenko et al., 2015).
def typeA(x):
    dimension = len(x)
    a = np.arange(1, dimension + 1)
    f = 1.0
    for i in range(dimension):
        f *= (abs(4.0 * x[i] - 2) + a[i]) / (1.0 + a[i])
    return [f]


dimension = 30
model = ot.PythonFunction(dimension, 1, typeA)
expectation = 1.0
inputDistribution = ot.JointDistribution([ot.Uniform(0.0, 1.0)] * dimension)
test_problem_type_A = {
    "expectation": expectation,
    "distribution": inputDistribution,
    "model": model,
    "name": "Type A",
}


# %%
# Type B
# See (Kucherenko et al., 2015).
def typeB(x):
    dimension = len(x)
    f = 1.0
    for i in range(1, dimension + 1):
        f *= (i - x[i - 1]) / (i - 0.5)
    return [f]


dimension = 30
model = ot.PythonFunction(dimension, 1, typeB)
expectation = 1.0
inputDistribution = ot.JointDistribution([ot.Uniform(0.0, 1.0)] * dimension)
test_problem_type_B = {
    "expectation": expectation,
    "distribution": inputDistribution,
    "model": model,
    "name": "Type B",
}


# %%
# Type C
# See (Kucherenko et al., 2015).
def typeC(x):
    dimension = len(x)
    f = 2**dimension * np.prod(x)
    return [f]


dimension = 30
model = ot.PythonFunction(dimension, 1, typeC)
expectation = 1.0
inputDistribution = ot.JointDistribution([ot.Uniform(0.0, 1.0)] * dimension)
test_problem_type_C = {
    "expectation": expectation,
    "distribution": inputDistribution,
    "model": model,
    "name": "Type C",
}

# %%
# See how influential the LOO is to Monte-Carlo sampling.
maximum_elapsed_time = 8.0
sample_size_min = 8  # Enforce a power 2 sample size
sample_size_factor = 2.0
iteration_maximum = 100
verbose = True

# %%
benchmark = [
    test_problem_ishigami,
    test_problem_sum_exp,
    test_problem_sum_unit,
    test_problem_product_exp,
    test_problem_type_A,
    test_problem_type_B,
    test_problem_type_C,
]
# for test_problem in benchmark:
for test_problem in benchmark:
    name = test_problem["name"]
    print(f"Testing {name}")
    graph = computeAndPlotCoefficientConvergence(
        test_problem,
        sample_size_min=sample_size_min,
        sample_size_factor=sample_size_factor,
        iteration_maximum=iteration_maximum,
        maximum_elapsed_time=maximum_elapsed_time,
        verbose=verbose,
    )
    graph.setTitle(f"{name}, " r"$n=2^e$")
    view = otv.View(graph, figure_kw={"figsize": (7.0, 4.5)})

# %%

Here are the results.