HPSR2001_minTransceivers(traff_trafficMatrix, phys, parameters)

% Algorithm for traffic grooming in optical networks to minimize the number of transceivers
% published in:
%
% [1]: Konda, V.R., Chow, T.Y.: Algorithm for Traffic Grooming
%      in Optical Networks to Minimize the Number of Transceivers.
%      IEEE Workshop on High Performance Switching and Routing, (2001),
%      pp. 218-221 
%
% Usage: [exitMsg exitFlag netState] = HPSR2001_minTransceivers(traff_trafficMatrix, phys, parameters)
%
% Abstract: This function solves the virtual topology desing problem by
% means of the algorithm proposed in the reference [1].
%
% This function solves the four classic subproblems into which the
% Virtual Topology Design Problem is possible to decompose:
%
% 1) Virtual Topology Subproblem
% 2) Lightpath Routing Subproblem
% 3) Wavelength Assignment Subproblem
% 4) Traffic Routing (over the Virtual Topology) Subproblem
%
% The algorithm proposed by Konda in [1] solves the subproblem 1: it
% selects the node pairs candidate to establish the lightpaths by
% minimizing the number of used transceivers. The subproblem 1 does not
% consider the lightpath routing between the previous node pairs
% (subproblem 2) neither the wavelength assignment (subproblem 3). The
% subproblem 2 and 3 are also jointly called RWA (Routing and Wavelength
% Assignmen) problem. The RWA is implemented by the function
% 'libraryRWA_kShortestRWA', which search the k shortest routes valid to
% establish a lightpath between a pair of node. (For more information see
% this function). Finally the subproblem 4) is solved by usinf the function
% 'libraryFR_optimalFlowAssignment', which implements a LP formulation to
% address the Optimal Flow Assignment Problem.
%
% Arguments:
% o	In: 
%   traff_trafficMatrix(NxN): Average traffic flow offered between node
%    pairs. The Traffic Matrix is a two-dimensional matrix with N (N:
%    number of nodes) rows and N columns. An entry(s,d) means the average
%    traffic flow from node 's' to node 'd', expressed in Gbps. The main
%    diagonal is full of 0s.
%
%  . phys: Phys Structure. More information about netState in section 
%    "Structure of phys variable" from Help.
%
%  . parameters: algorithm parameters. For this algorithm, it is the number
%    of 'k' shortest paths computed as candidate routes for a each
%    lightpath.).
%   
% o	Out:
%  . exitMsg: Exit Message. If the virtual topology design was successful
%    and all the traffic flows were routed, exitMsg='OK!!!'(exitFlag == 0); 
%    otherwise, exitMsg indicates 'No feasible solution found' (exitFlag == 1) 
%    or other situation (exitFlag >= 2).
%
%  . exitFlag: 
%           0, if it is possible to design the virtual topology and route
%              all the traffic flows over the virtual topology.
%           1, if it is not possible to find a feasible solution
%           2, if the flow routing is not optimal
%
%  . netState: NetState Structure. More information about netState in
%    section "Structure of netState variable" from Help. 
%
%


function [exitMsg exitFlag netState] = HPSR2001_minTransceivers(traff_trafficMatrix, phys, parameters)

exitMsg = 'OK';
exitFlag = 0;
netState = struct ('flowRoutingMatrix' , [] ,  'lightpathRoutingMatrix' , [] , 'lightpathTable' , [] , 'flowTable' , [] , 'numberOfOccupiedTWCs' , [],'numberOfOccupiedTxs' , [],'numberOfOccupiedRxs' , []);

numberOfNodes = max(max(phys.linkTable));
numberOfLinks = size(phys.linkTable,1);

lightpathTable=[];
lightpathRoutingMatrix=[];
netState.numberOfOccupiedTxs = zeros (numberOfNodes,1);
netState.numberOfOccupiedRxs = zeros (numberOfNodes,1);
netState.numberOfOccupiedTWCs = zeros (numberOfNodes,1); 
%%%%%%%%%%%%%%%%%%%%%%%%0. PARAMETERS PROCESSING %%%%%%%%%%%%%%%%%%%%%%%%%
%We define the default parameter values 
cellOfDefaultParameters = {['k'], ['1'], ['integer'], ['[1,inf)']};      
                       
