Index: /prextra/prextra/Contents.m
===================================================================
--- /prextra/prextra/Contents.m	(revision 2)
+++ /prextra/prextra/Contents.m	(revision 2)
@@ -0,0 +1,55 @@
+% PR Extra repository
+% Version 0.3 6-Jan-2005
+%
+%Data generation
+%---------------
+%genriver     Generate 'bending river'-like 2D test data
+%gendati      Generate random image windows
+%
+%Mappings
+%--------
+%unitm        Fixed unit mapping
+%unittm       Trainable unit mapping
+%
+%Classifiers
+%-----------
+%incsvc       Incremental support vector classifier
+%gatem        Gateway mapping
+%emparzenc    EM Parzen Classifier for semi-supervised learning
+%emc          EM algorithm for semi-supervised learning
+%linprogc     Linear Program Classifier (includes sparse solutions)
+%fsellpc      Low level routine for linprogc
+%lessc        Least Error in Sparse Subspaces classifier (David Tax)
+%adaboostc    Adaboost classifier
+%rbsvc        Automatic radial basis SVC
+%
+%auc          error under the curve estimator
+%
+%Image manipulation
+%------------------
+%impaint      Create and paint an label overlay on an image
+%impatch      Create and change a polygon label overlay on an image
+%neighbh      Compute the indices of neighboring pixels
+%
+%Feature Selection
+%-----------------
+%featsel2     Pairwise feature selection
+%
+%ilab: user-defined named identifiers
+%------------------------------------
+%enableilab   enable named identifiers in a dataset (ilab)
+%isvalidilab  validity check for named-identifier support in a dataset
+%getilablist  retrieve list of named identifiers
+%addilab      add named identifier to a dataset
+%getilab      retrieve identifier values given identifier names or indices
+%setilab      set the values for the per-example named-identifiers in a dataset
+%removeilab   remove named identifier columns
+%
+%
+%meta: named meta data in the dataset user field
+%-----------------------------------------------
+%enablemeta   enable the storage of key-value pairs in a user field
+%getmeta      retrieve metadata in a dataset by key names
+%setmeta      add/update key-value pair(s) to a dataset
+%removemeta   remove key-value pairs
+
Index: /prextra/prextra/add_obj_cl.m
===================================================================
--- /prextra/prextra/add_obj_cl.m	(revision 2)
+++ /prextra/prextra/add_obj_cl.m	(revision 2)
@@ -0,0 +1,275 @@
+%Auxiliary script for incsvc.
+%I especially did not make it a function to avoid a ton of arguments to
+%pass. Maybe it gives a small speedup...
+
+%fprintf('Add object number %d\n',c);
+% currently seen objects plus new object
+newsetD = (1:c)';
+
+% compute the kernel matrix entry for the new object:
+% (let us assume this is a row vector)
+K = feval(kernel,kpar,c,newsetD);
+Kc = (y(c)*y(newsetD)').*K;
+
+%  % first see if the object is accepted or not
+%  % (note that I could also have looked at the gradient, but then I
+%  % don't have a very honest 'distance to the boundary')
+%  this_label = K(c) - b;
+%  if ~isempty(setS)
+%    this_label = this_label - K(setS)*(y(setS).*alf(setS));
+%  end
+%  if ~isempty(setE)
+%    this_label = this_label - K(setE)*(y(setE).*alf(setE));
+%  end
+
+% Compute the gradient for the new object:
+grad(c,1) = y(c)*b - 1; %DXD here is a change for SVC
+if ~isempty(setS)
+	grad(c,1) = grad(c,1) + Kc(setS)*alf(setS);
+end
+if ~isempty(setE)
+	grad(c,1) = grad(c,1) + Kc(setE)*alf(setE);
+end
+
+% Right, here we go, can be add object c?:
+if grad(c)>0 % it is already classified ok,
+	if length(setR)>0
+		setR = [setR;c];   % put in the 'rest' set
+		Kr = [Kr; Kc(setS)];  % update Kr
+	else
+		setR = c;
+		Kr = Kc(setS);
+	end
+
+else         % we have to work
+
+	done = 0;  % to check if we are done
+	nrloops=0;
+	while ~done
+
+		% compute beta, not for the new object:
+		beta = zeros(length(setS)+1,1);
+		beta = -R*[y(c); Kc(setS)'];
+
+		if any(beta>1e100) & ~isempty(setS)
+			disp('Beta is too large: do you have identical objects with different labels??'); keyboard;
+		end
+
+		% compute gamma, also for the new object:
+		if isempty(setS)
+			% here something fishy is going on, when there is just a single
+			% object added. Then we cannot freely move alf, and we can only
+			% move b. In this case, b is moved until alf(c) enters setS
+			% (that means, the gradient becomes 0).
+			gamma = y(c)*y(newsetD);
+			duptonow = y(c)*b;
+		else
+			gamma = Kc';
+			if ~isempty(setE)
+				gamma(setE) = gamma(setE) + [y(setE) Ke]*beta;
+			end
+			if ~isempty(setR)
+				gamma(setR) = gamma(setR) + [y(setR) Kr]*beta;
+			end
+			gamma(c) = gamma(c) + [y(c) Kc(setS)]*beta;
+			gamma(setS) = 0;
+			duptonow = alf(c);
+		end
+
+      % now we have to see how large deltaAc can become...
+      % (1) check the own upper bound:
+      if isempty(setS)
+        deltaAcisC = inf; %because we're moving b, and not alf!
+      else
+        deltaAcisC = C;
+      end
+      % (2) check if own gradient becomes zero:
+      warning off;
+      deltaGcis0 = duptonow-grad(c)./gamma(c);
+      warning on;
+      % object moves the wrong way:
+      deltaGcis0(deltaGcis0<duptonow) = inf;
+      % (3) check upper bounds of the SVs:
+      deltaupperC = inf;
+      if ~isempty(setS)
+        warning off;
+        deltaupperC = duptonow + (C-alf(setS))./beta(2:end);
+        warning on;
+        % negative changes do not count:
+        deltaupperC(deltaupperC<duptonow) = inf;
+        % object moves the wrong way (or not at all):
+        deltaupperC(beta(2:end)<=0) = inf;
+        [deltaupperC,nrS_up] = min(deltaupperC);
+      end
+      % (4) check lower bounds of the SVs:
+      deltalowerC = inf;
+      if ~isempty(setS)
+        warning off;
+        deltalowerC = duptonow + -alf(setS)./beta(2:end);
+        warning on;
+        % negative changes do not count:
+        deltalowerC(deltalowerC<=duptonow) = inf;
+        % object moves the wrong way (or not at all)
+        deltalowerC(beta(2:end)>=0) = inf;
+        [deltalowerC,nrS_low] = min(deltalowerC);
+      end
+      % (5) check E gradients to become 0:
+      deltaGeis0 = inf;
+      if ~isempty(setE)
+        warning off; % divide by 0 is taken care of later...
+        deltaGeis0 = duptonow -grad(setE)./gamma(setE);
+        warning on;
+        deltaGeis0(deltaGeis0<duptonow) = inf;
+        deltaGeis0(gamma(setE)<=0) = inf;
+        [deltaGeis0,nrE_0] = min(deltaGeis0);
+      end
+      % (6) check R gradients to become 0:
+      deltaGris0 = inf;
+      if ~isempty(setR)
+        warning off; % divide by 0 is taken care of later...
+        deltaGris0 = duptonow -grad(setR)./gamma(setR);
+        warning on;
+        %deltaGris0(deltaGris0<duptonow) = inf;
+        deltaGris0(deltaGris0<=duptonow) = inf;
+        deltaGris0(gamma(setR)>=0) = inf;
+        [deltaGris0,nrG_0] = min(deltaGris0);
+      end
+
+      % which one is the most urgent one?
+      deltas = [deltaAcisC; deltaGcis0; deltaupperC; deltalowerC;...
+                deltaGeis0; deltaGris0];
+      [maxdelta,situation] = min(deltas);
+% fprintf('Situation %d (max_delta=%f)\n',situation,maxdelta);
+%
+      % update the parameters
+      if isempty(setS) % then we only change b:
+        %disp('setS is empty!');
+        b = y(c)*maxdelta;
+      else
+        alf(c) = maxdelta;
+        alf(setS) = alf(setS) + (maxdelta-duptonow)*beta(2:end);
+        b = b + (maxdelta-duptonow)*beta(1);
+%fprintf(' b: %f -> %f\n',bold,b);
+%alf'
+%fprintf('after: sum alf = %f\n',sum(alf));
+      end
+      grad = grad + (maxdelta-duptonow)*gamma;
+      
+      % do consistency check:
+      I = find(alf<tol);
+      if ~isempty(I)
+        J = find(alf(I)<-tol);
+%        if ~isempty(J)
+%          disp('One of the alpha''s became < 0!');
+%			 alf(I(J))
+%          keyboard;
+%		  end
+        alf(I) = 0;
+      end
+
+      % update the sets:
+      %fprintf('%d: situation %d  \n',c,situation);
+      switch situation
+      case 1   % object c goes to setE
+        alf(c) = C; % just to be sure...
+        if size(Ke,1)==0, Ke = []; end % make it really empty
+        Ke = [Ke; Kc(setS)];
+        setE = [setE; c];
+        done = 1;
+      case 2   % object c goes to setS
+        Ks = [Ks [y(c); Kc(setS)']; y(c) Kc([setS; c])];
+        Ke = [Ke Kc(setE)'];
+        Kr = [Kr Kc(setR)'];
+        % update R 
+        if isempty(setS) % compute it directly (to avoid the inf's)...
+          R = [-Kc(c) y(c); y(c) 0];
+        else
+          R = change_R(R,+c,beta,gamma(c));
+        end
+        setS = [setS;c];
+        done = 1;
+      case 3  % a support object hits upper bound
+        j = setS(nrS_up);             % number of the object
+        alf(j) = C;                   % just to be sure
+        if size(Ke,1)==0, Ke=[]; end  % make it really really empty
+        Ke = [Ke; Ks(nrS_up+1,2:end)];  % update Ke
+        setE = [setE;j];              % add to setE
+        Ks(nrS_up+1,:) = [];            % update all K's
+        Ks(:,nrS_up+1) = [];
+        Ke(:,nrS_up) = [];
+        if ~isempty(Kr), Kr(:,nrS_up) = []; end
+        setS(nrS_up) = [];            % remove from setS
+        R = change_R(R,-nrS_up,beta,gamma(j));
+      case 4  % a support object hits lower bound
+        j = setS(nrS_low);             % number of the object
+        alf(j) = 0;                    % just to be sure
+        if size(Kr,1)==0, Kr = []; end % make really empty
+        Kr = [Kr; Ks(nrS_low+1,2:end)];  % update Kr
+        setR = [setR;j];               % add to setE
+        Ks(nrS_low+1,:) = [];            % update all K's
+        Ks(:,nrS_low+1) = [];
+        if ~isempty(Ke), Ke(:,nrS_low) = []; end;
+        if ~isempty(Kr), Kr(:,nrS_low) = []; end;
+        setS(nrS_low) = [];            % remove from setS
+        R = change_R(R,-nrS_low,beta,gamma(j));
+      case 5  % an error becomes a support object
+        j = setE(nrE_0);              % number of the object
+        % adding to setS, means that all kernels have to be computed:
+        K = feval(kernel,kpar,j,newsetD);
+        Kj = (y(j)*y(newsetD)').*K;
+        % to update R, we have to have the beta of object j:
+        betaj = zeros(length(setS)+1,1);
+        betaj = -R*[y(j); Kj(setS)'];
+        Ks = [Ks; y(j) Kj(setS)];     % add row to Ks
+        Kr = [Kr Kj(setR)'];          % update Kr
+        Ke = [Ke Kj(setE)'];          % update Ke
+        Ke(nrE_0,:) = [];
+        setE(nrE_0) = [];             % update setE
+        setS = [setS;j];              % add to setS
+        Ks = [Ks [y(j); Kj(setS)']];   % and the extra column for Ks
+        if length(betaj)==1 % compute it directly (to avoid the inf's)...
+          R = [-Kj(j) y(j); y(j) 0];
+        else
+          % to update R, we also have to have the gamma of object j:
+          gammaj = Ks(end,:)*[betaj;1] ;
+          R = change_R(R,+j,betaj,gammaj);
+        end
+      case 6  % an other object becomes a support object
+        j = setR(nrG_0);              % number of the object
+        % adding to setS, means that all kernels have to be computed:
+        K = feval(kernel,kpar,j,newsetD);
+        Kj = (y(j)*y(newsetD)').*K;
+        % to update R, we have to have the beta of object j:
+        betaj = zeros(length(setS)+1,1);
+        betaj = -R*[y(j); Kj(setS)'];
+        Ks = [Ks; y(j) Kj(setS)];     % add row to Ks
+        Ks = [Ks [y(j); Kj([setS;j])']]; % and the extra column for Ks
+        Ke = [Ke Kj(setE)'];          % update Ke
+        Kr = [Kr Kj(setR)'];          % update Kr
+        Kr(nrG_0,:) = [];
+        setS = [setS;j];              % add to setS
+        setR(nrG_0) = [];             % update setR
+        if length(betaj)==1 % compute it directly (to avoid the inf's)...
+          R = [-Kj(j) y(j); y(j) 0];
+        else
+          % to update R, we also have to have the gamma of object j:
+          gammaj = Ks(end,:)*[betaj;1] ;
+          R = change_R(R,+j,betaj,gammaj);
+        end
+      end
+      %fprintf('end situation\n');
+      %alf,b
+		nrloops = nrloops+1;
+		%fprintf('nrloops = %d\n',nrloops);
+		if nrloops>50, done=1; end
+      %keyboard
+    end % changing alfa's till stable solution is found
+end % check for gradient<=0
+
+
+% now we can also compute the R, compute the output for the x(setS):
+% But fortunately, this value is
+%   R^2 = offset + b;
+% If we ignore the offset, we just have to use b!
+
+
Index: /prextra/prextra/addilab.m
===================================================================
--- /prextra/prextra/addilab.m	(revision 2)
+++ /prextra/prextra/addilab.m	(revision 2)
@@ -0,0 +1,62 @@
+%ADDILAB add named identifier to a dataset
+%
+%   B=ADDILAB(A,REF,LABS)
+%
+% INPUT
+%    A     input dataset with ilab structure (see MAKEILAB)
+%    REF   name, cell array with names or numerical indices of
+%          identifier columns
+%    LABS  identifiers to be set. If all identifiers are numerical, 
+%          matrix may be supplied, otherwise a cell array.
+%
+% OUTPUT
+%    B     output dataset
+%
+% DESCRIPTION
+% This function adds one or more named identifiers to the dateset
+% and sets their per-example values. To be used, the named-identifier 
+% facility must be enabled in a dataset (see ENABLEILAB).
+%
+% SEE ALSO
+% ENABLEILAB, REMOVEILAB, GETILAB, SETILAB, GETILABLIST
+
+% $Id: addilab.m,v 1.2 2005/05/02 13:05:46 pavel Exp $
+
+function a=addilab(a,name,labs)
+    
+    isvalidilab(a);
+
+    u=a.user;
+    
+    % is the name already defined?
+    ind=strcmp(u.ilab.names,name);
+    if any(ind)
+	error([name ' already defined in ilab!']);
+    else
+	
+	% check the validity of labs to be added
+	if length(labs)~=size(a,1)
+	    length(labs)
+	    error('wrong size of ident labels to be included!');
+	end
+
+%	labs=reshape(labs,size(a,1),1);
+	
+	if isnumeric(labs)
+	    labs=num2cell(labs);
+	end
+	
+	if ~iscell(labs)
+	    error('ident labels to include must be either cell array or numerical vector');
+	end
+	
+	% update ilab
+	u.ilab.names{end+1}=name;
+	a.user=u;
+	ident=a.ident;
+	ident(:,end+1)=labs;
+	a.ident=ident;
+	
+    end
+    
+    return
Index: /prextra/prextra/auc.m
===================================================================
--- /prextra/prextra/auc.m	(revision 2)
+++ /prextra/prextra/auc.m	(revision 2)
@@ -0,0 +1,54 @@
+%AUC Area-Under-the-Curve error estimator
+%
+%   E = AUC(A,W)
+%   E = AUC(A*W)
+%   E = A*W*AUC
+%
+% Returns the ROC area_under_the_curve error 
+% for the test set A on the classifier W.
+%
+% Two class problems only
+
+% 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 e = auc(a,w)
+
+	if nargin == 0 | isempty(a)
+		e = mapping('auc','fixed');
+		return;
+	end
+	
+	if nargin == 1, d = a*classc; end
+	if nargin == 2, d = a*w*classc; end
+	
+	m = size(d,1);
+	c = getsize(d,3);
+	if c ~= 2
+		error('Routine implemented for two-class problems only')
+	end
+	
+	nlab = getnlab(d);
+ 	s = +d; s = sort(s(:));
+	thr = [s' 1]; 
+	
+	e1 = []; e2 = [];
+
+	% NLAB_OUT will be one where B is larger than THR.
+	I = matchlablist(d.lablist,d.featlab); % Correct possible changes class orders
+	
+	nlab_out = (repmat(+d(:,I(1)),1,m*c+1) > repmat(thr,m,1));
+
+	% ERRS will be 1 where the numeric label is unequal to NLAB_OUT
+	% (i.e., where errors occur).
+	errs = (repmat((nlab==1),1,m*c+1) ~= nlab_out);
+			
+	% Split the cases where ERRS = 1 into cases for CLAS (E1) and all 
+	% other classes (E2).
+
+	e1 = mean(errs(find(nlab==1),:),1);
+	e2 = mean(errs(find(nlab~=1),:),1);
+	
+	J = [1:m*c];
+	e = (e1(J+1)-e1(J))*(e2(J) + e2(J+1))'/2;
Index: /prextra/prextra/auclpc.m
===================================================================
--- /prextra/prextra/auclpc.m	(revision 2)
+++ /prextra/prextra/auclpc.m	(revision 2)
@@ -0,0 +1,105 @@
+%AUCLPC Area Under the Curce Linear Programming Classifier
+%
+% [W,R] = AUCLPC(A,C,rtype, par, unitnorm, usematlab)
+%
+%   A  Dataset
+%   C  Trade off parameter, default C = 10
+%   W  Linear Classifier
+%   R  Indices of support objects
+%
+% For other parameters see AUCLPM.
+% The classifier optimizes the AUC for two-class problems using AUCLPM.
+% The best threshold for the training set A is found using
+% the class priors. Multi-class problems are solved using the
+% one-against-rest scheme of MCLASSC (MODE = 'single').
+%
+% REFERENCE
+% D.M.J. Tax and C.J. Veenman, Tuning the hyperparameter of an AUC-optimized
+%    classifier, Proc. BNAIC 2005, 2005.
+%
+% SEE ALSO
+% DATASETS, MAPPINGS, AUCLPM, MCLASSC
+
+% Copyright: D.M.J. Tax, 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,r] = auclpc(a,C, rtype, par, unitnorm, usematlab)
+
+if (nargin < 6)
+	usematlab = 0;
+end
+if (nargin < 5)
+	unitnorm = 1;
+end
+if (nargin < 4)
+	par = 0.25;
+end
+if (nargin < 3)
+	rtype = 'subk';
+end
+if (nargin < 2)
+	prwarning(3,'Lambda set to ten');
+	C = 10; 
+end
+
+if nargin < 1 | isempty(a)
+    w = mapping(mfilename,{C,rtype,par,unitnorm,usematlab});
+    w = setname(w,'AUC_LP');
+    return
+end
+
+[m,k,c] = getsize(a);
+if c > 2
+    %w = mclassc(a,feval(mfilename));
+    w = mclassc(a,auclpc);
+    w = setname(w,'AUC_LP');
+    N = [];
+    for j=1:c
+        N = [N w.data{j}.data{1}.data{1}.data];
+    end
+    r = unique(N);
+    return
+end
+
+v = auclpm(a,C,rtype, par, unitnorm, usematlab);
+nlab = getnlab(a);
+p = getprior(a);
+d = +(a*v);
+[dd,J] = sort(d);
+n1 = sum(nlab==1);
+n2 = sum(nlab==2);
+%e = n1-n2;
+%J1 = cumsum(nlab(J) == 1);
+%J2 = n2-cumsum(nlab(J) == 2);
+e1 = p(1)*cumsum(nlab(J) == 1)/n1 + p(2)*(1 - cumsum(nlab(J) == 2)/n2);
+e2 = p(2)*cumsum(nlab(J) == 2)/n2 + p(1)*(1 - cumsum(nlab(J) == 1)/n1);
+
+[min1,k1] = min(e1);
+[min2,k2] = min(e2);
+%if min1 < min2          % NO!! We should assume that AUCLPM gives a proper
+%    if k1 == length(d)  % and well defined class order: 1 -> -, 2 > +
+%        w0 = d(J(k1)) + realmin;
+%    else
+%        w0 = (d(J(k1))+d(J(k1+1)))/2;
+%    end
+%    w = v.data.u - v.data.v;
+%else
+    if k2 == length(d)
+        w0 = d(J(k2)) + realmin;
+    else
+        w0 = (d(J(k2))+d(J(k2+1)))/2;
+    end
+    w = v.data.v - v.data.u;
+%end
+%J = find(w~=0);
+prec = 1e-7;   % this is really choosing the significant guys
+ss = sum(abs(w));
+J  = find(abs(w) > ss*prec);
+
+w = featsel(k,J)*affine(w(J),w0,a(:,J),getlablist(a),k,c); 
+w = setout_conv(w,1);
+w = setname(w,'AUC_LP');
+r = J;
+    
+    
Index: /prextra/prextra/auclpm.m
===================================================================
--- /prextra/prextra/auclpm.m	(revision 2)
+++ /prextra/prextra/auclpm.m	(revision 2)
@@ -0,0 +1,231 @@
+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 = 1;
+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
+
+% define the correct name:
+if unitnorm
+	cl_name = sprintf('AUC-LP %s (%s) 1norm',rtype,num2str(par));
+else
+	cl_name = sprintf('AUC-LP %s (%s)',rtype,num2str(par));
+end
+if (nargin < 1) | (isempty(x))
+	w = mapping(mfilename,{C,rtype,par,unitnorm,usematlab});
+	w = setname(w,cl_name);
+	return
+end
+
+if ~ismapping(C)   % train the mapping
+
+	% 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
+%if nargout>1
+%	returnA = A;
+%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
+constr = orgA*(z(1:k)-z(k+1:2*k));
+% the number of satisfied constraints (=AUC:)
+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);
+			message(4,'Optimum C = %f\n',C);
+		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(:);
+%if nargout>1
+%	returnA = zeta;
+%end
+	else
+		error('Only a two-class classifier is implemented');
+	end
+	% now find out how sparse the result is:
+	%nr = sum(beta>1e-6);
+	rel = (abs(u-v)>0);
+	nr = sum(rel);
+	
+	% and store the results:
+	%W.wsc = wsc;
+	W.u = u; %the ultimate weights
+	W.v = v;
+	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,sprintf('%s %s',cl_name,rtype));
+	
+else
+	% Evaluate the classifier on new data:
+	W = getdata(C);
+	n = size(x,1);
+
+	% linear classifier:
+	out = x*(W.u-W.v);
+
+	% 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);
+	%w = setdat(x,[-out out],C);
+	%w = setdat(x,sigm([-out out]),C);
+
+end
+		
+return
Index: /prextra/prextra/change_R.m
===================================================================
--- /prextra/prextra/change_R.m	(revision 2)
+++ /prextra/prextra/change_R.m	(revision 2)
@@ -0,0 +1,41 @@
+function R = change_R(R,c,beta,gammac)
+% R = change_R(R,c,beta,gammac)
+% 
+% Auxiliary function for the incremental SVDD, see there.
+
+% Copyright: D.M.J. Tax, D.M.J.Tax@prtools.org
+% Faculty EWI, Delft University of Technology
+% P.O. Box 5031, 2600 GA Delft, The Netherlands
+
+% Note that n is one larger than set S, because of the
+% added parameter b at the beginning!
+n = size(R,1);
+
+% In some unfortunate cases the gamma_c can be very very small (like 0),
+% and in these cases we don't want to blow up the whole thing. Therefore
+% we define a small eps and define that as the smallest gamma_c
+% possible.
+smalleps = 1e-12;
+
+if c>0                          % we add object c to R
+	if abs(gammac)>smalleps
+		R = [R zeros(n,1); zeros(1,n+1)] + ...
+			([beta;1]/gammac)*[beta' 1];
+		% Improve the numerical stability (overflow) for large beta and
+		% gammac, by first dividing by gammac and then multiplying again
+		% by beta.
+		% Fix suggested by Mauro Del Rio
+	else
+		warning('dd_tools:change_R:DivideByZero','We are about to divide by 0.');
+		R = [R zeros(n,1); zeros(1,n+1)] + ...
+			([beta;1]/(sign(gammac)*smalleps))*[beta' 1];
+	end
+else                       % we remove object c from R
+	c = -c; % ok, get rid of the sign
+	Irm = [1:c,c+2:n];
+	if R(c+1,c+1)>smalleps  % whatever...
+		R = R(Irm,Irm) - R(Irm,c+1)*R(c+1,Irm)/R(c+1,c+1);
+	else
+		R = R(Irm,Irm) - R(Irm,c+1)*R(c+1,Irm)/smalleps;
+	end
+end
Index: /prextra/prextra/comparch2.m
===================================================================
--- /prextra/prextra/comparch2.m	(revision 2)
+++ /prextra/prextra/comparch2.m	(revision 2)
@@ -0,0 +1,60 @@
+
+function ee = comparch2(type,L,n,iter)
+
+if nargin < 4, iter = 25; end
+if nargin < 3 | isempty(n), n = 50; end
+if nargin < 2, L = []; end
+
+W = classfs;
+if ~isempty(L), W = W(L); end
+
+rand('state',0);
+ee.names = getname(W{1});
+for j=2:length(W)
+	ee.names = char(ee.names,getname(W{j}));
+end
+
+learnsizes = n;
+	
+e = zeros(length(W),length(learnsizes),iter);
+b = genarche(type,1000,0);
+for j=1:iter
+	for m = 1:length(learnsizes)
+		a = genarche(type,learnsizes(m),j);
+		for i=1:length(W)
+			e(i,m,j) = b*(a*W{i})*testc;
+			disp([j,m,i])
+		end
+	end
+end
+
+ee.error = e;
+ee.iter = iter;
+ee.learnsizes = learnsizes;
+
+disp(type)
+disp(' ')
+for m = 1:length(learnsizes)
+	em = mean(squeeze(e(:,m,:))');
+	es = std(squeeze(e(:,m,:))');
+	[eem,J] = sort(em);
+	for i=1:length(W)
+		fprintf(1,'%7.4f %8.5f   %s \n',em(J(i)),es(J(i))/sqrt(iter),ee.names(J(i),:));
+	end
+end
+
+if isempty(L) % full set, store result
+	save([type '_' num2str(iter) '_' num2str(learnsizes(1))],'ee');
+end
+
+%figure;
+%scatterd(genarche(type,200,0));
+%axis equal
+%fsave('scatter_plot',14,2,5,1.5);
+%ee.error = mean(e,3);
+%ee.xvalues = learnsizes;
+%ee.std = std(e,[],3)/sqrt(iter);
+%plotr(ee,[],[],[],'errorbar')
+%save(['arch2_' getmapping_file(W{n})],'ee','iter');
+
+			
Index: /prextra/prextra/comparche.m
===================================================================
--- /prextra/prextra/comparche.m	(revision 2)
+++ /prextra/prextra/comparche.m	(revision 2)
@@ -0,0 +1,33 @@
+
+W = {qdc,udc,ldc,nmc,treec,lmnc,parzenc,knnc,knnc([],1),loglc,svc,naivebc,rbnc,fdsc};
+
+ee.names = getname(W{1});
+for j=2:length(W)
+	ee.names = char(ee.names,getname(W{j}));
+end
+
+learnsizes = [2,5,7,10,15,20,30,50,70,100];
+
+iter = 25;
+for n = 1:8
+	
+e = zeros(length(W),length(learnsizes),iter);
+b = genarche(W{n},1000,0);
+for j=1:iter
+	for m = 1:length(learnsizes)
+		a = genarche(W{n},learnsizes(m),j);
+		for i=1:length(W)
+			e(i,m,j) = b*(a*W{i})*testc;
+			disp([n,j,m,i])
+		end
+	end
+end
+
+ee.error = mean(e,3);
+ee.xvalues = learnsizes;
+ee.std = std(e,[],3)/sqrt(iter);
+plotr(ee,[],[],[],'errorbar')
+save(['arche_' getmapping_file(W{n})],'ee','iter');
+
+end
+			
Index: /prextra/prextra/createA.m
===================================================================
--- /prextra/prextra/createA.m	(revision 2)
+++ /prextra/prextra/createA.m	(revision 2)
@@ -0,0 +1,194 @@
+function [A,Nxi,A2] = createA(X,y,rtype,par,seed)
+% [A,Nxi,A2] = CREATEA(X,Y,RTYPE,PAR,SEED)
+%
+% Create the data matrix containing all pairwise difference vectors in
+% data matrix X (with their corresponding labels Y, -1/+1).
+% Because the size of this data matrix can become huge (ALL pairwise
+% difference vectors is a lot!), you can subsample it by choosing an
+% appropriate RTYPE.
+%
+%  RTYPE 'full'   use all constraints
+%        'subs'   randomly subsample PAR constraints 
+%        'subk'   randomly subsample a fraction PAR of the constraints
+%        'knn'    use the PAR nearest neighbors in the other class
+%        'xval'   subsample and use remaining constraints to optimize C
+%        'xvalk'  subsample a fraction k*n and use remaining constraints
+%                 to optimize C
+%        'kmeans' use k-means clustering with k=PAR
+%        'randk'  subsample objects to get PAR*(Npos+Nneg) constraints
+%
+% The SEED is optional, it is the seed for the random sampling.
+%
+if nargin<5
+	seed = [];
+end
+% If a seed is defined, set it:
+if ~isempty(seed)
+	rand('state',seed);
+end
+
+A2 = [];
+%---create A for optauc
+
+k = size(X,2);
+
+% compute how many xi-s we expect:
+Ineg = find(y==-1);
+Ipos = find(y==+1);
+Nneg = length(Ineg);
+Npos = length(Ipos);
+
+% depending on the reduction type
+switch rtype
+case 'full'  % take all the possibilities:
+	Nxi = Nneg*Npos;
+	A = zeros(Nxi,k);
+	% run over all possibilities:
+	dummyk=0;
+	for i=1:Nneg
+		for j=1:Npos
+			dummyk = dummyk+1;
+			A(dummyk,:) = X(Ineg(i),:)-X(Ipos(j),:);
+		end
+	end
+case 'subk'  % subsample the possibilities, but now not a fixed number,
+	%but k times the number of training objects:
+	Nxi = ceil(par*size(X,1));
+	A = zeros(Nxi,k);
+	Ip = floor(Npos*rand(Nxi,1))+1; Ip = Ip(1:Nxi);
+	In = floor(Nneg*rand(Nxi,1))+1; In = In(1:Nxi);
+	for i=1:Nxi
+		diffx = X(Ineg(In(i)),:) - X(Ipos(Ip(i)),:);
+		A(i,:) = diffx;
+	end
+case 'subs'  % subsample the possibilities:
+	Nxi = par;
+	A = zeros(Nxi,k);
+	Ip = floor(Npos*rand(Nxi,1))+1; Ip = Ip(1:Nxi);
+	In = floor(Nneg*rand(Nxi,1))+1; In = In(1:Nxi);
+	for i=1:Nxi
+		diffx = X(Ineg(In(i)),:) - X(Ipos(Ip(i)),:);
+		A(i,:) = diffx;
+	end
+case 'knn'  % only use the k nearest neighbors
+	Nxi = ceil((Nneg+Npos)*par);
+	A = zeros(Nxi,k);
+	% first process all the neg. examples:
+	D = sqeucldistm(X(Ineg,:),X(Ipos,:));
+	[dummy,I] = sort(D,2);
+	dummyk = 0;
+	for i=1:Nneg
+		for j=1:par
+			thispos = I(i,j);
+			diffx = X(Ineg(i),:)-X(Ipos(thispos),:);
+			dummyk = dummyk+1;
+			A(dummyk,:) = diffx;
+		end
+	end
+	% then to all the pos. examples:
+	D = D';  % (no need to recompute D)
+	[dummy,I] = sort(D,2);
+	for i=1:Npos
+		for j=1:par
+			thispos = I(i,j);
+			diffx = -X(Ipos(i),:)+X(Ineg(thispos),:);
+			dummyk = dummyk+1;
+			A(dummyk,:) = diffx;
+		end
+	end
+case 'randk'  % randomly chosen objs such that you have k(Npos+Nneg)
+	           % constraints
+	q = sqrt(par*(Npos+Nneg)/(Npos*Nneg));
+	qpos = ceil(q*Npos); qneg = ceil(q*Nneg);
+	Nxi = qpos*qneg;
+	A = zeros(Nxi,k);
+
+	% first select the neg. examples:
+	I = randperm(Nneg); In = Ineg(I(1:qneg));
+	% then select the pos. examples:
+	I = randperm(Npos); Ip = Ipos(I(1:qpos));
+	% run over all possibilities:
+	dummyk=0;
+	for i=1:qneg
+		for j=1:qpos
+			dummyk = dummyk+1;
+			A(dummyk,:) = X(In(i),:)-X(Ip(j),:);
+		end
+	end
+case 'xval'  % take all the possibilities and use part for testing:
+	Nxi = Nneg*Npos;
+	A = zeros(Nxi,k);
+	% run over all possibilities:
+	dummyk=0;
+	for i=1:Nneg
+		for j=1:Npos
+			diffx = X(Ineg(i),:)-X(Ipos(j),:);
+			dummyk = dummyk+1;
+			A(dummyk,:) = diffx;
+		end
+	end
+	% get part of data for constraints, the rest for evalation:
+	I = randperm(Nxi);
+   if par>=size(A,1)
+		warning(sprintf('More constraints requested than available (%d and %d)',par,size(A,1)));
+		disp('Now using half for training and testing');
+		par = ceil(size(A,1)/2);
+	end
+	% if data is really really huge, then subsample more...
+	Mega=100000;
+	if length(I)-par>Mega
+		A2 = A(I((par+1):(par+Mega)),:);
+	else
+		A2 = A(I((par+1):end),:);
+	end
+	A = A(I(1:par),:);
+	Nxi = par;
+case 'xvalk'  % take all the possibilities and use part for testing:
+	par = par*size(X,1);
+	Nxi = Nneg*Npos;
+	A = zeros(Nxi,k);
+	% run over all possibilities:
+	dummyk=0;
+	for i=1:Nneg
+		for j=1:Npos
+			diffx = X(Ineg(i),:)-X(Ipos(j),:);
+			dummyk = dummyk+1;
+			A(dummyk,:) = diffx;
+		end
+	end
+	% get part of data for constraints, the rest for evalation:
+	I = randperm(Nxi);
+   if par>=size(A,1)
+		warning(sprintf('More constraints requested than available (%d and %d)',par,size(A,1)));
+		disp('Now using half for training and testing');
+		par = ceil(size(A,1)/2);
+	end
+	% if data is really really huge, then subsample more...
+	Mega=100000;
+	if length(I)-par>Mega
+		A2 = A(I((par+1):(par+Mega)),:);
+	else
+		A2 = A(I((par+1):end),:);
+	end
+	A = A(I(1:par),:);
+	Nxi = par;
+case 'kmeans'
+	wp = kmeans_dd(X(Ipos,:),0.1,par);
+	wn = kmeans_dd(X(Ineg,:),0.1,par);
+	Xp = wp.data.w;
+	Xn = wn.data.w;
+	Nxi = par*par;
+	A = zeros(Nxi,k);
+	dummyk=0;
+	for i=1:par
+		for j=1:par
+			diffx = Xn(i,:)-Xp(j,:);
+			dummyk = dummyk + 1;
+			A(dummyk,:) = diffx;
+		end
+	end
+otherwise
+	error(sprintf('Type %s is not defined',rtype));
+end
+
+return
Index: /prextra/prextra/delfigs.m
===================================================================
--- /prextra/prextra/delfigs.m	(revision 2)
+++ /prextra/prextra/delfigs.m	(revision 2)
@@ -0,0 +1,3 @@
+%DELFIGS Delete all figures
+
+delete(get(0,'children'))
Index: /prextra/prextra/ecoc.m
===================================================================
--- /prextra/prextra/ecoc.m	(revision 2)
+++ /prextra/prextra/ecoc.m	(revision 2)
@@ -0,0 +1,139 @@
+%ECOC Error-correcting output code
+% 
+%      W = ECOC(A,BASECL,CM,RULE)
+%
+% INPUT
+%    A       Dataset
+%    BASECL  Untrained base classifier
+%    CM      Coding matrix
+%    RULE    Combination rule
+%
+% OUTPUT
+%    W       ECOC classifier
+%
+% DESCRIPTION
+% Computation of the classifier using Error-correcting output codes to
+% create a multi-class classifier from a base classifier/dichotomizer.
+% The classes in dataset A are relabeled according to the coding matrix
+% CM (containing +1, -1 or 0). The default coding matrix is a 1-vs-rest
+% coding. For a three-class problem this becomes:
+%     CM = [ 1 -1 -1;
+%           -1  1 -1;
+%           -1 -1  1];
+%
+% For the evaluation of objects, the outputs of the dichotomizers are
+% combined using the combination rule RULE. Currently the following
+% rules are implemented:
+%     'none'     take the maximum output (this only works for 1-vs-all
+%                coding matrices)
+%     'hamming'  standard Hamming distance between the discretized
+%                classifier outcomes and the coding matrix. When a zero
+%                entry in the coding matrix appears, this entry is not
+%                used in the distance computation
+%
+% NOTE: the order of the classes as they are used in the coding matrix
+% is determined by the lablist of A.
+%
+% SEE ALSO
+%    MCLASSC, MAXC
+
+% Copyright: D.M.J. Tax, D.M.J.Tax@prtools.org
+% Faculty EWI, Delft University of Technology
+% P.O. Box 5031, 2600 GA Delft, The Netherlands
+  
+function W = ecoc(a,basecl,CM,rule)
+
+if nargin < 4 | isempty(rule), rule = 'hamming'; end
+if nargin < 3 | isempty(CM), CM = []; end
+if nargin < 2 | isempty(basecl), basecl = fisherc; end
+if nargin < 1 | isempty(a) 
+	W = mapping(mfilename,{w,L,rule});
+	W = setname(W,'ECOC');
+	return
+end
+
+if ~ismapping(basecl) | ~istrained(basecl)   %training
+
+	[n,k,c] = getsize(a);
+	nlab = getnlab(a);
+	% define the coding matrix if not given:
+	if isempty(CM)
+		CM = repmat(-1,c,c) + 2*eye(c);
+	else
+		% check if the coding matrix is OK
+		if size(CM,1)~=c
+			error('Coding matrix does not have correct number of rows.');
+		end
+		if (max(CM(:))>1)|(min(CM(:))<-1)
+			error('Entries in the coding matrix should be -1 or +1.');
+		end
+	end
+	l = size(CM,2);
+   w = cell(1,l);
+
+	% Train it:
+	for i=1:l
+		% create a new dataset
+		lab = nlab;
+		for j=1:c
+			lab(lab==j) = CM(j,i);
+		end
+		% should I remove the objects with 0 label?
+		b = setlabels(a,lab);
+		% train the base dichotomizer:
+		w{i} = b*basecl;
+	end
+
+	%and save all useful data:
+	W.CM = CM;
+	W.w = w;
+	W.rule = rule;
+	W = mapping(mfilename,'trained',W,getlablist(a),k,c);
+	W = setname(W,'ECOC');
+
+else                               %testing
+
+	% Extract the data:
+	W = getdata(basecl);
+	m = size(a,1);
+	[c,l] = size(W.CM);
+
+	% Apply all the base classifiers:
+	z = zeros(m,l);
+	for i=1:l
+		tmp = a*W.w{i};
+		id = findfeatlab(tmp,1);
+		z(:,i) = +tmp(:,id);
+	end
+	z = 2*z-1;  %DXD classifier output should be between -1 and +1!!!
+%[W.CM; z]
+	% and apply the combination rule:
+	out = zeros(m,c);
+	switch W.rule
+	case 'none'
+		if size(z,2)~=c
+			error('You have to combine when you don''t do one-vs-rest.');
+		end
+		out = z;
+	case 'hamming'
+		for i=1:m
+			zz = sign(z(i,:)); % discrete output
+			% hamming distance of object i to all code words:
+			% remove entries for which the coding matrix has a 0:
+			CM0 = (W.CM~=0);
+			d = sum(CM0.*abs(repmat(zz,c,1)-W.CM),2);
+%keyboard
+			[md,Icl] = min(d);
+			% class Icl won!
+			out(i,Icl) = 1;
+		end
+	otherwise
+		error('This combination rule is not implemented yet.');
+	end
+
+	% Store the outputs
+	W = setdat(a,out,basecl);
+end
+return
+
+
Index: /prextra/prextra/emc.m
===================================================================
--- /prextra/prextra/emc.m	(revision 2)
+++ /prextra/prextra/emc.m	(revision 2)
@@ -0,0 +1,88 @@
+%EMC EM Classifier using semi-supervised data
+%
+%   W = EMC(A,B,CLASSF,LABTYPE,FID)
+%   W = A*EMC([],B,CLASSF,LABTYPE,FID)
+%
+% INPUT
+%   A       Labeled dataset used for training
+%   B       Additional unlabeled dataset
+%   CLASSF  Untrained classifier (default QDC)
+%   LABTYPE Label type to be used (crisp (default) or soft)
+%   FID     File ID to write progress to (default [], see PRPROGRESS)
+%
+% OUTPUT
+%   W      Trained classifier
+%
+% DESCRIPTION
+% Using the EM algorithm the classifier CLASSF is used iteratively
+% on the joint dataset [A;B]. In each step the labels of A are reset
+% to their initial values. Initial labels in B are neglected.
+% Labels of LABTYPE 'soft' are not supported by all classifiers. 
+%
+% SEE ALSO
+% DATASETS, MAPPINGS, EMCLUST, PRPROGRESS
+
+% 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 = emc(a,b,classf,labtype,fid)
+	if nargin < 5, fid = []; end
+	if nargin < 4 | isempty(labtype), labtype = 'crisp'; end
+	if nargin < 3 | isempty(classf), classf = qdc; end	
+	if nargin < 2, b = []; end
+	if nargin < 1 | isempty(a)
+		w = mapping(mfilename,'untrained',{b,classf,labtype,fid});
+		w = setname(w,'EM CLassifier');
+		return
+	end
+
+	islabtype(a,'crisp','soft');
+	isvaldset(a,1,2); % at least 2 object per class, 2 classes
+	if isempty(b)
+		w = a*classf;
+		return
+	end
+	if size(a,2) ~= size(b,2)
+		error('Datasets should have same number of features')
+	end
+	
+	c = getsize(a,3);
+	epsilon = 1e-6;
+	change = 1;
+	nlab = getnlab(a);
+	lablist = getlablist(a);
+	a = setlabels(a,nlab); 
+	a = setlabtype(a,labtype);
+	switch labtype
+		case 'crisp'
+			lab = zeros(size(b,1),1);
+		case 'soft'
+			lab = zeros(size(b,1),c);
+	end
+	b = dataset(+b);
+	w = a*classf;
+	
+	prprogress(fid,'\nem_classifier optimization\n')		
+	while change > epsilon
+		d = b*w;
+		switch labtype
+			case 'crisp'
+				labb = d*labeld;
+				change = mean(lab ~= labb); 
+			case 'soft'
+				labb = d*classc;
+				change = mean(mean((+(labb-lab)).^2));
+			otherwise
+				error('Wrong LABTYPE given')
+		end
+		lab = labb;
+		b = setlabtype(b,labtype,lab);
+		c = [a; b];
+		w = c*classf;
+		prprogress(fid,'  change = %d\n', change)
+	end
+	
+	J = getlabels(w);
+	w = setlabels(w,lablist(J,:));
+	
Index: /prextra/prextra/emparzenc.m
===================================================================
--- /prextra/prextra/emparzenc.m	(revision 2)
+++ /prextra/prextra/emparzenc.m	(revision 2)
@@ -0,0 +1,94 @@
+%EMPARZENC EM-algorithm for semi-supervised learning by parzenc
+%
+%   [W,V] = EMPARZENC(A,B,N,FID)
+%   W = A*EMPARZENC([],B,N,FID)
+%
+% INPUT
+%   A       Labeled dataset used for training
+%   B       Additional unlabeled dataset
+%   N       Number of smoothing parameter steps (default 1)
+%   FID     File ID to write progress to (default [], see PRPROGRESS)
+%
+% OUTPUT
+%   W      Trained classifier, based on A and B
+%   V      Trained classifier based on A only
+%
+% DESCRIPTION
+% Using the EM algorithm the PARZENC classifier is used iteratively
+% on the joint dataset [A;B]. In EM each step the labels of A are reset
+% to their initial values. Initial labels in B are neglected. They
+% are iteratively updated as soft labels obtained by classifying B
+% by the actual W. The EM algorithm is run for a fixed smoothing
+% parameter of PARZENC. This is repeated for smaller smoothing 
+% parameters in N steps, using harmonic interpolation between HL and HU,
+% in which HL is the smoothing parameter estimate obtained from PARZENML
+% applied to A and HU the estimate obtained from PARZENML applied to B.
+% For N = 1, the average of HL and HU is used.
+%
+% SEE ALSO
+% DATASETS, MAPPINGS, EMCLUST, EMC, PARZENC, PARZENML, PRPROGRESS
+
+% 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,v] = emc(a,b,n,fid)
+	if nargin < 4, fid = []; end
+	if nargin < 3, n = 1; end
+	if nargin < 2, b = []; end
+	if nargin < 1 | isempty(a)
+		w = mapping(mfilename,'untrained',{b,classf,labtype,fid});
+		w = setname(w,'EMParzen CLassifier');
+		return
+	end
+
+	if size(a,2) ~= size(b,2)
+		error('Datasets should have same number of features')
+	end
+	
+	c = getsize(a,3);
+	epsilon = 1e-6;
+	nlab = getnlab(a);
+	lablist = getlablist(a);
+	a = setlabels(a,nlab); 
+	a = setlabtype(a,'soft');
+	%ws = scalem([+a; +b],'variance');
+	ws = unitm;
+	a = a*ws;
+	b = b*ws;
+	lab = zeros(size(b,1),c);
+	hl = parzenml(a);
+	b = dataset(+b);
+	hu = parzenml(b);
+	hl = max(hl,hu) * 1.05;
+	hu = min(hl,hu) * 0.95;
+	if n == 1
+		h = (hl+hu)/2;
+	else
+		dh = (log(hl) - log(hu))/(n-1);
+		h = exp([log(hl):-dh:log(hu)]);
+	end
+	c = a;
+	
+	first = 1;
+	for j=1:length(h)
+		hh = h(j);
+		prprogress(fid,['\nem_classifier optimization, h = ' num2str(hh) '\n'])	
+		change = 1;
+		while change > epsilon
+			w = parzenc(c,hh);
+			if first, v = w; first = 0; end
+			d = b*w;
+			labb = d*classc;
+			change = mean(mean((+(labb-lab)).^2));
+			lab = labb;
+			b = setlabtype(b,'soft',lab);
+			c = [a; b];
+			prprogress(fid,'  change = %d\n', change)
+		end
+	end
+	
+	J = getlabels(w);
+	w = ws*setlabels(w,lablist(J,:));
+	v = ws*setlabels(v,lablist(J,:));
+	
Index: /prextra/prextra/enableilab.m
===================================================================
--- /prextra/prextra/enableilab.m	(revision 2)
+++ /prextra/prextra/enableilab.m	(revision 2)
@@ -0,0 +1,37 @@
+%ENABLEILAB - enable named identifiers in a dataset (ilab)
+%
+%    B=ENABLEILAB(A)
+%
+% INPUT
+%    A   input dataset
+%
+% OUTPUT
+%    B   output dataset with initialized ilab structure
+%
+% DESCRIPTION
+% Enables the use of multiple named identifiers in a dataset. For each
+% example in a dataset, multiple identifiers may be defined. Each
+% identifier has a name. The names (identifier labels i.e. ilab) are stored
+% in a ilablist in a user field of a dataset. This function removes all the
+% identifiers stored in a dataset, add one column, called ident, and
+% populates it with unique indices.
+%
+% Technically, ENABLEILAB adds a field ilab to the user structure.  The
+% field names of the ilab stores the column names of the dataset ident
+% cell-array. Other fields of user structure than ilab are not altered.
+%
+% SEE ALSO
+% ADDILAB, REMOVEILAB, GETILAB, SETILAB, GETILABLIST
+
+function a=enableilab(a)
+
+    % reseting the identifiers
+    a.ident=num2cell((1:size(a,1))');
+    prwarning(3,'identifiers were reset in the dataset');
+    
+    ilab.names={'ident'};
+    u.ilab=ilab;
+    
+    a.user=u;
+    
+    return
Index: /prextra/prextra/enablemeta.m
===================================================================
--- /prextra/prextra/enablemeta.m	(revision 2)
+++ /prextra/prextra/enablemeta.m	(revision 2)
@@ -0,0 +1,32 @@
+%ENABLEMETA  enable the storage of key-value pairs in a user field
+%
+%    B=ENABLEMETA(A)
+%
+% INPUT
+%    A   input dataset
+%
+% OUTPUT
+%    B   output dataset with initialized meta structure
+%
+% DESCRIPTION
+% Enables the use of named key-value pairs in a dataset. These are
+% stored inside the user field of a dataset.
+%
+% Technically, ENABLEMETA adds a field meta to the user structure.  It
+% contains the cell-arrays 'keys' and 'values'. Other fields of the user
+% structure than meta are not altered.
+%
+% SEE ALSO
+% ADDMETA, REMOVEMETA, GETMETA, SETMETA
+
+function a=enablemeta(a)
+    
+    u=a.user;
+
+    meta.keys={};
+    meta.values={};
+    u.meta=meta;
+    
+    a.user=u;
+    
+    return    
Index: /prextra/prextra/expandd.m
===================================================================
--- /prextra/prextra/expandd.m	(revision 2)
+++ /prextra/prextra/expandd.m	(revision 2)
@@ -0,0 +1,40 @@
+%EXPANDD Expand integer vector to a matrix occurrence table
+% 
+%   A = EXPANDD(X,M)
+% 
+% INPUT
+%   X   Vector containing positive integers
+%   M   Dimensionality of the returned matrix A
+%
+% OUTPUT
+%   A   Matrix 
+%
+% DESCRIPTION
+% The vector X (e.g. numeric labels obtained from RENUMLAB) is expanded to 
+% a matrix A of the size [LENGTH(X) x M] such that A(i,j) = 1 if X(i) == j 
+% and A(i,j) = 0, otherwise.
+% As a result SUM(A) is a frequency table of j in X. MAX(A) is a 0/1 vector 
+% indicating the occurrence of some j in X. FIND(MAX(A)) gives all j that 
+% occur in X.
+% 
+% SEE ALSO 
+% RENUMLAB
+
+% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
+% Faculty of Applied Physics, Delft University of Technology
+% P.O. Box 5046, 2600 GA Delft, The Netherlands
+
+% $Id: expandd.m,v 1.1 2005/10/26 08:27:30 duin Exp $
+
+function a = expandd(x,m)
+
+	prtrace(mfilename);
+
+	x = x(:);
+	n = length(x);
+	if (nargin == 1) 
+		m = max(x); 
+	end
+	a = (x(:,ones(1,m)) == ones(n,1) * linspace(1,m,m));
+
+return;
Index: /prextra/prextra/featsel2.m
===================================================================
--- /prextra/prextra/featsel2.m	(revision 2)
+++ /prextra/prextra/featsel2.m	(revision 2)
@@ -0,0 +1,85 @@
+%FEATSEL2 Pairwise feature selection for classification
+% 
+% [W,R] = FEATSEL2(A,CRIT,K,T,FID)
+%
+% INPUT	
+%   A    Training dataset
+%   CRIT Name of the criterion or untrained mapping 
+%        (default: 'NN', i.e. the 1-Nearest Neighbor error)
+%   K    Number of features to select (default: K = 0, return optimal set)
+%   T    Tuning dataset (optional)
+%   FID  File ID to write progress to (default [], see PRPROGRESS)
+%
+% OUTPUT
+%   W    Output feature selection mapping
+%   R    Matrix with step-by-step results
+%
+% DESCRIPTION
+%
+% Best pairs of features are evaluated inidividually, i.e. independently
+% from other pairs. 
+% 
+% SEE ALSO 
+% MAPPINGS, DATASETS, FEATEVAL, FEATSELF, FEATSELLR, FEATSEL,
+% FEATSELO, FEATSELB, FEATSELI, FEATSELP, FEATSELM, PRPROGRESS
+
+% Copyright: A. Harol, 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,r] = featself(a,crit,ksel,t,fid)
+		 
+  prtrace(mfilename);
+	
+  if (nargin < 2) | isempty(crit)
+    prwarning(2,'no criterion specified, assuming NN');
+    crit = 'NN';
+  end
+  if (nargin < 3) | isempty(ksel)
+    ksel = 0;
+  end
+  if (nargin < 4)
+    prwarning(3,'no tuning set supplied (risk of overfit)');
+    t = [];
+  end
+	if (nargin < 5)
+		fid = [];
+	end
+	
+	if nargin == 0 | isempty(a)
+		% Create an empty mapping:
+		w = mapping(mfilename,{crit,ksel,t});
+	else	
+		[m,k] = size(a);
+		for j1=1:k
+			for j2=j1+1:k
+				J(j1,j2) = feateval(a(:,[j1,j2]),crit,t);
+			end
+		end
+		w = pairs_sort(J, ksel);
+		w = featsel(k,w);
+		prprogress(fid,'featself  finished\n')
+	end
+	w = setname(w,'Pairwise FeatSel');
+
+return
+
+function mn = pairs_sort(J, MAX_F)      
+% Artsiom Harol
+	if MAX_F == 0, MAX_F = size(J,1); end
+  ind2=find(J(:));
+  a=J(ind2);
+  [b,ind_sort]=sort(-a);
+  %ind_sort=(flipud(ind_sort'))';
+  ind2=ind2(ind_sort);
+  [m,n]=ind2sub(size(J),ind2);
+  
+  mn=[m(:)';n(:)'];
+ % mn=[m(1:MAX_F)';n(1:MAX_F)'];
+  [umn] = unique(mn);
+  for k=1:length(umn)
+    ind = find(umn(k) == mn);
+    mn(ind(2:end)) = [];
+  end  
+  mn = mn(1:MAX_F)';
+return;  
Index: /prextra/prextra/findilab.m
===================================================================
--- /prextra/prextra/findilab.m	(revision 2)
+++ /prextra/prextra/findilab.m	(revision 2)
@@ -0,0 +1,41 @@
+%FINDILAB find indices of examples having a given identifier
+%
+%   LABS=FINDILAB(A,NAME,VALUES)
+%
+% INPUT
+%    A        input dataset with ilab structure (see ENABLEILAB)
+%    NAME     name of identifier to search
+%    VALUES   value, a vector, or a cell array with identifier values 
+%
+% OUTPUT
+%    IND      indices of examples containing identifier VALUES
+%
+% DESCRIPTION
+% FINDILAB returns indices of examples that contain any of specified
+% values.
+%
+% SEE ALSO
+% ENABLEILAB, REMOVEILAB, findilab, SETILAB, findilabLIST
+
+% $Id: findilab.m,v 1.1 2005/05/03 15:00:21 pavel Exp $
+
+function ind=findilab(a,name,values)
+    
+    isvalidilab(a,name);
+    
+    t=getilab(a,name);
+    
+    ind=[];
+    
+    if isnumeric(values)
+	for i=1:length(values)
+	    ind=[ind; find(t==values(i))];
+	end
+    end
+    if iscell(values)
+	for i=1:length(values)
+	    ind=[ind; find(t==values{i})];
+	end
+    end
+    
+    return
Index: /prextra/prextra/fsellpc.m
===================================================================
--- /prextra/prextra/fsellpc.m	(revision 2)
+++ /prextra/prextra/fsellpc.m	(revision 2)
@@ -0,0 +1,286 @@
+%FSELLPC LP (linear programming) classifier
+%
+%   [W1,W2,W3] = fsellpc(A,BIAS,TYPE,PARAM)
+%
+% INPUT
+%   A     Dataset
+%   BIAS  YES or NO (optional; default: 1 (YES))
+%   TYPE  Type of a classifier
+%         'SIMPLE'   - the most simple formulation; all features used; PARAM = [];
+%         'STANDARD' - minimization of the training misclassification errors;
+%                      all features used; PARAM = [];
+%         'C-SPARSE' - feature selection (sparse solution); a formulation similar
+%                      to the LP_1 SVM; PARAM is a tradeoff parameter.
+%                      (optional; DEFAULT: 1).
+%         'MU-SPARSE' - feature selection (sparse solution); a formulation similar
+%                       to the LP_1 SVM, based on the paper of Graepel, Herbrich, Smola etc
+%                       'Classification on proximity data with LP-machines'.
+%                       PARAM is a tradeoff parameter, usually PARAM = 0.05 or 0.1.
+%                       It is an upper bound on the misclassfied training objects.
+%                       So, for well separable problems, PARAM = 0.01 or PARAM = 0.02.
+%                       (optional; DEFAULT: the LOO 1-NN error * 1.3).
+%  PARAM Parameter connected to the TYPE, as above
+%
+% OUTPUT
+%   W1  LP Classifier in the complete space
+%   W2  LP Classifier in a reduced space
+%   W3  Feature selection mapping; the indices of chosen features are in +W3.
+%
+% DEFAULTS
+%   BIAS  = 1
+%   TYPE  = 'STANDARD'
+%   PARAM = []
+%
+% DESCRIPTION
+% Classification problem on a N x M data A with LP-machines. Assume a two-class
+% problem. Let DLPC select J features. Then:
+% W1 is an M x 2 classifier in the original space, W2 is an J x 2 classifier
+% in the feature space defined by the J chosen features and W3 is an M x R feature
+% selection such that W1 = W3 * W2. Note that the indices of the selected features
+% can be retrieved by +W3.
+%
+% A linear classifier is built on A:
+%
+%     f(A(x,F)) = A(x,F) * w + w0,
+%
+% where F are the features and w are the weights. If BIAS is 1, then w0 is
+% also sought, otherwise it equals 0. This means that the hyperplane is
+% forced to go through the origin.
+%
+% For C-class problems, C classifiers are trained, one against all others.
+% In such a case, only W1 is returned and W3 in now NOT a feature selection,
+% but directly the indices of the selected features.
+%
+% DEFAULT:
+%   BIAS  = 1
+%   TYPE  = 'STANDARD'
+%   PARAM = 1
+%
+
+% Elzbieta Pekalska, Robert P.W. Duin, e.pekalska@tudelft.nl
+% Faculty of Electrical Engineering, Mathematics and Computer Science,
+% Delft University of Technology, The Netherlands.
+
+
+
+
+function [W1,W2,W3] = fsellpc(a,is_w0,type,par,usematlab,prec)
+
+if nargin < 6, prec = 1e-7; end
+if nargin < 5, usematlab = 0; end
+if nargin < 3 | isempty(type), type = 'standard'; end
+if nargin < 4 | isempty(par),
+  switch upper(type)
+    case {'MU-SPARSE'}
+      par = max(0.01,1.3*testk(a,1));     % upperbound error: 1.3 * loo 1-nn error
+    case {'C-SPARSE'}
+      par = 1;
+    case {'SIMPLE','STANDARD'},
+      par = [];
+    otherwise
+      disp(type)
+      error('Wrong type.')
+  end
+end
+if nargin < 2 | isempty(is_w0), is_w0 = 1; end
+if nargin < 1 | isempty(a)
+  W1 = mapping(mfilename,{is_w0,type,par,usematlab});
+  W1 = setname(W1,'FSELLPC');
+  W2 = [];
+  W3 = [];
+  return
+end
+
+
+
+if ~isdataset(a),
+  error('The first parameter should be a dataset.')
+end
+if ~isnumeric(is_w0) | (is_w0 ~= 0 & is_w0 ~= 1),
+  error('The second parameter should be 0 or 1.');
+end
+
+
+lab      = getnlab(a);
+lablist  = getlablist(a);
+[m,k,C]  = getsize(a);
+
+
+z = (is_w0 > 0);  % is the bias used or not?
+
+% This is the status of the optimization procedure.
+% For GLPK, this is the exit code; see GLPKMEX for details.
+% For Matlab LINPROG, if negative then no solution is found.
+
+status = 1;
+
+
+% If more than 2 classes, train the classifier one-against-all.
+if C > 2,
+
+% W1 = mclassc(a,mapping(mfilename,{is_w0,type,par,usematlab}));
+
+  W1 = [];
+  W2 = [];
+  W3 = [];
+  N  = [];
+  for i=1:C
+    mlab = 2 - (lab == i);
+    aa   = dataset(+a,mlab);
+    [v1,v2,v3]= fsellpc(aa,is_w0,type,par,usematlab);
+    j = +v3;
+    if isempty(v1),
+      W1 = [];
+      W2 = [];
+      W3 = [];
+      prwarning(1,'No solution found.');
+      return;
+    end
+    W1 = [W1,setlabels(v1(:,1),lablist(i,:))];
+    W2 = [W2;setlabels(v2(:,1),lablist(i,:))];
+    W3(j) = ones(length(j),1);
+    N = [N j];
+  end
+  [N1,N2,N3] = unique(N);
+  W3 = featsel(k,N1);
+  W2 = featsel(length(N1),N3)*W2;
+  return
+
+else
+
+  Y1 = 3 - 2 * lab;   % labels     +/-1
+  Y  = ones(k,1);
+
+  alpha(1:k+1,1) = 0;
+
+  aa = +a;
+  switch upper(type)
+    case {'SIMPLE'}
+      f = zeros(k+z,1);
+      b = -ones(m,1);
+      if is_w0,
+        A = -[(Y1*Y').* aa  Y1];
+      else
+        A = -[(Y1*Y').* aa];
+      end
+      [al,fval,status] = linprog(f,A,b);
+      alpha(1:k+z) = al;
+
+
+
+    case {'STANDARD'}
+      L  = ones(k,1)/k;
+
+      f  = [zeros(k+z,1); L];
+      lb = [-Inf .*ones(k+z,1); zeros(k,1)];
+      ub = Inf .* ones(2*k+z,1);
+      b  = -ones(m,1);
+      if is_w0,
+        A  = -[(Y1*Y').* aa  Y1  eye(m,k)];
+      else
+        A  = -[(Y1*Y').* aa  eye(m,k)];
+      end
+      [al,fval,ststus] = linprog(f,A,b,[],[],lb,ub);
+      alpha(1:k+z) = al(1:k+z);
+
+
+
+    case {'C-SPARSE'}
+      L  = ones(k,1);
+      ub = Inf .* ones(3*k+z,1);
+      lb = [zeros(2*k,1); -Inf.*ones(z,1); zeros(k,1)];
+      b  = -ones(m,1);
+      ay = (Y1*Y').* aa;
+      if is_w0,
+        f  = [ones(2*k,1); 0; par*L];
+        A  = -[ay  -ay  Y1 eye(m,k)];
+      else
+        f  = [ones(2*k,1); par*L];
+        A  = -[ay  -ay  eye(m,k)];
+      end
+      if (exist('glpkmex')>0) & (usematlab==0)
+        smin     = 1;  % solve minimum
+        ctype    = char(ones(m,1)*abs('U'));      % sign of inequalities
+        vartype  = char(ones(3*k+z,1)*abs('C'));  % continous variables
+%       lpsolver = 1;                             % Revised Simlex Method
+        lpsolver = 2;                             % Interior Point Method
+        params.msglev = 0; % no outputs
+        [sss,hostname] = unix('hostname');
+        hostname = hostname(1:end-1);
+        if strcmp(hostname,'saturnus') | strcmp(hostname,'polaris') | strcmp(hostname,'neptunus')
+          [al,fval,status] = glpkmex_redhat(smin,f,A,b,ctype,lb,ub,vartype,params,lpsolver);
+        else
+          [al,fval,status] = glpkmex(smin,f,A,b,ctype,lb,ub,vartype,params,lpsolver);
+        end
+      else
+        [al,fval,status] = linprog (f,A,b,[],[],lb,ub);
+      end
+      alpha(1:k) = al(1:k) - al(k+1:2*k);
+      if is_w0,
+        alpha(k+1) = al(2*k+1);
+      end
+
+
+    case {'MU-SPARSE'}
+      L   = ones(k,1)/k;
+      f   = [zeros(2*k+z,1); L; -par];
+      ub  = Inf .* ones(3*k+1+z,1);
+      lb  = [zeros(2*k,1); -Inf.*ones(z,1); zeros(k+1,1)];
+      Aeq = [ones(2*k,1); zeros(k+1+z,1)]';
+      beq = 1;
+      b   = zeros(m,1);
+      ay  = (Y1*Y').* aa;
+
+      if is_w0,
+        A  = -[ay  -ay  Y1  eye(m,k) -ones(m,1)];
+      else
+        A  = -[ay -ay eye(m,k) -ones(m,1)];
+      end
+
+      if (exist('glpkmex')>0) & (usematlab==0)
+        smin     = 1;  % solve minimum
+        ctype    = char([ones(m,1)*abs('U'); 'S']);  % sign of inequalities
+        vartype  = char(ones(3*k+1+z,1)*abs('C'));   % continous variables
+%       lpsolver = 1;                                % Revised Simlex Method
+        lpsolver = 2;                                % Interior Point Method
+        params.msglev = 0; % no outputs
+        [sss,hostname] = unix('hostname');
+        hostname = hostname(1:end-1);
+        if strcmp(hostname,'saturnus') | strcmp(hostname,'polaris') | strcmp(hostname,'neptunus')
+          [al,fval,status] = glpkmex_redhat(smin,f,[A; Aeq],[b; beq],ctype,lb,ub,vartype,params,lpsolver);
+        else
+          [al,fval,status] = glpkmex(smin,f,[A; Aeq],[b; beq],ctype,lb,ub,vartype,params,lpsolver);
+        end
+      else
+        [al,fval,status] = linprog(f,A,b,Aeq,beq,lb,ub);
+      end
+      alpha(1:k) = al(1:k) - al(k+1:2*k);
+      if is_w0,
+        alpha(k+1) = al(2*k+1);
+      end
+
+    otherwise
+      disp(type)
+      error ('Wrong type.');
+    end
+  end
+
+  % Choose features
+  ss = sum(abs(alpha(1:k)));
+  J  = find(abs(alpha(1:k)) > ss*prec);
+
+  if isempty(J) | (status <= 0) | (status > 181 | status == 150),
+    prwarning(1,'Fisher classifier is trained.');
+    W1 = fisherc(a);
+    W2 = W1;
+    W3 = featsel(k,[1:k]);
+  else
+    W3 = featsel(k,J);
+    w = [Y; 1] .* alpha(1:k+1);
+    W2 = affine(w(J),w(k+1),a(:,J),lablist,k,2);
+    W2 = cnormc(W2,a(:,J));
+    W1 = W3*W2;
+    W1 = setname(W1,'FSELLPC');
+    W2 = setname(W2,'FSELLPC');
+  end
+return;
Index: /prextra/prextra/gatem.m
===================================================================
--- /prextra/prextra/gatem.m	(revision 2)
+++ /prextra/prextra/gatem.m	(revision 2)
@@ -0,0 +1,120 @@
+%GATEM  Gate classifier
+%
+%       W = GATEM(A,WGATE,W)
+%
+% Make a gate classifier WGATE that splits the feature space in subsets,
+% train a classifier W{i} on each of the subsets and combine the results
+% again. Thus W is a cell-array (cell-vector) of untrained classifiers.
+% Both the gate classifier as the secondary classifiers are trained on
+% dataset A.  This gatem is like a mixture of experts, where each
+% classifier is focussing on a special area in the feature space.
+%
+% Example:  a = gendatb;    % generate a 2-class banana dataset
+%           w = gatem(a, ldc, {parzenc, qdc});
+%                           % split the feature space in two by an ldc
+%                           % train a parzenc in one half, and a qdc in
+%                           % the other half
+%           scatterd(a); plotc(w)
+%
+% In this version the ordering of the data that is passed to the
+% classifiers W{i} is based on the ordering as defined by
+% getlabels(A*WGATE).  Probably this should be made more flexible in a
+% future version.
+
+% Copyright: D.M.J. Tax, D.M.J.Tax@prtools.org
+% Faculty EWI, Delft University of Technology
+% P.O. Box 5031, 2600 GA Delft, The Netherlands
+
+function w = gatem(a,wgate,W)
+
+if nargin < 1 | isempty(a) 
+	% When no inputs are given, we are expected to return an empty
+	% mapping:
+	W = mapping(mfilename,{wgate,W});
+	W = setname(W,'Gate classifier');
+	return
+end
+
+if ~strcmp(getmapping_file(wgate),mfilename)  % training
+	% training consist of training the gate classifier, splitting the
+	% dataset A according to the output label, and training the
+	% individual classifiers on the smaller datasets.
+
+	% First train the gate:
+	if ~istrained(wgate)
+		wgate = a*wgate;
+	end
+	gatelabels = getlabels(wgate);
+	[k,c] = size(wgate);
+
+	% Do some checking
+	if length(W)~=c
+		error('The number of classifiers in W should be equal to the number of classes');
+	end
+
+	% Now map the data through the classifier:
+	out = a*wgate;
+	% and extract the indices of the sub-datasets:
+	[dummy,I] = max(+out,[],2);
+
+	% Train the secondary mappings:
+	for i=1:c
+		%Empty classifiers and trained classifiers don't have to be trained
+		if ~isempty(W{i}) & ~istrained(W{i}) 
+			% extract the datasets,
+			b = a(find(I==i),:);
+			% be careful, this can go wrong:
+			if isempty(b)
+				warning('One of the classifiers did not get any training data');
+			else
+				% and train the mapping
+				W{i} = b*W{i};
+			end
+		end
+		
+		% Does this classifier output the same as the gate??
+		if ~isempty(W{i})
+			labels{i} = matchlablist(getlabels(W{i}),gatelabels);
+		else
+			labels{i} = [];
+		end
+	end
+
+	% Now store everything:
+	V.wgate = wgate;
+	V.W = W;
+	V.labels = labels;
+	w = mapping(mfilename,'trained',V,gatelabels,k,c);
+	w = setname(w,'Gate classifier');
+
+else                  % evaluation
+	% get the data out:
+	V = getdata(wgate);
+	[n,k] = size(a);
+
+	% Apply the gate classifier:
+	out = +(a*V.wgate);
+	[dummy,lab] = max(+out,[],2);
+
+	% Now apply the secondary classifiers
+	for i=1:length(V.W)
+		if ~isempty(V.W{i})
+			% extract the datasets,
+			I = find(lab==i);
+			if ~isempty(I)  % so, if there is data to classify:
+				% compute the output:
+				tmpout = +(a(I,:)*V.W{i});
+				% store the output:
+				%out(I,:) = tmpout(:,V.labels{i});
+				out(I,V.labels{i}) = tmpout;
+			end
+		end
+	end
+
+	% Fill in the data, keeping all other fields in the dataset intact
+	w = setdata(a,out);
+	w = set(w,'featlab',getlabels(wgate),'featsize',getsize_out(wgate));
+
+end
+
+return
Index: /prextra/prextra/genarche.m
===================================================================
--- /prextra/prextra/genarche.m	(revision 2)
+++ /prextra/prextra/genarche.m	(revision 2)
@@ -0,0 +1,261 @@
+%GENARCHE Generate classifier archetypical data
+%
+%		A = GENARCHE(CLASSF,N,SEED,VERSION)
+%
+% INPUT
+%  CLASSF : Classifer (untrained or name in string)
+%  N      : Number of objects per class, default [50 50]
+%  SEED   : Initial state for random generator, default as it is
+%  VERSION: Version, default is most recent one
+%
+% OUTPUT
+%  A      : Dataset
+
+function a = genarche(classf,n,seed,version)
+
+if nargin < 4, version = 0; end
+if nargin < 3, seed = []; end
+if nargin < 2, n = 50; end
+
+if ~isempty(seed)
+	rand('state',seed);
+	randn('state',seed);
+end
+
+if ismapping(classf) & isuntrained(classf)
+	classname = getname(classf);
+elseif isstr(classf)
+	classname = classf;
+else
+		error('Classifier should be given by string or by untrained mapping')
+end
+
+switch classname
+	case {'qdc','QDA'}
+		a = gendath([n n]);
+		a = a * [4 1; -1 4] ./ sqrt(17);
+	case {'udc','UDA'}
+		a = gendath([n n]);
+	case {'ldc','fisherc','LDA'}
+		a = gendatd([n n]);
+	case {'nmc','Nearest-Mean'}
+		a = gendats([n n]);
+	case {'treec','Dec-Tree'}
+		m = genclass(n,[2/3 1/3]);
+		a = [rand(n,1)*3+2 rand(n,1)*5]; 
+		a = [a; rand(m(1),1)*2 rand(m(1),1)*5; rand(m(2),1)*5 rand(m(2),1)+5];
+		a = dataset(a,genlab([n n]),'prior',0);
+	case {'LM-Neural-Net'}
+		m = genclass(n,[2/3 1/3]);
+		a = [rand(n,1)*3+2 rand(n,1)*5]; 
+		a = [a; rand(m(1),1)*2 rand(m(1),1)*5; rand(m(2),1)*5 rand(m(2),1)+5];
+		a = a*[1 1; -1 1];
+		a = dataset(a,genlab([n n]),'prior',0);
+	case {'lmnc','neurc'}
+		state1 = rand('state');
+		state2 = randn('state');
+		rand('state',0); 
+		randn('state',0); 
+		a = gendatb([200,200],0.5); 
+		w = lmnc(a,5,40);              % find a feasible LM network with 5 hu
+		rand('state',state1); 
+		randn('state',state2); 
+        b = gauss(4*n,[10 10],[9 0; 0 9]);
+        b = exp(b/5);
+        b = b - repmat(mean(b)+[1 1],4*n,1);
+		%b = gendatb([2*n 2*n],1);    % generate a too large dataset
+		d = +(b*w);                    % classify
+		u = rand(4*n,1);               % reverse labels at random according posteriors
+		J = find(u < min(+d,[],2));
+		[dd,lab0] = max(+d,[],2);
+        lab = lab0;
+		lab(J) = 3-lab0(J);
+		L1 = find(lab == 1);
+		L2 = find(lab == 2);
+		a = dataset([b(L1(1:n),:);b(L2(1:n),:)],genlab([n n]),'prior',0);
+	case {'parzenc','Parzen'}
+		a = gendatb([n n],2);
+	case {'knnc','K-NN'}
+		a = gendatb([n n],2);
+		a = exp(a/7);
+	case {'nnc','1nnc','spirals','1-NN'}
+		x = rand(n,1)*4*pi+0.25*pi; 
+		a = [x.*x.*sin(x),x.*x.*cos(x)];
+		y = rand(n,1)*4*pi+0.25*pi; 
+		b = [-y.*y.*sin(y),-y.*y.*cos(y)];
+		a = dataset([a;b],genlab([n n]),'prior',0);
+	case {'svc_nu'}
+		% make a linearly separable problem, where one of the classes has
+		% a terribly long moon-shape tail
+		n = round(n/2);
+		S = cat(3,5*eye(2),7*eye(2));
+		extrablob = gauss([n n],[-35 35; 15 0],S);
+		a = [gendatd([n n],2,3.5); extrablob];
+    case {'asym-linear','asym-linear4'}
+        %a = gendatd([2*n 2*n],2,3);
+        b = gendats([n n]);
+        b1 = seldat(b,1) + repmat([0 10],n,1);
+        b2 = seldat(b,2) + repmat([0 -10],n,1);
+        a = gendat([a; b1; b2],[n n]);
+    case {'asym-linear2'}
+				a1 = [rand(n,1) rand(n,1)*10];
+				a2 = [rand(20*n,1)*0.1-1 rand(20*n,1)*10];
+				b1 = [rand(n,1) rand(n,1)*10];
+				b2 = [rand(2*n,1)*100+1 rand(2*n,1)*10];
+				%b3 = [rand(18*n,1)+100 rand(18*n,1)*10];
+				a = [a1;a2;b1;b2]*[1 1;-1 1];
+				a = [a; +gauss(18*n,[40,80],(10*eye(2)))];
+        a = dataset(a,genlab([21*n 21*n]));
+        a = setprior(a,0);
+				a = gendat(a,[n n]);
+    case {'asym-linear3'}
+				a1 = [rand(n,1)*10 rand(n,1)*100];
+				b1 = [rand(n,1)*10+9 rand(n,1)*100+80];
+				a1 = dataset([a1;b1]*[1 1; -1 1],genlab([n n]));
+				b2 = gauss(n,[-100 150],10*eye(2));
+				a2 = gauss(n,[-50 25],10*eye(2));
+				a2 = dataset([a2;b2],genlab([n n]));
+        a = setprior([a1;a2],0);
+				a = gendat(a,[n n]);
+    case {'SVM-1','asym-linear4'}
+				a1 = [[rand(20*n,1)*5 rand(20*n,1)*100]; +gauss(2*n,[8,95],4*eye(2))];
+				b1 = [[rand(20*n,1)*5+11 rand(20*n,1)*100]; +gauss(2*n,[8,5],4*eye(2))];
+				a = [a1;b1]*[2 0; 0 1];
+				%a1 = dataset([a1;b1]*[1 1; -1 1],genlab([n n]));
+				a = dataset(a*[1 1; -1 1],genlab([22*n 22*n]));
+				%b2 = gauss(n,[-20 80],10*eye(2));
+				%a2 = gauss(n,[-90 80],10*eye(2));
+				%a2 = dataset([a2;b2],genlab([n n]));
+        a = setprior(a,0);
+        %a = setprior([a1;a2],0);
+				a = gendat(a,[n n]);
+    case {'Logistic','asym-linear5'}
+				a1 = [gauss(n,[5,50],[1 0; 0 100]); +gauss(n,[-3,25],4*eye(2))];
+				b1 = [gauss(n,[7.5,50],[1 0; 0 100]); +gauss(n,[16,75],4*eye(2))];
+				%a1 = dataset([a1;b1]*[1 1; -1 1],genlab([n n]));
+				a = dataset([a1;b1]*[1 1; -1 1],genlab([2*n 2*n]));
+				%b2 = gauss(n,[-20 80],10*eye(2));
+				%a2 = gauss(n,[-90 80],10*eye(2));
+				%a2 = dataset([a2;b2],genlab([n n]));
+        a = setprior(a,0);
+        %a = setprior([a1;a2],0);
+				a = gendat(a,[n n]);
+    case('line-plane')
+        a = [rand(n,1) 0.25*((rand(n,1)*2).^2)];
+        a = [a; rand(n,2)*[1 0; 0 0.01]+repmat([0 -0.01],n,1)]*[1 1; -1 1];
+        a = dataset(a,genlab([n n]));
+	case {'rbnc','RB-Neural-Net'}
+		a = gendatb([n n]);
+	case 'River'
+		a = genriver([n n],1,0.3);
+	case {'diamonds','Dis-Rep-LC'}
+		a = rand(n,2);
+		m = genclass(n,0.25*ones(1,4));
+		d = (sqrt(2)-1)/2;
+		b = [];
+		b = [b; randxy(m(1),[-d -d],[0 1])];
+		b = [b; randxy(m(2),[-d 1],[1 d+1])];
+		b = [b; randxy(m(3),[1 0],[1+d 1+d])];
+		b = [b; randxy(m(4),[0 -d],[1+d 0])];
+		a = dataset([a;b],genlab([n n]),'prior',0);
+		a = a*[1 1; 1 -1];
+	case {'rbsvc','RB-SVM'}
+		a = circ(n,1,0.5);
+		m = genclass(n,[2/3 1/3]);
+		b = circ(m(1),sqrt(1.5),1);
+		c = circ(m(2),0.5);
+		a = dataset([a;b;c],genlab([n n]),'prior',0);
+	case 'circ2'
+		a = circ(n,1,0.5);
+		m = genclass(n,[2/3 1/3]);
+		b = circ(m(1),sqrt(1.5),1);
+		c = circ(m(2),0.5);
+		a = dataset([a;b;c],genlab([n n]),'prior',0);
+		a(:,2) = 3*a(:,2);
+		a = a*[1 1; 1 -1];
+	case 'chess4'
+		ma = genclass(n,ones(1,8)/8);
+		mb = genclass(n,ones(1,8)/8);
+		a = randxy(ma(1),[0 0],[1 1]);
+		a = [a; randxy(ma(2),[2 0],[3 1])];
+		a = [a; randxy(ma(3),[1 1],[2 2])];
+		a = [a; randxy(ma(4),[3 1],[4 2])];
+		a = [a; randxy(ma(5),[0 2],[1 3])];
+		a = [a; randxy(ma(6),[2 2],[3 3])];
+		a = [a; randxy(ma(7),[1 3],[2 4])];
+		a = [a; randxy(ma(8),[3 3],[4 4])];
+		b = randxy(mb(1),[1 0],[2 1]);
+		b = [b; randxy(mb(2),[3 0],[4 1])];
+		b = [b; randxy(mb(3),[0 1],[1 2])];
+		b = [b; randxy(mb(4),[2 1],[3 2])];
+		b = [b; randxy(mb(5),[1 2],[2 3])];
+		b = [b; randxy(mb(6),[3 2],[4 3])];
+		b = [b; randxy(mb(7),[0 3],[1 4])];
+		b = [b; randxy(mb(8),[2 3],[3 4])];
+		a = dataset([a;b],genlab([n n]),'prior',0);
+		a = a*[1 1; 1 -1];
+	case 'chess41'
+		ma = genclass(n,ones(1,8)/8);
+		mb = genclass(n,ones(1,8)/8);
+		a = randxy(ma(1),[0 0],[1 1]);
+		a = [a; randxy(ma(2),[2 0],[3 1])];
+		a = [a; randxy(ma(3),[1 1],[2 2])];
+		a = [a; randxy(ma(4),[3 1],[4 2])];
+		a = [a; randxy(ma(5),[0 2],[1 3])];
+		a = [a; randxy(ma(6),[2 2],[3 3])];
+		a = [a; randxy(ma(7),[1 3],[2 4])];
+		a = [a; randxy(ma(8),[3 3],[4 4])];
+		b = randxy(mb(1),[1 0],[2 1]);
+		b = [b; randxy(mb(2),[3 0],[4 1])];
+		b = [b; randxy(mb(3),[0 1],[1 2])];
+		b = [b; randxy(mb(4),[2 1],[3 2])];
+		b = [b; randxy(mb(5),[1 2],[2 3])];
+		b = [b; randxy(mb(6),[3 2],[4 3])];
+		b = [b; randxy(mb(7),[0 3],[1 4])];
+		b = [b; randxy(mb(8),[2 3],[3 4])];
+		a = dataset([a;b],genlab([n n]),'prior',0);
+		a(:,2) = 3*a(:,2);
+		a = a*[1 1; 1 -1];
+	case {'Naive-Bayes'}
+		a = [[randn(n,1)/6-0.5; randn(n,1)/6+0.5] 2*rand(2*n,1)-1];
+		b = [2*rand(2*n,1)-1 [randn(n,1)/6-0.5; randn(n,1)/6+0.5]];
+		a = dataset([a;b],genlab([2*n 2*n]));
+		a = setprior(a,0);
+		a = gendat(a,[n n]);
+	otherwise
+		error(sprintf('%s is not implemented',classname))
+end
+
+if ismapping(classf)
+	a = setname(a,[getname(classf) '']);
+else
+	a = setname(a,classf);
+end
+
+function r = randxy(n,x,y)
+
+if nargin < 1, n = 100; end
+if nargin < 2, x = [0 0]; end
+if nargin < 3, y = x + [1 1]; end
+
+r1 = rand(n,2) .* repmat(y-x,n,1);		
+r = r1 + repmat(x,n,1);	
+
+function x = circ(n,r1,r2)
+
+if nargin < 3, r2 = 0; end
+if nargin < 2, r1 = 1; end
+
+m = ceil(2*n*(r1*r1)/(r1*r1 - r2*r2));
+x = rand(m,2) - repmat([0.5 0.5],m,1);
+x = x*r1*2;
+d = sqrt(sum(x.*x,2));
+J = find(d < r1 & d > r2);
+x = x(J(1:n),:);
+
+function x = gausst(n,u,s,t)
+x = randn(n,1);
+x = t*x.*exp(abs(t)*x);
+x = x - mean(x) + u;
+x = s * x ./ std(x);
+
Index: /prextra/prextra/gendati.m
===================================================================
--- /prextra/prextra/gendati.m	(revision 2)
+++ /prextra/prextra/gendati.m	(revision 2)
@@ -0,0 +1,71 @@
+%GENDATI Create dataset from randomly selected windows in a given image
+%
+%		A = GENDATI(IMAGE,WSIZE,N,LABEL)
+%
+% INPUT
+%   IMAGE  - Image of any dimensionality
+%   WSIZE  - Vector with size of the window
+%   N      - Number windows to be generated
+%   LABEL  - Optional string or number with label for all objects
+%
+% OUTPUT
+%   A      - Dataset of N objects and prod(WSIZE) features
+%
+% DESCRIPTION
+%	Windows of the specified size are arbitrarily positioned in the
+% image, converted to a row vector and stored in the dataset A.
+%
+% If specified, all objects have label LABEL. Otherwise they are
+% unlabeled.
+%
+% SEE ALSO
+% DATASETS
+
+% 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 a = gendati(im,ws,m,label)
+
+if nargin < 4, label = []; end
+if nargin < 3, m = 100; end
+
+n = length(ws);
+imsize = size(im);
+if length(imsize) ~= n
+	error('Window size should have same number of components as image size')
+end
+
+if any(ws > imsize)
+	error('Window size should not be larger than image size')
+end
+
+we = prod(ws);
+n = length(ws);
+window_offset = [0:ws(1)-1]';
+for j=2:n
+	wo = window_offset;
+	for i=1:ws(j)-1
+		wo = wo + prod(imsize(1:j-1));
+		window_offset = cat(j,window_offset,wo);
+	end
+end
+window_offset = window_offset(:)';
+
+N = cell(1,n);
+for j=1:n
+	R = imsize(j)-ws(j)+1;
+	N(j) = {ceil(rand(m,1)*R)};
+end
+N = sub2ind(imsize,N{:});
+
+a = zeros(m,prod(ws));
+for i=1:m
+	a(i,:) = im(N(i)+window_offset);
+end
+
+a = dataset(a);
+a = setfeatsize(a,ws);
+if ~isempty(label)
+	a = setlabels(a,label);
+end
Index: /prextra/prextra/gendatw.m
===================================================================
--- /prextra/prextra/gendatw.m	(revision 2)
+++ /prextra/prextra/gendatw.m	(revision 2)
@@ -0,0 +1,44 @@
+%GENDATW Sample dataset by given weigths
+%
+%   B = GENDATW(A,V,N)
+%
+% INPUT
+%   A    Dataset
+%   V    Vector with weigths for each object in A
+%   N    Number of objects to be generated (default size A);
+%
+% OUTPUT
+%   B    Dataset
+%
+% DESCRIPTION
+% The dataset A is sampled using the weigths in V as a prior distribution.
+
+function b = gendatw(a,v,n)
+
+isdataset(a);
+if nargin < 3, n  = size(a,1); end
+
+v = v./sum(v);
+if any(v<0)
+    error('Weights should be positive');
+end
+
+mins = 0;
+nn = 0;
+while(mins < 2)
+	N = genclass(n,v);
+	L = [];
+	while any(N > 0)
+	    L = [L find(N > 0)];
+	    N = N-1;
+	end
+	b = a(L,:);
+	mins = min(classsizes(b));
+	nn = nn + 1;
+	if nn > 100
+		error('Problems with weighted subsampling. Please enlarge training set.')
+	end
+end
+
+return
+
Index: /prextra/prextra/genriver.m
===================================================================
--- /prextra/prextra/genriver.m	(revision 2)
+++ /prextra/prextra/genriver.m	(revision 2)
@@ -0,0 +1,50 @@
+%GENRIVER Generate 2D 2-class River data
+%
+%		A = GENRIVER(N,L,S)
+%
+% INPUT
+%   N    Number of objects in each of the classes (default: [50 50])
+%   L    Length of the river, default: L = 1 (i.e. 2 bends)
+%   S    Width of uniform class distribution on the river borders,
+%        default: S = 1;
+%
+% OUTPUT
+%   A    Generated dataset
+%
+% Generate a separable 'bending river'-like 2D dataset of two classes.
+% Class priors are equal.
+%
+% EXAMPLE
+%  repset = genriver([10 10],4);
+%  trainset = genriver([],4);
+%  w = proxm(repset,'d');
+%  v = w*fisherc(trainset*w);
+%  figure; scatterd(trainset); plotc(v)
+%
+% SEE ALSO
+% DATASETS, PRDATASETS
+
+% Copyright: A. Harol, 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 a = genriver(N,len,s,d)
+
+	prtrace(mfilename);
+	
+	if nargin < 4, d = 0; end
+	if nargin < 3, s = 1; end
+	if nargin < 2 | isempty(len), len = 1; end
+	if nargin < 1 | isempty(N), N = [50 50]; end
+
+	p = [0.5 0.5];
+	N = genclass(N,p);	
+	x1=len*2*pi*rand(1,N(1));
+	x2=sin(x1)+(s*rand(size(x1)))+d;
+	d1 = dataset([x1',x2'],ones(size(x1')));
+	x3=len*2*pi*rand(1,N(2));
+	x4 = sin(x3)+(-s+s*rand(size(x3)))-d;
+	d2 = dataset([x3',x4'],2*ones(size(x3')));
+	a = setprior([d1;d2],p);
+	
+return
Index: /prextra/prextra/getilab.m
===================================================================
--- /prextra/prextra/getilab.m	(revision 2)
+++ /prextra/prextra/getilab.m	(revision 2)
@@ -0,0 +1,77 @@
+%GETILAB retrieve identifier values given identifier names or indices
+%
+%   LABS=GETILAB(A,REF,J)
+%
+% INPUT
+%    A      input dataset with ilab structure (see ENABLEILAB)
+%    REF    name, cell array with names or numerical indices of
+%           identifier columns
+%    J      indices specifying examples to be considered (optional)
+%
+% OUTPUT
+%    VALUES identifier values. If all identifiers are numerical, matrix is
+%           returned, otherwise cell array.
+%
+% DESCRIPTION
+% GETILAB returns values of identifiers, specified by names or numerical
+% indices. This function requires that the dataset supports
+% named-identifiers (see ENABLEILAB). Optional vector J specifies subset
+% of examples to be considered.
+%
+% SEE ALSO
+% ENABLEILAB, REMOVEILAB, GETILAB, SETILAB, GETILABLIST
+
+% $Id: getilab.m,v 1.4 2005/05/03 13:00:39 serguei Exp $
+
+function labs=getilab(a,ref,J)
+
+  isvalidilab(a);                                                    
+                                                                     
+  if ~exist('ref','var') | (exist('ref','var') & isempty(ref))       
+    % retrieve all identifiers                                       
+    ref=1:length(a.user.ilab.names);                                 
+  end                                                                
+                                                                     
+  if ~exist('J','var')                                               
+    J=1:size(a,1);                                                   
+  end                                                                
+                                                                     
+  % get the column to be returned                                    
+                                                                     
+  col=[];                                                            
+  if isnumeric(ref)                                                  
+    col=ref;                                                         
+  else                                                               
+    if ischar(ref)                                                   
+      ref={ref}; % turn it into cell array                           
+    end                                                              
+    if ~iscell(ref)                                                  
+      error('numerical index, name or cell array of names expected') 
+    end                                                              
+    % we know ref is a cell array                                    
+    for i=1:length(ref)                                              
+      ind=strcmp(a.user.ilab.names,ref{i});                          
+      if sum(ind)==0                                                 
+        error(['identifier label ''' ref{i} ''' not defined!']);     
+      end                                                            
+      col=[col find(ind==1)];                                        
+    end                                                              
+  end                                                                
+                                                                     
+  if max(col)>size(a.user.ilab.names,2)                              
+    error('index exceeds number of defined ident labels');           
+  end                                                                
+                                                                     
+  % retrieve ident labels                                            
+  labs=a.ident(J,col);                                               
+  % if possible, convert the ilabs into numerical output             
+  try                                                                
+    labs = cell2mat(labs);                                           
+    % check if the cell array contained single objects               
+    if any(size(labs) ~= [length(J), length(col)]);                  
+      labs = a.ident(J,col);                                         
+    end                                                              
+  catch                                                              
+  end                                                                
+                                                                     
+  return                                                             
Index: /prextra/prextra/getilablist.m
===================================================================
--- /prextra/prextra/getilablist.m	(revision 2)
+++ /prextra/prextra/getilablist.m	(revision 2)
@@ -0,0 +1,27 @@
+%GETILABLIST retrieve list of named identifiers
+%
+%   S=GETILABLIST(A)
+%
+% INPUT
+%    A      input dataset with ilab structure (see ENABLEILAB)
+%
+% OUTPUT
+%    S      cell-array with identifier names.
+%
+% DESCRIPTION
+% GETILABLIST returns a list of identifiers defined in a dataset.  The
+% dataset must support the user-defined named identifiers (see
+% ENABLEILAB).
+%
+% SEE ALSO
+% ENABLEILAB, REMOVEILAB, GETILAB, SETILAB, GETILABLIST
+
+% $Id: getilablist.m,v 1.1 2005/05/02 13:05:46 pavel Exp $
+
+function il=getilablist(a)
+    
+    isvalidilab(a);
+    
+    il=a.user.ilab.names;
+    
+    return;
Index: /prextra/prextra/getmeta.m
===================================================================
--- /prextra/prextra/getmeta.m	(revision 2)
+++ /prextra/prextra/getmeta.m	(revision 2)
@@ -0,0 +1,59 @@
+%GETMETA retrieve metadata in a dataset by key names
+%
+%   [VALUES,KEYS]=GETMETA(A,KEYS)
+%   [VALUES,KEYS]=GETMETA(A)
+%
+% INPUT
+%    A       input dataset with meta structure (see MAKEMETA)
+%    KEYS    key (string) or a cell array of keys (optional)
+%
+% OUTPUT
+%    VALUES  values to be returned. Multiple values are returned in
+%            a cell array.
+%    KEYS    a cell-array with all keys
+%
+% DESCRIPTION
+% Retrieve values for the given keys. If no keys are specified,
+% all values are returned.
+%
+% SEE ALSO
+% REMOVEMETA, GETMETA, SETMETA
+
+% $Id: getmeta.m,v 1.2 2005/05/02 14:10:08 pavel Exp $
+
+function [out,keys]=getmeta(a,inkeys)
+    
+    isvalidmeta(a);
+
+    if ~exist('inkeys','var')
+	out=a.user.meta.values;
+	if nargout==2, keys=a.user.meta.keys; end
+	return
+    end
+    
+    u=a.user;
+    keys=u.meta.keys;
+    values=u.meta.values;    
+    
+    if ischar(inkeys)
+	inkeys={inkeys}; % turn it into cell array
+    end
+    if ~iscell(inkeys)
+	error('name or a cell array of names expected as keys')
+    end
+    % we know inkeys is a cell array
+    out={};
+    for i=1:length(inkeys)
+	ind=strcmp(a.user.meta.keys,inkeys{i});
+	if sum(ind)==0
+	    error(['key ''' inkeys{i} ''' not defined!']); 
+	end
+	ind=find(ind);
+	out{i}=values{ind};
+    end
+    
+    if length(out)==1
+	out=out{1};
+    end
+
+    return
Index: /prextra/prextra/hausdm.m
===================================================================
--- /prextra/prextra/hausdm.m	(revision 2)
+++ /prextra/prextra/hausdm.m	(revision 2)
@@ -0,0 +1,101 @@
+%HAUSDM Hausdorff distance between datasets of image blobs
+%
+%	[DH,DM] = hausdm(A,B,fid)
+%
+% Computes a Hausdorff distance matrix between the sets of binary images 
+% A and B, or datasets containing them as features. 
+% Progress is reported in fid (fid = 1: on the sreeen).
+% If A and B are image sets and 
+% size(A) is [may,max,na]
+% size(B) is [mby,mbx,nb]
+% then DH is the Hausdorff distance matrix, size(DH) = [na,nb], between 
+% these sets. DM is the Modified Hausdorff distance matrix.
+% Preferably na <= nb (faster).
+%
+% See M.-P. Dubuisson and A.K. Jain, Modified Hausdorff distance for object
+%     matching, Proceedings 12th IAPR International Conference on Pattern
+%     Recognition (Jerusalem, October 9-13, 1994), vol. 1, IEEE, Piscataway,
+%     NJ, USA,94CH3440-5, 1994, 566-568.
+%
+% $Id: hausdm.m,v 1.1 2005/04/25 05:54:36 duin Exp $
+
+% Copyright: R.P.W. Duin, duin@ph.tn.tudelft.nl
+% Faculty of Applied Physics, Delft University of Technology
+% P.O. Box 5046, 2600 GA Delft, The Netherlands
+
+function [dh,dm] = hausdm(A,B,fid)
+if nargin < 3, fid = 0; end
+
+if isdataset(A) & isdataset(B)
+	[dh,dm] = hausdm(data2im(A),data2im(B),fid);
+	dh = setdata(A,dh);
+	dm = setdata(A,dm);
+	return
+end
+
+[ma1,ma2,na] = size(A);
+[mb1,mb2,nb] = size(B);
+dh = zeros(na,nb);
+dm = zeros(na,nb);
+for i=1:na
+	a = A(:,:,i); 
+	J = find(any(a));
+	J = [min(J):max(J)];
+	K = find(any(a'));
+	K = [min(K):max(K)];
+	a = double(a(K,J));
+	if length(a(:)) > 0
+		a = bord(a,0);
+	end
+	ca = contourc(a,[0.5,0.5]);
+	J = find(ca(1,:) == 0.5);
+	ca(:,[J J+1]) =[];
+	ca = ca - repmat([1.5;1.5],1,size(ca,2));
+	ca = ca/max(ca(:));
+	ca = ca - repmat(max(ca,[],2)/2,1,size(ca,2));
+	for j = 1:nb
+		b = B(:,:,j);
+		J = find(any(b));
+		J = [min(J):max(J)];
+		K = find(any(b'));
+		K = [min(K):max(K)];
+		b = double(b(K,J));
+		if length(b(:)) > 0
+			b = bord(b,0);
+		end
+		cb = contourc(b,[0.5,0.5]);
+		J = find(cb(1,:) == 0.5);
+		cb(:,[J J+1]) =[];
+		cb = cb - repmat([1.5;1.5],1,size(cb,2));
+		cb = cb/max(cb(:));
+		cb = cb - repmat(max(cb,[],2)/2,1,size(cb,2));
+		dab = sqrt(distm(ca',cb'));
+		dh(i,j) = max(max(min(dab)),max(min(dab')));
+		dm(i,j) = max(mean(min(dab)),mean(min(dab')));
+%		if fid, disp([i,j,dh(i,j),dm(i,j)]); end
+		fprintf(fid,'%5d %5d %10.3f %8.3f \n',i,j,dh(i,j),dm(i,j));
+	end
+end
+	
+% 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);
+%ipcontr(0);
+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: /prextra/prextra/im_update.m
===================================================================
--- /prextra/prextra/im_update.m	(revision 2)
+++ /prextra/prextra/im_update.m	(revision 2)
@@ -0,0 +1,164 @@
+function im_update(arg1)
+% im_update(arg1)
+%
+% Auxiliary function for impaint.m It takes care for the callbacks when
+% the mouse button is clicked or released.
+%
+% See also  impaint.m
+
+    fig = gcbf;
+    allUD = get(fig,'userdata');
+    UD=allUD.impaint; % our portion of UD
+
+    switch lower(arg1)
+     case 'fill'
+      % Check if it is inside the image:
+      set(fig,'units','normalized');
+      pos = get(fig,'currentpoint');
+      if pos(1)<UD.pos(1) | pos(1)>UD.pos(3) | pos(2)<UD.pos(2) | pos(2)>UD.pos(4)
+	  return
+      end
+      % Compute the position:
+      pos(1) = round((pos(1)-UD.pos(1))*UD.W);
+      pos(2) = round(UD.sz(1) - (pos(2)-UD.pos(2))*UD.H);
+
+      switch get(fig,'SelectionType')
+       case 'normal'    % left click
+			%disp('left click');
+			t=floor((UD.brushsz-1)/2);
+			I = -t:t;
+			% just paint it foreground:
+			cdata = get(UD.msk_im,'cdata');
+			cdata(pos(2)+I,pos(1)+I,1) = UD.msk_fg(UD.currlab,1);
+			cdata(pos(2)+I,pos(1)+I,2) = UD.msk_fg(UD.currlab,2);
+			cdata(pos(2)+I,pos(1)+I,3) = UD.msk_fg(UD.currlab,3);
+			set(UD.msk_im,'cdata',cdata);
+			UD.labels(pos(2)+I,pos(1)+I) = UD.currlab;
+       case 'alt'       % right click
+			%disp('right click');
+			% just paint it background:
+			t=floor((UD.brushsz-1)/2);
+			I = -t:t;
+			cdata = get(UD.msk_im,'cdata');
+			cdata(pos(2)+I,pos(1)+I,1) = UD.msk_bg(1);
+			cdata(pos(2)+I,pos(1)+I,2) = UD.msk_bg(2);
+			cdata(pos(2)+I,pos(1)+I,3) = UD.msk_bg(3);
+			set(UD.msk_im,'cdata',cdata);
+			UD.labels(pos(2)+I,pos(1)+I) = 0;
+      end
+      % now I want to trace the cursor:
+      set(fig,'WindowButtonMotionFcn','im_update(''fill'')');
+      
+     case 'stopfill'
+      % stop tracing the cursor:
+      set(fig,'WindowButtonMotionFcn',[]);
+
+     case 'keypress-hidden'
+      t=get(fig,'currentchar');  
+      if strcmp(UD.state,'off') & strcmp(t,'t')
+
+	  UD.state='on';
+
+	  % show the overlay axes:
+	  set(UD.msk_im,'visible','on');
+
+	  % use our impaint handlers (previous are stored)
+	  set(fig,'WindowButtonDownFcn','im_update(''fill'')');
+	  set(fig,'WindowButtonUpFcn','im_update(''stopfill'')');
+	  set(fig,'KeyPressFcn','im_update(''keypress'')');	
+	  
+	  % change the window title
+	  set(fig,'Name',sprintf('imprint on, brush: %d, class: %d',UD.brushsz,UD.currlab));
+	  
+      else
+	  % invoke old handler
+	  eval(UD.oldhandlers.keypress);
+      end
+      
+     case 'keypress'
+      t=get(fig,'currentchar');
+      if double(t)==46
+	  % first allow decimal point
+	  UD.charbuffer=[UD.charbuffer '.'];
+      else
+	  t2=str2num(t); % is is a number?
+	  if ~isempty(t2)
+	      % it's a diggit -> add to the buffer
+	      UD.charbuffer=[UD.charbuffer t];
+	  else
+	      % it's a letter command: execute
+	      % get the current numerical content from buffer:
+	      num=str2num(UD.charbuffer);
+
+	      switch t
+	       case 't' % toggle the impaint on and off
+		if strcmp(UD.state,'on')
+		    UD.state='off';
+
+		    % hide the overlay axes:
+		    set(UD.msk_im,'visible','off');
+		    
+		    % supply original handlers
+		    set(fig,'WindowButtonDownFcn',UD.oldhandlers.btndown);
+		    set(fig,'WindowButtonUpFcn',UD.oldhandlers.btnup);
+		    
+		    % keypress handler must allow to return back
+		    set(fig,'KeyPressFcn','im_update(''keypress-hidden'')');	
+		    
+		    % change the window title
+		    set(fig,'Name','imprint off');
+
+		    set(fig,'UserData',[]);
+		    allUD.impaint=UD;
+		    set(fig,'UserData',allUD);
+		    return
+		    
+		end
+	       case 'b' % change the brush size
+		if isempty(num)
+		    fprintf(1,'impaint: use a numerical prefix followed by letter command!');
+		end
+		if num>0 
+		    UD.brushsz=num;
+		end
+	       case 'c' % class selector
+		if isempty(num)
+		    fprintf(1,'impaint: use a numerical prefix followed by letter command!');
+		end
+		if num>0 
+		    UD.currlab=num;
+		end
+	       case 'a' % alpha level
+		if isempty(num)
+		    fprintf(1,'impaint: use a numerical prefix followed by letter command!');
+		end
+		if num>=0 & num<=1
+		    msk2 = findobj(get(UD.msk_ax,'children'),'type','image');
+		    set(msk2,'alphadata',num);
+		end
+	       case 's' % save labels
+		def='lab';
+		if ~isempty(num), def=[def num2str(num)]; end
+		name=inputdlg('Enter label vector name','Assign labels to a workspace variable',1,{def});
+		if ~isempty(name)
+		    lab=UD.labels;
+		    lab=lab(1:UD.labsize(1),1:UD.labsize(2));
+		    assignin('base',name{1},lab);
+		end
+	       otherwise
+%		fprintf(1,'unknown command');
+	      end
+	      set(fig,'Name',sprintf('imprint on, brush: %d, class: %d',UD.brushsz,UD.currlab));	   
+	      UD.charbuffer=''; % clear the buffer
+	  end
+      end
+    end
+
+    set(fig,'UserData',[]);
+    allUD.impaint=UD;
+    set(fig,'UserData',allUD);
+
+    % DXD it does not seem to work....:
+    %set(UD.msk_im,'erasemode','normal');
+
+    return
Index: /prextra/prextra/impaint.m
===================================================================
--- /prextra/prextra/impaint.m	(revision 2)
+++ /prextra/prextra/impaint.m	(revision 2)
@@ -0,0 +1,218 @@
+function out = impaint(fig,varargin)
+% out = impaint(fig)
+%
+% Create an overlay on the current active image and allow for painting the
+% overlay.
+% Possible calls:
+%   impaint(h,'on');               start the overlay
+%   impaint(h,'off');              remove the whole overlay
+%   impaint(h,classnr,[r g b]);    define a new class with colour [r g b]
+% out = impaint(h,'get');          get the label image out
+%
+% You can leave out h, then the current active figure is used.
+%
+% There are also some interactive changes that you can make:
+%   t    toggles impaint on and off
+%  Nc    selects another class label (prefix this command with the
+%         class number N)
+%  Nb    sets the brush size to N (prefix this command with the brush
+%         size N)
+%   s    save the current labels in Matlab variable 'lab'
+%
+% See also  im_update
+
+
+% When nothing is defined:
+if nargin<2
+	varargin = {'on'};
+end
+% Special case, because an integer can be interpreted as a handle, we
+% assume that when a single input is given, it is actually a label
+if nargin==1
+	if isnumeric(fig)
+		varargin{1} = fig;
+		fig = gcf;
+	end
+end
+if nargin<1
+	fig = gcf;
+end
+% make sure that the first input is the figure handle:
+if ~ishandle(fig)
+	varargin = {fig varargin{:}};
+	fig = gcf;
+end
+
+% Check if we want to get/set the active class, or (re)define a class:
+if isnumeric(varargin{1})
+	allUD = get(fig,'userdata');
+	if ~isfield(allUD,'impaint')
+		error('This image does not have an overlay defined.');
+	end
+	UD=allUD.impaint; % our part of UD
+	if length(varargin)==1  % we want to set/get the active class
+		if varargin{1}>size(UD.msk_fg,1)
+			error('This class is not defined in the overlay.');
+		else
+			UD.currlab = varargin{1};
+		end
+		% and (maybe more important) get the labels to the output
+		if nargout>0
+			out = UD.labels;
+		end
+	else                    % we want to set the colours
+		% We can process a whole list of colours...
+		for i=1:2:(length(varargin)-1)
+			[nr,clr] = checkimpaintcolor(varargin{i},varargin{i+1});
+			UD.msk_fg(nr,:) = clr;
+			% Definition of the active label
+			UD.currlab = nr;
+		end
+	end
+	% store everthing in the figure:
+	set(fig,'UserData',[]);
+	allUD.impaint=UD;
+	set(fig,'UserData',allUD);
+	return;
+end
+
+
+% Now check all the 'string' possibilities:
+switch lower(varargin{1})
+case 'on'
+
+	% Set up variables:
+	allUD=get(fig,'UserData'); % existing user data
+	UD = [];   % Here I store the UserData
+
+	% Put the image on screen
+	UD.im_ax = findobj(get(fig,'children'),'type','axes');
+	UD.im_im = findobj(get(UD.im_ax,'children'),'type','image');
+	% get the size of the image
+	sz = size(get(UD.im_im,'cdata'));
+
+	%set(fig,'busyaction','cancel','DoubleBuffer','on');
+
+	% Label storage:
+	UD.labels = zeros(sz);  % background has label=0
+	UD.labsize=size(UD.labels);
+	
+	% Mask image
+	if length(sz)==2
+		mskim = zeros([sz 3]);
+	elseif length(sz)==3
+		mskim = zeros(sz);
+	else
+		error('I cannot make a mask image for this image data.');
+	end
+	%mskim(1:50, 1:140, 1) = 1; % just to show something...
+	% Definitions of the fore- and background colour
+	if length(varargin)==3
+	        [nr,clr] = checkimpaintcolor(varargin{2},varargin{3});
+		UD.msk_fg(nr,:) = clr;
+		% Definition of the active label
+		UD.currlab = nr;
+	else
+		UD.msk_fg = [1 0 0; 0 1 0; 0 0 1];   % standard: the first class is red.
+		% Definition of the active label
+		UD.currlab = 1;
+	end
+	UD.msk_bg = [0 0 0];
+	% Definition of the brush-size
+	UD.brushsz = 4;
+	% define character command buffer:
+	UD.charbuffer='';
+	UD.state='on'; % impaint is toggled on
+
+	% make a second axis for the overlay:
+	UD.pos = get(UD.im_ax,'position');
+	UD.msk_ax = axes('position', UD.pos);
+	set(fig,'currentaxes',UD.msk_ax);
+	UD.msk_im = imagesc(mskim);
+	set(UD.msk_ax,'visible','off');
+	msk2 = findobj(get(UD.msk_ax,'children'),'type','image');
+	set(msk2,'alphadata',0.5);
+
+	% set the window title:
+	set(fig,'Name',sprintf('imprint on, brush: %d, class: %d',UD.brushsz,UD.currlab));
+	
+	% to speedup coordinate transformations later:
+	UD.sz = sz;
+	UD.W = sz(2)./UD.pos(3);
+	UD.H = sz(1)./UD.pos(4);
+	UD.pos(3) = UD.pos(3)+UD.pos(1);
+	UD.pos(4) = UD.pos(4)+UD.pos(2);
+
+	% now we have to define a callback
+	oldh.btndown=get(fig,'WindowButtonDownFcn');
+	set(fig,'WindowButtonDownFcn','im_update(''fill'')');
+	oldh.btnup=get(fig,'WindowButtonUpFcn');
+	set(fig,'WindowButtonUpFcn','im_update(''stopfill'')');
+	oldh.keypress=get(fig,'KeyPressFcn');
+	set(fig,'KeyPressFcn','im_update(''keypress'')');	
+
+	% preserve previously set handlers
+	UD.oldhandlers=oldh;
+	
+	% store everthing in the figure:
+	set(fig,'UserData',[]);
+	allUD.impaint=UD; % set our substructure in UD
+	set(fig,'UserData',allUD);
+	
+	out=fig;
+	return
+	
+case 'get'
+	allUD = get(fig,'userdata');
+	if ~isfield(allUD,'impaint')
+		error('This image does not have an overlay defined.');
+	end
+	out = allUD.impaint.labels;
+	
+	return
+
+case 'off'
+	% Remove all the traces from the image:
+	allUD = get(fig,'userdata');
+	if ~isfield(allUD,'impaint')
+		error('This image does not have an overlay defined.');
+	end
+	UD=allUD.impaint;
+	delete(UD.msk_ax);
+
+	% restore old handlers
+	set(fig,'WindowButtonDownFcn',UD.oldhandlers.btndown);
+	set(fig,'WindowButtonUpFcn',UD.oldhandlers.btnup);
+	set(fig,'KeyPressFcn',UD.oldhandlers.keypress);	
+	
+	rmfield(allUD,'impaint');
+	set(fig,'userdata',allUD);
+
+	return
+
+end
+
+return
+
+function [nr,clr] = checkimpaintcolor(nr,clr)
+if ~isnumeric(nr)
+	error('The class number should be a scalar.');
+end
+if (length(nr)~=1)
+	error('The class number should be a scalar.');
+end
+if ~isnumeric(clr)
+	error('The color definition should be 1x3 matrix.');
+end
+if (length(clr)~=3)
+	error('The color definition should be 1x3 matrix.');
+end
+if any(clr<0)
+	error('Colour values should be larger than 0.');
+end
+if any(clr>1)
+	error('Colour values should be smaller than 1.');
+end
+% make it a nice row vector:
+clr = clr(:)';
+return
Index: /prextra/prextra/impaintbrush.m
===================================================================
--- /prextra/prextra/impaintbrush.m	(revision 2)
+++ /prextra/prextra/impaintbrush.m	(revision 2)
@@ -0,0 +1,20 @@
+function impaintbrush(h,sz)
+% impaintbrush(h,sz)
+%
+% Set the brush size as it is used in impaint.m.
+
+allUD = get(h,'UserData');
+
+if ~isfield(allUD,'impaint')
+	error('No overlay is defined.');
+end
+UD = allUD.impaint;
+
+UD.brushsz = sz;
+
+allUD.impaint = UD;
+set(h,'UserData',allUD);
+
+return
+
+
Index: /prextra/prextra/impatch.m
===================================================================
--- /prextra/prextra/impatch.m	(revision 2)
+++ /prextra/prextra/impatch.m	(revision 2)
@@ -0,0 +1,323 @@
+function lab = impatch(fig,action,varargin)
+%IMPATCH Create and change a polygon label overlay on an image
+%
+%      OUT = IMPATCH(H,ACTION,VARARGIN)
+%
+% Draw some patches as overlay over an image, given by the handle H.
+% Possible actions are:
+%
+%  Press left button to draw lines, right button to end the patch.
+% >> impatch(h,'label',2,[0 0 1])    defines label 2 in blue
+% >> impatch(h,'label',3,'g')        defines label 2 in blue
+% >> impatch(h,'label',1)            set active label to 1
+% >> out = impatch(h,'get')          get the labels out.
+%  Press ESC to cancel drawing and deleting the polygon.
+%
+% When H is not given, it will take the current active figure.
+%
+% Example (when you have dipimage in your path):
+%   im = laplace;
+%   impaint                % now click somewhere in the image,
+%                           % use leftclick to add more lines
+%                           % use rightclick to close the polygon
+%   impaint('label',2,'b') % make a second label, in blue
+%   impaint('label',3,[0 1 0]) %make third label in green
+%   out = impaint('get')   % get the label image out
+
+% Copyright: D.M.J. Tax, D.M.J.Tax@prtools.org
+% Faculty EWI, Delft University of Technology
+% P.O. Box 5031, 2600 GA Delft, The Netherlands
+
+if nargin<2
+	action = [];
+end
+if nargin<1
+	fig = gcbf;
+end
+% If the figure handle is not given, add it and move all other
+% arguments:
+if ischar(fig)
+	varargin = {action varargin{:}};
+	action = fig;
+	fig = gcbf;
+end
+if ~ishandle(fig) & isnumeric(fig) % we forgot to put 'label' ???
+	varargin = {fig action varargin{:}};
+	action = 'label';
+	fig = gcf;
+end
+
+
+% First check if the figure has already some patches:
+if isempty(fig)
+	fig = gcf;
+	UD = get(fig,'userdata');
+	if isempty(UD) | ~isfield(UD,'impatch')
+		% OK, now we have to initialize it...
+		% maybe for speedup??:
+		set(fig,'busyaction','cancel','doublebuffer','on');
+		UD.impatch = [];
+		% initialize
+		UD.impatch.mode = 'start';
+		UD.impatch.curr = 0;  % current polygon
+		UD.impatch.currlab = 1;  % current label
+		UD.impatch.clrs = [1 0 0];
+		UD.impatch.nrp = 0;   % no polygons defined
+		UD.impatch.h = [];
+		UD.impatch.hlab = [];
+		% take care of the position and size
+		% check if there is an axes object:
+		hh = get(fig,'children');
+		for i=length(hh)
+			isax(i) = strcmp('axes',get(hh(i),'type'));
+		end
+		if sum(isax)==0
+			error('I cannot find an axes!');
+		end
+		if sum(isax)>1
+			error('There are several axes defined:-(');
+		end
+		hax = hh(find(isax));
+		UD.impatch.units = get(hax,'units');
+		set(hax,'units','normalized');
+		UD.impatch.pos = get(hax,'position');
+		him = get(hax,'children');
+		% for recomputing the coordinates:
+		UD.impatch.sz = size(get(him,'cdata'));
+		UD.impatch.W = UD.impatch.sz(2)./UD.impatch.pos(3);
+		UD.impatch.H = UD.impatch.sz(1)./UD.impatch.pos(4);
+		UD.impatch.pos(3) = UD.impatch.pos(3)+UD.impatch.pos(1);
+		UD.impatch.pos(4) = UD.impatch.pos(4)+UD.impatch.pos(2);
+
+		% take care for the callback function:
+		%get(fig,'windowbuttondownfcn')  %DXD
+		UD.impatch.olddownfcn = get(fig,'windowbuttondownfcn');
+		UD.impatch.oldmotionfcn = get(fig,'windowbuttonmotionfcn');
+		UD.impatch.oldkeypressfcn = get(fig,'keypressfcn');
+		set(fig,'userdata',UD);
+		
+		set(fig,'WindowButtonDownFcn','impatch');
+		set(fig,'KeyPressFcn','impatch(''character'')');
+	end
+end
+
+% The special case that we get direct actions, instead of mouse events:
+if ~isempty(action)
+	switch action
+	case 'get'     % we want to get the label image out
+		% set up:
+		UD = get(gcf,'userdata');
+		% go over all the patches:
+		lab = zeros(UD.impatch.sz);
+		[X,Y] = meshgrid(1:UD.impatch.sz(2), 1:UD.impatch.sz(1));
+		for i=1:length(UD.impatch.h)
+			vx = get(UD.impatch.h(i),'xdata');
+			vy = get(UD.impatch.h(i),'ydata');
+			im = inpolygon(X,Y,vx,vy);
+			I = find(im);
+			lab(I) = UD.impatch.hlab(i);
+			%lab = im;
+		end
+		return
+	case 'label'   % add a new label or change active label
+		% set up:
+		UD = get(gcf,'userdata');
+		newnr = varargin{1};
+		if newnr>size(UD.impatch.clrs,1)  % add a new label
+			if length(varargin)<2
+				error('I am expecting now a colour definition');
+			end
+			if (newnr>size(UD.impatch.clrs,1)+1)
+				newnr=size(UD.impatch.clrs,1)+1;
+			end
+			% check if the color definition is acceptable:
+			newclr = checkcolour(varargin{2});
+			% and commit:
+			UD.impatch.currlab = newnr;
+			UD.impatch.clrs(newnr,:) = newclr;
+		else                             % change the active label
+			if newnr>0
+				UD.impatch.currlab = newnr;
+			end
+		end
+		% Save the results:
+		set(gcf,'userdata',UD);
+		return;
+	case 'off'  % remove all traces from the figure
+		UD = get(gcf,'userdata');
+		if ~isempty(UD) & isfield(UD,'impatch')
+			set(gcf,'windowbuttondownfcn',UD.impatch.olddownfcn);
+			set(gcf,'windowbuttonmotionfcn',UD.impatch.oldmotionfcn);
+			set(gcf,'keypressfcn',UD.impatch.oldkeypressfcn);
+			delete(UD.impatch.h);
+			UD = rmfield(UD,'impatch');
+			if size(fieldnames(UD),1)==0, UD = []; end
+			set(gcf,'userdata',UD);
+			return
+		end
+	case 'character'   % A user pressed a key
+		ch = double(get(gcf,'currentcharacter'));
+		if isempty(ch), return, end
+		UD = get(gcf,'userdata');
+		switch ch
+		case 27 % ESCAPE
+%			switch UD.impatch.mode
+%			case 'draw'  % cancel this complete line
+				UD.impatch.hlab(UD.impatch.curr) = [];
+				delete(UD.impatch.h(UD.impatch.curr));
+				UD.impatch.h(UD.impatch.curr) = [];
+				UD.impatch.nrp = UD.impatch.nrp-1;
+				UD.impatch.curr = 0;
+				UD.impatch.mode = 'start';
+				set(gcf,'userdata',UD);
+				set(fig,'windowbuttonmotionfcn',UD.impatch.oldmotionfcn);
+%			case 'move'  % cancel the movement of the node
+%			end
+		end
+		return
+
+	end
+end
+
+% Get the userdata
+UD = get(fig,'userdata');
+% Check if the cursor position is inside the image:
+set(fig,'units','normalized');
+pos = get(fig,'currentpoint');
+if pos(1)<UD.impatch.pos(1) | pos(1)>UD.impatch.pos(3) | ...
+	pos(2)<UD.impatch.pos(2) | pos(2)>UD.impatch.pos(4)
+	return
+end
+% Compute the position:
+pos(1) = round((pos(1)-UD.impatch.pos(1))*UD.impatch.W);
+pos(2) = round(UD.impatch.sz(1) - (pos(2)-UD.impatch.pos(2))*UD.impatch.H);
+
+% Now do the actions that need the position of the cursor:
+if strcmp(action,'move');
+	switch UD.impatch.mode
+	case 'draw'
+		h = UD.impatch.h(UD.impatch.curr);
+		xd = get(h,'xdata'); yd = get(h,'ydata');
+		xd(end) = pos(1);    yd(end) = pos(2);
+		set(h,'xdata',xd);   set(h,'ydata',yd);
+	case 'move' % we want to move a node...
+		minnode = UD.impatch.currnode;
+		hh = UD.impatch.h(UD.impatch.curr);
+		xd = get(hh,'xdata'); yd = get(hh,'ydata');
+		xd(minnode) = pos(1); yd(minnode) = pos(2);
+		set(hh,'xdata',xd);   set(hh,'ydata',yd);
+	end
+	% Save the results:
+	set(fig,'userdata',UD);
+	return
+end
+
+% Otherwise we have to check mouse buttons
+switch UD.impatch.mode
+case 'start'   % start a new line
+	switch get(fig,'selectiontype')
+	case 'alt'
+		if length(UD.impatch.h)==0
+			return;
+		end
+		% Find nearest node
+		mind = inf; minp = 0; minnode = 0;
+		for i=1:length(UD.impatch.h)
+			hh = UD.impatch.h(i);
+			coord = [get(hh,'xdata') get(hh,'ydata')];
+			N = size(coord,1);
+			diff = repmat(pos,N,1) - coord;
+			[diff,nodenr] = min(sum(diff.*diff,2));
+			if diff<mind
+				mind = diff;
+				minp = i;
+				minnode = nodenr;
+			end
+		end
+		% If we are close enough, we should move the node...
+		if mind<100  % magic distance
+			UD.impatch.curr = minp;
+			UD.impatch.currnode = minnode;
+			hh = UD.impatch.h(minp);
+			xd = get(hh,'xdata'); yd = get(hh,'ydata');
+			xd(minnode) = pos(1); yd(minnode) = pos(2);
+			set(hh,'xdata',xd);   set(hh,'ydata',yd);
+			UD.impatch.mode = 'move';
+			set(fig,'WindowButtonMotionFcn','impatch(''move'')');
+		end
+	case 'normal'
+		% Set the number of the current polygon
+		UD.impatch.nrp = UD.impatch.nrp+1;
+		UD.impatch.curr = UD.impatch.nrp;
+		% Define the colour for this line
+		thisclr = UD.impatch.clrs(UD.impatch.currlab,:);
+		UD.impatch.h(UD.impatch.curr) = line([pos(1) pos(1)],[pos(2) pos(2)]);
+		UD.impatch.hlab(UD.impatch.curr) = UD.impatch.currlab;
+		set(UD.impatch.h(UD.impatch.curr),'color',thisclr);
+		UD.impatch.mode = 'draw';
+		set(fig,'WindowButtonMotionFcn','impatch(''move'')');
+	end
+	% Save the results:
+	set(fig,'userdata',UD);
+	return
+case 'draw'
+	switch get(fig,'selectiontype')
+	case 'normal' % extend the line with a new node
+		h = UD.impatch.h(UD.impatch.curr);
+		xd = get(h,'xdata'); yd = get(h,'ydata');
+		xd(end+1) = pos(1);    yd(end+1) = pos(2);
+		set(h,'xdata',xd);   set(h,'ydata',yd);
+	case 'alt'    % end the line, and fill the polygon
+		h = UD.impatch.h(UD.impatch.curr);
+		xd = get(h,'xdata'); yd = get(h,'ydata');
+		% get rid of the last point, that was only glued to the pointer:
+		xd(end) = []; yd(end) = [];
+		delete(h);  % will this be OK??
+		myclr = UD.impatch.clrs(UD.impatch.hlab(UD.impatch.curr),:);
+		h = patch(xd,yd, myclr);
+		set(h,'facealpha',0.5);
+		UD.impatch.h(UD.impatch.curr) = h;
+		UD.impatch.curr = 0;
+		UD.impatch.mode = 'start';
+		set(fig,'windowbuttonmotionfcn',UD.impatch.oldmotionfcn);
+	end
+case 'move'
+	minnode = UD.impatch.currnode;
+	hh = UD.impatch.h(UD.impatch.curr);
+	xd = get(hh,'xdata'); yd = get(hh,'ydata');
+	xd(minnode) = pos(1); yd(minnode) = pos(2);
+	set(hh,'xdata',xd);   set(hh,'ydata',yd);
+	UD.impatch.mode = 'start';
+end
+
+% Save the results:
+set(fig,'userdata',UD);
+
+return
+
+function newclr = checkcolour(newclr)
+
+if ischar(newclr)
+	if length(newclr)>1
+		error('I am expecting a single character');
+	end
+	defclrchar = 'rbgymc';
+	defclrs = [1 0 0; 0 0 1; 0 1 0; 1 1 0; 1 0 1; 0 1 1];
+	newclr = find(defclrchar==newclr);
+	if isempty(newclr)
+		error('Not a legal colour definition');
+	end
+	newclr = defclrs(newclr,:);
+else
+	newclr = newclr(:)';
+	if length(newclr)~=3
+		error('The colour definition should have R G and B');
+	end
+	if any(newclr<0)
+		error('Colour values should be larger than 0.');
+	end
+	if any(newclr>1)
+		error('Colour values should be smaller than 1.');
+	end
+end
+return
Index: /prextra/prextra/inckernel.m
===================================================================
--- /prextra/prextra/inckernel.m	(revision 2)
+++ /prextra/prextra/inckernel.m	(revision 2)
@@ -0,0 +1,110 @@
+function K = inckernel(par,I,J);
+%INCKERNEL Kernel definition for incsvdd
+%
+%              K = INCKERNEL(PAR,I,J);
+%
+% Computation of the kernel function for the incremental SVDD. It is
+% assumed that there is a global variable X_incremental, containing the
+% objects. This is to avoid unnecessary overhead for very large
+% datasets. Therefore we will also not use the 'standard' ways to comute
+% the kernel (i.e. proxm).
+%
+% The kernel is defined by PAR.  We assume that PAR is a structure
+% containing:
+%  PAR.type   the kernel type
+%       'kernel'       | 'k': X_incremental(I,J)
+%       'polynomial'   | 'p': sign(xI*xJ'+1).*(xI*xJ'+1).^s
+%       'exponential'  | 'e': exp(-(||xI-xJ||)/s)
+%       'radial_basis' | 'r': exp(-(||xI-xJ||.^2)/(s*s))
+%       'sigmoid'      | 's': sigm((sign(xI*xJ').*(xI*xJ'))/s)
+%       'distance'     | 'd': ||xI-xJ||.^s
+% And the special case:
+%       'inline'       | 'i': use s as an inline function
+%  PAR.s      the free parameter,
+%             when more than 1 free parameter is required, then par.s
+%             can be a vector of length 2 or more.
+% The index vectors I and J indicate between which objects in
+% X_incremental the kernel should be computed.
+%
+% See also: incsvdd, proxm
+
+% Copyright: D. Tax, R.P.W. Duin, davidt@ph.tn.tudelft.nl
+% Faculty of Applied Physics, Delft University of Technology
+% P.O. Box 5046, 2600 GA Delft, The Netherlands
+
+if nargin<2
+	I = [];
+	J = [];
+end
+
+global X_incremental;
+if isempty(X_incremental)
+  error('No data matrix X_incremental defined');
+end
+if isa(X_incremental,'dataset');
+  error('Please make X_incremental a normal matlab array');
+end
+
+A = X_incremental(I,:);
+B = X_incremental(J,:);
+
+switch par.type
+case {'kernel' 'k'}
+	if isempty(I) & isempty(J)
+		K = X_incremental;
+	else
+		K = X_incremental(I,J);
+	end
+case {'polynomial' 'p'}
+	K = A*B'; 
+	[n,d] = size(A);
+	[m,d] = size(B);
+	if par.s ~= round(par.s)
+		K = K + ones(n,m);
+		K = sign(K).*abs(K).^par.s;
+	elseif par.s ~= 1
+		K = K + ones(n,m);
+		K = K.^par.s;
+	end
+case {'sigmoid' 's'}
+	K = A*B'; 
+	if length(par.s)>1
+		K = sigm(K/par.s(1) + par.s(2));   %DXD: I need this sometimes
+	else
+		K = sigm(K/par.s);
+	end
+case {'exponential' 'e'}
+	K = sqeucldistm(A,B);
+	J = find(K<0);
+	K(J) = zeros(size(J));
+	K = exp(-sqrt(K)/par.s);
+case {'radial_basis' 'r'}
+	K = sqeucldistm(A,B);
+	J = find(K<0);
+	K(J) = zeros(size(J));
+	K = exp(-K/(par.s*par.s));
+case {'inv_radial_basis' 'i'}
+	K = sqeucldistm(A,B);
+	J = find(K<0);
+	K(J) = zeros(size(J));
+	K = exp(sqrt(K)/par.s);
+case {'distance' 'd'}
+	K = sqeucldistm(A,B);
+	J = find(K<0);
+	K(J) = zeros(size(J));
+	if par.s ~= 2
+		K = K.^(par.s/2);
+	end
+case {'inline' 'i'}
+	error('Sorry, this does not work yet.');
+case {'combp1s2'}
+	K = sqeucldistm(A,B);
+	J = find(K<0);
+	K(J) = zeros(size(J));
+	K = exp(-K/(2*2));
+	K = par.s*A*B' + (1-par.s)*K; 
+otherwise
+	error(sprintf('Unknown proximity type: %s',par.type))
+end
+
+return
Index: /prextra/prextra/incsvc.m
===================================================================
--- /prextra/prextra/incsvc.m	(revision 2)
+++ /prextra/prextra/incsvc.m	(revision 2)
@@ -0,0 +1,161 @@
+function w = incsvc(a,ktype,par,C)
+%INCSVC Incremental support vector classifier
+%
+%     W = INCSVC(A,KTYPE,PAR,C)
+%
+% INPUT
+%   A      Dataset
+%   KTYPE  Type of the kernel (optional; default: 'p')
+%    PAR   Kernel parameter (optional; default: 1)
+%    C     Regularization parameter (optional; default: 1)
+%
+%  OUTPUT
+%    W     Mapping: Support Vector Classifier
+%
+%  DESCRIPTION
+%  Optimizes a support vector classifier for the dataset A by an
+%  incremental procedure to perform the quadratic programming. The
+%  classifier can be of one of the types as defined by PROXM. Default is
+%  linear (TYPE = 'p', PAR = 1). The kernel computation is done by
+%  INCKERNEL which is more lightweight than PROXM.  C < 1 allows for
+%  more class overlap. Default C = 1.
+%
+%  See also ADD_OBJ_CL, SVC, INCKERNEL
+
+% Copyright: D.M.J. Tax, D.M.J.Tax@prtools.org
+% Faculty EWI, Delft University of Technology
+% P.O. Box 5031, 2600 GA Delft, The Netherlands
+
+if nargin < 4 | isempty(C)
+        C = 1;
+        prwarning(3,'Regularization parameter C set to 1\n');
+end
+if nargin < 3 | isempty(par)
+        par = 1;
+        prwarning(3,'Kernel parameter par set to 1\n');
+end
+if nargin < 2 | isempty(ktype)
+        ktype = 'p';
+        prwarning(3,'Polynomial kernel type is used\n');
+end
+if nargin < 1 | isempty(a)
+        w = mapping(mfilename,{ktype,par,C});
+        w = setname(w,'Inc. Support Vector Classifier');
+        return;
+end
+
+if ~isa(ktype,'mapping') %training
+	%Basic sizes::
+	[N,k,c] = getsize(a);
+	%kernel definition:
+	kernel = 'inckernel';
+	kpar.type = ktype;
+	kpar.s = par;
+	%setting for displaying: (I know you can program it shorter, but this is
+	%more clear):
+	if N>=1000
+		dodisplay = 1;
+	else
+		dodisplay = 0;
+	end
+
+	if c==2
+		% randomize the data:
+		%I = randperm(N); a = a(I,:);
+		a = unrandomize(a);
+
+		% settings for the program:
+		global X_incremental;
+		X_incremental = +a;
+		y = 3 - 2*getnlab(a);
+
+		% here we go:
+		alf = zeros(N,1);   % weights
+		grad = [];          % the gradient of all seen objects
+		setR = [];          % 'rest' set
+		Kr = [];            % kernel matrix of the rest objects(RxS)
+		setD = [];
+		setS = [];
+		setE = [];
+		Ke = [];
+		Ks = 0;
+		b = 0;
+		R = inf;
+		tol = 1e-8;
+
+		% startup:
+		for c=1:N
+			if dodisplay & mod(c,100)==0
+				fprintf('%d/%d ',c,N);
+			end
+			add_obj_cl;
+		end
+
+		% make the classifier
+		W.kernel = kernel;
+		W.kpar = kpar;
+		J = [setS;setE];
+		W.J = J;
+		W.sv = X_incremental(J,:);
+		W.v = y(J).*alf(J,:);
+		W.b = b;
+		w = mapping(mfilename,'trained',W,getlablist(a),k,2);
+		w = setname(w,'Inc. Support Vector Classifier');
+
+	else   % multi-class classifier:
+		
+		w = mclassc(a,mapping(mfilename,{ktype,par,C}));
+		
+	end
+
+else   %execution
+	W = +ktype;
+	[n,d] = size(a);
+	laba = getlab(a);
+	orga = a;
+	a = +a;
+	global X_incremental;
+	X_incremental = [W.sv; zeros(1,d)];
+	nra = size(W.sv,1)+1; I = 1:nra;
+	out = repmat(W.b,n,1);
+	for i=1:n
+		X_incremental(nra,:) = a(i,:);
+		Ka = feval(W.kernel,W.kpar,nra,I);
+		out(i) = out(i) + Ka(1:(nra-1))*W.v;
+	end
+	newout = [out -out];
+	w = setdat(orga,newout,ktype);
+end
+
+return
+
+
+
+function a= unrandomize(a);
+% Unrandomize a dataset;
+
+[n,k,c]=getsize(a);
+if c~=2
+	error('I assume 2 classes');
+end
+
+nlab = getnlab(a);
+I1 = find(nlab==1);
+I2 = find(nlab==2);
+
+if length(I1)<length(I2)
+	cmin = length(I1);
+else
+	cmin = length(I2);
+	tmpI = I1;
+	I1 = I2;
+	I2 = tmpI;
+end
+
+J=[];
+J(1:2:2*cmin) = I1;
+J(2:2:2*cmin) = I2(1:cmin);
+J = [J,I2((cmin+1):end)'];
+a = a(J,:);
+
+return
Index: /prextra/prextra/isvalidilab.m
===================================================================
--- /prextra/prextra/isvalidilab.m	(revision 2)
+++ /prextra/prextra/isvalidilab.m	(revision 2)
@@ -0,0 +1,57 @@
+%ISVALIDILAB validity check for named-identifier support in a dataset
+%
+%   I=ISVALIDILAB(A,FIELDS)
+%
+% INPUT
+%   A       input dataset
+%   FIELDS  (opt) name or cell array with required identifier names 
+%
+% OUTPUT
+%   I       logical output: the dataset is/is not valid
+%
+% DESCRIPTION 
+% This function validates if the dataset contains named identifiers.  IN
+% order to enable the use of named identifiers in a dataset, use
+% ENABLEILAB.  If a name or a cell array with identifier names is supplied
+% in the FIELDS parameter, these are checked for existence.
+%
+% SEE ALSO
+% ENABLEILAB, REMOVEILAB, GETILAB, SETILAB, GETILABLIST
+
+% $Id: isvalidilab.m,v 1.2 2005/05/02 13:05:46 pavel Exp $
+
+function i=isvalidilab(a,fields)
+
+    if ~exist('fields','var'), fields=[]; end
+    if ischar(fields), fields={fields}; end
+    
+    isdataset(a);
+
+    i=0;
+    user=a.user;
+    if isfield(user,'ilab')
+	t=user.ilab;
+	if isstruct(t) & isfield(t,'names') & iscell(t.names)
+	    i=1;
+	end
+    end
+
+    if (nargout == 0) & (i == 0)
+	error(['Dataset with ilab structure in user field expected. Use help enableilab for details.'])
+    end
+    
+    % validate the required fields, if supplied
+    if ~isempty(fields)
+	for j=1:length(fields)
+	    if ~any(strcmp(user.ilab.names,fields{j}))
+		i=0;
+		if (nargout == 0) & (i == 0)
+		    error(['required field ''' fields{j} ''' is not present in a dataset']);
+		end
+	    end
+	end
+    end
+    
+    
+    return
+	
Index: /prextra/prextra/isvalidmeta.m
===================================================================
--- /prextra/prextra/isvalidmeta.m	(revision 2)
+++ /prextra/prextra/isvalidmeta.m	(revision 2)
@@ -0,0 +1,53 @@
+%ISVALIDMETA meta validity check for datasets
+%
+%   I=ISVALIDMETA(A,KEYS)
+%
+% INPUT
+%   A       input dataset
+%   KEYS  (opt) name or cell array with names of required meta keys
+%
+% OUTPUT
+%   I       logical output: the dataset is/is not valid
+%
+% DESCRIPTION
+% This function validates if the dataset contains meta capability i.e.
+% can store additional named sets of identifiers for each object.
+% The metalist containing the names of ident columns is stored in the
+% dataset user field. If a name or a cell array with column names is
+% supplied in KEYS parameter, these are checked for existence.
+
+function i=isvalidmeta(a,inkeys)
+
+    if ~exist('inkeys','var'), inkeys=[]; end
+    if ischar(inkeys), inkeys={inkeys}; end
+    
+    isdataset(a);
+
+    i=0;
+    user=a.user;
+    if isfield(user,'meta')
+	t=user.meta;
+	if isstruct(t) & isfield(t,'keys') & iscell(t.keys) & isfield(t,'values') & iscell(t.values)
+	    i=1;
+	end
+    end
+
+    if (nargout == 0) & (i == 0)
+	error(['Dataset with meta structure in user field expected. Use help enablemeta for details.'])
+    end
+    
+    % validate the required keys, if supplied
+    if ~isempty(inkeys)
+	for j=1:length(inkeys)
+	    if ~any(strcmp(user.meta.keys,inkeys{j}))
+		i=0;
+		if (nargout == 0) & (i == 0)
+		    error(['required key ''' inkeys{j} ''' is not present in a dataset']);
+		end
+	    end
+	end
+    end
+    
+    
+    return
+	
Index: /prextra/prextra/lasso.h
===================================================================
--- /prextra/prextra/lasso.h	(revision 2)
+++ /prextra/prextra/lasso.h	(revision 2)
@@ -0,0 +1,75 @@
+/* Copyright (C) 1998
+   Berwin A Turlach <bturlach@stats.adelaide.edu.au> */
+
+/* This library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Library General Public License
+   as published by the Free Software Foundation; either version 2 of
+   the License, or (at your option) any later version. */
+
+/* This library is distributed in the hope that it will be useful, but
+   WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Library General Public License for more details. */
+
+/* You should have received a copy of the GNU Library General Public
+   License along with this library; if not, write to the Free Software
+   Foundation, Inc., 59 Temple Place, Suite 330, Boston,
+   MA 02111-1307, USA. */
+/*#define S_Plus*/
+#define Matlab
+#ifndef BT_LASSO_H
+#define BT_LASSO_H
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <float.h>
+#if defined(S_Plus)
+#include <S.h>
+#if defined(SPLUS_VERSION) && SPLUS_VERSION >= 4000
+#  include <newredef.h>
+#endif
+#endif
+
+#if defined(Matlab)
+#   include "mex.h"
+#endif
+
+void lasso(double *x, long *pn, long *pm, double *pt,
+           double *beta, double *y, double *yhat1, double *r,
+           double *lagrangian, long *psuc,  long *pverb, long *pas_sub);
+void mult_lasso(double *x, long *pn, long *pm, double *pt, long *pl,
+                double *beta, double *y, double *yhat1, double *r,
+                double *lagrangian, long *psuc, long *pverb);
+void rs_lasso(double *x, double *y, long *pn, long *porder,
+           double *pt, double *beta, double *yhat1,  double *yhat2,
+           double *r, long *psuc, long *pverb, long *pas_sub);
+void rs_auto(double *x, double *y, long *pn, long *porder,
+           double *pt, double *beta, double *yhat1, double *yhat2,
+           double *r,  long *pcrit1, double *pprec,
+           long *psuc, long *pverb,
+           double *tgr, double *tgr_val, long *ptgr_len);
+#if !defined(S_Plus)
+#define Calloc(n,t)       (t *) malloc((n)*sizeof(t))
+#define Realloc(p,n,t)    (t *) realloc(p,(n)*sizeof(t))
+#define Free(p)           free(p), p=NULL
+#include <string.h>
+#define Memcpy(p,q,n)     memcpy(p,q,(n)*sizeof(*(p)))
+#endif
+
+#define TRUE 1
+#define FALSE 0
+/*#define RMAT(i,j) (1.0e-10+*(rmat+(j)*((j)+1)/2+(i))) //??*/
+#define RMAT(i,j) (*(rmat+(j)*((j)+1)/2+(i))) /*??*/
+#define QMAT(i,j) (*(qmat+(j)*q_nrow+(i)))
+#define RLAST(i)  (*(rmat+(r_ncol-1)*r_ncol/2+(i)))
+#define RLTEL     (*(rmat+r_ncol*(r_ncol+1)/2-1))
+#define RCOL(j)   (rmat+(j)*((j)+1)/2)
+#define QCOL(j)   (qmat+(j)*q_nrow)
+#define RLTCOL    (rmat+(r_ncol-1)*r_ncol/2)
+/*#define QR_CHUNK 10*/
+#if defined (S_Plus)
+#define ERRMSG(where,msg) errmsg(msg)
+#else
+#define ERRMSG(where,msg) errmsg(where,msg)
+#endif
+#endif
Index: /prextra/prextra/lasso.m.org
===================================================================
--- /prextra/prextra/lasso.m.org	(revision 2)
+++ /prextra/prextra/lasso.m.org	(revision 2)
@@ -0,0 +1,59 @@
+function [beta]=lasso(X,Y,lambda)
+% function [beta]=lasso(X,Y,lambda)
+% This function minimize 
+% f(beta)=1/N*(Y-X*beta)'*(Y-X*beta)+lambda*sum(abs(beta)) 
+% with respect to beta.
+%
+% X is a matrix of explanatory variables
+% Y is a vector with the response variables
+% lambda is the weight of the Lagrange factor
+% N is the number of observations
+
+      [n,p]=size(X);      
+      X=X/sqrt(n);
+      Y=Y/sqrt(n);
+      Xtwo=2*X'*X;
+      XY=2*X'*Y;
+
+%% Improvement suggested by Sijmen de Jong
+
+      if(p>n)
+        beta=X'*pinv( X * X') * Y;
+      else
+        beta=pinv(Xtwo)*XY;
+      end
+%%
+      notfound=1;
+
+      while(notfound==1)
+        notzero=find(abs(beta)>1e-5);
+        lnz=length(notzero);
+        therest=setdiff(1:p,notzero);
+        if(length(therest>0))
+          beta(therest)=0;
+        end;
+        if(lnz==0)
+          notfound=0;
+        else
+          g=-XY(notzero)+Xtwo(notzero,notzero)*beta(notzero)+lambda*sign(beta(notzero));
+          Z=sparse(1:lnz,1:lnz,0.5./abs(beta(notzero)),lnz,lnz);
+          if(lnz<n)
+            betastep=(Xtwo(notzero,notzero)+lambda*Z)\g;
+          else
+            Ainv=diag(1./diag(2*lambda*Z));
+            betastep=Ainv*g-Ainv*(2*X(:,notzero)'*((X(:,notzero)*Ainv*2*X(:,notzero)'+eye(n))\(X(:,notzero)*(Ainv*g))));
+          end
+          tempsing=find(betastep==Inf);
+          if((max(isinf(betastep))+max(isnan(betastep)))>0)
+            newbeta=ones(size(tempbeta))*NaN;
+            notfound=0;
+          else                                    
+            newbeta=beta(notzero)-betastep;
+            if(max(abs(betastep))<1E-9)
+              notfound=0;
+            end;
+          end
+          beta(notzero)=newbeta;
+        end
+      end
+
Index: /prextra/prextra/lassoc.m
===================================================================
--- /prextra/prextra/lassoc.m	(revision 2)
+++ /prextra/prextra/lassoc.m	(revision 2)
@@ -0,0 +1,65 @@
+function w = lassoc(x, lambda)
+% w = lassoc(x, lambda)
+%prtrace(mfilename);
+
+mustScale=0;
+
+if (nargin < 2)
+	%prwarning(3,'Lambda set to one');
+	lambda = 1;
+end
+if (nargin < 1) | (isempty(x))
+	w = mapping(mfilename,{lambda});
+	w = setname(w,'LASSO classifier');
+	return
+end
+
+if ~ismapping(lambda)   % train the mapping
+
+    % Unpack the dataset.
+    islabtype(x,'crisp');
+    %isvaldset(x,1,2); % at least 1 object per class, 2 classes
+    [n,k,c] = getsize(x);
+
+    % Is this necessary??
+    if mustScale
+        wsc = scalem(x,'variance');
+        x.data = x.data*wsc;
+    end
+
+    if c ~= 2  % two-class classifier:
+        error('Only a two-class classifier is implemented');
+    end
+
+    beta=-lasso(+x,3-2*getnlab(x),lambda);
+
+    % now find out how sparse the result is:
+    nr = sum(abs(beta)>1.0e-8);
+
+    % and store the results:
+    if mustScale
+        W.wsc = wsc;
+    end
+
+    W.beta = beta; % the ultimate weights
+    W.nr = nr;
+    w = mapping(mfilename,'trained',W,getlablist(x),size(x,2),c);
+    w = setname(w,'LASSO classifier');
+
+else
+    % Evaluate the classifier on new data:
+    W = getdata(lambda);
+    n = size(x,1);
+
+    % scaling and linear classifier:
+    if mustScale
+        x = x*W.wsc;
+    end
+    out = x*W.beta;
+
+    % and put it nicely in a prtools dataset:
+    w = setdat(x,sigm([-out out]),lambda);
+
+end
+
+return
Index: /prextra/prextra/lessc.m
===================================================================
--- /prextra/prextra/lessc.m	(revision 2)
+++ /prextra/prextra/lessc.m	(revision 2)
@@ -0,0 +1,137 @@
+function [w,I] = lessc(x, C, ftype, include_bias)
+%LESSC Least Error in Sparse Subspaces classifier
+%
+%     W = LESSC(X, C, FTYPE, INCLUDE_BIAS)
+%
+% Train a linear classifier which also performs feature selection.
+% In this version we do:
+%                min \sum_i w_i + C*delta_i
+%          s.t. forall_i   w^T f(x_i) > 1 - delta_i
+%                       sum_i |w_i| = 1
+% where f(x_i) is in principle free, but as a start we use the nearest
+% mean idea:
+%
+%         f(x_i) = (x-mu2).^2 - (x-mu1).^2
+% See for further definitions of f(x_i) lessfx.
+%
+% Dxd   15-3-2004
+prtrace(mfilename);
+
+if (nargin < 4)
+	% To include a bias term in the model, we extend the number of features
+	% by one:
+	include_bias = 0;
+end
+if (nargin < 3)
+	prwarning(3,'Use default function fx.');
+	ftype = 1; 
+end
+if length(ftype)>1
+	include_bias = ftype(2);
+	ftype = ftype(1);
+end
+if (nargin < 2)
+	prwarning(3,'C set to one');
+	C = 1; 
+end
+if (nargin < 1) | (isempty(x))
+	w = mapping(mfilename,{C,ftype,include_bias});
+	w = setname(w,'LESS classifier');
+	return
+end
+
+  	
+if ~ismapping(C)   % train the mapping
+
+	% Unpack the dataset.
+	islabtype(x,'crisp');
+	isvaldset(x,1,2); % at least 1 object per class, 2 classes
+	[m,k,c] = getsize(x); 
+
+	if c == 2   % two-class classifier
+
+		% get -1/+1 labels:
+		nlab = getnlab(x);
+		y = 2*nlab-3;
+
+		% train and apply the feature mapping:
+		par = lessfx(ftype,x);
+		f = lessfx(par,x);
+
+		if (include_bias)
+			f = [f ones(m,1)];
+			k = k+1;
+		end
+
+		% In the LP formulation, we define the free parameter vector as:
+		%  [delta; w]
+		% setup the constraints:
+		yf = -repmat(y,1,k).*f;
+		% standard version when we have Ax<b  and  Aeq x = b;
+		A = [-eye(m)    -(+yf)];
+		b = -ones(m,1);
+		%Aeq = [zeros(1,m) ones(1,k)]; beq = 1;
+		%if (include_bias), Aeq(1,end)=0; end
+		Aeq = []; beq = [];
+		% function to optimize:
+		c = [repmat(C,1,m) ones(1,k)];
+		%c = [ones(1,m) repmat(C,1,k)];
+		if (include_bias), c(end) = 0; end
+		% upper and lower bounds:
+		lb = zeros(m+k,1);
+		if (include_bias), lb(end) = -inf; end
+		ub = repmat(inf,m+k,1);
+
+		% optimize
+		if (exist('glpkmex')==3)
+			[out,dummy]=glpkmex(1,c',A,b,repmat('U',m,1),lb,[],repmat('C',m+k,1));
+		else
+			out = linprog(c,A,b,Aeq,beq,lb,ub);
+		end
+		w = out((m+1):end);
+
+		% find out how many features are relevant:
+		if (include_bias)
+			I = find(abs(w(1:(end-1)))>0);
+			nr = length(I);
+		else
+			I = find(abs(w)>0);
+			nr = length(I);
+		end
+
+		% Store the classifier
+		W.extend = include_bias;
+		W.par = par;
+		W.w = w;
+		W.nr = nr;
+		w = mapping(mfilename,'trained',W,getlablist(x),size(x,2),2);
+		w = setname(w,'LESS classifier');
+	
+	else   % multi-class classifier:
+		
+		%error('Multiclass not implemented yet');
+		w = mclassc(x,mapping(mfilename,{C,ftype,include_bias}));
+		v = w.data{1}.data{1}.data.w;
+		for i=2:length(w.data)
+			v = v + w.data{i}.data{1}.data.w;
+		end
+		I = find(abs(v)>0);
+		
+	end	
+else
+	% Evaluate the classifier on new data:
+	W = getdata(C);
+
+	% It is a simple linear classifier:
+	if (W.extend)
+		out = [lessfx(W.par,x) ones(size(x,1),1)]*W.w;
+	else
+		out = lessfx(W.par,x)*W.w;
+	end
+
+	% and put it nicely in a prtools dataset:
+	w = setdat(x,sigm([out -out]),C);
+
+end
+		
+return
Index: /prextra/prextra/lessfx.m
===================================================================
--- /prextra/prextra/lessfx.m	(revision 2)
+++ /prextra/prextra/lessfx.m	(revision 2)
@@ -0,0 +1,104 @@
+function out = lessfx(par,x)
+%LESSFX simple dataset mappings
+%
+%     PAR = LESSFX(TYPE,X)
+%
+% 'Train' or define a mapping of different types.
+% TYPE can be:
+%  1:    (x - m_2).^2 - (x - m_1).^2       normal nearest means
+%  2:  ((x-m_2).^2)/s1 - ((x-m_1).^2)/s2   weighted nearest means
+%  3:    (x-|m|_2).^2 - (x-|m|_1).^2       nearest medians
+%  4: exp(-((x-M1).^2)./S1) - exp(-((x-M2).^2)./S2);
+%  5:        "                 "     with medians instead of %  avergs
+%
+% New data is mapped using:
+%     Y = LESSFX(PAR,X)
+% where X in the input dataset, and PAR is obtained as above...
+%
+% This is used in LESS.M
+
+
+if ~isstruct(par)  % we train the parameters of the function
+	out = [];
+	out.type = par;
+	switch out.type
+	case 0
+		out.bla = [];
+	case 1   % basic version
+		out.u = +meancov(x);
+	case 2   % mean-var
+		[u,g] = meancov(x);
+		out.u = +u;
+		out.g(1,:) = diag(g(:,:,1)) + mean(diag(g(:,:,1)));
+		out.g(2,:) = diag(g(:,:,2)) + mean(diag(g(:,:,2)));
+	case 3   % median
+		out.u(1,:) = med(+seldat(x,1));
+		out.u(2,:) = med(+seldat(x,2));
+	case 4   % mean-var
+		[u,g] = meancov(x);
+		out.u = +u;
+		out.g(1,:) = diag(g(:,:,1)) + mean(diag(g(:,:,1)));
+		out.g(2,:) = diag(g(:,:,2)) + mean(diag(g(:,:,2)));
+	case 5   % median-MSD
+		[u,g] = meancov(x);
+		X1 = seldat(x,1);
+		X2 = seldat(x,2);
+		out.u(1,:) = med(+X1);
+		out.u(2,:) = med(+X2);
+		out.g(1,:) = medstd(+X1,out.u(1,:)) + mean(diag(g(:,:,1)));
+		out.g(2,:) = medstd(+X2,out.u(2,:)) + mean(diag(g(:,:,2)));
+	case 6 % median-MSD
+		X1 = seldat(x,1);
+		X2 = seldat(x,2);
+		out.u(1,:) = med(+X1);
+		out.u(2,:) = med(+X2);
+		out.g(1,:) = medstd(+X1,out.u(1,:));
+		out.g(2,:) = medstd(+X2,out.u(2,:));
+	otherwise
+		error('This function is not defined');
+	end
+else                 % we evaluate the function on new data:
+	[m,k] = size(x);
+	switch par.type
+	case 0
+		out = +x;
+	case 1
+		M1 = repmat(par.u(1,:),m,1);
+		M2 = repmat(par.u(2,:),m,1);
+		out = -(2*(+x).*(M2-M1) + M1.*M1 - M2.*M2);
+	case 2
+		M1 = repmat(par.u(1,:),m,1);
+		M2 = repmat(par.u(2,:),m,1);
+		S1 = repmat(par.g(1,:),m,1);
+		S2 = repmat(par.g(2,:),m,1);
+		out = ((x-M2).^2)./S2 - ((x-M1).^2)./S1;
+	case 3
+		M1 = repmat(par.u(1,:),m,1);
+		M2 = repmat(par.u(2,:),m,1);
+		out = -(2*(+x).*(M2-M1) + M1.*M1 - M2.*M2);
+	case 4
+		M1 = repmat(par.u(1,:),m,1);
+		M2 = repmat(par.u(2,:),m,1);
+		S1 = repmat(par.g(1,:),m,1);
+		S2 = repmat(par.g(2,:),m,1);
+		out = exp(-((x-M1).^2)./S1) - exp(-((x-M2).^2)./S2);
+	case 5
+		M1 = repmat(par.u(1,:),m,1);
+		M2 = repmat(par.u(2,:),m,1);
+		S1 = repmat(par.g(1,:),m,1);
+		S2 = repmat(par.g(2,:),m,1);
+		out = ((x-M2).^2)./S2 - ((x-M1).^2)./S1;
+	case 6
+		M1 = repmat(par.u(1,:),m,1);
+		M2 = repmat(par.u(2,:),m,1);
+		S1 = repmat(par.g(1,:),m,1);
+		S2 = repmat(par.g(2,:),m,1);
+		%out = ((x-M2).^2)./S2 - ((x-M1).^2)./S1;
+		out = (sigm(10-((x-M1).^2)./(1*S1)) - sigm(10-((x-M2).^2)./(1*S2)));
+	otherwise
+		error('This function is not defined');
+	end
+
+end
+
+return
Index: /prextra/prextra/linprogc.m
===================================================================
--- /prextra/prextra/linprogc.m	(revision 2)
+++ /prextra/prextra/linprogc.m	(revision 2)
@@ -0,0 +1,82 @@
+%LINPROGC Sparce linear programming classifier
+%
+%   W = LINPROGC(A,N,TYPE,PARAM)
+%   W = A*LINPROGC([],N,TYPE,PARAM)
+%
+% INPUT
+%   A     Dataset
+%   N     Degree of classifier, 1,2,3,  ...
+%   TYPE  Linear programming option, see below
+%   PARAM Additional parameter
+%
+% OUTPUT
+%   W     Classifier
+%
+% The parameter TYPE stands for the possible LP-machines:
+% 'SIMPLE'   -  the most simple formulation; no sparse solution, hence 
+%               all features are used.
+% 'STANDARD' -  minimization of the training misclassification errors; no 
+%               sparse solution, hence all features are used.
+% 'C-SPARSE' -  a sparse solution (feature selection); a formulation similar 
+%               to the LP_1 SVC, adjusted for the purpose of feature selection. 
+%               PARAM is a tradeoff parameter, as in the traditional SVC.
+%               Usually, PARAM = 1 or more. Default PARAM = 1.
+% 'NU-SPARSE' - a sparse solution (feature selection); a formulation similar 
+%               to the LP_1 SVC, based on the paper of Graepel, Herbrich, Smola etc
+%               'Classification on proximity data with LP-machines' adjusted
+%               for the purpose of feature selection. 
+%               PARAM is a tradeoff parameter, usually PARAM = 0.05 or 0.1. 
+%               It is an upper bound on the misclassfied training objects.      
+%               So, for well separable problems, PARAM = 0.01 or PARAM = 0.05.
+%               Default PARAM is the leave-one-out 1-NN error of A.
+%
+% DESCRIPTION
+% A linear programming solution, similar to the LP_1 SVC is found for the 
+% polynomial classifier (degree N, default N=1).
+% Sparse solutions based on TYPE = C-SPARSE and TYPE = NU-SPARSE reduce the
+% number of polynomial terms, thereby performing an automatic feature
+% selection.
+%
+% This routine is effectively a wrapper around FSELLPC
+%
+% SEE ALSO
+% MAPPINGS, DATASETS, POLYC, FSELLPC
+
+% Elzbieta Pekalska, Robert P.W. Duin, e.pekalska@ewi.tudelft.nl
+% Faculty of Electrical Engineering, Mathematics and Computer Science,
+% Delft University of Technology, The Netherlands.
+
+function w = linprogc(a,n,type,param)
+if nargin < 4, param =[]; end
+if nargin < 3, type = 'standard'; end
+if nargin < 2, n =  1; end
+name = [type '-LP'];
+if nargin < 1 | isempty(a)
+	w = mapping(mfilename,{n,type});
+	w = setname(w,name);
+	return
+end
+vs = scalem(a,'variance');
+b = a*vs;
+switch type
+	case {'simple', 'Simple', 'SIMPLE'}
+		type = 'Simple';
+		param = 1;
+	case {'standard', 'Standard', 'STANDARD'}
+		type = 'Standard';
+		param = 1;
+	case {'c-sparse', 'C-Sparse', 'c-Sparse'}
+		type = 'c-Sparse';
+		if isempty(param), param = 1; end
+	case {'nu-sparse','nu-Sparse','\nu-Sparse'}
+		type = 'nu-Sparse';
+		if isempty(param), param = testk(b,1); end
+	otherwise
+		error('Wrong Type')
+end
+vp = fsellpc([],1,type,param);
+w = polyc(b,vp,n,1); 
+w = vs*w;
+name = [type '-LP'];
+w = setname(w,name);
+return
Index: /prextra/prextra/lpsvc.m
===================================================================
--- /prextra/prextra/lpsvc.m	(revision 2)
+++ /prextra/prextra/lpsvc.m	(revision 2)
@@ -0,0 +1,77 @@
+%LPSVC Sparse linear programming classifier according to Mangasarian
+%
+%   [W,E] = LPSVC(A,N,PoARAM)
+%   W     = A*LPSVC([],N,PARAM)
+%
+% INPUT
+%   A     Dataset
+%   N     N-fold cross-validationocedure
+%         N = 1: no cross-validation
+%         N = size A: leave-one-out pr
+%   PARAM Linear programming option
+%         -1 Easy
+%         0  Hard (default)
+%         Otherwise: nu
+%
+% OUTPUT
+%   W     Classifier
+%   E     Internal error estimate
+%
+% DESCRIPTION
+% This routine calls LPSVM by Fung and Mangasarian. If you use it, please
+% refer to the below paper.
+%
+% REFERENCE
+% G.M. Fung and O.L. Mangasarian, A Feature Selection Newton Method for
+% Support Vector Machine Classification, Computational Optimization and
+% Aplications, vol. 28, 2004, 185-202.
+%
+% SEE ALSO
+% MAPPINGS, DATASETS, LINPROGC
+
+% Elzbieta Pekalska, Robert P.W. Duin, e.pekalska@ewi.tudelft.nl
+% Faculty of Electrical Engineering, Mathematics and Computer Science,
+% Delft University of Technology, The Netherlands.
+
+function [w,J] = lpsvm(a,n,param)
+if nargin < 3, param =0; end
+if nargin < 2, n =  1; end
+name = 'LPSVC';
+if nargin < 1 | isempty(a)
+	w = mapping(mfilename,{n,param});
+	w = setname(w,name);
+	return
+end
+islabtype(a,'crisp');
+isvaldset(a,1,2); % at least 1 object per class, 2 classes
+[m,k,c] = getsize(a);
+nlab = getnlab(a);	
+	% The SVC is basically a 2-class classifier. More classes are
+	% handled by mclassc.
+if c == 2   % two-class classifier
+	y = 3 - 2*nlab;
+	[v,v0,ea,et] = lpsvm(+a,y,n,param);
+	J = find(v~=0);
+	if isempty(J)
+		prwarning(1,'No features found. Fisher classifier is trained.')
+		wf = featsel(k,1:k);
+		w = wf*fisherc(a);
+	else
+		wf = featsel(k,J);
+		flab = getfeatlab(a);
+		w = wf*affine(v(J),-v0,flab(J,:),getlablist(a),length(J),c);
+		w = cnormc(w,a);
+	end
+	
+else
+	
+	[w,J] = mclassc(a,mapping(mfilename,{n,param}));
+		
+end
+
+
+
+
+
+
+
Index: /prextra/prextra/lpsvm.m
===================================================================
--- /prextra/prextra/lpsvm.m	(revision 2)
+++ /prextra/prextra/lpsvm.m	(revision 2)
@@ -0,0 +1,289 @@
+function [w,gamma,trainCorr, testCorr, cpu_time, nu]=lpsvm(A,d,k,nu,output,delta)
+% version 1.3
+% last revision: 07/07/03
+%===========================================================================
+% Usage: [w, gamma,trainCorr,testCorr,time,nu]=lpsvm(A,d,k,nu,output,delta);
+%
+% A and d are both required, everything else has a default
+% An example: [w gamma train test time nu] = lpsvm(A,d,10);
+%
+% Input parameters:
+%    A: Data points
+%    d: 1's or -1's
+%    k: way to divide the data set into test and training set
+%       if k = 0: simply run the algorithm without any correctness
+%         calculation, this is the default
+%       if k = 1: run the algorithm and calculate correctness on
+%         the whole data set
+%       if k = any value less than the # of rows in the data set:
+%         divide up the data set into test and training
+%         using k-fold method
+%       if k = # of rows in the data set: use the 'leave 1' method
+%
+%    output: 0 - no output, 1 - produce output, default is 0
+%    nu:             weighted parameter
+%                    -1 - easy estimation
+%                    0  - hard estimation
+%                    any other value - used as nu by the algorithm
+%                    default - 0
+%    delta:  default is 10^-3
+%===================================================================
+% Output parameters:
+%
+%       w:              the normal vector of the classifier
+%       gamma:          the threshold
+%       trainCorr:      training set correctness
+%       testCorr:       test set correctness
+%       cpu_time:       time elapsed
+%       nu:             estimated value (or specified value) of nu
+%==========================================================================
+
+if nargin<6
+delta=1e-3;
+end
+
+if nargin<5
+output=0;
+end
+
+if ((nargin<4)|(nu==0))
+     nu = EstNuLong(A,d);  % default is hard estimation
+elseif nu==-1  % easy estimation
+nu = EstNuShort(A,d);
+end
+
+if nargin<3
+k=0;
+end
+
+r=randperm(size(d,1));d=d(r,:);A=A(r,:);    % random permutation
+
+tic;
+
+trainCorr=0;
+testCorr=0;
+
+if k==0
+[w, gamma,iter] = core(A,d,nu,delta);
+cpu_time=toc;
+  if output==1
+fprintf(1,'\nNumber of Iterations: %d',iter);
+fprintf(1,'\nElapse time: %10.2f\n\n',cpu_time);
+  end
+  return
+end
+
+%if k==1 only training set correctness is calculated
+if k==1
+tic;
+[w, gamma,iter] = core(A,d,nu,delta);
+trainCorr = correctness(A,d,w,gamma);
+cpu_time = toc;
+  if output == 1
+fprintf(1,'\nTraining set correctness: %3.2f%%',trainCorr);
+fprintf(1,'\nNumber of Iterations: %d',iter);
+fprintf(1,'\nElapse time: %10.2f\n\n',cpu_time);
+  end
+  return
+end
+
+[sm sn]=size(A);
+accuIter = 0;
+
+cpu_time = 0;
+indx = [0:k];
+indx = floor(sm*indx/k);    %last row numbers for all 'segments'
+% split trainining set from test set
+for i = 1:k
+Ctest = []; dtest = [];Ctrain = []; dtrain = [];
+
+Ctest = A((indx(i)+1:indx(i+1)),:);
+dtest = d(indx(i)+1:indx(i+1));
+
+Ctrain = A(1:indx(i),:);
+Ctrain = [Ctrain;A(indx(i+1)+1:sm,:)];
+dtrain = [d(1:indx(i));d(indx(i+1)+1:sm,:)];
+tic;
+[w, gamma,iter] = core(Ctrain,dtrain,nu,delta);
+thisToc = toc;
+
+tmpTrainCorr = correctness(Ctrain,dtrain,w,gamma);
+tmpTestCorr = correctness(Ctest,dtest,w,gamma);
+
+ if output==1
+fprintf(1,'________________________________________________\n');
+fprintf(1,'Fold %d\n',i);
+fprintf(1,'Training set correctness: %3.2f%%\n',tmpTrainCorr);
+fprintf(1,'Testing set correctness: %3.2f%%\n',tmpTestCorr);
+fprintf(1,'Number of iterations: %d\n',iter);
+fprintf(1,'Elapse time: %10.2f\n',thisToc);
+end
+
+trainCorr = trainCorr + tmpTrainCorr;
+testCorr = testCorr + tmpTestCorr;
+accuIter = accuIter + iter; % accumulative iterations
+cpu_time = cpu_time + thisToc;
+
+end % end of for (looping through test sets)
+
+     trainCorr = trainCorr/k;
+     testCorr = testCorr/k;
+     cpu_time=cpu_time/k;
+
+if output == 1
+     fprintf(1,'==============================================');
+fprintf(1,'\nTraining set correctness: %3.2f%%',trainCorr);
+fprintf(1,'\nTesting set correctness: %3.2f%%',testCorr);
+fprintf(1,'\nAverage number of iterations: %d',accuIter/k);
+fprintf(1,'\nAverage cpu_time: %10.2f\n',cpu_time);
+end
+
+return;  % lpsvm function return
+
+%%%%%%%%%%% core calculation function %%%%%%%%%%%%%%%%%%%%%
+function [w, gamma, iter] = core(A,d,nu,delta);
+
+[m,n]=size(A);
+
+if m>=n
+[w,gamma,iter]=lpsvm_with_smw(A,d,nu,delta);
+else
+[w,gamma,iter]=lpsvm_without_smw(A,d,nu,delta);
+end
+
+return
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%   lpsvm when m>=n                                   %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [w,gamma,iter]=lpsvm_with_smw(A,d,nu,delta)
+%with SMW without armijo
+%parameters
+epsi=10^(-3);alpha=10^(2);tol=10^(-5);maxiter=100;
+[m,n]=size(A);
+en=ones(n,1);
+em=ones(m,1);
+
+% initial u
+u=ones(m,1);iter=0;
+epsi=epsi*em;nu=nu*em;
+diff=1;
+DA=spdiags(d,0,m,m)*A;
+while (diff>tol) & (iter<maxiter)
+    uold=u;
+    iter=iter+1;
+    du=d.*u;Adu=A'*du;
+    pp=max(Adu-en,0);np=max(-Adu-en,0);
+    dd=sum(du)*d;unu=max(u-nu,0);uu=max(-u,0);
+    %Gradient
+    g=-epsi+(d.*(A*pp))-(d.*(A*np))+dd+unu-alpha*uu;
+    %Hessian
+    E=spdiags(sqrt(sign(np)+sign(pp)),0,n,n);
+    H=[DA*E d];
+    f=1./(delta+sign(unu)+alpha*sign(uu));
+    F=spdiags(f,0,m,m);gg=f.*g;HT=H';
+    di=(eye(n+1)+HT*(F*H))\(HT*gg);
+    di=H*di;di=f.*di;di=-gg+di;u=u+di;
+    diff=norm(g);
+end
+
+w=1/epsi(1)*(pp-np);
+gamma=-(1/epsi(1))*sum(du);
+iter;
+return
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%   lpsvm when m<n                                    %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [w,gamma,iter]=lpsvm_without_smw(A,d,nu,delta)
+%without sherman and without armijo
+%parameters
+epsi=10^(-1);alpha=10^3;tol=10^(-3);maxiter=50;
+[m,n]=size(A);
+en=ones(n,1);
+em=ones(m,1);
+%initial u
+u=ones(m,1);iter=0;
+epsi=epsi*em;nu=nu*em;
+diff=1;
+DA=spdiags(d,0,m,m)*A;
+while (diff>tol) & (iter<maxiter)
+    uold=u;
+    iter=iter+1;
+    du=d.*u;Adu=A'*du;
+    pp=max(Adu-en,0);np=max(-Adu-en,0);
+    dd=sum(du)*d;unu=max(u-nu,0);uu=max(-u,0);
+    %Gradient
+    g=-epsi+(d.*(A*pp))-(d.*(A*np))+dd+unu-alpha*uu;
+    %Hessian
+    E=spdiags(sqrt(sign(np)+sign(pp)),0,n,n);
+    H=[DA*E d];
+    F=spdiags(delta+sign(unu)+alpha*sign(uu),0,m,m);
+    di=-((H*H'+F)\g);
+    u=u+di;
+    diff=norm(g);
+end
+du=d.*u;Adu=A'*du;
+pp=max(Adu-en,0);np=max(-Adu-en,0);
+w=1/epsi(1)*(pp-np);gamma=-(1/epsi(1))*sum(du);
+return
+
+%%%%%%%%%%%%%%%% correctness calculation %%%%%%%%%%%%%%%%
+
+function corr = correctness(AA,dd,w,gamma)
+
+p=sign(AA*w-gamma);
+corr=length(find(p==dd))/size(AA,1)*100;
+return
+
+%%%%%%%%%%%%%%EstNuLong%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% hard way to estimate nu if not specified by the user
+function value = EstNuLong(C,d)
+
+[m,n]=size(C);e=ones(m,1);
+H=[C -e];
+if m<201
+H2=H;d2=d;
+else
+r=rand(m,1);
+[s1,s2]=sort(r);
+H2=H(s2(1:200),:);
+d2=d(s2(1:200));
+end
+
+lamda=1;
+[vu,u]=eig(H2*H2');u=diag(u);p=length(u);
+yt=d2'*vu;
+lamdaO=lamda+1;
+
+cnt=0;
+while (abs(lamdaO-lamda)>10e-4) &(cnt<100)
+     cnt=cnt+1;
+     nu1=0;pr=0;ee=0;waw=0;
+     lamdaO=lamda;
+     for i=1:p
+     nu1= nu1 + lamda/(u(i)+lamda);
+pr= pr + u(i)/(u(i)+lamda)^2;
+ee= ee + u(i)*yt(i)^2/(u(i)+lamda)^3;
+waw= waw + lamda^2*yt(i)^2/(u(i)+lamda)^2;
+   end
+lamda=nu1*ee/(pr*waw);
+end
+
+value = lamda;
+if cnt==100
+    value=1;
+end
+
+return
+
+%%%%%%%%%%%%%%%%%EstNuShort%%%%%%%%%%%%%%%%%%%%%%%
+
+% easy way to estimate nu if not specified by the user
+function value = EstNuShort(C,d)
+
+value = 1/(sum(sum(C.^2))/size(C,2));
+return
Index: /prextra/prextra/mexlasso.c
===================================================================
--- /prextra/prextra/mexlasso.c	(revision 2)
+++ /prextra/prextra/mexlasso.c	(revision 2)
@@ -0,0 +1,1003 @@
+/* Copyright (C) 1998
+   Berwin A Turlach <bturlach@stats.adelaide.edu.au> */
+
+/* This library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Library General Public License
+   as published by the Free Software Foundation; either version 2 of
+   the License, or (at your option) any later version. */
+
+/* This library is distributed in the hope that it will be useful, but
+   WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Library General Public License for more details. */
+
+/* You should have received a copy of the GNU Library General Public
+   License along with this library; if not, write to the Free Software
+   Foundation, Inc., 59 Temple Place, Suite 330, Boston,
+   MA 02111-1307, USA. */
+
+#include "lasso.h"
+//#include "fortify.h"
+
+#define PRECISION FLT_EPSILON
+
+static void lasso_alloc(long n, long m);
+static void lasso_free(void);
+static void qr_init(int n);
+static void qr_incr(void);
+static void qr_free(void);
+static void qr_del(int l, int aug);
+static void qr_add(double *x, int swap);
+#if defined (S_Plus)
+static void errmsg(char* string);
+#else
+static void errmsg(char* where, char* string);
+#endif
+
+  static int QR_CHUNK = 10;
+
+  static double *xtr=NULL, *btmp=NULL, *qtr=NULL,
+                *rinvt_theta=NULL, *step=NULL, ytyd2=0.0;
+  static int *theta=NULL, *nz_x=NULL, num_nz_x=0;
+  static double *qmat=NULL, *rmat=NULL;
+  static int qr_max_size=0, r_ncol=0, q_nrow=0, q_use_row=0;
+  static char *no_dyn_mem_message="Cannot allocate dynamic memory";
+
+void lasso(double *x,           /* input: data */
+           long *pn,            /* input: rows */
+           long *pm,            /* input: columns */
+           double *pt,          /* input: threshold */
+           double *beta,        /* input/ouput: intial/final beta */
+           double *y,           /* input: labels */
+
+           double *yhat1,       /* output */
+           double *r,           /* output */
+           double *lagrangian,  /* output */
+           long *psuc,          /* output: succeed: 0/-1 */
+           long *pverb,         /* input: verbose: yes/no */
+           long *pas_sub)       /* input: no!? */
+{
+
+  double t = *pt, prec;
+  long n = *pn, m = *pm, verb = *pverb, as_sub = *pas_sub;
+  int not_solved;
+  double *x_elem=NULL, tmp, max_val, b_1norm;
+  int i, j, max_ind;
+  double p_obj, d_obj;
+  int num_iter=0, max_num;
+  double b_1norm_old, mu;
+  double rho_up, rho_low, rho;
+  int to_del;
+  double *q_elem, wtc, wtw;
+  double *r_elem;
+  int k;
+  int add;
+  if( !as_sub)
+    lasso_alloc(n,m);
+
+  prec = sqrt(PRECISION);
+
+  if( as_sub ){
+
+  b_1norm = 0.0;
+  for(j=0;j<num_nz_x;j++)
+    b_1norm += fabs(beta[nz_x[j]]);
+  }else{
+
+  b_1norm = 0.0;
+  num_nz_x = 0;
+  for(j=0; j<m; j++){
+    if(fabs(beta[j]) > prec){
+      b_1norm += fabs(beta[j]);
+      nz_x[num_nz_x] = j;
+      num_nz_x++;
+    }else
+      beta[j] = 0.0;
+  }
+  }
+
+  if( b_1norm > t){
+    if(verb){
+      mexPrintf("******************************\n");
+      mexPrintf("Rescaling beta from L1-norm %f to %f\n",b_1norm,t);
+      mexEvalString("drawnow;");
+
+    }
+    for(j=0; j<num_nz_x; j++)
+      beta[nz_x[j]] = beta[nz_x[j]] * t/b_1norm;
+    b_1norm = t;
+  }
+
+  for(i=0; i< n; i++)
+    yhat1[i] = 0.0;
+  for(j=0; j < num_nz_x; j++){
+    /* x_elem points to the first element in the column of X to which
+       the j-th entry in nz_x points  */
+    x_elem = x + nz_x[j]*n;
+    tmp = beta[nz_x[j]];
+    for(i=0; i < n; i++){
+      yhat1[i] += *x_elem*tmp;
+      x_elem++;           /* now we point to the next element in X */
+    }
+  }
+
+  /* calculate the residual vector */
+  for(i=0; i< n; i++)
+    r[i] = y[i]-yhat1[i];
+
+  /* multiply X^T with the residual vector */
+  x_elem = x;
+  for(j=0; j < m; j++){
+    tmp = 0.0;
+    for(i=0; i<n; i++){
+      tmp += *x_elem*r[i];
+      x_elem++;
+    }
+    xtr[j] = tmp;
+  }
+
+  max_val = fabs(xtr[0]);
+  max_ind = 0;
+  for(j=1; j<m; j++)
+    if( fabs(xtr[j]) > max_val ){
+      max_val = fabs(xtr[j]);
+      max_ind = j;
+    }
+
+  if( !as_sub ){
+
+  qr_add(y,TRUE);
+  ytyd2 = *rmat * *rmat/2.0;
+  for(j=0;j<num_nz_x;j++){
+    qr_add(x+nz_x[j]*n, TRUE);
+    if(fabs(beta[nz_x[j]])<prec)
+      theta[j] = xtr[nz_x[j]] < 0 ? -1 : 1;
+    else
+      theta[j] = beta[nz_x[j]] < 0 ? -1 : 1;
+
+  }
+  }
+  if( num_nz_x==0 ){
+    nz_x[0] = max_ind;
+    num_nz_x = 1;
+    if(verb){
+      mexPrintf("******************************\n");
+      mexPrintf("  -->\tAdding variable: %d\n",max_ind+1);
+      mexEvalString("drawnow;");
+    }
+    qr_add(x+max_ind*n, TRUE);
+    theta[0] = xtr[max_ind] < 0 ? -1 : 1;
+  }
+  *psuc=0;
+  if(verb){
+
+
+  /* Find out how many times [[max_val]] is attained */
+  tmp = (1.0-prec)*max_val;
+  if( tmp < prec ) tmp = 0.0;
+  max_num = 1;
+  for(j=0; j<m; j++){
+    if( tmp <= fabs(xtr[j]) && j!=max_ind){
+      /*we found another element equal to the (current)
+        maximal absolute value */
+      max_num++;
+    }
+  }
+
+  /* for the value of the dual objective function we need
+     to calculate the L2 norm of the vector of fitted values */
+  p_obj = 0.0;
+  d_obj = 0.0;
+  for(i=0;i<n;i++){
+    p_obj += r[i]*r[i];
+    d_obj += yhat1[i]*yhat1[i];
+  }
+  p_obj /= 2.0;
+  d_obj = ytyd2 - d_obj/2.0 - t*max_val ;
+  mexPrintf("******************************\n");
+  mexPrintf("\nIteration number: %d\n", num_iter);
+  mexPrintf("Value of primal object function      : %f\n", p_obj);
+  mexPrintf("Value of dual   object function      : %f\n", d_obj);
+  mexPrintf("L1 norm of current beta              : %f", b_1norm);
+  mexPrintf(" <= %f\n", t);
+  mexPrintf("Maximal absolute value in t(X)%%*%%r   : %e", max_val);
+  mexPrintf(" attained %d time(s)\n", max_num);
+  mexPrintf("Number of parameters allowed to vary : %d\n", num_nz_x);
+  mexEvalString("drawnow;");
+  num_iter++;
+  }
+  while(1){
+    do{
+
+  q_elem = qmat;
+
+  for(j=0;j<num_nz_x;j++){
+    tmp = 0.0;
+    for(i=0;i<n;i++,q_elem++) {
+      tmp += *q_elem*r[i];
+    }
+    qtr[j] = tmp;
+  }
+
+  if( b_1norm < (1.0-prec)*t )
+    mu = 0.0;
+  else{
+
+  /* z=R^{-T}\theta can be calculated by solving R^Tz=\theta */
+  r_elem = rmat;
+  for(j=0;j<num_nz_x;j++,r_elem++){
+    tmp = theta[j];
+    for(k=0;k<j;k++,r_elem++)
+      tmp -= *r_elem * rinvt_theta[k];
+    rinvt_theta[j] = tmp /  *r_elem;
+  }
+
+  wtc = 0.0;
+  wtw = 0.0;
+  for(j=0;j<num_nz_x;j++){
+    wtc += rinvt_theta[j]*qtr[j];
+    wtw += rinvt_theta[j]*rinvt_theta[j];
+  }
+  mu = (wtc-(t-b_1norm))/wtw;
+  }
+
+  if(mu<=0.0)
+    for(j=0;j<num_nz_x;j++)
+      step[j] = qtr[j];
+  else
+    for(j=0;j<num_nz_x;j++)
+      step[j] = qtr[j] - mu*rinvt_theta[j];
+
+  /* h=R^{-1}z can be calculated by solving Rh=z */
+  for(j=num_nz_x-1;j>=0;j--){
+    tmp = step[j];
+    for(k=num_nz_x-1;k>j;k--)
+      tmp -= RMAT(j,k) * step[k];
+    step[j] = tmp /  RMAT(j,j);
+  }
+  for(j=0;j<num_nz_x;j++){
+    btmp[j] = beta[nz_x[j]];
+    beta[nz_x[j]] += step[j];
+  }
+  b_1norm_old=b_1norm;
+
+  b_1norm = 0.0;
+  for(j=0;j<num_nz_x;j++) {
+    b_1norm += fabs(beta[nz_x[j]]);
+  }
+
+  not_solved=FALSE;
+  if( b_1norm > (1+prec)*t){
+    not_solved=TRUE;
+    if(b_1norm_old < (1.0-prec)*t){
+
+    if(verb) {
+      mexPrintf("  -->\tStepping onto the border of the L1 ball.\n");
+      mexEvalString("drawnow;");
+    }
+    rho_up  = rho = 1.0;
+    rho_low = 0.0;
+    while( fabs(t-b_1norm) > prec*t ){
+      if( b_1norm > t){
+        rho_up = rho;
+        rho = (rho+rho_low)/2.0;
+      }
+      if( b_1norm < t){
+        rho_low = rho;
+        rho = (rho+rho_up)/2.0;
+      }
+      if(rho < prec) break;
+      for(j=0; j<num_nz_x; j++) {
+        beta[nz_x[j]] = btmp[j]+rho*step[j];
+      }
+
+  b_1norm = 0.0;
+  for(j=0;j<num_nz_x;j++)
+    b_1norm += fabs(beta[nz_x[j]]);
+    }
+    for(j=0;j<num_nz_x;j++){
+      if(fabs(beta[nz_x[j]])<prec)
+        theta[j] = btmp[j] > 0 ? -1 : 1;
+      else
+        theta[j] = beta[nz_x[j]] < 0 ? -1 : 1;
+    }
+    }else{
+
+
+  rho = 1.0;
+  to_del = -1;
+  for(j=0; j<num_nz_x; j++){
+    if(fabs(step[j]) > prec){
+      tmp = -btmp[j]/(step[j]);
+      if( 0.0 < tmp && tmp < rho ){
+        rho = tmp;
+        to_del = j;
+      }
+    }
+  }
+  if(to_del < 0 ){
+    *psuc= -1;
+    goto EXIT_HERE;
+  }
+  for(j=0; j<num_nz_x; j++) {
+    beta[nz_x[j]] = btmp[j]+rho*step[j];
+  }
+  if(verb) {
+    mexPrintf("  -->\tRemoving variable: %d",nz_x[to_del]+1);
+    mexEvalString("drawnow;");
+  }
+  beta[nz_x[to_del]] = 0.0;
+  qr_del(to_del,TRUE);
+  for(j=to_del+1; j< num_nz_x; j++){
+    nz_x[j-1] = nz_x[j];
+    theta[j-1] = theta[j];
+  }
+  num_nz_x--;
+  
+  b_1norm = 0.0;
+  for(j=0;j<num_nz_x;j++)
+    b_1norm += fabs(beta[nz_x[j]]);
+  if(verb)
+    if(b_1norm < (1-prec)*t)
+      mexPrintf(", and stepping into the interior of the L1 ball\n");
+    else
+      mexPrintf("\n");
+    }
+    mexEvalString("drawnow;");
+  }
+
+  for(i=0; i< n; i++)
+    yhat1[i] = 0.0;
+  for(j=0; j < num_nz_x; j++){
+    /* x_elem points to the first element in the column of X to which
+       the j-th entry in nz_x points  */
+    x_elem = x + nz_x[j]*n;
+    tmp = beta[nz_x[j]];
+    for(i=0; i < n; i++){
+      yhat1[i] += *x_elem*tmp;
+      x_elem++;           /* now we point to the next element in X */
+    }
+  }
+
+  /* calculate the residual vector */
+  for(i=0; i< n; i++) {
+    r[i] = y[i]-yhat1[i];
+  }
+    }while(not_solved);
+    
+  
+  for(i=0; i< n; i++)
+    yhat1[i] = 0.0;
+  for(j=0; j < num_nz_x; j++){
+    /* x_elem points to the first element in the column of X to which
+       the j-th entry in nz_x points  */
+    x_elem = x + nz_x[j]*n;
+    tmp = beta[nz_x[j]];
+    for(i=0; i < n; i++){
+      yhat1[i] += *x_elem*tmp;
+      x_elem++;           /* now we point to the next element in X */
+    }
+  }
+
+  /* calculate the residual vector */
+  for(i=0; i< n; i++)
+    r[i] = y[i]-yhat1[i];
+
+
+  /* multiply X^T with the residual vector */
+  x_elem = x;
+  for(j=0; j < m; j++){
+    tmp = 0.0;
+    for(i=0; i<n; i++){
+      tmp += *x_elem*r[i];
+      x_elem++;
+    }
+    xtr[j] = tmp;
+  }
+  max_val = fabs(xtr[0]);
+  max_ind = 0;
+  for(j=1; j<m; j++)
+    if( fabs(xtr[j]) > max_val ){
+      max_val = fabs(xtr[j]);
+      max_ind = j;
+    }
+
+    if(verb){
+
+
+  /* Find out how many times [[max_val]] is attained */
+  tmp = (1.0-prec)*max_val;
+  if( tmp < prec ) tmp = 0.0;
+  max_num = 1;
+  for(j=0; j<m; j++){
+    if( tmp <= fabs(xtr[j]) && j!=max_ind){
+      /*we found another element equal to the (current)
+        maximal absolute value */
+      max_num++;
+    }
+  }
+
+  /* for the value of the dual objective function we need
+     to calculate the L2 norm of the vector of fitted values */
+  p_obj = 0.0;
+  d_obj = 0.0;
+  for(i=0;i<n;i++){
+    p_obj += r[i]*r[i];
+    d_obj += yhat1[i]*yhat1[i];
+  }
+  p_obj /= 2.0;
+  d_obj = ytyd2 - d_obj/2.0 - t*max_val ;
+  mexPrintf("******************************\n");
+  mexPrintf("\nIteration number: %d\n", num_iter);
+  mexPrintf("Value of primal object function      : %f\n", p_obj);
+  mexPrintf("Value of dual   object function      : %f\n", d_obj);
+  mexPrintf("L1 norm of current beta              : %f", b_1norm);
+  mexPrintf(" <= %f\n", t);
+  mexPrintf("Maximal absolute value in t(X)%%*%%r   : %e", max_val);
+  mexPrintf(" attained %d time(s)\n", max_num);
+  mexPrintf("Number of parameters allowed to vary : %d\n", num_nz_x);
+  mexEvalString("drawnow;");
+  num_iter++;
+    }
+    
+  add = TRUE;
+  for(i=0; i<num_nz_x; i++)
+    if(nz_x[i]==max_ind){
+      add = FALSE;
+      break;
+    }
+  if(add){
+    qr_add(x+max_ind*n,TRUE);
+    nz_x[num_nz_x] = max_ind;
+    theta[num_nz_x] = xtr[max_ind]<0 ? -1 : 1;
+    num_nz_x++;
+  }else{
+    *lagrangian = max_val;
+    break;
+  }
+
+  b_1norm = 0.0;
+  for(j=0;j<num_nz_x;j++)
+    b_1norm += fabs(beta[nz_x[j]]);
+  if(verb) {
+    mexPrintf("  -->\tAdding variable: %d\n",max_ind+1);
+    mexEvalString("drawnow;");
+  }
+  }
+  EXIT_HERE:
+  if( !as_sub)
+    lasso_free();
+}
+void mult_lasso(double *x, long *pn, long *pm, double *pt, long *pl,
+                double *beta, double *y, double *yhat1, double *r,
+                double *lagrangian, long *psuc, long *pverb)
+{
+
+  double prec;
+  long n = *pn, m = *pm, l = *pl, verb = *pverb, as_sub = TRUE, i, j;
+  lasso_alloc(n,m);
+
+qr_add(y,TRUE);
+ytyd2 = *rmat * *rmat/2.0;
+prec = sqrt(PRECISION);
+num_nz_x = 0;
+for(j=0; j<m; j++){
+  if(fabs(beta[j]) > prec){
+    qr_add(x+j*n, TRUE);
+    nz_x[num_nz_x] = j;
+    num_nz_x++;
+  }else
+    beta[j] = 0.0;
+}
+*psuc = 0;
+  for(i=0; i<l; i++){
+
+if(verb){
+  mexPrintf("\n\n++++++++++++++++++++++++++++++\n");
+  mexPrintf("Solving problem number %ld with bound %f\n", i+1, pt[i]);
+  mexPrintf("++++++++++++++++++++++++++++++\n");
+  mexEvalString("drawnow;");
+}
+if(i>0) Memcpy(beta,beta-m,m);
+lasso(x, pn, pm, pt+i, beta, y, yhat1, r, lagrangian, psuc, pverb, &as_sub);
+if( *psuc < 0 ){
+  goto EXIT_HERE;
+}
+
+beta += m;
+yhat1 += n;
+r += n;
+lagrangian++;
+  }
+  EXIT_HERE:
+  lasso_free();
+}
+static void lasso_alloc(long n, long m){
+
+#if defined(S_Plus)
+  if( nz_x != NULL || theta != NULL || xtr != NULL || btmp != NULL ||
+      qtr != NULL || rinvt_theta != NULL || step != NULL ||
+      num_nz_x != 0 || ytyd2 != 0.0){
+    MESSAGE "Possible memory corruption or memory leak.\n  We"
+      "advise to restart your S+ session"  WARNING(NULL_ENTRY);
+     lasso_free();
+  }
+#endif
+
+  QR_CHUNK = m;  // Added: CJV
+
+  nz_x = Calloc(m,int);
+  if( nz_x == NULL )
+    ERRMSG("lasso_alloc", no_dyn_mem_message);
+  theta = Calloc(m,int);
+  if( theta==NULL )
+    ERRMSG("lasso_alloc", no_dyn_mem_message);
+  xtr = Calloc(m,double);
+  if( xtr==NULL )
+    ERRMSG("lasso_alloc", no_dyn_mem_message);
+  btmp = Calloc(m,double);
+  if( btmp==NULL )
+    ERRMSG("lasso_alloc", no_dyn_mem_message);
+  qtr = Calloc(m,double);
+  if( qtr==NULL )
+    ERRMSG("lasso_alloc", no_dyn_mem_message);
+  rinvt_theta = Calloc(m,double);
+  if( rinvt_theta==NULL )
+    ERRMSG("lasso_alloc", no_dyn_mem_message);
+  step = Calloc(m,double);
+  if( step==NULL )
+    ERRMSG("lasso_alloc", no_dyn_mem_message);
+  qr_init(n);
+}
+static void lasso_free(void){
+  num_nz_x=0;
+  ytyd2 = 0.0;
+  Free(nz_x);
+  Free(theta);
+  Free(xtr);
+  Free(btmp);
+  Free(qtr);
+  Free(rinvt_theta);
+  Free(step);
+  qr_free();
+}
+static void qr_init(int n) {
+#if defined (S_Plus)
+  if(qr_max_size!=0 || r_ncol!=0 || q_nrow!=0 || q_use_row!=0 ||
+     qmat!=NULL || rmat!=NULL){
+    MESSAGE "Possible memory corruption or memory leak.\n  We"
+      "advise to restart your S+ session"  WARNING(NULL_ENTRY);
+    qr_free();
+  }
+#endif
+  qr_max_size = QR_CHUNK;
+  r_ncol = 0;
+  q_nrow = n;
+  qmat = Calloc(n*qr_max_size,double);
+  if(qmat==NULL)
+    ERRMSG("qr_init", no_dyn_mem_message);
+  rmat = Calloc(qr_max_size*(qr_max_size+1)/2,double);
+  if(rmat==NULL)
+    ERRMSG("qr_init", no_dyn_mem_message);
+}
+static void qr_incr(void) {
+  qr_max_size += QR_CHUNK;
+
+  /* reallocate R always */
+  rmat = Realloc(rmat,qr_max_size*(qr_max_size+1)/2,double);
+  if(rmat==NULL)
+    ERRMSG("qr_incr", no_dyn_mem_message);
+
+  /* reallocate Q only if necessary and only to maximal necessary size */
+  if( qr_max_size >= q_nrow){
+    if( qr_max_size-QR_CHUNK < q_nrow){
+      qmat = Realloc(qmat,q_nrow*q_nrow,double);
+      if(qmat==NULL)
+        ERRMSG("qr_incr", no_dyn_mem_message);
+    }
+  }else{
+    qmat = Realloc(qmat,q_nrow*qr_max_size,double);
+    if(qmat==NULL)
+      ERRMSG("qr_incr", no_dyn_mem_message);
+  }
+}
+static void qr_free(void)
+{
+  qr_max_size = 0;
+  r_ncol = 0;
+  q_nrow = 0;
+  q_use_row = 0;
+  Free(qmat);
+  Free(rmat);
+}
+static void qr_del(int l, int aug) {
+
+  double c, s, tau, nu, *col_k, *col_kp1, *a, *b, tmp;
+  int i, j, k, l0;
+
+  if( l<0 || l>=r_ncol )
+    ERRMSG("qr_del", "Invalid column number");
+  r_ncol--;
+  if( l==r_ncol)
+    if( aug )
+      ERRMSG("qr_del", "Trying to delete last column of augmented matrix");
+    else
+      return;
+
+  /* this is TRUE if $m\le n$ ($m-1\le n$ if the matrix is augmented) */
+  if (r_ncol < q_nrow ){
+    for(k=l; k<r_ncol; k++) /* Update the factorisation and be done */
+    {
+
+col_k = RCOL(k);  /* first element in column $k$ in R */
+Memcpy(col_k,col_k+k+1,k+1);
+a = col_k+k;
+b = a+k+2;
+tau = fabs(*a)+fabs(*b);
+if( tau == 0.0 ) continue; /* both elements are zero
+                              nothing to update */
+nu = tau*sqrt((*a/tau)*(*a/tau)+(*b/tau)*(*b/tau));
+c = *a/nu;
+s = *b/nu;
+*a = nu;
+  b += k+2;
+  a = b-1;
+  for(j=k+2;j<=r_ncol;j++, a+=j, b+=j){
+    tmp = c * *a + s * *b;
+    *b  = c * *b - s * *a;
+    *a  = tmp;
+  }
+  col_k = QCOL(k);
+  col_kp1 = col_k+q_nrow;
+  for(j=0;j<q_nrow;j++){
+    tmp = c*col_k[j]+s*col_kp1[j];
+    col_kp1[j] = c*col_kp1[j]-s*col_k[j];
+    col_k[j] = tmp;
+  }
+    }
+  }else
+    if( l < q_nrow ){
+    /* Update columns upto $m$ and than shift remaining columns */
+      for(k=l; k<q_nrow-1; k++){
+
+col_k = RCOL(k);  /* first element in column $k$ in R */
+Memcpy(col_k,col_k+k+1,k+1);
+a = col_k+k;
+b = a+k+2;
+tau = fabs(*a)+fabs(*b);
+if( tau == 0.0 ) continue; /* both elements are zero
+                              nothing to update */
+nu = tau*sqrt((*a/tau)*(*a/tau)+(*b/tau)*(*b/tau));
+c = *a/nu;
+s = *b/nu;
+*a = nu;
+  b += k+2;
+  a = b-1;
+  for(j=k+2;j<=r_ncol;j++, a+=j, b+=j){
+    tmp = c * *a + s * *b;
+    *b  = c * *b - s * *a;
+    *a  = tmp;
+  }
+  col_k = QCOL(k);
+  col_kp1 = col_k+q_nrow;
+  for(j=0;j<q_nrow;j++){
+    tmp = c*col_k[j]+s*col_kp1[j];
+    col_kp1[j] = c*col_kp1[j]-s*col_k[j];
+    col_k[j] = tmp;
+  }
+      }
+      l0 = q_nrow-1;
+      
+  col_k   = RCOL(l0); /* first element in column $l_0$ in R */
+  for(i=l0; i<r_ncol; i++, col_k +=i)
+     Memcpy(col_k,col_k+i+1,q_nrow);
+
+  /* [[rmat + (q_nrow-1)*q_nrow/2+q_nrow-1]] is last element in column
+     [[q_nrow]] in R */
+  a = rmat + (q_nrow-1)*(q_nrow+2)/2;
+  if( *a < 0 ){
+    /* [[j<r_ncol]] sufficient since we shifted the columns already */
+    for(j=l0;j<r_ncol;j++,a+=j) 
+      *a = - *a;
+    col_k = qmat+(q_nrow-1)*q_nrow;
+    for(i=0;i<q_nrow;i++,col_k++)
+      *col_k = -*col_k;
+  }
+    }else{ /* just shift last columns to the left */
+      l0 = l;
+      
+  col_k   = RCOL(l0); /* first element in column $l_0$ in R */
+  for(i=l0; i<r_ncol; i++, col_k +=i)
+     Memcpy(col_k,col_k+i+1,q_nrow);
+  
+  /* [[rmat + (q_nrow-1)*q_nrow/2+q_nrow-1]] is last element in column
+     [[q_nrow]] in R */
+  a = rmat + (q_nrow-1)*(q_nrow+2)/2;
+  if( *a < 0 ){
+    /* [[j<r_ncol]] sufficient since we shifted the columns already */
+    for(j=l0;j<r_ncol;j++,a+=j)
+      *a = - *a;
+    col_k = qmat+(q_nrow-1)*q_nrow;
+    for(i=0;i<q_nrow;i++,col_k++)
+      *col_k = -*col_k;
+  }
+    }
+}
+static void qr_add(double *x, int swap){
+
+  double tmp, norm_orig, norm_init, norm_last, *q_new_col;
+  int i;
+  double *col_l, *col_lm1, *q_elem, *r_new_col;
+  double norm_new;
+  int j;
+  double c, s, tau, nu, *a, *b;
+
+  if( r_ncol == qr_max_size )
+    qr_incr();
+
+norm_orig = 0.0;
+tmp = 0.0;
+for(i=0; i<q_nrow; i++) tmp += x[i]*x[i];
+norm_orig = sqrt(tmp);
+if(r_ncol<q_nrow){
+  q_new_col = QCOL(r_ncol); /* points to the first element of new column */
+  tmp = 0.0;
+  for(i=0; i<q_nrow; i++){
+    q_new_col[i] = x[i]/norm_orig;
+    tmp += q_new_col[i]*q_new_col[i];
+  }
+  if(r_ncol==0){
+    *rmat = norm_orig;
+    r_ncol++;
+    return;
+  }
+  norm_init = norm_last = sqrt(tmp);
+}
+r_new_col = RCOL(r_ncol);
+for(i=0; i<=r_ncol; i++) r_new_col[i] = 0.0;
+if( r_ncol >= q_nrow ){
+  /* Multiply new column with Q-transpose in [[qr_add]] */
+  q_elem = qmat;
+  for(i=0;i<q_nrow;i++){
+    tmp = 0.0;
+    for(j=0;j<q_nrow;j++,q_elem++)
+      tmp += *q_elem * x[j];
+    r_new_col[i] = tmp;
+  }
+  if( swap ){ /* swap last two columns */
+    col_l   = r_new_col;
+    col_lm1 = r_new_col-r_ncol;
+    for(i=0; i<q_nrow; i++,col_l++,col_lm1++){
+      tmp = *col_l;
+      *col_l = *col_lm1;
+      *col_lm1 = tmp;
+    }
+    col_lm1--;
+    col_l--;
+    if(q_nrow==r_ncol && *col_lm1<0){
+      /* if new R[$n$,$n$] is negative, than change signs */
+      *col_lm1 = - *col_lm1;
+      *col_l   = - *col_l;
+      q_new_col = qmat+(q_nrow-1)*q_nrow;
+      for(j=0;j<q_nrow;j++)
+        q_new_col[j] = -q_new_col[j];
+    }
+  }
+  r_ncol++;
+  return;
+}
+  while(1){
+
+q_elem = qmat;
+for(j=0; j<r_ncol;j++){
+  tmp = 0.0;
+  for(i=0;i<q_nrow;i++,q_elem++)
+    tmp += *q_elem * q_new_col[i];
+  r_new_col[j] += tmp;
+  q_elem -= q_nrow;
+  for(i=0;i<q_nrow;i++,q_elem++)
+    q_new_col[i] -= *q_elem*tmp;
+}
+tmp = 0.0;
+for(i=0;i<q_nrow;i++) tmp += q_new_col[i]*q_new_col[i];
+norm_new = sqrt(tmp);
+
+  if( norm_new >= norm_last/2.0) break;
+  if( norm_new > 0.1*norm_init*PRECISION )
+    norm_last = norm_new;
+  else{
+    norm_init = norm_last = 0.1*norm_last*PRECISION;
+    for(i=0;i<q_nrow;i++)
+      q_new_col[i] = 0.0;
+    if(q_use_row == q_nrow)
+      ERRMSG("qr_add","Cannot orthogonalise new column\n");
+
+    q_new_col[q_use_row] = norm_new = norm_last;
+    q_use_row++;
+   }
+  }
+  
+  for(i=0;i<q_nrow;i++)
+    q_new_col[i] /= norm_new;
+  for(i=0;i<r_ncol;i++)
+    r_new_col[i] *= norm_orig;
+  r_new_col[r_ncol] = norm_new*norm_orig;
+  if(swap){
+
+col_l = r_new_col; /* first element in last column in R */
+col_lm1 = r_new_col-r_ncol; /* first element in column before last in R */
+for(j=0;j<r_ncol;j++,col_l++,col_lm1++){
+  tmp = *col_l;
+  *col_l = *col_lm1;
+  *col_lm1 = tmp;
+}
+a = col_lm1-1;
+b = col_l;
+tau = fabs(*a)+fabs(*b);
+if( tau > 0.0 ){
+  /*Calculate Givens rotation*/
+  nu = tau*sqrt((*a/tau)*(*a/tau)+(*b/tau)*(*b/tau));
+  c = *a/nu;
+  s = *b/nu;
+
+  *a = nu;   /* Fix column before last in R */
+  *b = -s* *(b-1);  /* Fix last element of last column in R */
+  *(b-1) = c * *(b-1); /* Fix second last element of last column in R */
+  if( *b < 0){
+    *b = -*b;
+    
+  col_lm1 = QCOL(r_ncol-1);
+  col_l   = col_lm1+q_nrow;
+  for(j=0;j<q_nrow;j++){
+    tmp = c*col_lm1[j]+s*col_l[j];
+    col_l[j] = -c*col_l[j]+s*col_lm1[j];
+    col_lm1[j] = tmp;
+  }
+  
+  }else{
+    
+  col_lm1 = QCOL(r_ncol-1);
+  col_l   = col_lm1+q_nrow;
+  for(j=0;j<q_nrow;j++){
+    tmp = c*col_lm1[j]+s*col_l[j];
+    col_l[j] = c*col_l[j]-s*col_lm1[j];
+    col_lm1[j] = tmp;
+  }
+  }
+}
+
+  }
+  r_ncol++;
+
+}
+#if defined (S_Plus)
+static void errmsg(char *string){
+  PROBLEM "%s\n", string RECOVER(NULL_ENTRY);
+}
+#elif defined(Matlab)
+static void errmsg(char *where, char *string){
+  char str[1024];
+  sprintf(str, "Error in %s: %s\n", where, string);
+  mexErrMsgTxt(str);
+}
+#else
+static void errmsg(char *where, char *string){
+  fprintf(stderr, "Error in %s: %s\n", where, string);
+  exit(EXIT_FAILURE);
+}
+#endif
+
+int assertFail(const char *ex, const char *file, const char *func, const int line)
+{
+static bool ignore = false;
+
+    if (!ignore) {
+        char str[1024];
+        sprintf(str, "%s:%u: %s(): failed assertion `%s'\n",
+                     file, line, func, ex);
+        mexErrMsgTxt(str);
+    } /* if */
+    return(1);
+}
+
+static void fortifyPrintf(const char *s)
+{
+   mexPrintf(s);
+}
+
+void
+mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
+{
+   mxArray *array_ptr;
+   int i;
+   long dataSize;
+   long dimensions;
+   long succeed;
+   long verbose = 0;
+   long as_sub = 0;
+   double lagrangian;
+   double threshold;
+   double *mx;
+   double *y;
+   double *yhat1;
+   double *r;
+   double *beta;
+   int arg;
+
+#ifdef FORTIFY
+   Fortify_SetOutputFunc(fortifyPrintf);
+   Fortify_EnterScope();
+#endif
+
+   if (nrhs > 3) {
+      mexErrMsgTxt("Syntax error: too many arguments\n");
+      return;
+   } /* if */
+
+   if (nrhs < 3) {
+      mexErrMsgTxt("Syntax error: arguments missing\n");
+      return;
+   } /* if */
+
+   if (nlhs > 1) {
+      mexErrMsgTxt("Syntax error: too many return arguments\n");
+      return;
+   } /* if */
+
+   arg = 2;
+   array_ptr = (mxArray *) prhs[arg];
+   if (mxGetClassID(array_ptr) != mxDOUBLE_CLASS) {
+      mexErrMsgTxt("Threshold argument is not a number\n");
+      return;
+   } /* if */
+
+   threshold = mxGetScalar(array_ptr);
+
+   arg = 1;
+   array_ptr = (mxArray *) prhs[arg];
+   if (mxGetClassID(array_ptr) != mxDOUBLE_CLASS) {
+      mexErrMsgTxt("Label argument is not a number\n");
+      return;
+   } /* if */
+
+   y = mxGetPr(array_ptr);
+
+   arg = 0;
+   array_ptr = (mxArray *) prhs[arg];
+   if (mxGetClassID(array_ptr) != mxDOUBLE_CLASS) {
+      mexErrMsgTxt("Data set argument is not a double matrix\n");
+      return;
+   } /* if */
+
+   dataSize = mxGetM(array_ptr);
+   dimensions = mxGetN(array_ptr);
+
+   mx = mxGetPr(array_ptr);
+
+   arg = 0;
+   plhs[arg] = mxCreateDoubleMatrix(dimensions, 1, mxREAL);
+
+   yhat1 = new double[dataSize];
+   r = new double[dataSize];
+   beta = mxGetPr(plhs[arg]);
+
+   for (i=0; i<dimensions; i++) {
+       beta[i] = 0.0;
+   } /* for */
+
+   if (dimensions > 0) {
+       lasso(mx, &dataSize, &dimensions, &threshold, beta, y,
+             yhat1, r, &lagrangian, &succeed, &verbose, &as_sub);
+   } else {
+       mexPrintf("Dataset with 0 dimensions\n");
+   } /* if */
+
+   delete [] yhat1;
+   delete [] r;
+
+#ifdef  FORTIFY
+   Fortify_LeaveScope();
+#endif
+
+}
+
Index: /prextra/prextra/mexlasso.cc
===================================================================
--- /prextra/prextra/mexlasso.cc	(revision 2)
+++ /prextra/prextra/mexlasso.cc	(revision 2)
@@ -0,0 +1,1003 @@
+/* Copyright (C) 1998
+   Berwin A Turlach <bturlach@stats.adelaide.edu.au> */
+
+/* This library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Library General Public License
+   as published by the Free Software Foundation; either version 2 of
+   the License, or (at your option) any later version. */
+
+/* This library is distributed in the hope that it will be useful, but
+   WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Library General Public License for more details. */
+
+/* You should have received a copy of the GNU Library General Public
+   License along with this library; if not, write to the Free Software
+   Foundation, Inc., 59 Temple Place, Suite 330, Boston,
+   MA 02111-1307, USA. */
+
+#include "lasso.h"
+//#include "fortify.h"
+
+#define PRECISION FLT_EPSILON
+
+static void lasso_alloc(long n, long m);
+static void lasso_free(void);
+static void qr_init(int n);
+static void qr_incr(void);
+static void qr_free(void);
+static void qr_del(int l, int aug);
+static void qr_add(double *x, int swap);
+#if defined (S_Plus)
+static void errmsg(char* string);
+#else
+static void errmsg(char* where, char* string);
+#endif
+
+  static int QR_CHUNK = 10;
+
+  static double *xtr=NULL, *btmp=NULL, *qtr=NULL,
+                *rinvt_theta=NULL, *step=NULL, ytyd2=0.0;
+  static int *theta=NULL, *nz_x=NULL, num_nz_x=0;
+  static double *qmat=NULL, *rmat=NULL;
+  static int qr_max_size=0, r_ncol=0, q_nrow=0, q_use_row=0;
+  static char *no_dyn_mem_message="Cannot allocate dynamic memory";
+
+void lasso(double *x,           /* input: data */
+           long *pn,            /* input: rows */
+           long *pm,            /* input: columns */
+           double *pt,          /* input: threshold */
+           double *beta,        /* input/ouput: intial/final beta */
+           double *y,           /* input: labels */
+
+           double *yhat1,       /* output */
+           double *r,           /* output */
+           double *lagrangian,  /* output */
+           long *psuc,          /* output: succeed: 0/-1 */
+           long *pverb,         /* input: verbose: yes/no */
+           long *pas_sub)       /* input: no!? */
+{
+
+  double t = *pt, prec;
+  long n = *pn, m = *pm, verb = *pverb, as_sub = *pas_sub;
+  int not_solved;
+  double *x_elem=NULL, tmp, max_val, b_1norm;
+  int i, j, max_ind;
+  double p_obj, d_obj;
+  int num_iter=0, max_num;
+  double b_1norm_old, mu;
+  double rho_up, rho_low, rho;
+  int to_del;
+  double *q_elem, wtc, wtw;
+  double *r_elem;
+  int k;
+  int add;
+  if( !as_sub)
+    lasso_alloc(n,m);
+
+  prec = sqrt(PRECISION);
+
+  if( as_sub ){
+
+  b_1norm = 0.0;
+  for(j=0;j<num_nz_x;j++)
+    b_1norm += fabs(beta[nz_x[j]]);
+  }else{
+
+  b_1norm = 0.0;
+  num_nz_x = 0;
+  for(j=0; j<m; j++){
+    if(fabs(beta[j]) > prec){
+      b_1norm += fabs(beta[j]);
+      nz_x[num_nz_x] = j;
+      num_nz_x++;
+    }else
+      beta[j] = 0.0;
+  }
+  }
+
+  if( b_1norm > t){
+    if(verb){
+      mexPrintf("******************************\n");
+      mexPrintf("Rescaling beta from L1-norm %f to %f\n",b_1norm,t);
+      mexEvalString("drawnow;");
+
+    }
+    for(j=0; j<num_nz_x; j++)
+      beta[nz_x[j]] = beta[nz_x[j]] * t/b_1norm;
+    b_1norm = t;
+  }
+
+  for(i=0; i< n; i++)
+    yhat1[i] = 0.0;
+  for(j=0; j < num_nz_x; j++){
+    /* x_elem points to the first element in the column of X to which
+       the j-th entry in nz_x points  */
+    x_elem = x + nz_x[j]*n;
+    tmp = beta[nz_x[j]];
+    for(i=0; i < n; i++){
+      yhat1[i] += *x_elem*tmp;
+      x_elem++;           /* now we point to the next element in X */
+    }
+  }
+
+  /* calculate the residual vector */
+  for(i=0; i< n; i++)
+    r[i] = y[i]-yhat1[i];
+
+  /* multiply X^T with the residual vector */
+  x_elem = x;
+  for(j=0; j < m; j++){
+    tmp = 0.0;
+    for(i=0; i<n; i++){
+      tmp += *x_elem*r[i];
+      x_elem++;
+    }
+    xtr[j] = tmp;
+  }
+
+  max_val = fabs(xtr[0]);
+  max_ind = 0;
+  for(j=1; j<m; j++)
+    if( fabs(xtr[j]) > max_val ){
+      max_val = fabs(xtr[j]);
+      max_ind = j;
+    }
+
+  if( !as_sub ){
+
+  qr_add(y,TRUE);
+  ytyd2 = *rmat * *rmat/2.0;
+  for(j=0;j<num_nz_x;j++){
+    qr_add(x+nz_x[j]*n, TRUE);
+    if(fabs(beta[nz_x[j]])<prec)
+      theta[j] = xtr[nz_x[j]] < 0 ? -1 : 1;
+    else
+      theta[j] = beta[nz_x[j]] < 0 ? -1 : 1;
+
+  }
+  }
+  if( num_nz_x==0 ){
+    nz_x[0] = max_ind;
+    num_nz_x = 1;
+    if(verb){
+      mexPrintf("******************************\n");
+      mexPrintf("  -->\tAdding variable: %d\n",max_ind+1);
+      mexEvalString("drawnow;");
+    }
+    qr_add(x+max_ind*n, TRUE);
+    theta[0] = xtr[max_ind] < 0 ? -1 : 1;
+  }
+  *psuc=0;
+  if(verb){
+
+
+  /* Find out how many times [[max_val]] is attained */
+  tmp = (1.0-prec)*max_val;
+  if( tmp < prec ) tmp = 0.0;
+  max_num = 1;
+  for(j=0; j<m; j++){
+    if( tmp <= fabs(xtr[j]) && j!=max_ind){
+      /*we found another element equal to the (current)
+        maximal absolute value */
+      max_num++;
+    }
+  }
+
+  /* for the value of the dual objective function we need
+     to calculate the L2 norm of the vector of fitted values */
+  p_obj = 0.0;
+  d_obj = 0.0;
+  for(i=0;i<n;i++){
+    p_obj += r[i]*r[i];
+    d_obj += yhat1[i]*yhat1[i];
+  }
+  p_obj /= 2.0;
+  d_obj = ytyd2 - d_obj/2.0 - t*max_val ;
+  mexPrintf("******************************\n");
+  mexPrintf("\nIteration number: %d\n", num_iter);
+  mexPrintf("Value of primal object function      : %f\n", p_obj);
+  mexPrintf("Value of dual   object function      : %f\n", d_obj);
+  mexPrintf("L1 norm of current beta              : %f", b_1norm);
+  mexPrintf(" <= %f\n", t);
+  mexPrintf("Maximal absolute value in t(X)%%*%%r   : %e", max_val);
+  mexPrintf(" attained %d time(s)\n", max_num);
+  mexPrintf("Number of parameters allowed to vary : %d\n", num_nz_x);
+  mexEvalString("drawnow;");
+  num_iter++;
+  }
+  while(1){
+    do{
+
+  q_elem = qmat;
+
+  for(j=0;j<num_nz_x;j++){
+    tmp = 0.0;
+    for(i=0;i<n;i++,q_elem++) {
+      tmp += *q_elem*r[i];
+    }
+    qtr[j] = tmp;
+  }
+
+  if( b_1norm < (1.0-prec)*t )
+    mu = 0.0;
+  else{
+
+  /* z=R^{-T}\theta can be calculated by solving R^Tz=\theta */
+  r_elem = rmat;
+  for(j=0;j<num_nz_x;j++,r_elem++){
+    tmp = theta[j];
+    for(k=0;k<j;k++,r_elem++)
+      tmp -= *r_elem * rinvt_theta[k];
+    rinvt_theta[j] = tmp /  *r_elem;
+  }
+
+  wtc = 0.0;
+  wtw = 0.0;
+  for(j=0;j<num_nz_x;j++){
+    wtc += rinvt_theta[j]*qtr[j];
+    wtw += rinvt_theta[j]*rinvt_theta[j];
+  }
+  mu = (wtc-(t-b_1norm))/wtw;
+  }
+
+  if(mu<=0.0)
+    for(j=0;j<num_nz_x;j++)
+      step[j] = qtr[j];
+  else
+    for(j=0;j<num_nz_x;j++)
+      step[j] = qtr[j] - mu*rinvt_theta[j];
+
+  /* h=R^{-1}z can be calculated by solving Rh=z */
+  for(j=num_nz_x-1;j>=0;j--){
+    tmp = step[j];
+    for(k=num_nz_x-1;k>j;k--)
+      tmp -= RMAT(j,k) * step[k];
+    step[j] = tmp /  RMAT(j,j);
+  }
+  for(j=0;j<num_nz_x;j++){
+    btmp[j] = beta[nz_x[j]];
+    beta[nz_x[j]] += step[j];
+  }
+  b_1norm_old=b_1norm;
+
+  b_1norm = 0.0;
+  for(j=0;j<num_nz_x;j++) {
+    b_1norm += fabs(beta[nz_x[j]]);
+  }
+
+  not_solved=FALSE;
+  if( b_1norm > (1+prec)*t){
+    not_solved=TRUE;
+    if(b_1norm_old < (1.0-prec)*t){
+
+    if(verb) {
+      mexPrintf("  -->\tStepping onto the border of the L1 ball.\n");
+      mexEvalString("drawnow;");
+    }
+    rho_up  = rho = 1.0;
+    rho_low = 0.0;
+    while( fabs(t-b_1norm) > prec*t ){
+      if( b_1norm > t){
+        rho_up = rho;
+        rho = (rho+rho_low)/2.0;
+      }
+      if( b_1norm < t){
+        rho_low = rho;
+        rho = (rho+rho_up)/2.0;
+      }
+      if(rho < prec) break;
+      for(j=0; j<num_nz_x; j++) {
+        beta[nz_x[j]] = btmp[j]+rho*step[j];
+      }
+
+  b_1norm = 0.0;
+  for(j=0;j<num_nz_x;j++)
+    b_1norm += fabs(beta[nz_x[j]]);
+    }
+    for(j=0;j<num_nz_x;j++){
+      if(fabs(beta[nz_x[j]])<prec)
+        theta[j] = btmp[j] > 0 ? -1 : 1;
+      else
+        theta[j] = beta[nz_x[j]] < 0 ? -1 : 1;
+    }
+    }else{
+
+
+  rho = 1.0;
+  to_del = -1;
+  for(j=0; j<num_nz_x; j++){
+    if(fabs(step[j]) > prec){
+      tmp = -btmp[j]/(step[j]);
+      if( 0.0 < tmp && tmp < rho ){
+        rho = tmp;
+        to_del = j;
+      }
+    }
+  }
+  if(to_del < 0 ){
+    *psuc= -1;
+    goto EXIT_HERE;
+  }
+  for(j=0; j<num_nz_x; j++) {
+    beta[nz_x[j]] = btmp[j]+rho*step[j];
+  }
+  if(verb) {
+    mexPrintf("  -->\tRemoving variable: %d",nz_x[to_del]+1);
+    mexEvalString("drawnow;");
+  }
+  beta[nz_x[to_del]] = 0.0;
+  qr_del(to_del,TRUE);
+  for(j=to_del+1; j< num_nz_x; j++){
+    nz_x[j-1] = nz_x[j];
+    theta[j-1] = theta[j];
+  }
+  num_nz_x--;
+  
+  b_1norm = 0.0;
+  for(j=0;j<num_nz_x;j++)
+    b_1norm += fabs(beta[nz_x[j]]);
+  if(verb)
+    if(b_1norm < (1-prec)*t)
+      mexPrintf(", and stepping into the interior of the L1 ball\n");
+    else
+      mexPrintf("\n");
+    }
+    mexEvalString("drawnow;");
+  }
+
+  for(i=0; i< n; i++)
+    yhat1[i] = 0.0;
+  for(j=0; j < num_nz_x; j++){
+    /* x_elem points to the first element in the column of X to which
+       the j-th entry in nz_x points  */
+    x_elem = x + nz_x[j]*n;
+    tmp = beta[nz_x[j]];
+    for(i=0; i < n; i++){
+      yhat1[i] += *x_elem*tmp;
+      x_elem++;           /* now we point to the next element in X */
+    }
+  }
+
+  /* calculate the residual vector */
+  for(i=0; i< n; i++) {
+    r[i] = y[i]-yhat1[i];
+  }
+    }while(not_solved);
+    
+  
+  for(i=0; i< n; i++)
+    yhat1[i] = 0.0;
+  for(j=0; j < num_nz_x; j++){
+    /* x_elem points to the first element in the column of X to which
+       the j-th entry in nz_x points  */
+    x_elem = x + nz_x[j]*n;
+    tmp = beta[nz_x[j]];
+    for(i=0; i < n; i++){
+      yhat1[i] += *x_elem*tmp;
+      x_elem++;           /* now we point to the next element in X */
+    }
+  }
+
+  /* calculate the residual vector */
+  for(i=0; i< n; i++)
+    r[i] = y[i]-yhat1[i];
+
+
+  /* multiply X^T with the residual vector */
+  x_elem = x;
+  for(j=0; j < m; j++){
+    tmp = 0.0;
+    for(i=0; i<n; i++){
+      tmp += *x_elem*r[i];
+      x_elem++;
+    }
+    xtr[j] = tmp;
+  }
+  max_val = fabs(xtr[0]);
+  max_ind = 0;
+  for(j=1; j<m; j++)
+    if( fabs(xtr[j]) > max_val ){
+      max_val = fabs(xtr[j]);
+      max_ind = j;
+    }
+
+    if(verb){
+
+
+  /* Find out how many times [[max_val]] is attained */
+  tmp = (1.0-prec)*max_val;
+  if( tmp < prec ) tmp = 0.0;
+  max_num = 1;
+  for(j=0; j<m; j++){
+    if( tmp <= fabs(xtr[j]) && j!=max_ind){
+      /*we found another element equal to the (current)
+        maximal absolute value */
+      max_num++;
+    }
+  }
+
+  /* for the value of the dual objective function we need
+     to calculate the L2 norm of the vector of fitted values */
+  p_obj = 0.0;
+  d_obj = 0.0;
+  for(i=0;i<n;i++){
+    p_obj += r[i]*r[i];
+    d_obj += yhat1[i]*yhat1[i];
+  }
+  p_obj /= 2.0;
+  d_obj = ytyd2 - d_obj/2.0 - t*max_val ;
+  mexPrintf("******************************\n");
+  mexPrintf("\nIteration number: %d\n", num_iter);
+  mexPrintf("Value of primal object function      : %f\n", p_obj);
+  mexPrintf("Value of dual   object function      : %f\n", d_obj);
+  mexPrintf("L1 norm of current beta              : %f", b_1norm);
+  mexPrintf(" <= %f\n", t);
+  mexPrintf("Maximal absolute value in t(X)%%*%%r   : %e", max_val);
+  mexPrintf(" attained %d time(s)\n", max_num);
+  mexPrintf("Number of parameters allowed to vary : %d\n", num_nz_x);
+  mexEvalString("drawnow;");
+  num_iter++;
+    }
+    
+  add = TRUE;
+  for(i=0; i<num_nz_x; i++)
+    if(nz_x[i]==max_ind){
+      add = FALSE;
+      break;
+    }
+  if(add){
+    qr_add(x+max_ind*n,TRUE);
+    nz_x[num_nz_x] = max_ind;
+    theta[num_nz_x] = xtr[max_ind]<0 ? -1 : 1;
+    num_nz_x++;
+  }else{
+    *lagrangian = max_val;
+    break;
+  }
+
+  b_1norm = 0.0;
+  for(j=0;j<num_nz_x;j++)
+    b_1norm += fabs(beta[nz_x[j]]);
+  if(verb) {
+    mexPrintf("  -->\tAdding variable: %d\n",max_ind+1);
+    mexEvalString("drawnow;");
+  }
+  }
+  EXIT_HERE:
+  if( !as_sub)
+    lasso_free();
+}
+void mult_lasso(double *x, long *pn, long *pm, double *pt, long *pl,
+                double *beta, double *y, double *yhat1, double *r,
+                double *lagrangian, long *psuc, long *pverb)
+{
+
+  double prec;
+  long n = *pn, m = *pm, l = *pl, verb = *pverb, as_sub = TRUE, i, j;
+  lasso_alloc(n,m);
+
+qr_add(y,TRUE);
+ytyd2 = *rmat * *rmat/2.0;
+prec = sqrt(PRECISION);
+num_nz_x = 0;
+for(j=0; j<m; j++){
+  if(fabs(beta[j]) > prec){
+    qr_add(x+j*n, TRUE);
+    nz_x[num_nz_x] = j;
+    num_nz_x++;
+  }else
+    beta[j] = 0.0;
+}
+*psuc = 0;
+  for(i=0; i<l; i++){
+
+if(verb){
+  mexPrintf("\n\n++++++++++++++++++++++++++++++\n");
+  mexPrintf("Solving problem number %ld with bound %f\n", i+1, pt[i]);
+  mexPrintf("++++++++++++++++++++++++++++++\n");
+  mexEvalString("drawnow;");
+}
+if(i>0) Memcpy(beta,beta-m,m);
+lasso(x, pn, pm, pt+i, beta, y, yhat1, r, lagrangian, psuc, pverb, &as_sub);
+if( *psuc < 0 ){
+  goto EXIT_HERE;
+}
+
+beta += m;
+yhat1 += n;
+r += n;
+lagrangian++;
+  }
+  EXIT_HERE:
+  lasso_free();
+}
+static void lasso_alloc(long n, long m){
+
+#if defined(S_Plus)
+  if( nz_x != NULL || theta != NULL || xtr != NULL || btmp != NULL ||
+      qtr != NULL || rinvt_theta != NULL || step != NULL ||
+      num_nz_x != 0 || ytyd2 != 0.0){
+    MESSAGE "Possible memory corruption or memory leak.\n  We"
+      "advise to restart your S+ session"  WARNING(NULL_ENTRY);
+     lasso_free();
+  }
+#endif
+
+  QR_CHUNK = m;  // Added: CJV
+
+  nz_x = Calloc(m,int);
+  if( nz_x == NULL )
+    ERRMSG("lasso_alloc", no_dyn_mem_message);
+  theta = Calloc(m,int);
+  if( theta==NULL )
+    ERRMSG("lasso_alloc", no_dyn_mem_message);
+  xtr = Calloc(m,double);
+  if( xtr==NULL )
+    ERRMSG("lasso_alloc", no_dyn_mem_message);
+  btmp = Calloc(m,double);
+  if( btmp==NULL )
+    ERRMSG("lasso_alloc", no_dyn_mem_message);
+  qtr = Calloc(m,double);
+  if( qtr==NULL )
+    ERRMSG("lasso_alloc", no_dyn_mem_message);
+  rinvt_theta = Calloc(m,double);
+  if( rinvt_theta==NULL )
+    ERRMSG("lasso_alloc", no_dyn_mem_message);
+  step = Calloc(m,double);
+  if( step==NULL )
+    ERRMSG("lasso_alloc", no_dyn_mem_message);
+  qr_init(n);
+}
+static void lasso_free(void){
+  num_nz_x=0;
+  ytyd2 = 0.0;
+  Free(nz_x);
+  Free(theta);
+  Free(xtr);
+  Free(btmp);
+  Free(qtr);
+  Free(rinvt_theta);
+  Free(step);
+  qr_free();
+}
+static void qr_init(int n) {
+#if defined (S_Plus)
+  if(qr_max_size!=0 || r_ncol!=0 || q_nrow!=0 || q_use_row!=0 ||
+     qmat!=NULL || rmat!=NULL){
+    MESSAGE "Possible memory corruption or memory leak.\n  We"
+      "advise to restart your S+ session"  WARNING(NULL_ENTRY);
+    qr_free();
+  }
+#endif
+  qr_max_size = QR_CHUNK;
+  r_ncol = 0;
+  q_nrow = n;
+  qmat = Calloc(n*qr_max_size,double);
+  if(qmat==NULL)
+    ERRMSG("qr_init", no_dyn_mem_message);
+  rmat = Calloc(qr_max_size*(qr_max_size+1)/2,double);
+  if(rmat==NULL)
+    ERRMSG("qr_init", no_dyn_mem_message);
+}
+static void qr_incr(void) {
+  qr_max_size += QR_CHUNK;
+
+  /* reallocate R always */
+  rmat = Realloc(rmat,qr_max_size*(qr_max_size+1)/2,double);
+  if(rmat==NULL)
+    ERRMSG("qr_incr", no_dyn_mem_message);
+
+  /* reallocate Q only if necessary and only to maximal necessary size */
+  if( qr_max_size >= q_nrow){
+    if( qr_max_size-QR_CHUNK < q_nrow){
+      qmat = Realloc(qmat,q_nrow*q_nrow,double);
+      if(qmat==NULL)
+        ERRMSG("qr_incr", no_dyn_mem_message);
+    }
+  }else{
+    qmat = Realloc(qmat,q_nrow*qr_max_size,double);
+    if(qmat==NULL)
+      ERRMSG("qr_incr", no_dyn_mem_message);
+  }
+}
+static void qr_free(void)
+{
+  qr_max_size = 0;
+  r_ncol = 0;
+  q_nrow = 0;
+  q_use_row = 0;
+  Free(qmat);
+  Free(rmat);
+}
+static void qr_del(int l, int aug) {
+
+  double c, s, tau, nu, *col_k, *col_kp1, *a, *b, tmp;
+  int i, j, k, l0;
+
+  if( l<0 || l>=r_ncol )
+    ERRMSG("qr_del", "Invalid column number");
+  r_ncol--;
+  if( l==r_ncol)
+    if( aug )
+      ERRMSG("qr_del", "Trying to delete last column of augmented matrix");
+    else
+      return;
+
+  /* this is TRUE if $m\le n$ ($m-1\le n$ if the matrix is augmented) */
+  if (r_ncol < q_nrow ){
+    for(k=l; k<r_ncol; k++) /* Update the factorisation and be done */
+    {
+
+col_k = RCOL(k);  /* first element in column $k$ in R */
+Memcpy(col_k,col_k+k+1,k+1);
+a = col_k+k;
+b = a+k+2;
+tau = fabs(*a)+fabs(*b);
+if( tau == 0.0 ) continue; /* both elements are zero
+                              nothing to update */
+nu = tau*sqrt((*a/tau)*(*a/tau)+(*b/tau)*(*b/tau));
+c = *a/nu;
+s = *b/nu;
+*a = nu;
+  b += k+2;
+  a = b-1;
+  for(j=k+2;j<=r_ncol;j++, a+=j, b+=j){
+    tmp = c * *a + s * *b;
+    *b  = c * *b - s * *a;
+    *a  = tmp;
+  }
+  col_k = QCOL(k);
+  col_kp1 = col_k+q_nrow;
+  for(j=0;j<q_nrow;j++){
+    tmp = c*col_k[j]+s*col_kp1[j];
+    col_kp1[j] = c*col_kp1[j]-s*col_k[j];
+    col_k[j] = tmp;
+  }
+    }
+  }else
+    if( l < q_nrow ){
+    /* Update columns upto $m$ and than shift remaining columns */
+      for(k=l; k<q_nrow-1; k++){
+
+col_k = RCOL(k);  /* first element in column $k$ in R */
+Memcpy(col_k,col_k+k+1,k+1);
+a = col_k+k;
+b = a+k+2;
+tau = fabs(*a)+fabs(*b);
+if( tau == 0.0 ) continue; /* both elements are zero
+                              nothing to update */
+nu = tau*sqrt((*a/tau)*(*a/tau)+(*b/tau)*(*b/tau));
+c = *a/nu;
+s = *b/nu;
+*a = nu;
+  b += k+2;
+  a = b-1;
+  for(j=k+2;j<=r_ncol;j++, a+=j, b+=j){
+    tmp = c * *a + s * *b;
+    *b  = c * *b - s * *a;
+    *a  = tmp;
+  }
+  col_k = QCOL(k);
+  col_kp1 = col_k+q_nrow;
+  for(j=0;j<q_nrow;j++){
+    tmp = c*col_k[j]+s*col_kp1[j];
+    col_kp1[j] = c*col_kp1[j]-s*col_k[j];
+    col_k[j] = tmp;
+  }
+      }
+      l0 = q_nrow-1;
+      
+  col_k   = RCOL(l0); /* first element in column $l_0$ in R */
+  for(i=l0; i<r_ncol; i++, col_k +=i)
+     Memcpy(col_k,col_k+i+1,q_nrow);
+
+  /* [[rmat + (q_nrow-1)*q_nrow/2+q_nrow-1]] is last element in column
+     [[q_nrow]] in R */
+  a = rmat + (q_nrow-1)*(q_nrow+2)/2;
+  if( *a < 0 ){
+    /* [[j<r_ncol]] sufficient since we shifted the columns already */
+    for(j=l0;j<r_ncol;j++,a+=j) 
+      *a = - *a;
+    col_k = qmat+(q_nrow-1)*q_nrow;
+    for(i=0;i<q_nrow;i++,col_k++)
+      *col_k = -*col_k;
+  }
+    }else{ /* just shift last columns to the left */
+      l0 = l;
+      
+  col_k   = RCOL(l0); /* first element in column $l_0$ in R */
+  for(i=l0; i<r_ncol; i++, col_k +=i)
+     Memcpy(col_k,col_k+i+1,q_nrow);
+  
+  /* [[rmat + (q_nrow-1)*q_nrow/2+q_nrow-1]] is last element in column
+     [[q_nrow]] in R */
+  a = rmat + (q_nrow-1)*(q_nrow+2)/2;
+  if( *a < 0 ){
+    /* [[j<r_ncol]] sufficient since we shifted the columns already */
+    for(j=l0;j<r_ncol;j++,a+=j)
+      *a = - *a;
+    col_k = qmat+(q_nrow-1)*q_nrow;
+    for(i=0;i<q_nrow;i++,col_k++)
+      *col_k = -*col_k;
+  }
+    }
+}
+static void qr_add(double *x, int swap){
+
+  double tmp, norm_orig, norm_init, norm_last, *q_new_col;
+  int i;
+  double *col_l, *col_lm1, *q_elem, *r_new_col;
+  double norm_new;
+  int j;
+  double c, s, tau, nu, *a, *b;
+
+  if( r_ncol == qr_max_size )
+    qr_incr();
+
+norm_orig = 0.0;
+tmp = 0.0;
+for(i=0; i<q_nrow; i++) tmp += x[i]*x[i];
+norm_orig = sqrt(tmp);
+if(r_ncol<q_nrow){
+  q_new_col = QCOL(r_ncol); /* points to the first element of new column */
+  tmp = 0.0;
+  for(i=0; i<q_nrow; i++){
+    q_new_col[i] = x[i]/norm_orig;
+    tmp += q_new_col[i]*q_new_col[i];
+  }
+  if(r_ncol==0){
+    *rmat = norm_orig;
+    r_ncol++;
+    return;
+  }
+  norm_init = norm_last = sqrt(tmp);
+}
+r_new_col = RCOL(r_ncol);
+for(i=0; i<=r_ncol; i++) r_new_col[i] = 0.0;
+if( r_ncol >= q_nrow ){
+  /* Multiply new column with Q-transpose in [[qr_add]] */
+  q_elem = qmat;
+  for(i=0;i<q_nrow;i++){
+    tmp = 0.0;
+    for(j=0;j<q_nrow;j++,q_elem++)
+      tmp += *q_elem * x[j];
+    r_new_col[i] = tmp;
+  }
+  if( swap ){ /* swap last two columns */
+    col_l   = r_new_col;
+    col_lm1 = r_new_col-r_ncol;
+    for(i=0; i<q_nrow; i++,col_l++,col_lm1++){
+      tmp = *col_l;
+      *col_l = *col_lm1;
+      *col_lm1 = tmp;
+    }
+    col_lm1--;
+    col_l--;
+    if(q_nrow==r_ncol && *col_lm1<0){
+      /* if new R[$n$,$n$] is negative, than change signs */
+      *col_lm1 = - *col_lm1;
+      *col_l   = - *col_l;
+      q_new_col = qmat+(q_nrow-1)*q_nrow;
+      for(j=0;j<q_nrow;j++)
+        q_new_col[j] = -q_new_col[j];
+    }
+  }
+  r_ncol++;
+  return;
+}
+  while(1){
+
+q_elem = qmat;
+for(j=0; j<r_ncol;j++){
+  tmp = 0.0;
+  for(i=0;i<q_nrow;i++,q_elem++)
+    tmp += *q_elem * q_new_col[i];
+  r_new_col[j] += tmp;
+  q_elem -= q_nrow;
+  for(i=0;i<q_nrow;i++,q_elem++)
+    q_new_col[i] -= *q_elem*tmp;
+}
+tmp = 0.0;
+for(i=0;i<q_nrow;i++) tmp += q_new_col[i]*q_new_col[i];
+norm_new = sqrt(tmp);
+
+  if( norm_new >= norm_last/2.0) break;
+  if( norm_new > 0.1*norm_init*PRECISION )
+    norm_last = norm_new;
+  else{
+    norm_init = norm_last = 0.1*norm_last*PRECISION;
+    for(i=0;i<q_nrow;i++)
+      q_new_col[i] = 0.0;
+    if(q_use_row == q_nrow)
+      ERRMSG("qr_add","Cannot orthogonalise new column\n");
+
+    q_new_col[q_use_row] = norm_new = norm_last;
+    q_use_row++;
+   }
+  }
+  
+  for(i=0;i<q_nrow;i++)
+    q_new_col[i] /= norm_new;
+  for(i=0;i<r_ncol;i++)
+    r_new_col[i] *= norm_orig;
+  r_new_col[r_ncol] = norm_new*norm_orig;
+  if(swap){
+
+col_l = r_new_col; /* first element in last column in R */
+col_lm1 = r_new_col-r_ncol; /* first element in column before last in R */
+for(j=0;j<r_ncol;j++,col_l++,col_lm1++){
+  tmp = *col_l;
+  *col_l = *col_lm1;
+  *col_lm1 = tmp;
+}
+a = col_lm1-1;
+b = col_l;
+tau = fabs(*a)+fabs(*b);
+if( tau > 0.0 ){
+  /*Calculate Givens rotation*/
+  nu = tau*sqrt((*a/tau)*(*a/tau)+(*b/tau)*(*b/tau));
+  c = *a/nu;
+  s = *b/nu;
+
+  *a = nu;   /* Fix column before last in R */
+  *b = -s* *(b-1);  /* Fix last element of last column in R */
+  *(b-1) = c * *(b-1); /* Fix second last element of last column in R */
+  if( *b < 0){
+    *b = -*b;
+    
+  col_lm1 = QCOL(r_ncol-1);
+  col_l   = col_lm1+q_nrow;
+  for(j=0;j<q_nrow;j++){
+    tmp = c*col_lm1[j]+s*col_l[j];
+    col_l[j] = -c*col_l[j]+s*col_lm1[j];
+    col_lm1[j] = tmp;
+  }
+  
+  }else{
+    
+  col_lm1 = QCOL(r_ncol-1);
+  col_l   = col_lm1+q_nrow;
+  for(j=0;j<q_nrow;j++){
+    tmp = c*col_lm1[j]+s*col_l[j];
+    col_l[j] = c*col_l[j]-s*col_lm1[j];
+    col_lm1[j] = tmp;
+  }
+  }
+}
+
+  }
+  r_ncol++;
+
+}
+#if defined (S_Plus)
+static void errmsg(char *string){
+  PROBLEM "%s\n", string RECOVER(NULL_ENTRY);
+}
+#elif defined(Matlab)
+static void errmsg(char *where, char *string){
+  char str[1024];
+  sprintf(str, "Error in %s: %s\n", where, string);
+  mexErrMsgTxt(str);
+}
+#else
+static void errmsg(char *where, char *string){
+  fprintf(stderr, "Error in %s: %s\n", where, string);
+  exit(EXIT_FAILURE);
+}
+#endif
+
+int assertFail(const char *ex, const char *file, const char *func, const int line)
+{
+static bool ignore = false;
+
+    if (!ignore) {
+        char str[1024];
+        sprintf(str, "%s:%u: %s(): failed assertion `%s'\n",
+                     file, line, func, ex);
+        mexErrMsgTxt(str);
+    } /* if */
+    return(1);
+}
+
+static void fortifyPrintf(const char *s)
+{
+   mexPrintf(s);
+}
+
+void
+mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
+{
+   mxArray *array_ptr;
+   int i;
+   long dataSize;
+   long dimensions;
+   long succeed;
+   long verbose = 0;
+   long as_sub = 0;
+   double lagrangian;
+   double threshold;
+   double *mx;
+   double *y;
+   double *yhat1;
+   double *r;
+   double *beta;
+   int arg;
+
+#ifdef FORTIFY
+   Fortify_SetOutputFunc(fortifyPrintf);
+   Fortify_EnterScope();
+#endif
+
+   if (nrhs > 3) {
+      mexErrMsgTxt("Syntax error: too many arguments\n");
+      return;
+   } /* if */
+
+   if (nrhs < 3) {
+      mexErrMsgTxt("Syntax error: arguments missing\n");
+      return;
+   } /* if */
+
+   if (nlhs > 1) {
+      mexErrMsgTxt("Syntax error: too many return arguments\n");
+      return;
+   } /* if */
+
+   arg = 2;
+   array_ptr = (mxArray *) prhs[arg];
+   if (mxGetClassID(array_ptr) != mxDOUBLE_CLASS) {
+      mexErrMsgTxt("Threshold argument is not a number\n");
+      return;
+   } /* if */
+
+   threshold = mxGetScalar(array_ptr);
+
+   arg = 1;
+   array_ptr = (mxArray *) prhs[arg];
+   if (mxGetClassID(array_ptr) != mxDOUBLE_CLASS) {
+      mexErrMsgTxt("Label argument is not a number\n");
+      return;
+   } /* if */
+
+   y = mxGetPr(array_ptr);
+
+   arg = 0;
+   array_ptr = (mxArray *) prhs[arg];
+   if (mxGetClassID(array_ptr) != mxDOUBLE_CLASS) {
+      mexErrMsgTxt("Data set argument is not a double matrix\n");
+      return;
+   } /* if */
+
+   dataSize = mxGetM(array_ptr);
+   dimensions = mxGetN(array_ptr);
+
+   mx = mxGetPr(array_ptr);
+
+   arg = 0;
+   plhs[arg] = mxCreateDoubleMatrix(dimensions, 1, mxREAL);
+
+   yhat1 = new double[dataSize];
+   r = new double[dataSize];
+   beta = mxGetPr(plhs[arg]);
+
+   for (i=0; i<dimensions; i++) {
+       beta[i] = 0.0;
+   } /* for */
+
+   if (dimensions > 0) {
+       lasso(mx, &dataSize, &dimensions, &threshold, beta, y,
+             yhat1, r, &lagrangian, &succeed, &verbose, &as_sub);
+   } else {
+       mexPrintf("Dataset with 0 dimensions\n");
+   } /* if */
+
+   delete [] yhat1;
+   delete [] r;
+
+#ifdef  FORTIFY
+   Fortify_LeaveScope();
+#endif
+
+}
+
Index: /prextra/prextra/mk_lasso.m
===================================================================
--- /prextra/prextra/mk_lasso.m	(revision 2)
+++ /prextra/prextra/mk_lasso.m	(revision 2)
@@ -0,0 +1,4 @@
+clear functions
+
+ mex -O -v -output ./lasso mexlasso.cc
+
Index: /prextra/prextra/neighbh.m
===================================================================
--- /prextra/prextra/neighbh.m	(revision 2)
+++ /prextra/prextra/neighbh.m	(revision 2)
@@ -0,0 +1,186 @@
+function J = neighbh(a,nbtype)
+% J = neighbh(sz,nbtype)
+% J = neighbh(a,nbtype)
+% 
+% Compute the indices of the pixels that are in the neighborhood of a
+% pixel. This is done for all pixels in a whole image and returned in
+% the index array J. The pixels are enumerated in this standard order
+% (for an 3x5 image):
+%     1   4   7  10  13
+%     2   5   8  11  14
+%     3   6   9  12  15
+% When there is no neighbor (i.e. the pixel is on the boundary), index 0
+% is returned.
+% Then J becomes for a '4NN' neighborhood:
+%     0  0  2  4      % pixel 1 have neighbors 2 and 4...
+%     0  1  3  5      % pixel 3 has 3 neighbors..
+%     0  2  0  6
+%     1  0  5  7
+%     2  4  6  8
+%     and so on ...
+%
+% We can define neighborhood types nbtype:
+%   '4NN'   4-nearest neighborhood
+%   '8NN'   8-nearest neighborhood
+%   '9NN'   9-nearest neighborhood
+%  '25NN'  25-nearest neighborhood
+% It is also possible to give your own definition of the neighborhood.
+% You should give the indices w.r.t. the origin, for example:
+%   nbtype = [-3 0;
+%             +3 0;
+%              0 -1];
+%
+
+if isdataset(a)
+	if ~chckcoord(a)
+		error('Does not look like an image dataset');
+	end
+	sz = getobjsize(a);
+else
+	if ~(size(a)==[1 2])
+		error('I want to have the size of the image');
+	end
+	sz = a;
+end
+% The total number of pixels in the image:
+N = prod(sz);
+% The definitions of the predefined neighborhoods:
+if ischar(nbtype)
+	switch nbtype
+	case '4NN'
+		% left    up     down    right
+		pos = [-1 0; 0 -1; 0  1; 1 0];
+	case '8NN'
+		pos = [-1 -1; -1 0; -1 1; 0 -1; 0 1; 1 -1; 1 0; 1 1];
+	case '9NN'
+		pos = [-1 -1; -1 0; -1 1; 0 -1; 0 0; 0 1; 1 -1; 1 0; 1 1];
+	case '25NN'
+		pos = [-2 -2; -2 -1; -2  0; -2  1; -2  2; -1 -2; -1 -1;...
+			 -1  0; -1  1; -1  2; 0 -2; 0 -1; 0  0; 0  1; 0  2; ...
+			  1 -2; 1 -1; 1  0; 1  1; 1  2; 2 -2; 2 -1; 2  0;...
+			  2  1; 2  2];
+	case '4star'
+		%      left   up   down right
+		pos = [-2 0; 0 -2; 0  2; 2 0];
+	case '4cross'
+		%      upleft  upright downleft downright
+		pos = [-2 -2;   2 -2;   -2  2;   2 2];
+	otherwise
+		error('This neighborhood is not defined');
+	end
+else
+	if size(nbtype,2)~=2
+		error('The neighborhood definition should be a nx2 matrix.');
+	end
+	pos = nbtype;
+end
+% convert these neighborhood positions into changes in indices:
+dx = sz(1)*pos(:,1)' + pos(:,2)';
+% in horizontal limits (lower and upper):
+ckhorz = -pos(:,1)';
+I = find(ckhorz<0);
+ckhorz(I)=ckhorz(I)+sz(2)+1;
+% and in vertical limits (lower and upper):
+ckvert = -pos(:,2)';
+I = find(ckvert<0);
+ckvert(I)=ckvert(I)+sz(1)+1;
+
+% Setup the general parameters
+D = length(dx);
+
+% Define the pixels, and compute their coordinates:
+I = (1:N)';
+hcoord = floor((I-1)/sz(1))+1;
+vcoord = mod(I-1,sz(1))+1;
+
+% here the results will be stored:
+J = zeros(N,D);
+
+% and here we go:
+for i=1:D
+	% Here the main work of constructing the J matrix is performed,
+	% it is actually just these two lines. The trick is to shift the I
+	% vector up or down and store it alongside the already created matrix
+	% J.
+
+	if dx(i)>0   % look at pixels in the + direction
+
+		J(:,i) = [I((1+dx(i)):end); zeros(dx(i),1)];
+
+	else         % look at pixels in the - direction
+
+		J(:,i) = [zeros(-dx(i),1); I(1:end+dx(i))];
+
+	end
+	% Unfortunately, most of the work is actually in the checking of the
+	% bounds:
+	% Check for the horizontal bounds:
+	% look at the pixels left:
+	if pos(i,1)<0
+		if ckhorz(i)
+			excp = find(hcoord<=ckhorz(i));
+			J(excp,i) = 0;
+		end
+	else  % and pixels right
+		if ckhorz(i)
+			excp = find(hcoord>=ckhorz(i));
+			J(excp,i) = 0;
+		end
+	end
+	% Check for the vertical bounds:
+	% look at the pixels *above*:
+	if pos(i,2)<0
+		if ckvert(i) 
+			excp = find(vcoord<=ckvert(i));
+			J(excp,i) = 0;
+		end
+	else   % otherwise the pixels below:
+		if ckvert(i)
+			excp = find(vcoord>=ckvert(i));
+			J(excp,i) = 0;
+		end
+	end
+end
+
+return
+
+
+function out = chckcoord(a)
+% out = chckcoord(a)
+%
+% check if dataset a contains pixels regularly sampled from an image
+% It looks a the objsize and maybe the X- and Y-coordinates.
+
+out = 1;
+
+sz = getobjsize(a);
+
+% first check if we even have some 2D objsize
+if length(sz)==1
+	out = 0;
+	return
+end
+
+featlab = getfeatlab(a);
+% now check the possible X coordinate:
+xcoord = strmatch('X',featlab);
+if ~isempty(xcoord)
+	cx = +a(:,xcoord);
+	dx = diff(cx(1:sz(1)));
+	if std(dx)~=0
+		out = 0;
+		return
+	end
+end
+% now check the possible Y coordinate:
+ycoord = strmatch('Y',featlab);
+if ~isempty(ycoord)
+	cy = +a(:,ycoord);
+	dy = diff(cy(1:sz(1):end));
+	if std(dy)~=0
+		out = 0;
+		return
+	end
+end
+
+return
Index: /prextra/prextra/rbsvc.m
===================================================================
--- /prextra/prextra/rbsvc.m	(revision 2)
+++ /prextra/prextra/rbsvc.m	(revision 2)
@@ -0,0 +1,43 @@
+%RBSVC Automaric radial basis SVM
+%
+%   W = RBSVC(A,FID)
+
+function w = rbsvc(a,fid)
+
+if nargin < 2, fid = []; end
+if nargin < 1 | isempty(a)
+	w = mapping(mfilename,fid);
+	return
+end
+
+nu = max(testk(a,1),0.05);
+
+d = sqrt(+distm(a));
+sigmax = min(max(d)); % max: smallest furthest neighboe distance
+
+d = d + 1e100*eye(size(a,1));
+sigmin = max(min(d)); % min: largest nearest neighbor distance
+
+q = (1/sigmin - 1/sigmax)/9;
+sig = [1/sigmax:q:1/sigmin];
+sig = ones(1,length(sig))./[1/sigmax:q:1/sigmin];
+if isempty(sig)
+	sigopt = sigmin;
+else
+	w = [];
+	errmin = inf;
+	sigopt = 0;
+	prprogress(fid,'     ')
+	for j=1:length(sig)
+		s = sprintf('\b\b\b\b\b%5.0f',j);
+		prprogress(fid,s);
+		err = crossval(a,svc_nu([],'r',sig(j),nu),10,1,fid);
+		if err < errmin
+			errmin = err;
+			sigopt = sig(j);
+		end
+	end
+	prprogress(fid,'\b\b\b\b\b')
+end
+w = svc_nu(a,'r',sigopt,nu);
+
Index: /prextra/prextra/removeilab.m
===================================================================
--- /prextra/prextra/removeilab.m	(revision 2)
+++ /prextra/prextra/removeilab.m	(revision 2)
@@ -0,0 +1,73 @@
+%REMOVEILAB remove named identifier columns
+%
+%   B=REMOVEILAB(A,REF)
+%
+% INPUT
+%    A     input dataset with ilab structure (see MAKEILAB)
+%    REF   name, cell array with names or numerical indices of
+%          identifier columns
+% OUTPUT
+%    B     output dataset
+%
+% DESCRIPTION 
+% This cunction removes named identifiers (both names and per-example
+% values) from a ilab-enabled dataset (see ENABLEILAB).  If no identifier
+% name or index is given (REF is empty), the named-identifier capability is
+% removed from the dataset.
+%
+% SEE ALSO
+% ENABLEILAB, REMOVEILAB, GETILAB, SETILAB, GETILABLIST
+
+% $Id: removeilab.m,v 1.2 2005/05/02 13:05:46 pavel Exp $
+
+function a=removeilab(a,ref)
+    
+    isvalidilab(a);
+
+    if ~exist('ref','var')
+	u=a.user;
+	u=rmfield(u,'ilab');
+	a.user=u;
+	prwarning(1,'ilab capability was removed from a dataset');
+	return
+    end
+    
+    u=a.user;
+    
+    % get the columns to be set
+    col=[];
+    if isnumeric(ref)
+	col=ref;    
+    else
+	if ischar(ref)
+	    ref={ref}; % turn it into cell array
+	end
+	if ~iscell(ref)
+	    error('numerical index, name or cell array of names expected')
+	end
+	% we know ref is a cell array
+	for i=1:length(ref)
+	    ind=strcmp(a.user.ilab.names,ref{i});
+	    if sum(ind)==0
+		error(['identifier label ''' ref{i} ''' not defined!']); 
+	    end
+	    col=[col find(ind==1)];
+	end
+    end
+    
+    if max(col)>size(a.user.ilab.names,2)
+	error('index exceeds number of defined ident labels');
+    end
+
+    names=u.ilab.names;
+    ncol=sort(setdiff(1:length(a.user.ilab.names),col));
+    names=names(ncol);
+    u.ilab.names=names;
+    a.user=u;
+    
+    % update ident
+    ident=a.ident;
+    ident(:,col)=[];
+    a.ident=ident;
+    
+    return
Index: /prextra/prextra/removemeta.m
===================================================================
--- /prextra/prextra/removemeta.m	(revision 2)
+++ /prextra/prextra/removemeta.m	(revision 2)
@@ -0,0 +1,57 @@
+%REMOVEMETA remove key-value pairs
+%
+%   B=REMOVEMETA(A,KEYS)
+%
+% INPUT
+%    A     input dataset with meta structure (see ENABLEMETA)
+%    KEYS  key (string) or a cell array with keys
+%
+% OUTPUT
+%    B     output dataset
+%
+% DESCRIPTION
+% Remove the key-value pairs specified by the keys. If no keys are
+% specified, the meta facility is removed from the dataset.
+%
+% SEE ALSO
+% REMOVEMETA, GETMETA, SETMETA
+
+% $Id: removemeta.m,v 1.3 2005/05/02 14:12:10 pavel Exp $
+
+function a=removemeta(a,inkeys)
+    
+    isvalidmeta(a);
+
+    u=a.user;    
+    if ~exist('inkeys','var')
+	u=rmfield(u,'meta');
+	a.user=u;
+	prwarning(1,'meta capability was removed from a dataset');
+	return
+    end
+	
+    keys=u.meta.keys;
+    values=u.meta.values;    
+    
+    if ischar(inkeys)
+	inkeys={inkeys}; % turn it into cell array
+    end
+    if ~iscell(inkeys)
+	error('name or a cell array of names expected as keys')
+    end
+    % we know inkeys is a cell array
+    for i=1:length(inkeys)
+	ind=strcmp(a.user.meta.keys,inkeys{i});
+	if sum(ind)==0
+	    error(['key ''' inkeys{i} ''' not defined!']); 
+	end
+	ind=find(ind);
+	keys(ind)=[];
+	values(ind)=[];
+    end
+    
+    u.meta.keys=keys;
+    u.meta.values=values;
+    a.user=u;
+    
+    return
Index: /prextra/prextra/setilab.m
===================================================================
--- /prextra/prextra/setilab.m	(revision 2)
+++ /prextra/prextra/setilab.m	(revision 2)
@@ -0,0 +1,73 @@
+%SETILAB set the values for the per-example named-identifiers in a dataset
+%
+%   B=SETILAB(A,REF,VALUES,J)
+%
+% INPUT
+%    A       input dataset with ilab structure (see MAKEILAB)
+%    REF     name, cell array with names or numerical indices of
+%            identifier columns
+%    VALUES  identifier values to be set. If all values are numerical, 
+%            matrix may be supplied, otherwise a cell array.
+%    J       sample indices to be retrieved (optional)
+%
+% OUTPUT
+%    B       output dataset
+%
+% DESCRIPTION
+% SETILAB sets the values of identifiers, specified by names for all
+% examples in a dataset (or for a subset, given by example indices J).
+% VALUES may be a matrix, in case all values to be set are numerical, or
+% a cell array otherwise.
+%
+% SEE ALSO
+% ENABLEILAB, REMOVEILAB, GETILAB, SETILAB, GETILABLIST
+
+function a=setilab(a,ref,labs,J)
+    
+    isvalidilab(a);
+
+    if isempty(ref)
+	% replace all identifiers
+	ref=1:length(a.user.ilab.names);
+    end
+    
+    if ~exist('J','var')
+	J=1:size(a,1);
+    end
+    
+    % get the columns to be set
+    col=[];
+    if isnumeric(ref)
+	col=ref;    
+    else
+	if ischar(ref)
+	    ref={ref}; % turn it into cell array
+	end
+	if ~iscell(ref)
+	    error('numerical index, name or cell array of names expected')
+	end
+	% we know ref is a cell array
+	for i=1:length(ref)
+	    ind=strcmp(a.user.ilab.names,ref{i});
+	    if sum(ind)==0
+		error(['identifier label ''' ref{i} ''' not defined!']); 
+	    end
+	    col=[col find(ind==1)];
+	end
+    end
+    
+    if max(col)>size(a.user.ilab.names,2)
+	error('index exceeds number of defined ident labels');
+    end
+ 
+    if isnumeric(labs)
+	labs=num2cell(labs);
+    end
+    
+    % retrieve ident labels
+    temp=a.ident;
+    temp(J,col)=labs;
+    a.ident=temp;
+    
+    return
+    
Index: /prextra/prextra/setmeta.m
===================================================================
--- /prextra/prextra/setmeta.m	(revision 2)
+++ /prextra/prextra/setmeta.m	(revision 2)
@@ -0,0 +1,52 @@
+%SETMETA add/update key-value pair(s) to a dataset
+%
+%   B=SETMETA(A,KEYS,VALUES)
+%
+% INPUT
+%    A       input dataset with meta structure (see MAKEMETA)
+%    KEYS    string key or a cell-array with key names
+%    VALUES  key-specific values. If more keys are defined, VALUES
+%            must be stored in a cell-array
+%
+% OUTPUT
+%    B     output dataset
+%
+% DESCRIPTION
+% This function sets the key-value pair(s) in a user field of a dataset.
+%
+% SEE ALSO
+% REMOVEMETA, GETMETA, SETMETA
+
+function a=setmeta(a,inkeys,invalues)
+    
+    isvalidmeta(a);
+
+    u=a.user;
+
+    if ischar(inkeys), inkeys={inkeys}; end
+    if ~iscell(invalues), invalues={invalues}; end    
+
+    if ~iscell(inkeys) | ~iscell(invalues)
+	error('unsupported format of keys or values')
+    end
+    if length(inkeys)~=length(invalues)
+	error('keys and values to be set must have equal sizes');
+    end
+
+    % update meta in user
+    u=a.user;
+    for j=1:length(inkeys)
+	% is the key already defined?
+	ind=strcmp(u.meta.keys,inkeys{j});
+	if any(ind)
+	    k=find(ind);
+	    u.meta.keys{k}=inkeys{j};
+	    u.meta.values{k}=invalues{j};
+	else % not found, add a new key-value
+	    u.meta.keys{end+1}=inkeys{j};
+	    u.meta.values{end+1}=invalues{j};
+	end	
+    end
+    a.user=u;
+    
+    return
Index: /prextra/prextra/showfigs.m
===================================================================
--- /prextra/prextra/showfigs.m	(revision 2)
+++ /prextra/prextra/showfigs.m	(revision 2)
@@ -0,0 +1,19 @@
+%SHOWFIGS Show all figures on the screen
+
+function showfigs
+h = sort(get(0,'children'));
+n = length(h);
+k = ceil(sqrt(n));
+s = 0.9/k;   % screen stitch
+r = 0.96;     % image size reduction
+t = 0.055;    % top gap
+b = 0.005;    % border gap
+fig = 0;
+for i=1:k
+	for j=1:k
+		fig = fig+1;
+		if fig > n, break; end
+		set(h(fig),'units','normalized','position',[(j-1)*s+b,(1-t)-i*s,s*r,s*r]);
+		figure(h(fig));
+	end
+end
Index: /prextra/prextra/testmc.m
===================================================================
--- /prextra/prextra/testmc.m	(revision 2)
+++ /prextra/prextra/testmc.m	(revision 2)
@@ -0,0 +1,52 @@
+%TESTMC More classifier tests
+%
+% E = TESTMC(A,W,TYPE)
+% E = TESTMC(A*W,TYPE)
+
+function e = testmc(a,w,varargin)
+
+prtrace(mfilename);
+
+if nargin < 1
+	error('Insufficient parameters given')
+end
+
+isdataset(a);
+
+if nargin == 1
+	d = a;
+	type = {'conf'};
+elseif ismapping(w)
+	d = a*w;
+	if nargin < 3
+		type = {'conf'};
+	else
+		type = varargin;
+	end
+elseif isstr(w)
+	d = a;
+	type = {w};
+elseif iscell(w)
+	d = a;
+	type = w;
+else
+	error('Illegal call')
+end
+
+d = d*maxc;
+I = matchlablist(getlablist(d),getfeatlab(d));
+J = I(d.nlab); % J(i) is the column of the true label of object i
+e = zeros(1,length(type));
+m = size(d,1);
+for j=1:length(type)
+	switch(type{j})
+		case 'conf'
+            d = +normm(d,1);
+			dd = d([1:m]'+(J-1)*m);
+			e(j) = 1-mean(dd);
+		otherwise
+			error('Illegal type')
+	end
+end
+			
+			
Index: /prextra/prextra/unitm.m
===================================================================
--- /prextra/prextra/unitm.m	(revision 2)
+++ /prextra/prextra/unitm.m	(revision 2)
@@ -0,0 +1,36 @@
+%UNITM Fixed unit mapping
+% 
+%   W = UNITM
+%   B = A*UNITM
+%   B = UNITM(A)
+%
+% INPUT
+%   A   Array or dataset
+%
+% OUTPUT
+%   W   Unit mapping, if applied to dataset, it is returned unchanged
+%   B   A
+%
+% DESCRIPTION
+% This is a fixed unit mapping that maps any dataset on itself. There is
+% also a trainable unit mapping, named UNITTM.
+%
+% SEE ALSO
+% MAPPINGS, UNITTM
+
+% 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 = unitm (a,v)
+
+	prtrace(mfilename);
+
+	if (nargin == 0) | (isempty(a))
+		w = mapping(mfilename,'fixed');
+		w = setname(w,'Unit Mapping');
+	else
+		w = a;
+	end
+
+	return
Index: /prextra/prextra/unittm.m
===================================================================
--- /prextra/prextra/unittm.m	(revision 2)
+++ /prextra/prextra/unittm.m	(revision 2)
@@ -0,0 +1,43 @@
+%UNITTM Trainable unit mapping
+% 
+%   W = UNITTM(A)
+%   W = A*UNITTM
+%
+% INPUT
+%   A   Array or dataset
+%
+% OUTPUT
+%   W   Unit mapping, if applied to dataset, it is returned unchanged
+%
+% DESCRIPTION
+% This is a trainable unit mapping that maps any dataset on itself. 
+% The difference with the fixed mapping is UNITM is that UNITM(A) returns
+% A and that W = UNITTM(A) returns a mapping with input and output 
+% dimensionality the feature size of A. Consequently B*W return B, but
+% checks the feature size against A.
+%
+% SEE ALSO
+% MAPPINGS, UNITM
+
+% 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 = unittm (a,v)
+
+	prtrace(mfilename);
+
+	if (nargin == 0) | (isempty(a))     % untrained definition
+		w = mapping(mfilename,'untrained');
+		w = setname(w,'Unit Mapping');
+	elseif (nargin == 1) | isempty(v)   % training
+		isdataset(a);
+		[m,k] = size(a);
+		w = mapping(mfilename,'trained',[],getfeatlab(a),k,k);
+	elseif ismapping(v)                 % execution
+		w = a;                            % return input  
+	else
+		error('Un-pars-able call')
+	end
+
+	return
Index: /prextra/prextra/wvotec.m
===================================================================
--- /prextra/prextra/wvotec.m	(revision 2)
+++ /prextra/prextra/wvotec.m	(revision 2)
@@ -0,0 +1,112 @@
+%WVOTEC Weighted combiner (Adaboost weights)
+%
+%  W = WVOTEC(A,V)   compute weigths and store
+%  W = WVOTEC(V,U)   Construct weighted combiner using weights U
+%
+%  INPUT
+%    A      Labeled dataset
+%    V      Parallel or stacked set of trained classifiers
+%    U      Set of classifier weights
+%
+%  OUTPUT
+%    W      Combined classifier
+%
+% DESCRIPTION
+% The set of trained classifiers V is combined using weighted
+% majority voting. If given the weights U are used. If not
+% given, the weights are computed from the classification
+% results of the labeled dataset A using 0.5*log((1-E)/E)
+% if E is the classifier error. 
+%
+% SEE ALSO
+% MAPPINGS, DATASETS,
+
+% 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 = wvotec(a,v)
+
+prtrace(mfilename);
+
+if nargin < 1 | isempty(a)
+	w = mapping(mfilename,'untrained');
+elseif nargin < 2
+	error('Illegal call')
+elseif isdataset(a)                    % train or classify
+	if ~strcmp(v.mapping_file,mfilename) % training
+		if isparallel(v)                   % parallel combiner
+			n = 0;
+			e = zeros(1,length(v.data));
+			for j=1:length(v.data)
+				vv = v.data{j};
+				d = a(:,n+1:n+size(vv,1))*vv*classc;
+				e(j) = testmc(d);
+				n = n+size(vv,1);
+			end
+		elseif isstacked(v)                % stacked combiner
+			e = zeros(1,length(v.data));
+			for j=1:length(v.data)
+				vv = v.data{j};
+				e(j) = testmc(a,vv);
+			end
+		else
+			error('Classifier combination expected')
+		end
+		                  % Find weights								
+		L = find(e < 1-max(getprior(a))); % take classifier better than prior
+		alf = zeros(1,length(e));
+		alf(L) = log((1-e(L))./e(L))/2;
+		alf = alf/sum(alf);
+		                  % construct the classifier
+		[m,k,c] = getsize(a);
+		w = mapping(mfilename,'trained',{v,alf},getlabels(vv),k,c);
+		w = setname(w,'Weighted Voting');
+	else                                 % testing
+		alf = v.data{2};                   % get the weights
+		u = v.data{1};                     % get the set of classifiers
+		m = size(a,1);
+		dtot = zeros(m,size(v,2));
+		if isparallel(u)                   % parallel combiner
+			n = 0;
+			for j=1:length(u.data)           % weight them 
+				vv = u.data{j};
+				aa = a(:,n+1:n+size(vv,1));
+				d = a(:,n+1:n+size(vv,1))*vv;
+				[dd,jj] = max(+d,[],2);
+				dd = zeros(size(dtot));
+				dd([1:m]'+(jj-1)*m) = alf(j);
+				dtot = dtot + dd;
+				n = n+size(vv,1);
+			end
+		elseif isstacked(u)                % stacked combiner
+			for j=1:length(u.data)           % weight them
+				vv = u.data{j};
+				d = a*vv;
+				[dd,jj] = max(+d,[],2);
+				dd = zeros(size(dtot));
+				dd([1:m]'+(jj-1)*m) = alf(j);
+				dtot = dtot + dd;
+			end
+		else
+			error('Classifier combination expected')
+		end
+		w = setdat(d,dtot);
+  end
+  
+else                  % store classifier from given weights
+    
+  ismapping(a);
+  u = v;              % the weights
+  v = a;              % the combined classifier
+  n = length(v.data);
+  if length(u) ~= n
+    error('Wrong number of weights given')
+  end
+  [k,c] = getsize(v.data{1});
+	w = mapping(mfilename,'trained',{v,u},getlabels(v{1}),k,c);
+	w = setname(w,'Weighted Voting');
+end
+
+		
+			
