

rxn_groupnumbers = [];
clear function rxn_prob_undef

for i = 1:num_rxns
    eval(['rxn_groupnumbers = union(' rxn_list{i}{1} '.rxn_grp, rxn_groupnumbers);'])
end
    

num_rxn_grps = length(rxn_groupnumbers);
undef_rxn = 0;

for i = rxn_groupnumbers
    
    rxnnum_grp(i) = 0;
    rxn_grp{i}.comps = [];
    rxn_grp{i}.global_comp_gid = [];
    
end
    
for j = 1:num_rxns
    
    clear tmp;
    num_new_rxns = 0;
    tmp_depend = [];
    
    eval(['this_rxn = ' rxn_list{j}{1} ';'])
    
    group_num = this_rxn.rxn_grp;
    if ~isfield(this_rxn, 'type')
        this_rxn.type = 'undef';
    end
    rxntype = this_rxn.type;
    this_rxn.name = rxn_list{j}{1};

    for k = group_num       % provides an easy way to put the same reaction into multiple reaction groups, which means you really have multiple reactions.
           switch rxntype
            
            case {'uni-uni'}
                [tmp num_new_rxns rxn_grp{k}] = prepare_uni_rxn(this_rxn, rxn_grp{k}, cmpd_info.cmpd_name, 1);
                
            case {'uni-multi'}
                [tmp num_new_rxns rxn_grp{k}] = prepare_uni_rxn(this_rxn, rxn_grp{k}, cmpd_info.cmpd_name, 2);
                
%            case {'homo-dimerization'}
%                [tmp num_new_rxns rxn_grp{k}] = prepare_homo_rxn(this_rxn, rxn_grp{k}, cmpd_info.cmpd_name);
                
            case {'binding'}
                [tmp num_new_rxns rxn_grp{k}] = prepare_binding_rxn(this_rxn, rxn_grp{k}, cmpd_info.cmpd_name);
                
            case {'nth-order'}
	%            [tmp num_new_rxns rxn_grp{i}.comps] = prepare_bidirectional_rxn(filename,rxn_grp{i}.comps);
                [tmp num_new_rxns rxn_grp{k}] = prepare_nth_order_rxn(this_rxn, rxn_grp{k}, cmpd_info.cmpd_name);
                
            case {'transv1'}
                [tmp num_new_rxns rxn_grp{k}] = prepare_transv1_rxn(this_rxn, rxn_grp{k}, cmpd_info.cmpd_name);
                
            case {'generalized-mich-ment'}
                [tmp num_new_rxns rxn_grp{k}] = prepare_genmm_rxn(this_rxn, rxn_grp{k}, cmpd_info.cmpd_name);
            
            case {'undef'}
                [tmp num_new_rxns rxn_grp{k}] = prepare_undef_rxn(this_rxn, rxn_grp{k}, cmpd_info.cmpd_name,volume);
                undef_rxn = 1;
 
            case {'transcv2'}
                [tmp num_new_rxns rxn_grp{k} cmpd_info.cmpd_name cmpd_info.orig_cmpd_vals tmp_depend tmp_min_depend num_cmpds] = prepare_transcv2_rxn(this_rxn, rxn_grp{k}, cmpd_info.cmpd_name, cmpd_info.orig_cmpd_vals);
                cmpd_info.num_cmpds = num_cmpds;
                
            otherwise
                error(['Reaction type ' rxntype ' is not recognized.  Exiting program.'])
                   
            end
     
        for l = 1:num_new_rxns
            rxnnum_grp(k) = rxnnum_grp(k) + 1;
            tmp{l}.id = j;
            if num_new_rxns > 1
                tmp{l}.name = [rxn_list{j}{1} '_' num2str(l)];
            else
                tmp{l}.name = rxn_list{j}{1};
            end
            rxn_grp{k}.rxn{rxnnum_grp(k)} = tmp{l};
            if isempty(tmp_depend)
                rxn_grp{k}.rxn_dependency_matrix{rxnnum_grp(k)} = [];
                min_depend{k}(rxnnum_grp(k)) = 0;
            else
                rxn_grp{k}.rxn_dependency_matrix{rxnnum_grp(k)} = tmp_depend{l};
                min_depend{k}(rxnnum_grp(k)) = tmp_min_depend(l);
            end
        end
        
    end
end

