%MEANCOV Fixed mapping estimating means and covariances of multiclass data
% 
%   [U,G] = MEANCOV(A,N)
%   [U,G] = A*MEANCOV([],N)
%   [U,G] = A*MEANCOV(N)
% 
%  INPUT
%   A	  Dataset
%   N	  Normalization to be useduse 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, given as a dataset with all annotation of A
%   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 (<a href="http://37steps.com/prtools">PRTools Guide</a>) 
% 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(varargin)

	argin = shiftargin(varargin,'scalar');
  argin = setdefaults(argin,[],0);
  
  if mapping_task(argin,'definition')
    
    U = define_mapping(argin,'fixed');
    
  elseif mapping_task(argin,'training')			% Compute mapping.
  
    [a,n] = deal(argin{:});

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

  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