.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/algorithms/nd_interpolation.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_algorithms_nd_interpolation.py>`
        to download the full example code.

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

.. _sphx_glr_auto_examples_algorithms_nd_interpolation.py:


Use N-dimensional Histogram functionality and Interpolation
===========================================================

-  could be used for example to read and interpolate an lookup table or
   IRF.
-  In this example, we load a sample energy reconstruction lookup-table
   from a FITS file
-  In this case it is only in 2D cube (to keep the file size small):
   ``SIZE`` vs ``IMPACT-DISTANCE``, however the same method will work
   for any dimensionality

.. GENERATED FROM PYTHON SOURCE LINES 14-26

.. code-block:: Python
   :lineno-start: 15


    import matplotlib.pylab as plt
    import numpy as np
    from astropy.io import fits
    from scipy.interpolate import RegularGridInterpolator

    from ctapipe.utils import Histogram
    from ctapipe.utils.datasets import get_dataset_path

    # %matplotlib inline









.. GENERATED FROM PYTHON SOURCE LINES 27-35

load an example datacube
~~~~~~~~~~~~~~~~~~~~~~~~

(an energy table generated for a small subset of HESS simulations) to
use as a lookup table. Here we will use the ``Histogram`` class, which
automatically loads both the data cube and creates arrays for the
coordinates of each axis.


.. GENERATED FROM PYTHON SOURCE LINES 35-42

.. code-block:: Python
   :lineno-start: 36


    testfile = get_dataset_path("hess_ct001_energylut.fits.gz")
    energy_hdu = fits.open(testfile)["MEAN"]
    energy_table = Histogram.from_fits(energy_hdu)
    print(energy_table)






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

 .. code-block:: none


    Downloading hess_ct001_energylut.fits.gz:   0%|          | 0.00/29.9k [00:00<?, ?B/s]
    Downloading hess_ct001_energylut.fits.gz:  34%|███▍      | 10.2k/29.9k [00:00<00:00, 84.6kB/s]
    Downloading hess_ct001_energylut.fits.gz: 100%|██████████| 29.9k/29.9k [00:00<00:00, 241kB/s] 
    Histogram(name='Histogram', axes=['LSIZ' 'DIST'], nbins=[100 100], ranges=[[5.e-01 6.e+00]
     [0.e+00 2.e+03]])




.. GENERATED FROM PYTHON SOURCE LINES 43-50

construct an interpolator that we can use to get values at any point:
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Here we will use a ``RegularGridInterpolator``, since it is the most
appropriate for this type of data, but others are available (see the
SciPy documentation)


.. GENERATED FROM PYTHON SOURCE LINES 50-57

.. code-block:: Python
   :lineno-start: 51


    centers = [energy_table.bin_centers(ii) for ii in range(energy_table.ndims)]
    energy_lookup = RegularGridInterpolator(
        centers, energy_table.hist, bounds_error=False, fill_value=-100
    )









.. GENERATED FROM PYTHON SOURCE LINES 58-69

``energy_lookup`` is now just a continuous function of ``log(SIZE)``,
``DISTANCE`` in m.

Now plot some curves from the interpolator.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Note that the LUT we used is does not have very high statistics, so the
interpolation starts to be affected by noise at the high end. In a real
case, we would want to use a table that has been sanitized (smoothed and
extrapolated)


.. GENERATED FROM PYTHON SOURCE LINES 69-89

.. code-block:: Python
   :lineno-start: 70


    lsize = np.linspace(1.5, 5.0, 100)
    dists = np.linspace(50, 100, 5)

    plt.figure()
    plt.title("Variation of energy with size and impact distance")
    plt.xlabel("SIZE (P.E.)")
    plt.ylabel("ENERGY (TeV)")

    for dist in dists:
        plt.plot(
            10**lsize,
            10 ** energy_lookup((lsize, dist)),
            "+-",
            label="DIST={:.1f} m".format(dist),
        )

    plt.legend(loc="best")





.. image-sg:: /auto_examples/algorithms/images/sphx_glr_nd_interpolation_001.png
   :alt: Variation of energy with size and impact distance
   :srcset: /auto_examples/algorithms/images/sphx_glr_nd_interpolation_001.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


    <matplotlib.legend.Legend object at 0x79f08c49cee0>



.. GENERATED FROM PYTHON SOURCE LINES 90-93

Using the interpolator, reinterpolate the lookup table onto an
:math:`N \times N` grid (regardless of its original dimensions):


.. GENERATED FROM PYTHON SOURCE LINES 93-103

.. code-block:: Python
   :lineno-start: 94


    N = 300
    xmin, xmax = energy_table.bin_centers(0)[0], energy_table.bin_centers(0)[-1]
    ymin, ymax = energy_table.bin_centers(1)[0], energy_table.bin_centers(1)[-1]
    xx, yy = np.linspace(xmin, xmax, N), np.linspace(ymin, ymax, N)
    X, Y = np.meshgrid(xx, yy)
    points = list(zip(X.ravel(), Y.ravel()))
    E = energy_lookup(points).reshape((N, N))









.. GENERATED FROM PYTHON SOURCE LINES 104-107

Now, let’s plot the original table and the new one (E). The color bar
shows :math:`\log_{10}(E)` in TeV


.. GENERATED FROM PYTHON SOURCE LINES 107-139

.. code-block:: Python
   :lineno-start: 108


    plt.figure(figsize=(12, 5))
    plt.nipy_spectral()

    # the uninterpolated table
    plt.subplot(1, 2, 1)
    plt.xlim(1.5, 5)
    plt.ylim(0, 500)
    plt.xlabel("log10(SIZE)")
    plt.ylabel("Impact Dist (m)")
    plt.pcolormesh(
        energy_table.bin_centers(0), energy_table.bin_centers(1), energy_table.hist.T
    )
    plt.title(f"Raw table, uninterpolated {energy_table.hist.T.shape}")
    cb = plt.colorbar()
    cb.set_label(r"$\log_{10}(E/\mathrm{TeV})$")

    # the interpolated table
    plt.subplot(1, 2, 2)
    plt.pcolormesh(np.linspace(xmin, xmax, N), np.linspace(ymin, ymax, N), E)
    plt.xlim(1.5, 5)
    plt.ylim(0, 500)
    plt.xlabel("log10(SIZE)")
    plt.ylabel("Impact Dist(m)")
    plt.title("Interpolated to a ({0}, {0}) grid".format(N))
    cb = plt.colorbar()
    cb.set_label(r"$\log_{10}(E/\mathrm{TeV})$")

    plt.tight_layout()
    plt.show()





.. image-sg:: /auto_examples/algorithms/images/sphx_glr_nd_interpolation_002.png
   :alt: Raw table, uninterpolated (100, 100), Interpolated to a (300, 300) grid
   :srcset: /auto_examples/algorithms/images/sphx_glr_nd_interpolation_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 140-145

In the high-stats central region, we get a nice smooth interpolation
function. Of course we can see that there are a few more steps to take
before using this table: \*
- need to deal with cases where the table had low stats near the edges (smooth or extrapolate, or set bounds)
- may need to smooth the table even where there are sufficient stats, to avoid wiggles


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

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


.. _sphx_glr_download_auto_examples_algorithms_nd_interpolation.py:

.. only:: html

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

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

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

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

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

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: nd_interpolation.zip <nd_interpolation.zip>`


.. only:: html

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

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