if undef_rxn == 1
% Create a c-file with Kinetic Law for "undef" reactions
    temp_file = 'undef_prob';
    fid = fopen(temp_file, 'w');
    if fid == -1
            errordlg('Cannot open temporary undef_prob file','Open File Error!')
            return
    end
    
  %  Need to define variables in c-file only once
      def = zeros(num_cmpds,1);
 
    fprintf(fid,'#include \"runsim.h\"\n\n');
 
    fprintf(fid,'double rxn_prob_undef(\n  struct rxn *arxn,\n  struct cmpd_data *cmpds  )\n{\n');
    fprintf(fid,'/*  Propensity functions for reactions with defined Kinetic Law\n');
    fprintf(fid,'**   Parameters are assigned in prepare_undef_reactions\n');
    fprintf(fid,'**   Need to assign the number of compounds and kinetic constants  \n*/\n');
    fprintf(fid,'     double ai=0;\n');
 

    is = 1;
    out_str{is} = sprintf('     switch(arxn->id)\n\t{\n');
    is=is+1;

    
    for i=1:num_rxns
        clear rxn
        
        rxn = rxn_grp{1}.rxn{i};
        if strcmp(rxn.type,'undef')
        out_str{is} = sprintf('       case %u:\n',i-1);
        is=is+1;
        formula = strrep(rxn.formula, ' ', '');

 %  Translate a Kinetic Law formula       
     
 %  1. Extract all variables
     
      s = symvar(formula);
 
 %  Assign values for parameters and volume  
   n_kin_const = size(rxn.StoP_name,2);
        if size(s) > 0
            for i = 1:size(s)
                temp_str = [];
                found = 0;
                % Check whether it is a rate        
                 for j = 1:n_kin_const
                     if strcmp(s{i},rxn.StoP_name{j})
                                             % Replace rate with it's StoP name
                        num_str =sprintf('arxn->rates[%u].value',j-1); 
                        formula = strrep(formula, s{i},num_str);
                        found = 1;
                        break;
                     end
                 end
                  if found == 0
                    if strncmp(s{i}, 'vol',3)    % Check whether it is a system's volume
                        num_str =num2str(volume,6); 
                        formula = strrep(formula, s{i},num_str);
                    end
                  end
                  
            end
        end
    
            num_subs = size(rxn.substrates,1);
            for j=1:num_subs                     
                id = rxn.substrates(j,1);
                num_str = sprintf('cmpds[%u].value;',id-1);
               if ( def(id,1) == 0 )
                     fprintf(fid,'   int %s;\n',rxn.names{j});
                     def(id,1) = 1;
                end
           
                out_str{is} = sprintf('          %s = %s\n',rxn.names{j}, num_str);
                is=is+1;
            end
                out_str{is} = sprintf('          ai = %s;\n',formula);
                is=is+1;
                out_str{is} = sprintf('\t  break;\n\n');
               is=is+1;
        end
    end
    
    for si=1:(is-1)
       fprintf(fid,'%s',out_str{si});
    end  
    fprintf(fid,'       }    /*   End switch     */\n');  
    fprintf(fid,'   return ai; \n  }\n');  
    fclose(fid);

    movefile(temp_file, fullfile('bns_mat_c','rxn_prob_undef.c'));
     curr_dir = pwd; 
     cd ('bns_mat_c');
        mex_makefile
     cd (curr_dir);
end

% give warning if compound is not used in a reaction, possible mispelling
for j = 1:num_cmpds
    cmpd_used = 0;
    for i = rxn_groupnumbers
        this_cmpd_used = find(rxn_grp{i}.global_comp_gid == j);
        if ~isempty(this_cmpd_used)
            cmpd_used = 1;
            break
        end
    end
    if cmpd_used == 0
        warning(['Compound ' cmpd_info.cmpd_name{j} ' listed in initial compounds lists is not used in any reaction.'])
    end
end
   
   
for i = rxn_groupnumbers
    % create reaction dependency matrix for this reaction group
    rxn_grp{i}.rxn_dependency_matrix = create_dependency_matrix(rxn_grp{i}.rxn, rxn_grp{i}.rxn_dependency_matrix, min_depend{i});
end

if length(rxn_groupnumbers) == 1
    for i = 1:length(plot_info.monitor_cmpds_gid)
        idx_cmpd_grplist = find(plot_info.monitor_cmpds_gid(i) == rxn_grp{rxn_groupnumbers(1)}.global_comp_gid);
        if ~isempty(idx_cmpd_grplist)
            rxn_grp{rxn_groupnumbers(1)}.monitor_comps(i) = idx_cmpd_grplist;
        else
            error(['The compound ' cmpd_info.cmpd_name{plot_info.monitor_cmpds_gid(i)} ' listed as being monitored in storing_and_plotting.m is not used in any reaction'])
        end
    end
else
    rxn_grp{rxn_groupnumbers(1)}.monitor_comps = [];    % doesn't make sense to take snapshots of compound values during simulation of each reaction group if there
                                                        % are multiple reaction groups
end


