
% Firstly please download libsvm toolbox from:
%       http://www.csie.ntu.edu.tw/~cjlin/libsvm/
% Put the downloaded toolbox in the current work directory, and run the
% 'make.m' file contained in the .\libsvm-3.11\matlab folder to compile the
% toolbox. Note here that the SVM algorithm contained in this toolbox is
% written in C for efficiency, but the author has also designed the
% interface for Matlab to call the functions.

% You can also download a local copy from:
%       http://www.viplab.cs.nott.ac.uk/course/G53MLE/libsvm-3.11.zip
% If you download the local copy, then you can use the following command to
% unzip it:
%                unzip('libsvm-3.11.zip');
% The C function in the local copy has already been compiled, hence you can directly use it.

function lab5()

    svm_classification();
    pause;
    randomforest_classification();

end

function svm_classification()
    addpath('libsvm-3.11/matlab');
    load 'MNIST_subset.mat';

    %% Normalization
    for i=1:size(tr_feats,1)
        tr_feats(i,:)=tr_feats(i,:)/sqrt(sum(tr_feats(i,:).^2));
    end
    for i=1:size(te_feats,1)
        te_feats(i,:)=te_feats(i,:)/sqrt(sum(te_feats(i,:).^2));
    end

    %% trainning: train 10 1vs.all binary classifiers. You can also train 10
    %% perceptrons in the same manner.
    classid = unique(sort(tr_label));
    ntrain  = size(tr_feats,1);
    ndim    = size(tr_feats,2);
    for i = 1:length(classid)
        npos = sum(tr_label==classid(i));
        nneg = ntrain - npos;
        wpos = nneg/npos;
        class_labels = (tr_label==classid(i))+1;
        svmstr = sprintf('-t 1 -d 5 -r 1 -g 1 -c 10 -w2 %.2f -b 1',nneg/npos);  % you can change the parameters here to implement a kernel svm. For a detailed description of these parameters, please refer to 'svmtrain.c'
        model{i}.svm = svmtrain(class_labels,tr_feats,svmstr);
    end

    %% testing
    classid = unique(te_label);
    probs = zeros(length(te_label),length(classid));
    for i = 1:length(classid)
        fprintf('Classifying for Digit :  %i \n', classid(i));
        class_label = (te_label==classid(i)) + 1;
        [l,a,d] = svmpredict(class_label,te_feats,model{i}.svm,'-b 1');
        probs(:,i) = d(:,find(model{i}.svm.Label==2));
    end
    [dummy,predicted_label] = max(probs,[],2);
    num_correct = sum(predicted_label == te_label+1);
    acc = num_correct*100/length(te_label);
    fprintf('\n\nAccuracy on MNIST %.2f%% (%.2f%%) \n', acc, 100-acc);
end


