format rat

disp('Example 1, p. 265 (KB)')
c = [5; 6; 0; 0];
A = [10 3 1 0; 2 3 0 1];
b = [52; 18];
basicvars = [3 4]
slackvars = [3 4]
[tableau,x,basic,feasible,optimal] = checkbasic1(A,b,c,basicvars);
tableau

disp('incoming 2')
disp('outgoing 4')
basicvars = [3 2]
[tableau,x,basic,feasible,optimal] = checkbasic1(A,b,c,basicvars);
tableau
optimal

disp('incoming 1')
disp('outgoing 3')
basicvars = [1 2]
[tableau,x,basic,feasible,optimal] = checkbasic1(A,b,c,basicvars);
tableau
optimal

disp('noninteger optimal solution found')
disp('use eq. 1 to find cutting plane')

A = [[A, zeros(2,1)]; [floor(tableau(1,1:end-1))-tableau(1,1:end-1) 1]]
b = [b;floor(tableau(1,end))-tableau(1,end)]
c = [c;0]
A(3,:)=A(3,:)+(1/8)*A(1,:)+(7/8)*A(2,:)
b(3)=b(3)+(1/8)*b(1)+(7/8)*b(2) % replace the new constraints in terms of the original variables
%(think of a problem in standard form) 
%An alternative way for doing this would have been to use checkbasic1.m with basicvars=[3 4 5] and
%then to extract a suitable matrix from the tableau.


basicvars = [1 2 5]
slackvars = [slackvars 5]
[tableau,x,basic,feasible,optimal] = checkbasic1(A,b,c,basicvars);
tableau
feasible

disp('nonfeasible basic solution. Convert to dual problem')
[dualA,dualb,dualc] = dualproblem(A,b,c,slackvars)
dualbasicvars = setdiff(1:5,basicvars)
[dualtableau,x,basic,feasible,optimal] = checkbasic1(dualA,dualb,dualc,dualbasicvars);
dualtableau

disp('incoming 5')
disp('outgoing 3')
dualbasicvars = [5 4];
[dualtableau,x,basic,feasible,optimal] = checkbasic1(dualA,dualb,dualc,dualbasicvars);
dualtableau
optimal
disp('Noninteger optimal solution of dual problem found')
disp('Convert to the primal problem')
basicvars = setdiff(1:5,dualbasicvars)
[tableau,x,basic,feasible,optimal] = checkbasic1(A,b,c,basicvars);
tableau

disp('Introduce a cutting plane using the second equation.')
A = [[A, zeros(3,1)]; [floor(tableau(2,1:end-1))-tableau(2,1:end-1) 1]]
b = [b;floor(tableau(2,end))-tableau(2,end)]
c = [c;0]
A(4,:) = A(4,:)+A(2,:)+(1/3)*A(3,:)
b(4)=b(4)+b(2)+(1/3)*b(3)

basicvars=[basicvars 6]
slackvars=[slackvars 6]
[tableau,x,basic,feasible,optimal] = checkbasic1(A,b,c,basicvars);
tableau

feasible

disp('The problem is infeasible as expected') 
disp('Convert to the dual problem')
[dualA,dualb,dualc] = dualproblem(A,b,c,slackvars);
dualbasicvars = setdiff(1:6,basicvars)
[dualtableau,x,basic,feasible,optimal] = checkbasic1(dualA,dualb,dualc,dualbasicvars);
dualtableau
optimal 

disp('incoming 6 outgoing 4')
dualbasicvars = [6 5]
[dualtableau,x,basic,feasible,optimal] = checkbasic1(dualA,dualb,dualc,dualbasicvars);
dualtableau
optimal 

disp('incoming 3 outgoing 5')
dualbasicvars = [6 3]
[dualtableau,x,basic,feasible,optimal] = checkbasic1(dualA,dualb,dualc,dualbasicvars);
dualtableau
optimal 

basicvars = setdiff(1:6,dualbasicvars)


[tableau,x,basic,feasible,optimal] = checkbasic1(A,b,c,basicvars);
tableau
feasible
optimal
disp('A noninteger optimal solution has been found')
disp('Construct cutting plane using eq. 1')
A = [[A, zeros(4,1)]; [floor(tableau(1,1:end-1))-tableau(1,1:end-1) 1]]
b = [b;floor(tableau(1,end))-tableau(1,end)]
c = [c;0]
A(5,:)=A(5,:)+(4/31)*A(1,:)+(28/31)*A(4,:)
b(5)=b(5)+(4/31)*b(1)+(28/31)*b(4)
basicvars=[basicvars 7]
slackvars=[slackvars 7]
[tableau,x,basic,feasible,optimal] = checkbasic1(A,b,c,basicvars);
tableau
feasible

disp('The problem is infeasible as expected') 
disp('Convert to the dual problem')
[dualA,dualb,dualc] = dualproblem(A,b,c,slackvars);
dualbasicvars = setdiff(1:7,basicvars)
[dualtableau,x,basic,feasible,optimal] = checkbasic1(dualA,dualb,dualc,dualbasicvars);
dualtableau
optimal 
disp('incoming 7 outgoing 3')
dualbasicvars = [7 6]
[dualtableau,x,basic,feasible,optimal] = checkbasic1(dualA,dualb,dualc,dualbasicvars);
dualtableau
optimal 

basicvars = setdiff(1:7,dualbasicvars)
[tableau,x,basic,feasible,optimal] = checkbasic1(A,b,c,basicvars);
tableau
feasible
optimal

disp('A noninteger optimal solution has been found')
disp('Construct cutting plane using eq. 2')
A = [[A, zeros(5,1)]; [floor(tableau(2,1:end-1))-tableau(2,1:end-1) 1]]
b = [b;floor(tableau(2,end))-tableau(2,end)]
c = [c;0]
A(6,:)=A(6,:)+(1/4)*A(5,:)
b(6)=b(6)+(1/4)*b(5)
basicvars=[basicvars 8]
slackvars=[slackvars 8]
[tableau,x,basic,feasible,optimal] = checkbasic1(A,b,c,basicvars);
tableau
feasible

disp('The problem is infeasible as expected') 
disp('Convert to the dual problem')
[dualA,dualb,dualc] = dualproblem(A,b,c,slackvars);
dualbasicvars = setdiff(1:8,basicvars)
[dualtableau,x,basic,feasible,optimal] = checkbasic1(dualA,dualb,dualc,dualbasicvars);
dualtableau
optimal

disp('incoming 8 outgoing 7')
dualbasicvars = [6 8]
[dualtableau,x,basic,feasible,optimal] = checkbasic1(dualA,dualb,dualc,dualbasicvars);
dualtableau
optimal

basicvars = setdiff(1:8,dualbasicvars)
[tableau,x,basic,feasible,optimal] = checkbasic1(A,b,c,basicvars);
tableau
feasible
optimal
disp('An integer optimal solution has been found with value 39')