K shortest paths in a graph represented by a sparse matrix (Yen's algorithm)

by

 

01 Mar 2012 (Updated )

Determine the K shortest paths from node S to node T.

graphkshortestpaths.m
function [ DIST, PATH ] = graphkshortestpaths( G, S, T, K )
%
% [ DIST, PATH ] = graphkshortestpaths( G, S, T, K ) determines the K shortest paths from node S
% to node T. weights of the edges are all positive entries in the n-by-n adjacency matrix
% represented by the sparse matrix G. DIST are the K distances from S to T; PATH is a cell array
% with the K shortest paths themselves.
%
% the shortest path algorithm used is Dijkstra's algorithm (graphshortestpath).
%
% **Please note that the algorithm implemented here is an undirected version of Yen's algorithm**
%
% - Yen, JY. Finding the k shortest loopless paths in a network; Management Science 17(11):712-6.
%
% 03/01/2013: I would like to thank Oskar Blom Gransson for helping me find a bug in the previous version.

% find A^1
[ DIST( 1 ), PATH{1} ] = graphshortestpath( G, S, T );

candidate_paths = { }; % list of candidate paths
candidate_dists = [ ]; % distance of each candidate path

% find A^2 ... A^K
for k = 2:K
	k_G = G; % version of G used in this iteration (some edges will be removed)

	% for each node travelled in A^{k-1}
	for i = 1:( length( PATH{k-1} ) - 1 )
		i_node = PATH{k-1}( i );

		% iterate over all previous paths and examine if 1..i overlaps with A^{k-1}
		for j = 1:k-1
			if( length( PATH{j} ) >= i & ( all( PATH{j}( 1:i ) == PATH{k-1}( 1:i ) ) ) )
				% it does; remove the following edge that appears in A^j
				j_next_node = PATH{j}( i+1 );
				k_G( i_node, j_next_node ) = 0; k_G( j_next_node, i_node ) = 0;
			end
		end

		% calculate shortest path from i to T
		[ dist_i_t, path_i_t ] = graphshortestpath( k_G, i_node, T );
		% if path exists, concatenate with 1..i-1 and add to candidates list
		if( dist_i_t < Inf )
			path_1_i_t = [ PATH{k-1}( 1:i-1 ) path_i_t ];
			dist_1_i_t = graphpathdistance( G, path_1_i_t ); % we can safely use G- removed
									 % edges will not appear
			% add resulting path to candidates list
			candidate_paths{end+1} = path_1_i_t;
			candidate_dists( end+1 ) = dist_1_i_t;
		end
	end

	% no candidates; all shortest paths found
	if( isempty( candidate_dists ) )
		return
	end

	% take shortest path from candidates list as kth path
	[ y, i ] = sort( candidate_dists );
	DIST( k ) = candidate_dists( i( 1 ) );
	PATH{k} = candidate_paths{i( 1 )};
	% remove shortest path (and all of its copies) from the candidates list
	remove_indices = [];
	for idx = 1:length( candidate_paths )
		if( length( PATH{k} ) == length( candidate_paths{idx} ) & all( PATH{k} == candidate_paths{idx} ) )
			remove_indices( end+1 ) = idx;
		end
	end
	candidate_dists( remove_indices ) = [];
	candidate_paths( remove_indices ) = [];
end


	function DIST = graphpathdistance( G, PATH )
	%
	% DIST = graphpathdistance( PATH ) calculates the distance travelled by PATH in graph G.

	% convert PATH into edges
	edges = sub2ind( size( G ), PATH( 1:end-1 ), PATH( 2:end ) );
	% sum weights over edges
	DIST = full( sum( G( edges ) ) );


Contact us