function [Pz_d_final, Pw_z_final, max_acc, max_accval, LogL_final, nIter_final, elapse_final] ...
    = LTM(X, Y, K, W, options, Pz_d, Pw_z, classnum, output_period, labeled_ids, unlabeled_ids );

% Locally-consistent Topic Modeling (LTM) using EM
%
% Notation:
% X ... (mFea x nSmp) term-document matrix (observed data)
%       X(i,j) stores number of occurrences of word i in document j
%
%       mFea  ... number of words (vocabulary size)
%       nSmp  ... number of documents
% K     ... number of topics
% W     ... weight matrix of the affinity graph 
%
% options ... Structure holding all settings
%
% You only need to provide the above four inputs.
%
% Pz_d ... P(z|d)
% Pw_z ... P(w|z) corresponds to beta parameter in LDA
%
%
% References:
% [1] Deng Cai, Xuanhui Wang, Xiaofei He, "Probabilistic Dyadic Data
% Analysis with Local and Global Consistency", ICML 2009.
%
%
% This software is based on the implementation of pLSA from
%
% Peter Gehler
% Max Planck Institute for biological Cybernetics
% pgehler@tuebingen.mpg.de
% Feb 2006
% http://www.kyb.mpg.de/bs/people/pgehler/code/index.html
%
%   
%   version 1.0 --Nov/2008 
%   
%
%   Written by Deng Cai (dengcai AT gmail.com)
%
ZERO_OFFSET = 1e-200;
useC = true;

differror = 1e-7;
if isfield(options,'error')
    differror = options.error;
end

maxIter = [];
if isfield(options, 'maxIter')
    maxIter = options.maxIter;
end

nRepeat = 1;
if isfield(options,'nRepeat')
    nRepeat = options.nRepeat;
end

minIterOrig = 100;
if isfield(options,'minIter')
    minIterOrig = options.minIter;
end
minIter = minIterOrig-1;

meanFitRatio = 0.1;
if isfield(options,'meanFitRatio')
    meanFitRatio = options.meanFitRatio;
end

alpha = 1;
if isfield(options,'alpha')
    alpha = options.alpha;
end

lambdaB = 0;
if isfield(options,'lambdaB')
    lambdaB = options.lambdaB;
end

Verbosity = 0;
if isfield(options,'Verbosity')
    Verbosity = options.Verbosity;
end

if min(min(X)) < 0
    error('Input should be nonnegative!');
end

[mFea,nSmp]=size(X);
if ~exist('Pz_d','var')
    [Pz_d,Pw_z] = pLSA_init(X,K);
else
    nRepeat = 1;
end

Pd = sum(X)./sum(X(:));
Pd = full(Pd);
Pw_d = mex_Pw_d(X,Pw_z,Pz_d);

