function optx=plsimplex(d,p,c,m,a,lb,ub,t)
%
% Simplex code for the pseudo-linear hospital fee increase problem.
% (See paper for details)
% Input:
% d (scalar) number of departments
% p [1..d] number of procedures per departments
% nb (implicitely defined) is p[1]+p[2]+...+p[d] total number of procedures
% c [1..nb] Total charges 
% m [1..nb] Medicare charges 
% a [1..nb] Other charges
% lb [1..nb] Lower bound on x
% ub [1..nb] Upper bound on x
% t (scalar) Total charge increase
% Serge G. Kruk
% $Id: plsimplex.m,v 1.3 1999/03/03 15:10:55 sgkruk Exp $	
nb=sum(p);					  			% Number of procedures in total
optimal=0;					  			% Boolean for success
n=zeros(nb,1);				  			% This will hold gradient of objective
gamma=zeros(nb,1);			  			% Multiplier of upper bound constraints
lambda=zeros(nb,1);			  			% Multiplier of lower bound constraints
oldf=0;f=0;								% Objective function values
%
% Phase one: Find a feasible solution on a vertex or detect
% infeasibility.
%
x=lb;						  			% Try all at lower bound
if (c'*lb > t)				  			% Bad data, this has no solution
  disp('Infeasible program. Lower bound too high.');
  optx=[];
  return;
elseif (c'*ub < t)			  			% Bad data, this has no solution
  disp('Infeasible program. Upper bound too low.');
  optx=[];
  return;
else					   	  			% Feasible program
  for i=1:nb,				  			% Each variable in turn
	x(i)=ub(i);				  			% Move up to upper bound
	if (c'*x > t)			  			% If we moved too high, move back
	  x(i)=ub(i)-((c'*x-t)/c(i));k=i;break;
	end						  			% End of last variable moved up
  end						  			% End each variable
end							  			% End of feasible phase one
if (isempty(k))							% There is only one feasible point
  optx=x;return
end
x
%
% Phase two: Move from vertex to vertex until optimal solution is
% found.
%
%
while (~optimal)			  			% Until flag is changed
  %
  % Compute objective and gradient (will be stored in vector n[1..nb])
  %
  for i=1:d,				  			% For each department
	if (i==1), j=1; f=0; end	 
	v=[j:p(i)+j-1];			  			% Index of vars corr. to department i
	n(v)=a(v)+((c(v)'*x(v)).*m(v) - (m(v)'*x(v)).*c(v))/(c(v)'*x(v))^2;
	f=f+a(v)'*x(v)+(m(v)'*x(v))/(c(v)'*x(v));
	j=j+p(i);
  end
  if (oldf==f)
	disp('Degenerate case');
  elseif (oldf > f)
	disp('Objective is not pseudo-linear...');
	optx=x;								% Objective is NOT pseudo-linear.
	optimal=1;							% Say we are done
  end
  oldf=f;
  %
  % Compute multipliers
  %
  mu=-n(k)/c(k);			  			% Multiplier of equality constraint
  iu=find(x==ub);			  			% Index of upper bounded variables
  if (~isempty(iu)),
	gamma(iu)=zeros(size(iu)); 			% To satisfy complementarity
	lambda(iu)=n(iu)+mu.*c(iu);         % To satisfy null gradient Lagrangean 
  end
  il=find(x==lb);			  			% Index of lower bounded variables
  if (~isempty(il)),
	lambda(il)=zeros(size(il));			% To satisfy complementarity
	gamma(il)=-n(il)-mu.*c(il);         % To satisfy null gradient lagrangean 
  end
  gamma(k)=0;lambda(k)=0;				% Multipliers of slack constraints
  %
  % We have multipliers, check for optimality and improve if not attained
  %
  [value,cdropped]=min([lambda', gamma', 0]);
  if (value >= 0)			  			% All mutipliers non-negative
	disp('Optimal solution found');
	optx=x;optimal=1;		  			% We are done
  else						  			% Negative multiplier, we can improve
	%
	% Pick constraint to drop from most negative multiplier.
	% Find maximum step length to remain feasible and update variable
	% and index of slack variable.
	%
	oldx=x;								% Keep old solution around
	if (cdropped > nb)		  			% We must drop a variable upper bound
	  l=cdropped-nb;  					% Index of dropped constraint is l
	  if (((lb(k)-x(k))*c(k))/c(l) < x(l)-ub(l))
		x(k)=x(k)-c(l)*(ub(l)-x(l))/c(k); 
		x(l)=ub(l);			  			% Move to upper bound (exactly)
	  else					  
		x(l)=x(l)-c(k)*(lb(k)-x(k))/c(l); 
		x(k)=lb(k);			  			% Move to lower bound (exactly)
		k=l;		  		  			% To indicate that k is index of slack
	  end
	else
	  l=cdropped;	  					% Index of dropped constraint is l
	  if (((x(k)-ub(k))*c(k))/c(l) > lb(l)-x(l))
		x(l)=x(l)-c(k)*(ub(k)-x(k))/c(l); 
		x(k)=ub(k);			  			% Move to upper bound (exactly)
		k=l;		  		  			% To indicate that k is index of slack
	  else					  
		x(k)=x(k)-c(l)*(lb(l)-x(l))/c(k); 
		x(l)=lb(l);			  			% Move to lower bound (exactly)
	  end					  
	end
  end						  			% End of code on a negative multiplier
end							  			% End of while(!optimal)

