function [rxn_grp, cmpd_info] = read_sbml_model(model)
%  Read *.mat file that contain sbml-type model 
%  File is produced by SBMLToolkit 
%
%           Created  by Slava Chushak         
%                     version 1.0 - 12/06/2005

    
clear cmpd_info
clear function rxn_prob_undef

cmpd_info.cmpd_name = [];
cmpd_info.orig_cmpd_vals = [];
 
  rxn_grp{1}.comps = [];
  rxn_grp{1}.global_comp_gid = [];
    tmp_depend = [];
 
%  Prepare cmpd_info structure

     num_cmpds =  size(model.species,2);
     cmpd_info.num_cmpds = num_cmpds;
  
     
     if (model.SBML_level == 2)
         name_str = 'model.species(i).id';
     else
         name_str = 'model.species(i).name';
     end
     
     for i = 1:num_cmpds
         cmpd_info.cmpd_name{i} = eval(name_str);
         cmpd_info.orig_cmpd_vals(i,1) = model.species(i).initialAmount;
         cmpd_info.orig_cmpd_vals(i,2) = 0;
     end
     
  % Create list of parameters   
  n_param = size(model.parameter,2);
      if (model.SBML_level == 2)
         name_str = 'model.parameter(i).id';
     else
         name_str = 'model.parameter(i).name';
      end
   
  if n_param > 0   
      for i = 1:n_param
          parameters{i} = eval(name_str);
      end
  end
  
%  Prepare reactions
    num_rxn = size(model.reaction,2);
    
 % Create a c-file with Kinetic Law for "undef" reactions
 %                         Slava Chushak    01/20/2006


    temp_file = 'undef_prob';
    fid = fopen(temp_file, 'w');
    if fid == -1
            error('Cannot open temporary undef_prob file')
            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 and rates are assigned in prepare_undef_reactions\n');
    fprintf(fid,'**   Need to assign the number of compounds   \n*/\n');
    fprintf(fid,'     double ai=0;\n');

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

    for k = 1:num_rxn
        clear rxn
        
      num_subs = size(model.reaction(k).reactant,2);
      num_prods = size(model.reaction(k).product,2);
      out_str{is} = sprintf('       case %u:\n',k-1);
      is=is+1;

              for i = 1:num_subs
                  subs_name = model.reaction(k).reactant(i).species;
                  rxn_grp{1} = build_comps(subs_name, rxn_grp{1}, cmpd_info.cmpd_name);
                  comp_id = find_comp_id(subs_name,rxn_grp{1}.comps);
                  rxn.substrates(i,1) = comp_id;
                  cmpds_id = find_comp_id(subs_name,cmpd_info.cmpd_name);
                  if (model.species(cmpds_id).boundaryCondition ==1)
                      rxn.substrates(i,2) = 0;
                  else
                      rxn.substrates(i,2) = model.reaction(k).reactant(i).stoichiometry;
                  end
                  names{i} = subs_name;
              end

                for i = 1:num_prods
                    prod_name = model.reaction(k).product(i).species;
                    rxn_grp{1} = build_comps(prod_name, rxn_grp{1}, cmpd_info.cmpd_name);
                    comp_id = find_comp_id(prod_name,rxn_grp{1}.comps);
                    rxn.products(i,1) = comp_id;
                     cmpds_id = find_comp_id(prod_name,cmpd_info.cmpd_name);
                     if (model.species(cmpds_id).boundaryCondition ==1)
                         rxn.products(i,2) = 0;
                     else
                         rxn.products(i,2) = model.reaction(k).product(i).stoichiometry;
                     end
                end

           rxn.type = 'undef';
           rxn.id = k;
           rxn.name = ['rxn' num2str(k)];
