function [qnRefactored, nrefactorings, maxRep] = softwareRefactor(qn, scAlternatives, ...
					mergeAlternatives, maxDuplications, maxJobReassignments)
	% SOFTWAREREFACTOR
	%
	% This function combines the three types of software refactorings according
	% to the following algorithm:
	%
	% 1. Apply all speed-increasing substitutions
	% 2. Apply all speed-increasing merges
	% 3. If at least a modification has been done, go back to step 1
	% 4. Apply job reassignment
	% 5. Return the result
	%
	% INPUT
	% qn: queuing network
	% 
	% scAlternatives: array of cells in which each cell contains a matrix. Each
	% row of the internal matrix contains alternative rates for each class, columns are the classes.
	% The indexes of the array are the ones of the replaceable components.
	% components.
	%
	% mergeAlternatives: matrix of cells in which each cell contains a matrix. Each
	% row of the internal matrix contains alternative rates for each class, columns are the classes.
	%
	% The indexes of the matrix are the ones of the mergeable components.
	% components.
   % maxDuplications: maximum number of additional duplications for the current queuing
	% center (vector of qn.M size)
	% maxJobReassignments: max number of iterations for reassigning jobs
	%
	% OUTPUT
	% qnRefactored: new qn after applying refactorings
	% nrefactorings: number of refactorings that have been performed
	% maxRep: replicability (number of max replica) for the queueing centers of
	%			 the resulting queueing network
	
	nrefactorings = 0;
	currentRefactorings = 1;
	myqn = qn;
	selected = true(1,qn.M);
	
	while currentRefactorings > 0
		currentRefactorings = 0;
		myScAlt = scAlternatives(selected);
		myMgAlt = mergeAlternatives(selected,selected);
		
		[myqn, nref] = softwareRefactorSCReplacement(myqn, myScAlt);
		nrefactorings = nrefactorings + nref;
		currentRefactorings = currentRefactorings + nref;
		
		[myqn, nref] = softwareRefactorSCMerge(myqn, myMgAlt);
		nrefactorings = nrefactorings + nref;
		currentRefactorings = currentRefactorings + nref;
		
		% remove queueing centers that have been merged
		selected = false(1, qn.M);
		
		for ii=1:qn.M
			for jj=1:myqn.M
				if strcmp(qn.nodeNames{ii},myqn.nodeNames{jj})
					selected(ii) = true;
				end
			end
		end
	end
	
	myMaxDuplications = zeros(myqn.M, 1);
	
	for ii=1:qn.M
		currentLabel = qn.nodeNames{ii};
		duplications = maxDuplications(ii);
		jj=0;
		
		while duplications>0
			jj = mod(jj, myqn.M)+1;
			
			if strcmp(myqn.nodeNames{jj}, currentLabel)
				myMaxDuplications(jj) = myMaxDuplications(jj) + 1;
				duplications = duplications - 1;
			end
		end
	end
	
	[qnRefactored, nref] = softwareRefactorJobReassignment(myqn, myMaxDuplications, ...
		maxJobReassignments);
	
	nrefactorings = nrefactorings + nref;
	
	% CLEAN ME, recalculate maximum duplications (this contains duplicated code
	% that can be optimized
	myMaxDuplications = zeros(qnRefactored.M, 1);
	
	for ii=1:qn.M
		currentLabel = qn.nodeNames{ii};
		duplications = maxDuplications(ii);
		jj=0;
		
		while duplications>0
			jj = mod(jj, qnRefactored.M)+1;
			
			if strcmp(qnRefactored.nodeNames{jj}, currentLabel)
				myMaxDuplications(jj) = myMaxDuplications(jj) + 1;
				duplications = duplications - 1;
			end
		end
	end
	maxRep = myMaxDuplications;
	
end