% Solves: Maximize c^Tx +  d^Tw
%         s.t.     Ax + Dw = b
%                  B_1x       = d_1
%                  x >= 0, w >= 0
% using the Dantzig-Wolfe decomposition scheme
% The idea is to modify the revised simplex code "revised_simplex.m"
% on the course webpage to handle the block angular structure in the LP


% Kartik's MATLAB code for MA/OR/IE 505
% April 19, 2006

% Key in problem parameters

c = [1 0 0 0]';
d = [0 0]';
A = [1 3 0 0; 2 -3 0 0];
D = [1 0; 0 1];
B_1 = [1 0 1 0; 0 1 0 1];
b = [9/4 0]';
d_1 = [1 1]';

% Set the tolerance parameter
eps = 1e-6;

%%%%%%%%%%%%%%%%%%%
% The extreme points of
% X = {x| B_1x = d_1, x >= 0}
% are (0 0 1 1), (1 0 0 1), (0 1 1 0), and (1 1 0 0) respectively
% Let us start with the origin of X and the columns of D in the initial basis
% Use this information to construct A_B and c_B

A_B = [A*[0 0 1 1]'  D; 1 0 0];
c_B = [c'*[0 0 1 1]'; d];

btilde = [b; 1];

% Construct the first basic feasible solution for the refomulated problem
x_B = A_B\btilde

iter = 0;
x1 = [0 0 1 1]';

while 1 == 1,
iter = iter + 1;

% Update lower bounds as discussed in class
% The value of a basic feasible solution in every iteration provides this
% bound
if iter == 1,
  lb(iter) = c_B'*x_B;
else
  lb(iter) = max(lb(iter-1),c_B'*x_B);
end  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 2 :- Solve A_B^Ty = c_B
% Set up the objective function for the subproblem
% Subproblem is
% theta_1 = Min (c-A^Ty)^Tx 
%           s.t. B_1x = d_1, x >= 0
% Reduced costs are s_1 = theta_1 - z_1

y = A_B'\c_B;
z1 = y(3);
y = y(1:2);
s_1 = c-A'*y;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the subproblem using the revised simplex code on the course
% webpage.
% Use the previous primal optimal solution to the subproblem as the
% starting point

[obj1,x1,u1] = revised_simplex_nondeg(s_1,B_1,d_1,1e-6,x1);

% Update the upper bounds using the subproblem information
if (iter == 1),
  ub(iter) = b'*y + d_1'*u1;
else
  ub(iter) = min(ub(iter-1),b'*y + d_1'*u1);
end  

sN = [obj1-z1; d-D'*y];

[sNmax,k] = max(sN);

if sNmax <= eps,
   fprintf('We are done\n');
   fprintf('Number of iterations is %d\n',iter);
   fprintf('Optimal objective value is %f\n',c_B'*x_B);
   return;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step 3 :- Solve A_Bd = a_{N(k)}
% Find theta = Min_{i=1,...,m|d_i > 0} xB(i)/d(i)
% Let theta = xB(l)/d(l)
% x_{B(l)} is the leaving basic variable
% Also check for unboundedness if d <= 0

if (k == 1),
  aN = [A*x1; 1];
  temp = c'*x1;
elseif (k ==2),  
  aN = [D(:,1); 0]';
  temp = d(1);
else
  aN = [D(:,2); 0];
  temp = d(2);
end  

dcol = A_B\aN;
zz = find(dcol > eps)';

if (isempty(zz))
  error('System is unbounded\n');
end  

[theta,ii] = min(x_B(zz)./dcol(zz));
l= zz(ii(1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Update x_B, A_B, and c_B

A_B(:,l) = aN;
x_B = x_B - theta*dcol;
x_B(l) = theta;
c_B(l) = temp;

end; % while