function [log_pdf_val, hparm] = hmm_reest(hparm, event_data, istart, nsamp, NIT, ...
    data_wts, iplot);
% function [log_pdf_val, hparm] = hmm_reest(hparm, event_data, istart, nsamp, NIT, ...
%   [data_wts], [iplot]);
%
% This program updates the HMM parameters (re-estimates them) using Baum-Welch.
%
% INPUT:
%   hparm:
%   event_data:  This is expected to be a matrix (nFeat x nSamples).
%      However, if this is actually a feature data structure, then 
%      we use a hmm_accum_big() instead of hmm_accum().
%   istart:
%   nsamp:
%   NIT: number of iterations, zero for an evaluation only.  In the zero
%     case, data_wts is ignored.
%   data_wts: default to []
%   iplot:
%   crp_data:  If this optional argument is present, then mrhmm_accum()
%     is used instead of hmm_accum(). Also, istart and dsamp are ignored.
%     This is for class-specific HMM (allows a different feature set for
%     each state).
%
% OUTPUT:
%   log_pdf_val:
%   hparm:
%

%
% 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.
%


% Algorithm description:
% parameters governing training:
%    pout_increase_per_thousand_threshold and max_n_pout_no_increase work
%    together to define a stopping condition.  If max_n_pout_no_increase
%    iterations each result in an increase of pOut that is less than
%    pout_increase_per_thousand_threshold times the number of thousands of
%    samples (columns) in event_data, then training is aborted.
%    But, the count of "weak" iterations is halved whenever the number of
%    modes change during an iteration.
%
%    var_fac: used for expanding the variance as a type of annealing.
%      The input var_fac is a kind of annealing term which is normally 1,
%      but if higher will expand the variances of all the kernels by that factor.

global data_directory;

feat_info=[];
[dim, ntot] = size(event_data);

nevent = length(istart);

if(nargin < 6 || isempty(data_wts) ), data_wts=ones(nevent,1); end;
if(nargin < 7 || isempty(iplot))  , iplot=0; end;

if(iplot== -1), be_quiet=1; iplot=0; else; be_quiet=0; end;


% if older version of parameter without "powr", set field
if(~isfield(hparm,'powr')), hparm.powr = 2; end;
% The input "powr" is the power of the exponent in the PDF, normally
% equal to 2 (for Gaussian).  Higher numbers make the kernels roll off faster,
% with powr=infinity producing a flat elliptically-shaped table kernel.
% Smaller numbers produce heavier tails. Powr=1 is exponential kernel.


if(length(istart) ~= length(nsamp) || sum(nsamp) ~= ntot),
  error('length(istart) ~= length(nsamp) || sum(nsamp) ~= ntot');
end;

DEBUG = 0;

if NIT <= 0	% NIT = 0, eval only
    % Accumulate the record
    
    if(DEBUG),
        save tmp.mat  hparm  event_data  istart  nsamp
    end;
    pout = hmm_accum(hparm, event_data, istart, nsamp);
    
    % Write a dump file if there is a NaN
    if (any(isnan(pout(:)))),
        fprintf('bad pout value\n');
        dumpCrashReport(hparm, event_data, istart, nsamp, crp_data, feat_info);
    end
    log_pdf_val = pout;
    return
end

% do training
if(isfield(hparm,'bias')), BIAS = hparm.bias;  else,  BIAS=1; end;
if(isfield(hparm,'idiag')),   DIAG = hparm.idiag;  else,    DIAG=0; end;
plast = -inf;
N = hparm.N;

pout_increase_per_thousand_threshold = .1;
n_pout_no_increase = 0;
max_n_pout_no_increase = 10;
var_fac = 1.0;
Npdf=N;

% determine current n_modes, which may be less than originally selected
for ipdf=1:Npdf,
    current_n_modes_a(ipdf) = length(hparm.pdf(ipdf).modes);
