qmat.lagrange
=============

.. py:module:: qmat.lagrange

.. autoapi-nested-parse::

   Base module for Barycentric Lagrange Approximation, based on `[Berrut & Trefethen, 2004] <https://doi.org/10.1137/S0036144502417715>`_.
   Allows to easily build integration / interpolation / derivative matrices, from any list of node points.

   .. rubric:: Examples

   >>> # Base usage to generate a quadrature matrix
   >>> from qmat.lagrange import LagrangeApproximation, np
   >>>
   >>> grid = np.linspace(0, 1, num=5)
   >>> approx = LagrangeApproximation(grid)
   >>> Q = approx.getIntegrationMatrix([(0, tau) for tau in grid])
   >>>
   >>> # Interpolation
   >>> fGrid = np.linspace(0, 1, num=200)
   >>> u = np.exp(grid)
   >>> P = approx.getInterpolationMatrix(fGrid)
   >>> uFine = P @ u
   >>>
   >>> # Alternative interpolation using the object as a function
   >>> uFine = approx(fGrid)
   >>>
   >>> # Derivative
   >>> D = approx.getDerivativeMatrix()
   >>> du = D @ u



Classes
-------

.. autoapisummary::

   qmat.lagrange.LagrangeApproximation


Functions
---------

.. autoapisummary::

   qmat.lagrange.computeFejerRule
   qmat.lagrange.getSparseInterpolationMatrix


Module Contents
---------------

