%{
A quick multivariate decision tree that I whipped up.
Multivariate here implies splitting across many attributes
via a hyperplane instead of a single attribute. This demo only
deals with 2D data.

The splitting is done using a RANSAC approach by trying different
pairing of data points to find the best split. 
Not exactly optimal but easy to implement!
%}

% parameters
iters = 1000; % ransac iterations
tree_depth = 4;
min_pts_per_node = 5; % minimum points for positive/negative samples to be split

% synthetic data generation
pts = 1000;
radius = 10;
inner_radius = 5;

% generate the data
label = [];
X = [];

for i=1:pts
    x = -radius + 2*radius*rand();
    y = -radius + 2*radius*rand();
    r = x*x + y*y;

    X = [X; x y];

    angle = atan2(y,x)*180/pi;

    if r < inner_radius*inner_radius % circle data
    %if (angle > -135 && angle < -45) % triangle data
        label = [label; 1];
    else
        label = [label; 0];
    end
end

figure(1);
clf
hold on

% plot the samples
pos = X(find(label == 1), :);
neg = X(find(label == 0), :);
plot(pos(:,1), pos(:,2), "@12", neg(:,1), neg(:,2), "@13");
xlabel("x");
ylabel("y");

% data structure that gets passed around
root.X = X;
root.label = label;
root.index = 1;

to_split = [root];

% allocate memory for the decision tree parameters
% using a matrix to store the parameters
tree = zeros(2^tree_depth, 5); % column 1-3 is theta parameter for the splitting line, column 4,5 is the jump index for left/right node
next_free = 2; % next free index

% splitting process ...
% done in a Breadth first search manner
for s=1:tree_depth
    tmp_split = [];

    for i=1:size(to_split,1)
        X2 = to_split(i).X;
        label2 = to_split(i).label;

        pos_idx = find(label2 == 1);
        neg_idx = find(label2 == 0);

        if size(pos_idx,1) == 0 || size(neg_idx, 1) == 0
            continue;
        end

        pos = X2(pos_idx, :);
        neg = X2(neg_idx, :);

        best_gain = 0;
        best_theta = [];
        best_pos_side = [];
        best_neg_side = [];

        for j=1:iters
            idx1 = int32(size(pos,1)*rand()) + 1;
            idx2 = int32(size(neg,1)*rand()) + 1;

            idx1 = min(size(pos,1), idx1);
            idx2 = min(size(neg,1), idx2);

            p = pos(idx1,:);
            n = neg(idx2,:);

            [gain, pos_side, neg_side, theta] = SplitData(X2, label2, pos, neg, p, n);

            if gain > best_gain
                best_gain = gain;
                best_theta = theta;
                best_pos_side = pos_side;
                best_neg_side = neg_side;
            end
        end

        if s != tree_depth
            % check left node's purity
            pos_idx = find(label2(best_pos_side) == 1);
            neg_idx = find(label2(best_pos_side) == 0);

            if size(pos_idx,1) == 0 || size(neg_idx, 1) == 0
                printf("Pure node found in left\n");
            else
                if size(pos_idx,1) > min_pts_per_node && size(neg_idx,1) > min_pts_per_node
                    left.X = X2(best_pos_side, :);
                    left.label = label2(best_pos_side, :);
                    left.index = next_free++;

                    tree(to_split(i).index, 4) = left.index;
                    tmp_split = [tmp_split; left;];   
                end
            end

            % check right node's purity
            pos_idx = find(label2(best_neg_side) == 1);
            neg_idx = find(label2(best_neg_side) == 0);

            if size(pos_idx,1) == 0 || size(neg_idx, 1) == 0
                printf("Pure node found in right\n");
            else
                if size(pos_idx,1) > min_pts_per_node && size(neg_idx,1) > min_pts_per_node
                    right.X = X2(best_neg_side, :);
                    right.label = label2(best_neg_side, :);
                    right.index = next_free++;

                    tree(to_split(i).index, 5) = right.index;
                    tmp_split = [tmp_split; right;];
                end
            end
        end

        tree(to_split(i).index, 1:3) = best_theta;

        col = [0;1;0]; % colour
        PlotLine(best_theta, [-10 -10], [10 10], col, 2);
    end

    to_split = tmp_split;
end

%tree

% Classify and get accuracy
correct = 0;

for i=1:pts
    r = Classify(tree, X(i,:));

    if label(i) == 1 && r >= 0
        correct++;
    elseif label(i) == 0 && r < 0
        correct++;
    endif
end

accuracy = correct / pts;

str = sprintf("Accuracy: %.2f%%\n", accuracy*100);
title(str);

% Plot the boundary lines
spacing = 100;

u = linspace(-radius, radius, spacing);
v = u;
z = zeros(spacing, spacing);

for i=1:length(u)
    for j=1:length(v)
        r = Classify(tree, [u(i) v(j)]);
        z(i,j) = r;
    end
end

z = z'; % need to transpose to get correct visual result
contour(u, v, z, [0, 0], 'r', 'LineWidth', 2)
axis([-10,10,-10,10],"square");
hold off

