function [x,lambda] = nlp(x0,funobj,funcon,c,A,b1,b2,vlb,vub,itermx)
% NLP	Nonlinear Programming using MINOS.
%
% -----------------------------------------------------------------------
%
%	x = nlp(x0) solves the nonlinear programming problem:
%
%	minimize f(x) 	subject to: 	g(x) <= 0 
%	   x
%
%	where x0 is an n x 1 vector, initial starting point.  NOTE: x0
%	MUST be provided.
% 	
%	And the file 'funobj.m' is defined and contains:
%		function [f,gradf]=funobj(x) that computes the objective
%		function f and the gradient gradf at the point x.
%
%	And, the file 'funcon.m' is defined and contains:
%		function [g,gJac]=funcon(x) which computes the nonlinear
%		constraint value g and the constraint Jacobian gJac for
%		the nonlinear constraints at the point x.
%		NOTE: The matrix gJac MUST be a SPARSE MATLAB matrix.  This
%		      can be done by including the line: 
%				gJac = sparse(gJac);
% 		      after gJac has been computed.
%
% ------------------------------------------------------------------------
%
%	[x,lambda] = nlp(x0) returns the set of Lagrangian multipliers,
%	lambda, at the solution.
%
% ------------------------------------------------------------------------
%
%	x = nlp(x0,'funobj','funcon') solves the same nlp as above, the
%	only difference being that the computation of the objective 
%	function f and the gradient gradf at the point x is done by the
%	function [f,gradf] = funobj(x) defined in the file "funobj.m". 
%	The computation of the constraint value g and the Jacobian matrix 
%	gJac at x is done by the function [g,gJac] = funcon(x) defined in
%	the file "funcon.m".
%
% -------------------------------------------------------------------------
%
%	x = nlp(x0,'funobj','funcon',c) solves the following nlp:
%
%	minimize  f(x) + c'x	subject to:	g(x) <= 0
%	   x			
%
%	where c is a n x 1 vector of linear objective coefficients. And
%	f,gradf,g,gJac are computed by the functions [f,gradf] = funobj(x),
%	and [g,gJac] = funcon(x) in the files "funobj.m" and "funcon.m".
%
% -------------------------------------------------------------------------
%
%	x = nlp(x0,'funobj','funcon',c,A,b1,b2) solves the following nlp:
%
%	minimize  f(x) + c'x	subject to:	g(x) <= 0
%	   x				    b1 <= Ax <= b2
%
%	where A is an m x n matrix; b1 and b2 are m x 1 vectors defining
%	linear constraints and upper and lower bounds on the linear 
%	constraints.  For a description of 'funobj' and 'funcon' see
%	above.
%
% -------------------------------------------------------------------------
%
%	x = nlp(x0,'funobj','funcon',c,A,b1,b2,vlb,vub) solves the 
%	following nlp:
%
%	minimize  f(x) + c'x	subject to:	g(x) <= 0
%	    x				    b1 <= Ax <= b2
%					    vlb <= x <= vub
%
%	where vlb and vub are n x 1 vectors of upper and lower bounds on
%	the variables x.  For a description of 'funobj' and 'funcon' see
%	above.
%
% -------------------------------------------------------------------------
%
%	x = nlp(x0,'funobj','funcon',c,A,b1,b2,vlb,vub,itermx) sets the
%	maximum number of iterations to itermx.
%
% -------------------------------------------------------------------------
%
%	NOTE: if the problem has nonlinear constraints, try to provide an
%	      initial starting point x0 which is feasible!
%
%	NOTE: MINOS will create a print file and a summary file!
%
% -------------------------------------------------------------------------

% do the initial checks on the inputs
if ((nargin < 1) | (nargin > 10))
	error('Wrong number of input parameters');
end;

n = length(x0);
big = 1.0e20;
if (nargin < 10)
	itermx = 100000;
end 
if (nargin < 9)
	vub = big*ones(n,1);
else
	vub = vub(:);
	nvub = length(vub);
	if (nvub ~= n)
		error('vub MUST have the same number of elements as x0!')
	else
		vub = min(vub,big*ones(n,1));
	end % if
end % if
if (nargin < 8)
	vlb = -big*ones(n,1);
else
        vlb = vlb(:);
	nvlb = length(vlb);
	if (nvlb ~= n)
		error('vlb MUST have the same number of elements as x0!')
	else
		vlb = max(vlb,-big*ones(n,1));
	end % if
end % if
if ((nargin < 7)&(nargin >= 5))
	error('If linear constraints are present, they MUST be b1 <= Ax <= b2')
end % if ((nargin < 7)&(nargin >=5))
if (nargin >= 7)
	b2 = b2(:);
	m = length(b2);
	b2 = min(b2,big*ones(m,1));
