% Cluster 2D real-world graphs with spectral clustering and compare with k-means
% USI, ICS, Lugano
% Numerical Computing 

clear variables;
close all;
warning OFF;

addpath ../datasets
addpath ../datasets/Meshes

MATS = {'airfoil1', 'barth', 'grid2', '3elt'};

for i = 1:length(MATS)
    load(strcat(MATS{i}, '.mat'));

    % Specify the number of clusters
    K = 4;
    % Read graph
    W   = Problem.A;
    Pts = Problem.aux.coord;
    n   = size(Pts,1);
    % dummy var
    dummy_map = ones(n,1);

    figure;
    spy(W)
    xlabel(strcat(MATS{i}, ': adjacency matrix'))
    %matlab2tikz('showInfo', false, strcat('../../', MATS{i}, '_adj.tex'))
    
    %% 2a) Create the Laplacian matrix and plot the graph using the 2nd and 3rd eigenvectors
    [L,~] = CreateLapl(W);
    [V,~] = eigs(L, 4, 'smallestabs');

    % Plot and compare
    figure;
    subplot(2,2,1);
    gplot(W,Pts)
    title(strcat(MATS{i}, ': nodal coordinates'))
    subplot(2,2,2);
    gplot(W,V(:, 2:3))
    title(strcat(MATS{i}, ': eigenvector coordinates'))

    %% 2b) Cluster each graph in K = 4 clusters with the spectral and the 
    %      k-means method, and compare yourresults visually for each case.

    [D_spec,x_spec]     = kmeans_mod(V, K, n);
    [D_kmeans,x_kmeans] = kmeans_mod(Pts, K, n);

    % Compare and visualize
    subplot(2,2,3);
    gplotmap(W,Pts,x_spec)
    title(strcat(MATS{i}, ': spectral clusters'))
    subplot(2,2,4);
    gplotmap(W,Pts,x_kmeans)
    title(strcat(MATS{i}, ': k-means clusters'))
    matlab2tikz('showInfo', false, strcat('../../', MATS{i}, '_clu.tex'))

    cx = sum(x_spec == 1:4);
    ck = sum(x_kmeans == 1:4);
    fprintf('%10s spectral: %4d %4d %4d %4d k-means: %4d %4d %4d %4d\n', ...
        MATS{i}, cx(1), cx(2), cx(3), cx(4), ck(1), ck(2), ck(3), ck(4));  
    
    figure;
    subplot(1,2,1);
    title(strcat(MATS{i}, ': spectral histogram'));
    histogram(x_spec, [0:4] + 0.5);
    subplot(1,2,2);
    title(strcat(MATS{i}, ': k-means histogram'));
    histogram(x_kmeans, [0:4] + 0.5);
    matlab2tikz('showInfo', false, strcat('../../', MATS{i}, '_hist.tex'));
    
    %% 2c) Calculate the number of nodes per cluster
    [Spec_nodes,Kmeans_nodes] = ClusterMetrics(K,dummy_map,dummy_map);
end