%MATLAB M-FILE FOR READING REACTION NETWORKS
%-------------------------------------------
% Copyright 2005 Dion Vlachos, Kapil Mayawala, Abhijit Chatterjee
% Written by Kapil Mayawala, Department of Chemical Engineering, University of Delaware, Newark, DE 19716, USA
% Tested for MATLAB versions 6.1.0 (R12.1) and 7.0.0 (R14). 
%!!NOTE!!
%--------
%This program is free software; you can redistribute it and/or modify
%it under The terms of the GNU General Public License as published by 
%the Free Software Foundation; either version 2 of the License, or 
%(at your option) any later version.

%This program is distributed in the hope that it will be useful, 
%but WITHOUT ANY WARRANTY; without even the implied warranty of 
%MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
%GNU General Public License for more details.

%You should have received a copy of the GNU General Public
%License along with this program; if not, write to the Free Software
%Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
%USA

%***********************************************************************************************************

% Purpose: Given a reaction mechanism and reaction rate constants in a text file, it generates 
% 5 text files containing stoichiometric data, reaction rate constants, and related data needed for Binomial
% Tau-leap code. ( visit http://dion.che.udel.edu/multiscale/software.html )  

%************************************************************************************************************

% BEFORE YOU USE
% Note:
%      1. Reactions with only one or two species in each direction can be treated. 
%      2. Operators '+' and '=' and the numerical value of the reaction rate constant need to be 
%         separated by at least one space from the previous and the next entry.
%      3. Autocatalytic reactions 
%     		like A + B = 3B 
%	    cannot be entered.
%      4. A Species name cannot be two separate words - for example 
%			carbon + oxygen = carbon dioxide is incorrect
%		but, carbon + oxygen = carbon_dioxide is ok


%************************************************************************************************************

% How to use: When run, a dialog box appears where the user has to enter the input file
%             name, such as reactions.txt, which contains the input as explained below using an example. 

% Example: Consider a reaction mechanism involving 3 reactants (A, B and C), 4 species (A, B, C, and A2), and 3 
%          reactions. The reaction rate constants for the three reactions are 3.1, 5.1 and 2.7. 
%          A + B --> 2C    k=3.1
%          B --> 3A        k=5.1
%          2A --> A2       k=2.7  

%-----------------------------------------------------------------------------------------

% Input file (e.g., reactions.txt):
%          A + B = 2C 3.1
%          B = 3A 5.1
%          2A = A2 2.7

%------------------------------------------------------------------------------------------

% Output files:
%          1. reactants.txt: It assigns numbers to each reactant, i.e.,
%                            1  A
%                            2  B
%                            3  C
%                            4  A2
%          2. stoichiometry.txt: It contains the stoichiometric reaction matrix, i.e.,
%                             -1 -1 2 0
%                             3 -1 0 0
%                             -2 0 0 1
%          3. input_connectivity.txt: The connectivity is generated in the following format
%                             and shows which species are involved in each reaction : 
%                             Reaction_number Connectivity_1 Connectivity_2 Connectivity_3 Connectivity_4
%                i.e,         Reaction_number ReactantNumber1 ReactantNumber2 ReactantNumber3 ReactantNumber4                       
%                i.e.,        1 1 2 3 0
%                             2 2 1 0 0 
%                             3 1 4 0 0
%          4. input_stoic.txt: The stoichiometry of participating species in reactions is generated in the 
%                              following format: 
%                             Reaction_number   Stoich_1   Stoich_2    Stoich_3   Stoich_4
%                i.e.,        1 -1 -1 2 0
%                             2 -1 3 0 0
%                             3 -2 1 0 0 
%          5. input_dyna.txt: The reaction rate constants are generated in the following format: 
%                             Reaction_number   Kinetic_constant
%                i.e.,        1 3.1
%                             2 5.1
%                             3 2.7

%--------------------------------------------------------------------------------------------------

% Chatterjee, A., K. Mayawala, J.S. Edwards and D.G. Vlachos. "Time accelerated Monte Carlo simulations using the binomial t-leap method." Bioinformatics submitted.

%**********************************************************************************************************

clc;
clear;
helpdlg('This program generates files for Binomioal Tau-leap code'); 
prompt  = {'Enter reactions file name:'};
title   = 'Input reaction data';
lines = 1;
inputfile  = inputdlg(prompt,title,lines);
inputfile1=char(inputfile);

[U,B,C,D,P,Q,M,K]=textread(inputfile1,'%s %s %s %s %s %s %s %s');
n=size(U,1)

bb1=size(B,1);
B(bb1+1:n,:)=cellstr(' ');

bb1=size(C,1);
C(bb1+1:n,:)=cellstr(' ');

bb1=size(D,1);
D(bb1+1:n,:)=cellstr(' ');

bb1=size(P,1);
P(bb1+1:n,:)=cellstr(' ');

bb1=size(Q,1);
Q(bb1+1:n,:)=cellstr(' ');

bb1=size(M,1);
M(bb1+1:n,:)=cellstr(' ');

bb1=size(K,1);
K(bb1+1:n,:)=cellstr(' ');

for i=1:1:n
    if strcmp(B(i),'=')
        if (strcmp(D(i),'+') & strcmp(K(i),''))
            K(i)=Q(i);
            M(i)=P(i);
            P(i)=C(i);
            C(i)=cellstr(' ');
        end
    end

    if strcmp(B(i),'+')
        if (strcmp(D(i),'=') & strcmp(K(i),''))
            K(i)=Q(i);
            M(i)=cellstr(' ');
        end
    end

    if strcmp(B(i),'=')
        if strcmp(K(i),'')
            K(i)=D(i);
            M(i)=cellstr(' ');
            P(i)=C(i);
            C(i)=cellstr(' ');
        end
    end
end

K=char(K);
K=str2num(K);

%*************************************************
%coeff separation
U1=char(U);
C1=char(C);
P1=char(P);
M1=char(M);

U2=U1; C2=C1; P2=P1; M2=M1;
l1=length(U1(1,:));
l2=length(C1(1,:));
l3=length(P1(1,:));
l4=length(M1(1,:));
coeff=ones(n,4);

for i=1:1:n
    if l1>0
        d1=str2num(U1(i,1));
        d2=char(d1);
        if ~strcmp(d2,'')
            dummy1=0;
            for j=1:1:l1
                d1=str2num(U1(i,j));
                d2=char(d1);
                if strcmp(d2,'')
                    break
                end
                dummy1=dummy1+1;
            end

            coeff(i,1)=str2num(U1(i,1:dummy1));
            dummy=U1(i,:);
            U1(i,:)=' ';
            for j=1:1:l1-dummy1
                U1(i,j)=dummy(j+dummy1);
            end
        end
    end

    if l2>0
        d1=str2num(C1(i,1));
        d2=char(d1);
        if ~strcmp(d2,'')
            dummy1=0;
            for j=1:1:l2
                d1=str2num(C1(i,j));
                d2=char(d1);
                if strcmp(d2,'')
                    break
                end
                dummy1=dummy1+1;
            end

            coeff(i,2)=str2num(C1(i,1:dummy1));
            dummy=C1(i,:);
            C1(i,:)=' ';
            for j=1:1:l2-dummy1
                C1(i,j)=dummy(j+dummy1);
            end
        end
    end

    if l3>0
        d1=str2num(P1(i,1));
        d2=char(d1);
        if ~strcmp(d2,'')
            dummy1=0;
            for j=1:1:l3
                d1=str2num(P1(i,j));
                d2=char(d1);
                if strcmp(d2,'')
                    break
                end
                dummy1=dummy1+1;
            end
            coeff(i,3)=str2num(P1(i,1:dummy1));
            dummy=P1(i,:);
            P1(i,:)=' ';
            for j=1:1:l3-dummy1
                P1(i,j)=dummy(j+dummy1);
            end
        end
    end

    if l4>0
        d1=str2num(M1(i,1));
        d2=char(d1);
        if ~strcmp(d2,'')
            dummy1=0;
            for j=1:1:l4
                d1=str2num(M1(i,j));
                d2=char(d1);
                if strcmp(d2,'')
                    break
                end
                dummy1=dummy1+1;
            end
            coeff(i,4)=str2num(M1(i,1:dummy1));
            dummy=M1(i,:);
            M1(i,:)=' ';
            for j=1:1:l4-dummy1
                M1(i,j)=dummy(j+dummy1);
            end
        end
    end
end

U=cellstr(U1);
C=cellstr(C1);
P=cellstr(P1);
M=cellstr(M1);

%*************************************************


compds(1)=U(1);
n_compds=1;

for i=1:1:n
    flag=0;
    if l1>0
        if ~strcmp(U(i),'')
            for j=1:1:n_compds
                if strcmp(U(i),compds(j))
                    flag=1;
                end
            end
            if flag==0
                n_compds=n_compds+1;
                compds(n_compds)=U(i);
            end
        end
    end

    flag=0;
    if l2>0
        if ~strcmp(C(i),'')
            for j=1:1:n_compds
                if strcmp(C(i),compds(j))
                    flag=1;
                end
            end
            if flag==0
                n_compds=n_compds+1;
                compds(n_compds)=C(i);
            end
        end
    end

    flag=0;
    if l3>0
        if ~strcmp(P(i),'')
            for j=1:1:n_compds
                if strcmp(P(i),compds(j))
                    flag=1;
                end
            end
            if flag==0
                n_compds=n_compds+1;
                compds(n_compds)=P(i);
            end
        end
    end

    flag=0;
    if l4>0
        if ~strcmp(M(i),'')
            for j=1:1:n_compds
                if strcmp(M(i),compds(j))
                    flag=1;
                end
            end
            if flag==0
                n_compds=n_compds+1;
                compds(n_compds)=M(i);
            end
        end
    end

end


index1=1:1:n;
index1=index1';
index2 = NaN*ones(n,1);
index3 = NaN*ones(n,1);
index4 = NaN*ones(n,1);
index5 = NaN*ones(n,1);

for i=1:1:n
    index2(i)=0; index3(i)=0; index4(i)=0; index5(i)=0;
    for j=1:1:n_compds
        if l2>0
        if strcmp(C(i),compds(j))
            index2(i)=j;
        end
        end
        if l1>0
        if strcmp(U(i),compds(j))
            index3(i)=j;
        end
        end
        if l4>0
        if strcmp(M(i),compds(j))
            index4(i)=j;
        end
        end
        if l3>0
        if strcmp(P(i),compds(j))
            index5(i)=j;
        end
        end
    end
end
m = max([index3; index2; index5; index4]);

S=zeros(n,m);

for i=1:1:n
    if index3(i)==index2(i)
        if index3(i)>0
            S(i,index3(i))=-coeff(i,1)-coeff(i,2);
        end
    else
        if index3(i)>0
            S(i,index3(i))=-coeff(i,1);
        end
        if index2(i)>0
            S(i,index2(i))=-coeff(i,2);
        end
    end
    if index5(i)==index4(i)
        if index5(i)>0
            S(i,index5(i))=coeff(i,3)+coeff(i,4);
        end
    else
        if index5(i)>0
            S(i,index5(i))=coeff(i,3);
        end
        if index4(i)>0
            S(i,index4(i))=coeff(i,4);
        end
    end
end

for i=1:1:n
    if index2(i)==index3(i)
        index2(i)=0;
    end
    if index4(i)==index5(i)
        index4(i)=0;
    end
end

input_connective_data=[index1 index3 index2 index5 index4];
for i=1:1:n
    k1=0; k2=0; k3=0; k4=0;
    if index3(i)~=0
        k1=S(i,index3(i));
    end
    if index2(i)~=0
        k2=S(i,index2(i));
    end
    if index5(i)~=0
        k3=S(i,index5(i));
    end
    if index4(i)~=0
        k4=S(i,index4(i));
    end
    input_stoic_data(i,:)=[i k1 k2 k3 k4];
end

input_dyna_data=[index1 K];
for i=1:1:n
    for j=1:1:4
        if input_connective_data(i,j)==0
            input_connective_data(i,j)=input_connective_data(i,j+1);
            input_stoic_data(i,j)=input_stoic_data(i,j+1);
            input_connective_data(i,j+1)=0;
            input_stoic_data(i,j+1)=0;
        end
    end
end
Stoichiometry=S
input_connective_data
input_stoic_data
input_dyna_data

fid = fopen('reactants.txt','w');
for i=1:1:n_compds
    fprintf(fid,'%i',i);
    g=char(compds(i));
    fprintf(fid,'  %s\n',g);
end
fclose(fid)


dlmwrite('input_connectivity.txt',input_connective_data,' ');
dlmwrite('input_stoic.txt',input_stoic_data,' ');
dlmwrite('input_dyna.txt',input_dyna_data,' ');
dlmwrite('stoichiometry.txt',S,' ');

helpdlg('Output files have been generated: input_connectivity.txt, input_stoic.txt, input_dyna.txt, reactants.txt', 'Program Complete');