libraryFR_optimalFlowAssignment(traff_trafficMatrix, linkTable , linkCapacities , linkCostPerGbps, varargin)

%%% libraryFR_optimalFlowAssignment
%
% Usage: [exitFlag flowTable flowRoutingMatrix cost] =
% libraryFR_optimalFlowAssignment(traff_trafficMatrix, linkTable , linkCapacities , linkCostPerGbps)
%
% Abstract: This function solves the Traffic Flow Routing over a network
% topology. It implements an optimal Flow Assignment problem. 
% The problem is established as linear programming problem, 
% and solved with the linprog solver
%  The function provides the solution which minimizes a given cost function. 
% The cost function employed is the sum of the costs in each link. 
% The cost in each link is give by the product of the Gbps 
% carried and cost per Gbps in the link.
%
% 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.
%
% .  linkTable(M,2): M-by-2 integer matrix. Each row is a link in the graph. 
%   first column is the ingress node (1...N), second the egress node
%   (1...N)
%
% .  linkCapacities(M,1): Capacity in Gbps of each link. This is the
% maximum traffic that the link can carry.
%
% .  linkCostPerGbps(M,1): Cost of carrying 1 Gbps in each link. If the
% vector is constant, the average number of hops is minimized. If the
% vector contains a measure of link distance, the average distance traversed 
% by the carried traffic is minimized
%
% o	Out:
%  . exitFlag: 
%           0, if LINPROG converged to a solution X.(It means that all the flows
%           can be routed)
%           1, if maximum number of iterations reached.
%           2, if the traffic flow routing is not FEASIBLE.
%
%  . flowTable(F,3): F-by-3 matrix. First column de ingress node, second
%    the egress node, third the traffic offered in Gbps.
%
%  . flowRoutingMatrix (F,L): F-by-L integer matrix where F is the number of 
%    flows and L is the number of lightpaths. Each row is a flow 'f' and each 
%    column is a lightpath 'l'. If a flow 'f' uses a lightpath 'l', the
%    entry (f,l) is equal to the value of the flow 'f' carried by the
%    lightpath 'f'. If no lightpath is used by the flow 'f', the entry is
%    equal to '0'. If a row is a row of zeros, the the traffic flows was
%    not routed
%
%   cost: value of the cost function at optimum  

function [exitFlag flowTable flowRoutingMatrix cost] = libraryFR_optimalFlowAssignment(traff_trafficMatrix, linkTable , linkCapacities , linkCostPerGbps, varargin)

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

% Create the flowTable from the trafficMatrix
numberOfFlows = nnz(traff_trafficMatrix);
flowTable = zeros (numberOfFlows,3);
contFlow = 0;
for ingressNode=1:NUMBERNODES
    for egressNode=1:NUMBERNODES
        if (ingressNode ~= egressNode) && (traff_trafficMatrix(ingressNode,egressNode) ~= 0),
            contFlow = contFlow + 1;          
            flowTable (contFlow,:) = [ingressNode egressNode traff_trafficMatrix(ingressNode,egressNode)];            
        end
    end
end
NUMBERFLOWS = contFlow;
NUMBERDECISIONVARIABLES = NUMBERFLOWS * NUMBERLINKS;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Non-negativity constraints
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = -eye ( NUMBERDECISIONVARIABLES);
b = zeros ( NUMBERDECISIONVARIABLES , 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Flow conservation constraints (one per flow and node)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for nodeID=1:NUMBERNODES
    % inputLinks are the incoming edges
    inputLinks = find (linkTable (:,2) == nodeID);
    % outputLinks are the outgoing edges
    outputLinks = find (linkTable (:,1) == nodeID);
   
    for flowId=1:NUMBERFLOWS
        flowOriginID = flowTable (flowId , 1);
        flowDestID = flowTable (flowId , 2);    
        flowGbps = flowTable (flowId , 3);
      
        constraint = zeros (1 ,  NUMBERDECISIONVARIABLES);
        
        for contOutputLink=1:length (outputLinks)
            linkId = outputLinks (contOutputLink);     
            columnA = columnMatrixA (linkId , flowId , NUMBERFLOWS);
            constraint (columnA) = 1;
        end
        for contInputLink=1:length (inputLinks)
            linkId = inputLinks (contInputLink);
            columnA = columnMatrixA (linkId , flowId , NUMBERFLOWS);
            constraint (columnA) = -1;
        end            
            
      
        A = [A ; constraint ; -constraint];
        
        if (nodeID == flowOriginID)
            b  = [b ; flowGbps ; -flowGbps]; % node is ingress of the flow
        elseif (nodeID == flowDestID)            
            b = [b ; -flowGbps ; flowGbps]; % node is ingress of the flow
        else
            b = [b ; 0 ; 0]; % node is intermediate to the flow
        end
    end    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constraint carried traffic in each link limited by link capacity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for linkId = 1:NUMBERLINKS
    constraint = zeros (1, NUMBERDECISIONVARIABLES);
    for flowId=1:NUMBERFLOWS
        columnA = columnMatrixA (linkId , flowId , NUMBERFLOWS);
        constraint (columnA) = 1;
    end
    
    A = [A ; constraint];
    b = [b ; linkCapacities(linkId)]; 
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Cost function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
costFunction = zeros (1, NUMBERDECISIONVARIABLES);
for linkId = 1:NUMBERLINKS
    for flowId=1:NUMBERFLOWS
        columnA = columnMatrixA (linkId , flowId , NUMBERFLOWS);
        costFunction (columnA) = linkCostPerGbps (linkId);
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% LINPROG
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

options=optimset('display','off');

[flowRoutingVector,cost,exitFlagLP] = linprog (costFunction',A,b,[],[],[],[],[],options);
%  . exitFlagLP: The value returned by the linprog solver
%      1  LINPROG converged to a solution X. (Means the flows can be carried)
%       0  Maximum number of iterations reached.
%      -2  No feasible point found.
%      -3  Problem is unbounded.
%      -4  NaN value encountered during execution of algorithm.
%      -5  Both primal and dual problems are infeasible.
%      -7  Search direction became too small; no further progress can be
%      made. 
switch exitFlagLP
    case 1 
        exitFlag = 0;
    case 0
        exitFlag = 1;
    case -2
        exitFlag = 2;
        flowTable =[];
        flowRoutingMatrix=[];
        return
    case -3
        exitFlag = 2;
        flowTable =[];
        flowRoutingMatrix=[];
        return
    case -4
        exitFlag=2;
        flowTable =[];
        flowRoutingMatrix=[];
        return
    case -5
        exitFlag=2;
        flowTable =[];
        flowRoutingMatrix=[];
        return
    case -7
        exitFlag=2;
        flowTable =[];
        flowRoutingMatrix=[];
        return
end

for cont=1:length(flowRoutingVector)
    if (flowRoutingVector(cont)<1E-6)
        flowRoutingVector(cont) = 0;
    end
end
    

flowRoutingMatrix = zeros (NUMBERFLOWS , NUMBERLINKS);
for flowId=1:NUMBERFLOWS
    for linkId=1:NUMBERLINKS
        flowRoutingMatrix(flowId,linkId) = flowRoutingVector (columnMatrixA (linkId , flowId , NUMBERFLOWS));
    end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Receives a flowId (1,2,...) and a linkId (1,2,...) and return the index
% of the column in the constraints matrix for the decision variable
% associated to that flow and link
function columnID = columnMatrixA (linkId , flowId , NUMBERFLOWS)

    columnID = flowId  + NUMBERFLOWS * (linkId-1);