function [x,stat] = quad_prog_gsmo(H, f, a, b, LB, UB, x0, opt)
% solve m number of quadratic programming (QP) problems using
% vectorization by the generalized SMO algorithm.
% Note, 
% 1. all QPs have the same constratints 
% 2. a = ones(n, 1) 2012-12-18 for speed concern  
% 3. b = 1 2012-12-18
%
% Synopsis:
%  [x,stat] = libqp_gsmo(H,f,a,b,LB,UB) to do
%  [x,stat] = libqp_gsmo(H,f,a,b,LB,UB,x0) to do
%  [x,stat] = libqp_gsmo(H,f,a,b,LB,UB,[],opt) to do 
%  [x,stat] = libqp_gsmo(H,f,a,b,LB,UB,x0,opt) supported only now
% 
% Description:
%  This function implements the Generalized SMO algorithm which solves
%  the following convex quadratic programming task:
%
%   min 0.5*x'*H*x + f'*x  
%    x                                      
%
%   subject to   a'*x = b 
%                LB(i) <= x(i) <= UB(i)   for all i=1:n
%
% Reference:
%  S.-S. Keerthi, E.G.Gilbert. Convergence of a Generalized SMO Algorithm for SVM 
%   Classifier Design. Technical Report CD-00-01, Control Division, Dept. of Mechanical 
%   and Production Engineering, National University of Singapore, 2000. 
%   http://citeseer.ist.psu.edu/keerthi00convergence.html   
%
% Input:
%  H [m x (n*n)] m weighting matrix in lexicographical order
%  f [m x n] first order term (p)
%  a [m x n] Vector which must not contain zero entries.
%  b [1 x 1] Scalar.
%  LB [m x n] Lower bound; -inf is allowed.
%  UB [m x n] Upper bound; inf is allowed.
%
% Optional inputs: 
%  x0 [m x n] Initial solution.
%
%  options [struct] 
%    .TolKKT [1 x 1] Determines relaxed KKT conditions (default TolKKT=0.001);
%       it correspondes to $\tau$ in Keerthi's paper.
%    .verb [1 x 1] if > 0 then prints info every verb-th iterations (default 0)
%    .MaxIter [1 x 1] Maximal number of iterations (default inf).
%  
% Output:
%  x [m x n] Solution vector.
%  stat [struct]
%   .QP [m x 1] Primal objective value.
%   .exitflag [n x 1] Indicates which stopping condition was used:
%      nIter >= MaxIter                  ->  exitflag = 0
%      relaxed KKT conditions satisfied  ->  exitflag = 4  
%   .nIter [1x1] Number of iterations.
%
% Example: to fix for m QPs 
%  n=50; X = rand(1,n); H = diag(X); f = -10*rand(n,1);  
%  a=ones(1,n); b=1; tmp = rand(n,1); LB = tmp-1; UB = tmp+1;
%
%  nup = 1; L=n; ind = 1;
%  x0 = zeros(nup, L); 
%  xa = zeros(nup, 1); 
%  i = 0;
%  while any(sum(x0.* a(ind, :), 2)) ~= b
%    i = i + 1;
%    if i > length(a),
%      error('No feasible solution exists.');
%    end
%    x0(:, i) = min(UB(ind,i),max(LB(ind,i),(b-xa)./ a(ind,i)));
%    xa = xa + x0(:,i) .* a(ind, i);
% end
% tic; [x1,stat1] = libqp_gsmo(H,f,a,b,LB,UB,x0); fval1=stat1.QP, toc
% todo 2013-11-17: tic; [x2,stat2] = quad_prog_gsmo(X,f',a,b,LB',UB',x0); toc
%
% todo 2013-11-17: tic; [x3,stat3] = quad_prog_gsmo([X;X],[f,f]',[a;a], b, [LB,LB]', [UB,UB]', [x0;x0]); toc
%
% Yong Long, 2012-12-13


if nargin < 7, error('At least seven input arguments are required.'); end

if nargin < 8, opt = []; end
if ~isfield(opt,'TolKKT'), opt.TolKKT = 0.001; end
if ~isfield(opt,'MaxIter'), opt.MaxIter = inf; end
if ~isfield(opt,'verb'), opt.verb = 0; end

[m n] = size(x0);
[x, stat] = gsmo_solver(H, f, a, b, LB, UB, x0, opt);

% function value, only needed for debugging
x1 = repmat(x, [1 n]); %(x1 x2 x3 x1 x2 x3 ...)[m n*n]
x2 = reshape(repmat(x, [n 1]), [m n*n]); % (x1 x1... x2 x2....x3 x3...)
stat.QP = sum(x1 .* x2 .* H, 2); %[m 1]
stat.QP = 0.5 * stat.QP + sum(f.*x, 2);

%   
% gsmo_solver
%
% 
function [x, stat] = gsmo_solver(H, f, a, b, LB, UB, x0, opt)

[m n] = size(x0);

% Initialization 
stat.niter = 0;
stat.exitflag = 100;
x = x0;
xtmp = x;

Htmp = H; ftmp = f; LBtmp = LB; UBtmp = UB;
xInd = ones(m, 1) == 1;

% Main optimization loop  
while(stat.exitflag == 100)
    clear gradQPtmp
    stat.niter = stat.niter + 1;
    
    xtmp1 = repmat(xtmp, [1 n]); %[m n*n]
    gradQPtmp = reshape(sum(reshape(xtmp1 .* Htmp, [m n n]), 2), [m n]); 
    gradQPtmp = gradQPtmp + ftmp; % [m x n]
    
    
    F = gradQPtmp; % [m x n]
    
    % find the most violating pair of variables
    % ie. minimum F_i in I_up 
    % and maximum F_i in I_low
    I_0 = LBtmp < xtmp & xtmp < UBtmp;
    I_1 = xtmp == LBtmp;
    I_3 = xtmp == UBtmp;
    
    I_up = I_0 | I_1; 
    upind = I_up == 0;
    I_low = I_0 | I_3; 
    lowind = I_low == 0;
    
    Fup = F.*I_up; Fup(upind) = inf; %[m x n]
    Flow = F.*I_low; Flow(lowind) = -inf; %[m x n]
    
    [minF_up, u] = min(Fup, [], 2); %[m x 1]
    [maxF_low, v] = max(Flow, [], 2); %[m x 1]
  
    
    % check KKT conditions 
    updateInd = (maxF_low - minF_up) > opt.TolKKT;% elements needs updates
          
    xInd(xInd == 1) = updateInd;
    
    if any(updateInd ~= 0) 
        xtmp = xtmp(updateInd,:);
        ftmp = ftmp(updateInd, :);
        Htmp = Htmp(updateInd, :);
        LBtmp = LBtmp(updateInd, :);
        UBtmp = UBtmp(updateInd, :); 
        gradQPtmp = gradQPtmp(updateInd, :); 
        
        m = size(xtmp, 1);
        mm = 1:1:m;
        uu = (u(updateInd)' - 1) * m + mm;
        vv = (v(updateInd)' - 1) * m + mm;
        uu3 = ((u(updateInd)' - 1) * (n + 1) ) * m + mm;
        vv3 = ((v(updateInd)' - 1) * (n + 1) ) * m + mm;
        uv3 = ((u(updateInd)' - 1) * n + v(updateInd)' - 1) * m + mm;
        
        % SMO update of the most violating pair
        tau = (gradQPtmp(vv) - gradQPtmp(uu)) ./ (Htmp(uu3) + Htmp(vv3) - 2*Htmp(uv3)); % [1 x m]
    
        % ends of line segments
        tau_lb = max(LBtmp(uu) - xtmp(uu), xtmp(vv) - UBtmp(vv)); % a[u] = 1, [1 x m]
        tau_ub = min(UBtmp(uu) - xtmp(uu), xtmp(vv) - LBtmp(vv)); % a[v] = 1, [1 x m]
        tau = min(max(tau, tau_lb), tau_ub);
    
        xtmp(uu) = xtmp(uu) + tau;
        xtmp(vv) = xtmp(vv) - tau;
        x(xInd,:) = xtmp;
    else
        stat.exitflag = 4;
    end


    if stat.niter >= opt.MaxIter % stop when any element reaches MaxIter
        stat.exitflag = 4;
    end
    
    
end




   
 

