function [Pareto_Fmat, Pareto_xmat, NBImultiplrs] = NBImorethan2(x0, spac_1, nvars, n, ne, ni, NBImultip, VLB, VUB, tol)

%	[Pareto_Fmat, Pareto_xmat] = NBImorethan2(x0, spac_1, nvars, n, ne, ni, NBImultip, VLB, VUB, tol)
%
%	This program implements Normal-Boundary Intersection for finding Pareto optimal
%	points of multi-objective problems 
%	The user enters an integer spac_1, 1/spac_1 being the uniform spacing, i.e.,
%	the possible values of any w_i are
%		[0, 1/k, 2/k, 3/k, ... , (k-1)/k, 1]
%
%	THIS VERSION CAN TACKLE MORE THAN 2 OBJECTIVES.
%	Weights are generated and picked using the tree approach.
%	Uses the same spacing for all w components, the spacing being 1/spac_1.
%	This version doesn't start on points on CHIM+ but not on CHIM.

% 
%	USES QUASI-NORMAL
% 	Copyright (c)	Indraneel Das, Rice University 

%	Function files needed:
%	Fns.m
%	max_tFUN.m (which calls FUN.m and Fns.m)
%	FUN1.m (for finding shadowmin)
%		max_tFUN.m doesn't need to be written afresh each time
%	Other problem non-independent files written:
%		Generate_wts(n,k)
%		num_subahead(n,k,w1,w2)


 global nvars n ne ni shadowMin F_index qpp_pt w  Phi normal

more off

	if nargin<3, nvars = input('Enter the number of variables in your problem   '); end
	if nargin<4, n = input('Enter the number of objective functions in your problem   '); end
	if nargin<5, ne = input('Enter the number of equality constraints      '); end
	if nargin<6, ni = input('Enter the number of inequality constraints      '); end
	if ne > 0,
		options(13) = ne;
	end
	
%	createshootfiles;

	if nargin < 1, x0 = input('Enter xo   '); end
	if length(x0) ~= nvars,
		disp(' The xo you entered has incorrect dimensions')
		x0 = input('Enter xo correctly  ');
	end


	if nargin < 7, NBImultip = 0;  end
	if nargin < 8, VLB = []; end
	if nargin < 9, VUB = []; end
	if nargin < 10, tol = 1e-4; end
	options(1) = 1;
	options(2) = tol; % constr.m tolerance for x
	options(3) = tol; % constr.m tolerance for f
	options(14) = 150*nvars; % max # of iterations

	x0 = assert_col_vec(x0);

	Phi = zeros(n);

%**********************************************

flops(0) 

   find_shadow = input('Find shadowmin or read from file ? (1/0)  ');


   if (find_shadow == 1)

%	FIND SHADOW_MIN
	shadowMin = zeros(n,1);
	shadow_x = zeros( n, nvars); % x*'s in rows
	disp('Finding Shadow Minimum')
	xo = x0'
	for i = 1:n
		
		fprintf(' Use the above as starting point for minimizing F(%i) ?\n', i);
		reply = input('Enter 1 if yes, 0 if you want to enter one yourself  ');
	if reply == 0,
	 	fprintf('Enter the starting point for minimizing F(%i)\n', i);
	  	xo = input('  ');
	end


	   	F_index = i; % using F_index, not function_index
	   	[shadow_x(i,:), fiasco, shadowMin(i)] = constr_mo('FUN1', xo, options, VLB, VUB);
		xo = shadow_x(i,:) % possibly a better option for the next starting point
	end

   else % read shadowmin from the function file containing the shadowmin data
	[shadowMin, shadow_x] = shadow_data;
	
   end % finding shadowMin afresh

% 	FIND MATRIX PHI = [F(x1*) - F* | F(x2*)- F* | ... | F(xn*) - F*];

	for i = 1:n
		Phi(:,i) = Fns(shadow_x(i,:));
		Phi(i,i) = 0;
	end


	shadowMin		
	Phi
	shadow_x

