% ellipsoid algorithm solves LP
% x >= 0, Ax <= b, min < c,x >
% with b >= 0
% the last row of A has positive coefficients
n = 2; % dimension of space
m = 4; % number of additional inequality constraints
A = randn(m,n);
A(end,:) = abs(A(end,:));
b = abs(randn(m,1));
c = randn(1,n);
tic
xbest = zeros(n,1);
fbest = 0;
R = b(end)/min(A(end,:));
x = zeros(n,1);
P = R^2*eye(n);
detPold = R^(2*n);
fprintf("Iter. |      det(P) | viol. constr. | function val | best val.\n")
count = 0;
while true
    count = count + 1;
    d = b - A*x;
    [mi,minind] = min([d; x]);
    if mi >= 0 % no constraint violated
        fval = c*x;
        if fbest > fval
            xbest = x;
            fbest = fval;
        end
        g = c;
        minind = 0;
        fprintf("    %d | %5e |       %d | %5e | %5e\n",count,detPold,minind,fval,fbest);
    else % constraint minind violated most
        if minind <= m
            g = A(minind,:);
        else
            g = zeros(1,n);
            g(minind-m) = -1;
        end
        fprintf("    %d | %5e |       %d |               | %5e\n",count,detPold,minind,fbest);
    end
    Pg = P*g';
    gPg = g*Pg;
    x = x - 1/(n+1)*Pg/sqrt(gPg);
    P = n^2/(n^2-1)*(P - 2/(n+1)*(Pg*Pg')/(gPg));
    detP = det(P);
    if (detPold < detP) || (detP <= 0) || (gPg <= 0)
        break;
    else
        detPold = detP;
    end
end
ellipsoid_time = toc;
tic
cvx_begin
    variable x(n,1) nonnegative
    A*x <= b
    minimize( c*x )
cvx_end
cvx_time = toc;
fprintf("Best solution ellipsoid method:\n")
xbest,
fprintf("cvx solution:\n")
x,
fprintf("Ellipsoid method time: %5e sec\n",ellipsoid_time)
fprintf("cvx time: %5e sec\n",cvx_time)
