function [alphaStar, bStar, SVIndex] = yxcSVMtrain(X, Y, C, kernel, sigma)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If you have the Chinese book: 《數據挖掘中的新方法——支持向量機》, find the algorithm here:
% Classify C-SVC using algorithm 5.4.12 on page 193
%
% Otherwise, refer to this article: "Support Vector Machines: Hype or Hallelujah?" 
% Section 6: SUMMARY OF SVM METHOD on page 5
% http://www.sigkdd.org/explorations/issue2-2/bennett.pdf
%
% Here is also a very good tutorial: "A Tutorial on Support Vector Machines for Pattern Recognition"
% www.umiacs.umd.edu/~joseph/support-vector-machines4.pdf
% 
% A detailed description about iris flower data set:
% http://en.wikipedia.org/wiki/Iris_flower_data_set
%
% Author: Xuchen Yao. yaoxuchen@gmail.com
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Input:
% X: num x dim Data, X has num points, every point is described by dim vectors
% Y: n x 1 Data, {1, -1}
% C: penalty
% kernel: see yxcSVMkernel
% sigma: see yxcSVMkernel
% Output:
% alphaStar: Optimized alpha, see page 193
% bStar: bias, see page 193
% SVIndex: Index for support vectors

[num, dim] = size(X);
if dim ~= 2
    return;
end
if num ~= length(Y)
    return;
end
Y = Y(:);

%以下為純粹按照《數據挖掘中的新方法-支持向量機》P193頁算法，將其轉化為可用二次規劃工具包進行
%求解的形式，具體轉化方法參見P181證明中所設
H = (Y*Y').*yxcSVMkernel(X, X, kernel, sigma);
f = -ones(num, 1);
A = zeros(1, num);
b = 0;
Aeq = Y';
beq = 0;
lb = zeros(num, 1);
ub = C .* ones(num, 1);
x0 = zeros(num,1);
%qp_options = optimset('Display','off');
qp_options = optimset('MaxIter',10^3, 'LargeScale', 'off', 'Display','off');
[alphaStar,fval,exitflag,output] = quadprog(H,f,[],[],Aeq,beq,lb,ub,x0,qp_options);

nearZero = 10^-12;
% Assume the minor ones are all zero.
alphaStar(find(abs(alphaStar) < nearZero)) = 0;
alphaStar(find(alphaStar > C - nearZero)) = C;
% support vectors are those whose alpha value > 0
SVIndex = find(alphaStar > 0);
% support vectors on bound are those whose alpha value == max(alphaStar) == C
SVNotOnBoundIndex = find(alphaStar > 0 & alphaStar < max(alphaStar));
 
% bStar is an average value, not a random one. 
if ~isempty(SVNotOnBoundIndex)
    bStar = sum( Y(SVNotOnBoundIndex) -  H(SVNotOnBoundIndex, SVIndex) * alphaStar(SVIndex) .* Y(SVNotOnBoundIndex) )...
                 / length(SVNotOnBoundIndex);
else
    bStar = 0;
end