pywafo/wafo/padua.py

'''

All the software  contained in this library  is protected by copyright.
Permission  to use, copy, modify, and  distribute this software for any
purpose without fee is hereby granted, provided that this entire notice
is included  in all copies  of any software which is or includes a copy
or modification  of this software  and in all copies  of the supporting
documentation for such software.

THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
WARRANTY. IN NO EVENT, NEITHER  THE AUTHORS, NOR THE PUBLISHER, NOR ANY
MEMBER  OF THE EDITORIAL BOARD OF  THE JOURNAL  "NUMERICAL ALGORITHMS",
NOR ITS EDITOR-IN-CHIEF, BE  LIABLE FOR ANY ERROR  IN THE SOFTWARE, ANY
MISUSE  OF IT  OR ANY DAMAGE ARISING OUT OF ITS USE. THE ENTIRE RISK OF
USING THE SOFTWARE LIES WITH THE PARTY DOING SO.

ANY USE OF THE SOFTWARE  CONSTITUTES  ACCEPTANCE  OF THE TERMS  OF THE
ABOVE STATEMENT.


AUTHORS:
Per A Brodtkorb
Python code Based on matlab code written by:

Marco Caliari
University of Verona, Italy
E-mail: marco.caliari@univr.it

Stefano de Marchi, Alvise Sommariva, Marco Vianello
University of Padua, Italy
E-mail: demarchi@math.unipd.it, alvise@math.unipd.it,
        marcov@math.unipd.it

Reference
---------
Padua2DM: fast interpolation and cubature at the Padua points in Matlab/Octave
NUMERICAL ALGORITHMS, 56 (2011), PP. 45-60


Padua module
------------
In polynomial interpolation of two variables, the Padua points are the first
known example (and up to now the only one) of a unisolvent point set
(that is, the interpolating polynomial is unique) with minimal growth of their
Lebesgue constant, proven to be O(log2 n).
This module provides all the functions needed to perform interpolation and
cubature at the Padua points, together with the functions and the demos used
in the paper.

pdint.m                 : main function for interpolation at the Padua points
pdcub.m                 : main function for cubature at the Padua points
pdpts.m                 : function for the computation of the Padua points
padua_fit.m              : function for the computation of the interpolation
                          coefficients by FFT (recommended)
pdcfsMM.m               : function for the computation of the interpolation
                          coefficients by matrix multiplications
pdval.m                 : function for the evaluation of the interpolation
                          polynomial
pdwtsFFT.m              : function for the computation of the cubature
                          weights by FFT
pdwtsMM.m               : function for the computation of the cubature
                          weights by matrix multiplications (recommended)
funct.m                 : function containing some test functions
demo_pdint.m            : demo script for pdint
demo_cputime_pdint.m    : demo script for the computation of CPU time for
                          interpolation
demo_errors_pdint.m     : demo script for the comparison of interpolation with
                          coefficients computed by FFT or by matrix
                          multiplications
demo_pdcub              : demo script for pdcub
demo_cputime_pdcub.m    : demo script for the computation of CPU time for
                          cubature
demo_errors_pdcub.m     : demo script for the comparison of cubature with
                          weights computed by FFT or by matrix multiplications
demo_errors_pdcub_gl.m  : demo script for the comparison of different cubature
                          formulas
cubature_square.m       : function for the computation of some cubature
                          formulas for the square
omelyan_solovyan_rule.m : function for the computation of Omelyan-Solovyan
                          cubature points and weights
Contents.m              : Contents file for Matlab


'''
from __future__ import absolute_import, division
import numpy as np
from numpy.fft import fft
from wafo.dctpack import dct
from wafo.polynomial import map_from_interval, map_to_interval
# from scipy.fftpack.realtransforms import dct


