function [sorted_vertices, ...
          h_fes, h_bnd, h_fill, h_vert, h_int, h_max, g_labels] = ...
    plot_feasible(A, b, c, lower_b, upper_b, varargin)
% Plot the feasible region of a linear program
%
% file:      	plot_feasible.m, (c) Matthew Roughan, 2009
% created: 	Mon Jun 22 2009 
% author:  	Matthew Roughan 
% email:   	matthew.roughan@adelaide.edu.au
% 
% Plot the feasible region of the 2D linear program
%           maximize            f = c'*x
%           subject to       A x <= b
% on the region bounded by lower_b and upper_b. Note that the inputs
%     A  is a mx2 matrix, where there are m constraints
%     b  is a mx1 column vector, where there are m constraints
%     c  is a 2x1 column vector
%     lower_b, upper_b are 2x1 column vectors
%
% If the problem includes non-negativity bounds, i.e., 
%               x >= 0 
% which is quite common, then these need to be explicitly included into
% the constraints:  A x <= b, not implicitly through the bounds, which are mainly there to
% allow us to make a nice plot.
%
% The feasible region is filled with lines at a setable angle and separation, much as fill
% would do with a solid colour, though options allow you to fill with a color, or both.
%
% This code is primarily (for me) useful for teaching about linear programming, and only
% works for 2D problems (for the obvious reason that problems with more than two variables
% are hard to plot).
%
% Note that if you want to cross-hatch just call this routine twice, with a change in angle
% of 90 degrees.
% 
% There are a bunch of optional inputs in pairs, as for the plot command
%    'linecolor'         color of feasible region boundary lines (default black)
%    'backgroundcolor'   background color of feasible region (default no background)
%    'linewidth'         width of lines (default 1)
%    'linestyle'         style of the boundary lines
%    'filllinestyle'     style of the fill lines
%    'linesep'           separation of cross-hatch lines (default 1)
%                        if this is set to -1, don't put any fill lines in
%    'lineangle'         angle (in degrees) of cross-hatch lines
%                           the default is to put them along lines of equal values of the
%                           objective function
%    'hold'              hold_n=1 means keep previous plot
%                        hold_n=0 (default) means clear previous figure
%                        use this to plot multiple feasible regions on same plot, or to do
%                        crosshatching.
%
%    'plot_vertices'     if defined indicates that we should plot the vertices of the feasible region
%                        the value will determine the symbol to use to plot the vertices,
%                        e.g., 'rx', or 'o+'
%    'label_vertices'    if not 0, then label the vertices with their position (label_vertices=1)
%                        or value (label_vertices=1) or both (label_vertices=3)
%    'label_vertices_size'  size of vertex labels, default = 12
%    'label_vertices_color' size of vertex labels, default = 'k'
%    'label_vertices_prec'  precision of vertex labels (default = 2 decimal places)
%
%    'plot_intersections' if defined indicates that we should plot the intersection points of boundaries
%                        the value will determine the symbol to use to plot the vertices,
%                        e.g., 'rx', or 'o+'
%                        Note that vertices are also intersection points, and will be double
%                        plotted if used with the plot_vertices option
% 
%    'plot_max'          if not 0, then plot the location of the maximum, using the defined symbol
%
%    'extend_boundaries' if defined, plot each of the constraint lines past the feasible region
%                        using the linewidth and color defined above. The style of the line can
%                        be given as the value, e.g., ':'
%
% OUTPUTS:
%     sorted_vertices = a 2x(n+1) vector of the vertices of the feasible region
%                        note the last one is a repeat of the first to aid in plotting the region
%     
%
%     h_fes    = handles to boundary lines of the linear program
%     h_bnd    = handles to the boundary lines
%     h_fill   = handles to the fill lines
%     h_vert   = handles to the vertices (if plotted)
%     h_int    = handles to the intersection points (if plotted)
%     h_max    = a handle to the plotted maximum point 
%     g_labels = handles to vertex labels if used
% 
% version 0.1, July 22nd 2009
%
% TODO
%     speckled fill
%     testing for boundedness and autosizing the figure
%     plotting bounds with cross-hatch on feasible side
%     at present the code uses the optimization toolkit (through the linprog) function
%       which we use to solve the linear program, but as this is a 2D problem we should be able to
%       remove this dependence with a little work.
%     currently the intersections with upper and lower bounds for plotting are also included
%       in intersection set, and plotted
%     automate ability to put text description of the problem in the window
%

% check sizes of inputs
sA = size(A);
sb = size(b);
if (sA(1) ~= sb(1))
  error('incorrect input sizes\n');
end
if (sA(2) ~= 2)
  error('can only plot for 2 variables\n');
end
n = sA(1); % number of constraints
if (size(upper_b,1) ~=2 | size(upper_b,2) ~= 1)
  error('upper_b is the wrong size\n');
