 function [w,beta,residp,residn,alp,totalviolation,dualgap,flag]...
    = sep(Xp,Xn,penalty);

 % This subroutine is designed to calculate a linear discriminator
 % based on the training data Xp and Xn and the violation cost penalty.
 % It uses the DWD method of J.S. Marron and M.J. Todd,
 % "Distance Weighted Discrimination," TR 1339,
 % School of Operations Research and Industrial Engineering,
 % Cornell University, Ithaca, New York, available at
 % ftp://ftp.orie.cornell.edu/pub/techreps/TR1339.pdf.

 % The user will need to get the SDPT3 optimization package, 
 % available at http://www.math.nus.edu.sg/~mattohkc/sdpt3.html, 
 % and install it.

 % Input:  Matrices Xp and Xn, whose columns contain the training
 %            data for positive and negative instances respectively,
 %         and penalty, the cost for perturbing each residual.

 % Output: w and beta, such that x with w'x + beta positive (resp.,
 %            negative) is classified in the "+1" (resp. "-1") class.
 %         residp, the values of w'x + beta for x a column of Xp,
 %         residp, the values of w'x + beta for x a column of Xp.
 %         alp, the dual solution (usually, w = [Xp, - Xn]*alp
 %          normalized).
 %         totalviolation, the total amount added to the residuals.
 %         dualgap, the duality gap, a measure of the accuracy
 %            of the solutions found.
 %         flag, an indication of the success of the computation:
 %             0, success;
 %            -1, inaccurate solution;
 %            -2, problem infeasible or unbounded.

 %stime = cputime;
 flag = 0;

 % Find the dimensions of the data.

 [dp,np] = size(Xp);
 [dn,nn] = size(Xn);
 if (dn ~= dp), error('The dimensions are incompatible.'), end;
 d = dp;

 % Do the dimension reduction if in HDLSS setting.

 XpnY = [Xp, -Xn];
 n = np + nn;
 if (d > n),
    [Q,RpnY] = qr(XpnY,0);
    %qrtime = cputime - stime, stime = cputime;
    dnew = n;
 else,
    RpnY = XpnY;
    dnew = d;
 end;
 y = [ones(np,1); -ones(nn,1)];
 ym = y(2:n);

 % nv is the number of variables (eliminating beta),
 % nc the number of constraints.

 nv = 1 + dnew + 3*n + n;
 nc = 2*n;
 %nv,
 %nc,

 % Set up the block structure, constraint matrix, rhs, and cost vector.

 blk = cell(2,2);
 blk{1,1} = 'q';
 blk{1,2} = [dnew+1, 3*ones(1, n)];
 blk{2,1} = 'l';
 blk{2,2} = n;
 %blk,

 Avec = cell(2,1);
 A = zeros(nc,nv-n);
 col1 = RpnY(:,1);
 A(1:n-1,2:dnew+1) = (RpnY(:,2:n) - col1*ym')';
 A(1:n-1,dnew+5:3:dnew+1+3*n) = - speye(n-1,n-1);
 A(1:n-1,dnew+6:3:dnew+2+3*n) = + speye(n-1,n-1);
 A(1:n-1,dnew+2) = ym;
 A(1:n-1,dnew+3) = -ym;
 A(n,1) = 1;
 A(n+1:n+n,dnew+4:3:dnew+3+3*n) = speye(n,n);
 %A,
 Avec{1,1} = A';
 Avec{2,1} = [[-ym,speye(n-1,n-1)];zeros(1+n,n)]';

 b = [zeros(n-1,1);ones(1+n,1)];
 %b,

 C = cell(2,1);
 c = zeros(nv-n,1);
 c(dnew+2:3:dnew+1+3*n) = ones(n,1);
 c(dnew+3:3:dnew+2+3*n) = ones(n,1);
 %c,
 C{1,1} = c;
 C{2,1} = penalty*ones(n,1);
 %setuptime = cputime - stime, stime = cputime;

 % Solve the SOCP problem.

 parameters;
 %startup;
   global CACHE_SIZE   % cache size in kbytes 
   global LOOP_LEVEL   % loop unrolling level
   CACHE_SIZE = 256; 
   LOOP_LEVEL = 8;
   %  instead of messing things up by calling startup,
   %  better to use these lines instead
 OPTIONS.vers = 2;
 OPTIONS.maxit = 40;
 %OPTIONS.gaptol = 1e-10;
 [X0,lambda0,Z0] = infeaspt(blk,Avec,C,b);
 %startpttime = cputime - stime,
 [obj,X,lambda,Z,gaphist,infeashist,info] = ...
    sqlp(blk,Avec,C,b,X0,lambda0,Z0,OPTIONS);

 % If infeasible or unbounded, break.

 if (info(1) > 0), flag = -2; break, end;

 % Compute the normal vector w and constant term beta.

 X1 = X{1}; X2 = X{2};
 barw = X1(2:dnew+1);
 if (d>n),
    w = Q*barw;
 else,
    w = barw;
 end;
 %w17 = w(1:7)',
 beta = X1(dnew+2) - X1(dnew+3) - X2(1) - col1'*barw,
 normw = norm(w);
 if normw < 1 - 1e-3, normw, end;
 normwm1 = 0;
 if normw > 1 - 1e-3,
    w = w / normw;
    normwm1 = norm(w)-1;
    beta = beta / normw;
 end;

 % Compute the minimum of the supposedly positive
 % and the maximum of the supposedly negative residuals.
 % Refine the primal solution and print its objective value.

 residp = Xp'*w + beta*ones(np,1);
 residn = Xn'*w + beta*ones(nn,1);
 minresidp = min(residp), maxresidn = max(residn),
 res = XpnY'*w + beta*y;
 rsc = 1/sqrt(penalty);
 xi = rsc - res;
 xi = max(xi,0);
 totalviolation = sum(xi),
 minresidpmod = min(residp+xi(1:np));
 maxresidnmod = max(residn-xi(np+1:n));
 minxi = min(xi);
 maxxi = max(xi);
 resn = res + xi;
 rresn = 1./resn;
 format long
 primalobj = penalty*sum(xi) + sum(rresn),
 %devprimalobj = primalobj - obj(1);
 %if abs(devprimalobj) > 1e-3,
 %   devprimalobj,
 %   objxioff = penalty*(sum(xi)-sum(X2)),
 %end;

 % Compute the dual solution alp and print its objective value.

 alp = zeros(n,1);
 lambda1 = lambda(1:n-1);
 alp(1) = -ym'*lambda1;
 alp(2:n) = lambda1;
 alp = max(alp,0);
 sump = sum(alp(1:np));
 sumn = sum(alp(np+1:n));
 sum2 = (sump + sumn)/2;
 alp(1:np) = (sum2/sump)*alp(1:np);
 alp(np+1:n) = (sum2/sumn) * alp(np+1:n);
 maxalp = max(alp);
 if (maxalp > penalty | maxxi > 1e-3),
    alp = (penalty/maxalp) * alp;
 end;
 minalp = min(alp);
 p = RpnY*alp;
 eta = - norm(p);
 gamma = 2 * sqrt(alp);
 dualobj = eta + sum(gamma),
 %devdualobj = dualobj - obj(2);
 %if abs(devdualobj) > 1e-3,
 %   devdualobj,
 %end;
 format short

 % dualgap is the duality gap, a measure of the accuracy of the solution.

 dualgap = primalobj - dualobj,
 %if normw > 1 - 1e-3,
 %   wfromdual = -(XpnY*alp)/eta;
 %   normdifwprimwdual = norm(w-wfromdual);
 %   wd17 = wfromdual(1:7)';
 %end;
 if (dualgap > 1e-4),
    flag = -1;
 end;
 %arccoswe1 = w(1);
 %if (d>n),
 %   e1 = zeros(d,1); e1(1) = 1;
 %   pe1 = Q*(Q'*e1);
 %   npe1 = sqrt(pe1(1));
 %   arccoswpe1 = w'*pe1 / npe1,
 %   arccose1pe1 = npe1;
 %end;
 return

