function [W,H,err,iMSV] = sv_nmf(X,nu,ker_string,skip_norm,comp_h,V)

% svm_nmf - compute NMF on data matrix X such that X ~ WH using the single-class
% SVM technique, after (Essid, 2012)
%
% SYNTAX
% [W,H,err,iMSV] = sv_nmf(X,nu)
% [W,H,err,iMSV] = sv_nmf(X,nu,ker_string,skip_norm,comp_h,V)
%
% ARGUMENTS
% X - Data matrix where each column is an observation (data point) except
% when a custom kernel is used (with Libsvm option '-t 4'), where X is
% rather the Gram matrix as expected by Libsvm
%
% nu - Libsvm parameter nu (see Libsvm doc)
%     
% ker_string - Optional argument; string to specify the kernel to be used 
% with Libsvm (see Libsvm doc); default is '-t 0' i.e. linear kernel
%
% skip_norm - Optional argument; if set to zero the data will be first
% normalised to have unit euclidean norm; default is 0
%
% comp_h - Optional argument; if set to zero the activations H will not be
% computed; default is 1
%
% V - Optional argument; the original data matrix, only useful when using a
% custom kernel when X is the Gram matrix
%
% OUTPUT
% W - "Basis vector" matrix
% H - Activations matrix 
% err - Approximation error
% iMSV - Indices of basis vectors in columns of original data matrix X
%
% WARNING 
% This program has not yet been fully tested and optimised. A number of
% checks are yet to be implemented. For instance nothing is done here to
% make sure that the data is properly normalised (to unit-norm, see (Essid,
% 2012) for explanations), though this is done by default in the linear
% case (see Line 81).
% It requires the Matlab optimisation toolbox and libsvm-mat (expected to
% be found in the Matlab path)
%
% (Essid, 2012) S. Essid, "A SINGLE-CLASS SVM BASED ALGORITHM FOR COMPUTING AN 
% IDENTIFIABLE NMF," in IEEE International Conference on Acoustics, Speech 
% and Signal Processing, 2012.
%  
% Copyright 2012 Slim Essid, Telecom ParisTech
%

%  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/>.


%addpath('/libsvm-mat-3.0-1'); % of fast addtitive svms package...

% -- Set default values ---------------------------------------------------

if nargin<5
    comp_h = 1;
    if nargin<4
        skip_norm = 0;
        if nargin<3
            ker_string = '-t 0';
        end
    end
end

% -- Normalise data -------------------------------------------------------
% Use only with linear Kernel, improve by adding checks
if ~skip_norm    
    scl = sqrt(sum(X.^2,1));    
    X = X ./ repmat(scl,size(X,1),1);    
end

if nargin<6
    V = X;
end

% -- Find W using 1-class SVM ---------------------------------------------
fprintf('Finding W...');

if strfind(ker_string,'-t 4') % costum pre-computed kernel
    SVM = svmtrain(ones(size(X,1), 1), X, ...
        ['-m 1000 -s 2 ' ker_string ' -n '  num2str(nu)]);
