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

.. only:: html

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

        Click :ref:`here <sphx_glr_download_auto_examples_emulator_plot_ppgasp.py>`
        to download the full example code

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

.. _sphx_glr_auto_examples_emulator_plot_ppgasp.py:


PPGaSP: GP emulation for multi-output functions
===============================================

.. GENERATED FROM PYTHON SOURCE LINES 9-23

This example shows how to apply Gaussian process emulation to a multi-output
function using class :class:`.PPGaSP`.

The multi-output function that we are going to look at is the `DIAMOND`
(diplomatic and military operations in a non-warfighting domain) computer model.
It is used as a testbed to illustrate the `PP GaSP` emulator in the `R` package
`RobustGaSP`, see :cite:t:`Gu2019` for more detail.

The simulator has :math:`13` input parameters and :math:`5` outputs. Namely,
:math:`\mathbf{y}=f(\mathbf{x})` where :math:`\mathbf{x}=(x_1,\ldots,x_{13})^T`
and :math:`\mathbf{y}=(y_1,\ldots,y_5)^T`.

The training and testing data are provided in the folder '.../tests/data/'.
We first load the training and testing data. 

.. GENERATED FROM PYTHON SOURCE LINES 24-40

.. code-block:: default


    import numpy as np
    import os

    dir_data = os.path.abspath('../../../tests/data/')

    humanityX = np.genfromtxt(os.path.join(dir_data, 'humanityX.csv'), delimiter=',')
    humanityY = np.genfromtxt(os.path.join(dir_data, 'humanityY.csv'), delimiter=',')
    print(f"Number of training data points: ", humanityX.shape[0])
    print(f"Input dimension: ", humanityX.shape[1])
    print(f"Output dimension: ", humanityY.shape[1])

    humanityXt = np.genfromtxt(os.path.join(dir_data, 'humanityXt.csv'), delimiter=',')
    humanityYt = np.genfromtxt(os.path.join(dir_data, 'humanityYt.csv'), delimiter=',')
    print(f"Number of testing data points: ", humanityXt.shape[0])





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

 .. code-block:: none

    Number of training data points:  120
    Input dimension:  13
    Output dimension:  5
    Number of testing data points:  120




.. GENERATED FROM PYTHON SOURCE LINES 41-43

.. note:: You may need to modify ``dir_data`` according to where you save them
   on your local machine.

.. GENERATED FROM PYTHON SOURCE LINES 45-50

``humanityX`` and ``humanitY`` are the training data, corresponding to ``design``
and ``response`` respectively. ``humanityXt`` are testing input data, at which
we are going to make predictions once the emulator is trained. ``humanityYt``
are the true outputs at ``humanityXt``, which is then used to validate the
performance of the trained emulator.

.. GENERATED FROM PYTHON SOURCE LINES 52-54

To build a `PP GaSP` emulator for the above simulator,  first import class
:class:`.PPGaSP` by

.. GENERATED FROM PYTHON SOURCE LINES 55-58

.. code-block:: default


    from psimpy.emulator import PPGaSP








.. GENERATED FROM PYTHON SOURCE LINES 59-63

Then, create an instance of :class:`.PPGaSP`.  The parameter ``ndim``
(dimension of function input ``x``) must be specified. Optional parameters, such
as ``method``, ``kernel_type``, etc., can be set up if desired. Here, we leave
all the optional parameters to their default values.

.. GENERATED FROM PYTHON SOURCE LINES 64-68

.. code-block:: default


    emulator = PPGaSP(ndim=humanityX.shape[1])









.. GENERATED FROM PYTHON SOURCE LINES 69-71

Next, we train the `PP GaSP` emulator based on the training data using
:py:meth:`.PPGaSP.train` method.

.. GENERATED FROM PYTHON SOURCE LINES 72-75

.. code-block:: default


    emulator.train(design=humanityX, response=humanityY)





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

 .. code-block:: none

    The upper bounds of the range parameters are 299.331 299.54 297.1066 298.2833 298.9167 298.2838 299.6795 298.7972 299.2889 300.6043 300.4275 301.0958 301.0958 
    The initial values of range parameters are 5.98662 5.9908 5.942133 5.965665 5.978334 5.965676 5.99359 5.975944 5.985778 6.012085 6.008551 6.021916 6.021916 
    Start of the optimization  1  : 
    The number of iterations is  46 
     The value of the  marginal posterior  function is  -5279.534 
     Optimized range parameters are 25.46132 2.891428 5.30812 26.95519 291.1829 48.17274 84.56133 2.874385 39.12121 51.81304 0.5407651 1.728696 1.020605 
     Optimized nugget parameter is 0 
     Convergence:  TRUE 
    The initial values of range parameters are 12.37503 12.38367 12.28307 12.33172 12.3579 12.33174 12.38944 12.35296 12.37329 12.42767 12.42037 12.44799 12.44799 
    Start of the optimization  2  : 
    The number of iterations is  43 
     The value of the  marginal posterior  function is  -5279.534 
     Optimized range parameters are 25.46132 2.891428 5.308119 26.95519 291.1831 48.17271 84.56135 2.874385 39.1212 51.81301 0.5407651 1.728696 1.020605 
     Optimized nugget parameter is 0 
     Convergence:  TRUE 




.. GENERATED FROM PYTHON SOURCE LINES 76-79

With the trained emulator, we can make predictions for any arbitrary set of
input points using :py:meth:`PPGaSP.predict` method.
Here, we make predictions at testing input points ``humanityXt``.

.. GENERATED FROM PYTHON SOURCE LINES 80-83

.. code-block:: default


    predictions = emulator.predict(humanityXt)








.. GENERATED FROM PYTHON SOURCE LINES 84-86

We can validate the performance of the trained emulator based on the true outputs
``humanityYt``.

.. GENERATED FROM PYTHON SOURCE LINES 87-106

.. code-block:: default


    import matplotlib.pyplot as plt

    fig, ax = plt.subplots(5, 1, figsize=(6,15))

    for i in range(humanityY.shape[1]):

        ax[i].set_xlabel(f'Actual $y_{i+1}$')
        ax[i].set_ylabel(f'Emulator-predicted $y_{i+1}$')
        ax[i].set_xlim(np.min(humanityYt[:,i]),np.max(humanityYt[:,i]))
        ax[i].set_ylim(np.min(humanityYt[:,i]),np.max(humanityYt[:,i]))

        _ = ax[i].plot([np.min(humanityYt[:,i]),np.max(humanityYt[:,i])], [np.min(humanityYt[:,i]),np.max(humanityYt[:,i])])
        _ = ax[i].errorbar(humanityYt[:,i], predictions[:,i,0], predictions[:,i,3], fmt='.', linestyle='', label='prediction and std')
        _ = ax[i].legend()

    plt.tight_layout()





.. image-sg:: /auto_examples/emulator/images/sphx_glr_plot_ppgasp_001.png
   :alt: plot ppgasp
   :srcset: /auto_examples/emulator/images/sphx_glr_plot_ppgasp_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 107-109

We can also draw any number of samples at testing input ``humanityXt`` using 
:py:meth:`.PPGaSPsample()` method.

.. GENERATED FROM PYTHON SOURCE LINES 110-112

.. code-block:: default


    samples = emulator.sample(humanityXt, nsamples=10)
    print("Shape of samples: ", samples.shape)



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

 .. code-block:: none

    Shape of samples:  (120, 5, 10)





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

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


.. _sphx_glr_download_auto_examples_emulator_plot_ppgasp.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_ppgasp.py <plot_ppgasp.py>`

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

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


.. only:: html

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

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