class _ExampleFunctions(object):
    '''
    Computes test function in the points (x, y)

    Parameters
    ----------
    x,y : array-like
        evaluate the function in the points (x,y)
    i : scalar int (default 0)
        defining which test function to use. Options are
        0: franke
        1: half_sphere
        2: poly_degree20
        3: exp_fun1
        4: exp_fun100
        5: cos30
        6: constant
        7: exp_xy
        8: runge
        9: abs_cubed
        10: gauss
        11: exp_inv

    Returns
    -------
    z : array-like
        value of the function in the points (x,y)
    '''
    @staticmethod
    def franke(x, y):
        '''Franke function.

        The value of the definite integral on the square [-1,1] x [-1,1],
        obtained using a Padua Points cubature formula of degree 500, is
        2.1547794245591083e+000 with an estimated absolute error of 8.88e-016.

        The value of the definite integral on the square [0,1] x [0,1],
        obtained using a Padua Points cubature formula of degree 500, is
        4.06969589491556e-01 with an estimated absolute error of 8.88e-016.

        Maple: 0.40696958949155611906
        '''
        def _exp(x, y, loc, scale, p2=2):
            return np.exp(- (x-loc[0])**2/scale[0] - (y-loc[1])**p2/scale[1])
        # exp = np.exp
        x9, y9 = 9. * x, 9. * y
        return (3. / 4 * _exp(x9, y9, [2, 2], [4, 4]) +
                3. / 4 * _exp(x9, y9, [-1, -1], [49, 10], p2=1) +
                1. / 2 * _exp(x9, y9, [7, 3], [4, 4]) -
                1. / 5 * _exp(x9, y9, [4, 7], [1, 1]))

    @staticmethod
    def half_sphere(x, y):
        '''The value of the definite integral on the square [-1,1] x [-1,1],
        obtained using a Padua Points cubature formula of degree 2000, is
        3.9129044444568244e+000 with an estimated absolute error of 3.22e-010.
        '''
        return ((x - 0.5)**2 + (y - 0.5)**2)**(1. / 2)

    @staticmethod
    def poly_degree20(x, y):
        ''''Bivariate polynomial having moderate degree.
        The value of the definite integral on the square [-1,1] x
        [-1,1], up to machine precision, is 18157.16017316017 (see ref. 6).
        The value of the definite integral on the square [-1,1] x [-1,1],
        obtained using a Padua Points cubature formula of degree 500,
        is 1.8157160173160162e+004.

        2D modification of an example by L.N.Trefethen (see ref. 7), f(x)=x^20.
        '''
        return (x + y)**20

    @staticmethod
    def exp_fun1(x, y):
        ''' The value of the definite integral on the square [-1,1] x [-1,1],
        obtained using a Padua Points cubature formula of degree 2000,
        is 2.1234596326670683e+001 with an estimated absolute error of
        7.11e-015.
        '''
        return np.exp((x - 0.5)**2 + (y - 0.5)**2)

    @staticmethod
    def exp_fun100(x, y):
        '''The value of the definite integral on the square [-1,1] x [-1,1],
        obtained using a Padua Points cubature formula of degree 2000,
        is 3.1415926535849605e-002 with an estimated absolute error of
        3.47e-017.
        '''
        return np.exp(-100 * ((x - 0.5)**2 + (y - 0.5)**2))

    @staticmethod
    def cos30(x, y):
        ''' The value of the definite integral on the square [-1,1] x [-1,1],
        obtained using a Padua Points cubature formula of degree 500,
        is 4.3386955120336568e-003 with an estimated absolute error of
        2.95e-017.
        '''
        return np.cos(30 * (x + y))

    @staticmethod
    def constant(x, y):
        '''Constant.
        To test interpolation and cubature at degree 0.
        The value of the definite integral on the square [-1,1] x [-1,1]
        is 4.
        '''
        return np.ones(np.shape(x+y))

    @staticmethod
    def exp_xy(x, y):
        '''The value of the definite integral on the square [-1,1] x [-1,1]
        is up to machine precision is 5.524391382167263 (see ref. 6).
        2. The value of the definite integral on the square [-1,1] x [-1,1],
        obtained using a Padua Points cubature formula of degree 500,
        is 5.5243913821672628e+000 with an estimated absolute error of
        0.00e+000.
        2D modification of an example by L.N.Trefethen (see ref. 7),
        f(x)=exp(x).
        '''
        return np.exp(x + y)

    @staticmethod
    def runge(x, y):
        ''' Bivariate Runge function: as 1D complex function is analytic
        in a neighborhood of [-1; 1] but not throughout the complex plane.

        The value of the definite integral on the square [-1,1] x [-1,1],
        up to machine precision, is 0.597388947274307 (see ref. 6).
        The value of the definite integral on the square [-1,1] x [-1,1],
        obtained using a Padua Points cubature formula of degree 500,
        is 5.9738894727430725e-001 with an estimated absolute error of
        0.00e+000.

        2D modification of an example by L.N.Trefethen (see ref. 7),
        f(x)=1/(1+16*x^2).
        '''
        return 1. / (1 + 16 * (x**2 + y**2))

    @staticmethod
    def abs_cubed(x, y):
        '''Low regular function.
        The value of the definite integral on the square [-1,1] x [-1,1],
        up to machine precision, is 2.508723139534059 (see ref. 6).
        The value of the definite integral on the square [-1,1] x [-1,1],
        obtained using a Padua Points cubature formula of degree 500,
        is 2.5087231395340579e+000 with an estimated absolute error of
        0.00e+000.

        2D modification of an example by L.N.Trefethen (see ref. 7),
        f(x)=abs(x)^3.
        '''
        return (x**2 + y**2)**(3 / 2)

    @staticmethod
    def gauss(x, y):
        '''Bivariate gaussian: smooth function.
        The value of the definite integral on the square [-1,1] x [-1,1],
        up to machine precision, is 2.230985141404135 (see ref. 6).
        The value of the definite integral on the square [-1,1] x [-1,1],
        obtained using a Padua Points cubature formula of degree 500,
        is 2.2309851414041333e+000 with an estimated absolute error of
        2.66e-015.

        2D modification of an example by L.N.Trefethen (see ref. 7),
        f(x)=exp(-x^2).
        '''
        return np.exp(-x**2 - y**2)

    @staticmethod
    def exp_inv(x, y):
        '''Bivariate example stemming from a 1D C-infinity function.
        The value of the definite integral on the square [-1,1] x [-1,1],
        up to machine precision, is 0.853358758654305 (see ref. 6).
        The value of the definite integral on the square [-1,1] x [-1,1],
        obtained using a Padua Points cubature formula of degree 2000,
        is 8.5335875865430544e-001 with an estimated absolute error of
        3.11e-015.

        2D modification of an example by L.N.Trefethen (see ref. 7),
        f(x)=exp(-1/x^2).
        '''
        arg_z = (x**2 + y**2)
        # Avoid cases in which "arg_z=0", setting only in those instances
        # "arg_z=eps".
        arg_z = arg_z + (1 - np.abs(np.sign(arg_z))) * 1.e-100
        arg_z = 1. / arg_z
        return np.exp(-arg_z)

    def __call__(self, x, y, i=0):
        s = self
        test_function = [s.franke, s.half_sphere, s.poly_degree20, s.exp_fun1,
                         s.exp_fun100, s.cos30, s.constant, s.exp_xy, s.runge,
                         s.abs_cubed, s.gauss, s.exp_inv]
        return test_function[i](x, y)
