Index: distools/private/auclpm.m
===================================================================
--- distools/private/auclpm.m	(revision 13)
+++ distools/private/auclpm.m	(revision 13)
@@ -0,0 +1,234 @@
+function [w,returnA,r2] = auclpm(x, C, rtype, par, unitnorm, usematlab)
+%AUCLPM Find linear mapping with optimized AUC
+%
+%    W = AUCLPM(X, C, RTYPE, PAR)
+%
+% Optimize the AUC on dataset X and reg. param. C. This is done by
+% finding the weights W for which the ordering of the objects mapped
+% onto the line defined by W, is optimal. That means that objects from
+% class +1 is always mapped above objects from the -1 class. This
+% results in a constraint for each (+1,-1) pair of objects. The number
+% of constraints therefore become very large.  The AUC constraints can
+% be subsampled in different ways:
+%
+%		RTYPE     PAR
+%	  'full',     -   use all constraints
+%	  'subs',     N   subsample just N constraints
+%	  'subk',     k   subsample just k*#trainobj. constraints
+%	  'knn'       k   use only the k nearest neighbors
+%	  'xval'      N   subsample just N constraints and use the rest to
+%                 optimize C (this version can be very slow)
+%    'xvalk'     k   subsample k*#trainobj and use remaining constraints
+%                 to optimize C
+%    'kmeans'    k   use k-means clustering with k=PAR
+%    'randk'     subsample objects to get PAR*(Npos+Nneg) constraints
+%
+%    W = AUCLPM(X, C, RTYPE, PAR, UNITNORM)
+%
+% Finally, per default the difference vectors are normalized to unit
+% length. If you don't like that, set UNITNORM to 0.
+%
+% Default: RTYPE='subk'
+%          PAR  = 1.0;
+prtrace(mfilename);
+
+if (nargin < 6)
+	usematlab = 0;
+end
+if (nargin < 5)
+	unitnorm = 0;
+end
+if (nargin < 4)
+	par = 1.00;
+end
+if (nargin < 3)
+	rtype = 'subk';
+end
+if (nargin < 2)
+	prwarning(3,'Lambda set to ten');
+	C = 10; 
+end
+
+if (nargin < 1) | (isempty(x))
+	% just a check...
+	if ~isa(C,'double')
+		error('Please check your input parameters: C should be a double.');
+	end
+	w = mapping(mfilename,{C,rtype,par,unitnorm,usematlab});
+	w = setname(w,defauclpmname(C,rtype,par,unitnorm));
+	return
+end
+
+if ~ismapping(C)   % train the mapping
+
+	% just a check...
+	if ~isa(C,'double')
+		error('Please check your input parameters: C should be a double.');
+	end
+	% Unpack the dataset.
+	islabtype(x,'crisp');
+	isvaldset(x,1,2); % at least 1 object per class, 2 classes
+	[n,k,c] = getsize(x); 
+	% Check some values:
+	if par<=0
+		error('Parameter ''par'' should be larger than zero');
+	end
+
+	if c == 2  % two-class classifier:
+
+		labl = getlablist(x); dim = size(x,2);
+		% first create the target values (+1 and -1):
+		% make an exception for a one-class or mil dataset:
+		tnr = strmatch('target',labl);
+		if isempty(tnr) % no target class defined.
+			tnr = strmatch('positive',labl);
+			if isempty(tnr) % no positive class defined.
+				% we just take the first class as target class:
+				tnr = 1;
+			end
+		end
+		y = 2*(getnlab(x)==tnr) - 1;
+		tlab = labl(tnr,:);
+
+		% makes the mapping much faster:
+		X = +x; clear x;
+
+		%---create A for optauc
+        rstate = rand('state');
+		seed = 0;
+		[A,Nxi,Aval] = createA(X,y,rtype,par,seed);
+        rand('state',rstate);
+		if unitnorm
+			% normalize the length of A:
+			lA = sqrt(sum(A.*A,2));
+			lenn0 = find(lA~=0);  % when labels are flipped, terrible
+			                      % things can happen
+			A(lenn0,:) = A(lenn0,:)./repmat(lA(lenn0,:),1,size(A,2));
+			if ~isempty(Aval)
+				% also normalize the length of Aval:
+				lA = sqrt(sum(Aval.*Aval,2));
+				lenn0 = find(lA~=0);
+				Aval(lenn0,:) = Aval(lenn0,:)./repmat(lA(lenn0,:),1,size(Aval,2));
+			end
+		end
+orgA = A;
+		% negative should be present for the constraints:
+		A = [A -A];
+		% take also care for the xi:
+		A = [A -speye(Nxi)];
+		%A = [A -eye(Nxi)];
+		%---create f
+		% NO, do this later, maybe we want to optimize it!
+		%f = [ones(2*k,1); repmat(C,Nxi,1)];
+		%--- generate b
+		b = -ones(Nxi,1);   % no zeros, otherwise we get w=0
+		    % the constraint is changed here to <=-1
+		%---lower bound constraints
+		lb = zeros(2*k+Nxi,1);
+
+		% should we run over a range of Cs?
+		if ~isempty(Aval)
+			M = 25;
+			xtr = zeros(M,1);
+			xval = zeros(M,1);
+			C = logspace(-3,3,M);
+         % run over all the Cs
+			for i=1:length(C)
+				%---create f again:
+				f = [ones(2*k,1); repmat(C(i),Nxi,1)];
+				%---solve linear program
+				if (exist('glpk')>0) & ~usematlab
+					[z,fmin,status]=glpk(f,A,b,lb,[],repmat('U',Nxi,1),...
+						repmat('C',size(f,1),1),1);
+				elseif (exist('glpkmex')>0) & ~usematlab
+					[z,fmin,status]=glpkmex(1,f,A,b,repmat('U',Nxi,1),...
+						lb,[],repmat('C',size(f,1),1));
+				else
+					opts = optimset('Display','off','LargeScale','on','Diagnostics','off');
+					z = linprog(f,A,b,[],[],lb,[],[],opts);
+				end
+				constr = Aval*(z(1:k)-z(k+1:2*k));
+				% the number of satisfied constraints (=AUC:)
+				I = find(constr<-0);
+				if ~isempty(I)
+					xval(i) = length(I)/size(constr,1);
+				end
+				% the number of satisfied constraints (=AUC:)
+				constr = orgA*(z(1:k)-z(k+1:2*k));
+				I = find(constr<-0);
+				if ~isempty(I)
+					xtr(i) = length(I)/size(constr,1);
+				end
+			end
+if nargout>1
+	returnA = xval;
+	if nargout>2
+		r2 = xtr;
+	end
+end
+			[minxval,mini] = max(xval);
+			C = C(mini);
+		end
+		%---create f
+		f = [ones(2*k,1); repmat(C,Nxi,1)];
+		%---solve linear program
+		if (exist('glpkmex')>0) & ~usematlab
+			prwarning(7,'Use glpkmex');
+			param.msglev=0;
+			[z,fmin,status,xtra]=glpkmex(1,f,A,b,repmat('U',Nxi,1),...
+				lb,[],repmat('C',size(f,1),1),param);
+			alpha = []; %xtra.lambda;
+		else
+			[z,fmin,exitflag,outp,alpha] = linprog(f,A,b,[],[],lb);
+		end
+
+		%---extract parameters
+		u = z(1:k); u = u(:);
+		v = z(k+1:2*k); v = v(:);
+		zeta = z(2*k+1:2*k+Nxi); zeta = zeta(:);
+	else
+		error('Only a two-class classifier is implemented');
+	end
+	% now find out how sparse the result is:
+	rel = (abs(u-v)>1e-6);
+	nr = sum(rel);
+	if (nr==0)
+		error('None of the features is selected. Please make the C a bit larger.');
+	end
+	
+	% and store the results:
+	W.u = u-v; %the ultimate weights
+	W.alpha = alpha;
+	W.zeta = zeta;
+	W.nr = nr;
+	W.rel = rel;
+	W.C = C;
+	w = mapping(mfilename,'trained',W,tlab,dim,1);
+	w = setname(w,defauclpmname(C,rtype,par,unitnorm));
+	
+else
+	% Evaluate the classifier on new data:
+	W = getdata(C);
+	n = size(x,1);
+
+	% linear classifier:
+	out = x*W.u;
+
+	% and put it nicely in a prtools dataset:
+	% (I am not really sure what I should output, I decided to give a 1D
+	% output:)
+	w = setdat(x,out,C);
+
+end
+		
+return
+
+function cl_name = defauclpmname(C,rtype,par,unitnorm)
+% define the correct name:
+if unitnorm
+	cl_name = sprintf('AUClpm (C=%s, %s, k=%s) 1norm',num2str(C),rtype,num2str(par));
+else
+	cl_name = sprintf('AUClpm (C=%s, %s, k=%s)',num2str(C),rtype,num2str(par));
+end
+
+return
Index: distools/private/bord.m
===================================================================
--- distools/private/bord.m	(revision 13)
+++ distools/private/bord.m	(revision 13)
@@ -0,0 +1,22 @@
+% C = bord(A,n,m)
+% Puts a border of width m (default m=1) around image A
+% and gives it value n. If n = NaN: mirror image values.
+function C = bord(A,n,m);
+if nargin == 2; m=1; end
+[x,y] = size(A);
+if m > min(x,y)
+        mm = min(x,y);
+        C = bord(A,n,mm);
+        C = bord(C,n,m-mm);
+        return
+end
+if isnan(n)
+   C = [A(:,m:-1:1),A,A(:,y:-1:y-m+1)];
+   C = [C(m:-1:1,:);C;C(x:-1:x-m+1,:)];
+else
+   bx = ones(x,m)*n;
+   by = ones(m,y+2*m)*n;
+   C = [by;[bx,A,bx];by];
+end
+return
+
Index: distools/private/cuteigs.m
===================================================================
--- distools/private/cuteigs.m	(revision 13)
+++ distools/private/cuteigs.m	(revision 13)
@@ -0,0 +1,108 @@
+%CUTEIGS Select significant eigenvalues from a list
+%
+%   J = CUTEIGS(L,PREC)
+%
+% INPUT
+%   L     List of eigenvalues
+%   PREC  Precision parameter (optional; default: 0.0001)
+%
+% OUTPUT
+%   J   Index of selected eigenvalues
+%
+% DESCRIPTION
+% This is a low-level routine for SELEIGS, which serves PSEM, KPSEM, KPCA and
+% PEPCA. It makes use of PCHIP determing piecewise cubic Hermite interpolating
+% polynomial P, built on [1:length(LL) LL], where LL = -sort(-abs(L).^3) folloowed
+% by a normalization to [0,1]. A third power is used to emphasize the differences
+% better, as sometimes there is a very long tail of slowly dropping eigenvalues.
+% The number M of significant eigenvalues is determined such that the ratio
+% of the estimated derivative of P versus the largest derivative value is smaller
+% then the given precision, PREC * cumsum(LL).
+%
+% Note that M will differ between different samples and for different sizes of L.
+%
+% SEE ALSO
+% MAPPINGS, DATASETS, PSEM, KPSEM, KPCA, PEPCA
+%
+
+% Elzbieta Pekalska, e.pekalska@ewi.tudelft.nl
+% Faculty of Electrical Engineering, Mathematics and Computer Science,
+% Delft University of Technology, The Netherlands
+
+
+function [J,prec] = cuteigs(L,prec)
+
+if nargin < 2 | isempty(prec),
+  prec = 0.0001;
+end
+
+tol = 1e-14;
+
+n  = length(L);
+LL = L.^3;
+[lambda,J] = sort(-abs(LL)/sum(abs(LL)));
+lambda = -lambda;
+
+if n > 1,
+  M = findJ(lambda,prec);
+else
+  M = n;
+end
+JJ = J(1:M);
+
+I1 = find (L(JJ) >= 0);
+I2 = find (L(JJ) < 0);
+
+J  = [JJ(I1); JJ(I2)];
+
+
+
+
+
+
+
+function M = findJ (lambda,prec)
+
+n = length(lambda);
+if n >= 400,
+  Z = [1:80 81:4:n];
+elseif n >= 200,
+  Z = [1:50 51:2:n];
+else
+  Z =1:n;
+end
+
+if n > 3
+  % Hermite piecewise cubic interpolating polynomial
+  H = pchip(Z,lambda(Z));
+
+  % Derivative of the Hermite piecewise cubic polynomial
+  Hder = H;
+  Hder.order = 3;
+  Hder.coefs = zeros(H.pieces,3);
+  for i=1:Hder.pieces
+    pp = polyder(H.coefs(i,:));
+    Hder.coefs(i,4-length(pp):3) = pp;
+  end
+  pder= -ppval(1:n,Hder);   % Evaluated derivative
+else
+  pder= [lambda(1) -diff(lambda')];   % Approximated derivative
+end
+
+% pder(1) is an extrapolated value; pder(2) is used
+I  = find(pder/pder(2) < prec*cumsum(lambda'));
+II = diff(I);
+Q  = find(II > 1);
+if ~isempty(Q),
+  I = I(Q(1));
+end
+if ~isempty(I),
+  M = I(1)-1;
+else
+  if n <= 25,
+    M = n;
+  else
+    M = n-1;
+  end
+end
+if M < 1, M = 1; end
Index: distools/private/invsig.m
===================================================================
--- distools/private/invsig.m	(revision 13)
+++ distools/private/invsig.m	(revision 13)
@@ -0,0 +1,4 @@
+%
+function a = invsig(a)
+a = log(a+realmin) - log(1-a+realmin);
+return
Index: distools/private/seleigs.m
===================================================================
--- distools/private/seleigs.m	(revision 13)
+++ distools/private/seleigs.m	(revision 13)
@@ -0,0 +1,179 @@
+%SELEIGS Select eigenvalues from a list
+%
+%   J = SELEIGS(L,ALF)
+%
+% INPUT
+%   L   List of eigenvalues
+%   ALF Parameter determining the dimensionality and the eigenvalue-based mapping
+%       (0,1)   - fraction of the total (absolute value) preserved variance
+%        Inf    - no dimensionality reduction, keeping all dimensions (it's 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]
+%
+% OUTPUT
+%   J   Index of selected eigenvalues
+%
+% DESCRIPTION
+% This is a low-level routine for PSEM, KPSEM, KPCA and PEPCA. From a list of
+% eigenvalues it selects the ones according to ALF.
+%
+% SEE ALSO
+% MAPPINGS, DATASETS, PSEM, KPSEM, KPCA, PEPCA
+%
+
+% Elzbieta Pekalska, e.pekalska@ewi.tudelft.nl
+% Faculty of Electrical Engineering, Mathematics and Computer Science,
+% Delft University of Technology, The Netherlands
+
+
+function [J,sig] = seleigs(L,alf,prec)
+
+alfstr = lower(alf);
+
+if strcmp(alfstr,'cut'),
+  J = cuteigs(L,prec);
+else
+  % Check whether alf is a string
+  if ~isnumeric(alfstr),
+    zp = strfind(alfstr,'p');
+    zn = strfind(alfstr,'n');
+    if isempty(zp), zp = 0; end
+    if isempty(zn), zn = 0; end
+    ALF = 0;
+    if zp > 0,
+      if length(zp) > 1,
+        error('Wrong ALF.');
+      end
+      if zn > zp,
+        sstr = alfstr(1:zp-1);
+      elseif zn < zp,
+        sstr = alfstr(zn+1:zp-1);
+      else
+        error('Wrong ALF.');
+      end
+      if ~isempty(sstr),
+        ALF = str2num(sstr);
+      else
+        ALF = length(find(L>0));
+      end
+    end
+    if zn > 0,
+      if length(zn) > 1,
+        error('Wrong ALF.');
+      end
+      if zn > zp,
+        sstr = alfstr(zp+1:zn-1);
+      elseif zn < zp,
+        sstr = alfstr(1:zn-1);
+      else
+        error('Wrong ALF.');
+      end
+      if ~isempty(sstr),
+        ALF = [ALF str2num(sstr)];
+      else
+        ALF = [ALF length(find(L<0))];
+      end
+    else
+      ALF = [ALF 0];
+    end
+    alf = ALF;
+  end
+
+
+% alf is now the number of eigenvalues
+  if max(size(alf)) > 2 | (alf(1) == 0 & alf(2) == 0),
+    error('Wrong ALF.')
+  elseif max(size(alf)) == 2,
+    J = [];
+    Z  = find(L>0);
+    if alf(1) > 0,
+      if ~isempty(Z),
+        J = selnumalf(alf(1),L,Z)';
+      else
+        prwarning(1,'No positive dimensions are selected.');
+      end
+    end
+    Z  = find(L<0);
+    if alf(2) > 0,
+      if ~isempty(Z),
+        J = [J; selnumalf(alf(2),L,Z)'];
+      else
+        prwarning(1,'No negative dimensions are selected.');
+      end
+    end
+   else
+     J = selnumalf(alf,L,1:length(L));
+  end
+end
+
+
+if isempty(J),
+  error('Wrong choice of ALF; the subspace cannot be found.');
+end
+
+I1  = find(L(J) > 0);
+I2  = find(L(J) < 0);
+sig = [length(I1) length(I2)];
+
+[ll,K1] = sort(-L(J(I1)));
+[ll,K2] = sort(L(J(I2)));
+
+% J - index of selected eigenvalues; SORTED:
+% first decreasing positive ones folowed by increasing negative ones
+JJ = [];
+if sig(1) > 0
+  JJ = J(I1(K1));
+end
+if sig(2) > 0
+  JJ = [JJ; J(I2(K2))];
+end
+J = JJ;
+return
+
+
+
+
+
+function J = selnumalf(alf,L,I)
+
+if alf >= 1,
+  if alf == inf,
+    cs = cumsum(abs(L(I)));
+    nn = min(find(cs./cs(end) == 1));
+    J  = I(1:nn)';
+  elseif alf > length(I),
+    error('ALF is invalid.')
+  else
+    J = I(1:alf)';
+  end
+elseif alf > 0,                     % alf in (0,1), percentage
+  cs = cumsum(abs(L(I)));
+  nn = min(find (cs./cs(end) >= alf));
+  J  = I(1:nn)';
+elseif alf < 0,
+  I   = I(length(I):-1:1);
+  alf = abs(alf);
+  if alf >= 1,
+    if alf > length(I),
+      error('ALF is invalid.')
+    else
+      J  = I(1:alf)';
+    end
+  else                             % |alf| in (0,1), percentage
+    cs = cumsum(abs(L(I)));
+    nn = min(find (cs./cs(end) >= alf));
+    J  = I(1:nn)';
+  end
+else
+  error('Wrong ALF.');
+end
Index: distools/private/updname.m
===================================================================
--- distools/private/updname.m	(revision 13)
+++ distools/private/updname.m	(revision 13)
@@ -0,0 +1,39 @@
+%UPDNAME Update of a distance/similarity dataset
+%
+%   NAME = UPDNAME (NAME)
+%
+% A low level routine used in other functions
+
+% Copyright: Elzbieta Pekalska, e.pekalska@ewi.tudelft.nl
+% Faculty EWI, Delft University of Technology
+% P.O. Box 5031, 2600 GA Delft, The Netherlands
+
+
+
+function name = updname(name)
+
+str = lower(name);
+NAMES = {'distance matrix', 'similarity matrix', 'kernel matrix'};
+
+ok = 0;
+z  = 1;
+while ~ok & z<=length(NAMES)
+  NAME = NAMES{z};
+  m = length(NAME);
+  k = findstr(str,NAME);
+  if ~isempty(k),
+    if length(str) >=k+m+3 & strcmp(str(k+m+1:k+m+3),'for'),
+      name(k:k+m+3)=[];
+    else
+      name(k:k+m-1)=[];
+    end
+    if ~isempty(name) & strcmp(name(1), ' '),
+      name(1)= [];
+    end
+    ok = 1;
+  end
+  z = z+1;
+end
+if isempty(name),
+  name = 'Data';
+end
