%DEMO script for commGSP
%commGSP  is a toolbox for community-aware graph signal processing (GSP) for MATLAB.
%This script calls all other functions of the toolbox (beginning with "commGSP_")
%and shows results of community-aware GSP on a toy graph.
%The toolbox was developed by Miljan Petrovic, at Medical Image Processing
%Lab, EPFL. It follows the framework proposed in:
% M. Petrovic, R. Liegeois, T. A. W. Bolton, D. Van De Ville, "Community-aware
% Graph Signal Processing", Signal Processing Magazine: Graph Signal
% Processing: Foundations And Emerging Directions, Nov. 2020.
% DOI: 10.1109/MSP.2020.3018087

clc; close all; clear all;
%build adjacency matrix of a graph
A=[0 0 1 0 0 1; 0 0 1 0 0 1; 1 1 0 0 0 0 ; 1 1 0 0 0 0; ...
    1 0 0 1 0 0; 1 1 0 0 0 1];
A=A+A';
A=[A, zeros(size(A));  zeros(size(A)), A]; A(1,end)=1; A(end,1)=1;
N=size(A,1);
A=sparse(A);

%build different shift operators & eigendecompose
L = commGSP_build_matrix(A, 'laplacian');
Q = commGSP_build_matrix(A, 'modularity');
[UL, DL] = commGSP_eig(L, 'ascend');
[UA, DA] = commGSP_eig(A, 'descend');
[UQ, DQ] = commGSP_eig_implicit(A, 'modularity');

%graph fourier transform
x=sparse(1:N)'; cutoff=4;
x_hat_1 = commGSP_gft(A, x, 'modularity'); %GFT of a signal (back-end eigendecompositon)
x_hat_2 = commGSP_gft(A, x, UQ); %GFT of a signal (w/o back-end eigendecompositon)
Q_x = commGSP_signal_modularity(A, [x, rand(N,1)]);

%filtering
H=sparse(1:cutoff,1:cutoff,ones(cutoff,1),N,N);
y_1 = commGSP_filter_spectral(A, x, H, 'modularity');
p = commGSP_spectral2poly(A, H, 6, 'modularity',1);
y_2 = commGSP_filter_poly(A, x, p, 'modularity');

%optimal subsampling and reconstruction
BW=8; K=8;
S = commGSP_optimal_subgraph(A, 'modular', BW, K);
x_sub=S*x;
y_3 = commGSP_reconstruct(A, 'modular', BW, S, x_sub);

%denoising
ns=randn(N,1); ns=(ns-mean(ns))./std(ns); ns=ns.*std(x)./1;
pens=[0.1 1];
[y_4, mu] = commGSP_quad_regularization(A, x, x+ns, pens, 'modular');

%surrogate generation
Ns=3;
SG = commGSP_surrogate(A, x, 'modular', Ns);

%% plotting
g=graph(A); fsz=14; msz=10; lw=2;

display(['"Ramp" signal''s modularity measure on the graph is ', num2str(Q_x(1))])
display(['Random signal''s modularity measure on the graph is ', num2str(Q_x(2))])

figure('Name','Graph Matrices & Spectra')
pp=get(gcf,'Position'); pp(1:2)=0; pp(3)=pp(3).*2.2; set(gcf,'Position',pp);
subplot 231
    imagesc(A)
    title('Adjacency matrix')
    colorbar
    set(gca,'FontSize',fsz)
subplot 232
    imagesc(L)
    title('Laplacian matrix')
    colorbar
    set(gca,'FontSize',fsz)
subplot 233
    imagesc(Q)
    title('Modularity matrix')
    colorbar
    set(gca,'FontSize',fsz)
subplot 234
    plot(diag(DA),'k','LineWidth',lw)
    xlim([1 N])
    title('Adjacency eigenvalues')
    set(gca,'FontSize',fsz)
subplot 235
    plot(diag(DL),'b','LineWidth',lw)
    xlim([1 N])
    title('Laplacian eigenvalues')
    set(gca,'FontSize',fsz)
subplot 236
    plot(diag(DQ),'r','LineWidth',lw)
    xlim([1 N])
    title('Modularity eigenvalues')
    set(gca,'FontSize',fsz)

figure('Name','Graph Signals & Operator Results')
pp=get(gcf,'Position'); pp(3)=pp(3).*2.2; set(gcf,'Position',pp);
subplot 241
    plot(g, 'NodeLabel', repmat({''},N,1), 'NodeCData',x, 'MarkerSize',msz)
    title('Graph signal')
    colorbar
    set(gca,'FontSize',fsz)
subplot 242
    plot(g, 'NodeLabel', repmat({''},N,1), 'NodeCData',full(x_sub), 'MarkerSize',msz)
    title('Subsampled signal')
    colorbar
    set(gca,'FontSize',fsz)
subplot 243
    plot(g, 'NodeLabel', repmat({''},N,1), 'NodeCData',y_3, 'MarkerSize',msz)
    title('Reconstructed signal')
    colorbar
    set(gca,'FontSize',fsz)
subplot 244
    plot(g, 'NodeLabel', repmat({''},N,1), 'NodeCData',SG(:,1), 'MarkerSize',msz)
    title('Surrogate signal')
    colorbar
    set(gca,'FontSize',fsz)
subplot 245
    plot(g, 'NodeLabel', repmat({''},N,1), 'NodeCData',x+ns, 'MarkerSize',msz)
    title('Noisy signal')
    colorbar
    set(gca,'FontSize',fsz)
subplot 246
    plot(g, 'NodeLabel', repmat({''},N,1), 'NodeCData',y_4, 'MarkerSize',msz)
    title('Denoised signal')
    colorbar
    set(gca,'FontSize',fsz)
subplot 247
    plot(g, 'NodeLabel', repmat({''},N,1), 'NodeCData',y_1, 'MarkerSize',msz)
    title('Filtered signal (spectral filtering)')
    colorbar
    set(gca,'FontSize',fsz)
subplot 248
    plot(g, 'NodeLabel', repmat({''},N,1), 'NodeCData',y_2, 'MarkerSize',msz)
    title('Filtered signal (polynomial filtering)')
    colorbar
    set(gca,'FontSize',fsz)