function [ d, cost, avgRT, cdfRTindex, cdfRTvalue, evaluations, time ] = ... 
		findOptimalAllocationFull(qn, constraints, resources, seed)
% FINDOPTIMALALLOCATIONFULL
% Calculate an optimal allocation for the stations of a queuing network on a set
% of resources, subject to some constraints on the response time.
% INPUT
% qn: input queuing network
% constraints: response time constraints
% resources: types of available resources
% seed: random number generator seed
% OUTPUT
% d: allocation matrix
% cost: cost fo allocating the allocation matrix
% avgRT: mean response time
% cdfRTindex: response time percentiles indexes
% cdfRTvalue: response time percentiles values
% evaluations: number of evaluations
% time: time to computer the solution.

	time = cputime;
	delayNodes = qn.S'==-1;
	found = false;
	evaluations = 0;
	tolerance = 1;
	
	rng(seed);
	
	% Calculate feasible rates
	startingRates = findFeasibleRates(qn, constraints)
	for ii=1:numel(startingRates)
		startingRates(ii) = startingRates(ii)*(1+rand());
	end
	
	%startingRates
	
	% Calculate optimal rates
	[ rates, avgRT, cdfRTindex, cdfRTvalue, ~, ievaluations ] = findOptimalRates(qn, startingRates, constraints);
	evaluations = evaluations + ievaluations;

	%rates
	
	while ~found
		disp('Macroiteration STARTED');
		
		% Calculate optimal allocation
		[ d, cost ] = findOptimalAllocation(rates(~delayNodes), resources);
		
		% Check if the constraints are still fulfilled
		[constraintViolation, avgRT, cdfRTindex, cdfRTvalue] = ... 
			evaluateAllocatedQN(qn, rates(~delayNodes), d, ...
					constraints, resources);
		evaluations = evaluations+1;

		% dynamic increment :)
		increment = max(max((constraints.value() + ...
					constraints.distance(avgRT, cdfRTindex, cdfRTvalue))./ ...
							constraints.value()));
		% If not found, try to increase the rates by 1% on each resource at the time
		% and use the one that violates less the constraints.
		if constraintViolation<=tolerance
			found = true;
			continue;
		end
		
		disp(['Macroiteration ENDED with violation=' ... 
			num2str(constraintViolation) ' increment=' num2str(increment)]);
		bestConstraintViolation = constraintViolation;
		
		for ii=1:sum(~delayNodes)
			currentRates = rates(~delayNodes);
			currentRates(ii) = currentRates(ii)*1.1;
			[ currentD, ~ ] = ... 
				findOptimalAllocation(currentRates, resources);
			[ currentConstraintViolation, ~, ~, ... 
				~] = evaluateAllocatedQN(qn, currentRates, ... 
				currentD, constraints, resources);
			evaluations = evaluations+1;

			disp(['Microiteration ' num2str(ii) '/' sum(~delayNodes) ...
				' ENDED with violation=' num2str(currentConstraintViolation)]);
			
			if currentConstraintViolation<=bestConstraintViolation
				bestRates = currentRates;
				bestConstraintViolation = currentConstraintViolation;
			end
		end

		if bestConstraintViolation == constraintViolation
			rates(~delayNodes) = rates(~delayNodes).*increment;
		else
			rates(~delayNodes) = bestRates;
		end
	end

	%d = allocation.*repmat(rates(~delayNodes),1,sum(needed));
	time = cputime - time;
end

function [res, avgRT, cdfRTindex, cdfRTvalue] = ...
	evaluateAllocatedQN(qn, rates, d, constraints, resources)
% EVALUATEALLOCATEDQN
% INPUT
% qn: queuing network
% rates: rates allocated to the stations (nonDelayMx1)
% needed: number of needed resourches for each kind
% constraints for the queuing network
% OUTPUT
% res: negative if the network satisfies the constraints, positive if at least a
% constraint is violated. The actual positive value is equal to the gap of the
% most violating constraint.

	delayNodes = getDelayNodes(qn);
	currentRates = avgRates(qn);
	currentRates(~delayNodes,:) = currentRates(~delayNodes,:) .* repmat(rates,1,qn.K); 
	
	nqn = setRates(qn, currentRates);
	nqn = allocateQN(nqn, d, resources);
	[avgRT, cdfRTindex, cdfRTvalue] = evaluateRT(nqn);
	
	res = max(max(constraints.distance(avgRT, cdfRTindex, cdfRTvalue)));
	
	%constraintDistance = constraints.distance(avgRT, cdfRTindex, cdfRTvalue)
end

