function [ cls, TRANS, alpha, beta, loglvit] = stateinference( count, areg, ureg, mis, TRANS, alpha, beta, Nit)
%UNTITLED5 Summary of this function goes here
%   Detailed explanation goes here

sdx = int32(1:length(areg))';
rl = int32(histc(areg,ureg)); % region lengths
url = unique(rl); % unique region lengths
lurl = length(url); % number of different region lengths
furl = int32(histc(rl,url)); % frequencies of url.
count = double(count - 1);

% Compute log(count!)
tmp = cumsum([0; log((1:max(count)).')]);
dnorm = tmp(count+1);

% group regions by the region length to faciliate parallel computing
cbr = cell(lurl,5);
parfor i = 1:lurl
    cbr(i,:) = byreglength( count, dnorm, sdx, areg, ureg, rl, url(i) );
end

% number of locations covered by non-clustered reads.
misrd = sum(mis.tfreq,2);

% The Baum-Welch EM algorithm
[TRANS, alpha, beta, ~, ~] = mseqnbh(cbr(:,1:4), furl, lurl, misrd, TRANS, alpha, beta, Nit);

% The Viterbi algorithm
[loglvit, cls] = mseqvit(cbr(:,1:4), furl, lurl, TRANS, alpha, beta);

% to match locations with maximum a posteriori path. 
sdx = cat(1,cbr{:,5});
[~, sdx] = sort(sdx);
cls = cls(sdx,:);
end

% a subfunction for group regions by the region length
function crray = byreglength( count, dnorm, sdx, dreg, ureg, rl, url )
	idx = ismember(dreg,ureg(rl == url));
    brlc = count(idx);
    brld = dnorm(idx);
    brdx = sdx(idx);
    brlr = ureg(rl == url);
    crray = {url, brlr, brlc, brld, brdx};
end