.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/stochastic_processes/bispectral/bispectral_nd.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_stochastic_processes_bispectral_bispectral_nd.py>`
        to download the full example code or to run this example in your browser via Binder

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

.. _sphx_glr_auto_examples_stochastic_processes_bispectral_bispectral_nd.py:


Bispectral Representation Method - ND
=================================================================

In this example, the BiSpectral Representation Method is used to generate stochastic processes from a prescribed Power
Spectrum and associated Bispectrum. This example illustrates how to use the BSRM class for 'n' dimensional case and
compare the statistics of the generated stochastic processes with the expected values.

.. GENERATED FROM PYTHON SOURCE LINES 13-15

Import the necessary libraries. Here we import standard libraries such as numpy and matplotlib, but also need to
import the :class:`.BispectralRepresentation` class from the :py:mod:`stochastic_processes` module of UQpy.

.. GENERATED FROM PYTHON SOURCE LINES 18-26

.. code-block:: default


    from UQpy.stochastic_process import BispectralRepresentation
    import numpy as np
    from scipy.stats import skew
    import matplotlib.pyplot as plt
    from copy import deepcopy
    plt.style.use('seaborn')


.. GENERATED FROM PYTHON SOURCE LINES 27-28

The input parameters necessary for the generation of the stochastic processes are given below:

.. GENERATED FROM PYTHON SOURCE LINES 31-68

.. code-block:: default


    n_sim = 1000  # Num of samples
    dim = 2
    n = 1  # Num of dimensions

    T = 20  # Time(1 / T = dw)
    nt = 128  # Num.of Discretized Time
    F = 1 / T * nt / 2  # Frequency.(Hz)
    nf = 64  # Num of Discretized Freq.

    # # Generation of Input Data(Stationary)
    dt = T / nt
    t = np.linspace(0, T - dt, nt)
    df = F / nf
    f = np.linspace(0, F - df, nf)

    f_list = [f for _ in range(dim)]
    F_P = np.array(np.meshgrid(*f_list, indexing='ij'))
    P = 20 / np.sqrt(2 * np.pi) * np.exp(-1/2 * np.linalg.norm(F_P, axis=0) ** 2)

    F_B = np.meshgrid(*[*f_list, *f_list])
    b = 40 / (2 * np.pi) * np.exp(2 * (-1/2 * np.linalg.norm(F_B, axis=0) ** 2))
    B_Real = deepcopy(b)
    B_Imag = deepcopy(b)

    B_Real[0, :, :, :] = 0
    B_Real[:, 0, :, :] = 0
    B_Real[:, :, 0, :] = 0
    B_Real[:, :, :, 0] = 0
    B_Imag[0, :, :, :] = 0
    B_Imag[:, 0, :, :] = 0
    B_Imag[:, :, 0, :] = 0
    B_Imag[:, :, :, 0] = 0

    B_Complex = B_Real + 1j * B_Imag



.. GENERATED FROM PYTHON SOURCE LINES 69-70

Make sure that the input parameters are in order to prevent aliasing.

.. GENERATED FROM PYTHON SOURCE LINES 73-79

.. code-block:: default


    t_u = 2*np.pi/2/F

    if dt>t_u:
        print('Error')


.. GENERATED FROM PYTHON SOURCE LINES 80-81

Defining the Power Spectral Density Function.

.. GENERATED FROM PYTHON SOURCE LINES 84-103

.. code-block:: default


    BSRM_object = BispectralRepresentation(n_sim, P, B_Complex, [dt, dt], [df, df], [nt, nt], [nf, nf])
    samples = BSRM_object.samples

    from pylab import *
    t_list = [t for _ in range(dim)]
    T = np.array(np.meshgrid(*t_list, indexing='ij'))

    fig1 = plt.figure()
    plt.title('2d random field with a prescribed Power Spectrum and Bispectrum')
    pcm = pcolor(T[0], T[1], samples[0, 0], cmap='RdBu_r', vmin=-30, vmax=30)
    plt.colorbar(pcm, extend='both', orientation='vertical')
    plt.xlabel('$X_{1}$')
    plt.ylabel('$X_{2}$')
    plt.show()

    print('The mean of the samples is ', np.mean(samples), 'whereas the expected mean is 0.000')
    print('The variance of the samples is ', np.var(samples), 'whereas the expected variance is ', np.sum(P)*df*2)
    print('The skewness of the samples is ', np.mean(skew(samples, axis=0)), 'whereas the expected skewness is ',
          np.sum(B_Real)*df**2*6/(np.sum(P)*df*2)**(3/2))

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

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


.. _sphx_glr_download_auto_examples_stochastic_processes_bispectral_bispectral_nd.py:

.. only:: html

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


    .. container:: binder-badge

      .. image:: images/binder_badge_logo.svg
        :target: https://mybinder.org/v2/gh/SURGroup/UQpy/master?urlpath=lab/tree/notebooks/auto_examples/stochastic_processes/bispectral/bispectral_nd.ipynb
        :alt: Launch binder
        :width: 150 px



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

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

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

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


.. only:: html

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

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