Find cycles in an undirected graph

Asked by Sim on 7 Oct 2019

Latest activity Commented on by Matt J on 26 Dec 2019 at 17:25

104 views (last 30 days)

Hi, I need to find cycles in a graph, exactly as it was asked here (and apparently without fully clear/working solutions!):

find cycle in array https://ch.mathworks.com/matlabcentral/answers/425321-find-cycle-in-array
find a cycles in undirected graph https://ch.mathworks.com/matlabcentral/answers/421435-find-a-cycles-in-undirected-graph

Here my example/code:

clear all; clc; close all;
figure('Color',[1 1 1]);
s = [1 1 1 2 3 3 4 4 5 6 7 6 8 9 10 10 12 12 13 14 15 16 17 17 18 19 20 21 20 25];
t = [2 8 18 3 4 23 5 21 6 7 8 11 9 10 11 12 14 13 15 18 16 17 18 25 19 20 1 22 24 26];
G = graph(s,t);
x = [0.5 0 0 0 0.5 1 1.5 2 3 3 3 5.5 6 4 6 6 4 3 2 0.5 -1 -2 -1 1.5 4.5 4.5];
y = [0 0.5 1 1.5 2 2 1.5 1 1 1.5 2 1 0.5 0.5 0 -1 -1 -0.5 -1 -1 1 0.5 0.5 -0.5 -0.5 0];
h = plot(G,'XData',x,'YData',y,'linewidth',2,'MarkerSize',7);
nl = h.NodeLabel;
h.NodeLabel = '';
xd = get(h, 'XData');
yd = get(h, 'YData');
text(xd, yd, nl, 'FontSize',17, 'FontWeight','bold', ...
    'HorizontalAlignment','left', 'VerticalAlignment','top')
set(gca,'Fontsize',15,'FontWeight','Bold','LineWidth',2, 'box','on')
% Remove "branches"
xy = [x' y'];
while ~isempty(find(degree(G)==1))
degreeone = find(degree(G)==1);
G = rmnode(G,degreeone);
xy(degreeone,:) = [];
end

Here the corresponding Figure (after removal of "branches"):

My goal would be to find the following 5 cycles as output (i.e. lists of nodes composing each cycle):

1-2-3-4-5-6-7-8-1
6-7-8-9-10-11-6
1-8-9-10-12-14-18-1
1-18-19-20-1
12-13-15-16-17-18-14-12

Note 1:

This is the Sergii Iglin 's idea, that I found in https://ch.mathworks.com/matlabcentral/fileexchange/4266-grtheory-graph-theory-toolbox

This method is partially working for my purposes.

Unfortunately, the 2nd and 4th cycles are not what I needed/expected.

% Sergii Iglin
% https://iglin.org/All/GrMatlab/grCycleBasis.html
E = table2array(G.Edges);
Output_SI = grCycleBasis(E);
% [my part] From the Sergii Iglin's output to cycles nodes
for i = 1 : size(Output_SI,2)
    w = [];
    u = E(find(Output_SI(:,i)),:); % edges list
    w(1) = u(1,1);
    w(2) = u(1,2);
    u(1,:) = [];
    j = 2;
    while ~isempty(u)
        [ind,~] = find(u==w(j));
        [~,ind2] = ismember(u, u(ind,:), 'rows');
        g = u( ind2==1 ,:) ~= w(j);
        w(j+1) = u( ind2==1 , g);
        u( ind2==1 ,:) = [];
        j = j + 1;
    end
    cycles_SI{i} = w;
end
% Sergii Iglin's results
>> cycles_SI{:}
     1     2     3     4     5     6     7     8     1
     1     2     3     4     5     6    11    10     9     8     1
     1     8     9    10    12    14    18     1
     1     8     9    10    12    13    15    16    17    18     1
     1    18    19    20     1

Note 2:

This is the Christine Tobler 's idea that I found in https://ch.mathworks.com/matlabcentral/answers/353565-are-there-matlab-codes-to-compute-cycle-spaces-of-graphs

This method is partially working for my purposes.

Unfortunately, the 2nd and 4th cycles are not what I needed/expected.

% Christine Tobler
% https://ch.mathworks.com/matlabcentral/answers/353565-are-there-matlab-codes-to-compute-cycle-spaces-of-graphs
t = minspantree(G, 'Type', 'forest');
%  highlight(h,t)
nonTreeEdges = setdiff(G.Edges.EndNodes, t.Edges.EndNodes, 'rows');
cycles_CT = cell(size(nonTreeEdges, 1), 1);
for i = 1 : length(cycles_CT)
    src = nonTreeEdges(i, 1);
    tgt = nonTreeEdges(i, 2);
    cycles_CT{i} = [tgt shortestpath(t, src, tgt)];
