function [W, s, obj] = SymPLSA(A, r, max_iter, check_step, W0, s0)
% the symmetric version of Probabilistic Latent Semantic Indexing (PLSI)
% He Zhang, Sep. 2012
% For PLSI, see ref: T. Hofmann: Probabilistic latent semantic indexing, SIGIR, ACM 1999

n = size(A,1);
rand('seed', 0);
if ~exist('W0', 'var') || isempty(W0)
    W0 = rand(n,r);
end    
if ~exist('s0', 'var') || isempty(s0)
    s0 = ones(r,1)/r;
end    
W = W0;
s = s0;

W = bsxfun(@rdivide, W, sum(W)+eps);
s = s / sum(s);

for iter=1:max_iter
    Z = sp_factor_ratio(A, bsxfun(@times, W, s'), W');
    W = W .* (Z * W);
    W = bsxfun(@rdivide, W, sum(W)+eps);
    
    Z = sp_factor_ratio(A, bsxfun(@times, W, s'), W');
    s = s .* diag(W'*Z*W);
    s = s / sum(s);
    
    if mod(iter, check_step)==0
        Z = sp_factor_ratio(A, bsxfun(@times, W, s'), W');
        indnz = A>0;
        obj = sum(sum(A(indnz) .* log(Z(indnz)+eps)));
        fprintf('SymPLSA: iter=%d, obj=%6.5f\n', iter, full(obj));
    end
    
end

if nargout>2
    Z = sp_factor_ratio(A, bsxfun(@times, W, s'), W');
    indnz = A>0;
    obj = sum(sum(A(indnz) .* log(Z(indnz)+eps)));
    obj = full(obj);
end
