%SUBSC Subspace Classifier
%
%   W = SUBSC(A,F)
%   W = A*SUBSC([],F)
%   W = A*SUBSC(F)
%
% INPUT
%   A    Dataset
%   F    Desired model dimensionality or fraction of retained 
%              variance per class
%
% OUTPUT
%   W    Subspace classifier  
%
% DESCRIPTION
% Each class in the trainingset A is described by linear subspace of
% dimensionality F (F>=1), or such that at least a fraction F (F<1) of its 
% varianceis retained. This is realised by calling PCAM(AI,F) or for
% each subset AI of A (objects of class I). For each class a model is
% built that assumes that the distances of the objects to the class
% subspaces follow a one-dimensional distribution. 
%
% New objects are assigned to the class of the nearest subspace.
% Classification by D = B*W, in which W is a trained subspace classifier
% and B is a testset, returns a dataset D with one-dimensional densities
% for each of the classes in its columns.
%
% If F is NaN it is optimised by REGOPTC.
%
% REFERENCE
% E. Oja, The Subspace Methods of Pattern Recognition, Wiley, New York, 1984.
%
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% DATASETS, MAPPINGS, PCAM, FISHERC, FISHERM, GAUSSM, REGOPTC

% 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

function W = subsc(varargin)
  
	mapname = 'SubspaceC';
  argin = shiftargin(varargin,'scalar');
  argin = setdefaults(argin,[],1);
  
  if mapping_task(argin,'definition')       % Define mapping
    W = define_mapping(argin,'untrained',mapname);
    
  elseif mapping_task(argin,'training')			% Train a mapping.
  
    [A,N] = deal(argin{:});
	
    if isnan(N)    % optimize regularisation parameter
      defs = {1};
      parmin_max = [1,size(A,2)];
      W = regoptc(A,mfilename,{N},defs,[1],parmin_max,testc([],'soft'),0);

    else

    % handle training like A*subsc, A*subsc([],3), subsc(A)
    % PRTools takes care that they are all converted to subsc(A,N)

      islabtype(A,'crisp');     % allow crisp labels only
      isvaldfile(A,1,2);        % at least one object per class, two objects
      A = testdatasize(A,'features');
      A = setprior(A,getprior(A));
      [m,k,c] = getsize(A);     % size of the training set
      for j = 1:c               % run over all classes
        B = seldat(A,j);        % get the objects of a single class only
        u = mean(B);            % compute its mean
        B = B - repmat(u,size(B,1),1); % subtract mean
        v = pcam(B,N);          % compute PCA for this class
        v = v*v';               % trick: affine mappings in the original space
        B = B - B*v;            % differences of objects and their mappings
        s = mean(sum(B.*B,2));  % mean square error w.r.t. the subspace
        data(j).u = u;          % store mean
        data(j).w = v;          % store mapping
        data(j).s = s;          % store mean square distance
      end
                                % define trained mapping,
                                % store class labels and size
      W = prmapping(mfilename,'trained',data,getlablist(A),k,c);
      W = setname(W,mapname);
      
    end

  elseif mapping_task(argin,'trained execution') % Evaluation

    % handle evaluation of a trained subspace classifier W for a dataset A.
    % The command D = A*W is by PRTools translated into D = subsc(A,W)
    
    [A,W] = deal(argin{1:2}); % get the variables
     m = size(A,1);           % number of test objects
    [k,c] = size(W);          % mapping size: from K features to C classes
     d = zeros(m,c);          % output: C class densities for M objects

    for j=1:c                 % run over all classes
      u = W.data(j).u;        % class mean in training set
      v = W.data(j).w;        % mapping to subspace in original space
      s = W.data(j).s;        % mean square distance
      B = A - repmat(u,m,1);  % substract mean from test set
      B = B - B*v;            % differences objects and their mappings
      d(:,j) = sum(B.*B,2)/s; % convert to distance and normalise 
    end
    d = exp(-d/2)/sqrt(2*pi); % convert to normal density

    A = prdataset(A);           % make sure A is a dataset 
    d = setdata(A,d,getlabels(W)); % take data from D and use 
                              % class labels as given in W
                              % other information in A is preserved   
    W = d;                    % return result in output variable W

  else     

    error('Illegal call')     % this should not happen

  end
      
return