find_loops_vset

FIND_LOOPS_VSET Detect feedback loops from sets of variable values of
interest.

 LOOP_REP = FIND_LOOPS_VSET(FUNC, VSET, MAX_NUM_LOOPS) can be used as a wrapper for the
 functions numerical_jacobian_complex() and find_loops(). For the ODE function y'=FUNC(y),
 it returns the loop list at column vector of variable values y=VSET. 
 FUNC should be a function handle that depends on the variable values
 only; parameters and time have to be set to fixed values. 
 MAX_NUM_LOOPS is an optional input parameter that restricts the maximal
 number of reported loops (default: 1e6). LOOP_REP is a cell array of length
 1, LOOP_REP{1} returns the loop list as Matlab Table (formatted as in find_loops()).

0001 function[loop_rep,loop_rep_index,jac_rep,jac_rep_index]=find_loops_vset(func, vset, max_num_loops,compute_full_list)
0002 %FIND_LOOPS_VSET Detect feedback loops from sets of variable values of
0003 %interest.
0004 %
0005 % LOOP_REP = FIND_LOOPS_VSET(FUNC, VSET, MAX_NUM_LOOPS) can be used as a wrapper for the
0006 % functions numerical_jacobian_complex() and find_loops(). For the ODE function y'=FUNC(y),
0007 % it returns the loop list at column vector of variable values y=VSET.
0008 % FUNC should be a function handle that depends on the variable values
0009 % only; parameters and time have to be set to fixed values.
0010 % MAX_NUM_LOOPS is an optional input parameter that restricts the maximal
0011 % number of reported loops (default: 1e6). LOOP_REP is a cell array of length
0012 % 1, LOOP_REP{1} returns the loop list as Matlab Table (formatted as in find_loops()).
0013 
0014 % [LOOP_REP,LOOP_REP_INDEX]= FIND_LOOPS_VSET(FUNC, VSET)
0015 % Multiple sets of variable values can be supplied to the function which reports
0016 % loop lists corresponding to different classes of Jacobian matrices (those
0017 % with similar sign structure). VSET can be, for example, the transposed
0018 % solution vector when solving the ODE system with
0019 % one of the Matlab solvers, or concatenations of other sets of
0020 % variable values. For an ODE system of size _n_, the dimension of VSET should
0021 % be _nxk_ with _k_ being the number of different sets of variable values at
0022 % which loops should be detected (i.e. k column vectors of length n are concatenated).
0023 % LOOP_REP_INDEX has length k (one entry for each set of variable values in VSET)
0024 % and captures as an integer index which loop list corresponds to which set
0025 % of variable values. Please note that loop lists can be identical despite
0026 % belonging to different sign classes of Jacobians.
0027 %
0028 % [LOOP_REP,LOOP_REP_INDEX]= FIND_LOOPS_VSET(FUNC, VSET,MAX_NUM_LOOPS,COMPUTE_FULL_LIST)
0029 % Syntax as before. If the logical input COMPUTE_FULL_LIST is set to false (default: true),
0030 % the classes of Jacobians are further reduced and loop lists are not re-computed
0031 % for Jacobians that clearly do not allow for altered loop lists. This is
0032 % the case if no new regulation appear and
0033 % only signs of regulations are altered that are not member of any loop.
0034 % Loop lists can still be identical (e.g. if two sign
0035 % switches occur that are both affecting the same loops).
0036 %
0037 % [LOOP_REP,LOOP_REP_INDEX,JAC_REP,JAC_REP_INDEX]= FIND_LOOPS_VSET(_)
0038 % Syntax as before. JAC_REP returns the different detected signed Jacobians,
0039 % JAC_REP_INDEX gives the index of the signed Jacobian for each set of
0040 % variable values. Only if COMPUTE_FULL_LIST is set to false, LOOP_REP can
0041 % contain
0042 % fewer elements than JAC_REP, otherwise both have the same number of elements.
0043 %
0044 % Examples
0045 % Determine the loop list for a function at one set of variable values
0046 %   klin=[165,0.044,0.27,550,5000,78,4.4,5.1];
0047 %   knonlin=[0.3,2];
0048 %   loop_list = find_loops_vset(@(x)func_POSm4(0,x,klin,knonlin),[1 2 1 1])
0049 %
0050 % Determine the loop list along the solution of an ODE system, report at
0051 % most 100000 loops and do not compute the loop list for all Jacobians.
0052 %   sol=dlmread('li08_solution.tsv','delimiter','\t');
0053 %   [loop_rep,loop_rep_index,jac_rep,jac_rep_index]=find_loops_vset(...
0054 %       @(x)func(0,x), sol(2:end,:), 1e5, false)
0055 %
0056 % See also: FIND_LOOPS(), NUMERICAL_JACOBIAN_COMPLEX()
0057 
0058 
0059 
0060 if nargin<3
0061     max_num_loops=1e6;
0062 end
0063 if nargin<4
0064     compute_full_list=true;
0065 end
0066 
0067 %determine the initial Jacobian and loop list
0068 J_ini=numerical_jacobian_complex(func,vset(:,1));
0069 
0070 vset_count=size(vset,2); %how many different sets of variable values
0071 
0072 %determine all different Jacobian matrices (signs only)
0073 Jsig_tab=zeros(1,size(vset,2));
0074 Jsig_rep={sign(J_ini)};
0075 for i=1:vset_count
0076     J_temp=numerical_jacobian_complex(func,vset(:,i),1e-8);
0077     %check whether the new Jacobian is different from the Jacobians found earlier
0078     for j=1:length(Jsig_rep) 
0079         disttemp=max(max(abs(sign(J_temp)-Jsig_rep{j})));
0080         if disttemp==0 %if another Jacobian is the same
0081             Jsig_tab(i)=j;
0082             break
0083         end
0084     end
0085     if Jsig_tab(i)==0 %if no other Jacobian was the same before
0086         Jsig_rep{end+1}=sign(J_temp); %save the Jacobian
0087         Jsig_tab(i)=length(Jsig_rep);
0088     end         
0089 end
0090 jac_rep=Jsig_rep;
0091 jac_rep_index=Jsig_tab;
0092 
0093 %determine loops for the different Jacobians
0094 %loops should only differ if
0095 % - switch from zero no nonzero
0096 % - switch of an entry that is an edge of a loop
0097 
0098 loop_tab=zeros(1,size(vset,2));
0099 loop_tab(Jsig_tab==1)=1; %set the loop index according to the Jacobian that belongs
0100 loop_rep={find_loops(J_ini,max_num_loops)}; %J_ini equals Jsig_rep{1}
0101 
0102 if compute_full_list %if we should compute the full loop list
0103     
0104     for i=2:length(Jsig_rep) %check each Jacobian (if there are more than 1)
0105         J_temp=Jsig_rep{i};
0106         loop_rep{end+1}=find_loops(J_temp,max_num_loops);
0107     end
0108     loop_tab=jac_rep_index; %loop identity equals Jacobian identity
0109 
0110 else %if we should not compute the full loop list, more efficient (exclude Jacobians that cannot generate different loop lists)
0111   
0112     for i=2:length(Jsig_rep) %check each Jacobian (if there are more than 1)
0113         J_temp=Jsig_rep{i};
0114         %determine if there is a switch from zero to nonzero
0115         switch_to_nonzero=zeros(1,i-1);
0116         for j=1:(i-1) %compare with all earlier Jacobians
0117             switch_to_nonzero(j)=sum(sum(abs(J_temp(Jsig_rep{j}==0))))>0;
0118         end
0119         %determine if the sign switches of an edge that is contained in a loop
0120         change_in_a_loop=zeros(1,i-1);
0121         for j=1:(i-1) %compare with all earlier Jacobians
0122             if switch_to_nonzero(j)==0 %we only have to check this if we don't already have another change
0123                 [target_node_ind, source_node_ind]=find(not(J_temp-Jsig_rep{j}==0));
0124                 for k=1:length(target_node_ind)
0125                     change_in_a_loop(j)=not(isempty(find_edge(loop_rep{loop_tab(find(Jsig_tab==j,1))},... %the loop list corresponding to Jacobian j
0126                         source_node_ind(k),target_node_ind(k))));
0127                     if change_in_a_loop(j)==1 %in case there is a change in at least one edge
0128                         break
0129                     end
0130                 end
0131             end
0132         end
0133         %check if there is any earlier Jacobian that is neither changed from zero to
0134         %nonzero for J_temp, nor any change in an edge of a loop is observed --> this
0135         %Jacobian will lead to exactly the same loops as the Jacobian of
0136         %interest, J_temp
0137         loop_ind_temp=find(switch_to_nonzero==0 & change_in_a_loop==0,1); %find only the first Jacobian that has the same profile
0138         if isempty(loop_ind_temp) %in case no Jacobian coincides, we compute the loops anew
0139             loop_rep{end+1}=find_loops(J_temp,max_num_loops);
0140             loop_tab(Jsig_tab==i)=length(loop_rep);
0141         else %in case the Jacobian is similar to another Jacobian
0142             %(e.g. because the sign switch appeared in a nonzero entry that
0143             %does not generate any loop) we just use that Jacobian
0144             loop_tab(Jsig_tab==i)=loop_ind_temp;
0145         end   
0146     end
0147 end
0148 
0149 loop_rep_index=loop_tab;
0150 
0151 end
0152