Dijkstra's Algorithm - how to calculate the shortest path (having a state matrix)

Hello, I am slowly learning Matlab and now have assignment with shortest path using Dijkstra's algorithm.

I got the following code, which is saved as "dijkstra.m"

clear, clf
format short
%
% Number of nodes in network
Nodes=8;
%
% Creating the adjacency matrix
Distance_Matrix=[ 0 20  0  0  0  0 15  0 ; ...
	    20  0  8  9  0  0  0  0 ; ...
             0  8  0  6 15  0  0 10 ; ... 
	     0  9  6  0  7  0  0  0 ; ...
             0  0 15  7  0 22 18  0 ; ...
             0  0  0  0 22  0  0  0 ; ...
            15  0  0  0 18  0  0  0 ; ...
             0  0 10  0  0  0  0  0 ];
%
% Create the state matrix
State_Matrix=[ Inf 0  0 ; ...
	       Inf 0  0 ; ...
	       Inf 0  0 ; ...
	       Inf 0  0 ; ...
	       Inf 0  0 ; ...
	       Inf 0  0 ; ...
	       Inf 0  0 ; ...
	       Inf 0  0 ]; 
%
% Define the coordinates of the nodes.
Coord_Matrix=[ 0 5 ; 1 0 ; 5 1 ; 2 4 ; 3 6 ; 0 7 ; 4 8 ; 6 2 ];
%
% Drawing the network.
text(Coord_Matrix(1,1), Coord_Matrix(1,2), '1'),hold on
text(Coord_Matrix(2,1), Coord_Matrix(2,2), '2'),hold on
text(Coord_Matrix(3,1), Coord_Matrix(3,2), '3'),hold on
text(Coord_Matrix(4,1), Coord_Matrix(4,2), '4'),hold on
text(Coord_Matrix(5,1), Coord_Matrix(5,2), '5'),hold on
text(Coord_Matrix(6,1), Coord_Matrix(6,2), '6'),hold on
text(Coord_Matrix(7,1), Coord_Matrix(7,2), '7'),hold on
text(Coord_Matrix(8,1), Coord_Matrix(8,2), '8'),hold on
%
for i=1:Nodes
	for j=1:(i-1)
		if Distance_Matrix(i,j)~=0
			Graphx=[Coord_Matrix(i,1) Coord_Matrix(j,1)];
			Graphy=[Coord_Matrix(i,2) Coord_Matrix(j,2)];
			plot(Graphx,Graphy,':'),hold on
		end;
	end;
end;
%
% Ask for start and finish point
startPoint=input('Start point is:  ');
endPoint=input('End point is:  ');
%
axis('equal');
axis('off');

Now I need to compute the shortest path, while also to give the the first point and the end point, so that the program would calculate the shortest route and should display the total distance.

I have seen the video on YouTube and tried to insert it in the same script, however, the program ignores it and only draws the network and initial code and asks for start and end points, but the loops below are ignored:

% Shortest path 
function distances = dijkstra1(Distance_Matrix, State_Matrix)

distances = dijkstra1(Distance_Matrix,1);

% Start with a map of distances among Nodes 
%Initialize the distance to all points as infinity 
distances(1:Nodes) = Inf;

% Mark nodes as UNvisited 
visited(1:Nodes) = 0;

% Initialize the distance to the first point as zero. 
distances(State_Matrix) = 0; 

while sum(visited) < Nodes
        % Find all unvisited nodes
        possibles(1:Nodes) = Inf; 
        for index = 1:Nodes
            if visited(index) == 0
                possibles(index) = distances(index);
            end 
        end 
        % Then find the one with the smallest distance value 
        % Make this the current point, and its distance the current
        % distance. 
        [currentDistance, currentPoint] = min(possibles);

        % Given the distance to the current point, update the distance to
        % all its neighbors if the new distance will be smaller 
        for index = 1:Nodes
            newDistance = currentDistance + Distance_Matrix(currentPoint, index)
            if newDistance < distances(index)
                distances(index) = newDistance;
            end 
        end 

        % Mark the current node as visited 
        visited(currentPoint) = 1;
end
%  End
end

Could someone please tell me where is the mistake? Any hints?

Many thanks in advance.

[Torsten]

You must call the function "dijkstra1" from your main program ( and maybe change the function so that it only computes the distance between start and end point and not between all points ).