end
if (size(lower_b,1) ~=2 | size(lower_b,2) ~= 1)
  error('lower_b is the wrong size\n');
end

% set default output values
h_fes    = [];
h_bnd    = [];
h_fill   = [];
h_vert   = [];
h_int    = [];
h_max    = [];

% parse the variable input arguments
linecolor = 'k'; 
backgroundcolor = 0;
linewidth = 1;
linestyle = '-';
filllinestyle = '-';
linesep = 1;
% default line angle is so that lines will be iso-objective
if (abs(c(1)) < 1.0e-12)
  lineangle = 90; % degrees
else
  lineangle = atand(c(2)/c(1));
end
hold_n = 0;
plot_vertices = 0;
label_vertices = 0;
label_vertices_size = 12;
label_vertices_color = 'k';
label_vertices_prec = 2;
plot_intersections = 0;
extend_boundaries = 0;
plot_max = 0;
if (length(varargin) > 0)
  % process variable arguments
  for k = 1:2:length(varargin)
    if (ischar(varargin{k}))
      argument = char(varargin{k}); 
      value = varargin{k+1};
      switch argument
       case {'linecolor','lc'}
	linecolor = value;
       case {'backgroundcolor','bgc'}
	backgroundcolor = value;
       case {'linewidth','lw'}
	linewidth = value;
       case {'linestyle','ls'}
	 linestyle = value;
       case {'filllinestyle','fls'}
	 filllinestyle = value;
       case {'linesep','ls'}
	linesep = value;
       case {'lineangle','la'}
	lineangle = value;
       case {'hold'}
	hold_n = value;
       case {'plot_vertices', 'pv'}
	plot_vertices = value;
       case {'label_vertices', 'lv'}
	label_vertices = value;
       case {'label_vertices_size', 'lvs'}
	label_vertices_size = value;
       case {'label_vertices_color', 'lvc'}
	label_vertices_color = value;
       case {'label_vertices_prec', 'lvp'}
	label_vertices_prec = value;
       case {'plot_intersections', 'pi'}
	plot_intersections = value;
       case {'extend_boundaries', 'eb'}
	extend_boundaries = value;
       case {'plot_max', 'pm'}
	 plot_max = value;
       otherwise
	error('incorrect input parameters');
      end
    end
  end
end

% set up the plot region
if (hold_n==0)
  hold off
  plot(lower_b(1), lower_b(2));
end
diff_x = (upper_b(1)-lower_b(1))/20; 
diff_y = (upper_b(2)-lower_b(2))/20; 
set(gca, 'xlim', [lower_b(1)-diff_x upper_b(1)+diff_x]);
set(gca, 'ylim', [lower_b(2)-diff_y upper_b(2)+diff_y]);
hold on
epsilon = 1.0e-14;

% plot the boundary line lines of the linear program across the plot region
if ~(extend_boundaries == 0)
  for i=1:n
    % find end-points of the lines along the boundaries
    constraint = A(i,:);
    k = b(i);
    if (abs(constraint(1)) < epsilon)
      % fprintf('(abs(constraint(1) < epsilon): i = %d\n', i);
      if (lower_b(1)<=k/constraint(2) & k/constraint(2)<=upper_b(1))
	x = [lower_b(1) upper_b(1)];
	y = [1 1] * k/constraint(2);
      else 
	x = [];
	y = [];
      end
    elseif (abs(constraint(2)) < epsilon)
      % fprintf('(abs(constraint(2) < epsilon): i = %d\n', i);
      if (lower_b(2)<=k/constraint(1) & k/constraint(1)<=upper_b(2))
	x = [1 1] * k/constraint(1);
	y = [lower_b(2) upper_b(2)];
      else 
	x = [];
	y = [];
      end
    else
      % compute all four intersections, and decide which is feasible
      % horizontal intercept points
      y_1 = [lower_b(2) upper_b(2)];
      x_1 = (k - constraint(2)*y_1) / constraint(1);
      i_1 = find(lower_b(1)<x_1 & x_1<upper_b(1));
      % vertical intercept points
      x_2 = [lower_b(1) upper_b(1)];
      y_2 = (k - constraint(1)*x_2) / constraint(2);
      i_2 = find(lower_b(2)<y_2 & y_2<upper_b(2));
      
      x = [x_1(i_1) x_2(i_2)];
      y = [y_1(i_1) y_2(i_2)];
    end
    % plot the line
    h_fes(i) = plot(x, y, 'color', linecolor, 'linestyle', extend_boundaries);
  end
end

% find the vertices of the lines, plus the borders
Aex = [A; -eye(2); eye(2)];
bex = [b; -lower_b; upper_b];
m = n + 4;
intersections = [];
vertices = [];
for i=1:m-1
  for j=i+1:m
    M = [Aex(i,:); Aex(j,:)];
    bm = bex([i,j]);
    x = M \ bm;
    intersections = [intersections, x];
    % test for feasibility
    if (A*x <= b)
      vertices = [vertices, x];
    end
  end