example_functions = _ExampleFunctions()


def _find_m(n):
    ix = np.r_[1:(n + 1) * (n + 2):2]
    if np.mod(n, 2) == 0:
        n2 = n // 2
        offset = np.array([[0, 1] * n2 + [0, ]] * (n2 + 1))
        ix = ix - offset.ravel(order='F')
    return ix


def padua_points(n, domain=(-1, 1, -1, 1)):
    ''' Return Padua points

    Parameters
    ----------
    n : scalar integer
         interpolation degree
    domain : vector [a,b,c,d]
        defining the rectangle [a,b] x [c,d]. (default domain = (-1,1,-1,1))

    Returns
    -------
    pad : array of shape (2 x (n+1)*(n+2)/2) such that
        (pad[0,:], pad[1,: ]) defines the Padua points in the domain
        rectangle [a,b] x [c,d].
    or
     X1,Y1,X2,Y2 : arrays
         Two subgrids X1,Y1 and X2,Y2 defining the Padua points
    -------------------------------------------------------------------------------
    '''
    a, b, c, d = domain
    t0 = [np.pi] if n == 0 else np.linspace(0, np.pi, n + 1)
    t1 = np.linspace(0, np.pi, n + 2)
    zn = map_to_interval(np.cos(t0), a, b)
    zn1 = map_to_interval(np.cos(t1), c, d)

    Pad1, Pad2 = np.meshgrid(zn, zn1)
    ix = _find_m(n)
    return np.vstack((Pad1.ravel(order='F')[ix],
                      Pad2.ravel(order='F')[ix]))