% Number of nodes in network
Nodes=8;
%
% Creating the adjacency matrix
Distance_Matrix=[ 0 20  0  0  0  0 15  0 ; ...
    20  0  8  9  0  0  0  0 ; ...
    0  8  0  6 15  0  0 10 ; ... 
    0  9  6  0  7  0  0  0 ; ...
    0  0 15  7  0 22 18  0 ; ...
    0  0  0  0 22  0  0  0 ; ...
    15  0  0  0 18  0  0  0 ; ...
    0  0 10  0  0  0  0  0 ];
%
% Create the state matrix
State_Matrix=[ Inf 0  0 ; ...
    Inf 0  0 ; ...
    Inf 0  0 ; ...
    Inf 0  0 ; ...
    Inf 0  0 ; ...
    Inf 0  0 ; ...
    Inf 0  0 ; ...
    Inf 0  0 ]; 
%
% Define the coordinates of the nodes.
Coord_Matrix=[ 0 5 ; 1 0 ; 5 1 ; 2 4 ; 3 6 ; 0 7 ; 4 8 ; 6 2 ];
%
% Drawing the network.
text(Coord_Matrix(1,1), Coord_Matrix(1,2), '1'),hold on
text(Coord_Matrix(2,1), Coord_Matrix(2,2), '2'),hold on
text(Coord_Matrix(3,1), Coord_Matrix(3,2), '3'),hold on
text(Coord_Matrix(4,1), Coord_Matrix(4,2), '4'),hold on
text(Coord_Matrix(5,1), Coord_Matrix(5,2), '5'),hold on
text(Coord_Matrix(6,1), Coord_Matrix(6,2), '6'),hold on
text(Coord_Matrix(7,1), Coord_Matrix(7,2), '7'),hold on
text(Coord_Matrix(8,1), Coord_Matrix(8,2), '8'),hold on
%
for i=1:Nodes
    for j=1:(i-1)
        if Distance_Matrix(i,j)~=0
            Graphx=[Coord_Matrix(i,1) Coord_Matrix(j,1)];
            Graphy=[Coord_Matrix(i,2) Coord_Matrix(j,2)];
            plot(Graphx,Graphy,':'),hold on
        end;
    end;
end;
%
% Ask for start and finish point
startPoint=3;%input('Start point is:  ');
endPoint=7;%input('End point is:  ');
%
axis('equal');
axis('off');

distances = dijkstra1(Distance_Matrix, State_Matrix);

Array indices must be positive integers or logical values.

Error in solution>dijkstra1 (line 66)
distances(State_Matrix) = 0;

distances(startPoint,endPoint)

% Shortest path 
function distances = dijkstra1(Distance_Matrix, State_Matrix)
Nodes = size(Distance_Matrix,1);
% Start with a map of distances among Vertices 
%Initialize the distance to all points as infinity 
distances(1:Nodes) = Inf;
% Mark nodes as UNvisited 
visited(1:Nodes) = 0;
% Initialize the distance to the first point as zero. 
distances(State_Matrix) = 0; 
while sum(visited) < Vertices
    % Find all unvisited nodes
    possibles(1:Nodes) = Inf; 
    for index = 1:Nodes
        if visited(index) == 0
            possibles(index) = distances(index);
        end 
    end 
    % Then find the one with the smallest distance value 
    % Make this the current point, and its distance the current
    % distance. 
    [currentDistance, currentPoint] = min(possibles);
    % Given the distance to the current point, update the distance to
    % all its neighbors if the new distance will be smaller 
    for index = 1:Nodes
        newDistance = currentDistance + Distance_Matrix(currentPoint, index)
        if newDistance < distances(index)
            distances(index) = newDistance;
        end 
    end 
    % Mark the current node as visited 
    visited(currentPoint) = 1;
end
%  End
end

[Dan Kristen]

Ah thanks for response..

What is with this expression? Seems it's not finished - got error

distances(startPoint,endPoint)

ThemeCopy

Check for incorrect argument data type or missing argument in call to function 'distances'.
Error in dijkstra (line 80)
distances(startPoint,endPoint)

[Torsten]

I don't know the format in which the distances are saved. Usually it's a matrix where distances(i,j) is the distance between node i and node j. You will have to check your function for that.

[Dan Kristen]

oh what do you mean by format? Where to check that?

I wrote everything in the same script, then I ran the function in the separate file.

But I am wondering how should be written this expression:

distances(startPoint,endPoint)

[Dan Kristen]

https://uk.mathworks.com/help/matlab/ref/graph.distances.html

