function [exitFlag exitMsg netStruct] = fa_minimaxUtilization_xdp (netStruct, algorithmParameters, net2planParameters)

	% Algorithm for flow assignment problems which tries to minimize
	% the maximum utilization among all links within the network. It uses a
	% flow-path formulation implemented in CVX. For each demand, the set of
	% admissible paths is composed of k-shortest paths (loopless) between end nodes,
	% in number of hops. k is an input parameter to the algorithm
	%
	% This algorithm has no specific parameters.
	%
	% It requires CVX solver installed and running.
	%
	% PARAM: k | 10 | Number of admissible paths for each demand
	
	%% Sanity checks
	k = str2double(algorithmParameters.k);
	assert(isCvxInstalledAndRunning(), 'CVX solver must be installed and running');
	assert(k > 0, 'k should be positive');
	assert(k == floor(k) , 'k should be integer');

	assert(isfield(netStruct, 'linkTable'), 'A set of links is required');
	assert(isfield(netStruct, 'demandTable'), 'A set of demands is required');
	assert(isCvxInstalledAndRunning(), 'CVX solver must be installed and running');

	% Initialize some inputs
	E = size (netStruct.linkTable , 1); 
	D = size (netStruct.demandTable , 1); 
	u_e = netStruct.linkCapacityInErlangs;
	h_d = netStruct.offeredTrafficInErlangs;

	% Construct the path list
	[netStruct.pathList delta_dp] = constructPathList_kShortest (netStruct.linkTable , netStruct.linkLengthInKm , netStruct.demandTable , k);
		
	% Initialize some variables from the path list
	P = numel (netStruct.pathList);

	% delta_pe (PxE) = number of times a path p traverses a link e 
	delta_pe = seqLinksPerPath2linkOccupancyPerPath(netStruct.pathList, E);

	cvx_begin

		cvx_quiet(true);

		variables rhoMax; % maximum utilization in the network
		variables x_dp(D,P); % fraction of traffic h_d of demand d carried by path p 

		minimize( rhoMax )

		subject to

		% 1. If the path is not admissible for demand d, it does not carry traffic of d
		% x_dp = 0 for all p not in P_d, for all d
		% a. vectorial form (all constraints in one single line of code)
		x_dp <= delta_dp;
		% b. other simple form for the same
		% x_dp (find (delta_dp == 0) ) == 0;
		% c. equivalent expression with for-loops (slower)
%			 for d=1:D
%				 for p=1:P
%					 if (delta_dp(d,p) == 0)
%						 x_dp(d,p) == 0;
%					 end
%				 end
%			 end

		% 2. The sum of the traffic carried by the admissible paths = 100%
		% sum(p in P_d) x_dp = 1, forall d
		% a. vectorial form (all constraints in one single line of code)
		sum (x_dp .* delta_dp , 2) == 1;
		% b. equivalent expression with for-loops (slower)
%			 for d=1:D
%				 P_d = find (delta_dp (d,:) == 1);
%				 if isempty(P_d), continue; end
%				 sum(x_dp(d,P_d)) == 1; % sums in all the coordinates given by P_d
%			 end
			
		% 3. Constraints to limit the utilization in the links
		% y_e <= u_e * rhoMax, forall e
		% a. vectorial form (all constraints in one single line of code)
		h_d * x_dp * delta_pe <= u_e * rhoMax;
		% b. equivalent expression with for-loops (slower)
%		 for e=1:E
%			 trafInLink = 0;
%			 for d=1:D
%				 P_de = intersect (find(delta_dp(d,:) == 1) , find (delta_pe(:,e) == 1))
%				 if isempty(P_de), continue; end
%				 trafInLink = trafInLink + h_d(d) * sum(x_dp(d,P_de));
%			 end
%			 trafInLink <= u_e(e) * rhoMax;
%		 end

		% 4. Traffic can't be more than capacity in any link
		rhoMax <= 1; 


		% 5. Non-negativity constraints (vectorial form)
		x_dp >= 0;
		% rhoMax >= 0; % this is redundant, so can be removed

	cvx_end

	if strcmp(cvx_status, 'Solved') || strcmp(cvx_status, 'Inaccurate/Solved')

		exitFlag = 0;
		exitMsg = 'OK!';

		x_dp(x_dp > 1-net2planParameters.precisionFactor) = 1;
		x_dp(x_dp < net2planParameters.precisionFactor) = 0;

		netStruct.routingMatrix = x_dp;

	else

		exitFlag = -1;
		exitMsg = cvx_status;
		netStruct = [];

	end
		

		

end

function  [pathList delta_dp] = constructPathList_kShortest (linkTable , costPerLink , demandTable , k)
		
	D = size(demandTable,1);
	pathList = cell (1,0); % initialize as empty array
	delta_dp = sparse (D,0); % initialize as 0-columns matrix
	for d=1:D
		cellArraySequenceOfLinks = kShortestPath(linkTable , costPerLink , demandTable(d,1) ,  demandTable(d,2) , k);			

		for cont_p=1:numel(cellArraySequenceOfLinks)
			% if empty, then there is no path => go to next demand
			if isempty(cellArraySequenceOfLinks{cont_p}), continue; end 
			% add the path at the end of the pathList
			pathList{end+1} = cellArraySequenceOfLinks{cont_p};
			% update delta_dp. This automatically adds one column to the
			% right. The number of columns of delta_dp is equal to the
			% number of paths in pathList
			delta_dp (d,end+1) = 1;
		end
	end
end


