function [ d, cost ] = findOptimalAllocation(rates, resources)
%FINDOPTIMALALLOCATION Allocate rates to resources
% INPUT
%   rates: rate needed for each station (vector of size M=number of stations)
%   resources: available resources (instance of resources class)
% OUTPUT
%	 d: allocation matrix
%   cost: total cost for renting the resources
	
	% Calculate needed resources: needed
	% objective is: minimize the sum needed*resourcecost
	% constraint is: sum of needed*resourcerate > rates
	% needed are non-negative integer value.
	totalRate = sum(rates);
	[needed, cost] = intlinprog(resources.cost,1:resources.size, ...
		-resources.ecu,-totalRate,[],[],zeros(resources.size,1),[]);
	needed = round(needed)';
	
	
	% Find the allocation matrix with the following properties:
	
	% Minimize the number of partitions of each station.
	% Maximize the size of each partition.
	
	m = length(rates);
	r = resources.size;
	s = sum(needed);
	allocation = zeros(m,s);
	remainingRates = rates;
	remainingNeeded(s)=0;
	kk=0;
	for ii=1:r
		remainingNeeded(kk + (1:needed(ii))) = resources.ecu(ii);
		kk = kk + needed(ii);
	end
	
	while any(remainingRates>0)
		highestRate = find(remainingRates==max(remainingRates));
		highestRate = highestRate(1);
		highestNeeded = find(remainingNeeded==max(remainingNeeded));
		highestNeeded = highestNeeded(1);
		transfer = min([remainingRates(highestRate) remainingNeeded(highestNeeded)]);
		remainingRates(highestRate) = remainingRates(highestRate) - transfer;
		remainingNeeded(highestNeeded) = remainingNeeded(highestNeeded) - transfer;
		allocation(highestRate,highestNeeded) = transfer/rates(highestRate);
	end
	
	d = allocation.*repmat(rates,1,sum(needed));
	
end
