example.m

%%%% In this script we provide an example of how to calculate the SVD 
%%%% of the Reaction-Diffusion Stoichiometry matrix as discussed in
%%%% the manuscript:
%%%%
%%%% J. Wentz, D. Bortz, "Analytical singular value decomposition for a 
%%%% class of stoichiometry matrices" (2021).

addpath('helperFunctions')

%%  Define stoichiometry matrix
%Note that this is the complete stoichoiemtry matrix that contains all 
%reactions (in both subregion 1 and subregion 2). In the next code block we
%select the relevent columns of this matrix for each of the subregions.

Sr = [0 0 -1;
      -1 -1 0; 
      -1 0 1;        
      2 1 0;
      0 -1 0;
      0 1 0];

%% Set up (Definitions)

%%% USER DEFINED PARAMETERS %%%
n = 8;                  %Total number of spatial compartments
n1 = 2;                 %Number of compartments to left of barrier
BC = 'Mixed-ND';        %Boundary condition (Neumann, Mixed-ND, or Dirichlet)
gamma = .1;             %Rate of reactions to diffusion
M_p = [1,2,4,6];        %Indices of species that diffuse across barrier
M_m = [3,5];            %Indices of species that cannot diffuse across barrier
R1 = [1,3];             %Reactions that occur in first subregion
R2 = [1,2,3];           %Reactions that occur in second subregion

%%% AUTOMATICALLY DEFINED PARAMETERS %%%
Sr1 = Sr(:,R1);
Sr2 = Sr(:,R2);
n2 = n - n1;            %Number of compartments to the right of barrier
m_plus = length(M_p);   %Number of spcies that can diffuse across barrier
m_minus = length(M_m);  %Number of species that can't diffuse across barrier
m=size(Sr1,1);
p1=size(Sr1,2);
p2=size(Sr2,2);

%% MAIN RESULT: Calculate the approximate Singular Value Decomposition
%Here we provide an example of how use the function svdWithBarrier to
%calculate the SVD. The order of the singular vectors/values
%matches the ordering in the paper. Note that V_null_a is the null space 
%as defined in Thm 2.1 and V_null_b is the null space as defined in Prop
%5.4.

[U, U_null, V, V_null_a, V_null_b, Sigma] = svdWithBarrier(gamma,Sr1,Sr2,n,n1,BC,M_p,M_m);

%% Calculate the full RD stoichiometry matrix
%Here we provide code for directly calculating the RD stoichiometry matrix.

%Get diffusion stoichiometry matrix
Sd = stoicDiffusion(n,BC);

%Define D_plus, D_minus, and H
D_p_vec = zeros(m,1);
D_p_vec(M_p)=1;
D_p = diag(D_p_vec);

D_m_vec = zeros(m,1);
D_m_vec(M_m)=1;
D_m = diag(D_m_vec);

H = zeros(n,n+1);
H(n1,n1+1)=-1;
H(n1+1,n1+1)=1;

%Calculate the stoichiometry matrix
S = [gamma*blkdiag(kron(eye(n1),Sr1),kron(eye(n-n1),Sr2)), kron(Sd,D_p) + kron(Sd-H,D_m)];


%% TESTING - Make sure the expected result is obtained as gamma -> 0

% %Should approach 0 as gamma goes to 0
% max(max(S - U*Sigma*V'))
% 
% max(max(S'*U_null))  %Always orthonormal null space! (for all gamma)
% max(max(S*V_null_a)) %Approaches zero as gamma -> 0
% max(max(S*V_null_b)) %Always null space (not always orthonormal)
% 
% %Should always be orthogonal as gamma -> 0
% max(max(eye(size(S,2))-[V V_null_a]'*[V V_null_a]))
% max(max(eye(size(S,2))-[V V_null_b]'*[V V_null_b]))
% 
% %Should always by orthogonal
% max(max(eye(size(S,1))-[U U_null]'*[U U_null]))
% max(max(eye(rank(null(S)))-V_null_a'*V_null_a))


%% Reproduce error analysis discussed in Section 5.6

% for gamma = [.1,.01,.001]
%     [U, U_null, V, V_null_a, V_null_b, Sigma] = svdWithBarrier(gamma,Sr1,Sr2,n,n1,BC,M_p,M_m); 
%     S = [gamma*blkdiag(kron(eye(n1),Sr1),kron(eye(n-n1),Sr2)), kron(Sd,D_p) + kron(Sd-H,D_m)];       
%     [U2, Sigma2, V2] = svd(S);
%     q = rank(S);
%     SigmaErr = zeros(q,1);
%     Uerr = zeros(q,1);
%     Verr = zeros(q,1);
%     evec_pert_sorted = [];
% 
%     for i = 1:q
%         tol = 1*10^-13;
%         ii = find(ismembertol(diag(Sigma),real(Sigma2(i,i)),tol));
%         while isempty(ii)
%             disp('Empty, increasing tol...')
%             tol = tol + tol*.01;
%             ii = find(ismembertol(diag(Sigma),Sigma2(i,i),tol));    
%         end
%         if length(ii) > 1
%             disp('Found more than one.')
%         else
%             k = find(abs(U(:,ii))>10^-6,1);
%             SigmaErr(i) = abs(Sigma(ii,ii)-Sigma2(i,i));
%             Uerr(i) = sqrt(sum((sign(U(k,ii))*U(:,ii)-sign(U2(k,i))*U2(:,i)).^2));
%             Verr(i) = sqrt(sum((sign(U(k,ii))*V(:,ii)-sign(U2(k,i))*V2(:,i)).^2));
%         end
%     end
% 
%     Fig1 = figure(1);
%     set(Fig1,'defaulttextinterpreter','latex');
%     subplot(3,1,1)
%     semilogy(1:q,SigmaErr)
%     hold on
%     xlabel('Index of singular value ($i$)')
%     ylabel('$|\sigma^{(i)}-\sigma_{num}^{(i)}|$')
% 
%     subplot(3,1,2)
%     semilogy(1:q,Uerr)
%     hold on
%     xlabel('Index of left singular vector ($i$)')
%     ylabel('$\|u^{(i)}-u_{num}^{(i)}\|_2$')
% 
%     subplot(3,1,3)
%     semilogy(1:q,Verr)
%     hold on
%     xlabel('Index of right singular vector ($i$)')
%     ylabel('$\|v^{(i)}-v_{num}^{(i)}\|_2$')
% end
% 
%     
% lgd = legend('0.1','0.01','0.001','Orientation','horizontal','Position',[.47,.96,.1,.01]);
% set(Fig1,'Units','Inches');
% Fig1.Position = [1 1 5 7];
% pos = get(Fig1,'Position');
% set(Fig1,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[pos(3), pos(4)])
% print(Fig1,'fig_error','-dpdf','-r0')