function cycles = find_cycles(conn,node,cycles_through_node)
% give histogram of cycle length
% cycles is a list where each row gives
% [cycle_length #cycles_with_that_length]
% if cycles_through_node==1, only list cycles through node
% Note: lists cycles twice. If cycles go through node, the cycle is
% counted for each outgoing neighbor. If the cycle does not go through
% node, it is counted on branch where cycle is found twice.
% Note: this assumes that node has more than one connection!!
lc = length(conn);
% Make a spanning tree centered around node
% where each vertex in the tree is placed so that
% its shortest distance in graph is its distance to node.
% Keep previous vertex
% Keep address for each vertex based on branch points.
% Keep distance from each vertex to closest branch point.
% Flag for branch when just AFTER branch node to tell
% which direction to go from tree
E_tree = []; % edges in tree
E_cycle = []; % edges that make a cycle
% 1st_address1 1st_address2 cell1 cell2
rank_tree = struct('prev',cell(lc,1),'address',cell(lc,1),...
'dist',cell(lc,1),'branch',cell(lc,1));
S = []; % searched
R= node; %reached
rank_tree(node) = struct('prev',0,'address',[],....
'dist',0,'branch',0);
while ~isempty(R)
%fprintf(1,'current cell: %d\n',R(1));
nc = length(conn{R(1)});
% decide if we are at branching point,
% make note to add address if we are
if nc-length(intersect([R S],conn{R(1)}))>=2
branch_node = 1;
else
branch_node = 0;
end
for k=1:length(conn{R(1)})
if isempty(find([R S]==conn{R(1)}(k),1))
%fprintf(1,'not reached yet: %d\n',conn{R(1)}(k));
R = [R conn{R(1)}(k)];
rank_tree(conn{R(1)}(k)).prev = R(1);
% is R(1) part of the address,
% i.e. branch of previous branch point?
if rank_tree(R(1)).branch
rank_tree(conn{R(1)}(k)).dist = 1;
rank_tree(conn{R(1)}(k)).address = ...
[rank_tree(R(1)).address R(1)];
else
rank_tree(conn{R(1)}(k)).dist = ...
rank_tree(R(1)).dist + 1;
rank_tree(conn{R(1)}(k)).address = ...
rank_tree(R(1)).address;
end
if branch_node
rank_tree(conn{R(1)}(k)).branch = 1;
else
rank_tree(conn{R(1)}(k)).branch = 0;
end
elseif conn{R(1)}(k) ~= rank_tree(R(1)).prev
%fprintf(1,'already in tree: %d\n',conn{R(1)}(k));
% Neighbor is already in tree, so add edge to E_cycle if it's not
% already there. Add edges listing lower node first so that
% they're easier to find.
if R(1)<=conn{R(1)}(k)
%fprintf(1,'%d <= %d\n',R(1),conn{R(1)}(k));
if isempty(E_cycle) || ...
~any(E_cycle(:,3)==R(1) & E_cycle(:,4)==conn{R(1)}(k))
% edge not in E_cycle yet
E_cycle = [E_cycle; branch(R(1),rank_tree)...
branch(conn{R(1)}(k),rank_tree)...
R(1) conn{R(1)}(k)];
% fprintf(1,'1 - added cycle edge: %d %d %d %d\n',...
% branch(R(1),rank_tree),...
% branch(conn{R(1)}(k),rank_tree),...
% R(1), conn{R(1)}(k));
end
else % R(1) > conn{R(1)}(k)
%fprintf(1,'%d > %d\n',R(1),conn{R(1)}(k));
if isempty(E_cycle) || ...
~any(E_cycle(:,4)==R(1) & E_cycle(:,3)==conn{R(1)}(k))
% edge not in E_cycle yet
E_cycle = [E_cycle; ...
branch(conn{R(1)}(k),rank_tree)...
branch(R(1),rank_tree)...
conn{R(1)}(k) R(1)];
% fprintf(1,'2 - added cycle edge: %d %d %d %d\n',...
% branch(conn{R(1)}(k),rank_tree),...
% branch(R(1),rank_tree),...
% conn{R(1)}(k),R(1));
end
end
end
end
S = [S R(1)];
R = R(2:length(R));
end
% for i=1:length(conn)
% fprintf(1,'%d: address=[',i)
% for j=1:length(rank_tree(i).address)
% fprintf(1,'%d ',rank_tree(i).address(j))
% end
% fprintf(1,'], dist=%d\n',rank_tree(i).dist);
% end
% Cycles going through node are made from combinations
% of edges in E_cycle. First, compute cycles with
% one extra edge
%fprintf(1,'length E_cycle = %d\n',length(E_cycle(:,1)));
n_branches = length(conn{node});
sE = size(E_cycle);
cycles = cell(n_branches,1);
%E = E_cycle;
for i=1:sE(1)
dist = zeros(4,1);
if ~cycles_through_node && E_cycle(i,1)==E_cycle(i,2)
% cycle is on same branch, doesn't go through node
%fprintf(1,'edge is on single branch: %d %d %d %d\n',E_cycle(i,1),...
% E_cycle(i,2),E_cycle(i,3),E_cycle(i,4));
merge_ind = length(intersect(rank_tree(E_cycle(i,3)).address,...
rank_tree(E_cycle(i,4)).address)) + 1;
for j=3:4
if merge_ind<=length(rank_tree(E_cycle(i,j)).address)
dist(j) = rank_tree(E_cycle(i,j)).dist;
else
dist(j) = 0;
end
%fprintf(1,'node dist = %d\n',dist(j));
for k=(merge_ind+1):length(rank_tree(E_cycle(i,j)).address)
stop = rank_tree(E_cycle(i,j)).address(k);
dist(j) = dist(j) + rank_tree(stop).dist;
%fprintf(1,'%d dist = %d\n',stop,rank_tree(stop).dist);
end
end
cycle_length = sum(dist) + 3;
elseif E_cycle(i,1) ~= E_cycle(i,2)
for j=3:4
dist(j) = rank_tree(E_cycle(i,j)).dist;
for k=1:length(rank_tree(E_cycle(i,j)).address);
stop = rank_tree(E_cycle(i,j)).address(k);
dist(j) = dist(j)+rank_tree(stop).dist;
end
end
cycle_length = sum(dist)+1;
else % cycles_through_node==1
continue;
end
for j=1:2
%fprintf('starting branch for E_cycle is %d\n',E_cycle(i,j))
branch_ind = find(conn{node}==E_cycle(i,j),1);
sc = size(cycles{branch_ind});
if sc(1)>0
cycle_ind = find(cycles{branch_ind}(:,1)...
==cycle_length,1);
% cycle_ind could be empty if no cycles of this length yet
else
cycle_ind = [];
end
if isempty(cycle_ind)
cycles{branch_ind} = [cycles{branch_ind};...
cycle_length 1];
else
cycles{branch_ind}(cycle_ind,2) = ...
cycles{branch_ind}(cycle_ind,2)+1;
end
end %j=1:2
end %i=1:sE(1)
% subfunctions --------------------------------------------
function b=branch(cell,rank_tree)
% give branch of cell based on address
if isempty(rank_tree(cell).address)
b = cell;
else
b=rank_tree(cell).address(1);
end