Index: /distools/Contents.m
===================================================================
--- /distools/Contents.m	(revision 17)
+++ /distools/Contents.m	(revision 18)
@@ -53,5 +53,6 @@
 %  Routines supporting in learning from dissimilarity matrices
 %  -----------------------------------------------------------------------
-%  CROSSVALD    Cross-validation error for dissimilarity representations
+%  CROSSVALD    Cross-validation error for dissimilarity data
+%  CLEVALD      Classifier evaluation (learning curve) for dissimilarity data
 %  DISSPACES    Compute various spaces out of a dissimilarity matrix
 %  GENDDAT      Generate random training and test sets for dissimilarity data
Index: /distools/clevald.m
===================================================================
--- /distools/clevald.m	(revision 18)
+++ /distools/clevald.m	(revision 18)
@@ -0,0 +1,237 @@
+%CLEVALD Classifier evaluation (learning curve) for dissimilarity data
+%
+%   E = CLEVALD(D,CLASSF,TRAINSIZES,REPSIZE,NREPS,T)
+%
+% INPUT
+%   D          Square dissimilarity dataset
+%   CLASSF     Classifiers to be evaluated (cell array)
+%   TRAINSIZE  Vector of class sizes, used to generate subsets of D
+%              (default [2,3,5,7,10,15,20,30,50,70,100])
+%   REPSIZE    Representation set size per class (>=1), or fraction (<1)
+%              (default total, training set)
+%   NREPS      Number of repetitions (default 1)
+%   T          Test dataset (default [], use remaining samples in A)
+%
+% OUTPUT
+%   E          Error structure (see PLOTE) containing training and test
+%              errors
+%
+% DESCRIPTION
+% Generates at random, for all class sizes defined in TRAINSIZES, training
+% sets out of the dissimilarity dataset D. The representation set is either
+% equal to the training set (REPSIZE = []), or a fraction of it (REPSIZE  <1)
+% or a random subset of it of a given size (REPSIZE>1). This set is used
+% for training the untrained classifiers CLASSF. The resulting trained
+% classifiers are tested on the training objects and on the left-over test 
+% objects, or, if supplied, the testset T. This procedure is then repeated 
+% NREPS times.
+%
+% The returned structure E contains several fields for annotating the plot
+% produced by PLOTE. They may be changed by the users. Removal of the field
+% 'apperror' (RMFIELD(E,'apperror')) suppresses the draw of the error
+% curves for the training set.
+%
+% Training set generation is done "with replacement" and such that for each
+% run the larger training sets include the smaller ones and that for all
+% classifiers the same training sets are used.
+%
+% This function uses the RAND random generator and thereby reproduces 
+% if its seed is reset (see RAND).
+% If CLASSF uses RANDN, its seed should be reset as well.
+%
+% SEE ALSO
+% MAPPINGS, DATASETS, CLEVAL, TESTC, PLOTE
+
+% R.P.W. Duin, r.p.w.duin@prtools.org
+% Faculty EWI, Delft University of Technology
+% P.O. Box 5031, 2600 GA Delft, The Netherlands
+
+function e = cleval(a,classf,learnsizes,repsize,nreps,t)
+
+	prtrace(mfilename);
+
+  if (nargin < 6) 
+    t = [];
+  end;
+  if (nargin < 5) | isempty(nreps);
+    nreps = 1;
+  end;
+  if (nargin < 4) 
+    repsize = [];
+  end
+  if (nargin < 3) | isempty(learnsizes);
+    learnsizes = [2,3,5,7,10,15,20,30,50,70,100];
+  end;
+  if ~iscell(classf), classf = {classf}; end
+
+  % Assert that all is right.
+  isdataset(a); issquare(a); ismapping(classf{1}); 
+  if (~isempty(t)), isdataset(t); end
+
+  % Remove requested class sizes that are larger than the size of the
+  % smallest class.
+
+  mc = classsizes(a); [m,k,c] = getsize(a);
+  toolarge = find(learnsizes >= min(mc));
+  if (~isempty(toolarge))
+    prwarning(2,['training set class sizes ' num2str(learnsizes(toolarge)) ...
+                 ' larger than the minimal class size in A; remove them']);
+    learnsizes(toolarge) = [];
+  end
+  learnsizes = learnsizes(:)';
+
+  % Fill the error structure.
+
+  nw = length(classf(:));
+  datname = getname(a);
+
+  e.n       = nreps;
+  e.error   = zeros(nw,length(learnsizes));
+  e.std     = zeros(nw,length(learnsizes));
+  e.apperror   = zeros(nw,length(learnsizes));
+  e.appstd     = zeros(nw,length(learnsizes));
+  e.xvalues = learnsizes(:)';
+  e.xlabel = 'Training set size per class';
+  e.names   = [];
+  if (nreps > 1)
+    e.ylabel= ['Averaged error (' num2str(nreps) ' experiments)'];
+  elseif (nreps == 1)
+    e.ylabel = 'Error';
+  else
+    error('Number of repetitions NREPS should be >= 1.');
+  end;
+  if (~isempty(datname))
+    if isempty(repsize)
+      e.title = [datname ', Rep. Set = Train Set'];
+    elseif repsize < 1
+      e.title = [datname ', Rep. size = ' num2str(repsize) ' Train size'];
+    else
+      e.title = [datname ', Rep. size = ' num2str(repsize) ' per class'];
+    end
+  end
+  if (learnsizes(end)/learnsizes(1) > 20)
+    e.plot = 'semilogx';        % If range too large, use a log-plot for X.
+  end
+
+  % Report progress.
+	
+	s1 = sprintf('cleval: %i classifiers: ',nw);
+	prwaitbar(nw,s1);
+
+  % Store the seed, to reset the random generator later for different
+  % classifiers.
+
+	seed = rand('state');
+
+  % Loop over all classifiers (with index WI).
+
+  for wi = 1:nw
+		
+    if (~isuntrained(classf{wi}))
+      error('Classifiers should be untrained.')
+    end
+    name = getname(classf{wi});
+    e.names = char(e.names,name);
+    prwaitbar(nw,wi,[s1 name]);
+
+    % E1 will contain the error estimates.
+
+    e1 = zeros(nreps,length(learnsizes));
+    e0 = zeros(nreps,length(learnsizes));
+
+    % Take care that classifiers use same training set.
+
+    rand('state',seed); seed2 = seed;
+
+		% For NREPS repetitions...
+		
+		s2 = sprintf('cleval: %i repetitions: ',nreps);
+		prwaitbar(nreps,s2);
+
+		for i = 1:nreps
+	
+			prwaitbar(nreps,i,[s2 int2str(i)]);
+      % Store the randomly permuted indices of samples of class CI to use in
+      % this training set in JR(CI,:).
+
+			JR = zeros(c,max(learnsizes));
+			
+			for ci = 1:c
+
+				JC = findnlab(a,ci);
+
+        % Necessary for reproducable training sets: set the seed and store
+        % it after generation, so that next time we will use the previous one.
+				rand('state',seed2);
+
+				JD = JC(randperm(mc(ci)));
+				JR(ci,:) = JD(1:max(learnsizes))';
+				seed2 = rand('state'); 
+			end
+
+			li = 0;										% Index of training set.
+			
+			nlearns = length(learnsizes);
+			s3 = sprintf('cleval: %i sizes: ',nlearns);
+			prwaitbar(nreps,s3);
+			
+			for j = 1:nlearns 
+				
+				nj = learnsizes(j);
+				
+				prwaitbar(nlearns,j,[s3 int2str(j) ' (' int2str(nj) ')']);
+				li = li + 1;
+
+        % J will contain the indices for this training set.
+
+        J = [];
+        R = [];
+        for ci = 1:c
+          J = [J;JR(ci,1:nj)'];
+          if isempty(repsize)
+            R = [R JR(ci,1:nj)];
+          elseif repsize < 1
+            R = [R JR(ci,1:ceil(repsize*nj))];
+          else
+            R = [R JR(ci,1:min(nj,repsize))];
+          end
+            
+        end;
+
+				w = a(J,R)*classf{wi}; 					% Use right classifier.
+				e0(i,li) = a(J,R)*w*testc;
+				if (isempty(t))
+  				Jt = ones(m,1);
+					Jt(J) = zeros(size(J));
+					Jt = find(Jt); 								% Don't use training set for testing.
+					e1(i,li) = a(Jt,R)*w*testc;
+				else
+					e1(i,li) = t(:,R)*w*testc;
+				end
+
+			end
+			prwaitbar(0);
+
+		end
+		prwaitbar(0);
+
+    % Calculate average error and standard deviation for this classifier
+    % (or set the latter to zero if there's been just 1 repetition).
+
+		e.error(wi,:) = mean(e1,1);
+		e.apperror(wi,:) = mean(e0,1);
+		if (nreps == 1)
+			e.std(wi,:) = zeros(1,size(e.std,2));
+			e.appstd(wi,:) = zeros(1,size(e.appstd,2));
+		else
+			e.std(wi,:) = std(e1)/sqrt(nreps);
+			e.appstd(wi,:) = std(e0)/sqrt(nreps);
+		end
+	end
+	prwaitbar(0);
+
+	% The first element is the empty string [], remove it.
+	e.names(1,:) = [];
+
+return
+
Index: /distools/disnorm.m
===================================================================
--- /distools/disnorm.m	(revision 17)
+++ /distools/disnorm.m	(revision 18)
@@ -1,10 +1,10 @@
 %DISNORM Normalization of a dissimilarity matrix
 %
-% 	V = disnorm(D,OPT)
+% 	V = DISNORM(D,OPT)
 %   F = E*V
 %
 % INPUT
 %   D 	 NxN dissimilarity matrix or dataset, which sets the norm
-%   E    Matrix to be normalized
+%   E    Matrix to be normalized, e.g. D itself
 % 	OPT  'max' : maximum dissimilarity is set to 1 by global rescaling
 % 		   'mean': average dissimilarity is set to 1 by global rescaling (default)
@@ -17,4 +17,10 @@
 %   OPT = 'mean'
 %
+% DESCRIPTION
+% Operation on dissimilarity matrices, like the computation of classifiers
+% in dissimilarity space, may depend on the scaling of the dissimilarities
+% (a single scalar for the entire matrix). This routine computes a scaling
+% for a giving matrix, e.g. a training set and applies it to other
+% matrices, e.g. the same training set or based on a test set.
 
 % Copyright: Elzbieta Pekalska, ela.pekalska@googlemail.com
@@ -41,5 +47,5 @@
 %DEFINE mapping
 if isstr(opt)
-	discheck(D);
+%	discheck(D);
 	opt = lower(opt);
 	if strcmp(opt,'mean')
Index: /distools/distools_long.m
===================================================================
--- /distools/distools_long.m	(revision 18)
+++ /distools/distools_long.m	(revision 18)
@@ -0,0 +1,137 @@
+%%   DisTools Table of Contents
+%%  23-Nov-2009_17:27
+%%
+%%  This Matlab toolbox for the analysis of dissimilarity data works only 
+%%  if also the pattern recognition toolbox PRTools is available.
+%%  See http://prtools.org
+%%
+%%    E. Pekalska,  ela.pekalska@googlemail.com, University of Manchester
+%%    R.P.W. Duin,  r.duin@ieee.org, Delft University of Technology
+%%
+%% Characterization of dissimilarity matrices
+%% ------------------------------------------
+%% CHECKEUCL    Check whether a square dissimilarity matrix has a Euclidean behavior
+%% CHECKTR      Check whether a square dissimilarity matrix obeys triangle inequality
+%% CHARDMAT     Fiand several characteristic of (dis)similarity data
+%% CORRTR       Correct a square dissimilarity matrix to obey the triangle inequality
+%% DISCHECK     Dissimilarity matrix check
+%% DISNORM      Normalization of a dissimilarity matrix
+%% DISSTAT      Basic statistics of the dissimilarity matrix
+%  GOFCL        Goodness of clusters/classes separability vs compactness for dissimilarity data
+%  INTRDIM      Estimate Intrinsic dimension from dissimilarity data
+%% ISSQUARE     Check whether a matrix is square
+%% ISSYM        Check whether a matrix is symmetric
+%% ASYMMETRY    Compute asymmetry of dissimilarity matrix
+%% NEF          Negative Eigen Fraction
+%% NNE          Leave-one-out Nearest Neighbor error on a dissimilarity matrix
+%% NNERR        Exact expected NN error from a dissimilarity matrix
+%  NNERROR      Exact expected NN error from a dissimilarity matrix (2)
+%  VAT          Visual Assessment of cluster Tendency for dissimilarity matrices	
+% 
+%% 
+%% Dissimilarity Measures
+%% -----------------------------------------------
+%  BINDISTM     Dissimilarity matrix between binary vectors
+%  BLURDISTM    Blurred Euclidean distance matrix between blobs
+%   BLOBBOX     Find box around a binary blob and resample
+%  CORRDISTM    Distance matrix based on correlations
+%% COSDISTM     Distance matrix based on inner products
+%  DPROCRUSTDM  Distance matrix between datasets based on extended Procrustes problem
+%  EDGEDISTM    Distance matrix between images based on their edges
+%  EDITDISTM    Edit distance matrix between strings
+%% EUDISTM      Euclidean distance matrix
+%  EXPDISTM     Exponential-type of distance matrix
+%  FLPDISTM     Fast computation of the lp (p > 0)  distance matrix
+%% HAMDISTM     Hamming distance matrix between binary vectors
+%% HAUSDM       Hausdorff and modified Hausdorff distance between datasets of image blobs
+%  JACSIMDISTM  Jaccard-like distance matrix based on similarities
+%  LPDISTM      l_p (p > 0) distance matrix
+%  QDISTM       Distance matrix for quantitative variables
+%  RANKDISTM    Distance matrix between two data sets based on ranking
+%  SAMDISTM     Distance matrix based on Spectral Angular Mapper (SAM)
+%  STRKERM      String Kernel Matrix by Lodhi et al
+% 
+%% 
+%% Transformations and projections
+%% -----------------------------------------------
+%% DISSIMT      Fixed DISsimilarity-SIMilarity transformation
+%  KCENTERM     Kernel weighted centering mapping (also for a similarity matrix) 
+%% MAKESYM      Make a matrix symmetric
+%  PROXXM       Proximity mapping
+%  SIGMOID      Element-wise sigmoid tranformation of a matrix
+%  FASTMAPD     FastMap; inear projection of Euclidean distances
+%% PE_EM        Pseudo-Euclidean embedding (includes Classical Scaling as a special case)	
+%  SPHEM        Spherical Embedding
+% 
+%% 
+%% Classification in Pseudo-Euclidean Space and indefinite kernels
+%% -----------------------------------------------
+%  PE_AFFINE
+%% SETSIG       Set PE signature for mappings or datasets
+%% GETSIG       Set PE signature for mappings or datasets
+%% ISPE_DATASET Test dataset for PE signature setting 
+%% ISPE_EM      Test mapping for PE signature setting
+%% PE_DISTM     Square pseudo-Euclidean distance between two datasets
+%% PE_KERNELM   Compute kernel in PE space
+%  PE_LIBSVC    Libsvc for PE spaces
+%% PE_MTIMES    Matrix multiplication (inner product) in PE space
+%% PE_PARZENC   Parzen classifier in PE space
+%% PE_KNNC      KNN classifier in PE space
+%% PE_NMC       Nearest mean classifier in PE space
+%% PE_EM        Pseudo-Euclidean linear embedding of dissimilarities
+%  KEM          Kernel embedding
+%% PLOTSPECTRUM Plot spectrum of eigenvalues
+%  PSPCA        Pseudo-Euclidean Principal Component Analysis
+%
+%  Indefinte kernel routines
+%  -------------------------
+%  IKFD         Indefinite Kernel Fisher discriminant
+%  IKPCA        Indefinite Kernel PCA
+%
+%% 
+%% Routines supporting in learning from dissimilarity matrices
+%% -----------------------------------------------------------------------
+%% CROSSVALD    Cross-validation error for dissimilarity representations
+%  CLEVALD      Classifier evaluation (Learning curve)
+%% DISSPACES    Compute various spaces out of a dissimilarity matrix
+%% GENDDAT      Generate random training and test sets for dissimilarity data
+%% GENREP       Generate a representation set
+%% GENREPI      Generate indices for representation, learning and testing sets
+%% SELCDAT      Select Class Subset from a Square Dissimilarity Dataset
+%% PROTSELFD    Forward prototype selection   
+%  AUCDLPC      AUC-LP classifier on dissimilarity data
+%% DLPC         LP-classifier on dissimilarity (proximity) data
+%  DRSSCC       Dissimilarity-based random subspace combining classifier 
+%% KNNDC        K-Nearest Neighbor classifier for dissimilarity matrices
+%% PARZENDDC    Parzen classifier for dissimilarity matrices
+%  KFD          Kernel Fisher Discriminant
+%  KSVC         Kernel Support Vector classifier on a kernel matrix
+% 	KSVO        Kernel Support Vector Optimizer
+%  KSVC_NU      Kernel Support Vector classifier on a kernel matrix; nu-version 
+% 	KSVO_NU     Kernel Support Vector Optimizer; nu-version 
+%  MCLASSDC     Multi-Class Dissimilarity-based Classifier from Two-Class Discriminants
+%% TESTKD       Test k-NN classifier for dissimilarity data
+%% TESTPD       Test Parzen classifier for dissimilarity data
+%  TQDC         Trade-off Quadratic Discriminant (Regularized Bayes Normal Classifier)
+% 
+% 
+%  Graphs and distances
+%  -----------------------------------------------
+%  DISTGRAPH    Computes distances in a graph
+%  DMSTSPM      Finds the shortest paths along K minimum spanning trees
+%  DSPATH       Single shortest path in a (dissimilarity) graph
+%  DSPATHS      All shortest paths in a (dissimilarity) Graph
+%  GRAPHPATH    Compute shortest paths in a graph
+%  KMST         Finds K minimum spanning trees based on a distance matrix
+%  MSTPLOT      Plot minimum spanning trees
+%  NHGRAPH      Find a neighborhood graph and its shortest paths
+%  PLOTGRAPH    Plot a 2D graph
+%%
+%% EXAMPLES
+%% --------
+%% CROSSVALD_EX Crossvalidation of several classifiers
+
+%  Superfluous / outdated but still available
+%  -------------------------------------------------
+%  KPCA, AUGPSEM, PSEM
+%  
Index: /distools/kem.m
===================================================================
--- /distools/kem.m	(revision 18)
+++ /distools/kem.m	(revision 18)
@@ -0,0 +1,212 @@
+%KEM Kernel Embedding
+%
+%     W = KEM(D,ALF,P)
+%                 OR
+%     W = KEM(W,ALF)
+%
+% INPUT
+%   D   NxN kernel matrix (dataset)
+%   W   Trained linear embedding into a pseudo-Euclidean space
+%   ALF Parameter determining the dimensionality and the mapping (optional, default: Inf)
+%       (0,1)   - Fraction of the total (absolute value) preserved variance
+%        Inf    - No dimensionality reduction, keeping all dimensions (it's VERY noisy)
+%       'p'     - Projection into a Euclidean space based on positive eigenvalues only
+%       'PARp'  - Projection into a Euclidean space based on the PAR fraction of
+%                 positive eigenvalues; e.g. ALF = '0.9p'
+%       'n'     - Projection into a Euclidean space based on negative eigenvalues only
+%       'PARn'  - Projection into a (negative) Euclidean space based on the PAR fraction
+%                 of negative eigenvalues; e.g. ALF = '0.7n'
+%       'P1pP2n'- Projection into a Euclidean space based on the P1 positive eigenvalues
+%                 and P2 negative eigenvalues; e.g. ALF = '0.7p0.1n', ALF = '7p2n'
+%       1 .. N  - Number of dimensions in total
+%       [P1 P2] - P1 dimensions or preserved fraction of variance in the positive subspace
+%                 and P2 dimensions or preserved fraction of variance in the negative
+%                 subspace; e.g. ALF = [5 10], ALF = [0.9 0.1]
+%   P   Integer between 0 and N specifying which object is mapped at the origin;
+%       0 stands for the mean; (optional, default: 0)
+%
+% OUTPUT
+%   W   Linear embedding into a pseudo-Euclidean space
+%
+% DEFAULT
+%   P   = 0
+%   ALF = INF
+%
+% DESCRIPTION
+% Linear mapping W onto an M-dimensional Pseudo-Euclidean (PE) subspace from a 
+% symmetric, square kernel matrix D such that the similarities are preserved.
+% M is determined by ALF. E.g., the subspace is found such that at least a fraction
+% ALF of the total variance is preserved for ALF in (0,1). The resulting X is found
+% by D*W. The signature of the obtained PE space (numbers of positive and negative
+% directions) can be found by PE_SIG(W). The spectrum of the obtained space
+% can be found by PE_SPEC(W). 
+%
+% A trained mapping can be reduced further by:   W = KEM(W,ALF)
+% The signature of the obtained PE space can be found by PE_SIG(W)
+% The spectrum of 
+%
+% SEE ALSO
+% MAPPINGS, DATASETS, PE_EM, AUGPSEM, PCA, PE_PCA, PE_SPEC, GETSIG, SETSIG
+%
+% LITERATURE
+% 1. L. Goldfarb, A unified approach to pattern recognition, Pattern Recognition, vol.17,
+% 575-582, 1984.
+% 2. E. Pekalska, P. Paclik, and R.P.W. Duin, A Generalized Kernel Approach to
+% Dissimilarity-based Classification, Journal of Machine Learning Research,
+% vol.2, no.2, 175-211, 2002.
+
+% Copyright: Elzbieta Pekalska, ela.pekalska@googlemail.com
+% Faculty EWI, Delft University of Technology and
+% School of Computer Science, University of Manchester
+
+% This routine is adapted from PE_EM by R.P.W. Duin. 
+
+function [W,sig,L,Q] = kem(d,alf,pzero,prec)
+
+% PREC is the precision parameter used for the automatic
+% selection (heuristic) of the number of dominant eigenvalues.
+% This happens when SELEIGS is called with the parameter 'CUT'.
+
+if nargin < 4, prec = [];  end
+if nargin < 3, pzero = []; end
+if nargin < 2  alf = [];   end
+
+if nargin < 1 | isempty(d),
+  W = mapping(mfilename,{alf,pzero,prec});
+  W = setname(W,'PE embedding');
+  return
+end
+
+if isempty(prec),  prec = 1e-4; end
+if isempty(pzero), pzero = 0;   end
+if isempty(alf),   alf = inf;   end
+
+if (isdataset(d) | isa(d,'double'))
+  if ismapping(alf)
+    % APPLY THE MAPPING
+    [m,n] = size(d);
+    %d     = d.^2; % OK for kernels???
+
+    Q   = getdata(alf,'evec');  % Eigenvectors
+    me  = getdata(alf,'mdis');  % Vector of the average squared original dissimilarities
+    p   = getdata(alf,'mean');  % p=0 -> the mean of the embedded configuration lies at 0, 
+                                % otherwise, it lies at pzero
+    L   = getdata(alf,'eval');  % Eigenvalues
+
+    % Project new data depending on p
+    % (whether the mean or other object lies at the origin)
+    if p == 0,
+      H = -repmat(1,n,n)/n;
+      H(1:n+1:end) = H(1:n+1:end) + 1;     % H = eye(n) - ones(n,n)/n
+      W = -0.5 * (d - me(ones(m,1),:)) * H * Q * diag(sqrt(abs(L))./L);
+    else
+      W =  0.5 * (d(:,p) * ones(1,n) + me(ones(m,1),:) - d) * Q * diag(sqrt(abs(L))./L);
+    end
+
+    % Store signature in the USER field
+    if isdataset(W),
+      W = setname(W,['Projected ' updname(W.name)]);
+      W = setsig(W,getdata(alf,'sig'));
+    end
+    return
+  end
+end
+
+
+
+% REDUCE A TRAINED MAPPING
+if ismapping(d)
+  data = getdata(d);
+  if iscell(data) % check for old type of mapping
+    ispsem(d);
+    dat.evec = data{1};
+    dat.mdis = data{2};
+    dat.mean = data{3};
+    dat.eval = data{4};
+    dat.sig  = data{5};
+    dat.prec = data{6};
+    d = setmapping_file(d,mfilename);
+    d = setdata(d,dat);
+    if nargin < 2
+      W = d;
+      return % conversion only
+    else 
+      data = dat;
+    end
+  end
+  Q   = data.evec;  % Eigenvectors
+  L   = data.eval;  % Eigenvalues
+  m    = size(Q,1);
+
+  [ll,K] = sort(-abs(L));
+  L = L(K);
+  Q = Q(:,K);
+  [J,sig] = seleigs(L,alf,getdata(d,'prec'));
+  data.evec = Q(:,J);        % Eigenvectors
+  data.eval = L(J);          % Eigenvalues
+ 
+  W = mapping(mfilename,'trained',data,[],m,length(J));
+  W = setname(W,'PE embedding');
+  return
+end
+
+
+
+% TRAIN THE MAPPING
+% Tolerance value used in comparisons
+% if mean(+d(:)) < 1,
+%   tol = 1e-12;
+% else
+%   tol = 1e-10;
+% end
+% 
+[n,m] = size(d);
+% if ~issym(d,tol),
+%   prwarning(1,'Matrix should be symmetric. It is made symmetric by averaging.')
+%   d = 0.5*(d+d');
+% end
+
+if pzero > n,
+  error('Wrong third parameter.');
+end
+
+%d = (+d).^2;
+
+% if pzero == 0,
+%   % Project the data such that the mean lies at the origin
+%   H = -repmat(1/n,n,n);
+%   H(1:n+1:end) = H(1:n+1:end) + 1; % H = eye(n) - ones(n,n)/n
+%   H = -0.5 * H * d * H;            % H is now the matrix of inner products
+% else
+%   % Project the data such that pzero's object lies at the origin
+%   H = 0.5 * (d(:,pzero) * ones(1,n) + ones(n,1) * d(:,pzero)' - d);
+% end
+H = d;
+H = 0.5*(H+H');                    % Make sure H is symmetric
+
+[Q,L]  = preig(H);
+Q      = real(Q);
+l      = diag(real(L));
+[lm,Z] = sort(-abs(l));
+Q      = Q(:,Z);
+l      = l(Z);                     % Eigenvalues are sorted by decreasing absolute value
+
+[J,sig] = seleigs(l,alf,prec);     % J is the index of the selected eigenvalues
+data.eval = l(J);                          % Eigenvalues
+data.evec = Q(:,J);                        % Eigenvectors
+data.mdis = mean(+d,2)';
+data.mean = pzero;
+data.sig  = sig;
+data.prec = prec;
+
+%A  = Q * diag(sqrt(abs(L)));      % Data in a pseudo-Euclidean space
+
+% Determine the mapping depending on pzero
+if pzero == 0,
+  W = mapping(mfilename,'trained',data,[],m,sum(sig));
+else
+  data.mdis = +d(:,pzero)';
+  W = mapping(mfilename,'trained',data,[],m,sum(sig));
+end
+W = setname(W,'PE embedding');
+return
Index: /distools/parzend_map.m
===================================================================
--- /distools/parzend_map.m	(revision 18)
+++ /distools/parzend_map.m	(revision 18)
@@ -0,0 +1,47 @@
+%PARZEN_MAP Map a dissimilarity dataset on a Parzen densities based classifier
+% 
+% 	F = PARZEN_MAP(D,W)
+%
+% INPUT
+%   D   Dissimilarity dataset
+%   W   Trained Parzen dissimilarity classifier mapping
+%
+% OUTPUT
+%   F   Mapped dataset
+%
+% DESCRIPTION 
+% 
+% SEE ALSO
+% MAPPINGS, DATASETS, PARZENDDC
+
+% Copyright: R.P.W. Duin, r.p.w.duin@prtools.org
+% Faculty EWI, Delft University of Technology
+% P.O. Box 5031, 2600 GA Delft, The Netherlands
+
+function f = parzend_map(d,w)
+	prtrace(mfilename);
+
+	% If no mapping is supplied, train one.
+	if (nargin < 2)
+		w = parzenddc(d); 
+	end
+
+	pars = getdata(w);
+	I = getdata(w,'objects');
+	p = getdata(w,'prior');
+	nlab = getdata(w,'nlab');
+	h = getdata(w,'smoothpar');
+	q = getdata(w,'weights');
+	[k,c] = getsize(w);
+	[m,k] = size(d);
+	e = +d(:,I);
+	
+	f = zeros(m,c);
+	for j=1:c
+		J = find(nlab==j);
+		f(:,j) = p(j)*mean(repmat(q(J)',m,1).*exp(-(e(:,J).^2)./repmat((2.*h(J).*h(J))',m,1)),2);
+	end
+	
+	f = setdata(d,f,getlabels(w));
+	
+return
Index: /distools/parzenddc.m
===================================================================
--- /distools/parzenddc.m	(revision 18)
+++ /distools/parzenddc.m	(revision 18)
@@ -0,0 +1,152 @@
+%PARZENDDC Parzen density based classifier for dissimilarity data
+% 
+%  [W,H] = PARZENDDC(D)
+%  W     = PARZENDDC(D,H)
+% 
+% INPUT
+%   D   Dissimilarity dataset
+%   H   Smoothing parameters (optional; default: estimated from D)
+%
+% OUTPUT
+%  W    Trained Parzen classifier
+%  H    Smoothing parameters, estimated from the data
+%
+% DESCRIPTION
+%
+% SEE ALSO
+% DATASETS, MAPPINGS, PARZEND_MAP
+ 
+% Copyright: R.P.W. Duin, r.p.w.duin@prtools.org
+% Faculty EWI, Delft University of Technology
+% P.O. Box 5031, 2600 GA Delft, The Netherlands
+
+function [W,h] = parzenddc(d,h)
+  
+  prtrace(mfilename);
+  
+  if nargin < 2, h = []; end
+  
+  % No input arguments: return an untrained mapping.
+  if nargin == 0 | isempty(d)
+    W = mapping(mfilename,h); 
+    W = setname(W,'Parzen DisRep Classifier');
+    return; 
+  end
+  
+  if nargin < 2 | isempty(h);
+    %prwarning(5,'Smoothing parameters not specified, estimated from the data.');
+    %error('Smoothing parameter estimation not yet implemented')
+    h = [];
+  end
+
+  islabtype(d,'crisp','soft');
+  isvaldfile(d,2,2); % at least 2 objects per class, 2 classes
+  d = testdatasize(d);
+  d = testdatasize(d,'objects');
+  
+  [m,k,c] = getsize(d);
+  nlab = getnlab(d);
+
+  if ~isempty(h)       % Take user settings for smoothing parameters.
+    
+    nlab = getnlab(d);
+    if length(h) == 1
+      hh = repmat(h,m,1);
+    elseif length(h) == c
+      hh = zeros(m,1);
+      for j=1:c
+        J = find(nlab==j);
+        hh(J) = h(j);
+      end
+    else
+      error('Smoothing parameter has wrong size')
+    end
+
+  else   % Estimate smoothing parameters (copied and adapted from parzenc)
+  
+    %error('Not yet implemented')
+
+    % compute all object distances
+    D = +d;
+    dim = intrdim(D);
+    hinit = max(D(:)); % for sure a too high value
+    D = D.^2 + diag(inf*ones(1,m));
+    % find object weights q
+    q = classsizes(d);
+  
+    % find for each object its class freqency
+    of = q(nlab);
+  
+    % find object weights q
+    p = getprior(d);
+    %a = setprior(a,p);
+    q = p(nlab)./q(nlab);
+  
+    % initialise
+    h = hinit;
+    %h = 2.5
+    L = -inf;
+    Ln = 0;
+    z = 0.01^(1/dim); % initial step size
+
+    % iterate
+  
+    len1 = prprogress([],'parzenddc: error optimization smoothing parameter: ');
+    len2 = prprogress([],' %6.4f  %6.4f ',0,0);
+	  iter = 0;
+    prwaitbar(100,'parzenc: Optimizing smoothing parameter',m > 100);
+    while abs(Ln-L) > 0.0001 & z < 1
+      % In L we store the best performance estimate found so far.
+      % Ln is the actual performance (for the actual h)
+      % If Ln > L we improve the bound L, and so we rest it.
+		  iter = iter+1;
+		  prwaitbar(100,100-100*exp(-iter/10));
+      if Ln > L, L = Ln; end
+
+      r = -0.5/(h^2);
+      F = q(ones(1,m),:)'.*exp(D*r);           % density contributions
+      FS = sum(F)*((m-1)/m); IFS = find(FS>0); % joint density distribution
+      if islabtype(d,'crisp');
+        G = sum(F .* (nlab(:,ones(1,m)) == nlab(:,ones(1,m))'));
+      else
+        G = zeros(1,m);
+        for j=1:c
+          G = G + sum(F .* (d.targets(:,j) * d.targets(:,j)'));
+        end
+      end
+      G = G.*(of-1)./of;                       % true-class densities
+      % performance estimate, neglect zeros
+    
+      en = max(p)*ones(1,m);
+      en(IFS) = (G(IFS))./FS(IFS);
+      Ln = exp(sum(log(en+realmin))/m);
+  
+      closemess([],len2);
+      len2 = prprogress([],' %6.4f  %6.4f ',h,Ln);
+
+      if Ln < L            % compute next estimate
+        z = sqrt(z);       % adjust stepsize up (recall: 0 < z < 1)
+        h = h / z;         % if we don't improve, increase h (approach h_opt from below)
+      else
+        z = z*z;
+        h = h * z;         % if we improve, decrease h (approach h_opt from above)
+      end
+      %disp([Ln,L,z,h])
+    end
+    closemess([],len1+len2);
+    hh = h*ones(m,1);
+    prwaitbar(0);
+
+  end
+
+  par.objects = [1:m]';
+  par.prior   = getprior(d);
+  par.nlab    = getnlab(d);
+  par.smoothpar = hh;
+  par.weights = ones(m,1);
+    
+  W = mapping('parzend_map','trained',par,getlablist(d),k,c);
+  W = setname(W,'Parzen DisRep Classifier');
+  W = setcost(W,d);
+
+return;
