%OPAM Orthogonal projection approach: selects "pure" features
%  
%
%  U = OPAM
%  U = OPAM([], OPT)
%
%  W = A*U
%  [W, SEL_VAR_DIS, DIS] = MAP(A, U)
%  [W, SEL_VAR_DIS, DIS] = OPAM(A, OPT)
%
%  C = B*W
%  C = MAP(B, W)
%  C = OPAM(B, W)
%
%
% INPUT 
%   A             [M -by- K] training dataset
%   B             [N -by- K] test dataset
%
%   U             [K -by- ?] untrained OPAM mapping
%   W             [K -by- K_SEL] trained OPAM feature selection mapping;
%                 indices of selected variables stored in W.DATA.IND  
%
%   OPT.          optional OPA settings
%    MAX_NUM      the max number of variables to select, INF means selecting as many variables as possible,a
%                 until the rank of dataset is exausted; at least one variable will be always selected; by delault = INF
%    EPS          stopping criteion; algorithm stops ERR < EPS, such features are not included into the set; by default = 1e-2 
%    VERBOSE      force/disable messages: PRPROGRESS(VERBOSE,['%s\n'],MSG) will be run;
%                 by default = []
%
% OUTPUT 
%   U               [K -by- ?] untrained OPAM mapping
%   W               [K -by- K_SEL] trained OPAM feature selection mapping;
%                   indices of selected variables stored in W.DATA.IND  
%   SEL_VAR_IND     [K_SEL -by- 1] vector of selected variables indices 
%   SEL_VAR_DIS     [K_SEL -by- 2] two vectors of selected variables disssimilarities 
%                   in the first column are absolute dissims, and in the second are relative dissim
%   DIS             [K_SEL -by- K -by- 2] two matrices which rows are dissimmilarity spectra; 
%                   each row contains dissim of all variables at particular selection step;
%                   the first matrix contains abs dissim, the second one containes relative dissim
%   ERR             [K_SEL -by- 1] relative squared error of the data projection to the subspace spanned by the select features
%   C               [N -by- K_SEL] dataset, result of mapping execution on dataset B
%
% DESRIPTION
% Implementation of the OPA algorithm as described in
% "Orthogonal Projection Approach Applied to Peak Purity Asessment"  
%  by F. C. Snachez, J. Toft, B. van den Bogaert, D. L. Massart, Anla. Chem. 68, 1, 79-85, 1996 
% (applied to the transposed data matrix)
%
% Absolute dissim measure  abs_diss_(n+1) = ||dx||^2*det(Y_n^T Y_n)
% Relative dissim measure  rel_diss_(n+1) = ||dx||^2
%
% Here, Y_n is the set of n selected features (columns), which are normalized to have a unit length, dx is the part of 
% a feature x which is (part) orthogonal to span(Y). If n=0 (diss(1)) the dissimilarity to the mean feature is measured.
%
% SEE ALSO
% SIMPLISMAM, VARIMAXFM, SIMLISMAPS, OPAPS, VARIMAXOM, ALS, OPA

% Copyright: S.Verzakov, s.verzakov@ewi.tudelft.nl 
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: opam.m,v 1.7 2007/01/31 21:52:19 serguei Exp $

function [w, sel_var_dis, dis, err] = opam(a,w)

% No dataset given: return untrained mapping.
if (nargin < 1) | (isempty(a))
  if nargin < 2
    w = [];
  end  
  w = mapping(mfilename,'untrained',{w});
  w = setname(w,'OPA Mapping');
  return
end

isdataset(a);           % Assert that A is a dataset.

% training
if nargin < 2 | ~isa(w,'mapping') 
  if nargin < 2
    w = [];  
  end  

  [m,k] = getsize(a);
  
  [z, sel_var_ind, sel_var_dis, dis, err] = opa((+a)',w);
  dis = permute(dis,[2 1 3]);

  w_data.ind = sel_var_ind;

  %labels = [num2str([1:length(sel_var_ind)]','%-1d'),  num2str(sel_var_ind(:),'(%-1d)')];
  labels = a.featlab(sel_var_ind,:);
  
  % Save all useful data
  w = mapping(mfilename,'trained',w_data,labels,size(a,2),length(sel_var_ind));
  w = setname(w,'OPA Feature Selection');
  
% applying
else 
  w_data = +w; % Unpack the mapping.
  %w = a*cmapm(w.size_in,w_data.ind);
  a = a(:,w_data.ind);
  %a.featlab = w.labels;
  units = specunits(a);
  a = remove_spectra_info(a);
  if units
    user = a.user;
    user.type = 'spectra';
    user.units = units;
    a.user = user;
  end  
  a.name = [a.name ': selected channels'];
  w = a;
end

return