% Save kinetic constants in case we need to modify them
           n_kin_const = size(model.reaction(1,k).kineticLaw.parameter,2);
           if n_kin_const > 0
               for j = 1:n_kin_const
                    rxn.StoP_name{j} = model.reaction(1,k).kineticLaw.parameter(j).name;
                    rxn.StoP(j) = model.reaction(1,k).kineticLaw.parameter(j).value;
               end
           end
           rxn_grp{1}.rxn_dependency_matrix{k} = [];
           min_depend{1}(k) = 0;

% Write out Kinetic Law formula into the file

           formula = model.reaction(k).kineticLaw.formula;
           rxn.formula = formula;
 %  Translate a Kinetic Law formula       
     
 %  1. Extract all variables
     
      s = symvar(formula);

 %  Assign values for parameters, constants 
        if size(s) > 0          
            for i = 1:size(s)
                temp_str = [];
                    found_parameter = 0;
                    if n_param >0
                       for j = 1:n_param  %  Check whether it is a parameter
                            if strcmp(s{i},parameters{j})
                                                 % Replace parameter with it's value
                              num_str =num2str(model.parameter(j).value); 
                              formula = strrep(formula, s{i},num_str);
                              rxn.formula = strrep(rxn.formula, s{i},num_str);
                              found_parameter = 1;
                              break;
                            end
                       end
                    end
               if found_parameter == 0               
                        found_rate = 0;
                        % Check whether it is a rate
                  if n_kin_const > 0          
                        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_rate = 1;
                              break;
                            end
                        end
                  end
                        if ( found_rate == 0) && (model.SBML_level == 2)
                            % For Level 2 - check Ids
                            for j = 1:size(model.reaction(1,k).kineticLaw.parameter,2)
                                if strcmp(s{i},model.reaction(1,k).kineticLaw.parameter(j).id)
                                                              % Replace rate with it's value
                                num_str =num2str(model.reaction(1,k).kineticLaw.parameter(j).value); 
                                formula = strrep(formula, s{i},num_str);
                                found_rate = 1;
                                break;
                                end
                            end
                        end
                        if found_rate == 0
                            found_subs = 0;
                            % Check whether it is a substrate
                              for j=1:num_subs
                                if strcmp(s{i},names{j}) 
                                   id = rxn.substrates(j,1);                                
                                   num_str = sprintf('cmpds[%u].value',id-1);
                                   if ( def(id,1) == 0 )
                                       fprintf(fid,'   int %s;\n',s{i});
                                       def(id,1) = 1;
                                   end  
                                  out_str{is} = sprintf('       %s = %s;\n',names{j}, num_str);
                                  is=is+1;
                                   found_subs = 1;
                                break;
                                end
                              end
                           if found_subs == 0
                            % Final check in the list of all substrates
                                   rxn_grp{1} = build_comps(s{i}, rxn_grp{1}, cmpd_info.cmpd_name);
                                   id = find_comp_id(s{i},rxn_grp{1}.comps);
                                   num_str = sprintf('cmpds[%u].value;',id-1);
                                  if ( def(id,1) == 0 )
                                       fprintf(fid,'   int %s;\n',s{i});
                                       def(id,1) = 1;
                                  end
                                   out_str{is} = sprintf('          %s = %s\n',s{i}, num_str);
                                   is=is+1;
                                   found_subs = 1;  
                            end

                            if found_subs == 0
                               fclose(fid);
                               error(['Cannot found constant ' s{i} ' in the reaction ' model.reaction(1,k).name ]);
                            end
                        end
               end
                            
            end
            out_str{is} = sprintf('       ai = %s;\n',formula);
            is=is+1;
            out_str{is} = sprintf('       \n\tbreak;\n');
            is=is+1;
        end
        rxn_grp{1}.rxn{k} = rxn;
    end
    
 %       fprintf(fid,'     mexPrintf("   rxn_id= %%d;\\n", rxn_id);\n');
    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'));
    
        % create reaction dependency matrix for this reaction group
   rxn_grp{1}.rxn_dependency_matrix = create_dependency_matrix(rxn_grp{1}.rxn, rxn_grp{1}.rxn_dependency_matrix, min_depend{1});



    
