function [ w, as, alp, ep, lobs ] = ...
    bmem( muput, w, alp, ep, Nit )
%UNTITLED4 Summary of this function goes here
%   Detailed explanation goes here

% muput = {eno, edx, mut, mrlen, mis, tnzs, cidx, adx, ldx, smtd};

eno = muput{1};
edx = muput{2};
mut = muput{3};
mrlen = muput{4};
mis = muput{5};
tnzs = muput{6};
cidx = muput{7};
adx = muput{8};
ldx = muput{9};

s = zeros(length(eno),1); % class membership.

theta = zeros(tnzs + 1,1);
thetac = ones(tnzs + 1,1); % theta^c
thetamat = zeros(tnzs + 1, Nit );

mscr = zeros(mrlen,mrlen,2);
mhes = zeros(mrlen,mrlen,2);

lobs = zeros(Nit,1);

sn = 0.001; %stop if |lobs - lobsc| < sn
mnri = 10; %maximum iteration within the M step

mismv = cat(3,mis.mfreq{:});
mistv = zeros(mrlen,mrlen,2);
mistv(:,:,1) = mis.tfreq(:,1) * ones(1,mrlen);
mistv(:,:,2) = mis.tfreq(:,2) * ones(1,mrlen);


for nit = 1:Nit
    %EXPECTATION STEP
    % compute density of each position for each mixture hands
    loe  =  likhdobs(mut(eno), alp, ep, cide, true );
    lou  =  likhdobs(mut(~eno), alp, ep, cidu, false );
    %         wpm  = w*mexplike(edx, loe.m ,tec, teo);
    %         wpe  = (1-w)*mexplike(ft .idx, loe,tec, teo);
    
    % compute density of each location for each mixture hands
    wpm = w*accumarray(edx,loe.m,[], @prod);
    wpe = (1-w)*accumarray(edx,loe.e,[], @prod);
    
    % estimated posterior probability
    as  = wpm ./(wpm  + wpe );
    
    s(eno) = as(edx);
    
    mpe = ilogit(ep);
    temp = mismv(ldx) .* ep(ldx) + log(1-mpe(ldx)) .* mistv(ldx);
    lobs(nit) = sum(log(lou.e)) + sum(log(wpm  + wpe )) + sum(temp);
    
    % stopping rule
    thetamat(1,nit) = alp;
    thetamat(2:tnzs+1,nit) = ep(ldx);
    if nit > 1
        if norm(lobs(nit) - lobs(nit-1)) < sn
            fprintf('Estimates coverge to local maximums.')
            break;
        end
    end
    
    fprintf('# iteration: %d, L_obs = %f \n', nit, lobs(nit));
    
    %MAXIMIZATION STEP
    % Update w (mixture proportion)
    w = mean(as);
    
    % mixture proportion parts in the complete likelihood evaluation.
    %     wlc = as;
    %     wlc = sum(wlc)*log(w) + sum(1-wlc)*log(1-w);
    
    nri = 0; % Newton Raphson iteration
    while (norm(theta - thetac) > 0.1 || nri == 0)
        nri = nri + 1;
        if nri > mnri
            error('Too long time to increase the Q function.');
        end
        
        p.m = ilogit(alp + ep(cidx)); % true mutation probs.
        p.e = ilogit(ep(cidx)); % error mutation probs.
        
        %socre function: d E(L_com|data) / d theta
        scr = zeros(tnzs + 1,1);
        scrv = [s.*(mut - p.m), (1-s).*(mut - p.e)];
        scr(1) = sum(scrv(:,1)); % d E(L_com|data) / d alpha
        
        hes = zeros(tnzs + 1);
        hesv = -[s.*(1 - p.m).*p.m, (1-s).*(1 - p.e).*p.e];
        
        hes(1,1) = sum(hesv(:,1)); % d^2 E(L_com|data) / d alpha^2
        
        % used non-highly enriched locations to construct socre / hessian
        for k = 1:2
            mscr(:,:,k) = mis.mfreq{k} - (mis.tfreq(:,k) * ones(1,mrlen)) .* ilogit(ep(:,:,k));
            mhes(:,:,k) = - (mis.tfreq(:,k) * ones(1,mrlen)) .* (1 - ilogit(ep(:,:,k))).*ilogit(ep(:,:,k));
        end
        
        % faster summations of score function / hessian
        scr(2:(tnzs+1)) = accumarray(adx,sum(scrv,2)) + mscr(ldx);
        shesv.a = accumarray(adx,hesv(:,1));
        shesv.e = accumarray(adx,hesv(:,2));
        hes(1,2:(tnzs+1)) = shesv.a;
        hes(2:(tnzs+1),2:(tnzs+1)) = diag(shesv.a + shesv.e +  mhes(ldx));
        hes(:,1) = hes(1,:);
        
        % reshape (alp, ep) to obatain theta
        thetac(1) = alp;
        thetac(2:tnzs+1) = ep(ldx);
        
        % Newton-Raphson update for theta
        theta = thetac - linsolve(hes,scr); %theta = theta_c - Hessian * Score
        
        % reshape theta to obatain (alp, ep)
        alp = theta(1);
        ep(ldx) = theta(2:tnzs+1);
    end
end

end

function lobs = likhdobs(m, alp, ep, idx, lobsm )
% Observed Likelihood function on each Mixture Hands
% No product has been taken
% m is logical (mutation or not)
% lobs.m = likelihoods for true mutations
% lobs.e = likelihoods for error mutations

if lobsm
    lobs.m = (ilogit(alp + ep(idx)).^m).*(1-ilogit(alp + ep(idx))).^(1-m);
end
lobs.e = (ilogit(ep(idx)).^m).*(1-ilogit(ep(idx))).^(1-m);

end