% Parameters selection for c and gamma
% pick the kernel type
% pick a set of parameter (c,g)
% n-fold cross validation
% report the across-n-fold average accuracy acc(c,g) as a function of c and g
clear
clc
close all
% addpath to the libsvm toolbox
addpath('../libsvm-3.12/matlab');
%%
% % % % make a dataset
% % % dirData = './data';
% % % Nc = 4;
% % % Ns = 50;
% % % h = 15;
% % % r = 3;
% % % [data, label, run] = generateSpiralDataWithLabels(Nc,Ns,h,r);
% % % save(fullfile(dirData,'spiral_Nc4_test'),'data','label','run');
%% read the data set Spiral
dirData = './data';
load(fullfile(dirData,'spiral_Nc10_train'));
rawTrainData = data(:,1:2);
rawTrainLabel = label;
NTrain = size(rawTrainData,1);
[sortedTrainLabel, permIndex] = sortrows(rawTrainLabel);
sortedTrainData = rawTrainData(permIndex,:);
load(fullfile(dirData,'spiral_Nc10_test'));
rawTestData = data(:,1:2);
rawTestLabel = label;
NTest = size(rawTestData,1);
[sortedTestLabel, permIndex] = sortrows(rawTestLabel);
sortedTestData = rawTestData(permIndex,:);
% combine the data together just to fit my format
totalData = [sortedTrainData; sortedTestData];
totalLabel = [sortedTrainLabel; sortedTestLabel];
figure;
subplot(1,2,1); imagesc(totalLabel); title('class label');
subplot(1,2,2); imagesc(totalData); title('features');
[N D] = size(totalData);
labelList = unique(totalLabel(:));
NClass = length(labelList);
% #######################
% Determine the train and test index
% #######################
trainIndex = zeros(N,1); trainIndex(1:NTrain) = 1;
testIndex = zeros(N,1); testIndex( (NTrain+1):N) = 1;
trainData = totalData(trainIndex==1,:);
trainLabel = totalLabel(trainIndex==1,:);
testData = totalData(testIndex==1,:);
testLabel = totalLabel(testIndex==1,:);
%%
% #######################
% Automatic Cross Validation
% Parameter selection using n-fold cross validation
% #######################
stepSize = 5;
bestLog2c = 0;
bestLog2g = log2(1/D);
epsilon = 0.005;
bestcv = 0;
Ncv = 3; % Ncv-fold cross validation cross validation
deltacv = 10^6;
Nlimit = 1000;
cnt = 1;
breakLoop = 0;
while abs(deltacv) > epsilon && cnt < Nlimit
bestcv_prev = bestcv;
prevStepSize = stepSize;
stepSize = prevStepSize/2;
log2c_list = bestLog2c-prevStepSize: stepSize: bestLog2c+prevStepSize;
log2g_list = bestLog2g-prevStepSize: stepSize: bestLog2g+prevStepSize;
numLog2c = length(log2c_list);
numLog2g = length(log2g_list);
cvMatrix = zeros(numLog2c,numLog2g);
for i = 1:numLog2c
log2c = log2c_list(i);
for j = 1:numLog2g
log2g = log2g_list(j);
% % % % With some kernel
% % % cmd = ['-q -c ', num2str(2^log2c), ' -g ', num2str(2^log2g),' -t 2'];
% % % cv = get_cv_ac(trainLabel, [(1:NTrain)' trainData*trainData'], cmd, Ncv);
% With some precal kernel
cmd = ['-q -c ', num2str(2^log2c), ' -g ', num2str(2^log2g)];
cv = get_cv_ac(trainLabel, trainData, cmd, Ncv);
if (cv >= bestcv),
bestcv = cv; bestLog2c = log2c; bestLog2g = log2g;
bestc = 2^bestLog2c; bestg = 2^bestLog2g;
end
disp(['So far, cnt=',num2str(cnt),' the best parameters, yielding Accuracy=',num2str(bestcv*100),'%, are: C=',num2str(bestc),', gamma=',num2str(bestg)]);
% Break out of the loop when the cnt is up to the condition
if cnt > Nlimit, breakLoop = 1; break; end
cnt = cnt + 1;
end
if breakLoop == 1, break; end
end
if breakLoop == 1, break; end
deltacv = bestcv - bestcv_prev;
end
disp(['The best parameters, yielding Accuracy=',num2str(bestcv*100),'%, are: C=',num2str(bestc),', gamma=',num2str(bestg)]);
%%
% % % % #######################
% % % % Semi-manual cross validation
% % % % Parameter selection using n-fold cross validation
% % % % #######################
% % % % === traditional manual cv =====
% % % bestcv = 0;
% % % for log2c = -1:1:5,
% % % for log2g = -8:1:3,
% % % cmd = ['-q -c ', num2str(2^log2c), ' -g ', num2str(2^log2g)];
% % % cv = get_cv_ac(trainLabel, trainData, cmd, 3);
% % % if (cv >= bestcv),
% % % bestcv = cv; bestc = 2^log2c; bestg = 2^log2g;
% % % end
% % % %fprintf('%g %g %g (best c=%g, g=%g, rate=%g)\n', log2c, log2g, cv, bestc, bestg, bestcv);
% % % end
% % % end
% % % disp(['The best parameters, yielding Accuracy=',num2str(bestcv*100),'%, are: C=',num2str(bestc),', gamma=',num2str(bestg)]);
% % % % ===============================
%%
% #######################
% Train the SVM in one-vs-rest (OVR) mode
% #######################
% % % % With specific kernel
% % % bestParam = ['-q -c ', num2str(bestc), ', -g ', num2str(bestg),' -t 2'];
% % % model = ovrtrainBot(trainLabel, [(1:NTrain)' trainData*trainData'], bestParam);
% Without specific kernel
bestParam = ['-q -c ', num2str(bestc), ', -g ', num2str(bestg)];
model = ovrtrainBot(trainLabel, trainData, bestParam);
% bestParam = ['-q -c 8 -g 0.0625'];
% #######################
% Classify samples using OVR model
% #######################
% % With specific kernel
% [predict_label, accuracy, decis_values] = ovrpredictBot(testLabel, [(1:NTest)' testData*trainData'], model);
% With specific kernel
[predict_label, accuracy, decis_values] = ovrpredictBot(testLabel, testData, model);
[decis_value_winner, label_out] = max(decis_values,[],2);
% #######################
% Make confusion matrix
% #######################
[confusionMatrix,order] = confusionmat(testLabel,label_out);
% Note: For confusionMatrix
% column: predicted class label
% row: ground-truth class label
% But we need the conventional confusion matrix which has
% column: actual
% row: predicted
figure; imagesc(confusionMatrix');
xlabel('actual class label');
ylabel('predicted class label');
totalAccuracy = trace(confusionMatrix)/NTest;
disp(['Total accuracy from the SVM: ',num2str(totalAccuracy*100),'%']);
%%
% #######################
% Plot the results
% #######################
figure;
subplot(1,3,2); imagesc(predict_label); title('predicted labels'); xlabel('class k vs rest'); ylabel('observations'); colorbar;
subplot(1,3,1); imagesc(decis_values); title('decision values'); xlabel('class k vs rest'); ylabel('observations'); colorbar;
subplot(1,3,3); imagesc(label_out); title('output labels'); xlabel('class k vs rest'); ylabel('observations'); colorbar;
% plot the true label for the test set
tmp = min(exp(zscore(decis_value_winner)),10);
tmp = tmp-min(tmp(:))+1;
tmp = tmp/max(tmp);
patchSize = 50*tmp;
colorList = generateColorList(NClass);
colorPlot = colorList(testLabel,:);
figure;
scatter(testData(:,1),testData(:,2),patchSize, colorPlot,'filled'); hold on;
% plot the predicted labels for the test set
patchSize = patchSize/2;
colorPlot = colorList(label_out,:);
scatter(testData(:,1),testData(:,2),patchSize, colorPlot,'filled');
%%
% #######################
% Plot the decision boundary
% #######################
% Generate data to cover the domain
minData = min(totalData,[],1);
maxData = max(totalData,[],1);
stepSizePlot = (maxData-minData)/50;
[xI yI] = meshgrid(minData(1):stepSizePlot(1):maxData(1),minData(2):stepSizePlot(2):maxData(2));
% #######################
% Classify samples using OVR model
% #######################
fakeData = [xI(:) yI(:)];
% % % % With specific kernel
% % % [pdl, acc, dcsv] = ovrpredictBot(xI(:)*0,[(1:size(fakeData))' fakeData*trainData'], model);
% Without specific kernel
[pdl, acc, dcsv] = ovrpredictBot(xI(:)*0,fakeData, model);
% Note: when the ground-truth labels of testData are unknown, simply put
% any random number to the testLabel
[dcsv_winner, label_domain] = max(dcsv,[],2);
% plot the result
tmp = min(exp(zscore(dcsv_winner)),10);
tmp = tmp-min(tmp(:))+1;
tmp = tmp/max(tmp);
patchSize = 50*tmp;
colorList = generateColorList(NClass);
colorPlot = colorList(label_domain,:);
figure;
scatter(xI(:),yI(:),patchSize, colorPlot,'filled');
