function [x,T,B,A1,flg] = ILPsolve(A,b,c)
%
% This is a combined Phase I/Phase II calculation for solving ILP,
% integer linear programs. 
% We assume non-degeneracy in the entire computation.
% Written by Ming Gu for Math 170
% March 2007
%
% 
% Prepare for initial Phase I computation.
%
flg       = 0;
[m,n]     = size(A);
mask      = find(b<0);
A(mask,:) = - A(mask,:);
b(mask)   = - b(mask);
B         = (n+1:n+m)';
x         = [zeros(n,1);b];
c1        = [zeros(n,1);ones(m,1)];
y         = ones(m,1);
T         = [b eye(m); c1'*x,y'];
A1        = [A eye(m)];
%
% Starting Simplex Method
%
f = ['Starting Initial Phase I Simplex Iteration... '];
disp(f);
[ITER1,B,T,flg1] = ILPRevisedSimplex(A1,B,c1,T);
f = ['Phase I took ',num2str(ITER1),' Iterations. '];
disp(f);
if (flg1 ~= 0)
    f = ['Initial Phase I Failed, flg =  ',num2str(flg1)];
    disp(f);
    flg = flg1;
    return;
end
if (max(B)>n)
    f = ['Initial Phase I Detected Primal Infeasibility, max(B) = ',num2str(max(B))];
    disp(f);
    flg = 3;
    return;
end
x   = zeros(n,1);
x(B)= T(1:m,1);
y   = T(1:m,2:m+1)'*c(B);
T(end,:) = [c'*x,y'];
f = ['Starting Phase II/Phase I ILP Iteration... '];
disp(f);
ILP     = 1;
ILPITER = 0;
A1      = A;
b1      = b;
c1      = c;
xtol    = 1e-12;
bsave   = b;
while (ILP == 1)
    ILPITER = ILPITER + 1;
    x       = zeros(size(A1,2),1);
    x(B)    = A1(:,B)\bsave;
    T       =[x(B),inv(A1(:,B));x'*c1,c1(B)'/A1(:,B)];
    [ITER,B,T,flg2] = ILPRevisedSimplex(A1,B,c1,T);
    f = ['ILPITER = ',num2str(ILPITER), ' Phase II took ',num2str(ITER),' Iterations. '];
    disp(f);
    if (flg2 ~= 0)
        f = ['Phase II Failed, flg =  ',num2str(flg2)];
        disp(f);
        flg = flg2;
        return;
    end
    [m1,n1]=size(A1);
    x     = zeros(n1,1);
    x(B)  = T(1:m1,1);
    xb    = x(B);
    xfrac = xb - floor(xb);
    xfrac(xfrac>1-xtol) = 0;
    if (max(xfrac)< xtol)
        f = ['Integer Solution found in ',num2str(ILPITER), ' Iteations'];
        disp(f);
        x = round(x);
        x = x(1:n);
        ILP = 1;
        return;
    end
    [Y,I]               = min(abs(xfrac-1/2));
    beta                = zeros(1,size(T,1)-1);
    beta(I)             = 1;
    beta                = beta/A1(:,B);
    alpha               = beta * A1;
    alpha               = alpha - floor(alpha);
    alpha(alpha>1-xtol) = 0;
    bfrac               = xfrac(I);
    bsave               =[bfrac;bsave];
    Tsave = T;
    T                   = [bfrac 1, zeros(1,size(T,2)-1); 
                           T(1:end-1,1), zeros(size(T,1)-1,1) T(1:end-1,2:end); ...  
                           bfrac,1, zeros(1,size(T,2)-1)];
    Bsave               = B;
    B                   = [1;B+1];
    A1                  = [1 alpha, -1 ; zeros(size(A1,1),1) A1 zeros(size(A1,1),1)];
    c11                 = [1; zeros(size(A1,2)-1,1)];
    [ITER1,B,T,flg1] = ILPRevisedSimplex(A1,B,c11,T);
    y = T(end,2:end)';
    x = zeros(size(A1,2),1);
    x(B) = T(1:end-1,1);
    f = ['ILPITER = ',num2str(ILPITER), ' Phase I took ',num2str(ITER1),' Iterations. '];
    disp(f);
    if (flg1 ~= 0)
        f = ['Phase I Failed, flg =  ',num2str(flg1)];
        disp(f);
        flg = flg1;
        return;
    end
    if (min(B)==1)
        f = ['Phase I Detected Primal Infeasibility, min(B) = ',num2str(min(B))];
        flg = 1;
        disp(f);
        return;
    end
    A1  = A1(:,2:end);
    B   = B - 1;
    [m1,n1]  = size(A1);
    x   = zeros(n1,1);
    x(B)= T(1:m1,1);
    c1  = zeros(n1,1);
    c1(1:n) = c;
    y   = T(1:m1,2:m1+1)'*c1(B);
    T(end,:) = [c1'*x,y'];
end