def error_estimate(C0f):
    ''' Return interpolation error estimate from Padua coefficients
    '''
    n = C0f.shape[1]
    C0f2 = np.fliplr(C0f)
    errest = sum(np.abs(np.diag(C0f2)))
    if (n >= 1):
        errest = errest + sum(np.abs(np.diag(C0f2, -1)))
        if (n >= 2):
            errest = errest + sum(np.abs(np.diag(C0f2, -2)))
    return 2 * errest


def padua_fit(Pad, fun, *args):
    '''
    Computes the Chebyshevs coefficients

    so that f(x, y) can be approximated by:

           f(x, y) = sum cjk*Tjk(x, y)

    Parameters
    ----------
    Pad : array-like
        Padua points, as computed  with padua_points function.
    fun : function to be interpolated in the form
        fun(x, y, *args), where *args are optional arguments for fun.

    Returns
    -------
    coefficents: coefficient matrix
    abs_err   : interpolation error estimate

    Example
    ------
    >>> import numpy as np
    >>> import wafo.padua as wp
    >>> domain = [0, 1, 0, 1]
    >>> a, b, c, d = domain
    >>> points = wp.padua_points(21, domain)
    >>> C0f, abs_error = wp.padua_fit(points, wp.example_functions.franke)
    >>> x1 = np.linspace(a, b, 100)
    >>> x2 = np.linspace(c, d, 101)
    >>> val = wp.padua_val(x1, x2, C0f, domain, use_meshgrid=True)
    >>> X1, X2 = np.meshgrid(x1, x2)
    >>> true_val = wp.example_functions.franke(X1, X2)

    >>> np.allclose(val, true_val, atol=10*abs_error)
    True
    >>> np.allclose(np.abs(val-true_val).max(), 0.0073174614275738296)
    True
    >>> np.allclose(abs_error, 0.0022486904061664046)
    True

    import matplotlib.pyplot as  plt
    plt.contour(x1, x2, val)

    '''

    N = np.shape(Pad)[1]
    # recover the degree n from N = (n+1)(n+2)/2
    n = int(round(-3 + np.sqrt(1 + 8 * N)) / 2)
    C0f = fun(Pad[0], Pad[1], *args)
    if (n > 0):
        ix = _find_m(n)
        GfT = np.zeros((n + 2) * (n + 1))
        GfT[ix] = C0f * 2 / (n * (n + 1))
        GfT = GfT.reshape(n + 1, n + 2)
        GfT = GfT.T
        GfT[0] = GfT[0] / 2
        GfT[n + 1] = GfT[n + 1] / 2
        GfT[:, 0] = GfT[:, 0] / 2
        GfT[:, n] = GfT[:, n] / 2
        Gf = GfT.T
        # compute the interpolation coefficient matrix C0f by FFT
        Gfhat = np.real(fft(Gf, 2 * n, axis=0))
        Gfhathat = np.real(fft(Gfhat[:n + 1, :], 2 * (n + 1), axis=1))
        C0f = 2 * Gfhathat[:, 0:n + 1]
        C0f[0] = C0f[0, :] / np.sqrt(2)
        C0f[:, 0] = C0f[:, 0] / np.sqrt(2)
        C0f = np.fliplr(np.triu(np.fliplr(C0f)))
        C0f[n] = C0f[n] / 2

    return C0f, error_estimate(C0f)


