This was my work in developing the plotregion function in Julia



In [4]:

    
# Pkg.add("QuickHull")
using QuickHull
using Combinatorics









    



INFO: Recompiling stale cache file /home/juser/.julia/lib/v0.5/FixedSizeArrays.ji for module FixedSizeArrays.



In [10]:

    
collect(combinations(1:5,3))









    Out[10]:





10-element Array{Array{Int64,1},1}:
 [1,2,3]
 [1,2,4]
 [1,2,5]
 [1,3,4]
 [1,3,5]
 [1,4,5]
 [2,3,4]
 [2,3,5]
 [2,4,5]
 [3,4,5]



In [11]:

    
setdiff(1:5,[1,2,3])









    Out[11]:





2-element Array{Int64,1}:
 4
 5

%% Define the LP
% minimize    -x1 - 2*x2
% subject to  -2*x1 + x2 <= 2
%             -x1 + 2*x2 <= 7
%              x1        <= 3
%              x1, x2 >= 0
%
% Except, we need this in standard form:
% minimize f'*x subject to. A*x <= b

c = [-1 -2];
A = [-2 1; -1 2; 1 0];
b = [2; 7; 3];
lb = [0; 0];
ub = []; % No upper bound

%% Find the solution using linprog
x = linprog(c,A,b,[],[],lb);

%% Define the full LP with all the slacks.
c = [-1 -2 0 0 0];
A1 = [-2 1; -1 2; 1 0];
b = [2; 7; 3];
A = [A1 eye(3)];
%x = linprog(c,[],[],A,b,zeros(5,1));



In [13]:

    
function basic_feasible_point(A::Matrix,b::Vector,set::Vector)
    m,n = size(A)
    @assert length(set) == m "need more indices to define a BFP"
    binds = set # basic variable indices
    ninds = setdiff(1:size(A,1),binds) # non-basic
    B = A[:,binds]
    N = A[:,ninds]
    #cb = c[binds]
    #cn = c[ninds]
    
    xb = B\b
    x = zeros(eltype(xb),n)
    x[binds] = xb
    x[ninds] = zero(eltype(xb))    
    
    if any(xb .< 0)
        return (:Infeasible, x)
    else
        #lam = B'\cb
        #sn = cn - N'*lam
        return (:Feasible, x)
    end
end
A1 = [-2.0 1; -1 2; 1 0]
b = [2.0; 7; 3]
A = [A1 eye(3)]
basic_feasible_point(A,b,[1,2,3])









    



WARNING: Method definition basic_feasible_point(Array{T<:Any, 2}, Array{T<:Any, 1}, Array{T<:Any, 1}) in module Main at In[12]:2 overwritten at In[13]:2.






    



UndefVarError: c not defined

 in basic_feasible_point(::Array{Float64,2}, ::Array{Float64,1}, ::Array{Int64,1}) at ./In[13]:8



In [ ]:



In [ ]:

    
"""
Plot the feasible polytope 
Ax = b, x >= 0
for the first two components of x.
"""
function plotregion(A::Matrix,b::Vector)
    m,n = size(A)
    for inds in combinations(1:n,m)
        bfp=basic_feasible_point(A,b,inds)

%% Look at all vertices
nchoosek(5,3)

sets = nchoosek(1:5,3);

xlim([-10 10]);
ylim([0,8]);

for i=1:size(sets,1)

    
        lam = B'\cb';

        sn = cn' - N'*lam;
       
        
        x = zeros(5,1);
        x(bi) = xb;
        x(ni) = 0;
        
        fprintf('set %i, %i, %i is   feasible  %s  \n', SB(1), SB(2), SB(3), num2str(x'));
        
        sn

        plot(x(1),x(2),'*','MarkerSize',25);
    end
end