function [lambda1] = NBImultipliers (x, lowbd, upbd, n_hat, n, nvars, ne, ni)

% Returns the least square multiplier estimates for the NBI subproblem
% Note that constraints are linear in t, hence t need not be entered 
% (working with gradients).
% Tries to find non-negative multipliers corresponding to the constraints 
% Phi w + t n_hat = F(x) (i.e. lambda1 >= 0)
%
% Problem to be solved:
% 	min_{multipliers} 	norm(grad_of_lagrangian(x)*mult)^2
%		s.t.     mu >= 0, lambda1 <= 0

global L nvars

gtol = 1e-5;
xtol = 1e-4;
residtol = 1e-4;

[gradF, gradh, gradg] = GradFhg(x); % Central Difference gradients
[fdummy, h, g] = FUN(x); % g needed to identify active constraints

L = [gradF, gradh];

for i = 1:ni
	% append gradients of active constraints to L
	if abs(g(i)) < gtol % constraint is binding
		L = [L, gradg(:,i)];
	end
end

% Take care of multipliers corresponding to active bounds

if ((lowbd ~= [])&(upbd ~= []))
    for i = 1:nvars
	% Inf is the same as 1e6 or more in magnitude
	
	% Check if lower bd is active
	if ((lowbd ~= [])&((lowbd(i) > -1e6)&( abs(x(i) - lowbd(i)) < xtol )))
		e_i = zeros(nvars,1);
		e_i(i) = -1;
		L = [L, e_i];
		elseif  ((upbd~=[])&((upbd(i) < 1e6)&( abs(x(i) - upbd(i)) < xtol )))
			e_i = zeros(nvars,1);
			e_i(i) = 1;
			L = [L, e_i];
	end 	

    end 
end

num_multipliers = size(L,2) % number of columns of L with nontrivial multipliers

last_row = [n_hat', zeros(1,ne), zeros(1, (num_multipliers - ne - n))];
		
L = [L; last_row];

% Problem:
% Find least squares solution to Ly = rhs (equiv to gradient of lagrangian = 0)
%			s.t. y(ne+1:num_multipliers) >= 0
%	where rhs = [zeros(nvars,1); 1]

vlb = [-1e10*ones(1,n+ne), zeros(1, (num_multipliers-ne-n))]';
vub = [zeros(1,n), 1e10*ones(1, (num_multipliers-n))]';


% Start at least square soln of Ly = rhs (w/o bounds on multipliers, that is)

y = L\[zeros(nvars,1);1];
Residual = norm(L*y - [zeros(nvars,1);1])

% vital_multips = y

% it is possible that the above solution falls within the right bounds and residual is small
% if not, minimize norm of gradient of lagrangian subject to bounds

if ~(( sum(y(1:n) > vub(1:n)) == 0)&(sum(y(n+ne+1:num_multipliers) < vlb(n+ne+1:num_multipliers)) == 0))|(Residual > residtol)

	y0 = y;
	disp('Minimizing norm of gradient of lagrangian to find multipliers') 
	options(1) = 0;
	options(3) = 0.1*residtol; % Stopping tolerance on f
	[y, fiasco, f] = constr_mo('NBImultplrFUN', y, options, vlb, vub);
	newResidual = sqrt(f);
	
	if (Residual <= residtol)&(newResidual > residtol)
		disp('Failed to find KKT multipliers, KKT point under suspicion') 
		disp(' ');	 
		% Relax nonpositivity constraints on lambda1 and rerun optimization
		vub(1:n) = 1e10*ones(n,1);
		disp('Relaxing bounds on lambda1 and minimizing ||gradl|| again')
		[y, fiasco, f] = constr_mo('NBImultplrFUN', y0, options, vlb, vub);
		newResidual = sqrt(f);
		
	end

	if (newResidual > residtol)
		% This cannot be a KKT point
		disp(' Your multipliers might be incorrect')
		lambda1 = y0(1:n);
		pause(2)
	else
		Residual = newResidual;
		lambda1 = y(1:n);
	end

else
	lambda1 = y(1:n);
end 
 

