%MEANCOV Estimation of the means and covariances from multiclass data
% 
%   [U,G] = MEANCOV(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, 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.7 2010/02/08 15:31:47 duin Exp $

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

		% 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 < 1
		U = prmapping(mfilename,'fixed');
		return
	end

	if (nargin < 2 || isempty(n))
		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 (isdouble(a))			% A is a matrix: compute mean and covariances
		U = mean(a);							% 	in the usual way.
		if nargout > 1
			G = prcov(a,n); 
		end
    
  elseif (isdatafile(a))
    
		if nargout  < 2
    	U = meancov_datafile(a,n);
		else
      [U,G] = meancov_datafile(a,n);
    end
    
  elseif (isdataset(a))
    
		[m,k,c] = getsize(a);
		if nargout == 2
			G = zeros(k,k,c);
		end
  	if (islabtype(a,'crisp'))
			
			if (c==0)  % special solution if all data is unlabeled
				U = mean(+a);
				if nargout > 1
					G = prcov(+a,n); 
				end
      else
        Jall = findclasses(a);
				for i = 1:c     
					%J = findnlab(a,i);
          J = Jall{i};
					if isempty(J)
						U(i,:) = NaN(1,k);
						if nargout > 1
							G(:,:,i) = NaN(k,k);
						end
					else
						U(i,:) = mean(a(J,:),1);
						if (nargout > 1)
							G(:,:,i) = covm(a(J,:),n);	
						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 = mean(+a);
				if nargout > 1
					G = prcov(+a,n);
				end
			else
				U = zeros(c,k);
				for i = 1:c
					g = problab(:,i);
					if mean(g) == 0 % all targets are zero: zero mean, zero cov.
						G(:,:,i) = zeros(k,k);
					else
						% Calculate relative weights for the means.
						g = g/mean(g); 
						U(i,:) = mean(a.*repmat(g,1,k));	% Weighted mean vectors	
						if (nargout > 1)
							u  = mean(a.*repmat(sqrt(g),1,k));
							% needed to weight cov terms properly
							G(:,:,i) = covm(a.*repmat(sqrt(g),1,k),1) - U(i,:)'*U(i,:) + u'*u;
							% Re-normalise by M-1 if requested.
							if (n == 0)
								G(:,:,i) = m*G(:,:,i)/(m-1);
							end
						end
					end
				end
			end
			labu = getlablist(a);
  	else
			% Default action.
  		U = mean(a);
			if nargout > 1
  			G = covm(a,n);
			end
			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,a.prior); 
		end
    
  else
    error('Illegal datatype')
    
	end

return



%MEANCOV Datafile overload

function [u,g] = meancov_datafile(a,n)

	if nargin < 2, n = []; end
	
  [m,k,c] = getsize(a);
  k = 0; % just for detecting first loop
	next = 1;
	while next > 0
		[b,next] = readdatafile(a,next);
    if k == 0
      k = size(b,2);
      if c == 0
        u = zeros(1,k);
      else
        u = zeros(c,k);
      end
      if nargout > 1
        if c == 0
          g = zeros(k,k);
        else
          g = zeros(k,k,c);
        end
      end
    end
		%compute for each call to readdatafile contributions to mean and cov
		bb = + b;
		if c == 0
			u = u + sum(bb,1);
		else
			nlab = getnlab(b);
			for j=1:size(b,1)
				if nlab(j) > 0
					u(nlab(j),:) = u(nlab(j),:) + bb(j,:);
					if nargout > 1
						g(:,:,nlab(j)) = g(:,:,nlab(j)) + bb(j,:)'*bb(j,:);
					end
				end
			end
		end
	end

	f = classsizes(a);
	u = u ./ repmat(f',1,k);

	if nargout == 2
		for j=1:c
			g(:,:,j) = g(:,:,j)/f(j) - u(j,:)'*u(j,:);
			if isempty(n) || n == 0
				g(:,:,j) = (f(j)/(f(j)-1))*g(:,:,j);
			end
		end
	end
	u = prdataset(u,getlablist(b));
	u.prior = b.prior;
return
