.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/plot_multifidelity_quadrature.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_plot_multifidelity_quadrature.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_plot_multifidelity_quadrature.py:


Multi-fidelity Quadrature
-------------------------
The following provides an example of how to use multi-fidelity quadrature, e.g. multilevel Monte Carlo, control variates to estimate the mean of a high-fidelity model from an ensemble of related models of varying cost and accuracy. A set of detailed tutorials on this subject can be found in the tutorials section, e.g. :ref:`sphx_glr_auto_tutorials_multi_fidelity_plot_monte_carlo.py`.

.. GENERATED FROM PYTHON SOURCE LINES 7-8

Load the necessary modules

.. GENERATED FROM PYTHON SOURCE LINES 8-15

.. code-block:: default

    import numpy as np
    from pyapprox.benchmarks.benchmarks import setup_benchmark
    from pyapprox import multifidelity
    from pyapprox import interface
    # set seed for reproducibility
    np.random.seed(1)








.. GENERATED FROM PYTHON SOURCE LINES 16-17

First define an ensemble of models using :py:class:`~pyapprox.benchmarks.setup_benchmark`, see  :mod:`pyapprox.benchmarks`.

.. GENERATED FROM PYTHON SOURCE LINES 17-22

.. code-block:: default

    benchmark = setup_benchmark(
        "tunable_model_ensemble", theta1=np.pi/2*.95, shifts=[.1, .2])
    model_ensemble = interface.ModelEnsemble(benchmark.funs)
    hf_mean = benchmark.mean[0]








.. GENERATED FROM PYTHON SOURCE LINES 23-24

Initialize a multifidelity estimator. This requires an estimate of the covariance between the models and the model costs and the random variable representing the model inputs

.. GENERATED FROM PYTHON SOURCE LINES 24-40

.. code-block:: default


    # generate pilot samples to estimate correlation
    npilot_samples = int(1e2)
    # The models are trivial to evaluate so make up model costs
    model_costs = 10.**(-np.arange(3))


    stat_name = "mean"
    cov = multifidelity.estimate_model_ensemble_covariance(
        npilot_samples, benchmark.variable.rvs, model_ensemble,
        model_ensemble.nmodels)[0]
    stat = multifidelity.multioutput_stats[stat_name](benchmark.nqoi)
    stat.set_pilot_quantities(cov)
    est_name = "mlblue"
    est = multifidelity.get_estimator(est_name, stat, model_costs)








.. GENERATED FROM PYTHON SOURCE LINES 41-42

Define a target cost and determine the optimal number of samples to allocate to each model

.. GENERATED FROM PYTHON SOURCE LINES 42-57

.. code-block:: default

    target_cost = 1000
    est.allocate_samples(target_cost)
    args = [benchmark.variable] if est_name == "mlblue" else []
    samples_per_model = est.generate_samples_per_model(
        benchmark.variable.rvs)
    best_models = [benchmark.funs[idx] for idx in est._best_model_indices]
    values_per_model = [
        fun(samples) for fun, samples in zip(best_models, samples_per_model)]
    mf_mean = est(values_per_model)

    print("Multi-fidelity mean", mf_mean)
    print("Exact high-fidelity mean", hf_mean)
    print("Multi-fidelity estimator variance",
          est._covariance_from_npartition_samples(est._rounded_npartition_samples))





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Multi-fidelity mean tensor(-0.0121)
    Exact high-fidelity mean [0.]
    Multi-fidelity estimator variance tensor([[0.0004]])




.. GENERATED FROM PYTHON SOURCE LINES 58-67

Excercises
^^^^^^^^^^
Compare the multi-fidelity mean to the single-fidelity means using only one model

Increase the target cost

Change the correlation between the models by varying the theta1 argument to setup benchmarks

Change the estimator (via est_name). Names of the available estimators can be printed via

.. GENERATED FROM PYTHON SOURCE LINES 67-73

.. code-block:: default


    print(multifidelity.factory.multioutput_estimators.keys())

    #Change the statistic computed (via stat_name). Names of the implemented statistics can be printed via

    print(multifidelity.factory.multioutput_stats.keys())




.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    dict_keys(['cv', 'gmf', 'gis', 'grd', 'mfmc', 'mlmc', 'mc', 'mlblue'])
    dict_keys(['mean', 'variance', 'mean_variance'])





.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 0 minutes  0.013 seconds)


.. _sphx_glr_download_auto_examples_plot_multifidelity_quadrature.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example




    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_multifidelity_quadrature.py <plot_multifidelity_quadrature.py>`

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_multifidelity_quadrature.ipynb <plot_multifidelity_quadrature.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_