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)');