function [model, loo_err, loo_pred] = kls_train_multi(X,Y,model,bias)
% KLS_TRAIN_MULTI    Train a Kernel Least Square for Multiclass (One-vs-all)
% Input: X training data
%   [MODEL]                    = KLS_TRAIN(X,Y,MODEL)
%   [MODEL, LOO_ERR]           = KLS_TRAIN(X,Y,MODEL)
%   [MODEL, LOO_ERR, LOO_PRED] = KLS_TRAIN(X,Y,MODEL)
% Input: K training matrix
%   [MODEL]                    = KLS_TRAIN(K,Y,MODEL)
%   [MODEL, LOO_ERR]           = KLS_TRAIN(K,Y,MODEL)
%   [MODEL, LOO_ERR, LOO_PRED] = KLS_TRAIN(K,Y,MODEL)
%
%    This file is part of the DOGMA library for Matlab.
%    Copyright (C) 2009, Francesco Orabona
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.
%
%    Contact the author: bremen79 [at] gmail.com

if nargin<4
  bias = 1; 
end
n = length(Y);

tmp=model;
model=cell(1,tmp.n_cla);
for i=1:tmp.n_cla
    model{i}=tmp;
    if isfield(tmp, 'ker')==0
       model{i}.SV=X;
       model{i}.S=1:n;
    end
    model{i}.SV=X;
end

if isfield(model{1}, 'ker')==0
  model{1}.K=X;
else
  model{1}.K=feval(model{1}.ker,X,1:n,1:n,model{1}.kerparam);
end

if bias
    model{1}.K=[model{1}.K ones(n,1); ones(1,n) 0];

    if numel(model{1}.C)==1
        id=1/model{1}.C*eye(n+1); 
        id(end,end)=0;
    else
        id=diag(1./model{1}.C); 
        id(end+1,end+1)=0;
    end
    G=pinv(model{1}.K+id);
    dd=diag(G);

    ytmp=zeros(size(Y));
    loo_pred=zeros(numel(Y),model{1}.n_cla);
    for i=1:model{1}.n_cla
        idx_pos=find(Y==i);
        idx_neg=find(Y~=i);
        ytmp(idx_pos')=1;
        ytmp(idx_neg')=-1;
        x=G*[ytmp'; 0];

        model{i}.beta = x(1:end-1);
        model{i}.b = x(end);

        loo_pred(:,i)=ytmp'-model{i}.beta./dd(1:end-1);
    end
else
    if numel(model{1}.C)==1
        id=1/model{1}.C*eye(n); 
    else
        id=diag(1./model{1}.C); 
    end
    G=pinv(model{1}.K+id);
    dd=diag(G);
 
    ytmp=zeros(size(Y));
    loo_pred=zeros(numel(Y),numel(model{1}.n_cla));
    for i=1:model{1}.n_cla
        idx_pos=find(Y==i);
        idx_neg=find(Y~=i);
        ytmp(idx_pos')=1;
        ytmp(idx_neg')=-1;
        x=G*ytmp';

        model{i}.beta = x;
        model{i}.b = 0;

        loo_pred(:,i)=ytmp'-model{i}.beta./dd;
    end
end

[tmp,pred]=max(loo_pred,[],2);
loo_err=1-numel(find(pred==Y'))/numel(Y);
loo_pred=pred;
