function K = covPS(hyp, x, z, i)

% To be used with the GPML package: www.gaussianprocess.org/gpml/code
%
% PS covariance of the form:
%
% k(x,y) = cos(sin(2*pi*||x-y||/p) * exp( t*cos(2*pi*||x-y||/p) - t)
%
% where the hyperparameters are:
%
% hyp = [ log(p)    p gives the period of the spikes
%         t ]       t gives the shape of the valleys
%
% to see how the function looks try:
% x=linspace(1,500,500); p=10; t=10; y = cos(sin(x*2*pi/p)).*exp(t*cos(x*2*pi/p)-t); plot(y);                                   
%
% Copyright (c) by Daniel Preotiuc-Pietro, 2013-04-25.
%
% See also COVFUNCTIONS.M.

if nargin<2, K = '2'; return; end                  % report number of parameters
if nargin<3, z = []; end                                   % make sure, z exists
xeqz = numel(z)==0; dg = strcmp(z,'diag') && numel(z)>0;        % determine mode

n = size(x,1);
p = log(hyp(1));
t = hyp(2);

% precompute distances
if dg                                                               % vector kxx
  R = zeros(size(x,1),1);
else
  if xeqz                                                 % symmetric matrix Kxx
    R = sqrt(sq_dist(x'));
  else                                                   % cross covariances Kxz
    R = sqrt(sq_dist(x',z'));
  end
end

R=2*pi*R/exp(p);

if nargin<4                                                        % covariances
    K = cos(sin(R)).*exp(t*cos(R)-t);
else                                                               % derivatives
  if i==1
    K = sin(sin(R)).*cos(R).*R.*exp(t*(cos(R)-1))+cos(sin(R)).*exp(t*(cos(R)-1))*t.*R.*sin(R);
  elseif i==2
    K = cos(sin(R)).*exp(t*(cos(R)-1)).*(cos(R)-1);
  else
    error('Unknown hyperparameter')
  end
end

end
