.. _sphx_glr_examples_plot_ellipsoids.py:


===================
Bounding Ellipsoids
===================

Illustration of determining bounding ellipsoids in 3 dimensions.

For complex iso-likelihood surfaces (ones not well-approximated by a
single ellipsoid), Nestle will bound the points in multiple ellipsoids
when using the ``'multi'`` method. Here, we draw points uniformly
from a torus and use Nestle to determine an optimal set of ellipses
bounding the points.



.. code-block:: python


    import math

    import numpy as np
    from matplotlib import pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D

    import nestle

    # seed RNG for reproducible results
    np.random.seed(0)


    def rand_torus(ro, ri, npoints):
        """Generate points within a torus wth major radius `ro` and minor radius
        `ri` via rejection sampling"""
        out = np.empty((npoints, 3), dtype=np.float64)
        i = 0
        while i < npoints:
            # generate point within box
            x = np.random.uniform(-1., 1., size=3)
            x[0:2] *= ro + ri
            x[2] *= ri
        
            r = math.sqrt(x[0]**2 + x[1]**2) - ro
            if (r**2 + x[2]**2 < ri**2):
                out[i, :] = x
                i += 1
            
        return out


    def plot_ellipsoid_3d(ell, ax):
        """Plot the 3-d Ellipsoid ell on the Axes3D ax."""
    
        # points on unit sphere
        u = np.linspace(0.0, 2.0 * np.pi, 100)
        v = np.linspace(0.0, np.pi, 100)
        z = np.outer(np.cos(u), np.sin(v))
        y = np.outer(np.sin(u), np.sin(v))
        x = np.outer(np.ones_like(u), np.cos(v))
    
        # transform points to ellipsoid
        for i in range(len(x)):
            for j in range(len(x)):
                x[i,j], y[i,j], z[i,j] = ell.ctr + np.dot(ell.axes,
                                                          [x[i,j],y[i,j],z[i,j]])

        ax.plot_wireframe(x, y, z,  rstride=4, cstride=4, color='#2980b9', alpha=0.2)


    # Generate points within a torus
    npoints = 1000
    R = 1.0
    r = 0.1
    points = rand_torus(R, r, npoints)
    torus_vol = 2. * math.pi**2 * R * r**2

    # Determine bounding ellipsoids
    pointvol = torus_vol / npoints
    ells = nestle.bounding_ellipsoids(points, pointvol=pointvol)

    # plot
    fig = plt.figure(figsize=(10., 10.))
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(points[:, 0], points[:, 1], points[:, 2], c='k', marker='.')
    for ell in ells:
        plot_ellipsoid_3d(ell, ax)

    ax.set_xlim(-1., 1.)
    ax.set_ylim(-1., 1.)
    ax.set_zlim(-1., 1.)
    ax.set_aspect('equal')
    fig.tight_layout()
    plt.show()





.. image:: /examples/images/sphx_glr_plot_ellipsoids_001.png
    :align: center




Single ellipsoid
----------------

Nestle correctly identifies a cloud of points as belonging to a single
ellipsoid.



.. code-block:: python


    # Generate samples drawn from a single ellipsoid
    npoints = 100
    A = np.array([[0.25, 1., 0.5],
                  [1., 0.25, 0.5],
                  [0.5, 1., 0.25]])
    ell_gen = nestle.Ellipsoid([0., 0., 0.], np.dot(A.T, A))
    points = ell_gen.samples(npoints)
    pointvol = ell_gen.vol / npoints

    # Find bounding ellipsoid(s)
    ells = nestle.bounding_ellipsoids(points, pointvol)

    # plot
    fig = plt.figure(figsize=(10., 10.))
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(points[:, 0], points[:, 1], points[:, 2], c='k', marker='.')
    for ell in ells:
        plot_ellipsoid_3d(ell, ax)

    fig.tight_layout()
    plt.show()



.. image:: /examples/images/sphx_glr_plot_ellipsoids_002.png
    :align: center




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



.. container:: sphx-glr-footer


  .. container:: sphx-glr-download

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



  .. container:: sphx-glr-download

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

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

    `Generated by Sphinx-Gallery <http://sphinx-gallery.readthedocs.io>`_