DCol = full(sum(W,2));
D = spdiags(DCol,0,speye(size(W,1)));
L = D - W;
if isfield(options,'NormW') && options.NormW
    D_mhalf = DCol.^-.5;
    
    tmpD_mhalf = repmat(D_mhalf,1,nSmp);
    L = (tmpD_mhalf.*L).*tmpD_mhalf';
    clear D_mhalf tmpD_mhalf;
    
    L = max(L, L');
end
L = alpha*L;
dLen = full(sum(X,1));
OmegaL = (1-lambdaB)*spdiags(dLen',0,speye(size(W,1)))+L;
[R,p] = chol(OmegaL);
if p > 0
    error('L not positive definite!');
end

%fprintf( '%d: ', 0 );
%[ acc1, acc2 ] = compute_accuracies( Pz_d, N, Rd, Y, length( labeled_id ), labeled_id, unlabeled_id, classnum );

N = size( X, 2 );
max_accval = zeros( length(labeled_ids), 1 );
max_acc = zeros( length(labeled_ids), 1 );

tryNo = 0;
selectInit = 1;
nIter = 0;
LogL = [];
while tryNo < nRepeat
    tmp_T = cputime;
    tryNo = tryNo+1;
    maxErr = 1;
    i = 0;
    while(maxErr > differror)
        [Pw_z,Pz_d] = mex_EMstep(X,Pw_d,Pw_z,Pz_d);
        Pz_d = Pz_d.*repmat(dLen,K,1);
        for k = 1:K
            Pz_d(k,:) = (R\(R'\Pz_d(k,:)'))';
        end
        Pw_d = mex_Pw_d(X,Pw_z,Pz_d);
        
        %i = i + 1;
        %fprintf( '%d: ', i );
        %[ acc1, acc2 ] = compute_accuracies( Pz_d, N, Rd, Y, length( labeled_id ), labeled_id, unlabeled_id, classnum );

        nIter = nIter + 1;
        if nIter > minIter
            if selectInit
                newLogL = mex_logL(X,Pw_d,Pd);
                newLogL = newLogL - sum(sum((log(Pz_d + ZERO_OFFSET)*L).*Pz_d));
                LogL = [LogL newLogL]; %#ok<AGROW>
                maxErr = 0;
            else
                newLogL = mex_logL(X,Pw_d,Pd);
                newLogL = newLogL - sum(sum((log(Pz_d + ZERO_OFFSET)*L).*Pz_d));
                LogL = [LogL newLogL]; %#ok<AGROW>
                meanFit = meanFitRatio*meanFit + (1-meanFitRatio)*newLogL;
                maxErr = (meanFit-newLogL)/meanFit;
                
                if ~isempty(maxIter)
                    if nIter >= maxIter
                        maxErr = 0;
                    end
                end
            end
        else
            newLogL = mex_logL(X,Pw_d,Pd);
            newLogL = newLogL - sum(sum((log(Pz_d + ZERO_OFFSET)*L).*Pz_d));
            LogL = [LogL newLogL]; %#ok<AGROW>
        end
        if Verbosity
            if length(LogL) > 1
                disp(['tryNo: ',num2str(tryNo),' Iteration: ',num2str(nIter),' LogL: ',num2str(LogL(end)),' deltaLogL: ',num2str(LogL(end)-LogL(end-1)),' maxErr:',num2str(maxErr)]);
            else
                disp(['tryNo: ',num2str(tryNo),' Iteration: ',num2str(nIter),' LogL: ',num2str(LogL(end)),' maxErr:',num2str(maxErr)]);
            end
        end
        
        if mod( nIter, 10 ) == 0 && output_period > 0
            for k = 1 : length( labeled_ids )
                [ acc1, acc2, acc3, acc4, accval ] = compute_accuracies( Pz_d, N, K, Y, labeled_ids{k}, unlabeled_ids{k}, classnum );
                if accval >= max_accval(k)
                    max_accval(k) = accval;
                    max_acc(k) = acc2;
                end
            end
        end
        
        if mod( nIter, output_period ) == 0 && output_period > 0
            D = diag( sum( W ) );
            L = D - W;    
            residual = compute_residual( X, Pw_z, Pz_d, X > 0 );        
            dist_ratio = full(compute_distance_ratio( Pz_d, L, N ));
            fprintf( 1, '%d : %f %f\n', nIter, residual/N, dist_ratio );
            for k = 1 : length( labeled_ids )
                [ acc1, acc2, acc3, acc4, accval ] = compute_accuracies( Pz_d, N, K, Y, labeled_ids{k}, unlabeled_ids{k}, classnum );
                fprintf( 1, '  %d : %f %f %f\n', length(labeled_ids{k})/classnum, acc2, acc4, accval );
            end
        end
        
    end
    
    elapse = cputime - tmp_T;
    
    if tryNo == 1
        Pz_d_final = Pz_d;
        Pw_z_final = Pw_z;
        nIter_final = nIter;
        elapse_final = elapse;
        LogL_final = LogL;
        if useC
            Pw_d_final = Pw_d;
        else
            Pz_dw_final = Pz_dw;
        end
    else
        if LogL(end) > LogL_final(end)
            Pz_d_final = Pz_d;
            Pw_z_final = Pw_z;
            nIter_final = nIter;
            LogL_final = LogL;
            if selectInit
                elapse_final = elapse;
            else
                elapse_final = elapse_final+elapse;
            end
            if useC
                Pw_d_final = Pw_d;
            else
                Pz_dw_final = Pz_dw;
            end
        end
    end
    
    if selectInit
        if tryNo < nRepeat
            %re-start
            [Pz_d,Pw_z] = pLSA_init(X,K);
            if useC
                Pw_d = mex_Pw_d(X,Pw_z,Pz_d);
            else
                Pz_dw = [];
            end
            LogL = [];
            nIter = 0;
        else
            tryNo = tryNo - 1;
            minIter = 0;
            selectInit = 0;
            Pz_d = Pz_d_final;
            Pw_z= Pw_z_final;
            LogL = LogL_final;
            nIter = nIter_final;
            meanFit = LogL(end)*10;
            if useC
                Pw_d = Pw_d_final;
            else
                Pz_dw = Pz_dw_final;
            end
        end
    end
end


% [Pd,Pz_d,Pw_z,Pq,Pz_q] = pLSA_init(X,K,Xtest)
%
% initialize the probability distributions
function [Pz_d,Pw_z,Pz_q] = pLSA_init(X,K,Xtest)

[mFea,nSmp] = size(X);

Pz_d = rand(K,nSmp);
Pz_d = Pz_d ./ repmat(sum(Pz_d,1),K,1);

% random assignment
Pw_z = rand(mFea,K);
Pw_z = Pw_z./repmat(sum(Pw_z,1),mFea,1);

if nargin > 2
    nTest = size(Xtest,2);
    Pz_q = rand(K,nTest);
    Pz_q = Pz_q ./ repmat(sum(Pz_q),K,1);
end


