% Benchmark for recursively partitioning meshes, based on various
% bisection approaches
%
% D.P & O.S for Numerical Computing in USI



% add necessary paths
addpaths_GP;
nlevels_a = 3;
nlevels_b = 4;

%addpath /Users/maggicl/Git/matlab2tikz/src/;

fprintf('       *********************************************\n')
fprintf('       ***  Recursive graph bisection benchmark  ***\n');
fprintf('       *********************************************\n')

% load cases
cases = {
    'airfoil1.mat';
    '3elt.mat';
    'barth4.mat';
    'mesh3e1.mat';
    'crack.mat';
    };

nc = length(cases);
maxlen = 0;
for c = 1:nc
    if length(cases{c}) > maxlen
        maxlen = length(cases{c});
    end
end

for c = 1:nc
    fprintf('.');
    sparse_matrices(c) = load(cases{c});
end


fprintf('\n\n Report Cases         Nodes     Edges\n');
fprintf(repmat('-', 1, 40));
fprintf('\n');
for c = 1:nc
    spacers  = repmat('.', 1, maxlen+3-length(cases{c}));
    [params] = Initialize_case(sparse_matrices(c));
    fprintf('%s %s %10d %10d\n', cases{c}, spacers,params.numberOfVertices,params.numberOfEdges);
end

%% Create results table
fprintf('\n%7s %16s %20s %16s %16s\n','Bisection','Spectral','Metis 5.0.2','Coordinate','Inertial');
fprintf('%10s %10d %6d %10d %6d %10d %6d %10d %6d\n','Partitions',8,16,8,16,8,16,8,16);
fprintf(repmat('-', 1, 100));
fprintf('\n');


for c = 1:nc
    spacers = repmat('.', 1, maxlen+3-length(cases{c}));
    fprintf('%s %s', cases{c}, spacers);
    sparse_matrix = load(cases{c});
    
    if c == 5
        g = 1;
    else
        g = 0;
    end
    

    % Recursively bisect the loaded graphs in 8 and 16 subgraphs.
    % Steps
    % 1. Initialize the problem
    [params] = Initialize_case(sparse_matrices(c));
    A      = params.Adj;
    coords = params.coords;
    
    W = A;
    
    % 2. Recursive routines
    % i. Spectral  
    sp8 = runtest(W, coords, @(a, b) bisection_spectral(a, b, 0), ...
         3, 'Spectral_8', g);
    sp16 = runtest(W, coords, @(a, b) bisection_spectral(a, b, 0), ...
         4, 'Spectral', g);
   
    % ii. Metis
    m8 = runtest(W, coords, @(a, b) bisection_metis(a, b, 0), ...
         3, 'METIS_8', g);
    m16 = runtest(W, coords, @(a, b) bisection_metis(a, b, 0), ...
         4, 'METIS', g);
    
    % iii. Coordinate 
    c8 = runtest(W, coords, @(a, b) bisection_coordinate(a, b, 0), ...
         3, 'Coordinate_8', g);
    c16 = runtest(W, coords, @(a, b) bisection_coordinate(a, b, 0), ...
         4, 'Coordinate', g);
    
    % iv. Inertial
    i8 = runtest(W, coords, @(a, b) bisection_inertial(a, b, 0), ...
         3, 'Inertial_8', g);
    i16 = runtest(W, coords, @(a, b) bisection_inertial(a, b, 0), ...
         4, 'Inertial', g);
    
    close;
    for ignored = 1:8
    %    close;
    end
    
    
    % 3. Calculate number of cut edges
    % 4. Visualize the partitioning result
    
    fprintf('%6d %6d %10d %6d %10d %6d %10d %6d\n',sp8,sp16,...
    m8,m16,c8,c16,i8,i16);
    
end

%Bisection         Spectral       Metis 5.0.2    Coordinate       Inertial
%Partitions          8     16       8     16       8     16       8     16
%--------------------------------------------------------------------------
%airfoil1.mat ...   327    578     320    563     516    819     578    903
%3elt.mat .......   372    671     395    651     733   1168     880   1342
%barth4.mat .....   505    758     405    689     875   1306     888   1348
%mesh3e1.mat ....    75    124     75    117      75     122      76    116
%crack.mat ......   804   1303     784   1290     1343  1860    1061   1618

function n = runtest(A, xy, algorithm, depth, t, show_graphs)
    g = 0;
    if show_graphs == 1
        g = 1;
        %title(strcat(t, ' bisection for crack with n=16'));
        hold on;
    end
    n = recurse(A, xy, algorithm, 0.5, depth, g);
    if show_graphs == 1
        hold off;
        %matlab2tikz('showInfo', false, sprintf('../../ex4_%s.tex', t));
        figure;
    end
end

function n = recurse(A, xy, algorithm, color, depth, show_g) 
    warning('off','MATLAB:eigs:IllConditionedA')
    [partA, partB] = algorithm(A, xy);
    warning('on', 'MATLAB:eigs:IllConditionedA')
    [n, ~] = cutsize(A, partA);
    A1 = A(partA, partA);
    xy1 = xy(partA, :);
    A2 = A(partB, partB);
    xy2 = xy(partB, :);
  
    depth = depth - 1;
    if depth > 0
        n = n + recurse(A1, xy1, algorithm, color * 1.5, depth, show_g);
        n = n + recurse(A2, xy2, algorithm, color * 0.5, depth, show_g);
    elseif show_g == 1
        gplot(A1, xy1);
        gplot(A2, xy2);
    end
    
    if show_g == 1
        A12 = sparse(size(A));
        A12(partA, partB) = A(partA, partB);
        A12(partB, partA) = A(partB, partA);
        gplot(A12, xy);
    end 
end