find_loops_noscc

FIND_LOOPS_NOSCC Compute feedback loops from Jacobian or adjacency matrix.

 loop_tab = FIND_LOOPS_NOSCC(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, and the direction of the regulation (inhibition vs.
 activation) and 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. It defaults to 1e6.
 Note that the number of reported loops might be larger if the number of
 self-loops is larger.

 In contrast to FIND_LOOPS(), FIND_LOOPS_NOSCC does not use decomposition 
 into strongly connected components and exclusively relies on Johnson's 
 algorithm as adapted from find_elem_circuits by cmaes on git.

 Example
 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_noscc(J)
   loop_t.loop{1} % delivers the first loop

 $Author: kabaum $  $Date: 2020/06/08 10:58 $ $Revision: 0.1 $
 Copyright: GNU GPLv3 K. Baum 2020

 See also: FIND_LOOPS

0001 function [loop_tab] = find_loops_noscc(A_in,max_num_loops)
0002 %FIND_LOOPS_NOSCC Compute feedback loops from Jacobian or adjacency matrix.
0003 %
0004 % loop_tab = FIND_LOOPS_NOSCC(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, and the direction of the regulation (inhibition vs.
0019 % activation) and 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. It defaults to 1e6.
0022 % Note that the number of reported loops might be larger if the number of
0023 % self-loops is larger.
0024 %
0025 % In contrast to FIND_LOOPS(), FIND_LOOPS_NOSCC does not use decomposition
0026 % into strongly connected components and exclusively relies on Johnson's
0027 % algorithm as adapted from find_elem_circuits by cmaes on git.
0028 %
0029 % Example
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_noscc(J)
0034 %   loop_t.loop{1} % delivers the first loop
0035 %
0036 % $Author: kabaum $  $Date: 2020/06/08 10:58 $ $Revision: 0.1 $
0037 % Copyright: GNU GPLv3 K. Baum 2020
0038 %
0039 % See also: FIND_LOOPS
0040 
0041 
0042 %set maximal number of reported loops (plus self-loops)
0043 if nargin < 2
0044     max_num_loops = 1e6;
0045 end
0046 
0047 %transpose (required for correct direction detection) and turn into sparse
0048 %representation
0049 %i.e. entry A_in(2,3) determines how species 2 is regulated by species 3
0050 if ~issparse(A_in)
0051     A = sparse(A_in');
0052 else
0053     A = A_in';
0054 end
0055 
0056 %preparatory functions for the matrix
0057 
0058 %determine self-loops and remove entries on the diagonal
0059     function[self_loop_tab,A_red] = determine_and_remove_self_loops(C)
0060         %determine self-loops, output as table
0061         self_loops=sign(diag(C));
0062         self_loops_ind=find(~(self_loops==0));
0063         self_loops=self_loops(self_loops_ind);
0064         A_red=C-diag(diag(C)); %remove diagonal entries
0065         self_loop_tab=table(... %write the input into result table
0066             mat2cell(cat(2,self_loops_ind,self_loops_ind),ones(length(self_loops_ind),1)),...
0067             ones(length(self_loops_ind),1),...
0068             self_loops,...
0069             'VariableNames',{'loop','length','sign'});
0070     end
0071 
0072 %detect nonviable nodes (i.e. sinks or sources of the network) and set the
0073 %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; assumes that the diagonal entries have been set to zero
0077         nonviable_nodes=find(any(C)+any(C')<2);
0078         A_red=C;
0079         A_red(nonviable_nodes,:)=0; %set nonviable rows to zero
0080         A_red(:,nonviable_nodes)=0; %set nonviable columns to zero
0081     end
0082 
0083 %preprocess A
0084 [self_loop_tab,A] = determine_and_remove_self_loops(A); %determine self-loops and eradicate diag entries
0085 [A]=set_nonviable_nodes_to_zero(A);
0086 
0087 %functions for Johnson's algorithm (cmaes!) of simple circuit detection
0088 %require global variables stack, blocked, Blist, cycles
0089 
0090 n = size(A,1);
0091 Blist = cell(n,1); %this has something to do with blocking nodes
0092 blocked = false(1,n); %whether nodes are blocked or not
0093 stack=[]; %important for keeping track
0094 cycles = {}; %save the cycles in there, the function circuit alters this object
0095 numcycles=size(self_loop_tab,1); %count number of detected FBLs
0096 
0097     function unblock(u)
0098         blocked(u) = false;
0099         for w=Blist{u}
0100             if blocked(w)
0101                 unblock(w)
0102             end
0103         end
0104         Blist{u} = [];
0105     end
0106 
0107     function f = circuit(v, s, C,num_loops_max) %this function seems to detect paths from v to s in matrix C, saves them to "cycle"
0108         f = false;
0109         
0110         stack(end+1) = v;
0111         blocked(v) = true;
0112         
0113         for w=find(C(v,:)) %determine non-zero elements in row v of C
0114             if w == s
0115                 cycles{end+1} = [stack s];
0116                 numcycles=numcycles+1;
0117                 f = true;
0118             elseif ~blocked(w)
0119                 if circuit(w, s, C,num_loops_max) % if there is a path from w to s
0120                     f = true;
0121                 end
0122             end
0123             if numcycles>=num_loops_max %leave function if we have enough loops detected
0124                 return
0125             end
0126         end
0127              
0128         if f
0129             unblock(v)
0130         else
0131             for w = find(C(v,:))
0132                 if ~ismember(v, Blist{w})
0133                     Bnode = Blist{w};
0134                     Blist{w} = [Bnode v];
0135                 end
0136             end
0137         end
0138         
0139         stack(end) = [];
0140     end
0141    
0142 %go through every index
0143 s = 2; %start at node 2; we might go for starting at node n and revert the direction?
0144     while s <= n %go from s=2 up to n =length of A
0145     
0146         % Subgraph of G induced by {1,...,s,}
0147         F = A(1:s,1:s);
0148         blocked(1:s) = false;
0149         Blist(1:s) = cell(s,1);
0150         circuit(s, s, F,max_num_loops);
0151         s = s + 1; %go one species more
0152         if numcycles>=max_num_loops %leave the loop if we have enough FBLs detected
0153             break 
0154         end
0155 
0156     end
0157 
0158 %now we have all cycles or the maximally allowed number
0159 %start preparing the table
0160 
0161 %determine cycle signs
0162 
0163 loop_length=cellfun(@length,cycles)-1;
0164 function[sign_val]=loop_sign_func(cycinds, C)
0165     %determine sign from indices
0166     sign_val= mod(sum(sign(C(sub2ind(size(C),cycinds(1:(end-1)),cycinds(2:end))))<1),2)...
0167         *(-2)+1;
0168 end
0169 loop_sign=cellfun(@(x)loop_sign_func(x,A),cycles);
0170 
0171 %cast the loop table and add self-loops such that at most max_num_loops are
0172 %returned
0173 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),:));
0174 end
0175         
0176