end % if
if (nargin >= 6)
        b1 = b1(:);
	m = length(b1);
	b1 = max(b1,-big*ones(m,1));
end % if
if (nargin < 5)
	A = [];b1 = [];b2 = [];
end % if (nargin < 5)
if (nargin < 4)
	c = [];
end % if (nargin < 4)
if (nargin < 3)
	funcon = 'funcon';
end % if (nargin < 3)
if (nargin < 2)
	funobj = 'funobj';
end % if (nargin < 2)


% check to make sure A (if provided) does indeed have n columns
[m,n1] = size(A);
if (~isempty(A))
	m2 = length(b1);
	m3 = length(b2);
	if (m2 ~= m)
		error('b1 and b2 must have the same number of elements')
	end % if (m2~=m)
	if (n1~=n)
		error('A MUST have n columns!');
	end % if (n1~=n)
	if (m2~=m)
		error('A and b1,b2 MUST have the same number of rows!');
	end % if (m2~=m)
end % if (A~=[])

% randomly perturb the initial point to determine nonzero structure of Jacobian
xp = x0 + 1e-5*rand(n,1);
[g,gJac] = feval(funcon,xp);
m1 = length(g);
% check to see if gJac is sparse
if (~issparse(gJac))
	error('gJac MUST be SPARSE matrix.  Add the line: gJac = sparse(gJac) to <funcon>.m file');
end % if (~issparse(gJac))

[f,df] = feval(funobj,xp);
if (length(df)~=n)
	error('gradf MUST have the same number of elements as x0!')
end % if (length(df)~=n)

nnobj = n; % number of nonlinear variables
[nncon,nnjac] = size(gJac); % number of nonlinear constraints (nncon)
			    % number of nonlinear jacobian vars (nnjac)

if ((m1 > 0)&(nnjac ~= n))
	error('Jacobian matrix must have n columns.')
end % if ((m1 > 0)&(nnjac~=n))

% set up the constraint matrix a and append the linear coefficients onto
% the bottom!
c = c(:);
if (m == 0 & m1 == 0 & isempty(c))
  c = zeros(n,1);
end
a = sparse([gJac;A;c']);
ne = nnz(a);
if (~isempty(c))
	iobj = (m+m1+1);
else
	iobj = 0; % there are no linear terms in objective
end % if (c~= [])

% to modify this code to put upper and lower bounds on the nonlinear
% constraints g(x), say  r <= g(x) <= t, where r and t are column m1-vectors
% then uncomment the following lines and comment the next set of "if" 
% statements, and change the calling sequence to include r and t:
% if (A==[])
%	bl = [vlb;-t];
%	bu = [vub;-r];
% end
% if ((c==[])&(A~=[]))
%	bl = [vlb;-t;-b2];
%	bu = [vub;-r;-b1];
% end 
% if ((c~=[]))
%	bl = [vlb;-t;-b2;-big];
%	bu = [vub;-r;-b1;big];
% end 

% set up bl,bu;
if (isempty(A)) % there are no linear constraints!
	bl = [vlb;zeros(m1,1)];
	bu = [vub;big*ones(m1,1)];
end % if (A == [])
if ((isempty(c))&(~isempty(A)))
	bl = [vlb;zeros(m1,1);-b2];
	bu = [vub;big*ones(m1,1);-b1];
end % if ((c==[])&(A~=[]))
if ((~isempty(c)))
	bl = [vlb;zeros(m1,1);-b2;-big];
	bu = [vub;big*ones(m1,1);-b1;big];
end % if ((c~=[]))

% call nlpmex!
[x,lambda,inform] = minosm(funobj,funcon,a,x0,bl,bu,ne,nncon,nnobj,nnjac,iobj,itermx,gJac);

if (inform == 1)
	disp('MINOS: inform = 1, Problem Infeasible')
	disp('	     NOTE: Due to the presence of nonlinear constraints,')
	disp('		   MINOS may determine that the problem is')
	disp('		   infeasible even though it is feasible.  Try')
	disp('		   optimizing from a feasible starting point.')
	disp('		   Refer to MINOS manual.')
end % if (inform == 1)
if (inform == 2)
	disp('MINOS: inform = 2, Problem Unbounded')
	disp('	     NOTE: Due to presence of nonlinear objective,')
	disp('		   MINOS may determine that the problem is')
	disp('		   unbounded even though it is not unbounded.')
	disp('		   Refer to MINOS manual.')
end % if (inform == 2)
if (inform > 2 | inform < 0)
	disp(['MINOS: inform = ' num2str(inform)])
	disp('	     NOTE: Refer to MINOS manual.')
end % if (inform > 2)

x = x(1:n);
lambda = lambda(1:(m+m1)); 

return

% !
