function [Theta,logprob] = learnThetaAtFixedPhi(alphas,X,Phi,reltol,maxiter)
% Find topic probabilities for each document separately keeping word
% probabilities fixed

% Inputs:
%   alphas     - Dirichlet prior, K-by-1, double
%   X          - W-by-N sparse double matrix for W words and N documents
%   Phi        - W-by-K double matrix of word probabilities per topic for W
%                words and K topics
%   reltol     - Relative tolerance on the likelihood value for each
%                individual document, double scalar
%   maxiter    - Max number of iterations, double integer
%
% Outputs:
%   Theta      - K-by-N matrix of topic probabilities per document for K
%                topics and N documents, double
%   logprob    - N-by-1 vector of unconditional log-probabilities for
%                documents computed using input Phi and estimated Theta,
%                double

%   Copyright 2017 The MathWorks, Inc.

N = size(X,2);
K = size(Phi,2);

nj = sum(X,1);

alpha0 = sum(alphas);

if alpha0==0
    Theta = NaN(K,N,'like',Phi);
    logprob = NaN(N,1,'like',Phi);
else
    Theta = repmat(alphas/alpha0,1,N);
    logprob = zeros(N,1,'like',Phi);
end

pwz2 = sum(Phi,2)./K;

for n=1:N
    x = X(:,n);
    
    if all(x==0)
        continue;
    end
    
    phi_times_theta = pwz2;
    theta = ones(K,1,'like',Phi)./K;
    
    iter = 0;
    
    relDeltaObjective = Inf;
    
    prev_logprob = Inf;
    
    while relDeltaObjective>reltol && iter<maxiter
        iter = iter + 1;
        
        Pzdw = (theta'.*Phi)./phi_times_theta; % T-by-K
        Pzdw(phi_times_theta==0,:) = 0;
        Pzdw = x.*Pzdw;
        
        if alpha0 > 0
            theta = (sum(Pzdw,1)' + alphas)/(nj(n) + alpha0);
        else
            theta = sum(Pzdw,1)'./nj(n);
        end
        
        phi_times_theta = Phi*theta;
        
        nlogp = x.*log(phi_times_theta);
        nlogp(phi_times_theta==0 & x==0) = 0;
        
        logprob(n) = sum(nlogp);
        
        deltaObj          = logprob(n) - prev_logprob;
        prev_logprob      = logprob(n);
        relDeltaObjective = abs(deltaObj/logprob(n));
    end
    
    Theta(:,n) = theta;
    
%     if verbose>0 && mod(n,100)==0
%         fprintf('%i documents done.\n',n);
%     end
end

end
