function ypred = algoHM(Xtrain, ytrain, Xtest, C, B)

if(~exist('C','var'))
    C = 0.183; % 0.525 ± 0.006
end

if(~exist('B','var'))
    B = 0.0884; % 0.525 ± 0.006
end

gamma = 5000;

params = sprintf('-s 1 -c %f -B %f -q',C,B);

classes = size(ytrain,2);

Ztrain = bsxfun(@rdivide,Xtrain,sqrt(sum(Xtrain.*Xtrain,2)));
Ztest = bsxfun(@rdivide,Xtest,sqrt(sum(Xtest.*Xtest,2)));

clear Xtrain Xtest;

S0 = cell(classes,1);
S1 = cell(classes,1);

firstClass0 = zeros(1,classes);

% --- Train ---

ytrain = ytrain*2-1;
weights = zeros(size(Ztrain,1),classes);
for i=1:classes
    S0{i} = train(full(ytrain(:,i)), Ztrain, params);
    firstClass0(i) = S0{i}.Label(1);
    [~,~,weights(:,i)] = predict(zeros(size(Ztrain,1),1), Ztrain, S0{i});
end
weights=bsxfun(@times,weights,firstClass0);

ypred = +(weights>0);
iZero = 1:size(ypred,1);
iZero = iZero(~any(ypred,2));
for i=iZero
    [~,j] = max(weights(i,:));
    ypred(i,j) = 1;
end

Etrain = [Ztrain gamma*ypred];
Etrain = bsxfun(@rdivide,Etrain,sqrt(sum(Etrain.*Etrain,2)));
for i=1:classes
    S1{i} = train(full(ytrain(:,i)), Etrain, params);
end

clear C params ypred ydec *train i j idx;

% --- Test ---

weights = zeros(size(Ztest,1),classes);
for i=1:classes
    [~,~,weights(:,i)] = predict(zeros(size(Ztest,1),1), Ztest, S0{i});
end
weights=bsxfun(@times,weights,firstClass0);

ypred = +(weights>0);
iZero = 1:size(ypred,1);
iZero = iZero(~any(ypred,2));
for i=iZero
    [~,j] = max(weights(i,:));
    ypred(i,j) = 1;
end

Etest = [Ztest gamma*ypred];
Etest = bsxfun(@rdivide,Etest,sqrt(sum(Etest.*Etest,2)));
ypred = zeros(size(Etest,1),classes);
for i=1:classes
    ypred(:,i) = predict(zeros(size(Etest,1),1), Etest, S1{i});
end
ypred = +(ypred>0);

%fprintf('%.4f %.4f ',C,B);

end