end
% Christine Tobler's results
>> cycles_CT{:}
     8     7     6     5     4     3     2     1     8
    11    10     9     8     1     2     3     4     5     6    11
    18    14    12    10     9     8     1    18
    18    17    16    15    13    12    10     9     8     1    18
    20    19    18     1    20

Note 3:

Methods from Sergii Iglin and Christine Tobler give the same result!

Note 4:

The ideas / FileExchange submissions

Count all cycles in simple undirected graph version 1.2.0.0 (5.43 KB) by Jeff Howbert
Count Loops in a Graph version 1.1.0.0 (167 KB) by Joseph Kirk

kindly suggested here

https://ch.mathworks.com/matlabcentral/answers/425321-find-cycle-in-array

are not working for my case...

Any further idea / suggestion ?

Thanks a lot!

4 Comments

Show 1 older comment

Matt J on 8 Oct 2019

Because the ulterior motive is to partition the graph into non-overlapping polygons.

Sim on 9 Oct 2019

I found also this idea/method "Polygons From Set of Line segments", explained in

http://ac.aua.am/sergey_yeranosyan/Web/algorithm%20project/

and based on Ferreira A., Fonseca M.J., Jorge J.A. (2003) Polygon detection from a set of lines

http://web.ist.utl.pt/alfredo.ferreira/publications/12EPCG-PolygonDetection.pdf

The basic idea can be summarised as: "Finding small polygons is the same as searching for a Minimum Cycle Basis"

Screen Shot 2019-10-09 at 2.40.33 PM.png

% Pseudocode of the algorithm
% Minimum basis cycle [3]
% Initialize empty sets F,P
% P=All-Pairs-Shortest-Paths(G)
% for each v ∈ V
% for each (x,y) ∈ E
% if Px,v ∩ Pv,y = {v}
% C=Px,v ∪ Pv,y ∪ (x,y)
% add C to F
% Order-By-Length
% return Select-Cycles(F)
% Polygons from cycles [4]
% Initialize empty set P
% for each cycle C in F
% p=new polygon
% for each vertex v ∈ V
% add vertex v to p
% add polygon p to set P
% return P

Sim on 9 Oct 2019

@Steven Lord : Kindly Matt J gave the answer... In this moment, I would just need to know the "cycles" corresponding to the polygons which compose my network...

2 Answers

Link

Answer by Matt J on 8 Oct 2019

Edited by Matt J on 8 Oct 2019

Accepted Answer

polyregions.m

Because this sounds like a generally useful thing, I cooked up the attached polyregions class to do the partitioning that you described. It uses graph theoretic functions only.

Here is its application to the data example that you provided. Each partitioned polygon is contained in the polyshape array, pgon.

s = [1 1 1 2 3 3 4 4 5 6 7 6 8 9 10 10 12 12 13 14 15 16 17 17 18 19 20 21 20 25];
t = [2 8 18 3 4 23 5 21 6 7 8 11 9 10 11 12 14 13 15 18 16 17 18 25 19 20 1 22 24 26];
G = graph(s,t);
x = [0.5 0 0 0 0.5 1 1.5 2 3 3 3 5.5 6 4 6 6 4 3 2 0.5 -1 -2 -1 1.5 4.5 4.5];
y = [0 0.5 1 1.5 2 2 1.5 1 1 1.5 2 1 0.5 0.5 0 -1 -1 -0.5 -1 -1 1 0.5 0.5 -0.5 -0.5 0];
obj=polyregions(G,x,y);
pgon=polyshape(obj);
plot(obj);
hold on
plot(pgon);
hold off

33 Comments

Show 30 older comments

Sim on 10 Oct 2019

Hi Matt, a question about the removal of "branches" in your code:

function obj = prune(obj)
%Remove dangling branches

....In case I have a small loop inside one of the dangling "branches" that we need to remove, like the loop 25-26-27, would you have some suggestion to workaround it ?

The code corresponding to that Figure, with the looped-branch is the following:

s = [1 1 1 2 3 3 4 4 5 6 7 6 8 9 10 10 12 12 13 14 15 16 17 17 18 19 20 21 20 25 26 27 27];
t = [2 8 18 3 4 23 5 21 6 7 8 11 9 10 11 12 14 13 15 18 16 17 18 25 19 20 1 22 24 26 27 25 28];
G = graph(s,t);
x = [0.5 0 0 0 0.5 1 1.5 2 3 3 3 5.5 6 4 6 6 4 3 2 0.5 -1 -2 -1 1.5 4.5 4.5 5 5];
y = [0 0.5 1 1.5 2 2 1.5 1 1 1.5 2 1 0.5 0.5 0 -1 -1 -0.5 -1 -1 1 0.5 0.5 -0.5 -0.5 0 0 0.5];
xy = [x' y']
h = plot(G,'XData',xy(:,1),'YData',xy(:,2),'linewidth',2,'MarkerSize',7);
axis([-3 7 -1.5 2.5]);
nl = h.NodeLabel;
h.NodeLabel = '';
xd = get(h, 'XData');
yd = get(h, 'YData');
text(xd, yd, nl, 'FontSize',17, 'FontWeight','bold', ...
    'HorizontalAlignment','left', 'VerticalAlignment','top')
set(gca,'Fontsize',15,'FontWeight','Bold','LineWidth',2, 'box','on')

If I apply your code...

% Matt J Polygons/Cycles
obj=polyregions(G,x,y);
pgon=polyshape(obj);
plot(obj);
hold on
plot(pgon);
hold off
[pgon,nodeIndices]=polyshape(obj);
nodeIndices{:}

Inside the function

function [P,E]=polyshortest(obj,arg1,arg2)

I get this error:

nodeA =
    17
nodeB =
    21
P =
     []
Index exceeds array bounds.
Error in polyregions/polyshortest (line 90)
            E=G.findedge(P,[P(2:end),P(1)])
Error in polyregions/polyshape (line 118)
                [P,E]=polyshortest(obj,E0); % F4              
Error in GraphCycles (line 22)
pgon=polyshape(obj);

Thanks a lot!

Sim on 11 Oct 2019

Hi Matt, your solution perfectly works on this network..! Thanks a lot! :)

...However, when I try to apply this revised prune() method to my real case (a larger network), I get this error:

Index exceeds array bounds.
Error in polyregions (line 15)
              w=sqrt((x(s)-x(t)).^2 + (y(s)-y(t)).^2);
Error in polyregions/prune (line 44)
            obj=polyregions(graph(sk,tk),xp,yp);
Error in polyregions/polyshape (line 85)
            obj=prune(obj);
Error in floods (line 98)
pgon=polyshape(obj);

At a first glance, it looks like that the components of the "w" array (i.e. x(s), x(t), y(s), y(t)) cannot be computed due to

Index exceeds array bounds.

Trying to understand the reasons of this error, I have noticed that the x and y arrays, which contain the nodes coordinates, reduce their sizes when the function prune() is called. This occurs via the xp and yp arrays, that are initialised as x and y arrays, at the beginnning of the prune() function. Indeed, when we remove the "dangling branches", the xp and yp arrays are reduced:

xp(tips)=[];
yp(tips)=[];

The xp and yp arrays are then "exported" to the new "obj", at the end of the prune() function

obj=polyregions(graph(sk,tk),xp,yp);

Therefore, in the following/remaining part of the code, the updated graph (without dangling branches) contains the reduced x and y arrays.

Now, if we change the last line of the prune() function as

obj=polyregions(graph(sk,tk),obj.x,obj.y);

where obj.x and obj.y have the original x and y arrays, the code runs further, and almost at the very end, where pgon is calculated

U=union(pgon) %#ok<*LTARG>
if ~U.NumHoles, return; end
hgon=holes(U)
pgon=[pgon,hgon]
for i=1:numel(hgon)
    [~,P]=ismember(hgon(i).Vertices,[x(:),y(:)],'rows')
    nodeIndices{end+1}=P.' %#ok<*AGROW>
end

....and exactly there I get this error:

U = 
  polyshape with properties:
      Vertices: [349×2 double]
    NumRegions: 1
      NumHoles: 15
hgon = 
  15×1 polyshape array with properties:
    Vertices
    NumRegions
    NumHoles
Error using horzcat
Dimensions of arrays being concatenated are not consistent.
Error in polyregions/polyshape (line 117)
            pgon=[pgon,hgon]
Error in floods (line 98)
pgon=polyshape(obj);   

Am I doing something wrong?

Thanks a lot!

Sim on 12 Oct 2019

Thanks... By using this line of code