end

% find and fill the polygon defined by the vertices
k = convhull(vertices(1,:), vertices(2,:));
sorted_vertices(1,:) = vertices(1,k)';
sorted_vertices(2,:) = vertices(2,k)';
if (backgroundcolor == 0)
  h_bnd = plot(sorted_vertices(1,:), sorted_vertices(2,:),  ...
	       'color', linecolor, 'linewidth', linewidth, 'linestyle', linestyle);
else
  h_bnd = fill(sorted_vertices(1,:), sorted_vertices(2,:), backgroundcolor, ...
	       'edgecolor', linecolor, 'linewidth', linewidth, 'linestyle', linestyle);
end

% plot the intersection points and vertices
if ~(plot_intersections == 0)
  h_int = plot(intersections(1,:), intersections(2, :), plot_intersections);
end

if ~(plot_vertices == 0)
  h_vert = plot(sorted_vertices(1,:), sorted_vertices(2, :), plot_vertices);

  if (label_vertices == 1)
    % write the position of the vertex
    format_str = sprintf('(%%.%.0ff,%%.%.0ff)', label_vertices_prec, label_vertices_prec);
    for i=1:size(sorted_vertices, 2)
      temp = 10^label_vertices_prec;
      x = round(temp*sorted_vertices(1,i))/temp;
      y = round(temp*sorted_vertices(2,i))/temp;
      g_labels(i) = text(x+diff_x/2, y, sprintf(format_str, x, y), ...
			 'fontsize', label_vertices_size, 'color', label_vertices_color);
    end
  elseif (label_vertices == 2)
    % write the value of the vertex
    format_str = sprintf('f = %%.%.0ff', label_vertices_prec);
    for i=1:size(sorted_vertices, 2)
      temp = 10^label_vertices_prec;
      x = sorted_vertices(1,i);
      y = sorted_vertices(2,i);
      f = round(temp* (c' * [x; y]))/temp;
      g_labels(i) = text(x+diff_x/2, y, sprintf(format_str, f), ...
			 'fontsize', label_vertices_size, 'color', label_vertices_color);
    end
  elseif (label_vertices == 3)
    % write the position and value
    % write the position of the vertex
    format_str = sprintf('(%%.%.0ff,%%.%.0ff), f=%%.%.0ff', ...
			 label_vertices_prec, label_vertices_prec, label_vertices_prec);
    for i=1:size(sorted_vertices, 2)
      temp = 10^label_vertices_prec;
      x = sorted_vertices(1,i);
      y = sorted_vertices(2,i);
      f = round(temp* (c' * [x; y]))/temp;
      x = round(temp*x)/temp;
      y = round(temp*y)/temp;
      g_labels(i) = text(x+diff_x/2, y, sprintf(format_str, x, y, f), ...
			 'fontsize', label_vertices_size, 'color', label_vertices_color);
    end
  end

end

% plot the max if needed
if ~(plot_max == 0)
  x_max = linprog(-c, A, b, [], [], lower_b, upper_b);
  h_max = plot(x_max(1), x_max(2), plot_max);
end

% do the fill lines
if (linesep > 0)
  %   find min and max via linprog to see start and end points of lines
  %      add in implicit boundary constraints when do it
  c_line = [cosd(lineangle), sind(lineangle)];
  x_min = linprog(c_line, A, b, [], [], lower_b, upper_b);
  % plot(x_min(1), x_min(2), 'b*');
  
  x_max = linprog(-c_line, A, b, [], [], lower_b, upper_b);
  % plot(x_max(1), x_max(2), 'b*');
  
  distance = sqrt( sum((x_max-x_min).^2)  );
  
  %   draw lines clipping them at the boundaries
  for i=0:linesep:distance
    counter = ceil(i/linesep) + 1;
    
    % sigma = x_min + i*(x_max-x_min)/distance;
    sigma = x_min + i*c_line';
    % plot(sigma(1), sigma(2), '+', 'color', linecolor);
    
    % for each line, consider its intersepts with each possible 
    % boundary, and check whether they are feasible
    edge_vertices = [];
    for j=1:m
      M = [Aex(j,:); c_line];
      d = [bex(j); sum(c_line.*sigma')];
      x = M \ d;
      % test for feasibility
      if (A*x <= b)
	edge_vertices = [edge_vertices, x];
      end
    end
    edge_vertices = round(100*edge_vertices)'/100;
    edge_vertices = unique(edge_vertices, 'rows')';
    if (size(edge_vertices,2) >= 2)
      h_fill(counter) = plot(edge_vertices(1,:), edge_vertices(2,:), ...
			     'color', linecolor, 'linestyle', filllinestyle, 'linewidth', linewidth);
    end
  end
end