%We process the string of the algorithm parameters                       
[exitFlagOfProcessParam exitMsgOfProcessParam userParam] = ...
    processAlgorithmParameters (parameters, cellOfDefaultParameters);   
if (exitFlagOfProcessParam ~= 0), error(['>>> Bad Algorithm Parameters:', exitMsgOfProcessParam]); end

NUMBERLINKS = size (phys.linkTable,1);
NUMBERNODES = size (traff_trafficMatrix,1);

%1) Virtual Topology Selection Subproblem: Select the node pairs candidate to
%establish the lightpaths.

%%% 1.a) Create the initial structure of lightpaths according to [1]: All the
%%% node pairs connected by the sufficient lightpaths to route all the
%%% traffic demand.
netState_lightpathTable = zeros (0,2);
capacitiesMatrix = zeros (NUMBERNODES,NUMBERNODES);
lightpathSurplusCapacityVector = zeros (0,0);

for ingressNode=1:NUMBERNODES
    for egressNode=1:NUMBERNODES
        offeredTraffic = traff_trafficMatrix(ingressNode,egressNode);        
        pendingOfferedTraffic = offeredTraffic;
        numberExistingLPsBeforeAdding = size (netState_lightpathTable,1);
        for numberLPs=1:ceil(offeredTraffic/phys.lightpathCapacity)
            netState_lightpathTable = [netState_lightpathTable ; ingressNode egressNode];
            carriedTrafficInThisLP = min (pendingOfferedTraffic,phys.lightpathCapacity);
            lightpathSurplusCapacityVector = [lightpathSurplusCapacityVector ; (phys.lightpathCapacity-carriedTrafficInThisLP)];
            pendingOfferedTraffic = pendingOfferedTraffic - carriedTrafficInThisLP;
            capacitiesMatrix (ingressNode , egressNode) = capacitiesMatrix (ingressNode , egressNode) + phys.lightpathCapacity;
        end
    end
end

weMustStop = 0; 

%%% 1.b) We tear down a lightpath per iteration and we try to route the
%%% traffic carried by the removed LP throug the remaining virtual
%%% topology according to [1].
while (weMustStop ~= 1) && (sum(lightpathSurplusCapacityVector) >=  phys.lightpathCapacity)
    
    minimumSurplusCapacityConsumedTotal = inf;
    
    lpIDWithMinimumLoadTotal = [];

    for ingressNode=1:NUMBERNODES
        for egressNode=1:NUMBERNODES

            if (ingressNode==egressNode)
                continue
            end
            % Obtain the set of LPs between these node pairs
            setOfLPsBetweenThisNodePair_1 = find (netState_lightpathTable(:,1) == ingressNode);
            setOfLPsBetweenThisNodePair_2 = find (netState_lightpathTable(:,2) == egressNode);
            setOfLPsBetweenThisNodePair = intersect (setOfLPsBetweenThisNodePair_1,setOfLPsBetweenThisNodePair_2);
            if(isempty(setOfLPsBetweenThisNodePair))
                continue
            end
            % Obtain the LP between the node pair least loaded
            [maximumSurplusCapacityThisIteration,lightpathPositionInSet] = max (lightpathSurplusCapacityVector (setOfLPsBetweenThisNodePair));
            minimumCarriedTrafficThisIteration = phys.lightpathCapacity - maximumSurplusCapacityThisIteration;
            lpIDWithMinimumLoadThisIteration = setOfLPsBetweenThisNodePair (lightpathPositionInSet);
            % Try to route the capacity of the LP candidate in the surplus
            % topology
            linkCostVector = ones (size (netState_lightpathTable,1),1);
            lightpathSurplusCapacityVector (lpIDWithMinimumLoadThisIteration) = 0;
            [exitFlag_MCF , flowRoutingVectorThisIteration , minimumSurplusCapacityConsumedThisIteration] = libraryGraph_minCostFlow (ingressNode , egressNode , minimumCarriedTrafficThisIteration , netState_lightpathTable , linkCostVector , lightpathSurplusCapacityVector);

            lightpathSurplusCapacityVector (lpIDWithMinimumLoadThisIteration) = phys.lightpathCapacity -minimumCarriedTrafficThisIteration;

            if (exitFlag_MCF == 0) && (minimumSurplusCapacityConsumedThisIteration < minimumSurplusCapacityConsumedTotal)
                flowRoutingVectorTotal = flowRoutingVectorThisIteration;
                lpIDWithMinimumLoadTotal = lpIDWithMinimumLoadThisIteration;
                minimumSurplusCapacityConsumedTotal = minimumSurplusCapacityConsumedThisIteration;
           end
        end
    end
    if (isempty (lpIDWithMinimumLoadTotal))
        weMustStop = 1;
    else
        netState_lightpathTable (lpIDWithMinimumLoadTotal,:) = [];
        lightpathSurplusCapacityVector = lightpathSurplusCapacityVector - flowRoutingVectorTotal';
        lightpathSurplusCapacityVector (lpIDWithMinimumLoadTotal) = [];
    end
