find_loops

FIND_LOOPS Compute feedback loops from Jacobian or adjacency matrix.

 loop_tab = FIND_LOOPS(A_in, max_num_loops) returns a Matlab
 table with three columns, 'loop', 'length', 'sign', that represents 
 feedback loops as derived from matrix A_in.
 'loop' is a cell array returning the regulations, e.g. the output 
 1 2 3 1 means that variable 1 regulates variable 2, variable 2 regulates 
 variable 3 and variable 3 regulates variable 1.
 'length' is a positive integer indicating the number of variables forming
 the loop, i.e. the loop length.
 'sign' is a value in {-1,1} and indicates whether the feedback loop is 
 negative or positive, i.e. whether the number of negative interactions
 forming the loop is uneven or even.
 A_in is a (Jacobian) nxn-matrix that captures the regulations of n 
 variables on each other. Thereby, entry a_ij represents the regulation of 
 variable i by variable j. a_ij should be non-zero only if variable j 
 regulates variable i. The effect of the regulation (inhibition vs.
 activation) should be encoded in the sign of a_ij. 
 max_num_loops is an optional parameter indicating the maximal number of
 feedback loops that are detected and returned. To limit runtime, it 
 defaults to 1e6.
 
 FIND_LOOPS uses the decomposition into strongly connected components 
 by Tarjan's algorithm (function tarjan.m from Matlab File Exchange server)
 and furthermore relies on Johnson's algorithm as adapted from 
 find_elem_circuits.m by cmaes on github.

 Examples
 System of 4 variables with 4 self-loops (loops of length 1), all
 negative, 1 positive loop of length 3 and one negative loop of length 4
   J = [-1 0 0 -1; 1 -1 0 1; 0 1 -1 0; 0 0 1 -1];
   loop_t = find_loops(J)
   loop_t.loop{1} % delivers the first loop
 Limiting the output to 4 loops:
   loop_t = find_loops(J,4)

 $Author: kabaum $  $Date: 2020/04/06 09:46 $ $Revision: 0.1 $
 Copyright: GNU GPLv3 K. Baum 2020

 See also: FIND_LOOPS_NOSCC(), find_loops_vset()

