function [x,r,hist,nverts] = ppr_push_queue(P,alpha,v,tol,maxsteps)
% PPR_PUSH_QUEUE Implement a push strategy for PPR ranking
%
% [x,r,hist,nverts] = ppr_push_queue(P,alpha,b,tol,maxsteps) specifies
%   P - the random walk transition matrix
%   alpha - the value of alpha in the Katz matrix
%   v - the teleportation vector.  See below for the interpretation.
%   tol - a stopping tolerance (default = 1e-7);
%   maxsteps - a limit on the total number of steps (default = 10*size(A,1))
% and returns
%   x - the ppr ranking vector given by the solution of
%     (I-alpha*P')x = (1-alpha)*v
%   r - the residual vector 
%   hist - a history of the residuals at each step and the vertex pushed
%   nverts - the total number of unique vertices explored
%   
% Solve (I-alpha*P')*x = (1-alpha)*v by exploiting the graph nature of
% P to only access vertex out-degrees.
%
% If v has length==N, then it is a full right hand side, otherwise, it 
% is a list of vertices and we set v = (1-alpha)/|b| for all these vertices and
% 0 elsewhere.

% David F. Gleich
% Purdue University, 2011

% History
% -------
% :2010-11-03: Modified from katz_push.m from local Katz codes

if ~exist('tol','var') || isempty(tol), tol = 1e-7; end
if ~exist('maxsteps','var') || isempty(maxsteps), maxsteps = 10*size(A,1); end

% Algorithm: We store the non-zero elements of the residual in a queue and
% push elements from the queue.
% convert graph to csr
Pt = P'; % get fast access to columns of Pt;
N=size(P,1); % NUMBER of verties
x = zeros(N,1); % x is the ranking vector
r = zeros(N,1); % r is the residual vector
if length(v) == N, r=(1-alpha)*b(:); 
else r(v) = (1-alpha)/numel(v);
end

% queue structure
q=zeros(N,1); qsize=0; qhead=1; qtail=0; qind=false(N,1);

% setup book keeping for the history
nedges = 0;
hist = zeros(maxsteps,4); % nedges, sumresid, time, vertex
dt = tic; % start the timer
visited = zeros(N,1);
sumresid = 0;

for w=find(r)' % add elements from the current residual to the queue
    sumresid = sumresid + r(w); 
    assert(qsize < N); qtail=mod(qtail,N)+1; q(qtail)=w; 
    qsize=qsize+1; qind(w)=1;
end

step = 1;
for step=1:maxsteps
    assert(qsize>0);
    v=q(qhead); qhead=mod(qhead,N)+1; qind(v) = false; qsize=qsize-1; % pop queue head
    visited(v) = 1;
    x(v) = x(v) + r(v); % increment v's pagerank    
    val = alpha*r(v);% set the value to push out to others and kill our residual
    sumresid = sumresid - r(v); % decrement the residual sum
    r(v) = 0;
    % for each vertex adjacent to v, push to it!
    [neighs,~,probs] = find(Pt(:,v)); % find neighboring elements
    for ni=1:length(neighs)
        w = neighs(ni); pw = probs(ni);        
        r(w)=r(w)+val*pw; % increase the residual
        sumresid = sumresid + val*pw; % increase its sum
        if ~qind(w) && r(w)>tol % add w to the queue
            assert(qsize < N); qtail=mod(qtail,N)+1; q(qtail)=w; 
            qsize=qsize+1; qind(w)=1;
        end
    end
    nedges = nedges + numel(neighs);
    
    % store history
    hist(step,1) = nedges;
    hist(step,2) = sumresid;
    hist(step,3) = toc(dt);
    hist(step,4) = v;
    
    if qsize == 0, break; end % stop iterating!
end

hist = hist(1:step,:); % truncate hist
nverts = sum(visited);