function [ results, w, alp, ep, lobs ] = findbinding( mfreq,tfreq,mis,zs,dat,agg,cls,nhs,w,alp,Nit )
% the preparation step to implement the EM alorithm for the mixture model.
% the data is sliced to implement parallel compution efficiently.

mrlen = length(tfreq);

azs = agg.zs;
ard = agg.rd;
slc = agg.slc;

cls = (cls == nhs);

dzs = dat.zs;
slo = dat.slc(~dzs);

enidx = (ard ~= azs & cls);

enc = slc(enidx);
% eno = 1 if a position in a read mapped to highly enriched location
[eno, edx] = ismember(slo,enc); 
edx = int32(edx(eno)); % highly-enriched genomic locations ``edx''
[~, edx] = histc(edx, unique(edx));
edx = int32(edx); % read counts in highly-enriched genomic locations 
dmut = dat.mut(~dzs); % mutation or no-mutation
dpos = int32(dat.pos(~dzs)); % postions in reads
drlen = int32(dat.rlen(~dzs)); % lengths
dstrn = int32(dat.strn(~dzs) + 1); % stands

mis.mfreq{1} = double(mis.mfreq{1});
mis.mfreq{2} = double(mis.mfreq{2});
mis.tfreq = double(mis.tfreq);

mfreq{1} = double(mfreq{1});
mfreq{2} = double(mfreq{2});
tfreq = double(tfreq);

tnzs = sum(sum(cellfun(@(x) sum(~x), zs)));

ub = inf;
lb = -inf;

% INITIALIZE the heterogeneous effects
ep = zeros(mrlen,mrlen,2);
for k = 1:2
    ep(:,:,k) = mfreq{k} ./ (tfreq(:,k) * ones(1,mrlen));
end

ep = logit(ep);
ep(ep > ub) = ub;
ep(ep < lb) = lb;
ep(isnan(ep)) = -inf;

sudx = maptdx(tnzs,zs,mrlen); % an one-to-one map: theta -> ep
ldx = sub2ind(size(ep), sudx(2:end,1), sudx(2:end,2), sudx(2:end,3));
cidx = sub2ind(size(ep), drlen, dpos, dstrn);

epdx = cell(tnzs,1);
parfor k = 1:tnzs
    epdx{k} = int32(find(drlen == sudx(k+1,1) & dpos == sudx(k+1,2) & dstrn == sudx(k+1,3)));
end

adx = cidx; % SUB to be used for accumarray
for k = 1:tnzs
    adx(epdx{k}) = k;
end

% input variables
muput = {eno, edx, dmut, mrlen, mis, tnzs, cidx, adx, ldx};

% the EM for the mixture model
[ w, as, alp, ep, lobs ] = ...
      bmem( muput, w, alp, ep, Nit );

% to generate the matrix for results
[sas, sdx] = sort(as,'descend');

results = {sas, agg.chr(sdx),agg.reg(sdx),agg.loc(sdx),agg.rd(sdx),agg.mut(sdx),agg.zs(sdx),agg.strn(sdx)};

end

% function to obtain the matrix of the heterogeneous effets ``ep'' from theta
function tdx = maptdx( tnzs, zs, mrlen )
% an one-to-one map: theta -> ep
tdx = zeros(tnzs + 1,3,'int32');
k = 1;
for n = 1:2
    for i = 1:mrlen
        if any(~zs{n,i})
            jdx = 1:i;
            jdx = jdx(~zs{n,i});
            for j = jdx
                k = k + 1;
                tdx(k,:) = [i, j, n];
            end
        end
    end
end

end