Daniel@0: function [D,P] = dijk(A,s,t)
Daniel@0: %DIJK Shortest paths from nodes 's' to nodes 't' using Dijkstra algorithm.
Daniel@0: % [D,p] = dijk(A,s,t)
Daniel@0: %     A = n x n node-node weighted adjacency matrix of arc lengths
Daniel@0: %         (Note: A(i,j) = 0   => Arc (i,j) does not exist;
Daniel@0: %                A(i,j) = NaN => Arc (i,j) exists with 0 weight)
Daniel@0: %     s = FROM node indices
Daniel@0: %       = [] (default), paths from all nodes
Daniel@0: %     t = TO node indices
Daniel@0: %       = [] (default), paths to all nodes
Daniel@0: %     D = |s| x |t| matrix of shortest path distances from 's' to 't'
Daniel@0: %       = [D(i,j)], where D(i,j) = distance from node 'i' to node 'j' 
Daniel@0: %     P = |s| x n matrix of predecessor indices, where P(i,j) is the
Daniel@0: %         index of the predecessor to node 'j' on the path from 's(i)' to 'j'
Daniel@0: %         (use PRED2PATH to convert P to paths)
Daniel@0: %       = path from 's' to 't', if |s| = |t| = 1
Daniel@0: %
Daniel@0: %  (If A is a triangular matrix, then computationally intensive node
Daniel@0: %   selection step not needed since graph is acyclic (triangularity is a 
Daniel@0: %   sufficient, but not a necessary, condition for a graph to be acyclic)
Daniel@0: %   and A can have non-negative elements)
Daniel@0: %
Daniel@0: %  (If |s| >> |t|, then DIJK is faster if DIJK(A',t,s) used, where D is now
Daniel@0: %   transposed and P now represents successor indices)
Daniel@0: %
Daniel@0: %  (Based on Fig. 4.6 in Ahuja, Magnanti, and Orlin, Network Flows,
Daniel@0: %   Prentice-Hall, 1993, p. 109.)
Daniel@0: 
Daniel@0: % Copyright (c) 1998-2001 by Michael G. Kay
Daniel@0: % Matlog Version 5 22-Aug-2001
Daniel@0: 
Daniel@0: % Input Error Checking ******************************************************
Daniel@0: error(nargchk(1,3,nargin));
Daniel@0: 
Daniel@0: [n,cA] = size(A);
Daniel@0: 
Daniel@0: if nargin < 2 | isempty(s), s = (1:n)'; else s = s(:); end
Daniel@0: if nargin < 3 | isempty(t), t = (1:n)'; else t = t(:); end
Daniel@0: 
Daniel@0: if ~any(any(tril(A) ~= 0))       % A is upper triangular
Daniel@0:    isAcyclic = 1;
Daniel@0: elseif ~any(any(triu(A) ~= 0))   % A is lower triangular
Daniel@0:    isAcyclic = 2;
Daniel@0: else                             % Graph may not be acyclic
Daniel@0:    isAcyclic = 0;
Daniel@0: end
Daniel@0: 
Daniel@0: if n ~= cA
Daniel@0:    error('A must be a square matrix');
Daniel@0: elseif ~isAcyclic & any(any(A < 0))
Daniel@0:    error('A must be non-negative');
Daniel@0: elseif any(s < 1 | s > n)
Daniel@0:    error(['''s'' must be an integer between 1 and ',num2str(n)]);
Daniel@0: elseif any(t < 1 | t > n)
Daniel@0:    error(['''t'' must be an integer between 1 and ',num2str(n)]);
Daniel@0: end
Daniel@0: % End (Input Error Checking) ************************************************
Daniel@0: 
Daniel@0: A = A';    % Use transpose to speed-up FIND for sparse A
Daniel@0: 
Daniel@0: D = zeros(length(s),length(t));
Daniel@0: if nargout > 1, P = zeros(length(s),n); end
Daniel@0: 
Daniel@0: for i = 1:length(s)
Daniel@0:    j = s(i);
Daniel@0:    
Daniel@0:    Di = Inf*ones(n,1); Di(j) = 0;
Daniel@0:    
Daniel@0:    isLab = logical(zeros(length(t),1));
Daniel@0:    if isAcyclic ==  1
Daniel@0:       nLab = j - 1;
Daniel@0:    elseif isAcyclic == 2
Daniel@0:       nLab = n - j;
Daniel@0:    else
Daniel@0:       nLab = 0;
Daniel@0:       UnLab = 1:n;
Daniel@0:       isUnLab = logical(ones(n,1));
Daniel@0:    end
Daniel@0:    
Daniel@0:    while nLab < n & ~all(isLab)
Daniel@0:       if isAcyclic
Daniel@0:          Dj = Di(j);
Daniel@0:       else	% Node selection
Daniel@0:          [Dj,jj] = min(Di(isUnLab));
Daniel@0:          j = UnLab(jj);
Daniel@0:          UnLab(jj) = [];
Daniel@0:          isUnLab(j) = 0;
Daniel@0:       end
Daniel@0:       
Daniel@0:       nLab = nLab + 1;
Daniel@0:       if length(t) < n, isLab = isLab | (j == t); end
Daniel@0:       
Daniel@0:       [jA,kA,Aj] = find(A(:,j));
Daniel@0:       Aj(isnan(Aj)) = 0;
Daniel@0:             
Daniel@0:       if isempty(Aj), Dk = Inf; else Dk = Dj + Aj; end
Daniel@0:       
Daniel@0:       if nargout > 1, P(i,jA(Dk < Di(jA))) = j; end
Daniel@0:       Di(jA) = min(Di(jA),Dk);
Daniel@0:       
Daniel@0:       if isAcyclic == 1       % Increment node index for upper triangular A
Daniel@0:          j = j + 1;
Daniel@0:       elseif isAcyclic == 2   % Decrement node index for lower triangular A
Daniel@0:          j = j - 1;
Daniel@0:       end
Daniel@0:    end
Daniel@0:    D(i,:) = Di(t)';
Daniel@0: end
Daniel@0: 
Daniel@0: if nargout > 1 & length(s) == 1 & length(t) == 1
Daniel@0:    P = pred2path(P,s,t);
Daniel@0: end