function varargout = bns_mpi(varargin)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Use MPI interface - version 2
%% Master processor sends packet of jobs to each worker and
%% perform ssimulations itself. This approach reduces communication
%% and allows each processor to write event_log file
%%                         Slava Ch.     02/15/06
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


    tic
% Initialize MPI
MPI_Init;

global MPI_COMM_WORLD
% Create communicator
comm = MPI_COMM_WORLD;

% Get size and rank
comm_size = MPI_Comm_size(comm);
my_rank = MPI_Comm_rank(comm);

%Set the master node
master = 0;

% Create base message tegs
 coefs_tag = 10000;
 input_tag = 20000;
 output_tag = 30000;
% Master's work
if (my_rank == master)
    
  code_dir= pwd;

clear_stuff
default_directories
time_info.randseed = randseed;

rxns_dir = fullfile(models_dir,specific_model_dir);
param_dir = fullfile(models_dir,specific_model_dir,parameters_dir);

if exist(models_dir) == 0
    disp([models_dir '  was not found.  Fix models_dir entry in default_directories.m.']);
    return
end
if exist(rxns_dir) == 0
    disp([rxns_dir '  was not found.  Fix rxns_dir entry in default_directories.m.']);
    return
end
if exist(param_dir) == 0
    disp([param_dir '  was not found.  Fix param_dir entry in default_directories.m.']);
    return
end
if ~isempty(output_directory) & exist(output_directory) == 0
    disp([output_directory '  was not found.  Fix entry in default_directories.m.']);
    return
end

if isempty(output_directory)
    output_directory = 'output';
%    output_dir = fullfile(models_dir,specific_model_dir,parameters_dir,plot_info.output_directory);
    output_dir = fullfile(param_dir,output_directory);
else
    output_dir = output_directory;
end

output_file = fullfile(output_dir, output_filename);

 
     other_constants

%    Modified to work with SBML files    -    Slava Chushak
% Look for *.mat files in the rxns_dir
    file_list = dir(param_dir);
% If exist *.mat file, check whether it is an SBML model fiel
    use_sbml = 0;
    clear sbml_model
    for i = 1:size(file_list,1)
        file_name = file_list(i).name;
        if ~isempty(strfind(file_name,'.xml'))
            file_name = fullfile(param_dir,file_name);
            sbml_model = TranslateSBML(file_name);
            if isfield(sbml_model,'typecode')
                use_sbml = 1;
                break;
            end
        end
    end
    if use_sbml ==1
        [rxn_grp cmpd_info] = read_sbml_model(sbml_model);
        for i = 1:length(plot_info.monitor_cmpds)
            found_gid = 0;
         for j = 1:cmpd_info.num_cmpds
             if strcmp(plot_info.monitor_cmpds{i},cmpd_info.cmpd_name{j})
                  plot_info.monitor_cmpds_gid(i) = j;
                  found_gid = 1;
                  break;
             end
         end
         % Additional check for SBML Level2
         if (found_gid == 0) && (sbml_model.SBML_level == 2)
             for j = 1:cmpd_info.num_cmpds
                if strcmp(plot_info.monitor_cmpds{i},sbml_model.species(1,j).name)
                      plot_info.monitor_cmpds_gid(i) = j;
                      found_gid = 1;
                      break;
                end
             end
         end
         if found_gid == 0
             error([ 'Cannot find gid for plot_info.monitor_cmpds - '  plot_info.monitor_cmpds{i}]);
         end
  end

    else

     
    read_cmpds
    read_reactions
    read_reaction_constants
    prepare_rxn_grps

end
      mex gateway_erg.c rxn_prob_undef.c
      init_plots
      plot_info.event_log_flag = event_log_plot_flag;

    save(output_file, 'randseed', 'plot_info', 'time_info', 'cmpd_info')


