%GPR Trainable mapping for Gaussian Process regression
%
%     W = GPR(A,KERNEL,S_noise)
%     W = A*GPR([],KERNEL,S_noise)
%     W = A*GPR(KERNEL,S_noise)
%
%INPUT
%  A        Dataset used for training
%  KERNEL   Untrained mapping to compute kernel by A*(A*KERNEL)
%           during training, or B*(A*KERNEL) during evaluation with
%           dataset B
%  S_noise  Standard deviation of the noise 
%
%OUTPUT
%  W        Mapping: Gaussian Process regression
%
%DESCRIPTION
%Fit a Gaussian Process regressor on dataset A. For a nonlinear regressor,
%define kernel mapping KERNEL. For kernel definitions, have a look at
%proxm.m.
%
%SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
%DATASETS, MAPPINGS, SVMR, PROXM, LINEARR, TESTR, PLOTR

% Copyright: D.M.J. Tax, D.M.J.Tax@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

function y = gpr(varargin)

	argin = shiftargin(varargin,'prmapping',1);
	argin = shiftargin(argin,'scalar',2);
  argin = setdefaults(argin,[],proxm([],'p',1),1);
  
  if mapping_task(argin,'definition')
    y = define_mapping(argin,'untrained');
    y = setname(y,'Gaussian Proc. regression');
    
  elseif mapping_task(argin,'training')			% Train a mapping.
  
    [x,kernel,s_noise] = deal(argin{:});
    [n,d] = size(x);
    W.X = [+x ones(n,1)];
  % train:
    W.K = W.X*kernel;
    L = chol(+(W.X*W.K) + s_noise*s_noise*eye(n));
    W.w = L\(L'\gettargets(x));
  % store:
    W.kernel = kernel;
    y = prmapping(mfilename,'trained',W,1,d,1);
    y = setname(y,'Gaussian Proc. regression');
    
  else                                      % Evaluation
    
    [x,v] = deal(argin{1:2});
    W = getdata(v);
    out = [+x ones(size(x,1),1)]*W.K*W.w;
    y = setdat(x,out);
	
  end