you mean here, though it shows error again:

distances(Adj_Matrix) = distances(startPoint,endPoint)

Check for incorrect argument data type or missing argument in call to function 'distances'.

[Walter Roberson]

that function needs to be passed a graph() or digraph()

[Dan Kristen]

Like this?

total_distance = graph(distances(startPoint,endPoint));

[Torsten]

Walter Roberson means that a graph of your network has to be passed to "dijkstra1".

Where did you get this code from ? Because the variables "Nodes" and "Vertices" are undefined when you call it like you do.

[Dan Kristen]

I watched this video:

https://www.youtube.com/watch?v=otJAMukL3Iw

But yeah I forgot to change vertices to nodes. Corrected it.

[Torsten]

Nodes was also not defined - just added the definition in the first line.

Take one of the code

https://de.mathworks.com/matlabcentral/fileexchange/36140-dijkstra-algorithm

https://de.mathworks.com/matlabcentral/fileexchange/46883-dijkstra-algorithm

from the file exchange. They seem to be exactly what you are looking for. And they are documented ...

[Dan Kristen]

I just tried, it does not work. many errors, as well as outdated functions.

[Torsten]

Here is the documentation of inputs and outputs:

https://gjacopo.github.io/imtools/graph/dijk.html

A=[ 0 20  0  0  0  0 15  0 ; ...
    20  0  8  9  0  0  0  0 ; ...
    0  8  0  6 15  0  0 10 ; ...
    0  9  6  0  7  0  0  0 ; ...
    0  0 15  7  0 22 18  0 ; ...
    0  0  0  0 22  0  0  0 ; ...
    15  0  0  0 18  0  0  0 ; ...
    0  0 10  0  0  0  0  0 ];
start_verts = 1;
end_verts = 8;
[dist, paths] = dijk(A, start_verts, end_verts)

dist = 38

paths = 4×1

1 2 3 8

function [dist, varargout] = dijk(W, start_verts, end_verts)

[dist, P] = dijk_(W, start_verts, end_verts);
if nargout==2
    varargout{1} = pred2path_(P, start_verts, end_verts);
end

end % end of dijk

%--------------------------------------------------------------------------
function [D, P] = dijk_(A,s,t)
% Copyright (c) 1994-2003 by Michael G. Kay
% Matlog Version 7 02-Sep-2003

[n,cA] = size(A);

if nargin < 2 || isempty(s), s = (1:n)'; else s = s(:); end
if nargin < 3 || isempty(t), t = (1:n)'; else t = t(:); end

if ~any(any(tril(A) ~= 0))			% A is upper triangular
    isAcyclic = 1;
elseif ~any(any(triu(A) ~= 0))	% A is lower triangular
    isAcyclic = 2;
else										% Graph may not be acyclic
    isAcyclic = 0;
end

if n ~= cA
    error('A must be a square matrix');
elseif ~isAcyclic && any(any(A < 0))
    error('A must be non-negative');
