function model = gpCovGradsTest(model)

% GPCOVGRADSTEST Test the gradients of the likelihood wrt the covariance.
%
%	Description:
%
%	MODEL = GPCOVGRADSTEST(MODEL) tests the gradients of the covariance
%	to ensure they are correct.
%	 Returns:
%	  MODEL - the model that was tested.
%	 Arguments:
%	  MODEL - the model to be tested.
%	
%
%	See also
%	GPCREATE, GPCOVGRADS


%	Copyright (c) 2006, 2009 Neil D. Lawrence
% 	gpCovGradsTest.m CVS version 1.1
% 	gpCovGradsTest.m SVN version 178
% 	last update 2009-01-13T14:14:27.000000Z

% WARNING --- this isn't testing g_Lambda in gpCovGrads
  
changeVal = 1e-6;
switch model.approx
 case 'ftc'
  
 case {'dtc', 'dtcvar', 'fitc', 'pitc'}
  for i =1 :size(model.K_uu, 1)
    for j=1:i
      origK = model.K_uu(i, j);
      model.K_uu(i, j) = origK + changeVal;
      model.K_uu(j, i) = model.K_uu(i, j);
      [model.invK_uu, U] = pdinv(model.K_uu);
      model.logDetK_uu = logdet(model.K_uu, U);
      model = gpUpdateAD(model);
      objPlus = gpLogLikelihood(model);
      model.K_uu(i, j) = origK - changeVal;
      model.K_uu(j, i) = model.K_uu(i, j);
      [model.invK_uu, U] = pdinv(model.K_uu);
      model.logDetK_uu = logdet(model.K_uu, U);
      model = gpUpdateAD(model);
      objMinus = gpLogLikelihood(model);
      diffsK_uu(i, j) = (objPlus - objMinus)/(2*changeVal);
      diffsK_uu(j, i) = diffsK_uu(i, j);
      model.K_uu(i, j) = origK;
      model.K_uu(j, i) = origK;
      [model.invK_uu, U] = pdinv(model.K_uu);
      model.logDetK_uu = logdet(model.K_uu, U);
      model = gpUpdateAD(model);
    end
  end
  for i=1:size(model.K_uf, 1)
    for j=1:size(model.K_uf, 2)
      origK = model.K_uf(i, j);
      model.K_uf(i, j) = origK + changeVal;
      model = gpUpdateAD(model);
      objPlus = gpLogLikelihood(model);
      model.K_uf(i, j) = origK - changeVal;
      model = gpUpdateAD(model);
      objMinus = gpLogLikelihood(model);
      diffsK_uf(i, j) = (objPlus - objMinus)/(2*changeVal);
      model.K_uf(i, j) = origK;
      model = gpUpdateAD(model);
    end
  end
  
  [gK_uu, gK_uf, g_Lambda] = gpCovGrads(model, model.m);
  
  gK_uuMaxDiff = max(max(abs(2*(gK_uu-diag(diag(gK_uu))) ...
                             + diag(diag(gK_uu)) ...
                             - diffsK_uu)));
  gK_ufMaxDiff = max(max(abs(gK_uf - diffsK_uf)));
  
  fprintf('K_uu grad max diff %2.4f\n', gK_uuMaxDiff);
  if gK_uuMaxDiff > 1e-4
    disp(2*(gK_uu-diag(diag(gK_uu))) ...
         + diag(diag(gK_uu)) ...
         - diffsK_uu);
  end
  fprintf('K_uf grad max diff %2.4f\n', gK_ufMaxDiff);
  if gK_ufMaxDiff > 1e-4
    disp(gK_uf - diffsK_uf)
  end
end

