%  Paper reference:
%  M. Katsurai and S. Ono, "TrendNets: Mapping Research Trends from Dynamic Co-Word Networks via Sparse Representation," 
%  submitted to Scientometrics, 2019.

function [Scell, num_edges] = TrendNetsAdj(Wcell, lambda)

% [Scell, num_edges] = TrendNetsAdj(Wcell, 0.0001)
%    Apply Fast Sparse-Smooth Matrix Decomposition to a given matrix, resulting in sparce matrices for TrendNets

%    ------------------------------------------
%    Wcell (cell): a time series of adjacency matrices W_t (t=1, ..., T) constructed using Eq, (1). Specifically,
%    Wcell{1} = W for the first period
%    Wcell{2} = W for the second period ....
%    The size of each W is n_words * n_words.
%    The size of Wcell is 1 * T.
%  
%    lambda (double): a parameter that controls the sparsity of TrendNets.
%    ------------------------------------------
%    Scell (cell): a time series of sparse matrices S_t (t=1, ..., T), which are adjacency matrices of TrendNets.
%    The size of Scell is 1 * T.
%  
%    num_edges (double): the total number of edges in the resulting TrendNets
%    

%% Set parameters
stepsize = 0.25;
maxIter = 20;
stopcri = 0.000001;

if ~iscell(Wcell)
    error('Wcell should be cell.');
end

num_years = size(Wcell, 2);
if num_years <= 1
    error('W cell should contain multiple adjacency matrices.')
end

num_words = size(Wcell{1}, 1);
comb_num = num_words * (num_words-1) / 2;
c = zeros(comb_num, 2);
k = 1;
for i=1:1:(num_words-1)
    for j=(i+1):1:num_words
        c(k, :) = [i, j];
        k = k + 1;
    end
end
si = sub2ind([num_words, num_words], c(:,1), c(:,2)); 
M = sparse(num_years, comb_num);

for y=1:1:num_years
    M(y, :) = Wcell{y}(si);
end

M = M';
[rowsize, colsize] = size(M);

%% initialization
X = sparse(rowsize, colsize); % sparse matrix
Y = X;
t = 1;
D = @(z) z(:,[2:end, end-1])-z;
Dt = @(z) [-z(:,1), z(:,1:end-3) - z(:,2:end-2), z(:,end-2) + z(:,end) - z(:,end-1), z(:,end-1) - z(:,end)];
errors = {};
tic
for i = 1:maxIter  
    Xpre = X;
    X = Y + stepsize*Dt(D(M-Y));
    Z = spfun(@(a)(a-lambda*stepsize), abs(X));
    Z(Z<0) = 0;
    X = sign(X).*Z;
    
    tpre = t;
    t = 0.5*(1 + sqrt(1+4*(t^2)));
    
    Y = X + ((tpre-1)/t)*(X-Xpre);
    
    % stopping condition
    e = X - Xpre;
    e = norm(e(:),2);
    disp(['i = ', num2str(i), ' error = ', num2str(e,4)])
    errors{i} = e;
    if e < stopcri
        break;
    end
end
toc

X = X';
tmp = find(X>0);
num_edges = size(tmp,1);
Scell = {};
for i=1:1:size(X, 1)
    S = sparse(num_words, num_words);
    S(si) = X(i, :);
    S = S + S'; 
    Scell{i} = S;
end

