function [signed_alphas,b,f,w] = svmtrain_cplex(x_trn,y_trn,K,parms)
%
% Solve an SVM using CPLEX (Train + test)
% Inputs: 
%    parms: 
%      C       - miscalssification penalty 
%      K       - precalculated kernel; (use kertyp='givn');
% 
% Outputs: 
%      w - separating hyper plane 
%      alphas - largange multipliers (support vectors coeffs)
%      f - the value of the target fucntion at the optimum
% 
% (C) gal Chechik 2005
  
% Check args 
  y_trn = y_trn(:);
  n_trn = size(x_trn,1);
  if(n_trn ~= size(y_trn,1)), error('Inconsistent sizes'); end;
  if(any(isnan(x_trn(:)))), 
    error('\n\n==> Input contains NaN values\n\n\n');
  end
  if size(K,1) ~= n_trn
    error('\n\n==> Inconsistent kernel matrix %d %d n_trn=%d\n\n\n',size(K),n_trn);
  end
  
  % Train SVM
  [w,alphas,b,f]=train_dual(x_trn,y_trn,K,parms);
  signed_alphas = alphas .* y_trn;
  
  return
end
% ======================================  
function [wopt,alphas,b,fopt]=train_dual(x,y,K,parms)
%
% Interface the dual to cplex 

% Init
  [parms,c]         = take_from_struct(parms,'C',[]);
  [parms,cplex_source]  = take_from_struct(parms,'cplex_source','ethz');
  if isempty(K), error('empty K'); end;
  
  % The variables vector is the alphas (a_1,...,a_n)
  [n,d] = size(x);  
  if any(size(K) ~= [n n]), error('Inconsistent dimension'); end
    
  switch (cplex_source)
   case 'ethz', [alphas,fopt] = call_qpcplex_ethz(x,y,K,parms);
   case 'dmus', [alphas,fopt] = call_qpcplex_dmus(x,y,K,parms);    
  end
  alphas_orig = alphas;
  maxlogalpha = max(log10(alphas_orig));
  
  svs = find(log10(alphas)  > maxlogalpha-6);
  nosvs = (log10(alphas) <= maxlogalpha-6);
  alphas(nosvs)=0;
  wopt = zeros(1,d);
  for i=1:n
    wopt = wopt + alphas(i) .*y(i) .* x(i,:);
  end

  % Compute the threshold. Eq (7.32) p.203 
  bj = zeros(1,length(svs));
  for jj=1:length(svs)
    j = svs(jj); s=0;
    for i=1:n
      s = s + y(i)*alphas(i)*K(j,i);
    end
    bj(jj) = y(j) - s; 
  end
  if(isempty(bj)), 
    b = 0 ; 
  else
    b = mean(bj);
  end
  
  
return
end

% =====================================================
function [alphas,fopt] = call_qpcplex_ethz(x,y,K,parms)
%
% QP_CPLEX ETHZ 
%     [min|max] .5 x'Qx+C'x
%     s.t.  Ax ["="|"<="|">="] b
%           {x <= UB}
%           {x >= LB}
  
  [parms,c]         = take_from_struct(parms,'C', []);
  [parms,verbose]   = take_from_struct(parms,'verbose',1);  
  [parms,save_mps]  = take_from_struct(parms,'save_mps',0);
  n = size(x,1);
  
  opt_type = -1; % maximize

  % Target function \sum_{i=1:d} w_i^2 + c \sum_{i=1:n} xi_i
  Q = zeros(n,n);
  for i=1:n
    for j=1:n
      Q(i,j) = y(i) * y(j) * K(i,j);
    end
  end
  C = ones(n,1);

  % Correct assymetry due to numerical errors
  Q = (Q + Q')/2;
  
  % sum (alphai yi)=0
  A = y';
  b = 0;
  ctype = 'E';
  rangevar = 0;
  
  % 0<alpha_i <c
  lb = ones(n,1);
  lb(1:n) = 0;  
  ub = ones(n,1);
  ub(1:n) = c;    

  try
    [alphas,fopt] = qp_cplex(opt_type,-Q,C,A,b, ctype,lb,ub,rangevar, ...
			     verbose, save_mps);
  catch
    % Fix PSD's due to numerical errors 
    [v,d] = eig(Q);dd=diag(d);
    if(min(dd)<0),
      fprintf('fix negative eigvals\n');
      Q2 = Q - eye(size(Q))*min(dd);
    end

    try 
      [alphas,fopt] = qp_cplex(opt_type,-Q2,C,A,b, ctype,lb,ub,rangevar, ...
			       verbose, save_mps);
    catch
      % Fix PSD's due to numerical errors 
      [v,d] = eig(Q2);dd=diag(d);
      if(min(dd)<0),
	fprintf('fix negative eigvals AGAIN!\n');
	Q3 = Q2 - eye(size(Q2))*min(dd);
      end

    [alphas,fopt] = qp_cplex(opt_type,-Q3,C,A,b, ctype,lb,ub,rangevar, ...
			     verbose, save_mps);
    end % Second try/catch 

  end % First try/catch 


  if(max(alphas)<=0) 
    fprintf('\n\nOptimization failed: all alphas are zeros\n\n');
    % input('press enter to continue');
  end;

return
end


% =====================================================
function [alphas,fopt] = call_qpcplex_dmus(x,y,K,parms)
%
%   QP(H,f,A,b) solves the quadratic programming problem:
% 
%         min 0.5*x'Hx + f'x   subject to:  Ax <= b
%          x
%
% Note: we maximize for alpha, but the dual 
% is       sum ai  -sum_{ij} ai yi K() yj aj
% so we give +H to the optimizer 
% 
% Note: This funciton uses the interface to CPLEX written by 
% David R. Musicant

  [parms,c]         = take_from_struct(parms,'C', []);
  [parms,save_mps]  = take_from_struct(parms,'save_mps',0);%#ok
  n = size(x,1);

  % Target function \sum_{i=1:d} w_i^2 + c \sum_{i=1:n} xi_i
  H = zeros(n,n);
  for i=1:n
    for j=1:n
      H(i,j) = y(i) * y(j) * K(i,j);
    end
  end
  f = ones(n,1);
  
  % sum (alphai yi)=0
  A = y';
  b = 0;
  
  % 0<alpha_i <c
  lb = ones(n,1);
  lb(1:n) = 0;  
  ub = ones(n,1);
  ub(1:n) = c;    

  % Note: we maximize over alphas, but the dual is       
  %           sum ai  -sum_{ij} ai yi K(xi,xj) yj aj
  % so we give +H to the optimizer   
  n_equality=1;
  X0=ones(1,n)/n;
  [fopt,alphas] = qp_cplex_dmus(H,f,A,b,lb,ub,X0,n_equality);

  return
end