# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
# Licensed under the BSD 3-clause license (see LICENSE.txt)


from kernpart import Kernpart
import numpy as np

class PeriodicSpikes(Kernpart):
    """
    periodic spikes kernel

    .. math::

       k(r) = cos(sin(2*pi*r/p) * exp( s*cos(2*pi*r/p) - s) where r^2 = (x-y)^2
       
    :param input_dim: the number of input dimensions
    :type input_dim: int (input_dim=1 is the only value currently supported)
    :param shape: the shape parameter s, should be >1
    :type shape: float 
    :param period: the period p
    :type period: float
    :rtype: Kernpart object

    """
    def __init__(self,input_dim,shape=1.,period=1.):
        assert input_dim == 1, "For this kernel we assume input_dim=1"
        self.input_dim = input_dim
        self.num_params = 2
        self.name = 'ps'
        self.shape = shape
        self.period = period

    def _get_params(self):
        return np.hstack((self.shape,self.period))

    def _set_params(self,x):
        self.shape = x[0]
        self.period = x[1]

    def _get_param_names(self):
        return ['shape','period']

    def K(self,X,X2,target):
        if X2 is None: X2 = X
        pcomp = 2.*np.pi*(X-X2.T)/self.period
        target += np.cos(np.sin(pcomp))*np.exp(self.shape*(np.cos(pcomp)-1))

    def Kdiag(self,X,target):
        target += 1

    def dK_dtheta(self,dL_dK,X,X2,target):
        if X2 is None: X2 = X
	pcomp=2.*np.pi*(X-X2.T)/self.period
        dshape = np.cos(np.sin(pcomp))*np.exp(self.shape*(np.cos(pcomp)-1))*(np.cos(pcomp)-1)
        dper = sin(sin(pcomp))*cos(pcomp)*pcomp*np.exp(self.shape*(cos(pcomp)-1))+cos(sin(pcomp))*np.exp(self.shape*(np.cos(pcomp)-1))*self.shape*pcomp*np.sin(pcomp);

        target[0] += np.sum(dshape*dL_dK)
        target[1] += np.sum(dper*dL_dK)

    def dKdiag_dtheta(self,dL_dKdiag,X,target):
# !!! these are not implemented
        target[0] += np.sum(dL_dKdiag)
        # here self.lengthscale and self.power have no influence on Kdiag so target[1:] are unchanged

    def dK_dX(self,dL_dK,X,X2,target):
# !!! these are not implemented
        """derivative of the covariance matrix with respect to X."""
#        if X2 is None:
#            dist2 = np.square((X-X.T)/self.lengthscale)
#            dX = -2.*self.variance*self.power * (X-X.T)/self.lengthscale**2 *  (1 + dist2/2./self.lengthscale)**(-self.power-1)
#        else:
#            dist2 = np.square((X-X2.T)/self.lengthscale)
#            dX = -self.variance*self.power * (X-X2.T)/self.lengthscale**2 *  (1 + dist2/2./self.lengthscale)**(-self.power-1)
        target += np.sum(dL_dK*dX,1)[:,np.newaxis]

    def dKdiag_dX(self,dL_dKdiag,X,target):
        pass