tic

clear_stuff
other_constants
collect_rundata

if isempty(randseed)
    temclock = sum(fix(clock));
    randseed = temclock;
end
run_info.randseed = randseed; 

partname = [output_filename '*'];
delete(fullfile(output_dir, partname ));

exec_file = fullfile(output_dir, [output_filename '_timesexecuted']);
fid_exec = fopen(exec_file,'wt');

sbml_check

if use_sbml ~= 1
        read_cmpds
        read_reactions
        read_reaction_constants
        prepare_rxn_grps 
        
        StoP = [];
else
     curr_dir = pwd;
     cd ('bns_mat_c');
        mex_makefile
     cd (curr_dir);
    % Save initial kinetic constants
              for i_rxn = 1: num_rxns
               rxn = rxn_grp{1}.rxn{i_rxn};
               StoP{i_rxn} = rxn.StoP;
              end
end
 
num_cmpds = cmpd_info.num_cmpds;
orig_cmpd_vals = cmpd_info.orig_cmpd_vals;
num_grps = length(rxn_grp);

run_id = 0;

    for i_setnum = 1:numsets
      if use_paramset
         thisset = pset(i_setnum);
         [cmpd_info, rxn_grp] = update_parameters(thisset, orig_cmpd_vals, num_cmpds, cmpd_info, scaling_factor, use_sbml, num_rxns, rxn_grp, StoP);
        
        if isfield(thisset,'newrateparms')
            if use_sbml ~= 1  
            % Not sbml
                read_reactions
                read_reaction_constants
                prepare_rxn_grps   
            end
        end % newrateparms
       
        if isfield(thisset,'numruns')
            num_repeats = thisset.numruns;
        end
      end
    
        init_plots
        plot_info.mpi_save_param = 0;
        
        rxn_grp{1}.cmpd_nums = create_cmpd_matrix(rxn_grp{1}.global_comp_gid, cmpd_info);
        [modrxn, moddepend, mod_monitor] = prepare_matlab_indices_for_c(rxn_grp{1}.rxn, rxn_grp{1}.rxn_dependency_matrix, rxn_grp{1}.monitor_comps);
        write_bns_textfile(modrxn,rxn_grp{1}.cmpd_nums(:,1), mod_monitor, run_info, plot_info, dir_info, numsets);
        movefile([output_filename '_constants.txt'], output_dir);
        
     save_rundata(output_dir, output_filename, plot_info, 'plot_info', run_info, 'run_info', dir_info, 'dir_info', cmpd_info, 'cmpd_info', rxn_grp, 'rxn_grp');

        for i_repeat_num = 1:num_repeats

            run_id = run_id + 1;

            if run_id == 1  % only seed the random number generator for first simulation
                run_info.use_randseed = 1;
            else
                run_info.use_randseed = 0;
            end

            rxn_grp = runsim(cmpd_info, plot_info, run_info, dir_info, rxn_grp, run_id);
  
    % save to disk the number of times each reaction was executed 
       for i = 1:num_grps
         save_timesexecuted(fid_exec, i, i_setnum, run_id, rxn_grp{i}.times_rxns_executed);
       end

        end
    end


if plot_interval ~= -2
    bns_plot(output_dir, output_filename);
end

fclose(fid_exec);

real_time_elapsed = toc
