libraryGraph_kshortestPath (k,linkTable,ingressNode,egressNode)

% libraryGraph_kshortestPath
%
% Usage: [cellOfKPaths] = libraryGraph_kshortestPath
% (k, linkTable, ingressNode, egressNode)
%
% Abstract: This function computes the "K" different shortest paths between
% a pair of nodes (ingress node, egress node). A shortest path is defined
% as the path with the minimum total weight (or cost), where the total
% weight is the sum of the weights of the links which the path consist of. 
% If there are only "X" shortest paths between ingress node and egress
% node, where "X" < "K" (K is the demanded number of shortest paths), the
% function returns the "X" found shortest paths, and empty matrices up to
% reach the number "K". The implementation of this function is based on the
% K-shortest path algorithm proposed in :
% [1]	E.Q. Martins, An algorithm to ranking paths that may contains
% cycles, European Journal of Operational Research, vol. 18, pp. 123-130,
% 1984
%
% Arguments:
% o	In: 
%    k: Integer number. Demanded number of shortest paths between ingress
%     node and egress node.
%
%    linkTable(L,3): Matrix of L physical links by 3. Each row is a
%     physical link 'l', where the first and second columns of each row are
%     the nodeID (integer number) of the origin node 'x' and the nodeID
%     (integer number) of the destination node 'y' of this link 'l'
%     respectively. The third column is the link weight of the link 'l'.
%
%    ingressNode: Origin Node of the path.
%
%   . egressNode: Destination Node of the path.
%   
% o	Out:
%    cellOfKPaths(K,1): K-by-3 cell of vectors. Each row is one of the "K"
%     shortest paths. The rows (shortest paths) are sorted in decreasing
%     weight (or cost). Each row consist of the three columns:
%       - 1st) sequenceOfSPFLinkIds: Row vector with the sorted sequence of
%         IDs of the links which the shortest path between ingress and
%         egress consist of. The links IDs are sorted from the ingress node
%         to the egress node.
%       - 2nd) sequenceOfSPFNodeIds: Row vector with the sorted sequence of
%         IDs of the nodes which the shortest path  between ingress
%         and egress consist of. The nodes IDs are sorted from the ingress
%         node to the egress node.
%       - 3rd) totalCost: It is the total weight of the shortest path found
%         between ingress and egress. It is computed as the sum of the
%         weights of the links which the shortest path consist of.
%
%


