function [ys, times] = vfdiags(u, kernels, widths, spiopts)
%VFDIAGS computes the output of a Volterra filter in the diagonal 
%coordinate representation for a given input signal.
%
%   YS = VFDIAGS(U, KRS, WS) returns an L x P array whose columns are the
%   outputs of the homogeneous subsystems defined by the P kernels whose 
%   diagonal coordinate representations are stored in the cell array 
%   KERS. WS(p) specifies the width of the pth-order kernel representation,
%   that is, the maximum distance around the main kernel diagonal up to
%   which all diagonals are kept.
%
%   [YS, TIMES] = VFILTER(U, KRS, WS) returns also a vector TIMES
%   containing the time spent on the computation of each homogeneous
%   subsystem output.
%
%   YS = VFDIAGS(U, KRS, WS, SPIOPTS) allows applying a sparse-interpolated
%   strategy where the kernel diagonals are sparsified and the output is
%   filtered by a reconstruction filter SPIOPTS.g, which mitigates the
%   effect of zeroing kernel coefficients. In this case, SPIOPTS.L contains 
%   the interpolation factor, and only the nonzero coefficients are stored 
%   in the diagonals contained by the cell array kernels.
%
%   The diagonal coordinate representation of the kernels must be computed
%   from the vectorized Volterra kernels by means of the function kvdiags.
%
%   See also VFDIAGSP, KVDIAGS
%
%   CBVK Toolbox
%   JHMG, PMSB 2016

if nargin < 4
    spiopts = [];
end

% Model order and memories
P = length(kernels);

% Output signals and measured times
yps = cell(1,P);
times = zeros(1,P);

% Loops over pth-order subsystems
inttimer = tic;
yps{1} = convsi(kernels{1},u,spiopts);
times(1) = toc(inttimer);
for p=2:P
    inttimer = tic;
    yps{p} = vfdiagsp(u, kernels{p}, p, widths(p), spiopts);
    times(p) = toc(inttimer);
end

% Arranges outputs in array with correct size
MM = max(cellfun(@length,yps));
ys = zeros(MM,P);
for p=1:P
    ys(:,p) = [yps{p}(:); zeros(MM-length(yps{p}(:)))];
end




%memories = cellfun(@(x) length(x{1}),kernels);
%
% Temp fix
% memories=length(kernels{1})*ones(1,P);

% Default width = memory
% if nargin < 3
%     widths = memories;
% end;    
% 
% u = u(:);
% 
% Nu = length(u);
% N = max(memories);
% Ny = Nu+N-1;
% 
% % Initialization
% ys = zeros(Ny,P);
% times = zeros(P,1);
% 
% inttimer = tic;
% 
% % First order output
% %ylin = conv(kernels{1}{1},u);
% ylin = conv(kernels{1},u);
% ys(:,1) = [ylin; zeros(Ny-length(ylin),1)];
% times(1) = toc(inttimer);
% 
% if P==1
%     return;
% end;
% 
% % Array of input powers
% inppow{N} = zeros(1,Nu);
% for n=1:N                                         
%     inppow{n} = [u(n:end); zeros(n-1,1)]; 
% end;
% 
% for p = 2:P
%     inttimer = tic;
%     
%     W = widths(p);
%            
%     wp = zeros(1,p-1);
%     ndiag = 1;
%     while (~isempty(wp)) && (ndiag <= length(kernels{p}))
%                
%         h = kernels{p}{ndiag};
%         
%         inds = [wp(end)-wp(1:end-1) wp(end)]+1;
%         us = u;
%         for ii = inds
%             us = us.*inppow{ii};
%         end;
%         
%         h = [zeros(wp(end),1); h(:)];
%         yparcial = conv(h,us);
% 
%         ys(:,p) = ys(:,p) + [yparcial zeros(1,Ny-length(yparcial))];
%         
%         wp = nextinds(wp,W);
%         ndiag = ndiag+1;
%     end;
%     times(p) = toc(inttimer);
% end;
% 
