%PCAKLM Principal Component Analysis/Karhunen-Loeve Mapping
%       (PCA or MCA of overall/mean covariance matrix)
% 
%   [W,FRAC] = PCAKLM(TYPE,A,N)
%   [W,N]    = PCAKLM(TYPE,A,FRAC)
%
% INPUT
%  A           Dataset
%  TYPE        Type of mapping: 'pca' or 'klm'. Default: 'pca'.
%  N or FRAC   Number of dimensions (>= 1) or fraction of variance (< 1) 
%              to retain; if > 0, perform PCA; otherwise MCA.
%              Default: N = inf.
%
% OUTPUT
%  W           Affine Karhunen-Loeve mapping
%  FRAC or N   Fraction of variance or number of dimensions retained.
%
% DESCRIPTION
% Performs a principal component analysis (PCA) or minor component analysis
% (MCA) on the overall or mean class covariance matrix (weighted by the
% class prior probabilities). It finds a rotation of the dataset A to an
% N-dimensional linear subspace such that at least (for PCA) or at most (for
% MCA) a fraction FRAC of the total variance is preserved.
%
% PCA is applied when N (or FRAC) >= 0; MCA when N (or FRAC) < 0. If N is 
% given (abs(N) >= 1), FRAC is optimised. If FRAC is given (abs(FRAC) < 1), 
% N is optimised. 
%
% Objects in a new dataset B can be mapped by B*W, W*B or by A*KLM([],N)*B.
% Default (N = inf): the features are decorrelated and ordered, but no 
% feature reduction is performed.
%
% ALTERNATIVE
%
%   V = PCAKLM(A,TYPE,0)
% 
% Returns the cumulative fraction of the explained variance. V(N) is the 
% cumulative fraction of the explained variance by using N eigenvectors.
%
% This function should not be called directly, only trough PCA or KLM.
% Use FISHERM for optimizing the linear class separability (LDA).
% 
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% MAPPINGS, DATASETS, PCLDC, KLLDC, PCAM, KLM, FISHERM

% Copyright: R.P.W. Duin, r.p.w.duin@37steps.com
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: pcaklm.m,v 1.15 2010/02/08 15:29:48 duin Exp $

function [w,truefrac] = pcaklm (type,varargin)

  truefrac = [];
    
  argin = shiftargin(varargin,'scalar'); % allows for call like pcam(0.99)
  [a,frac] = setdefaults(argin,[],inf);
  type = setdefaults(type,'pcam');

  if (strcmp(type,'pcam'))
    mapname = 'PCA';
  elseif (strcmp(type,'klm'))
    mapname = 'Karhunen-Loeve Mapping';
  else
    error('Unknown type specified');
  end

  %DXD Make the name a bit more informative:
  if isfinite(frac)
    if (frac<1)
      mapname = [mapname sprintf(' ret. %4.1f%% var',100*frac)];
    else
      mapname = [mapname sprintf(' to %dD',frac)];
    end
  end

  % Empty dataset: return straightaway.
  if (nargin < 2) || (isempty(a))
    w = prmapping(type,frac);
    w = setname(w,mapname);
    return
  end
  
  %nodatafile(a);
  if ~isdataset(a)
    if isa(a,'double')
      a = prdataset(a,1);   % make sure we have a dataset
    else
      error('nodatafile','Data should be given in a dataset or as doubles')
    end
  end
  
  islabtype(a,'crisp','soft');
  %isvaldfile(a,1);       % at least 1 object per class
  a = setfeatdom(a,[]);  % get rid of domain testing

  [m,k,c] = getsize(a);
  p = getprior(a);
  a = setprior(a,p);  % make class frequencies our prior

  % If FRAC < 0, perform minor component analysis (MCA) instead of 
  % principal component analysis.
  mca = (frac < 0); frac = abs(frac);

  % Shift mean of data to origin.
  b = a*scalem(a); 
  
  % If there are less samples M than features K, first perform a lossless
  % projection to the (M-1) dimensional space spanned by the samples.
  if (m <= k)
    testdatasize(b,'objects');
    u = reducm(b); 
    b = b*u;
    k = size(b,2);
    frac = min(frac,k);
  else
    testdatasize(b,'features');
    u = [];
  end

  % Calculate overall or average class prior-weighted covariance matrix and
  % find eigenvectors F. 

  if (strcmp(type,'pcam'))
    if (c==0) || ~islabtype(a,'crisp')  % we have unlabeled data!
      G = prcov(+b); % use all
    else
      bb = [];
      classsiz = classsizes(b);
      L = find(classsiz > 0);
      U = zeros(1,k);
      Q = zeros(k,k);
      for j = L
        bj = +(seldat(b,j)*remclass);
        U = U + p(j)*mean(bj,1);
        Q = Q + p(j)*(bj'*bj)/classsiz(j);
      end
      G = Q - U'*U;
    end
  else
    %DXD For high dimensional dataset with many classes, we cannot
    %store all individual cov. matrices in memory (like in the next
    %line), but we have to compute them one by one:
    %[U,GG] = meancov(b,1);
    G = zeros(k,k);
    L = classuse(b);
    for i = L
      [U,GG] = meancov(seldat(b,i)*remclass,1);
      G = G + p(i)*GG;
    end
  end
  [F,V] = preig(G); % overdone if reducm has been called

  % v = V(I) contains the sorted eigenvalues:
  % descending for PCAM, ascending for MCA.
  if (mca)
    [v,I] = sort(diag(V));
  else
    [v,I] = sort(-diag(V)); v = -v;
  end
  
  if (frac == inf)           % Return all dimensions, decorrelated and ordered.
    n = k; truefrac = k;
  elseif (frac == 0)         % Just return cumulative retained variance.
    w = cumsum(v)/sum(v);
    return
  elseif (frac >= 1)         % Return FRAC dimensions.
    n = abs(frac); 
    if (n > k),
      error('illegal dimensionality requested');
    end
    I = I(1:n);
    sv = sum(v); 
    if (sv ~= 0),
      truefrac = cumsum(v(1:n))/sv;
    else,
      truefrac = 0;
    end;
  elseif (frac > 0)    % Return the N dimensions that retain at least (PCA)
                       % or at most (MCA) FRAC variance.
    J = find(cumsum(v)/sum(v) > frac);
    if (mca), n = J(1)-1; else, n = J(1); end;
    truefrac = n; I = I(1:n);
  end

  % If needed, apply pre-calculated projection to (M-1) dimensional subspace.
  if (~isempty(u))
    rot = u.data.rot*F(:,I); 
    off = u.data.offset*F(:,I);
  else
    rot = F(:,I); 
    off = -mean(a*F(:,I));
  end

  % Construct affine mapping.
  
  w = affine(rot,off,a);
  w = setdata(w,v,'eigenvalues');
  w = setname(w,mapname);
    
return