def paduavals2coefs(f):
    m = len(f)
    n = int(round(-1.5 + np.sqrt(.25 + 2 * m)))
    x = padua_points(n)
    idx = _find_m(n)
    w = 0 * x[0] + 1. / (n * (n + 1))
    idx1 = np.all(np.abs(x) == 1, axis=0)
    w[idx1] = .5 * w[idx1]
    idx2 = np.all(np.abs(x) != 1, axis=0)
    w[idx2] = 2 * w[idx2]

    G = np.zeros(idx.max() + 1)
    G[idx] = 4 * w * f

    use_dct = 100 < n
    if use_dct:
        C = np.rot90(dct(dct(G.T).T))  # , axis=1)
    else:
        t1 = np.r_[0:n + 1].reshape(-1, 1)
        Tn1 = np.cos(t1 * t1.T * np.pi / n)
        t2 = np.r_[0:n + 2].reshape(-1, 1)
        Tn2 = np.cos(t2 * t2.T * np.pi / (n + 1))
        C = np.dot(Tn2, np.dot(G, Tn1))

    C[0] = .5 * C[0]
    C[:, 1] = .5 * C[:, 1]
    C[0, -1] = .5 * C[0, -1]
    del C[-1]

    # Take upper-left triangular part:
    return np.fliplr(np.triu(np.fliplr(C)))
    # C = triu(C(:,end:-1:1));
    # C = C(:,end:-1:1);


# TODO: padua_fit2 does not work correctly yet.
def padua_fit2(Pad, fun, *args):
    # N = np.shape(Pad)[1]
    # recover the degree n from N = (n+1)(n+2)/2
    # _n = int(round(-3 + np.sqrt(1 + 8 * N)) / 2)
    C0f = fun(Pad[0], Pad[1], *args)
    return paduavals2coefs(C0f)


def _compute_moments(n):
    k = np.r_[0:n:2]
    mom = 2 * np.sqrt(2) / (1 - k ** 2)
    mom[0] = 2
    return mom


def padua_cubature(coefficients, domain=(-1, 1, -1, 1)):
    '''
    Compute the integral through the coefficient matrix.
    '''
    n = coefficients.shape[1]
    mom = _compute_moments(n)
    M1, M2 = np.meshgrid(mom, mom)
    M = M1 * M2
    C0fM = coefficients[0:n:2, 0:n:2] * M
    a, b, c, d = domain
    integral = (b - a) * (d - c) * C0fM.sum() / 4
    return integral


def padua_val(X, Y, coefficients, domain=(-1, 1, -1, 1), use_meshgrid=False):
    '''
    Evaluate polynomial in padua form at X, Y.

     Evaluate the interpolation polynomial defined through its coefficient
     matrix coefficients at the target points X(:,1),X(:,2) or at the
     meshgrid(X1,X2)

    Parameters
    ----------
    X, Y: array-like
        evaluation points.
    coefficients : array-like
         coefficient matrix
    domain : a vector [a,b,c,d]
         defining the rectangle [a,b] x [c,d]
    use_meshgrid: bool
        If True interpolate at the points meshgrid(X, Y)

    Returns
    -------
     fxy : array-like
         evaluation of the interpolation polynomial at the target points
    '''
    def transform(tn, x, a, b):
        xn = map_from_interval(x, a, b).clip(min=-1, max=1).reshape(1, -1)
        tx = np.cos(tn * np.arccos(xn)) * np.sqrt(2)
        tx[0] = 1
        return tx

    X, Y = np.atleast_1d(X, Y)
    original_shape = X.shape
    a, b, c, d = domain
    n = np.shape(coefficients)[1]

    tn = np.r_[0:n][:, None]
    tx1 = transform(tn, X.ravel(), a, b)
    tx2 = transform(tn, Y.ravel(), c, d)

    if use_meshgrid:  # eval on meshgrid points
        return np.dot(tx1.T, np.dot(coefficients, tx2)).T
    # scattered points
    val = np.sum(np.dot(tx1.T, coefficients) * tx2.T, axis=1)
    return val.reshape(original_shape)


if __name__ == '__main__':
    from wafo.testing import test_docstrings
    test_docstrings(__file__)