clear; % clear variables from memory

% read in the dataset
nsample = 100;
fid = fopen('1d.noise', 'r'); % open dataset for read
%fid = fopen('1d.clean', 'r');
[Dataset, count] = fscanf(fid, '%f %f', [2, nsample]); % read in examples (x,y) as a 2 * nsample matrix
fclose(fid); % close file id
Dataset = Dataset'; % matrix transpose, one of the matrix operations in MatLab
X = Dataset(:,1); % first column of matrix Dataset is assigned to X
Y = Dataset(:,2); % second column of matrix Dataset is assigned to Y

% the trade-off weights we want to investigate
C = [0.01, 0.1, 1, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, inf];
Margin = []; % margin; initialized as null
nSV = []; % number of support vector;
nMis = []; % number of misclassification;
Err = []; % training errors;

for n = 1 : max(size(C)),

% construct Hessian matrix; Hessian matrix is the Q matrix in our slides; also called Kernel matrix
H = zeros(nsample, nsample); % initialize H; set H to a nsample * nsample zero matrix
for i = 1 : nsample,
for j = 1 : nsample,
H(i,j) = ; % !!! please write your answer here !!!
end
end
H = H+1e-10*eye(size(H)); % add 1e-10 to the main diagonal of H; a trick to make H stable

F = -ones(nsample,1); % F' * Alpha corresponds to sigma_i(Alpha_i) in object function

% set up equality constraints
A = Y'; % corresponds to sigma_i(Alpha_i * Y_i) = 0
b = 0;

% set up upper and lower bounds for alpha: LB <= Alpha <= UB
UB = ; % !!! please write your answer here !!! note that UB and LB are nsample * 1 matrixes
LB = ; % !!! please write your answer here !!!

% starting point of alpha
Alpha0 = zeros(nsample, 1);

% optimizing alpha with quadratic programming
Alpha = quadprog(H, F, [], [], A, b, LB, UB, Alpha0);
% Alpha = qp(H, F, A, b, LB, UB, Alpha0, 1);

% tolerance for support vector detection; we will ignore the alphas less than tol
tol = 0.0001;

% calculate weight
w = 0;
for i = 1 : nsample,
w = w + Alpha(i) * Y(i) * X(i);
end

% calculate bias
bias = 0;
b1 = 0;
b2 = 0;
for i = 1 : nsample,
if (Alpha(i) > tol & Alpha(i) < C(n) - tol),
b1 = b1 + X(i) * w - Y(i);
b2 = b2 - 1;
end
end

if b2 ~= 0,
bias = b1 / b2;
else % unlikely
b1 = 0;
for i = 1 : nsample,
if Alpha(i) < tol,
b1 = b1 + X(i) * w - Y(i);
b2 = b2 - 1;
end
end

if b2 ~= 0,
bias = b1 / b2;
else % even unlikelier
b1 = 0;
for i = 1 : nsample,
b1 = b1 + X(i) * w - Y(i);
b2 = b2 - 1;
end
if b2 ~= 0,
bias = b1 / b2;
end
end
end

% margin = 2 / ||w||
Margin = [Margin, 2 / abs(w)]; % the operation A = [A, v] appends v to matrix A

% number of support vectors
nSV = [nSV, size(find(Alpha > tol), 1)];

% calculate # of misclassification and training error
m = 0;
e = 0;
for i = 1 : nsample,
predict = w * X(i) + bias; % Y = w * X + b
if predict >= 0 & Y(i) < 0,
m = m + 1;
end
if predict < 0 & Y(i) >= 0,
m = m + 1;
end
if Alpha(i) > tol, % consider support vectors only; why?
e = e + 1 - predict * Y(i);
end
end
nMis = [nMis, m];
Err = [Err, e];
end

% plot C_margin, C_trainingerror, C_misclassification, C_nsupportvector
% please use your code to make better plots instead of ours

Z = zeros(size(C));
for i = 1 : size(C, 2)
Z(i) = i;
end

figure
plot(Z, Margin);
title('Margin');
xlabel('C(i)');

figure
plot(Z, Err);
title('Training Error');
xlabel('C(i)');

figure
plot(Z, nMis);
title('# of Misclassification');
xlabel('C(i)');

figure
plot(Z, nSV);
title('# of Support Vector');
xlabel('C(i)');