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

%	[Pareto_Fmat, Pareto_xmat, NBImultiplrs] = NBIgeneral5(x0, W, nvars, n, ne, ni, NBImultip, VLB, VUB, tol)
%
%	This program implements the shooting method for finding Pareto optimal
%	points of multi-objective problems 
%	Either the user is asked for weights w, each representing a starting point on 
%	the Quasi Pareto Plane, OR, the user can enter all the desired combinations 
% 	of weights in W, each **row** of W corresponding to a weighting strategy.
%	This version works only for  (n = 2), since the
%	weight-picking strategy is hard-wired to n = 2. CAN HANDLE nvars > n.
% 
%	This simply finds the optimal point for a given w starting from the optimal
%	solution of a nearby w.
%	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


 global nvars n ne ni shadowMin F_index  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

	if nargin<2,
		Wmat = input('Enter W with weights in rows  ');
	end
	Wmat = Wmat'; % now the weights are *** in columns ***	
	num_wts = size(Wmat,2);

%	normalize all weights
	for i = 1:num_wts
		wt = Wmat(:,i);

		Wmat(:,i) = wt/sum(wt); 

	end
%	sort in ascending order of w1: n = 2 hardwire
	Wmat(1,:) = sort(Wmat(1,:)); % sorting first row
	Wmat(2,:) = 1 - Wmat(1,:); % entering 1 - w(1) on second row
%	if [1 0] and [0 1] are in the set of weights, delete them
	
	if (Wmat(1,1) == 0), Wmat = Wmat(:,2:num_wts); end
	num_wts = size(Wmat,2); 
	if (Wmat(1,num_wts) == 1), Wmat = Wmat(:,1:(num_wts-1)); end 
%	update num_wts
	num_wts = size(Wmat,2);


	disp(Wmat');
%	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, which is t in the NBI subproblem
	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 % 

% 	FIND MATRIX PHI = [F(x1*) - F* | F(x2*)- 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

	w = ones(n,1)/n; 

	%% 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) = 150*(nvars + 1);


% Starting point for first NBI subproblem is the minimizer of f_2(x)
x_t = [shadow_x(n,:), 1e-3];
if NBImultip == 1, NBImultiplrs = [0;1]; end
W = [0;1];
Pareto_Fmat = (Phi(:,2)+shadowMin)';
Pareto_xmat = shadow_x(2,:);
Flop = [Flop, 0];

	for k = 1:num_wts

	% reset flops
  	 flops(0)

		w = Wmat(:,k)
	   
		[x_trial, fiasco] = constr_mo('max_tFUN',x_t, options, VLB, VUB);
	   
		if fiasco == 0,

			x = (x_trial(1:nvars))';
			xstar = (x_trial(1:nvars))';
			x_t = x_trial; % next starting point
			
		        if NBImultip  == 1,
				% Computing lambda using NBImultipliers.m
				lambda1 = NBImultipliers(x,VLB,VUB,normal,n,nvars,ne,ni); % lambda1 in column	
				NBImultiplrs = [NBImultiplrs, lambda1/sum(lambda1)];
			end

			F_Pareto = (Fns(x) + shadowMin)'
			x_Pareto = x';

			Pareto_Fmat = [Pareto_Fmat; F_Pareto]; % Pareto opt vectors in rows
			Pareto_xmat = [Pareto_xmat; x_Pareto]; % Pareto opt solutions in rows
			W = [W,w];


    			Flop = [Flop, flops]; %row

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

	end %for	


% Flop report

% tag on [1, 0] to the report



	W = [W, [1; 0]];
	NBImultiplrs = [NBImultiplrs, [1;0]];
	Pareto_Fmat = [Pareto_Fmat; (Phi(:,1)+shadowMin)'];
	Pareto_xmat = [Pareto_xmat; shadow_x(1,:)];
	Flop = [Flop, 0];


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


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

if (NBImultip == 1)
	% Report multipliers
	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


format short


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

% Plot if possible
  if size(Pareto_Fmat,2) == 2,
	plot(Pareto_Fmat(:,1), Pareto_Fmat(:,2), 'r*', Pareto_Fmat(:,1), Pareto_Fmat(:,2), 'b--')
	title('Points obtained using NBI')
	xlabel('f_1')
	ylabel('f_2')
	% 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 new 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
			