%	Check to make sure that QPP is n-1 dimensional
	if rcond(Phi) < 1e-8,
		disp(' Phi matrix singular, aborting.')
		Phi
		break
	end

	%% Finding quasi-normal instead of normal direction
	normal = -Phi*ones(n,1);
	normal = normal/norm(normal)

	reply = input('Waiting, hit return to continue  ...');

%*************** Making normal point in the right direction ****************


	if (sum(normal>0) > sum(normal<0)) ,
		gamma = -1;
	else
		gamma = 1;
	end

	b = gamma*normal;
	initial_flops = flops;

options(4) = 5e-7; % tol on constraint violation
options(13) = ne + n; % no. of equalities
options(14) = 100*(nvars + min(n,2) + 1);

% Starting point for first NBI subproblem is the minimizer of f_1(x)

x_t = [shadow_x(1,:), 0];
Pareto_Fmat = [];
Pareto_xmat = [];
StartingPts = [];

NBImultiplrs = [];
pi_vec = -pi*ones(1,nvars+1);


%	GENERATE WEIGHTS

	Wmat = Generate_wts(n,spac_1);
	num_wts = size(Wmat, 1);

%-----------------------------------------------------
	% initialize w_prev
	w_prev = Wmat(1,:);

	num_fiascos = 0;

  for k = 1:num_wts

	% reset flops
  	 flops(0)
	start_cputime = cputime;
	start_realtime = clock;

	%************  Pick weight from tree  *********************
		
	w_unsc = (Wmat(k,:))'; % primed to make w a column
	w = w_unsc/spac_1; % scaling to make it a convex combination 

	%******************************************************************

	% Solve problem only if it is not minimizing one of the individual objectives
	indiv_fn_index = find(w == 1);

	if indiv_fn_index ~= []
		% w has a 1 in indiv_fn_index th component, zero in rest
		% Just read in solution from shadow data

		W = [W, w];
		w_prev = w_unsc;

		Pareto_Fmat = [Pareto_Fmat; (Phi(:,i)+shadowMin)'];

		Pareto_xmat = [Pareto_xmat; shadow_x(i,:), 0];
		StartingPts = [StartingPts; shadow_x(i,:), 0];
		Flop = [Flop, 0];

		if NBImultip == 1, NBImultiplrs = [NBImultiplrs, w]; end

	else
		% FIND NEAREST SUBPROBLEM

		% Compare w_unsc and w_prev beyond 2nd element to find out the nearest subproblem,
		%		i.e. the starting x
		
		if (sum(w_unsc(3:n) - w_prev(3:n)) ~= 0)

		%	Number of rows in Pareto_xmat so far = k - num_fiascos -1
	
			

			num_wts_ahead = num_subahead(n, spac_1, w_prev, w_unsc);

			startx_t = StartingPts(k - num_wts_ahead, :);

			if abs(sum(startx_t - pi_vec)) < 1e-8
					% The solution of this subproblem was a failure
					% Pick nearest vertex of the simplex as the starting point

				[maximum_component, nearestVertexIndex] = max(w_unsc);
				startx_t = [shadow_x(nearestVertexIndex,:), mean(Pareto_xmat(:,nvars+1))];

					% uses average of all the current t's as the starting t


				fprintf(' Vertex %i picked as starting point \n\n', nearestVertexIndex);
				pause(1)
			end 
		else
			% if  w_unsc and w_prev don't differ beyond second element
			% Pick previous subproblem as nearest

			startx_t = StartingPts(k - 1, :);

			if (abs(sum(startx_t - pi_vec)) == 0)  
					% startx_t exactly pi_vec for a failed problem
					% The solution of this subproblem was a failure
					% Pick nearest vertex of the simplex as the starting point

				[maximum_component, nearestVertexIndex] = max(w_unsc);

				startx_t = [shadow_x(nearestVertexIndex,:), mean(Pareto_xmat(:,nvars+1))];

					% uses average of all the current t's as the starting t

				fprintf(' Vertex %i picked as starting point, previous subproblem failed\n\n', nearestVertexIndex);
				pause(1)
			end
			
		end

		% SOLVE NBI SUBPROBLEM
		w
		[x_t, fiasco] = constr_mo('max_tFUN', startx_t, options, VLB, VUB);
		


		if fiasco == 0
			F_Pareto = (Fns(x_t(1:nvars)) + shadowMin)'

			Pareto_Fmat = [Pareto_Fmat; F_Pareto]; % Pareto opt vectors in rows
			Pareto_xmat = [Pareto_xmat; x_t]; % Pareto opt solutions with t* in rows
			W = [W,w];
			w_prev = w_unsc;
			
    			Flop = [Flop, flops]; %row
			StartingPts = [StartingPts; x_t];

			if NBImultip  == 1,
				% Computing lambda using NBImultipliers.m
			lambda1 = NBImultipliers(x_t(1:nvars),VLB,VUB,normal,n,nvars,ne,ni) 
			%lambda1 in column

				NBImultiplrs = [NBImultiplrs, lambda1/sum(lambda1)];
	
			end

   		else
			num_fiascos = num_fiascos + 1;
			giveup = input('Give up ?  (0/1)  ');
			if giveup == 1, break; end	

			StartingPts = [StartingPts; pi_vec];

			% It is being falsely assumed that pi_vec is an impossible solution (a bug)

   		end % if fiasco == 0	

	end % if w is not one of the canonical vectors

  end %for	
	
  total_cputime = cputime - start_cputime;

% Flop report


disp('Report on number of floating point operations')
fprintf(' Number of flops in initial calculations: %i\n\n', initial_flops); 
Flop_report = [W; Flop/1000]';
net_flops = initial_flops + sum(Flop);
disp('	 beta	     #flops(in K)     ')
	
%format rat
disp(Flop_report)
net_flops

fprintf('\n\n Number of failures in solving NBI subproblems = %i\n\n',   num_fiascos);
fprintf('Number of NBI subproblems solved successfully = %i\n\n',   size(Pareto_xmat,1) );
fprintf('Net CPU time (secs) = %10.4f\n', total_cputime);
fprintf('Total elapsed time (secs) = %10.4f\n\n', etime(clock, start_realtime));

pause = input(' Waiting, hit return to continue.  ');


format short


if (NBImultip == 1)
	% Report multipliers
	num_wts = size(W)
	num_multip = size(NBImultiplrs)	
	Multip = [W;NBImultiplrs]; 

	disp('Corresponding weights in convex combinations problem')
	disp('    	   beta          	        w  ')
	disp(Multip')
	NBImultiplrs = NBImultiplrs'; % the output matrix has weights in rows
	pause = input(' Waiting, hit return to continue.  ');
	disp(' ')
end

%****************************************

% Plot if possible
  if size(Pareto_Fmat,2) == 3,
	plot3(Pareto_Fmat(:,1), Pareto_Fmat(:,2), Pareto_Fmat(:,3),'r*')
	%title('Pareto points obtained using NBImorethan2')
	xlabel('f_1')
	ylabel('f_2')
	zlabel('f_3')
	
	% Offer print options

save_plot = input('Save plot to a .ps file ? (1 for yes, anything else for no)  ');
if (save_plot == 1)
	enter_title = input('Want to enter a title for the plot ? (1 for yes) ');
	if (enter_title == 1)
		plot_title = input('Enter title in single quotes  ');
		eval(['title(' 'plot_title' ')']);
	end
	
	axes_relabel;

	filename = input('Enter name of file with .ps extension (in single quotes)   ');
	eval(['print -depsc ' filename]);
end

  end


  if size(Pareto_xmat,2) == 2,
  	reply = input('Want a plot of Pareto optimal x ? (0/1)   ');
  	if reply == 1,  plot(Pareto_xmat(:,1), Pareto_xmat(:,2), 'b*');	end
  % Offer print options

save_plot = input('Save plot to a .ps file ? (1 for yes, anything else for no)  ');
if (save_plot == 1)
	enter_title = input('Want to enter a title for the plot ? (1 for yes) ');
	if (enter_title == 1)
		plot_title = input('Enter title in single quotes  ');
		eval(['title(' 'plot_title' ')']);
	end

	axes_relabel;

	filename = input('Enter name of file with .ps extension (in single quotes)   ');
	eval(['print -depsc ' filename]);
end

  end
			