0001 function [loop_tab] = find_loops(A_in,max_num_loops)
0002 %FIND_LOOPS Compute feedback loops from Jacobian or adjacency matrix.
0003 %
0004 % loop_tab = FIND_LOOPS(A_in, max_num_loops) returns a Matlab
0005 % table with three columns, 'loop', 'length', 'sign', that represents
0006 % feedback loops as derived from matrix A_in.
0007 % 'loop' is a cell array returning the regulations, e.g. the output
0008 % 1 2 3 1 means that variable 1 regulates variable 2, variable 2 regulates
0009 % variable 3 and variable 3 regulates variable 1.
0010 % 'length' is a positive integer indicating the number of variables forming
0011 % the loop, i.e. the loop length.
0012 % 'sign' is a value in {-1,1} and indicates whether the feedback loop is
0013 % negative or positive, i.e. whether the number of negative interactions
0014 % forming the loop is uneven or even.
0015 % A_in is a (Jacobian) nxn-matrix that captures the regulations of n
0016 % variables on each other. Thereby, entry a_ij represents the regulation of
0017 % variable i by variable j. a_ij should be non-zero only if variable j
0018 % regulates variable i. The effect of the regulation (inhibition vs.
0019 % activation) should be encoded in the sign of a_ij.
0020 % max_num_loops is an optional parameter indicating the maximal number of
0021 % feedback loops that are detected and returned. To limit runtime, it
0022 % defaults to 1e6.
0023 %
0024 % FIND_LOOPS uses the decomposition into strongly connected components
0025 % by Tarjan's algorithm (function tarjan.m from Matlab File Exchange server)
0026 % and furthermore relies on Johnson's algorithm as adapted from
0027 % find_elem_circuits.m by cmaes on github.
0028 %
0029 % Examples
0030 % System of 4 variables with 4 self-loops (loops of length 1), all
0031 % negative, 1 positive loop of length 3 and one negative loop of length 4
0032 %   J = [-1 0 0 -1; 1 -1 0 1; 0 1 -1 0; 0 0 1 -1];
0033 %   loop_t = find_loops(J)
0034 %   loop_t.loop{1} % delivers the first loop
0035 % Limiting the output to 4 loops:
0036 %   loop_t = find_loops(J,4)
0037 %
0038 % $Author: kabaum $  $Date: 2020/04/06 09:46 $ $Revision: 0.1 $
0039 % Copyright: GNU GPLv3 K. Baum 2020
0040 %
0041 % See also: FIND_LOOPS_NOSCC(), find_loops_vset()
0042 
0043 
0044 %set maximal number of reported loops (plus self-loops)
0045 if nargin < 2
0046     max_num_loops = 1e6;
0047 end
0048 
0049 %transpose (required for correct direction detection) and turn into sparse
0050 %representation
0051 %i.e. entry A_in(2,3) determines how species 2 is regulated by species 3
0052 if ~issparse(A_in)
0053     A = sparse(A_in'); %this command alters the format of how the matrix is represented
0054 else
0055     A=A_in';
0056 end
0057 
0058 %preparatory functions for the matrix
0059 
0060     function[self_loop_tab,A_red] = determine_and_remove_self_loops(C)
0061         %determine self-loops, output as table
0062         self_loops=sign(diag(C));
0063         self_loops_ind=find(~(self_loops==0));
0064         self_loops=self_loops(self_loops_ind);
0065         A_red=C-diag(diag(C)); %remove diagonal entries
0066         self_loop_tab=table(... %write the input into result table
0067             mat2cell(cat(2,self_loops_ind,self_loops_ind),ones(length(self_loops_ind),1)),...
0068             ones(length(self_loops_ind),1),...
0069             self_loops,...
0070             'VariableNames',{'loop','length','sign'});
0071     end
0072 
0073 %detect nonviable_nodes and set the matrix rows & columns to zero
0074     function[A_red]=set_nonviable_nodes_to_zero(C)
0075         %determine those variables that have no ingoing or no outgoing
0076         %edges!
0077         %assume that the diagonal entries have been set to zero already!
0078         nonviable_nodes=find(any(C)+any(C')<2);
0079         A_red=C;
0080         A_red(nonviable_nodes,:)=0; %set nonviable rows to zero
0081         A_red(:,nonviable_nodes)=0; %set nonviable columns to zero
0082     end
0083 
0084 %preprocess A
0085 [self_loop_tab,A] = determine_and_remove_self_loops(A); %determine self-loops and eradicate diag entries
0086 [A]=set_nonviable_nodes_to_zero(A);
0087 
0088 %now advance with determining strongly connected components
0089 %function from Matlab Exchange, see licence
0090 [comp_ids]=tarjan(A); %
0091 
0092 %size of largest connected component will be used to define the blocking
0093 %structures
0094 max_comp_size=length(comp_ids{1}); %they are sorted by descending size, this is therefore the maximal size
0095 
0096 %functions for Johnson's algorithm (cmaes!) of simple circuit detection
0097 %require global variables stack, blocked, Blist, cycles
0098 n = max_comp_size;
0099 Blist = cell(n,1); %this has something to do with blocking nodes
0100 blocked = false(1,n); %whether nodes are blocked or not
0101 stack=[]; %important for keeping track
0102 cycles = {}; %save the cycles in there, the function circuit alters this object
0103 numcycles=size(self_loop_tab,1); %keep track of found number of cycles, add detected self loops!
0104 
0105     function unblock(u)
0106         blocked(u) = false;
0107         for w=Blist{u}
0108             if blocked(w)
0109                 unblock(w)
0110             end
0111         end
0112         Blist{u} = [];
0113     end
0114 
0115     function f = circuit(v, s, C,ind_transfer,num_loops_max) %this function seems to detect paths from v to s in matrix C, saves them to "cycle"
0116         f = false;
0117         
0118         stack(end+1) = ind_transfer(v);
0119         blocked(v) = true;
0120         
0121         for w=find(C(v,:)) %determine non-zero elements in row v of C
0122             if w == s
0123                 cycles{end+1} = [stack ind_transfer(s)];
0124                 numcycles=numcycles+1; %keep track of number of cycles
0125                 f = true;
0126             elseif ~blocked(w)
0127                 if circuit(w, s, C, ind_transfer,num_loops_max) % if there is a path from w to s
0128                     f = true;
0129                 end
0130             end
0131             if numcycles>=num_loops_max %leave function if we have enough FBLs detected
0132                 return
0133             end
0134         end
0135         
0136         if f
0137             unblock(v)
0138         else
0139             for w = find(C(v,:))
0140                 if ~ismember(v, Blist{w})
0141                     Bnode = Blist{w};
0142                     Blist{w} = [Bnode v];
0143                 end
0144             end
0145         end
0146         
0147         stack(end) = [];
0148     end
0149 
0150 
0151 %go through every component of length larger than 1 - they are sorted in
0152 %descending order of their sizes!
0153 i=1;
0154 while i<=length(comp_ids)&&length(comp_ids{i})>1  
0155     
0156     A_temp=A(comp_ids{i},comp_ids{i}); %submatrix of the strongest connected component variables only
0157     
0158     s = 2; %start at node 2; we might go for starting at node n and revert the direction?
0159     while s <= size(A_temp,1) %go from s=2 up to length of A_temp
0160         
0161         %we could also go and find connected components for each of these
0162         %subgraphs! try in terms of runtime...
0163     
0164         % Subgraph of G induced by {1,...,s,}
0165         F = A_temp(1:s,1:s);
0166         blocked(1:s) = false;
0167         Blist(1:s) = cell(s,1);
0168         circuit(s, s, F,comp_ids{i},max_num_loops);
0169         if numcycles >= max_num_loops %leave the loop if enough cycles have been detected
0170             break
0171         end
0172         s = s + 1; %go one species more
0173 
0174     end
0175     if numcycles >= max_num_loops %leave the loop if enough cycles have been detected
0176         break
0177     end
0178 i=i+1; %go one connected component further
0179 end
0180 
0181 %now we have all cycles or the maximally allowed number
0182 %start preparing the table
0183 
0184 %determine cycle signs
0185 
0186 loop_length=cellfun(@length,cycles)-1;
0187 function[sign_val]=loop_sign_func(cycinds, C)
0188     %determine sign from indices
0189     sign_val= mod(sum(sign(C(sub2ind(size(C),cycinds(1:(end-1)),cycinds(2:end))))<1),2)...
0190         *(-2)+1;
0191 end
0192 loop_sign=cellfun(@(x)loop_sign_func(x,A),cycles);
0193 
0194 %cast the loop table and add self-loops such that at most max_num_loops are
0195 %returned
0196 loop_tab=cat(1,table(cycles',loop_length',loop_sign','VariableNames',{'loop','length','sign'}),self_loop_tab(1:min(height(self_loop_tab),max_num_loops-numcycles),:));
0197 
0198 end
0199         
0200