""" Kernel class.

It provides Gram matrix computation and incomplete Cholesky decomposition
capabilities.

"""

from scipy.spatial.distance import pdist, cdist, squareform
from numpy import sum, sqrt, exp, dot, ones, array, zeros, argmax, \
                  hstack, newaxis, copy, argsort


class Kernel(object):
    """ Kernel class """

    def __init__(self, par=None):
        """ Initialization.

        Parameters
        ----------
        par : dictionary, optional
              Name of the kernel and its parameters (default is
              {'name': 'RBF','sigma': 1}). The name of the kernel comes
              from 'RBF', 'exponential', 'Cauchy', 'student', 'Matern3p2',
              'Matern5p2', 'polynomial', 'ratquadr' (rational quadratic),
              'invmquadr' (inverse multiquadr).

        Examples
        --------
        >>> from ite.cost.x_kernel import Kernel
        >>> k1 = Kernel({'name': 'RBF','sigma': 1})
        >>> k2 = Kernel({'name': 'exponential','sigma': 1})
        >>> k3 = Kernel({'name': 'Cauchy','sigma': 1})
        >>> k4 = Kernel({'name': 'student','d': 1})
        >>> k5 = Kernel({'name': 'Matern3p2','l': 1})
        >>> k6 = Kernel({'name': 'Matern5p2','l': 1})
        >>> k7 = Kernel({'name': 'polynomial','exponent': 2,'c': 1})
        >>> k8 = Kernel({'name': 'ratquadr','c': 1})
        >>> k9 = Kernel({'name': 'invmquadr','c': 1})

        from numpy.random import rand
        num_of_samples, dim = 5, 2
        y1, y2 = rand(num_of_samples, dim), rand(num_of_samples+1, dim)
        y1b = rand(num_of_samples, dim)
        k1.gram_matrix1(y1)
        k1.gram_matrix2(y1, y2)
        k1.sum(y1,y1b)
        k1.gram_matrix_diagonal(y1)

        """

        # if par is None:
        #     par = {'name': 'RBF', 'sigma': 0.01}
        if par is None:
            par = {'name': 'RBF', 'sigma': 1}
        # if par is None:
        #    par = {'name': 'exponential', 'sigma': 1}
        # if par is None:
        #     par = {'name': 'Cauchy', 'sigma': 1}
        # if par is None:
        #     par = {'name': 'student', 'd': 1}
        # if par is None:
        #     par = {'name': 'Matern3p2', 'l': 1}
        # if par is None:
        #     par = {'name': 'Matern5p2', 'l': 1}
        # if par is None:
        #     par = {'name': 'polynomial','exponent': 2, 'c': 1}
        # if par is None:
        #     par = {'name': 'polynomial', 'exponent': 3, 'c': 1}
        # if par is None:
        #     par = {'name': 'ratquadr', 'c': 1}
        # if par is None:
        #     par = {'name': 'invmquadr', 'c': 1}

        # name:
        name = par['name']
        self.name = name

        # other attributes:
        if name == 'RBF' or name == 'exponential' or name == 'Cauchy':
            self.sigma = par['sigma']
        elif name == 'student':
            self.d = par['d']
        elif name == 'Matern3p2' or name == 'Matern5p2':
            self.l = par['l']
        elif name == 'polynomial':
            self.c = par['c']
            self.exponent = par['exponent']
        elif name == 'ratquadr' or name == 'invmquadr':
            self.c = par['c']
        else:
            raise Exception('kernel=?')

    def __str__(self):
        """ String representation of the kernel.

        Examples
        --------
        print(kernel)

        """

        return ''.join((self.__class__.__name__, ' -> ',
                        str(self.__dict__)))

    def gram_matrix1(self, y):
        """  Compute the Gram matrix = [k(y[i,:],y[j,:])]; i, j: running.

        Parameters
        ----------
        y : (number of samples, dimension)-ndarray
            One row of y corresponds to one sample.

        Returns
        -------
        g : ndarray.
            Gram matrix of y.

        Examples
        --------
        g = k.gram_matrix1(y)

        """

        if self.name == 'RBF':
            sigma = self.sigma
            g = squareform(pdist(y))
            g = exp(-g ** 2 / (2 * sigma ** 2))
        elif self.name == 'exponential':
            sigma = self.sigma
            g = squareform(pdist(y))
            g = exp(-g / (2 * sigma ** 2))
        elif self.name == 'Cauchy':
            sigma = self.sigma
            g = squareform(pdist(y))
            g = 1 / (1 + g ** 2 / sigma ** 2)
        elif self.name == 'student':
            d = self.d
            g = squareform(pdist(y))
            g = 1 / (1 + g ** d)
        elif self.name == 'Matern3p2':
            l = self.l