nodeIndices{end+1}=lp(P.');

I was getting the following error

Array indices must be positive integers or logical values.
Error in jigsaw/polyshape (line 137)
             nodeIndices{end+1}=lp(P.');
Error in program (line 94)
[pgon,cycles]=polyshape(obj)

....Then, I checked the previous line

[~,P]=ismember(hgon(i).Vertices,[x(:),y(:)],'rows');

and the three quantites inside it, i.e.

P
hgon(i).Vertices
[x(:),y(:)]

As example, I got something like this:

P =
         188
        1873
         191
           0
           0
        1902
         245
         281
           0
           0
           0
         177
         178
hgon(1).Vertices = 
 
                719206.761                262076.829
                720011.091                260159.216
                720306.696                259519.618
           720812.14899999           258050.23099999
           722903.21299998                256225.431
                726046.209                256652.186
                727673.923                255337.524
                732122.269                253113.636
           736764.74299995           256159.99799996
           728781.39999995           267902.73300001
                725492.548           269719.18599999
                722128.621                271490.111
                723849.229                268561.544
[x(:),y(:)] = 
                704726.255                    281756
                711700.953                265115.835
                712542.533                264359.824
                713410.689                263107.846
                713769.936                 262368.76
                714125.573                261718.406
                715803.539                261306.605
                717089.203                260113.503
                717635.691                259625.701
                694183.174                 253518.56
                etc...                    etc...

Then, I thought to replace this line

[~,P]=ismember(hgon(i).Vertices,[x(:),y(:)],'rows');

with (ismember with a tolerance)

[~,P]=ismembertol(hgon(i).Vertices,[x(:),y(:)],0.02,'ByRows',true)

and I did not get any error anymore... But I am just wondering why the coordinates of some "vertices" are not perfectly overlapping during the operation of "ismember", since all of those coordinates come from the same initial dataset...

In summary, I modified the very last loop of the code from this

for i=1:numel(hgon)
    [~,P]=ismember(hgon(i).Vertices,[x(:),y(:)],'rows');
    nodeIndices{end+1}=lp(P.'); %#ok<*AGROW>
end

to this:

for i=1:numel(hgon)
    [~,P]=ismembertol(hgon(i).Vertices,[x(:),y(:)],0.02,'ByRows',true)
    nodeIndices{end+1}=lp(P.'); %#ok<*AGROW>
end

and the last error disappeared....However, I am not sure if the usage of "ismembertol" is fully correct, instead of the previous "ismember"

Sim on 14 Oct 2019

jigsaw.m

Many thanks Matt :) I used what you kindly suggested (something similar)

dbstop in jigsaw if warning

in my list of commands... i.e. in this way...

dbstop in jigsaw if warning
obj=jigsaw(G,x,y);
obj.tol
[pgon,cycles]=polyshape(obj);
plot(obj)
hold on
plot(pgon);
hold off

....and I found out something, that I explaine here below...

Let's start from this warning

Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected
results. Input data has been modified to create a well-defined polyshape. 

which disappeares if, inside the

function [pgon,labelIndices]=polyshape(obj)

I replace this line

pgon(end+1)=polyshape(Vall(P,:),'Simplify',true);

by this one

pgon(end+1)=polyshape(Vall(P,:),'Simplify',false);

In your opinion, is this replacement a problem/issue for the calculation of polygons ?

Note: in general, I can see that the warnings disappear if I use "Simplify, false"

polyshape(..,'Simplify',false);

Let's go now to the "network squeezed towards the right side of the Figure window" and to the other warnings, which occur all together at once, as follows

Warning: Output contains an empty polyshape due to duplicate vertices, intersections, or other inconsistencies. 
Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected
results. Input data has been modified to create a well-defined polyshape. 
Warning: Boundaries with less than 3 points were removed. 

....All these warnings disappear and the picture is restored to the center as previously, if I replace the last loop of the code

for i=1:numel(hgon)
    P= knnsearch(Vall,hgon(i).Vertices); 
    pgon(end+1)=polyshape(Vall(P,:),'Simplify',true);
    labelIndices{end+1}=lp(P.'); %#ok<*AGROW>
end

with that one we used in a previous version of jigsaw.m

pgon=[pgon(:).',hgon(:).'];
for i=1:numel(hgon)
    [~,P]=ismembertol(hgon(i).Vertices,[x(:),y(:)],0.02,'ByRows',true);
    labelIndices{end+1}=lp(P.');
end  

By using the jigsaw.m code with the above mentioned modifications (see attached file too) I checked the "visual results" for both your network and mine.... At a first glance, it looks like that both networks give correct visual results, i.e. all polygons/holes are detected... However, I am aware that "ismembertol" is not the best, while the usage of "knnsearch" is a more robust solution.. but employing "knnsearch", I got many warnings and a squeezed picture which I would not know how to fix for my network/case...

Just for information, please see below the code I used - that you kindly provided - and the corresponding Figure I got with the attached (slightly) modified jigsaw.m

N=1000;
x = randn(N,1)*100;
y = randn(N,1)*100;
DT=delaunayTriangulation(x,y);
EdgeTable=table(DT.edges,'VariableNames',{'EndNodes'});
G=graph(EdgeTable);
% dbstop in jigsaw if warning
obj=jigsaw(G,x,y);
obj.tol
tic;
[pgon,cycles]=polyshape(obj);
plot(obj)
toc %Elapsed time is 16.075788 seconds.
hold on
plot(pgon);
hold off

Screen Shot 2019-10-15 at 12.35.34 AM.png

Sim on 9 Dec 2019

Cool!! I am sure it will be very useful to many people :) ..Hopefully I will cite you soon :)

Tarik Chowdhury on 24 Dec 2019 at 17:56

Hi, I have also used the same code for my data.

s = [1 1 2 2 2 3 4 4 4 5 6 6 6 7 7 9 9 10 12 13];
t = [2 5 3 4 5 4 5 7 9 6 11 12 13 8 9 10 14 11 13 14];
G = graph(s,t);
x = [15 13 14.5 11.5  11.5 8 9.5 20  7.0 5.5 5.9 5.5 4.5 2.7];
y = [03 04 5.50 06.5  02.5 2 9.5 20 6.5 5.0 3.5 1.0 3.0 6.0];
obj=polyregions(G,x,y);
pgon=polyshape(obj);
plot(obj);
hold on
plot(pgon);
hold off

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

I am getting the loops as following:

1 5 2

2 4 3

4 2 5

4 8 7

5 4 8 9 10 6

6 12 11

8 9 10 6 12 13

at the 4th loop why node * is coming also the node 8 is creating problem with some other loops too.

Matt J on 26 Dec 2019 at 17:25

@Tarik,

It looks like you are using an early, non-perfected version of the code. Here is what I get using spatialgraph2D from the link mentioned above,

s = [1 1 2 2 2 3 4 4 4 5 6 6 6 7 7 9 9 10 12 13];
t = [2 5 3 4 5 4 5 7 9 6 11 12 13 8 9 10 14 11 13 14];
G = graph(s,t);
x = [15 13 14.5 11.5  11.5 8 9.5 20  7.0 5.5 5.9 5.5 4.5 2.7];
y = [03 04 5.50 06.5  02.5 2 9.5 20 6.5 5.0 3.5 1.0 3.0 6.0];
obj=spatialgraph2D(G,x,y);
[pgon,loops]=polyshape(obj);
plot(obj);
hold on
plot(pgon);
hold off

>> loops{:}
ans =
     1     5     2
ans =
     2     4     3
ans =
     4     2     5
ans =
     4     9     7
ans =
     5     4     9    10    11     6
ans =
     6    13    12
ans =
     9    10    11     6    13    14

Link

Answer by darova on 7 Oct 2019

Just use for loop and cells since you already know indices of each polygon

s = [1 1 1 2 3 3 4 4 5 6 6 6 8 9 10 10 12 12 13 14 15 16 17 17 18 19 20 21 20 25];
t = [2 8 18 3 4 23 5 21 6 7 8 11 9 10 11 12 14 13 15 18 16 17 18 25 19 20 1 22 24 26];
x = [0.5 0 0 0 0.5 1 1.5 2 3 3 3 5.5 6 4 6 6 4 3 2 0.5 -1 -2 -1 1.5 4.5 4.5];
y = [0 0.5 1 1.5 2 2 1.5 1 1 1.5 2 1 0.5 0.5 0 -1 -1 -0.5 -1 -1 1 0.5 0.5 -0.5 -0.5 0];
ind = {{1 2 3 4 5 6 7 8 1}
    {6 7 8 9 10 11 6}
    {1 8 9 10 12 14 18 1}
    {1 18 19 20 1}
    {12 13 15 16 17 18 14 12}};
cla
% plot([x(s);x(t)],[y(s);y(t)],'.b')
hold on
for i = 1:length(ind)
    ix = cell2mat(ind{i});
    plot(x(ix),y(ix),'color',rand(1,3))
end
hold off
axis equal

5 Comments

Show 2 older comments

Sim on 7 Oct 2019

Thanks Darova!

darova on 8 Oct 2019

Untitled2.m

Here is an example. Any ideas how to remove branches?

See the attached script

Sim on 9 Oct 2019

Hi Darova, to remove branches I used this:

% Remove branches
xy = [x' y'];
while ~isempty(find(degree(G)==1))
degreeone = find(degree(G)==1);
G = rmnode(G,degreeone);
xy(degreeone,:) = [];
end

About your idea/proposal, it looks cool and promising! Also, a nice animation!

However, is there any way to extrapolate the nodes composing each "cycle"/"polygon" that you were able to isolate ?

Thanks a lot for your efforts!

You are now following this question