%PARZENDC Trainable classifier based on Parzen density estimation
% 
%  [W,H] = PARZENDC(A,H)
%  [W,H] = A*PARZENDC([],H)
%  [W,H] = A*PARZENDC(H)
% 
% INPUT
%  A   Dataset
%  H   Smoothing parameters (optional; default: estimated from A for each class)
%
% OUTPUT
%  W   Trained Parzen classifier
%  H   Smoothing parameters, estimated from the data
%
% DESCRIPTION
% For each of the classes in the dataset A, a Parzen density is estimated
% using PARZENML. For each class, a feature normalisation on variance is
% included in the procedure. As a result, the Parzen density estimate uses
% different smoothing parameters for each class and each feature.
%
% If a set of smoothing parameters H is specified, no learning is performed, 
% only the classifier W is produced. H should have the size of [C x K] if 
% A has C classes and K features. If the size of H is [1 x K] or [C x 1], 
% or [1 x 1], then identical values are assumed for all the classes and/or
% features.
%
% The densities for the points of a dataset B can be found by D = B*W.
% D is an [M x C] dataset, if B has M objects.
% 
% EXAMPLES
% See PREX_DENSITY for densities and PREX_PARZEN for differences between
% PARZENC, PARZENDC and PARZENM.
%
%
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% DATASETS, MAPPINGS, PARZENC, PARZENM, PARZENML, PREX_DENSITY
 
% 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,h] = parzendc(varargin)
  
	mapname = 'ParzenD';
  argin = shiftargin(varargin,'scalar');
  argin = setdefaults(argin,[],[]);
  
  if mapping_task(argin,'definition')
    
    W = define_mapping(argin,'untrained',mapname);
    
  elseif mapping_task(argin,'training')			% Train a mapping.
  
    [a,h] = deal(argin{:});
    islabtype(a,'crisp','soft');
  
    % remove too small classes, escape in case no two classes are left
    [a,m,k,c,lablist,L,W] = cleandset(a,1); 
    if ~isempty(W), h = NaN; return; end
      
    a = testdatasize(a);
    a = testdatasize(a,'objects');

    if ~isempty(h)       % Take user settings for smoothing parameters.

      if size(h,1) == 1, h = repmat(h,c,1); end
      if size(h,2) == 1, h = repmat(h,1,k); end
      if any(size(h) ~= [c,k])
        error('Array with smoothing parameters has a wrong size.');
      end

    else   % Estimate smoothing parameters

      % Scale A such that its mean is shifted to the origin and 
      % the variances of all features are scaled to 1. 
      ws = scalem(a,'variance');
      b = a*ws;  
      % SCALE is basically [1/mean(A) 1/STD(A)] based on the properties of SCALEM.
      scale = ws.data.rot;				
      if (size(scale,1) ~= 1) % formally ws.data.rot stores a rotation matrix 
        scale = diag(scale)'; % extract the diagonal if it does,
      end                     % otherwise we already have it

      h = zeros(c,k);
      s = sprintf('parzendc: smoothing per class ');
      prwaitbar(c,s);
      for j=classuse(b,2)
        prwaitbar(c,j,[s int2str(j)]);
        bb = seldat(b,j)*remclass; 	% BB consists of the j-th class only.
        h(j,:) = repmat(parzenml(bb),1,k)./scale;
      end
      prwaitbar(0);
    end

    W = prmapping('parzen_map','trained',{a,h,getprior(a)},getlablist(a),k,c);
    W = setname(W,mapname);
    W = allclass(W,lablist,L);   % complete classifier with missing classes
    W = setcost(W,a);
    
  end

return;
