clear;
clc;
%读取数据
%Read data
trainImages = loadMNISTImages('train-images.idx3-ubyte');
trainLabels = loadMNISTLabels('train-labels.idx1-ubyte');
testImages = loadMNISTImages('t10k-images.idx3-ubyte');
testLabels = loadMNISTLabels('t10k-labels.idx1-ubyte');
%PCA降维
%PCA dimensionality reduction
square = double(trainImages * trainImages');
[eigVect, eigVal] = eig(square);
eigVal = diag(eigVal);  
eigVal = cumsum(eigVal) / sum(eigVal);
k = find(eigVal >= 1e-2, 1);
pca = eigVect(:, k:end);
trainImages = pca' * trainImages;
testImages = pca' * testImages;
%选取部分样本便于测试
%Select some samples for testing
trainImages = trainImages(:,1:500)';
trainLabels = trainLabels(1:500,1);
testImages = testImages';
%获得样本规模
%Obtain sample size
[m,n] = size(trainImages);
[r,s] = size(testImages);
%调节sig1和C
%Adjust sig1 and C
%当前最佳sig1.sig2,C为12.3,16.75,30
%Current best sig1 Sig2, C = 12.3,16.75,30
tic;
kernels = zeros(m,r);
for i = 1 : r
    kernels(:,i) = Kernel(trainImages,testImages(i,:),12.3);
end
     
% omegas = zeros(n,10);
b = zeros(10,1);
alphas = zeros(m,10);
C = 30;
dis = 0.001;
tolerance = 10;
vote = zeros(r,10);
testResults = zeros(r,10);
Labels = zeros(r,1);
error = 0;
%OvR
fprintf('计算alpha...\n');
for i = 1 : 10
    one = i - 1;
    [b(i),alphas(:,i)] = kernelSVM(trainImages,trainLabels,C,dis,tolerance,one);
    fprintf('第%d个alpha计算完毕\n',i);
%     for j = 1 : 10
%         omegas(:,j) = ((alphas(:,i) .* trainLabels)' * trainImages)';
%     end
end
fprintf('投票中...\n');
%根据alpha和b计算出结果y
for i = 1 : 10
    tmpLabels = trainLabels;
    y = zeros(r,10);
    for j = 1 : m
        if tmpLabels(j) == (i - 1)
            tmpLabels(j) = 1;
        else
            tmpLabels(j) = -1;
        end
    end
    for j = 1 : r
        y(j,i) = (alphas(:,i) .* tmpLabels)' * kernels(:,j) + b(i);
%         y = omegas(:,i)' * testImages(j,:)' + b(i);
%        fprintf('y的值为%d\n',y);
        vote(j,i) = sign(y(j,i));
    end
    clear tmpLabels;
end
%投票得到待测样本类别
for i = 1 : r
    for j = 1 : 10
        if vote(i,j) == 1
            testResults(i,j) = testResults(i,j) + 1;
        else
            testResults(i,:) = testResults(i,:) + 1;
            testResults(i,j) = testResults(i,j) - 1;
        end
    end
end
fprintf('计算结果...');
for i = 1 : r
    [maxCount,pos] = max(testResults(i,:));
    Labels(i) = pos - 1;
end

for i = 1 : r
    if Labels(i) ~= testLabels(i)
        error = error + 1;
%         fprintf('错误%d',Labels(i));
%         fprintf('正确%d\n',testLabels(i));
    end
end

fprintf('正确率为%d\n',1-error/r);
toc;

% save data;
% system('shutdown -s');