else
    SVM = svmtrain(ones(size(X',1), 1), X', ...
        ['-m 1000 -s 2 ' ker_string ' -n '  num2str(nu)]);
end

fprintf('%d support vectors\n', SVM.totalSV);

% -- Get W columns margin support vectors
W = []; H = []; err = [];
iNMSV = []; % outlier indices in X
iMSV = []; % margin support vector indices in X

for ns = 1:SVM.totalSV    
    if SVM.sv_coef(ns) == 1, % outlier
        if strfind(ker_string, '-t 4')
            iNMSV = [iNMSV full(SVM.SVs(ns))]; % for custom kernel case
        else
            iNMSV = [iNMSV getSvRelIndex(full(SVM.SVs(ns,:)), X')];
        end
        continue,
    end
        
    if strfind(ker_string, '-t 4')
        iMSV = [iMSV  full(SVM.SVs(ns))];
        W = [W V(:,full(SVM.SVs(ns)))];
    else
        x = full(SVM.SVs(ns,:));
        W = [W x'];
        iMSV = [iMSV getSvRelIndex(x, X')];
    end
end
fprintf('%d basis vectors  (margin support vectors)\n',size(W,2));
K = length(iMSV);

% -- Compute H  -----------------------------------------------------------
if comp_h
    fprintf('Computing H...');
  
    H = zeros(size(W,2),size(V,2));
    err = zeros(1,size(V,2));
    
    if strfind(ker_string, '-t 0') % linear case        
        % nonnegative LS
        for i = 1:size(V,2)
            if mod(i,1000)==0, fprintf('.'); end
            [H(:,i),err(i)] = lsqnonneg(W,V(:,i));
        end
        fprintf('\n');
    
    else % kernel-based case
        
        % compute Gram matrix on support vectors (needed in optimiation)
        % to be improved by getting it out of SVM optimisation !!        
        if strfind(ker_string,'-t 4') % gram matrix given
            Gw = X(iMSV,iMSV+1);
        else
            Gw = zeros(K,K);
            for k = 1:K
                for l = k:K
                    Gw(k,l) = eval_ker(W(:,k), W(:,l), ker_string);
                end
            end
             Gw = Gw + triu(Gw,1)'; % symmetrise W
        end
        
        % solve optimisation problem for H using the Matlab opt toolbox
        % to be improved using a faster opt technique !!
        options = optimset('LargeScale','off','GradObj', 'on', 'Hessian', 'off', ...
            'DerivativeCheck', 'off', 'Diagnostics', 'off','Display','off',...
            'MaxFunEvals',10000,'TolX',1e-6,'TolFun',1e-10,'Algorithm','active-set');
        
        for i = 1:size(V,2)            
            vi = V(:,i);            
            h0 = abs(randn(K,1)) + ones(K,1);
            
            h_objfun = @(hi)prm_h_objfun(hi, W, Gw, vi, ker_string, X, iMSV, i);
            
            [H(:,i),fval,exitflag,output,lambda,grad] = ...
                fmincon(h_objfun, h0, [], [], [], [], zeros(K,1), inf, ...
                [], options );
        end        
    end
    
    if ~skip_norm 
        H = H*diag(scl);
    end
end


function iSV = getSvRelIndex(SV, X)

obs_score = sum(X==repmat(SV,size(X,1), 1), 2);
[valeur, iSV] = max(full(obs_score));
if (valeur~=length(SV)),
    error('--> Support vector not found in database!');
    %iSV = 0;
end

function [f, g, Hs] = prm_h_objfun(hi, W, Gw, vi, ker_string, X, iMSV, i)

K = size(W,2);

% intermediate var K(vi, vsk)
k_ik = zeros(K,1);
for k = 1:K
    if strfind(ker_string,'-t 4')
        k_ik(k) = X(iMSV(k), i);
    else
        k_ik(k) =  eval_ker(W(:,k), vi, ker_string);
    end
end

% objective function
f = 1 + hi'*Gw*hi - 2 * sum( hi.*k_ik, 1);

if nargout> 1 % gradient
    g = zeros(K,1);
    for k = 1:K         
        g(k) = 2*( sum( hi.*Gw(:,k), 1 ) - k_ik(k) );    
    end
    if nargout > 2 % Hessian
        Hs = 2*Gw;
    end
end

    
function kxy = eval_ker(x, y, ker_string)

switch ker_string
    case '-t 5' % intersection kernel
        kxy = sum( min([x,y], [], 2), 1);
    case '-t 6' % chi-squared kernel
        kxy = sum( 2*(x.*y)./(x+y), 1);
    case '-t 0' % linear (for validation purposes)
        kxy = x'*y;
    otherwise
        ik = strfind(ker_string,'-t 2');
        if ik % Gaussian
            q = str2double(ker_string(ik+7:end));            
            kxy = exp(- 2*q *(1-x'*y)); % supposing unit-norm vectors
        else
            error('Objective function for H not defined with this kernel');
        end
end