function [w,theta] = msmt_msmtlp(A,B)
% MSMT_MSMTlp solves linear program used in MSM-T
%
%	[w,theta] = msmt_msmtlp(A,B) solves the following lp:
%
%	For matrices A,B the following lp is solved:
%		
%	min		(1/m)ey + (1/k)ez
%	w,theta,y,z
%
%		s.t.	Aw - e(theta) + y >= e
%		       -Bw + e(theta) + z >= e
%					y >= 0
%					z >= 0
%				
%				w,theta: unrestricted
%
%	[w,theta,y,z] are elements of the solution set 
%
%	A represents the m n-dimensional elements of the point set A.
%	B represents the k n-dimensional elements of the point set B.
%	A and B are disjoint.
%	The lp determines the "best" seperating plane, wx = theta, 
%	between the two sets.
%
%	NOTE: e is a vector of 1's of appropriate dimension
%
%	Solve the dual to the lp above is much more efficient.  By 
%	observing the transpose of the constraint matrix C, (which is
%	the constraint matrix in the dual), we can eliminate many 
%	constraints and put that information in the form of upper and
%	lower bounds on the dual variables.
%
%	primal:	min	c'x
%	          s.t.	C(i)x = b(i) 		for i in M
%			C(i)x >= b(i) 		for i in N
%			x(j) >= 0		for j in O
%			x(j):unrestricted	for j in P
%
%	dual:	max	b'u
%		  s.t.	u(i):unrestricted	for i in M
%			u(i) >= 0		for i in N
%			C'(j) >= c(j)		for j in O
%			C'(j) = c(j)		for j in P

% NOTE: in this problem:
%	x(1:n) = w
%	x(n+1) = theta
%	x(n+2:n+m+1) = y
%	x(n+m+k+1) = z

% obtain the m,k,n:
[m,n] = size(A);
k = size(B,1);

% the vector c will be all zeros:
% in the dual problem, this becomes the right-hand-side vector, in our
% case, this vector will be 0.  The information from the non-zero terms
% of c will be taken care of as upper and lower bounds on the dual 
% variable, u.
%c = zeros(n+1,1);

% construct constraint matrix C:
%C(1:m,1:n) = A;
%C(1:m,n+1) = -ones(m,1);
%C(m+1:m+k,1:n) = -B;
%C(m+1:m+k,n+1) = ones(k,1);

% construct vector b:
%b = ones(m+k,1);

% since this dual is a maximization problem, and function lp solves
% minimization problems, must take the negative of the dual coeffiecient

% construct vlb,vub:
%vlb = zeros(m+k,1);
%vub = zeros(m+k,1);
%vub(1:m)= (1/m)*ones(m,1);
%vub(m+1:m+k) = (1/k)*ones(k,1);


% construct x0:
%x0 = zeros(m+k,1);

% solve this lp
%[obj,u,lambda,status,cstat,it] = cplex(-b,C',c,vlb,vub);
%[x,lambda] = lp(-b,C',c,vlb,vub);

% obtain the primal solution x:
%w = lambda(1:n);
%theta = lambda(n+1);

%c = [ zeros(1,n+1), (1/m)*ones(1,m), (1/k)*ones(1,k)]';
%C = [ -A, ones(m,1), -eye(m), zeros(m,k);
%       B, -ones(k,1), zeros(k,m), -eye(k);
%       zeros(m,n), zeros(m,1), -eye(m), zeros(m,k);
%       zeros(k,n), zeros(k,1), zeros(k,m), -eye(k) ];

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

%x = lp(c,C,b)
%w = -x(1:n);
%theta = -x(n+1);

% dual
c = [ -ones(1,m+k) ]';
C = [ A', -B'; -ones(1,m), ones(1,k)];
b = zeros(n+1,1);
vlb = zeros(m+k,1);
vub = [ (1/m)*ones(m,1); (1/k)*ones(k,1)];
u0 = zeros(m+k,1);

[x,lambda] = lp(c,C,b,vlb,vub,u0,n+1);

w = lambda(1:n);
theta = lambda(n+1);

return