.. py:function:: computeFejerRule(n)

   Compute a Fejer rule of the first kind, using DFT `[Waldvogel, 2006] <https://link.springer.com/article/10.1007/s10543-006-0045-4>`_.
   Inspired from quadpy (https://github.com/nschloe/quadpy @Nico_Schlömer)

   :Parameters: **n** (*int*) -- Number of points for the quadrature rule.

   :returns: * **nodes** (*np.1darray(n)*) -- The nodes of the quadrature rule
             * **weights** (*np.1darray(n)*) -- The weights of the quadrature rule.


.. py:class:: LagrangeApproximation(points, weightComputation='AUTO', scaleWeights=False, scaleRef='MAX', duplicates='USE_LEFT', fValues=None)



   Class approximating any function on a given set of points using barycentric
   Lagrange interpolation.

   Let note :math:`(t_j)_{0\leq j<n}` the set of points, then any scalar
   function :math:`f` can be approximated by the barycentric formula :

   .. math::
       p(x) =
       \frac{\displaystyle \sum_{j=0}^{n-1}\frac{w_j}{x-x_j}f_j}
       {\displaystyle \sum_{j=0}^{n-1}\frac{w_j}{x-x_j}},

   where :math:`f_j=f(t_j)` and

   .. math::
       w_j = \frac{1}{\prod_{k\neq j}(x_j-x_k)}

   are the barycentric weights.
   The theory and implementation is inspired from [1]_.

   :Parameters: * **points** (*list, tuple or np.1darray*) -- The given interpolation points, no specific scaling, but must be
                  ordered in increasing order.
                * **weightComputation** (*str, optional*) -- Algorithm used to compute the barycentric weights. Can be :

                  - 'FAST' : uses the analytic formula (unstable for large number of points)
                  - 'STABLE' : uses logarithmic difference and scaling of the weights
                  - 'CHEBFUN' : uses the same approach as in the chebfun package

                  The default is 'AUTO' : it tries the 'FAST' algorithm, and if an
                  overflow is detected, it switches to the 'STABLE' algorithm.
                * **scaleRef** (*str, optional*) -- Scaling used in the 'STABLE' algorithm for weight computation.
                  Can be :

                  - 'ZERO' : scaling based on the weight for the value closest to :math:`t=0`.
                  - 'MAX' : scaling based on the maximum weight value.

                  The default is 'MAX'.
                * **duplicates** (*str*) -- Which strategy to use in case of duplicated values within the interpolation
                  points. Can be :

                  - 'USE_LEFT' : uses the first value from the left in the values vector
                  - 'USE_RIGHT' : uses the first value from the right in the values vector

                  The default is 'USE_LEFT'.
                * **fValues** (*list, tuple or np.1darray*) -- Function values to be used when evaluating the LagrangeApproximation as a function

   :ivar points: The interpolating points
   :vartype points: np.1darray
   :ivar weights: The associated barycentric weights
   :vartype weights: np.1darray
   :ivar nPoints: The number of points, can also be retrieve with `n` (legacy alias)
   :vartype nPoints: int (property)
   :ivar uniquePoints: The unique interpolating points.
                       When there is no duplicates, points == uniquePoints.
   :vartype uniquePoints: np.1darray
   :ivar nUniquePoints: The number of unique points
   :vartype nUniquePoints: int (property)
   :ivar duplicates: The strategy used when there is duplicated interpolation points

   :vartype duplicates: str

   .. rubric:: References

   .. [1] Berrut, J. P., & Trefethen, L. N. (2004).
       "Barycentric Lagrange interpolation." SIAM review, 46(3), 501-517.
       URL: https://doi.org/10.1137/S0036144502417715


   .. py:attribute:: weights


   .. py:attribute:: weightComputation
      :value: 'AUTO'



   .. py:attribute:: fValues
      :value: None



   .. py:property:: nPoints
      :type: int


      The number of points


   .. py:attribute:: n

      Legacy alias for nPoints


   .. py:property:: nUniquePoints
      :type: int


      The number of unique points


   .. py:property:: hasDuplicates
      :type: bool


      Wether the points have duplicates or not


   .. py:method:: getInterpolationMatrix(times, duplicates=True)

      Compute the interpolation matrix for a given set of discrete "time"
      points.

      For instance, if we note :math:`\vec{f}` the vector containing the
      :math:`f_j=f(t_j)` values, and :math:`(\tau_m)_{0\leq m<M}`
      the "time" points where to interpolate.
      Then :math:`I[\vec{f}]`, the vector containing the interpolated
      :math:`f(\tau_m)` values, can be obtained using :

      .. math::
          I[\vec{f}] = P_{Inter} \vec{f},

      where :math:`P_{Inter}` is the interpolation matrix returned by this
      method.

      :Parameters: * **times** (*list-like or np.1darray*) -- The discrete "time" points where to interpolate the function.
                   * **duplicates** (*bool*) -- Wether or not take into account duplicates in the points.
                     This has no impact if all interpolating points are distinct.
                     Default is True.

      :returns: **P** -- The interpolation matrix, with :math:`M` rows (size of the **times**
                parameter) and :math:`n` columns.
      :rtype: np.2darray(M, n)



   .. py:method:: getIntegrationMatrix(intervals, numQuad='FEJER', duplicates=True)

      Compute the integration matrix for a given set of intervals.

      For instance, if we note :math:`\vec{f}` the vector containing the
      :math:`f_j=f(t_j)` values, and
      :math:`(\tau_{m,left}, \tau_{m,right})_{0\leq m<M}` the different
      interval where the function should be integrated using the barycentric
      interpolant polynomial.
      Then :math:`\Delta[\vec{f}]`, the vector containing the approximations
      of

      .. math::
          \int_{\tau_{m,left}}^{\tau_{m,right}} f(t)dt,

      can be obtained using :

      .. math::
          \Delta[\vec{f}] = P_{Integ} \vec{f},

      where :math:`P_{Integ}` is the interpolation matrix returned by this
      method.

      :Parameters: * **intervals** (*list of pairs*) -- A list of all integration intervals.
                   * **numQuad** (*str, optional*) -- Quadrature rule used to integrate the interpolant barycentric
                     polynomial. Can be :

                     - 'LEGENDRE_NUMPY' : Gauss-Legendre rule from Numpy
                     - 'LEGENDRE_SCIPY' : Gauss-Legendre rule from Scipy
                     - 'FEJER' : internally implemented Fejer-I rule

                     The default is 'FEJER'.
                   * **duplicates** (*bool*) -- Wether or not take into account duplicates in the points.
                     This has no impact if all interpolating points are distinct.
                     Default is True.

      :returns: **Q** -- The integration matrix, with :math:`M` rows (number of intervals)
                and :math:`n` columns.
      :rtype: np.2darray(M, n)



   .. py:method:: getDerivativeMatrix(order=1, duplicates=True)

      Generate derivative matrix of first or second order (or both) based on
      the Lagrange interpolant.
      The first order differentiation matrix :math:`D^{(1)}` approximates

      .. math::
          D^{(1)} u \simeq \frac{du}{dx}

      on the interpolation points. The formula is :

      .. math::
          D^{(1)}_{ij} = \frac{w_j/w_i}{x_i-x_j}

      for :math:`i \neq j` and

      .. math::
          D^{(1)}_{jj} = -\sum_{i \neq j} D^{(1)}_{ij}`

      The second order differentiation matrix :math:`D^{(2)}` approximates

      .. math::
          D^{(2)} u \simeq \frac{d^2u}{dx^2}

      on the interpolation points. The formula is :

      .. math::
          D^{(1)}_{ij} = -2\frac{w_j/w_i}{x_i-x_j}\left[
              \frac{1}{x_i-x_j} + \sum_{k \neq i}\frac{w_k/w_i}{x_i-x_k}
              \right]

      for :math:`i \neq j` and

      .. math::
          D^{(2)}_{jj} = -\sum_{i \neq j} D^{(2)}_{ij}

      ⚠️ If you want a derivative matrix with many points (~1000 or more),
      favor the use of `weightComputation="STABLE"` when initializing
      the `LagrangeApproximation` object. If not, some (very small) weights
      could be approximated by zeros, which would make the computation
      of the derivative matrices fail ...

      .. note::

         There is a typo in the formula for :math:`D^{(2)}` given in the paper
         of Berrut and Trefethen. The formula above is the correct one.

      :Parameters: * **order** (*int or str, optional*) -- The order of the derivative matrix, use "ALL" to retrieve both.
                     The default is 1.
                   * **duplicates** (*bool*) -- Wether or not take into account duplicates in the points.
                     This has no impact if all interpolating points are distinct.
                     Default is True.

      :returns: **D** -- Derivative matrix. If order="ALL", return a tuple containing all
                derivative matrices in increasing derivative order.
      :rtype: np.2darray or tuple of np.2darray



   .. py:method:: getDerivationMatrix(*args, **kwargs)


.. py:function:: getSparseInterpolationMatrix(inPoints, outPoints, order)

   Get a sparse interpolation matrix from `inPoints` to `outPoints` of order
   `order` using barycentric Lagrange interpolation.

   The matrix will have `order` entries per line, and tends to be banded when
   both `inPoints` and `outPoints` are equispaced and cover the same interval.

   :Parameters: * **inPoints** (*np.1darray*) -- The points you want to interpolate from
                * **outPoints** (*np.1darray*) -- The points you want to interpolate to
                * **order** (*int*) -- Order of the interpolation

   :returns: **A** -- Sparse interpolation matrix
   :rtype: scipy.sparse.csc_matrix(len(outPoints), len(inPoints))