%MEANCOV Estimation of the means and covariances from multiclass data
% 
%   [U,G] = meancov_new(A,N)
% 
%	INPUT
% 	A		Dataset
%   N		Normalization to use for calculating covariances: by M, the number
%					of samples in A (N = 1) or by M-1 (default, unbiased, N = 0).
%
% OUTPUT
% 	U		Mean vectors
% 	G		Covariance matrices
%
% DESCRIPTION  
%	Computation of a set of mean vectors U and a set of covariance matrices G
%	of the C classes in the dataset A. The covariance matrices are stored as a
%	3-dimensional matrix G of the size K x K x C, the class mean vectors as a
%	labeled dataset U of the size C x K.
%
% The use of soft labels or target labels is supported.
% 
%	SEE ALSO 
%	DATASETS, GAUSS, NBAYESC, DISTMAHA

% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
% Faculty of Applied Sciences, Delft University of Technology
% P.O. Box 5046, 2600 GA Delft, The Netherlands

% $Id: meancov.m,v 1.17 2005/10/19 08:32:11 duin Exp $

function [U,G] = meancov_new(a,n)

	prtrace(mfilename);

	% N determines whether the covariances are normalized by M (N = 1) or by 
	% M-1 (unbiased, N = 0), where M is the number of objects.

	if (nargin < 2)
		prwarning(4,'normalisation not specified, assuming by M-1');
		n = 0;
	end

	if (n ~= 1) && (n ~= 0)
		error('Second parameter should be either 0 or 1.')
    end

	if (~isa(a,'prdataset'))			% A is a matrix: compute mean and covariances
        
        DirectionalFeatures = [];
        
        if( isempty(DirectionalFeatures) )
            U = nanmean(a);
            G = nancov(a,n);
        else
            
            iCantFeatures = size(a,2);
            U = nan(1, iCantFeatures);
            NotDirectionalFeatures = setdiff(1:iCantFeatures, DirectionalFeatures);

            if( ~isempty(NotDirectionalFeatures) )
                U(NotDirectionalFeatures) = nanmean(a(:,NotDirectionalFeatures));							% 	in the usual way.
            end
            
            U(DirectionalFeatures) = angle(nanmean(exp(1i*a(:,DirectionalFeatures))));
            
            bNotNaN = ~any(isnan( a(:, DirectionalFeatures ) ),2);
            iCantElementos = sum(bNotNaN);
            
            if(iCantElementos < 2)
                
                G = zeros(iCantFeatures);
                
            else
                iAux = nan(iCantElementos, iCantFeatures);

                if( ~isempty(NotDirectionalFeatures) )
                    iAux(:,NotDirectionalFeatures) = a(bNotNaN,NotDirectionalFeatures) - repmat(U(NotDirectionalFeatures), iCantElementos,1);
                end

                for jj = DirectionalFeatures                                     
                    iAux1 = a(bNotNaN, jj );
                    iAux2 = repmat(U(jj), iCantElementos,1);
                    iAux3 = [ iAux1 - iAux2, 2*pi + iAux1 - iAux2, iAux1 - 2*pi - iAux2 ];
                    [dummy, iMinIndex] = min(abs(iAux3), [], 2);
                    iMinIndex = sub2ind(size(iAux3),1:size(iAux3,1),iMinIndex');
                    iAux(:,jj) = iAux3(iMinIndex);
                end 

                if(n == 1)
                    G = 1/iCantElementos*(iAux'*iAux);
                else
                    G = 1/(iCantElementos-1)*(iAux'*iAux);
                end
            end

        end
        
    else
        
		[m,k,c] = getsize(a);

        cFeaturesDomain = getfeatdom(a);
        
        DirectionalFeatures = [];
        
        for jj = 1:length(cFeaturesDomain)
            if(~isempty(cFeaturesDomain{jj}))
                DirectionalFeatures = [DirectionalFeatures jj];
            end
        end

        if (islabtype(a,'crisp'))
			
			if (c==0)
                
                aa = +a;
                
                if( isempty(DirectionalFeatures) )
                    U = nanmean(aa);
                    G = nancov(aa,n);
                else

                    iCantFeatures = k;
                    U = nan(1, iCantFeatures);
                    NotDirectionalFeatures = setdiff(1:iCantFeatures, DirectionalFeatures);

                    if( ~isempty(NotDirectionalFeatures) )
                        U(NotDirectionalFeatures) = nanmean(aa(:,NotDirectionalFeatures));							% 	in the usual way.
                    end
                    
                    U(DirectionalFeatures) = angle(nanmean(exp(1i*aa(:,DirectionalFeatures))));

                    bNotNaN = ~any(isnan( aa(:, DirectionalFeatures ) ),2);
                    iCantElementos = sum(bNotNaN);
                    iAux = nan(iCantElementos, iCantFeatures);
                    
                    if( ~isempty(NotDirectionalFeatures) )
                        iAux(:,NotDirectionalFeatures) = aa(bNotNaN,NotDirectionalFeatures) - repmat(U(NotDirectionalFeatures), iCantElementos,1);
                    end


                    for jj = DirectionalFeatures                                     
                        iAux1 = aa(bNotNaN, jj );
                        iAux2 = repmat(U(jj), iCantElementos,1);
                        iAux3 = [ iAux1 - iAux2, 2*pi + iAux1 - iAux2, iAux1 - 2*pi - iAux2 ];
                        [dummy, iMinIndex] = min(abs(iAux3), [], 2);
                        iMinIndex = sub2ind(size(iAux3),1:size(iAux3,1),iMinIndex');
                        iAux(:,jj) = iAux3(iMinIndex);
                    end 

%                     iAux1 = aa(:,DirectionalFeatures) - repmat(U(DirectionalFeatures), iCantElementos,1);
%                     boolAbs = abs( iAux1 ) > pi;
%                     SignX = sign( iAux1 );
% 
%                     bAux = boolAbs & SignX < 0;
% 
%                     aa( bAux, DirectionalFeatures ) = aa( bAux, DirectionalFeatures ) + 2*pi;
%                     aa( bAux, DirectionalFeatures ) = aa( bAux, DirectionalFeatures ) - 2*pi;
% 
%                     iAux(:,DirectionalFeatures) = rem( aa(:, DirectionalFeatures ) - repmat(U(DirectionalFeatures), iCantElementos,1) , 2*pi);

                    if(n == 1)
                        G = 1/iCantElementos*(iAux'*iAux);
                    else
                        G = 1/(iCantElementos-1)*(iAux'*iAux);
                    end

                end
                
            else

                U = nan(c,k);
                G = nan(k,k,c);
                
                if( isempty(DirectionalFeatures) )
                    
                    for ii = 1:c     
                        J = findnlab(a,ii);
                        if( ~isempty(J) )
                            aa = a(J,:);
                            bNotNaN = ~any(isnan( aa ),2);
                            U(ii,:) = nanmean( aa(bNotNaN,:) ,1);
                            if (nargout > 1 && sum(bNotNaN) > 1 )
                                G(:,:,ii) = covm(aa(bNotNaN,:),n);	
                            else
                                G(:,:,ii) = 0;
                            end
                        end
                    end
                    
                else
                
                    NotDirectionalFeatures = setdiff(1:k, DirectionalFeatures);
                    aa = +a;

                    U = nan(c, k);
                    
                    for ii = 1:c     
                        J = findnlab(a,ii);
                        if( ~isempty(J) )

                            if( ~isempty(NotDirectionalFeatures) )
                                U(ii,NotDirectionalFeatures) = nanmean(aa(J,NotDirectionalFeatures));							% 	in the usual way.
                            end
                            
                            U(ii,DirectionalFeatures) = angle(nanmean(exp(1i*aa(J,DirectionalFeatures))));
                            
                            bNotNaN = ~any(isnan( aa(J,:) ),2);
                            iCantElementos = sum(bNotNaN);
                            
                            if(iCantElementos < 2)
                                G(:,:,ii) = zeros(k);
                            else
                                
                                if (nargout > 1)

                                    iAux = nan(iCantElementos, k);

                                    if( ~isempty(NotDirectionalFeatures) )
                                        iAux(:,NotDirectionalFeatures) = aa(J(bNotNaN),NotDirectionalFeatures) - repmat(U(ii,NotDirectionalFeatures,:), iCantElementos,1);
                                    end

                                    for jj = DirectionalFeatures
                                        iAux1 = aa(J(bNotNaN), jj );
                                        iAux2 = repmat(U(ii,jj), iCantElementos,1);
                                        iAux3 = [ iAux1 - iAux2, 2*pi + iAux1 - iAux2, iAux1 - 2*pi - iAux2 ];
                                        [dummy, iMinIndex] = min(abs(iAux3), [], 2);
                                        iMinIndex = sub2ind(size(iAux3),1:size(iAux3,1),iMinIndex');
                                        iAux(:,jj) = iAux3(iMinIndex);
                                    end

                                    if(n == 1)
                                        G(:,:,ii) = 1/iCantElementos*(iAux'*iAux);
                                    else
                                        G(:,:,ii) = 1/(iCantElementos-1)*(iAux'*iAux);
                                    end
                                end
                            end
                        end
                    end
                end
			end
			labu = getlablist(a);
  	elseif (islabtype(a,'soft'))
  		problab = gettargets(a);
			% Here we also have to be careful for unlabeled data
			if (c==0)
				prwarning(2,'The dataset has soft labels but no targets defined: using targets 1');
				U = nanmean(+a);
				G = nancov(+a,n);
			else
				U = zeros(c,k);
				for ii = 1:c

					% Calculate relative weights for the means.
					g = problab(:,ii); nn = sum(g); g = g/mean(g); 

					U(ii,:) = mean(a.*repmat(g,1,k));	% Weighted mean vectors	

					if (nargout > 1)

						u  = mean(a.*repmat(sqrt(g),1,k));

						% this appears to be needed to weight cov terms properly
						G(:,:,ii) = covm(a.*repmat(sqrt(g),1,k),1) - U(ii,:)'*U(ii,:) + u'*u;

						% Re-normalise by M-1 if requested.
						if (n == 0)
							G(:,:,ii) = m*G(:,:,ii)/(m-1);
						end
					end
				end
            end
			labu = getlablist(a);
        else
			% Default action.
            U = mean(a);
            G = covm(a,n);
			labu = [];
        end

		% Add attributes of A to U.
		U = prdataset(U,labu,'featlab',getfeatlab(a), ...
                                'featsize',getfeatsize(a));
		if (~islabtype(a,'targets'))
			p = getprior(a);
			U = setprior(U,p); 
		end

	end;

return