end

%2) Routing and Wavelength Assignment Subproblem: Determines the physical
%fibre links which each lightpath consists of, and assigns a wavelength to
%each fibre link.
% We must calculate the lightpathRoutingMatrix of the LPs to be created
% (the remaining)
serialNumber=1;
for lpId = 1:size(netState_lightpathTable,1)
    
    lpIngressNode = netState_lightpathTable(lpId,1);
    lpEgressNode = netState_lightpathTable(lpId,2);
    if(netState.numberOfOccupiedTxs(lpIngressNode)<phys.numberTxPerNode(lpIngressNode) & netState.numberOfOccupiedRxs(lpEgressNode)<phys.numberRxPerNode(lpEgressNode) )
        
        [exitFlagRWA sequenceOfSPFLinkIds sequenceOfSPFNodeIds sequenceOfWavelengths occupiedTWCsPerNode newLightpathRoutingEntry] = ...
        libraryRWA_kShortestRWA(userParam.k, lpIngressNode, lpEgressNode, phys, netState);

        % update the netState if the LP is set up
        if (exitFlagRWA == 0)
            netState.lightpathTable=[netState.lightpathTable;serialNumber lpIngressNode lpEgressNode];
            serialNumber=serialNumber+1;
            netState.lightpathRoutingMatrix = [netState.lightpathRoutingMatrix;newLightpathRoutingEntry];
            netState.numberOfOccupiedTxs(lpIngressNode) = netState.numberOfOccupiedTxs(lpIngressNode) + 1;
            netState.numberOfOccupiedRxs(lpEgressNode) = netState.numberOfOccupiedRxs(lpEgressNode) + 1;  
            netState.numberOfOccupiedTWCs = netState.numberOfOccupiedTWCs + occupiedTWCsPerNode;
        else
            exitMsg = 'Not possible to establish a lightpath'; % Not possible to establish a lightpath;
            exitFlag = 1; 
            return;
        end
    end
end

%3) Traffic Flow Routing (over the virtual topology) Subproblem.
% We must calculate the flowRoutingMatrix on top of the virtual topology.
% We make it using the optimalFA. If everythign is OK, the capacities MUST
% be enough to carry the flows

[exitFlag_FR flowTable netState.flowRoutingMatrix cost] = libraryFR_optimalFlowAssignment(traff_trafficMatrix, netState.lightpathTable(:,2:3) , phys.lightpathCapacity*ones(1,size(netState.lightpathTable,1)) , ones(1,size(netState.lightpathTable,1)) , '');
netState.flowTable = [(1:size(flowTable,1))' flowTable -1*ones(size(flowTable,1),3)];% Add the serial number of the flows at fist column.
    
switch exitFlag_FR,
    case 0,
        exitFlag = 0;
        exitMsg = 'OK';
    case 1,
        exitFlag = 2;
        exitMsg = 'Some traffic flow was not routed or no optimal flow routing found';
    case 2
        exitFlag = 1;
        exitMsg = 'No feasible solution found';
    otherwise
        exitFlag = 1;
        exitMsg = 'No feasible solution found';
end