function randomforest_classification()
    %% Random forest classifier
    clear all;
    load 'MNIST_subset.mat';
    for i=1:size(tr_feats,1)
        tr_feats(i,:)=tr_feats(i,:)/sqrt(sum(tr_feats(i,:).^2));
    end
    for i=1:size(te_feats,1)
        te_feats(i,:)=te_feats(i,:)/sqrt(sum(te_feats(i,:).^2));
    end

    if ~exist('tmp_data','dir')
        mkdir('tmp_data');
    end
    
    label=zeros(size(tr_feats,1),10);
    for i=1:size(tr_feats,1)
        label(i,tr_label(i)+1)=1;
    end
    
    %% training
    rand('state',0);
    num_of_trees = 100;
    for i=1:num_of_trees     % generate 100 trees
        FH_generate_supervised_tree(tr_feats',i,label);
        fprintf('Finished %d / %d \n',i,num_of_trees);
    end
    
    %% testing
    for i=1:num_of_trees
        fprintf('Parsing tree %d / %d \n',i,num_of_trees);
        model=load(['tmp_data/tree' int2str(i) '.mat']);
        neighbors{i}=FH_parse_tree(te_feats',model);
    end

    catenated_neighbors=neighbors{1};
    for i=1:length(catenated_neighbors)
        for j=2:length(neighbors)
            catenated_neighbors{i}=[catenated_neighbors{i} neighbors{j}{i}];
        end
    end

    neighbors=catenated_neighbors;

    whole_prediction_probability=zeros(10,10000);
    for i=1:length(neighbors)    
        prediction_probability=zeros(10,1);
        current_neighbor=cell2mat(neighbors(i));
        for j=1:length(current_neighbor)
            prediction_probability(tr_label(current_neighbor(j))+1)=prediction_probability(tr_label(current_neighbor(j))+1)+1;
        end
        whole_prediction_probability(:,i)=prediction_probability;
    end
    [dummy,predicted_label] = max(whole_prediction_probability);
    num_correct = sum(predicted_label' == te_label+1);
    acc = num_correct*100/length(te_label);
    fprintf('\n\nAccuracy on MNIST %.2f%% (%.2f%%) \n', acc, 100-acc);    
end















function FH_generate_supervised_tree(feature,type,label)

    global Node Max_depth Min_leaf_sample

    Max_depth=10;   % adjustable parameters
    Min_leaf_sample=4;

    Node=struct('isleaf',[],'depth',[],'global_index',[],'thresh',[],'dim',[]);

    depth=1;
    node_index=1;

    if nargin<5
        global_index=1:size(feature,2);
    end

    if ~exist(['tmp_data/tree' int2str(type) '.mat'])
        grow(feature,depth,node_index,global_index,label);
        save(['tmp_data/tree' int2str(type) '.mat'],'Node');
    end
end



function grow(feature,depth,node_index,global_index,label)
    global Node Max_depth Min_leaf_sample
    Node(node_index)=struct('isleaf',[],'depth',[],'global_index',[],'thresh',[],'dim',[]);

    rand_index=randperm(size(feature,1));
    score=zeros(size(feature,1),1);
    threshold=zeros(size(feature,1),1);
    for i=rand_index(1:50)
        [score(i),threshold(i)] = compute_entropy_and_thresh(feature(i,:),label);
    end
    [dummy,index]=max(score);
    threshold=threshold(index);
    
    Node(node_index).thresh=threshold;
    Node(node_index).dim = index;
    Node(node_index).depth=depth;
    Node(node_index).global_index=global_index;

    left_index=find(feature(index,:)-threshold<=eps);        % to avoid the strange problem caused by precision
    right_index=setdiff(1:size(feature,2),left_index);
    
    if length(left_index)>Min_leaf_sample && depth<Max_depth
        grow(feature(:,left_index),depth+1,node_index*2,global_index(left_index),label(left_index,:));
    else
        makeleaf(depth+1,node_index*2,global_index(left_index));
    end    
    if length(right_index)>Min_leaf_sample && depth<Max_depth
        grow(feature(:,right_index),depth+1,node_index*2+1,global_index(right_index),label(right_index,:));
    else
        makeleaf(depth+1,node_index*2+1,global_index(right_index));
    end
end

function makeleaf(depth,node_index,global_index)
    global Node;
    Node(node_index).depth=depth;
    Node(node_index).isleaf='true';
    Node(node_index).global_index=global_index;
end

function [returned_score,returned_threshold] = compute_entropy_and_thresh(feature,label)  % traditional kd-tree criteria, the weak node test is to pick an arbitrary dimension and use a threshold function

    assert(size(feature,1)==1);
    
    returned_score=-100;
    returned_threshold=9999;
    
    if min(feature)==max(feature)    % impossible to distinguish samples from this feature
        return;
    end
    threshold_vector=min(feature)+(max(feature)-min(feature))*rand(1,15);
        
    for threshold=threshold_vector
        p_left=sum(label(feature<=threshold,:),1);
        p_right=sum(label(feature>threshold,:),1);
        if isempty(p_left) || isempty(p_right) || sum(p_left)==0 || sum(p_right)==0
            score=0;break;
        end
        p_left=p_left/sum(p_left);
        p_right=p_right/sum(p_right);

        entropy_left=zeros(size(p_left));
        index=find(p_left~=0);
        entropy_left(index)=p_left(index).*log2(p_left(index));
        entropy_left=-sum(entropy_left);

        entropy_right=zeros(size(p_right));
        index=p_right~=0;
        entropy_right(index)=p_right(index).*log2(p_right(index));
        entropy_right=-sum(entropy_right);

        ratio=size(find(feature<=threshold),2)/size(feature,2);
        if ratio==0 || ratio==1
            score=0;break;
        end
        score=entropy_left*ratio+entropy_right*(1-ratio);
        score=1/score;
        
        if score>returned_score
            returned_score=score;
            returned_threshold=threshold;
        end
    end
end


function neighbors=FH_parse_tree(feature,model)

    Node=model.Node;
    for i=1:size(feature,2) 
        current_feature=feature(:,i);
        index=1;
        while ~strcmp(Node(index).isleaf,'true')
            if current_feature(Node(index).dim)<=Node(index).thresh
                index=index*2;
            else
                index=index*2+1;
            end
        end
        neighbors{i}=Node(index).global_index;    
    end
end