function [cellOfKPaths] = libraryGraph_kshortestPath (k,linkTable,ingressNode,egressNode)
    
    cellOfKPaths = cell(k,3);

    % first we compute for k=1
    [spfk1_sequenceOfLinkIds , spfk1_sequenceOfNodeIds , totalCost] = libraryGraph_shortestPath (linkTable , ingressNode , egressNode);
       
    if (length(spfk1_sequenceOfLinkIds) == 0)
        disp ('There is no path');
        return
    end

    cellOfKPaths{1,1} = spfk1_sequenceOfLinkIds;
    cellOfKPaths{1,2} = spfk1_sequenceOfNodeIds;
    cellOfKPaths{1,3} = totalCost;

    if k == 1
        return;
    end

    originalGraph_maxNodeNumber = max(max(linkTable(:,1:2)));
    originalGraph_maxLinkNumber = size(linkTable,1);

    % coord i is the orginal node which this node is mirror of
    % coord i belongs to [1...maxNodeNumberK1]
    % there are as many coordinates as different nodes

    previousSPF_sequenceOfLinkIds = spfk1_sequenceOfLinkIds;
    previousSPF_sequenceOfNodeIds = spfk1_sequenceOfNodeIds;
    previousGraph_physLinkTable = linkTable;
    previousGraph_nodeMirrors = -1 * ones (1,originalGraph_maxNodeNumber); % the nodes that are mirror nodes from itself, are put to -1
    previousGraph_linkMirrors = -1 * ones (1,originalGraph_maxLinkNumber); % the links that are mirror nodes from itself, are put to -1

    for currentK=2:k  
        
        [previousSPF_intermediateLinkIdsPathToRemove,previousSPF_intermediateNodeIdsPathToRemove]  = sequenceOfIntermediateNodesAndLinks (previousSPF_sequenceOfLinkIds , previousSPF_sequenceOfNodeIds);
        
        currentGraph_physLinkTable = previousGraph_physLinkTable;
        currentGraph_nodeMirrors = previousGraph_nodeMirrors;% vector. The component i is the id of the node which the node i is mirror of. It is -1 if the node i is no mirror from anyone
        currentGraph_linkMirrors = previousGraph_linkMirrors;% vector. The component i is the id of the link which the link i is mirror of. It is -1 if the link i is no mirror from anyone
                                                             
        % step 1 and 2: create the mirror nodes to the intermediate nodes, and the mirror intermediate links 
        % step 3: 
        for contIntermediateNode=1:length(previousSPF_intermediateNodeIdsPathToRemove)
 
            originalNodeId = previousSPF_intermediateNodeIdsPathToRemove (contIntermediateNode);
            mirrorNodeId = length (currentGraph_nodeMirrors) + 1;
            currentGraph_nodeMirrors = [currentGraph_nodeMirrors originalNodeId];

            if contIntermediateNode > 1
                % mirror intermediate links of the path
                previousGraph_linkIdToMirror = previousSPF_intermediateLinkIdsPathToRemove (contIntermediateNode - 1);            
                currentGraph_sourceNodeLinkToAdd = length(currentGraph_nodeMirrors)-1;
                currentGraph_targetNodeLinkToAdd = length(currentGraph_nodeMirrors);
                currentGraph_weightLinkToAdd = currentGraph_physLinkTable (previousGraph_linkIdToMirror , 3);

                currentGraph_physLinkTable = [currentGraph_physLinkTable ; currentGraph_sourceNodeLinkToAdd currentGraph_targetNodeLinkToAdd currentGraph_weightLinkToAdd];
                currentGraph_linkMirrors = [currentGraph_linkMirrors previousGraph_linkIdToMirror];
            end

            % step 3: change the arcs that direct to the original internal nodes,
            % to direct to their mirrors (except the arcs of the path)
            linkIdsTargetedToOriginalIntermediateNodeId = find (previousGraph_physLinkTable(:,2) == originalNodeId);
            for contLinkIdToChange=1:length(linkIdsTargetedToOriginalIntermediateNodeId)
                linkIdToChange=linkIdsTargetedToOriginalIntermediateNodeId(contLinkIdToChange);
                %if the node does not belong to the path to delete, then we
                %link it.
                if (isempty (find (previousSPF_sequenceOfLinkIds==linkIdToChange)))
                    currentGraph_physLinkTable (linkIdToChange,2) = mirrorNodeId;
                end
            end
        end

       % step 4: move an arc
       if (length (previousSPF_intermediateLinkIdsPathToRemove) == 0)
           % If the path to delete was a single-hop path, there is no
           % intermediate links neither intermediate nodes. For this
           % reason, we must create a imaginary mirror node so that the
           % step 4 can works correct.
           currentGraph_nodeMirrors = [currentGraph_nodeMirrors -1];
       end
        lastLinkIdOfPath = previousSPF_sequenceOfLinkIds (end);
        currentGraph_physLinkTable (lastLinkIdOfPath , 1) = length (currentGraph_nodeMirrors); % enlazo ultimo link del camino a borrar, tal que ahora provenga del ultimo nodo espejo creado
       
        % calculate the spf in the new graph
        [spfkcurrent_sequenceOfLinkIds , spfkcurrent_sequenceOfNodeIds , totalCost] = libraryGraph_shortestPath (currentGraph_physLinkTable , ingressNode , egressNode);    

        if isempty (spfkcurrent_sequenceOfLinkIds)
            return;
        end
               
        % convert the mirror node ids found in the links and nodes, to their
        % original values
        spfToReturn_sequenceOfLinkIds = convertMirrorIds (spfkcurrent_sequenceOfLinkIds , currentGraph_linkMirrors , originalGraph_maxLinkNumber);
        spfToReturn_sequenceOfNodeIds = convertMirrorIds (spfkcurrent_sequenceOfNodeIds , currentGraph_nodeMirrors , originalGraph_maxNodeNumber);

        % add the new spf path to the cell
        cellOfKPaths{currentK,1} = spfToReturn_sequenceOfLinkIds;
        cellOfKPaths{currentK,2} = spfToReturn_sequenceOfNodeIds;
        cellOfKPaths{currentK,3} = totalCost;

        exitFlag = checkPath (spfToReturn_sequenceOfLinkIds , spfToReturn_sequenceOfNodeIds , linkTable , ingressNode , egressNode);
        if (exitFlag == -1) ||(cellOfKPaths{currentK-1,3}>totalCost)
            disp ('Error: checkpath: no son un camino de origen a destino o el coste es mejor y no peor');
            spfToReturn_sequenceOfLinkIds
            spfToReturn_sequenceOfNodeIds
            pause
        end
        

        previousSPF_sequenceOfLinkIds = spfkcurrent_sequenceOfLinkIds;
        previousSPF_sequenceOfNodeIds = spfkcurrent_sequenceOfNodeIds;
        
        previousGraph_physLinkTable = currentGraph_physLinkTable;
        previousGraph_nodeMirrors = currentGraph_nodeMirrors;
        previousGraph_linkMirrors = currentGraph_linkMirrors;
    end
end


function [convertedSequence] = convertMirrorIds (sequence , mirrorTable , originalGraph_maxNodeNumber)
    convertedSequence = sequence;
    for contId=1:length (sequence)
        while (convertedSequence (contId) > originalGraph_maxNodeNumber)
            convertedSequence (contId) = mirrorTable (convertedSequence (contId));
        end
    end
end


function [intermLinkIdSequence,intermNodeIdSequence] = sequenceOfIntermediateNodesAndLinks (sequenceOfLinkIds , sequenceOfNodeIds)
   if length(sequenceOfNodeIds) <= 2
        intermNodeIdSequence = [];
        intermLinkIdSequence = [];
        return
   elseif length(sequenceOfNodeIds) == 3
        intermNodeIdSequence = sequenceOfNodeIds (2);
        intermLinkIdSequence = [];
   else
        intermNodeIdSequence = sequenceOfNodeIds (2:end-1);
        intermLinkIdSequence = sequenceOfLinkIds (2:end-1);
   end       
end

function [exitFlag] = checkPath (sequenceLinkIds , sequenceNodeIds , linkTable , ingressNode , egressNode)
exitFlag = -1;
if (sequenceNodeIds(1) ~= ingressNode) || (sequenceNodeIds(end) ~= egressNode)
    return;
end

for contLinkId=1:length(sequenceLinkIds)
    linkId = sequenceLinkIds (contLinkId);
    linkSourceNode = linkTable(linkId,1);
    linkTargetNode = linkTable(linkId,2);
    % the origin node and the destinatio node are the nodes that
    % corresponds to sequenceNodeId
    if (linkSourceNode ~= sequenceNodeIds (contLinkId)) || (linkTargetNode ~= sequenceNodeIds (contLinkId+1))
        return;
    end
end

exitFlag = 0;

end