classdef (Sealed) ldaModel < matlab.mixin.CustomDisplay
%

%   Copyright 2017-2022 The MathWorks, Inc.

    properties(GetAccess=public,SetAccess=?tldaModel)
        FitInfo
        TopicConcentration
        CorpusTopicProbabilities
        WordConcentration
    end

    properties(GetAccess=?tldaModel,SetAccess=?tldaModel)
        FitAlpha0                      % logical scalar
        FitAlphas                      % logical scalar
        ModelParameters                % record model parameters
        Nkj                            % K-by-N
        Nwk                            % W-by-K
        PrivStochasticTheta            % K-by-N

        % The names for these two private properties follow the
        % dimensionality convention: word for the first dimension goes
        % first etc.
        PrivWordTopicProbabilities     % W-by-K
        PrivTopicDocumentProbabilities % K-by-N

        Wordlist                       % wordlist behind Vocabulary

        PrivTopicOrder                 % toggle to sort or not
        PrivTopicOrderVector           % sort order for topics
    end

    properties(Dependent)
        TopicOrder
    end

    properties(Dependent,GetAccess=public,SetAccess=?tldaModel)
        DocumentTopicProbabilities
        NumTopics
        TopicWordProbabilities
        Vocabulary
    end

    methods
        function thetaT = get.DocumentTopicProbabilities(mdl)
            thetaT = mdl.PrivTopicDocumentProbabilities';
        end

        function numTopics = get.NumTopics(mdl)
            numTopics = numel(mdl.CorpusTopicProbabilities);
        end

        function phi = get.TopicWordProbabilities(mdl)
            phi = mdl.PrivWordTopicProbabilities;
        end

        function s = get.Vocabulary(mdl)
            wl = mdl.Wordlist;
            s = wl.Vocabulary;
        end

        function order = get.TopicOrder(mdl)
            order = mdl.PrivTopicOrder;
        end

        function mdl = set.TopicOrder(mdl, order)
            vals = {'initial-fit-probability','unordered'};
            order = validatestring(order,vals);
            if ~isequal(order, mdl.PrivTopicOrder)
                perm = mdl.PrivTopicOrderVector;
                if strcmp(order,vals{1})
                    mdl = sortTopics(mdl, perm);
                else
                    mdl = unsortTopics(mdl, perm);
                end
                mdl.PrivTopicOrder = order;
            end
        end
    end


    % The following notation is used in this code:
    %   W       - Dictionary size
    %   N       - Number of documents
    %   K       - Number of topics
    %   Phi     - W-by-K matrix of word probabilities per topic
    %   Nwk     - W-by-K matrix of word counts per topic
    %   nk      - 1-by-K vector of aggregated word counts per topic
    %   Theta   - K-by-N matrix of topic probabilities per document
    %   Nkj     - K-by-N matrix of topic counts per document
    %   nj      - 1-by-N vector of word counts per document


    methods(Access=protected)
        function group = getPropertyGroups(~)
            plist = {'NumTopics' 'WordConcentration' 'TopicConcentration' ...
                     'CorpusTopicProbabilities' ...
                     'DocumentTopicProbabilities' 'TopicWordProbabilities' ...
                     'Vocabulary' 'TopicOrder' 'FitInfo'};
            group = matlab.mixin.util.PropertyGroup(plist);
        end
    end


    methods(Access=?tldaModel)
        function mdl = ldaModel()
        end

        function mdl = initializeTopicOrder(mdl)
            mdl.PrivTopicOrder = 'unordered';
            alphas = mdl.CorpusTopicProbabilities;
            [~,mdl.PrivTopicOrderVector] = sort(alphas,'descend');
        end

        function mdl = sortTopics(mdl,order)
            mdl.CorpusTopicProbabilities = mdl.CorpusTopicProbabilities(order);
            mdl.PrivWordTopicProbabilities = mdl.PrivWordTopicProbabilities(:,order);
            mdl.PrivTopicDocumentProbabilities = mdl.PrivTopicDocumentProbabilities(order,:);
            mdl.Nkj = mdl.Nkj(order,:);
            mdl.Nwk = mdl.Nwk(:,order);
        end

        function mdl = unsortTopics(mdl,order)
            mdl.CorpusTopicProbabilities(order) = mdl.CorpusTopicProbabilities;
            mdl.PrivWordTopicProbabilities(:,order) = mdl.PrivWordTopicProbabilities;
            mdl.PrivTopicDocumentProbabilities(order,:) = mdl.PrivTopicDocumentProbabilities;
            mdl.Nkj(order,:) = mdl.Nkj;
            mdl.Nwk(:,order) = mdl.Nwk;
        end

        function [score,ranked] = wordScores(mdl,scaling)
            [W,K] = size(mdl.PrivWordTopicProbabilities);

            if strcmp(scaling,'inversemean')
                Phi = mdl.PrivWordTopicProbabilities';

                iszero = Phi==0;

                logPhi = log(Phi);
                logPhi = logPhi - mean(logPhi,1);

                Phi = Phi.*logPhi;
                Phi(:,any(iszero,1)) = NaN;
                Phi(iszero) = 0;

                score = Phi';
            else
                score = mdl.PrivWordTopicProbabilities;
            end

            if nargout > 1
                ranked = zeros(W,K);
                for k=1:K
                    [~,ranked(:,k)] = sort(score(:,k),'descend');
                end
            end
        end
    end


    methods(Static,Hidden)
        function mdl = fit(fun,X,K,varargin)
            [varargin{:}] = convertStringsToChars(varargin{:});

            validateattributes(X,{'bagOfWords','bagOfNgrams','numeric'},{});

            p = inputParser;
            p.CaseSensitive      = false;
            p.FunctionName       = fun;
            p.PartialMatching    = true;

            p.addRequired('FunctionName');
            p.addRequired('InputData');
            p.addRequired('NumTopics');

            p.addParameter('LogLikelihoodTolerance',1e-4);
            p.addParameter('Verbose',1);
            p.addParameter('IterationLimit',[]);
            p.addParameter('InitialTopicConcentration',[]);
            p.addParameter('WordConcentration',[]);
            p.addParameter('Solver','cgs');
            p.addParameter('FitTopicConcentration',[]);
            p.addParameter('FitTopicProbabilities',[]);
            p.addParameter('ValidationData',[]);
            p.addParameter('ValidationFrequency',[]);
            p.addParameter('Stream',[]); % RandStream object
            p.addParameter('DataPassLimit',[]);
            p.addParameter('MiniBatchLimit',[]);
            p.addParameter('MiniBatchSize',[]);
            p.addParameter('LearnRateDecay',[]);
            p.addParameter('DocumentsIn','rows');
            order_vals = {'initial-fit-probability','unordered'};
            p.addParameter('TopicOrder',order_vals{1},@(s)any(validatestring(s,order_vals)));

            p.parse(fun,X,K,varargin{:});

            reltol      = p.Results.LogLikelihoodTolerance;
            verbose     = p.Results.Verbose;
            maxiter     = p.Results.IterationLimit;
            alpha0      = p.Results.InitialTopicConcentration;
            beta0       = p.Results.WordConcentration;
            solver      = p.Results.Solver;
            dofitalpha0 = p.Results.FitTopicConcentration;
            dofitalphas = p.Results.FitTopicProbabilities;
            Xval        = p.Results.ValidationData;
            freqval     = p.Results.ValidationFrequency;
            rsh         = p.Results.Stream;
            maxpass     = p.Results.DataPassLimit;
            maxbatch    = p.Results.MiniBatchLimit;
            batchsize   = p.Results.MiniBatchSize;
            decay       = p.Results.LearnRateDecay;
            docsIn      = p.Results.DocumentsIn;
            topicOrder  = p.Results.TopicOrder;
            words       = [];

            usemex = true;
            K = cast(K,'double');

            if isa(X,'bagOfWords') || isa(X,'bagOfNgrams')
                if isa(X,'bagOfWords') && ~isscalar(X)
                    error(message('textanalytics:ldaModel:NonscalarBagOfWords'));
                elseif isa(X,'bagOfNgrams') && ~isscalar(X)
                    error(message('textanalytics:ldaModel:NonscalarBagOfNgrams'));
                end
                words = X.Wordlist;
                X = X.getCountsTransposed;
            else
                try
                    X = clustering.ldautils.prepareData(fun,'Train',X,true,docsIn);
                catch ME
                    throw(ME)
                end
            end
            if isempty(words)
                words = textanalytics.internal.wordlist(string(1:size(X,1)));
            end
            vocab = words.Vocabulary;
            vocab = vocab(:);

            if isempty(vocab)
                error(message('textanalytics:ldaModel:EmptyTrainingData'));
            end

            [W,N] = size(X);

            if ~isempty(Xval)
                if isa(Xval,'bagOfWords') || isa(Xval,'bagOfNgrams')
                    if isa(Xval,'bagOfWords') && ~isscalar(Xval)
                        error(message('textanalytics:ldaModel:NonscalarBagOfWordsValidation'));
                    elseif isa(Xval,'bagOfNgrams') && ~isscalar(Xval)
                        error(message('textanalytics:ldaModel:NonscalarBagOfNgramsValidation'));
                    end
                    if isa(Xval,'bagOfWords')
                        Xval = Xval.setVocabulary(words);
                    else
                        Xval = Xval.setNgrams(words);
                    end

                    Xval = Xval.getCountsTransposed;
                else
                    try
                        Xval = clustering.ldautils.prepareData(...
                            fun,'Validation',Xval,false,docsIn);
                    catch ME
                        throw(ME)
                    end
                end

                if size(Xval,1)~=W
                    error(message('textanalytics:ldaModel:ValidationDataSize',W));
                end
            end

            try
                [eff_solver,dosc,dotensor,maxiter,alpha0,beta0,...
                 dofitalpha0,dofitalphas,usemex,batchsize,decay,freqval] = ...
                 clustering.ldautils.checkInputs(...
                    fun,X,Xval,K,solver,usemex,maxiter,maxpass,maxbatch,...
                    alpha0,beta0,dofitalpha0,dofitalphas,batchsize,decay,freqval);
            catch ME
                throw(ME)
            end

            validateattributes(reltol,{'single','double'},...
                               {'scalar','real'},fun,'LogLikelihoodTolerance');

            validateattributes(verbose,{'numeric'},...
                               {'scalar','real','integer','>=',0,'<=',1},fun,'Verbose');

            % Initialize alphas
            alphas = alpha0*ones(K,1)/K;

            % Get random stream
            if isempty(rsh)
                rsh = RandStream.getGlobalStream;
            end

            % Initialize Phi and Theta
            switch eff_solver
              case 'cgs'
                % Initialize topic assignments. Record rows, columns
                % and topic values (z) as 3 vectors of the same length.
                tstart = tic;
                [rowidx,colidx,z,Nwk,Nkj] = ...
                    clustering.ldautils.sample_topics(X,K,rsh);
                telapsed = toc(tstart);

                if verbose > 0
                    fprintf('%s\n',...
                            message('textanalytics:ldaModel:ReportInitialSamplingTime',...
                                    sprintf('%g',telapsed)).getString());
                end

                Phi = clustering.ldautils.nwk2phi(Nwk,beta0);
                Theta = clustering.ldautils.nkj2theta(Nkj,alphas);

              case {'avb' 'cvb0'}
                % Initial topic assignments
                Theta = rand(rsh,K,N,class(X));
                Theta = Theta./sum(Theta,1);

                % Initial word probabilities
                Phi = rand(rsh,W,K,class(X));
                Phi = Phi./sum(Phi,1);
            end

            % Initialize parameters for stochastic learning
            Theta_sc = [];
            scinfo = [];
            if dosc
                scinfo                            = struct;
                scinfo.MiniBatchSize              = batchsize;
                scinfo.LearnRateDecay             = decay;
                scinfo.PermutedIndices            = randperm(rsh,N);
                scinfo.LastProcessedDocumentIndex = 0;
            end

            % Solve
            switch eff_solver
              case {'avb' 'cvb0'}
                [alphas,Theta,Phi,fitted_Nkj,fitted_Nwk,...
                 fitinfo,scinfo,Theta_sc] = ...
                 clustering.ldautils.solve_em(...
                    eff_solver,dofitalpha0,dofitalphas,alphas,beta0,...
                    X,Xval,Theta,Phi,[],[],true,reltol,0,maxiter,verbose,...
                    usemex,dosc,scinfo,freqval,rsh,[],[],[],[]);

              case 'cgs'
                [alphas,Theta,Phi,fitted_Nkj,fitted_Nwk,fitinfo,~,~] = ...
                    clustering.ldautils.solve_em(...
                    eff_solver,dofitalpha0,dofitalphas,alphas,beta0,...
                    X,Xval,Theta,Phi,Nkj,Nwk,true,reltol,0,maxiter,verbose,...
                    usemex,dosc,[],freqval,rsh,rowidx,colidx,z,[]);
            end

            [alphas,Phi,Theta,fitted_Nwk,fitted_Nkj] = ...
                zeroOutBadTopics(alphas,Phi,Theta,fitted_Nwk,fitted_Nkj);

            mdl = ldaModel;

            mp.MiniBatchSize                   = batchsize;
            mp.FitTopicConcentration           = dofitalpha0;
            mp.FitTopicProbabilities           = dofitalphas;
            mp.LearnRateDecay                  = decay;
            mp.InitialTopicConcentration       = alpha0;
            mp.IterationLimit                  = maxiter;
            mp.LogLikelihoodTolerance          = reltol;

            mdl.ModelParameters                = mp;

            mdl.PrivTopicDocumentProbabilities = Theta;
            mdl.PrivWordTopicProbabilities     = Phi;

            if ~dotensor
                mdl.Nkj                        = fitted_Nkj;
                mdl.Nwk                        = fitted_Nwk;
            end

            alpha0 = sum(alphas);
            if alpha0 > 0
                mdl.TopicConcentration         = alpha0;
                mdl.CorpusTopicProbabilities   = alphas'/alpha0;
            else
                mdl.TopicConcentration         = 0;
                mdl.CorpusTopicProbabilities   = ones(1,K)/K;
            end

            mdl.WordConcentration              = beta0;

            mdl.FitInfo                        = fitinfo;

            if dosc
                mdl.PrivStochasticTheta        = Theta_sc;
                mdl.FitInfo.StochasticInfo     = scinfo;
            end

            mdl.Wordlist                       = words;

            mdl = initializeTopicOrder(mdl);
            mdl.TopicOrder                     = topicOrder;
        end

        function mdl = loadobj(mdl)
            if isempty(mdl.PrivTopicOrder)
                % loading a MAT file pre-18b
                mdl = initializeTopicOrder(mdl);
            end
        end
    end


    methods
        function mdl = resume(mdl,X,varargin)
            [varargin{:}] = convertStringsToChars(varargin{:});

            validateattributes(X,{'bagOfWords','bagOfNgrams','numeric'},{});

            solver = mdl.FitInfo.Solver;

            if strcmp(solver,'cgs')
                error(message('textanalytics:ldaModel:ResumeForCGSNotSupported'));
            end

            old_order = mdl.PrivTopicOrder;
            mdl.TopicOrder = 'unordered';

            alpha0 = mdl.TopicConcentration;
            alphas = alpha0*mdl.CorpusTopicProbabilities';
            K = numel(alphas);

            beta0 = mdl.WordConcentration;

            dofitalpha0 = mdl.ModelParameters.FitTopicConcentration;
            dofitalphas = mdl.ModelParameters.FitTopicProbabilities;

            p = inputParser;
            p.CaseSensitive      = false;
            p.FunctionName       = 'resume';
            p.PartialMatching    = true;

            p.addRequired('Model');
            p.addRequired('InputData');

            p.addParameter('LogLikelihoodTolerance',1e-4);
            p.addParameter('Verbose',1);
            p.addParameter('IterationLimit',[]);
            p.addParameter('DataPassLimit',[]);
            p.addParameter('MiniBatchLimit',[]);
            p.addParameter('FitTopicConcentration',dofitalpha0);
            p.addParameter('FitTopicProbabilities',dofitalphas);
            p.addParameter('ValidationData',[]);
            p.addParameter('ValidationFrequency',[]);
            p.addParameter('Stream',[]); % RandStream object
            p.addParameter('DocumentsIn','rows');

            p.parse(mdl,X,varargin{:});

            reltol      = p.Results.LogLikelihoodTolerance;
            verbose     = p.Results.Verbose;
            maxiter     = p.Results.IterationLimit;
            maxpass     = p.Results.DataPassLimit;
            maxbatch    = p.Results.MiniBatchLimit;
            dofitalpha0 = p.Results.FitTopicConcentration;
            dofitalphas = p.Results.FitTopicProbabilities;
            Xval        = p.Results.ValidationData;
            freqval     = p.Results.ValidationFrequency;
            rsh         = p.Results.Stream;
            docsIn      = p.Results.DocumentsIn;

            usemex = true;

            if isa(X,'bagOfWords') || isa(X,'bagOfNgrams')
                if isa(X,'bagOfWords') && ~isscalar(X)
                    error(message('textanalytics:ldaModel:NonscalarBagOfWords'));
                elseif isa(X,'bagOfNgrams') && ~isscalar(X)
                    error(message('textanalytics:ldaModel:NonscalarBagOfNgrams'));
                end
                vocab = mdl.Vocabulary(:);

                if isa(X,'bagOfWords')
                    new_vocab = X.Vocabulary(:);
                else
                    new_vocab = X.Wordlist.Vocabulary';
                end

                if numel(vocab)~=numel(new_vocab)
                    error(message('textanalytics:ldaModel:TrainingVocabularySizeForResume',numel(vocab)));
                end

                if ~all(vocab==new_vocab)
                    error(message('textanalytics:ldaModel:TrainingVocabularyWordMismatchForResume'));
                end

                X = X.getCountsTransposed;
            else
                try
                    X = clustering.ldautils.prepareData('resume','Train',X,true,docsIn);
                catch ME
                    throw(ME)
                end
            end

            [W,N] = size(X);

            if ~isempty(Xval)
                if isa(Xval,'bagOfWords') || isa(Xval,'bagOfNgrams')
                    if isa(Xval,'bagOfWords') && ~isscalar(Xval)
                        error(message('textanalytics:ldaModel:NonscalarBagOfWordsValidation'));
                    elseif isa(Xval,'bagOfNgrams') && ~isscalar(Xval)
                        error(message('textanalytics:ldaModel:NonscalarBagOfNgramsValidation'));
                    end
                    if isa(Xval,'bagOfWords')
                        Xval = Xval.setVocabulary(mdl.Wordlist);
                    else
                        Xval = Xval.setNgrams(mdl.Wordlist);
                    end

                    Xval = Xval.getCountsTransposed;
                    docsIn = 'columns';
                end
                try
                    Xval = clustering.ldautils.prepareData(...
                        'resume','Validation',Xval,false,docsIn);
                catch ME
                    throw(ME)
                end

                if size(Xval,1)~=W
                    error(message('textanalytics:ldaModel:ValidationDataSize',W));
                end
            end

            try
                [eff_solver,dosc,~,maxiter,~,~,dofitalpha0,dofitalphas,usemex,~,~,freqval] = ...
                    clustering.ldautils.checkInputs(...
                    'resume',X,Xval,K,solver,usemex,maxiter,maxpass,maxbatch,...
                    alpha0,beta0,dofitalpha0,dofitalphas,[],[],freqval);
            catch ME
                throw(ME)
            end

            Phi = mdl.PrivWordTopicProbabilities;

            if W~=size(Phi,1)
                error(message('textanalytics:ldaModel:TrainingDataNumWordsMismatch',size(Phi,1)));
            end

            scinfo = [];
            if dosc
                scinfo = mdl.FitInfo.StochasticInfo;
                Theta = mdl.PrivStochasticTheta;
            else
                Theta = mdl.PrivTopicDocumentProbabilities;
            end

            Nexp = size(Theta,2);
            if N~=Nexp
                error(message('textanalytics:ldaModel:TrainingDataNumDocumentsMismatch',Nexp));
            end

            dohist = ~isempty(mdl.FitInfo.History);

            first_iter = numel(mdl.FitInfo.History.Perplexity);

            if isempty(rsh)
                rsh = RandStream.getGlobalStream;
            end

            validateattributes(reltol,{'single','double'},...
                               {'scalar','real'},'resume','LogLikelihoodTolerance');

            validateattributes(verbose,{'numeric'},...
                               {'scalar','real','integer','>=',0,'<=',1},'resume','Verbose');

            [alphas,Theta,Phi,fitted_Nkj,fitted_Nwk,fitinfo,scinfo,Theta_sc] = ...
                clustering.ldautils.solve_em(...
                eff_solver,dofitalpha0,dofitalphas,alphas,mdl.WordConcentration,...
                X,Xval,Theta,Phi,mdl.Nkj,mdl.Nwk,...
                dohist,reltol,first_iter,maxiter,verbose,usemex,...
                dosc,scinfo,freqval,rsh,[],[],[],mdl.FitInfo.History);

            [alphas,Phi,Theta,fitted_Nwk,fitted_Nkj] = ...
                zeroOutBadTopics(alphas,Phi,Theta,fitted_Nwk,fitted_Nkj);

            mp                                 = mdl.ModelParameters;

            mp.FitTopicConcentration           = dofitalpha0;
            mp.FitTopicProbabilities           = dofitalphas;
            mp.IterationLimit                  = maxiter;
            mp.LogLikelihoodTolerance          = reltol;

            mdl.ModelParameters                = mp;

            mdl.PrivTopicDocumentProbabilities = Theta;
            mdl.PrivWordTopicProbabilities     = Phi;

            mdl.Nkj                            = fitted_Nkj;
            mdl.Nwk                            = fitted_Nwk;

            alpha0 = sum(alphas);
            if alpha0 > 0
                mdl.TopicConcentration         = alpha0;
                mdl.CorpusTopicProbabilities   = alphas'/alpha0;
            else
                mdl.TopicConcentration         = 0;
                mdl.CorpusTopicProbabilities   = ones(1,K)/K;
            end

            mdl.FitInfo                        = fitinfo;

            if dosc
                mdl.PrivStochasticTheta       = Theta_sc;
                mdl.FitInfo.StochasticInfo    = scinfo;
            end

            mdl.TopicOrder = old_order;
        end


        function Theta = transform(mdl,documents,varargin)
            [~,Theta] = predict(mdl,documents,varargin{:});
        end


        function [cidx,Theta] = predict(mdl,documents,varargin)
            [varargin{:}] = convertStringsToChars(varargin{:});

            p = inputParser;
            p.CaseSensitive      = false;
            p.FunctionName       = 'predict';
            p.PartialMatching    = true;

            p.addRequired('Model');
            p.addRequired('InputData');

            p.addParameter('LogLikelihoodTolerance',1e-4);
            p.addParameter('IterationLimit',100);
            p.addParameter('DocumentsIn','rows');

            p.parse(mdl,documents,varargin{:});

            reltol      = p.Results.LogLikelihoodTolerance;
            maxiter     = p.Results.IterationLimit;
            docsIn      = p.Results.DocumentsIn;

            if isa(documents,'bagOfWords')
                if ~isscalar(documents)
                    error(message('textanalytics:ldaModel:NonscalarBagOfWords'));
                end
                documents = setVocabulary(documents, mdl.Wordlist);
                documents = documents.getCountsTransposed;
                docsIn = 'columns';
            end

            if isa(documents,'bagOfNgrams')
                if ~isscalar(documents)
                    error(message('textanalytics:ldaModel:NonscalarBagOfNgrams'));
                end
                [~,wl] = ...
                    textanalytics.internal.encodeNgramsWordlist(mdl.Wordlist,documents);
                documents = setNgrams(documents, wl);
                documents = documents.getCountsTransposed;
                docsIn = 'columns';
            end

            if ~isnumeric(documents)
                textanalytics.internal.validateStringOrDocument(documents,2,nargin);
                if isa(documents,'tokenizedDocument') && ~iscolumn(documents)
                    error(message('textanalytics:ldaModel:TokenizedDocumentMustBeColumnVector'));
                end
                documents = textanalytics.internal.encodeNgramsWordlist(mdl.Wordlist,documents);
            end

            try
                documents = clustering.ldautils.prepareData(...
                    'predict','Prediction',documents,false,docsIn);
            catch ME
                throw(ME)
            end

            W = size(documents,1);

            Phi = mdl.PrivWordTopicProbabilities;
            if W~=size(Phi,1)
                error(message('textanalytics:ldaModel:PredictionDataNumWordsMismatch',size(Phi,1)));
            end

            validateattributes(maxiter,{'numeric'},...
                               {'scalar','real','>=',1,'integer'},'predict','IterationLimit');
            maxiter = double(maxiter);

            validateattributes(reltol,{'single','double'},...
                               {'scalar','real','>=',0},'predict','LogLikelihoodTolerance');
            reltol = double(reltol);

            alphas = mdl.CorpusTopicProbabilities'*mdl.TopicConcentration;

            Theta = clustering.ldautils.learn_theta_at_fixed_phi(...
                alphas,Phi,documents,reltol,maxiter);

            sumTheta = sum(Theta,1);
            Theta = Theta./sumTheta;
            bad = isnan(sumTheta);
            [~,cidx] = max(Theta,[],1);
            cidx(bad) = NaN;

            cidx = cidx';
            Theta = Theta';
        end


        function [logprob,ppl] = logp(mdl,documents,varargin)
            [varargin{:}] = convertStringsToChars(varargin{:});

            p = inputParser;
            p.CaseSensitive      = false;
            p.FunctionName       = 'logp';
            p.PartialMatching    = true;

            p.addRequired('Model');
            p.addRequired('InputData');

            p.addParameter('NumSamples',1000);
            p.addParameter('Stream',[]); % RandStream object
            p.addParameter('DocumentsIn','rows');

            p.parse(mdl,documents,varargin{:});

            nsample     = p.Results.NumSamples;
            rsh         = p.Results.Stream;
            docsIn      = p.Results.DocumentsIn;

            if isa(documents,'bagOfWords')
                if ~isscalar(documents)
                    error(message('textanalytics:ldaModel:NonscalarBagOfWords'));
                end
                documents = setVocabulary(documents, mdl.Wordlist);
                documents = documents.getCountsTransposed;
                docsIn = 'columns';
            end

            if isa(documents,'bagOfNgrams')
                if ~isscalar(documents)
                    error(message('textanalytics:ldaModel:NonscalarBagOfNgrams'));
                end
                [~,wl] = ...
                    textanalytics.internal.encodeNgramsWordlist(mdl.Wordlist,documents);
                documents = setNgrams(documents, wl);
                documents = documents.getCountsTransposed;
                docsIn = 'columns';
            end

            if ~isnumeric(documents)
                textanalytics.internal.validateStringOrDocument(documents,2,nargin);
                if isa(documents,'tokenizedDocument') && ~iscolumn(documents)
                    error(message('textanalytics:ldaModel:TokenizedDocumentMustBeColumnVector'));
                end
                [documents,~] = textanalytics.internal.encodeNgramsWordlist(mdl.Wordlist,documents);
            end

            try
                documents = clustering.ldautils.prepareData('logp','Prediction',documents,false,docsIn);
            catch ME
                throw(ME)
            end

            W = size(documents,1);

            Phi = mdl.PrivWordTopicProbabilities;

            if W~=size(Phi,1)
                error(message('textanalytics:ldaModel:PredictionDataNumWordsMismatch',size(Phi,1)));
            end

            validateattributes(nsample,{'numeric'},...
                               {'scalar','real','integer','>=',1},'logp','NumSamples');
            nsample = double(nsample);

            alphas = mdl.CorpusTopicProbabilities'*mdl.TopicConcentration;
            if sum(alphas)==0
                error(message('textanalytics:ldaModel:ZeroTopicConc'));
            end

            if isempty(rsh)
                rsh = RandStream.getGlobalStream;
            end

            % Pzw is a K-by-W matrix of topic probabilities per word
            Pzw = alphas.*Phi';
            Pzw = Pzw./sum(Pzw,1);

            logprob = clustering.ldautils.isip(alphas,Phi,documents,nsample,Pzw,rsh);

            if nargout > 1
                num_words = full(sum(documents(:)));
                ppl = exp(-nansum(logprob)/num_words);
            end
        end


        function T = topkwords(mdl,K,topic,varargin)
            [varargin{:}] = convertStringsToChars(varargin{:});

            p = inputParser;
            p.CaseSensitive      = false;
            p.FunctionName       = 'topkwords';
            p.PartialMatching    = true;

            p.addRequired('Model');
            p.addRequired('K');
            p.addRequired('TopicIndex');

            p.addParameter('Scaling','none');

            p.parse(mdl,K,topic,varargin{:});

            numWords = K;
            scaling  = p.Results.Scaling;

            numTopics = numel(mdl.CorpusTopicProbabilities);

            validateattributes(numWords,{'numeric'},...
                               {'scalar','real','>=',1},'topkwords','K',2);
            if numWords~=Inf
                validateattributes(numWords,{'numeric'},...
                                   {'integer'},'topkwords','K',2);
            end
            validateattributes(topic,{'numeric'},...
                               {'scalar','real','integer','>=',1,'<=',numTopics},'topkwords','TopicIndex',3);
            scaling = validatestring(scaling,{'none' 'inversemean'},'topkwords','Scaling');

            dict = mdl.Wordlist;
            numWords = min(numWords,dict.NumWords);

            [scores,ranked] = wordScores(mdl,scaling);

            idx = ranked(1:numWords,topic);

            dict_idx = dict(idx);
            dict_idx = dict_idx(:);

            T = table(dict_idx, scores(idx,topic), ...
                      'VariableNames',{'Word' 'Score'});
        end

        function h = wordcloud(mdl,topic,varargin)
            [varargin{:}] = convertStringsToChars(varargin{:});

            T = topkwords(mdl,Inf,topic);
            h = wordcloud(T.Word,T.Score,varargin{:});
        end
    end
end


function [alphas,Phi,Theta,fitted_Nwk,fitted_Nkj] = ...
        zeroOutBadTopics(alphas,Phi,Theta,fitted_Nwk,fitted_Nkj)
    bad = any(isnan(Phi),1)' | alphas==0 | isnan(alphas);
    if any(bad)
        warning(message('textanalytics:ldaModel:TopicsWithZeroProbs',sum(bad)));
        alphas(bad)       = 0;
        Phi(:,bad)        = 0;
        Theta(bad,:)      = 0;
        fitted_Nkj(bad,:) = 0;
        fitted_Nwk(:,bad) = 0;
    end
end
