function [ rates, avgRT, cdfRTindex, cdfRTvalue, sumOfRates, evaluations ] = findOptimalRates(qn, startingRates, constraints)
%FINDOPTIMALRATES Calculate the optimum rates
% INPUT:
% qn: queuing network
% startingRates: initial rates
% constraints: constraints on the response time
% OUTPUT:
% rates: new optimal rates
% avgRT: average respose time
% cdfRTindex: indexes of the response time distribution
% cdfRTvalue: values of the response time distribution
% sumOfRates: sum of the rates for non-delay nodes
% evaluations: number of evaluations of the fluid function

	% Options
	maxIterations = 100;
	maxError = 1e-4;

	
	% Initializations
	delayNodes = getDelayNodes(qn);
	fixedRates(qn.M,1) = false;
	visits = avgVisits(qn);
	evaluations = 0;
	
	% STEP 1: Find a feasible solution (not necessarily optimal) and fix the
	%			 rates of the delay nodes.
	rates = startingRates;
	initQNrates = avgRates(qn);
	rates(delayNodes) = 1;
	fixedRates(delayNodes,1) = true;
	
	% STEP 2: Scale all the non-fixed rates with maxError under maxIterations iterations
	%         using bisection.
	minRates = zeros(qn.M,1);
	maxRates = rates;
	iterations = 0;
	
	while iterations < maxIterations && any(~fixedRates)
		iterations = iterations + 1;
		rates(~fixedRates) = (maxRates(~fixedRates) + minRates(~fixedRates))/2;
		currentRates = initQNrates.*repmat(rates,1,qn.K);
		qn = setRates(qn, currentRates);
		[avgRT, cdfRTindex, cdfRTvalue] = evaluateRT(qn);
		evaluations = evaluations + 1;
		% minRates
		% maxRates
		% sumOfRates = sum(rates(~delayNodes))
		
		if constraints.verify(avgRT, cdfRTindex, cdfRTvalue)
			% Decrease the maximum value of all the rates
			%disp(['DECREASED all maxRates it=' num2str(iterations)]);
			maxRates = rates;
		else
			% Increase the minimum value of the bottleneck rates
			demands = avgDemands(qn, visits);
			maxDemands = max(demands(~delayNodes,:));
			bottleneck = false(qn.M, qn.K);
			bottleneck(~delayNodes, :) = demands(~delayNodes,:) == repmat(maxDemands, sum(~delayNodes),1);
			minRates(any(bottleneck,2)) = rates(any(bottleneck,2));
			%disp(['INCREASED minRate of node ' num2str(find(bottleneck)') ' it=' num2str(iterations)]);
		end
		
		% Calculate error
		error = max(max((maxRates(~fixedRates) - minRates(~fixedRates))/2));
		
		if error<=maxError
			% STEP 3: Fix the rate of the bottleneck resource
			distance = max(constraints.distance(avgRT, cdfRTindex, cdfRTvalue));
			bottleneckClass = find(distance == max(distance));
			bottleneckClass = bottleneckClass(1);
			demandsBottleneckClass = bottleneck(:, bottleneckClass);
			demandsBottleneckClass(fixedRates) = 0;
			rateToFix = demandsBottleneckClass == max(demandsBottleneckClass);			
			fixedRates(rateToFix) = true;
			% disp(['Fixed rate of node '  num2str(find(rateToFix)') ]);
		end
	end
	
	% STEP 4: Perform a last evaluation to compute output values
	rates = maxRates;
	qn = setRates(qn, initQNrates.*repmat(rates,1,qn.K));
	[avgRT, cdfRTindex, cdfRTvalue] = evaluateRT(qn);
	evaluations = evaluations + 1;
	sumOfRates = sum(rates(~delayNodes));
	
end