% Create data to be broadcasted
     data.cmpd_info = cmpd_info;
     data.plot_info = plot_info;
     data.time_info = time_info;
     data.rxn_grp = rxn_grp;
     data.volume = volume;
     data.stoch_algor = stoch_algor;
     
     MPI_Bcast(master, coefs_tag, comm, data);
      
     N_runs = time_info.num_runs;
     
  %    MAIN LOOP
   % Master processor sends number of runs for each of
   % the processors and also performs runs:

   % Calculate number of runs for each of the processors
 
   if (comm_size < N_runs)

      n_runs_i = floor(N_runs/comm_size);
        add_i  = mod(N_runs, comm_size);

        if add_i ~= 0
  % Send n+1 number of runs add_i workers
            n_runs_i1 = n_runs_i + 1;
            for i=1:add_i
              dest = mod((i - 1), (comm_size - 1)) + 1;
              dest_tag = input_tag + i;

              MPI_Send(dest, dest_tag, comm, n_runs_i1);
            end
        end
        n_procs = comm_size;
   else
        n_procs = N_runs;
        n_runs_1 = 1;
        add_i = 0;
   end

 % Send n_runs_i number of runs to rest of the workers

           for i= add_i+1 : n_procs
              dest = mod((i - 1), (comm_size - 1)) + 1;
              dest_tag = input_tag + i;

              MPI_Send(dest, dest_tag, comm, n_runs_i);
            end

    % End first part of master's work
    % Since workers are busy with simulations,
    % run simulations on master processor

       for i=1:n_runs_i

            % Do computations

%        [rxn_grp, data_time{i}, data_cmpds{i}] = runsim(cmpd_info, plot_info, time_info, rxn_grp, volume, stoch_algor);
         [rxn_grp, data_time{i}, data_cmpds{i}, event_log] = runsim(cmpd_info, plot_info, time_info, rxn_grp, volume, stoch_algor);

        if event_log_plot_flag == 1
            save_event_log(i, event_log, output_file);
        end

       end

    save(output_file, '-append', 'data_time', 'data_cmpds')

    % Receive results
        n_received = n_runs_i; % Results from master's simulations

        while n_received < N_runs

           [message_ranks, message_tags] = MPI_Probe( '*', '*', comm);
      % If message_ranks is not empty - receive the output
             if ~isempty(message_ranks)

      % Sort the received messages
                [m_tags_sorted, m_tags_sorted_idx] = sort(message_tags);

                dest = message_ranks(m_tags_sorted_idx(1));
                master_tag = message_tags(m_tags_sorted_idx(1));
                mesg_num = master_tag - output_tag;

                output = MPI_Recv( dest, master_tag, comm);

                i_runs = output.i_runs;

             for i=1:i_runs
                data_time{n_received + i} = output.data_time{i};
                data_cmpds{n_received + i} = output.data_cmpds{i};
             end

                n_received = n_received + i_runs;

             end
        end

          save(output_file, '-append', 'data_time', 'data_cmpds')

 else

   % Workers
   % Receive input data
     data = MPI_Recv( master, coefs_tag, comm);
            cmpd_info = data.cmpd_info;
            plot_info = data.plot_info;
            time_info = data.time_info;
            rxn_grp = data.rxn_grp;
            volume = data.volume;
            stoch_algor = data.stoch_algor;

   % Receive the input data and run simulations
   for i=1:(comm_size-1)
  % Compute who the destination is for this message.
   dest = mod((i - 1),(comm_size - 1)) + 1;

  % Check if this destination is me.
     if (my_rank == dest)
    % Compute tags.
       dest_tag = input_tag + i;
       master_tag = output_tag + i;

    % Receive input.
       i_runs =  MPI_Recv(master,dest_tag,comm);

       for j=1:i_runs

            % Do computations

%        [rxn_grp, data_time{j}, data_cmpds{j}] = runsim(cmpd_info, plot_info, time_info, rxn_grp, volume, stoch_algor);
         [rxn_grp, data_time{j}, data_cmpds{j}, event_log] = runsim(cmpd_info, plot_info, time_info, rxn_grp, volume, stoch_algor);
 
       end
            % Send results back to the master
              out.data_time = data_time;
              out.data_cmpds = data_cmpds;
              out.i_runs = i_runs; 
            MPI_Send(master, master_tag, comm, out);

     end

   end
end

MPI_Finalize;
disp(' SUCCESS');
      real_time_elapsed = toc

