function parm = gmix_norm(parm,newmean,newvar,atot, bias,diag_cov,usedmean);
%function parm = gmix_norm(parm,newmean,newvar,atot, [bias=0],[diag_cov=0],[usedmean=1]);
%
% Subrouting to complete Gaussian mixture update
% (see gmix_accum.m)
%
% Inputs:
%        parm     Input/Output parameters, 
%        newmean  accumulated mean estimates
%        newvar   accumulated covariance estimates
%        atot     total weights used to normalize the above
%        bias     (optional) bias=1 for BIAS method,
%                   0 for CONSTRAINT method (default)
%        usedmean   set to 1 if mean was subtracted during variance accumulation
%                 default: usedmean=1;
%
%  Dr. Paul M. Baggenstoss
%  Naval Undersea Warfare Center
%  Newport, RI
%  p.m.baggenstoss@ieee.org
%  Date              Reason
%  ------------      --------------
%  Nov 21, 1999      Initial release
%  Dec 7,  1999      MATLAB 5 Upgrade
%  Mar 3,  2000      Fixed bug - not normalizing min_std 
%                     and use eig() in place of chol()
%  Mar 13,  2000     Fixed bug - could divide by 0 if atot(k)=0;
%  Feb 16,  2009     added diagonal covariance matrix support
%  Apr 16,  2009     added usedmean argument
%  May  1,  2009     added FRACTION method 
%  Nov 22, 2017      replaced code segment to condition covariance
%                             with call to condx_chol
%
% 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.
%


global data_directory;

nmode = length(parm.modes);
DIM = length(parm.features);
min_std = [parm.features.min_std]';

if(nargin < 5), bias=0; end;
if(nargin < 6), diag_cov=0; end;
if(nargin < 7), usedmean=1; end;

assert(all(size(newmean) == [DIM,nmode]));
assert(all(size(newvar) == [DIM*DIM,nmode]));
assert(all(size(atot) == [nmode,1]));

if(sum(atot)==0), 
  fprintf('gmix_norm: warning sum(atot)==0\n'); 
  return; 
  end;

if(diag_cov),
  tmpvar=zeros(DIM,1);
else,
  tmpvar=zeros(DIM,DIM);
end;
for k=1:nmode,
    
    if (atot(k) > 0) ,
        % correct variance estimate - needed because 
        % 'old' mean was used in 'newvar' accumulator
        mu = newmean(:,k)/atot(k);  
        % mu0 is the mean that was used during variance accumulation
        mu0 = parm.modes(k).mean * usedmean;
        parm.modes(k).mean=mu;
        if(diag_cov),
           tv=newvar(:,k);
           tmpvar(:)=diag(reshape(tv,DIM,DIM))/atot(k);
        else,
           tmpvar(:)=newvar(:,k)/atot(k);
        end;
        if(diag_cov),
           tmpvar  = tmpvar  - (mu0-mu).^2;
           std_est = sqrt(tmpvar);
        else,
           tmpvar  = tmpvar  - (mu0-mu) * (mu0-mu)';
           std_est = sqrt(diag(tmpvar));
        end;

        % implement FRACTION method
        % each negative min_std entry is interpreted as a fraction of std_est
        min_std2 = min_std;
        ii=find(min_std<0);
        min_std2(ii) = (-min_std(ii)) .* std_est(ii);

        if(diag_cov),
            switch(bias),
                case 0, % CONSTRAINT method
                    tmpvar=max(tmpvar,min_std2(:).^2);
                case 1, % BIAS
                    tmpvar=tmpvar+min_std2(:).^2;
            end;
            Rout=diag(sqrt(tmpvar));
        else,
            
            % tmpvar=condx_chol(tmpvar,0,bias,min_std);
            tmpvar=condx_cov(tmpvar,bias,min_std2);
            
            Rout = chol(tmpvar);
        end;
        parm.modes(k).cholesky_covar = Rout;
    end;
    
    wts(k) = atot(k)/sum(atot);
    %CAS: what's this?  shouldn't we complain up front for negative weights?
    parm.modes(k).weight = wts(k) * (wts(k) > 0); 
end;
