cluster_gen.m

clear;
addpath('../ManifoldClassification');

% parameters
n = 5000; % number of points
k = 5; % number of nearest neighbors
dim = 3;

% create data
ppm = [0.5*n 0.5*n];
centers = [5 12; 12 12;];
stdev = 2;
[data2,labels]=makegaussmixnd(centers,stdev,ppm);
labels = labels -1;

% convert to manifold
t = data2(:,1);
u = data2(:,2);
data3 = [t.*cos(t) u t.*sin(t) labels'];
data3 = data3(randperm(size(data3,1)),:);
labels = data3(:,end);
coords = data3(:,1:end-1);
% clean up workspace
clearvars data2 t u stdev;

% k nearest neighbors
d_euc = squareform( pdist(coords, 'euclidean' ) ); 
[~, nn_idx] = sort(d_euc);
nn_idx = nn_idx(2:k+1,:); % indicies
% clean up workspace
clearvars d_euc;


% adjacency matrix & diversity scores
adj = zeros(n,n); % adjacency matrix
pre_div = zeros(k,n); % preliminary diversity scores
div_scores = zeros(n,n);
for c = 1:n
    col = nn_idx(:,c);
    for cc = 1:k
        idx = col(cc);
        adj(idx,c) = 1;
        pre_div(cc,c) = ~(labels(idx) == labels(c));
    end
    div_scores(c,:) = adj(c,:)*(sum(pre_div(:,c))/k);
end
% based on k nearest neighbors and their k nearest neighbours
div = zeros(n,1); % diversity scores [0,1] where 1 is most heterogenous
for c = 1:n
    div(c) = sum(div_scores(:,c))/(k);
end
% clean up workspace
clearvars pre_div div_scores idx col c cc;

div_sort = sortrows(div);
max_div = div_sort(end); % maximum diversity score
st = std(div(div>0)); % standard deviation of diversity scores
% range of scores based on standard deviation
range = [max_div-(2*st) max_div]; 

% % Plot diversity scores & threshold
% plot(1:n, div_sort, '-')
% title('Diversity Scores')
% hold on
% plot([1 n], [range(1) range(1)], 'r-')
% hold off;
% % clean up workspace
% clearvars max_div div_sort st;

% find point cloud
method = zeros(n,1);
idx_range = find(div > range(1));
method(idx_range) = 1;

% find convex hull
cloud_nn = unique(reshape(nn_idx(:,idx_range)',[],1));
cloud = coords(cloud_nn,:);
ch = delaunayn(cloud); % delaunay triangulation
% find points in hull
in_hull = ~isnan(tsearchn(cloud, ch, coords));
cloud_idx = find(in_hull);

method(cloud_idx) = 1;
% % plot manifold with clusters
% plot_labelled(coords, 'a', labels, 'gb', 'Manifold');
% % plot method separation
% plot_labelled(coords, 'a', method, 'kr', 'nn Point Cloud');

full_data = [coords labels method];
reg_idx = full_data(:,5) == 1;
knn_idx = full_data(:,5) == 0;
reg_coords = coords(reg_idx, :);
reg_labels = labels(reg_idx);
knn_coords = coords(knn_idx, :);
knn_labels = labels(knn_idx);

clearvars reg_idx knn_idx full_data method idx_range div max_div;

% create models
knn_model = fitcknn(knn_coords, knn_labels, 'NumNeighbors', 5, ...
    'Distance', 'euclidean');
reg_model = glmfit(reg_coords, reg_labels, 'binomial');


% k-fold evaluation
parts = 10;
[kn, ~] = size(knn_coords);
rn = n - kn;
kt = floor(kn/parts);
rt = floor(rn/parts);
kfold_mse = zeros(2, 10);
for c=1:parts
    kr = [(c-1)*kt c*kt]; % range of indicies for knn test data
    rr = [(c-1)*rt c*rt]; % range of indicies for regression test data
    if c==1 % sets training interval to all points after test data
        knn_model = fitcknn(knn_coords((kr(2)+1):kn,:), ...
            knn_labels((kr(2)+1):kn,:), 'NumNeighbors', 5, ...
            'Distance', 'euclidean');
        reg_model = glmfit(reg_coords((rr(2)+1):rn,:), ...
            reg_labels((rr(2)+1):rn,:), 'binomial');
    elseif c==parts % sets training interval to all points before test data
        knn_model = fitcknn(knn_coords(1:kr(1),:), knn_labels(1:kr(1),:), ...
            'NumNeighbors', 5, 'Distance', 'euclidean');
        reg_model = glmfit(reg_coords(1:rr(1),:), reg_labels(1:rr(1)), ...
            'binomial');
        % compensating for rounding
        kr(2) = kn;
        rr(2) = rn;
    else % sets training interval to points surrounding test interval
        knn_model = fitcknn(knn_coords([1:kr(1) (kr(2)+1):kn],:), ...
            knn_labels([1:kr(1) (kr(2)+1):kn],:), ...
            'NumNeighbors', 5, 'Distance', 'euclidean');
        reg_model = glmfit(reg_coords([1:rr(1) (rr(2)+1):rn],:), ...
            reg_labels([1:rr(1) (rr(2)+1):rn]), 'binomial');
    end
    % knn evaluation
    knn_pred = predict(knn_model, knn_coords((kr(1)+1):kr(2), :));
    knn_comp = [knn_pred knn_labels((kr(1)+1):kr(2))];
    kfold_mse(1,c) = sum((knn_comp(:,1)-knn_comp(:,2)).^2)/kn;
    % regression evaluation
    reg_pred = glmval(reg_model, reg_coords((rr(1)+1):rr(2),:), 'logit');
    reg_comp = [reg_pred reg_labels((rr(1)+1):rr(2))];
    kfold_mse(2,c) = sum((reg_comp(:,1)-reg_comp(:,2)).^2)/rn;
end
% averaged MSE across all evaluations
avg_mse = [sum(kfold_mse(1,:))/parts; sum(kfold_mse(2,:))/parts];
model_mse = sum(avg_mse)/2;

clearvars kr rr kt rt kn rn parts;

% near close error
% quantifying uncertainty ^ ?
% F1 graph with inflection point