function [pout,newmean,newvar,atot,A_hat,Pi_hat] = hmm_accum(hparm, ...
    event_data, starts, nsamp, newmean, newvar, atot, ...
    A_hat, Pi_hat, rec_wts, powr, var_fac);
%function [pout,newmean,newvar,atot,A_hat,Pi_hat] = hmm_accum(hparm, ...
%    event_data, starts, nsamp, newmean, newvar, atot, ...
%    A_hat, Pi_hat, rec_wts, powr, var_fac);
%
% 
% This function is used by hmm_reest to accumulate the contribution of 
% one data record.
% CAS: It looks like it actually accumulates the contribution of all data
% records in the parallel vectors, istart and nsamp.
% 
%
% data_wts is optional and defaults to uniform weighting.
% powr is optional and defaults to 2 (Gaussian). HOWEVER: powr IS UNUSED!
% var_fac is optional and defaults to 1 (no variance expansion).
%
%

%
% NAVAL UNDERSEA WARFARE CENTER
% class-specific@npt.nuwc.navy.mil
% http://www.npt.nuwc.navy.mil/csf
% Paul M. Baggenstoss
% 401-832-8240
% Do not disseminate.
%


% NOTE: This function has a corresponding mex implementation.

% CAS: mex code outputs pout as a column vector, now we do the same here.
    
nrec = length(starts);
[dim, ntot] = size(event_data);   
if(nargin <= 11 || isempty(var_fac) ) , var_fac=1; end;
if(nargin <= 10 || isempty(powr) ), powr=2; end;
if(nargin <= 9 || isempty(rec_wts) ), rec_wts = ones(nrec,1); end;

if(nargin <= 4)
    % zero accumulators
    N = hparm.N;
    A_hat = zeros(N,N);
    Pi_hat = zeros(N,1);
    for i = 1:hparm.N
        nmode = length(hparm.pdf(i).modes);
        newmean{i} = zeros(dim,nmode);
        newvar{i} = zeros(dim*dim,nmode);
        atot{i} = zeros(nmode,1);
    end
end;

if(length(nsamp) ~= nrec), error('mismatch nsamp vs. starts'); end;
if(ntot ~= sum(nsamp)), error('ntot ~= sum(nsamp)'); end;
if(starts(1) ~= 1), error('starts(1) ~= 1'); end;

pout=zeros(1,nrec);
for irec=1:nrec,
       idx = starts(irec) + [0:nsamp(irec)-1];
       x = event_data(:,idx);
       [newmean,newvar,atot,A_hat,Pi_hat,pout2,gamma,l] = ...
                hmm_accum_rec(hparm, x, newmean, newvar, atot, A_hat, Pi_hat, pout(irec),rec_wts(irec));
       pout(irec)=pout2;
end;
pout = pout(:);
return


function   [newmean,newvar,atot,A_hat,Pi_hat,pout,gamma,l] = ...
    hmm_accum_rec(hparm,x,newmean, newvar,atot,A_hat,Pi_hat,pout,wt);

BAD = -9e37;
N=hparm.N;
[dim,nsamp]=size(x);

[l, pk] = hmm_get_probs(hparm, x);

lmax=max(l,[],2);
igood=find(lmax > BAD/2);
gamma=l;

bout=ones(N,1);
[alphas,alognorm,betas,blognorm,tgamma,A_hat] = ...
    baumwelch(l(igood,:),lmax(igood),hparm.A,hparm.Pi,A_hat,bout);
gamma(igood,:)=tgamma;

Pi_hat = Pi_hat + gamma(1,:)' * wt;
pout = pout + (log(sum(alphas(end,:)))+alognorm(end))*wt;

for istate=1:N,
    gparm = hparm.pdf(istate);
    nmode = length(gparm.modes);
    wts = [gparm.modes.weight]';
    datawts = pk{istate} .* repmat(gamma(:,istate),1,nmode) *wt;
    for imode=1:nmode,
        atot{istate}(imode)= atot{istate}(imode)+sum(datawts(:,imode));
        newmean{istate}(:,imode) =  newmean{istate}(:,imode) ...
            + (datawts(:,imode)'*x')';
        y = x' - repmat(gparm.modes(imode).mean',nsamp,1);
        y = y .* repmat(sqrt( datawts(:,imode)),1,dim);
        tmp = y'*y;
        newvar{istate}(:,imode) =  newvar{istate}(:,imode) + tmp(:);
    end;
end;
return
