# Source code for compmech.interpolate

    r"""
    Interpolate (:mod:`compmech.interpolate`)
    ==================================================

    .. currentmodule:: compmech.interpolate

    This module includes some interpolation utilities that will be used in other
    modules.

    .. autofunction:: inv_weighted

    .. autofunction:: interp

    """
    from collections import Iterable

    import numpy as np

    from compmech.logger import msg, warn

    def inv_weighted(data, mesh, num_sub, col, ncp=5, power_parameter=2):
        r"""Interpolates the values taken at one group of points into
        another using an inverse-weighted algorithm

        In the inverse-weighted algorithm a number of `n_{CP}` measured points
        of the input parameter ``data`` that are closest to a given node in
        the input parameter ``mesh`` are found and the imperfection value of
        this node (represented by the normal displacement `{w_0}_{node}`) is
        calculated as follows:

        .. math::
            {w_0}_{node} = \frac{\sum_{i}^{n_{CP}}{{w_0}_i\frac{1}{w_i}}}
                                {\sum_{i}^{n_{CP}}{\frac{1}{w_i}}}

        where `w_i` is the imperfection at each measured point, calculated as:

        .. math::
            w_i = \left[(x_{node}-x_i)^2+(y_{node}-y_i)^2+(z_{node}-y_i)^2
                  \right]^p

        with `p` being a power parameter that when increased will increase the
        relative influence of a closest point.

        Parameters
        ----------
        data : numpy.ndarray, shape (N, ndim+1)
            The data or an array containing the imperfection file. The values
            to be interpolated must be in the last column.
        mesh : numpy.ndarray, shape (M, ndim)
            The new coordinates where the values will be interpolated to.
        num_sub : int
            The number of sub-sets used during the interpolation. The points
            are divided in sub-sets to increase the algorithm's efficiency.
        col : int
            The index of the column to be used in order to divide the data
            in sub-sets. Note that the first column index is ``0``.
        ncp : int, optional
            Number of closest points used in the inverse-weighted interpolation.
        power_parameter : float, optional
            Power of inverse weighted interpolation function.

        Returns
        -------
        ans : numpy.ndarray
            A 1-D array with the interpolated values. The size of this array
            is ``mesh.shape[0]``.

        """
        if mesh.shape[1] != data.shape[1]-1:
            raise ValueError('Invalid input: mesh.shape[1] != data.shape[1]')

        msg('Interpolating... ')
        num_sub = int(num_sub)
        mesh_size = mesh.shape[0]

        # memory control
        mem_limit = 1024*1024*1024*8*2    # 2 GB
        mem_entries = int(mem_limit / 64) # if float64 is used
        sec_size = int(mesh_size/num_sub)
        while sec_size**2*10 > mem_entries:
            num_sub +=1
            sec_size = int(mesh_size/num_sub)
            if sec_size**2*10 <= mem_entries:
                warn('New num_sub: {0}'.format(int(mesh_size/float(sec_size))))
                break

        mesh_seq = np.arange(mesh.shape[0])

        mesh_argsort = np.argsort(mesh[:, col])
        mesh_seq = mesh_seq[mesh_argsort]
        back_argsort = np.argsort(mesh_seq)

        mesh = np.asarray(mesh[mesh_argsort], order='F')

        length = mesh[:, col].max() - mesh[:, col].min()

        data = np.asarray(data[np.argsort(data[:, col])], order='F')

        ans = np.zeros(mesh.shape[0], dtype=mesh.dtype)

        # max_num_limits defines how many times the log will print
        # "processed ... out of ... entries"
        max_num_limits = 10
        for den in range(max_num_limits, 0, -1):
            if num_sub % den == 0:
                limit = int(num_sub/den)
                break

        for i in range(num_sub+1):
            i_inf = sec_size*i
            i_sup = sec_size*(i+1)

            if i % limit == 0:
                msg('\t processed {0:7d} out of {1:7d} entries'.format(
                      min(i_sup, mesh_size), mesh_size))
            sub_mesh = mesh[i_inf : i_sup]
            if not np.any(sub_mesh):
                continue
            inf = sub_mesh[:, col].min()
            sup = sub_mesh[:, col].max()

            tol = 0.03
            if i == 0 or i == num_sub:
                tol = 0.06

            while True:
                cond1 = data[:, col] >= inf - tol*length
                cond2 = data[:, col] <= sup + tol*length
                cond = np.all(np.array((cond1, cond2)), axis=0)
                sub_data = data[cond]
                if not np.any(sub_data):
                    tol += 0.01
                else:
                    break

            dist = np.subtract.outer(sub_mesh[:, 0], sub_data[:, 0])**2
            for j in range(1, sub_mesh.shape[1]):
                dist += np.subtract.outer(sub_mesh[:, j], sub_data[:, j])**2
            asort = np.argsort(dist, axis=1)
            lenn = sub_mesh.shape[0]
            lenp = sub_data.shape[0]
            asort_mesh = asort + np.meshgrid(np.arange(lenn)*lenp,
                                             np.arange(lenp))[0].transpose()
            # getting the distance of the closest points
            dist_cp = np.take(dist, asort_mesh[:, :ncp])
            # avoiding division by zero
            dist_cp[(dist_cp==0)] == 1.e-12
            # fetching the imperfection of the sub-data
            imp = sub_data[:, -1]
            # taking only the imperfection of the closest points
            imp_cp = np.take(imp, asort[:, :ncp])
            # weight calculation
            total_weight = np.sum(1./(dist_cp**power_parameter), axis=1)
            weight = 1./(dist_cp**power_parameter)
            # computing the new imp
            imp_new = np.sum(imp_cp*weight, axis=1)/total_weight
            # updating the answer array
            ans[i_inf : i_sup] = imp_new

        ans = ans[back_argsort]

        msg('Interpolation completed!')

        return ans

    def interp(x, xp, fp, left=None, right=None, period=None):
        """
        One-dimensional linear interpolation

        Returns the one-dimensional piecewise linear interpolant to a function
        with given values at discrete data-points.

        .. note:: This function has been incorporated in NumPy >= 1.10.0 and will be soon
                  removed from here.

        Parameters
        ----------
        x : array_like
            The x-coordinates of the interpolated values.
        xp : 1-D sequence of floats
            The x-coordinates of the data points, must be increasing if argument
            ``period`` is not specified. Otherwise, ``xp`` is internally sorted
            after normalizing the periodic boundaries with ``xp = xp % period``.
        fp : 1-D sequence of floats
            The y-coordinates of the data points, same length as ``xp``.
        left : float, optional
            Value to return for ``x < xp[0]``, default is ``fp[0]``.
        right : float, optional
            Value to return for ``x > xp[-1]``, default is ``fp[-1]``.
        period : float, optional
            A period for the x-coordinates. This parameter allows the proper
            interpolation of angular x-coordinates. Parameters ``left`` and
            ``right`` are ignored if ``period`` is specified.

        Returns
        -------
        y : {float, ndarray}
            The interpolated values, same shape as ``x``.

        Raises
        ------
        ValueError
            If ``xp`` and ``fp`` have different length
            If ``xp`` or ``fp`` are not 1-D sequences
            If ``period==0``

        Notes
        -----
        Does not check that the x-coordinate sequence ``xp`` is increasing.
        If ``xp`` is not increasing, the results are nonsense.
        A simple check for increasing is::

            np.all(np.diff(xp) > 0)

        Examples
        --------
        >>> xp = [1, 2, 3]
        >>> fp = [3, 2, 0]
        >>> interp(2.5, xp, fp)
        1.0
        >>> interp([0, 1, 1.5, 2.72, 3.14], xp, fp)
        array([ 3. ,  3. ,  2.5 ,  0.56,  0. ])
        >>> UNDEF = -99.0
        >>> interp(3.14, xp, fp, right=UNDEF)
        -99.0

        Plot an interpolant to the sine function:

        >>> x = np.linspace(0, 2*np.pi, 10)
        >>> y = np.sin(x)
        >>> xvals = np.linspace(0, 2*np.pi, 50)
        >>> yinterp = interp(xvals, x, y)
        >>> import matplotlib.pyplot as plt
        >>> plt.plot(x, y, 'o')
        [<matplotlib.lines.Line2D object at 0x...>]
        >>> plt.plot(xvals, yinterp, '-x')
        [<matplotlib.lines.Line2D object at 0x...>]
        >>> plt.show()

        Interpolation with periodic x-coordinates:

        >>> x = [-180, -170, -185, 185, -10, -5, 0, 365]
        >>> xp = [190, -190, 350, -350]
        >>> fp = [5, 10, 3, 4]
        >>> interp(x, xp, fp, period=360)
        array([7.5, 5., 8.75, 6.25, 3., 3.25, 3.5, 3.75])

        """
        if period is None:
            return np.interp(x, xp, fp, left, right)
        else:
            if period==0:
                raise ValueError('Argument `period` must be a non-zero value')
            period = abs(period)
            if not isinstance(x, Iterable):
                x = [x]
            x = np.asarray(x)
            xp = np.asarray(xp)
            fp = np.asarray(fp)
            if xp.ndim != 1 or fp.ndim != 1:
                raise ValueError('Data points must be 1-D sequences')
            if xp.shape[0] != fp.shape[0]:
                raise ValueError('Inputs `xp` and `fp` must have the same shape')
            # eliminating discontinuity between periods
            x = x % period
            xp = xp % period
            asort_xp = np.argsort(xp)
            xp = xp[asort_xp]
            fp = fp[asort_xp]
            xp = np.concatenate((xp[-1:]-period, xp, xp[0:1]+period))
            fp = np.concatenate((fp[-1:], fp, fp[0:1]))
            return np.interp(x, xp, fp)

    if __name__=='__main__':
        a = np.array([[1.1, 1.2, 10],
                      [1.2, 1.2, 10],
                      [1.3, 1.3, 10],
                      [1.4, 1.3, 10],
                      [1.5, 1.3, 10],
                      [2.6, 2.3, 5],
                      [2.7, 2.3, 5],
                      [2.6, 2.1, 5],
                      [2.7, 2.1, 5],
                      [2.8, 2.2, 5],
                      [2.8, 2.2, 5],
                      [5.6, 5.3, 20],
                      [5.7, 5.3, 20],
                      [5.6, 5.1, 20],
                      [5.7, 5.1, 20],
                      [5.8, 5.2, 20],
                      [5.8, 5.2, 20]])

        b = np.array([[1., 1.],
                      [2., 2.],
                      [4., 4.],
                      [5., 5.]])

        print inv_weighted(a, b, num_sub=1, col=1, ncp=10)