function [logl, cls] = mseqvit(cbr, furl, lurl, TRANS, alpha, beta)

% Check input arguments
narginchk(6, 6);

% Data length
Td = cat(1,cbr{:,1});
ccnt = cbr(:,3);
cnorm = cbr(:,4);

% Number of mixture components
Q = nbh_chk(TRANS, alpha, beta);
alpha = reshape(alpha, 1, Q);
beta  = reshape(beta, 1, Q);

% HMM transition probabilities in log
TRANS = log(TRANS + realmin);

% ldens = cell(lurl,1);
ccls =  cell(lurl,1);
logl = zeros(lurl,1);

%%% 0: Compute density values on a log scale
parfor d = 1:lurl
    [logl(d), ccls{d}] = dpeq(clfemis(Td(d), ccnt{d}, cnorm{d}, furl(d), Q, alpha, beta), Q, TRANS, Td(d), furl(d));
end

logl = sum(logl);
cls = cat(1,ccls{:});
% Dynamic programming recursion

end

function ldens = clfemis(Td, count, dnorm, nreg, Q, alpha, beta )
% evalutate log emission probabilites
% to be used as a cell function.
T = Td * nreg;

ldens = ones(T,1)*(alpha.*log(beta./(1+beta)) - gammaln(alpha)) ...
    - count*log(1+beta) + gammaln(count*ones(1,Q) + ones(T,1)*alpha) ...
    - dnorm*ones(1,Q);

end

function [logl, cls] = dpeq( ldens, Q, TRANS, T, nreg )
% Dynamic programming recursion

nel = T*nreg;

% Partial loglikelihood array and bactracking array
PLOGL = zeros(nel, Q);
BCKTR = zeros(nel, Q, 'int32');

if nreg > 1
    bdx = (0:(nreg-1)).'*T;
    mdx = repmat(1:Q,nreg,1);
    mdx = mdx(:);
    
    t = 1;
    idx = bdx + t;
    PLOGL(idx,:) = ldens(idx,:) + ones(nreg,1) * TRANS(1,:);
    
    for t = 2:T
        idx = bdx + t;
        [temp1, temp2] = max(repmat(PLOGL(idx-1,:)',[1,1,Q]) + reshape(TRANS(:,mdx),Q,nreg,Q));
        PLOGL(idx,:) = squeeze(temp1) + ldens(idx,:);
        BCKTR(idx-1,:) = squeeze(temp2);
    end
    
    PLOGL(idx,:) = ones(nreg,1) * TRANS(:,1)' + PLOGL(idx,:);
    
    % Backtracking
    cls = zeros(nel,1,'int32');
    [logl, cls(idx)] = max(PLOGL(idx,:),[],2);
    
    for t = (T-1):-1:1
        idx = bdx + t;
        cls(idx) = BCKTR(sub2ind([nel Q], idx, cls(idx+1)));
    end
else
    PLOGL(1,:) = ldens(1,:) + TRANS(1,:);
    for t = 2:T
        [PLOGL(t,:), BCKTR(t-1,:)] = max((PLOGL(t-1,:)' * ones(1,Q)) + TRANS);
        PLOGL(t,:) = PLOGL(t,:) + ldens(t,:);
    end
    
    t = T;
    [PLOGL(t,:), BCKTR(t-1,:)] = max((PLOGL(t-1,:)' * ones(1,Q)) + TRANS);
    PLOGL(t,:) = TRANS(:,1)' + PLOGL(t,:) + ldens(t,:);
    
    % Backtracking
    cls = zeros(T,1,'int32');
    [logl, cls(T)] = max(PLOGL(T,:));
    for t = (T-1):-1:1
        cls(t) = BCKTR(t, cls(t+1));
    end
end

logl = sum(logl);

end