elseif any(s < 1 | s > n)
    error(['''s'' must be an integer between 1 and ',num2str(n)]);
elseif any(t < 1 | t > n)
    error(['''t'' must be an integer between 1 and ',num2str(n)]);
end

A = A';		% Use transpose to speed-up FIND for sparse A

D = zeros(length(s),length(t));
P = zeros(length(s),n);

for i = 1:length(s)
    j = s(i);
    
    Di = Inf*ones(n,1); Di(j) = 0;
    
    isLab = false(length(t),1);
    if isAcyclic ==  1
        nLab = j - 1;
    elseif isAcyclic == 2
        nLab = n - j;
    else
        nLab = 0;
        UnLab = 1:n;
        isUnLab = true(n,1);
    end
    
    while nLab < n && ~all(isLab)
        if isAcyclic
            Dj = Di(j);
        else	% Node selection
            [Dj,jj] = min(Di(isUnLab));
            j = UnLab(jj);
            UnLab(jj) = [];
            isUnLab(j) = 0;
        end
        
        nLab = nLab + 1;
        if length(t) < n, isLab = isLab | (j == t); end
        
        [jA,~,Aj] = find(A(:,j));
        Aj(isnan(Aj)) = 0;
        
        if isempty(Aj), Dk = Inf; else Dk = Dj + Aj; end
        
        P(i,jA(Dk < Di(jA))) = j;
        Di(jA) = min(Di(jA),Dk);
        
        if isAcyclic == 1			% Increment node index for upper triangular A
            j = j + 1;
        elseif isAcyclic == 2	% Decrement node index for lower triangular A
            j = j - 1;
        end
        
        %disp( num2str( nLab ));
    end
    D(i,:) = Di(t)';
end

end

%--------------------------------------------------------------------------
function rte = pred2path_(P,s,t)
% Copyright (c) 1994-2006 by Michael G. Kay
% Matlog Version 9 13-Jan-2006 (http://www.ie.ncsu.edu/kay/matlog)

[~,n] = size(P);

if nargin < 2 || isempty(s), s = (1:n)'; else s = s(:); end
if nargin < 3 || isempty(t), t = (1:n)'; else t = t(:); end

if any(P < 0 | P > n)
    error(['Elements of P must be integers between 1 and ',num2str(n)]);
elseif any(s < 1 | s > n)
    error(['"s" must be an integer between 1 and ',num2str(n)]);
elseif any(t < 1 | t > n)
    error(['"t" must be an integer between 1 and ',num2str(n)]);
end

rte = cell(length(s),length(t));

[~,idxs] = find(P==0);

for i = 1:length(s)
    %    if rP == 1
    %       si = 1;
    %    else
    %       si = s(i);
    %       if si < 1 | si > rP
    %          error('Invalid P matrix.')
    %       end
    %    end
    si = find(idxs == s(i));
    for j = 1:length(t)
        tj = t(j);
        if tj == s(i)
            r = tj;
        elseif P(si,tj) == 0
            r = [];
        else
            r = tj;
            while tj ~= 0
                if tj < 1 || tj > n
                    error('Invalid element of P matrix found.')
                end
                r = [P(si,tj); r];                                          %#ok
                tj = P(si,tj);
            end
            r(1) = [];
        end
        rte{i,j} = r;
    end
end

if length(s) == 1 && length(t) == 1
    rte = rte{:};
end

%rte = t;
while 0%t ~= s
    if t < 1 || t > n || round(t) ~= t
        error('Invalid "pred" element found prior to reaching "s"');
    end
    rte = [P(t) rte];                                                   %#ok
    t = P(t);
end
end

[Dan Kristen]

Thanks Torsten!

seems challenging even just to understand this code.

I got the errors though: W - is adjacency matrix?

how it is then expressed in the rest of the code? because later in the code it is expressed as A?

 Not enough input arguments.
Error in dijk1 (line 2)
[dist, P] = dijk_(W, start_verts, end_verts);
 
>> dijk1
Not enough input arguments.

Error in dijk1 (line 5)
[n,cA] = size(A);
 

[Torsten]

how it is then expressed in the rest of the code? because later in the code it is expressed as A?

In each function, the variables names can change for the input data from those of the calling function.

I don't know which changes you made to the input matrix A. As you can see above, the code works with the cost matrix you defined.

All information about the code is on the page I gave you the link to:

Syntax

[dist, paths] = DIJK(A, start_verts, end_verts);

Inputs

A : (n,n) node-node weighted adjacency matrix of arc lengths.

A(i,j) = 0 => Arc (i,j) does not exist;

A(i,j) = NaN => Arc (i,j) exists with null weight.

start_verts : FROM node indices; default: s=[], ie. paths from all nodes.

end_verts : TO node indices; default: s=[], ie. paths to all nodes.

Outputs

D : matrix of size (length(s),length(t)) storing the shortest path distances from s to t: D(i,j) is the distance from node i to node j.

P : shortest paths.

Remarks

• If A is a triangular matrix, then computationally intensive node selection step not needed since graph is acyclic (triangularity is a sufficient, but not a necessary, condition for a graph to be acyclic) and A can have non-negative elements.
• If length(s)>>length(t), then DIJK is faster if DIJK(A',t,s) used, where D is now transposed and P represents successor indices.

References

[AMO93] R. Ahuja, T. Magnanti and J. Orlin: "Network Flows: Theory, Algorithms, and Applications", Prentice-Hall, 1993.

[DIJK] Source code made available by Michael G. Kay (see copyright) at http://web.mit.edu/cocosci/isomap/code/dijk.m http://cbi.nyu.edu/svn/mrTools/trunk/mrLoadRet/ROI/pred2path.m

Acknowledgment

Copyright (c) 1994-2006 by Michael G. Kay Matlog Version 9 13-Jan-2006 (http://www.ie.ncsu.edu/kay/matlog)