end
state_wts = ones(N,1) / N;
for it = 1:NIT,
    
    % zero accumulators
    A_hat = zeros(N,N);
    Pi_hat = zeros(N,1);

    for istate = 1:N
            % nmodes can change during iteration, but
            % current_n_modes_a(istate) is kept accurate, so it could be used
            % here.
            nmode = length(hparm.pdf(istate).modes);
            newmean{istate} = zeros(dim,nmode);
            newvar{istate} = zeros(dim*dim,nmode);
            atot{istate} = zeros(nmode,1);
    end;

    % Perform a step of EM (baumwelch algo)
    if(DEBUG),
        fprintf('HMM_REEST: Debug mode: saving tmp.mat\n');
        save tmp.mat atot A_hat Pi_hat data_wts hparm newmean newvar
    end;
    [pout, newmean, newvar, atot, A_hat, Pi_hat] = ...
        hmm_accum(hparm, event_data, istart, nsamp, newmean, newvar, ...
        atot, A_hat, Pi_hat, data_wts, hparm.powr, var_fac);
    if(DEBUG),
        % Store Q-function so we can later check for increase
        q1=sum(pout(:) .* data_wts(:));
        %for kk=1:3, fprintf('HMMREEST1: %f %f\n',pout(kk), data_wts(kk)); end;
    end;
    
    % check for errors
    ib=0;
    for istate=1:N, ib = ib + any(isnan(atot{istate})); end;
    if(ib || any(isnan(pout(:)))),
        fprintf('hmm_reest: ib || isnan(pout)\n'); 
        keyboard
    end;
   
    newvar = fixWeirdNumericalErrors(newvar);
    
    % determine relative weight of each state
    for istate=1:N,
        state_wts(istate) = sum(atot{istate});
    end;
    if any(isnan(state_wts))'
        'nan in state_wts'
         keyboard
    end
    state_wts = state_wts/sum(state_wts);

    % calculate the effective number of samples assigned to each state
    %for istate=1:N,
    %       neff(istate) = state_wts(istate)*ntot;
    %       % static states inherit m
    %       if(isfield(hparm,'static') & hparm.static(istate) & istate<N & hparm.static(istate+1)==0),
    %            neff(istate) = neff(istate) + state_wts(istate+1)*ntot;
    %       end;
    %end;
        
     min_state_wt = 1e-10;
    
    %TODO: something like: hparm = hmm_prune(hparm, min_state_wt);
    % Prune very small states
    if(any(state_wts < min_state_wt)),
        fprintf('Warning: pruning some states\n');
        igood = find(state_wts >= min_state_wt);
        
        A_hat = A_hat(igood, igood);
        Pi_hat = Pi_hat(igood);
        atot = atot(igood);
        current_n_modes_a = current_n_modes_a(igood);
        newmean = newmean(igood);
        newvar = newvar(igood);
        state_wts = state_wts(igood);
        N = length(igood);

        if(isfield(hparm,'pdf_to_featureset_index')),
        else,
              Npdf=N;
        end;
        
        hparm.pdf = hparm.pdf(igood);
        hparm.A = hparm.A(igood,igood); % will be overwritten
        hparm.Pi = hparm.Pi(igood); % will be overwritten
        hparm.N = N;
        if(isfield(hparm,'mask')),
            hparm.mask = hparm.mask(igood,igood);
        end;
        if(isfield(hparm,'pi_mask')),
            hparm.pi_mask = hparm.pi_mask(igood);
        end;
    end;
    
    % Normalize A_hat and Pi_hat
    if(isfield(hparm,'mask')),
        A_hat = A_hat .* hparm.mask;
    end;
    asum=sum(A_hat, 2);  % Sum the rows
    asum(find(asum==0)) = 1;
    
    hparm.A = A_hat ./ repmat(asum,1,N);
    if(isfield(hparm,'pi_mask')),
        Pi_hat = Pi_hat .* hparm.pi_mask;
    end;
    hparm.Pi = Pi_hat / sum(Pi_hat);

    % Update Gaussian mixtures
    for ipdf=1:Npdf,
        gparm = hparm.pdf(ipdf);
        gparm = gmix_norm(gparm, newmean{ipdf}, ...
            newvar{ipdf}, atot{ipdf}, BIAS,DIAG);
        
        % Prune the very very small weights
        gparm = gmix_deflate(gparm, 1e-5, 1e-5);

        % determine effective number of (segment) samples for this state
        numEffSamples = state_wts(ipdf) *  ntot;
        [gparm, didPrune] = gmix_prune(gparm, numEffSamples);
        if didPrune  
            fprintf('hmm: pruned state %d down to %d modes.\n',...
                ipdf, length(gparm.modes));
        end
        hparm.pdf(ipdf) = gparm;
    end;

    % check whether any of the n_modes changed
    n_modes_changed = 0;
    for ipdf = 1:Npdf,
        if current_n_modes_a(ipdf) ~= length(hparm.pdf(ipdf).modes) 
            n_modes_changed = 1;
            current_n_modes_a(ipdf) = length(hparm.pdf(ipdf).modes);
        end;
    end;
    if( n_modes_changed ),
        n_pout_no_increase = ceil(n_pout_no_increase / 2);  % maybe this should be a parameter
    end;
    
    sumPout = sum(pout(:) .* data_wts(:));
    if(it > 1),
        if(~be_quiet),
            fprintf('Iteration %g, Pout=%g: delta=%g', it, sumPout, sumPout-plast);
         end;
    else,
        %fprintf('Iteration %g, Pout=%g', it, sumPout);
    end;
    if(DEBUG),
        pout2 = hmm_accum(hparm, event_data, istart, nsamp);
        q2=sum(pout2(:) .* data_wts(:));
        %for kk=1:3, fprintf('HMMREEST2: %f %f\n',pout2(kk), data_wts(kk)); end;
        if(q2 < q1), 
            fprintf('HMMREEST: Q function decrease!\n');
            keyboard
        else,
            fprintf('HMMREEST: Q function increase %f\n',q2-q1);
        end;
    end;

    if( sumPout - plast < pout_increase_per_thousand_threshold*(ntot/1000)),
        n_pout_no_increase = n_pout_no_increase + 1;
        
        if( n_pout_no_increase >= max_n_pout_no_increase ),
            if(~be_quiet),
                fprintf('  n_no_incr=%d\n',  n_pout_no_increase);
            end;
            break;  % breaks out of for it=1:NIT loop
        end;
    end;

    plast = sumPout;
    if(it > 1 && ~be_quiet),
      fprintf('  n_no_incr=%d\n',  n_pout_no_increase);
    end;
    
end; % it

log_pdf_val = pout;


