%
% stochastic_sim
%
% chemical equations in create_model:
%
function stochastic_sim

%%%
%%% Definitions of global variables shared with the other files
%%%

% Variables
global num_reactions;
global num_species;
global num_external_species;
global species_names;
global external_species_names;
global max_cells;
global max_alive_cell_index;
global num_alive_cells;
global add_species_after_cell_division_index;
global add_species_after_cell_division_amount;
global num_species_to_add_after_cell_division;
global reaction_reactants;
global reaction_reactants_index;
global reaction_reactants_num;
global reaction_products;
global reaction_rates;
global num_reactions;
global max_reactions;
global max_reactants;
global initial_alive_cells;
global external_reaction_reactants;
global external_reaction_products;
global external_reaction_rates;
global num_external_reactions;
global KILLER_THRESHOLD;
global DIVIDE_THRESHOLD;
global num_external_input_reactions;
global external_input_reactions;
global reaction_has_external_inputs;
global curr_time;
global cell_population;


% Various simulation definitions
global REGRESSION;         % Monte Carlo Simulation Times
global END_TIME;
global RECORD_STEP;
global PLASMID_COPIES;
global CUR_SPECIES_VALUES;
global EXTERNAL_CUR_SPECIES_VALUES;
global REC_SPECIES;
global EXTERNAL_REC_SPECIES;
global AVALS;
global EXTERNAL_AVALS;
global CELL_INTERNAL_AVALS;
global CELL_EXTERNAL_INPUT_AVALS;
global REC_time;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialization & Starting The Simulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

setup_sim_params;  % user defined simulation parameters

%
% other variables
%
species_names          = {};
external_species_names = {};

num_reactions                = 0;
num_external_reactions       = 0;
num_external_input_reactions = 0;
external_input_reactions     = [];
reaction_reactants_index     = zeros(max_reactions,max_reactants);
reaction_reactants_num       = [];

num_species_to_add_after_cell_division = 0;

reaction_reactants = {''};
reaction_products  = {''};
reaction_rates     = {''};
species_names      = {};

external_reaction_reactants = {''};
external_reaction_products  = {''};
external_reaction_rates     = {''};
external_species_names      = {};

%
% call another file to create the model  -- customize if needed
%
create_model_pop_ctrl;

Alive_index  = lookup_species('Alive');
Volume_index = lookup_species('Volume');
Killer_index = lookup_species('Killer');

num_species          = size(species_names,2);
num_external_species = size(external_species_names,2);
AVALS                = zeros(max_cells,num_reactions);
EXTERNAL_AVALS       = zeros(num_external_reactions,1);
REC_POPULATION       = zeros(END_TIME/RECORD_STEP + 1,1);
CELL_INTERNAL_AVALS       = -1 .* ones(max_cells,1);
CELL_EXTERNAL_INPUT_AVALS = -1 .* ones(max_cells,1);

% arrays to store current simulation values
REC_SPECIES = zeros(max_cells,num_species,END_TIME/RECORD_STEP + 1);
EXTERNAL_REC_SPECIES = zeros(num_external_species,END_TIME/RECORD_STEP + 1);

% array to store statistics about multiple simulations
sum_REC_SPECIES = zeros(max_cells,num_species,END_TIME/RECORD_STEP + 1);
EXTERNAL_sum_REC_SPECIES = zeros(max_cells,num_species,END_TIME/RECORD_STEP + 1);

% Setup Time Record Field
index = 1;
REC_time = zeros(size(0 : RECORD_STEP : END_TIME, 2),1);
for i = 0 : RECORD_STEP : END_TIME,
    REC_time(index) = i;
    index = index + 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stochastic Simulation Core Loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1 : REGRESSION
    
    % setup initial species values in the cells and external
    CUR_SPECIES_VALUES = zeros(max_cells,num_species);
    EXTERNAL_CUR_SPECIES_VALUES = zeros(num_species);
    
    for c=1:initial_alive_cells
        CUR_SPECIES_VALUES(c,lookup_species('P_lux'))  = 20;
        CUR_SPECIES_VALUES(c,lookup_species('P_luxR')) = 20;
        CUR_SPECIES_VALUES(c,Alive_index) = 1;
    end
    
    % initialize internal variables
    sum_avals     = 0.0;
    index         = 1;
    curr_time     = 0.0;
    max_alive_cell_index = initial_alive_cells;
    num_alive_cells      = initial_alive_cells;
    
    % display population during simulation
    figure(99);
    clf;
    display_period = 5;
    next_time_to_display = display_period;
    
    %
    % main loop to simulate the system of reactions
    %
    while( curr_time < END_TIME )
        
        % CUR_SPECIES_VALUES
        
        % random numbers used to determine which next reaction and when
        r1 = rand;    % when
        r2 = rand;    % which
        
        if( (r1 ~= 0) && (r2 ~= 0) )
            %
            % compute AVALS for each reaction and sum of AVALS
            %
            sum_avals = 0;
            % first, look at purely external reactions
            for er=1:num_external_reactions
                EXTERNAL_AVALS(er) = external_reaction_rates{er};
                for s=1:size(external_reaction_reactants{er,2},2)
                    s_index = external_reaction_reactants{er,2}(s);
                    if (s_index > 0)
                        error('Internal error: reactant index for external reaction is >0');
                    else
                        rv = EXTERNAL_CUR_SPECIES_VALUES(-s_index);
                    end
                    EXTERNAL_AVALS(er)= EXTERNAL_AVALS(er) * rv;
                end
                sum_avals = sum_avals + EXTERNAL_AVALS(er);
            end
            
            % next, examine all reactions that involve cells
            for c=1:max_alive_cell_index
                if (CELL_INTERNAL_AVALS(c) == -1)
                    % sometimes need to recompute all AVALS for all reactions 
                    % involving this cell
                    CELL_INTERNAL_AVALS(c)       = 0;
                    CELL_EXTERNAL_INPUT_AVALS(c) = 0;
                    for r=1:num_reactions
                        % compute AVAL for cell c, reaction r
                        AVALS(c,r) = reaction_rates{r};
                        for s=1:reaction_reactants_num(r)
                            s_index = reaction_reactants_index(r,s);
                            if (s_index > 0)
                                rv = CUR_SPECIES_VALUES(c,s_index);
                            else
                                rv = EXTERNAL_CUR_SPECIES_VALUES(-s_index);
                            end
                            AVALS(c,r)= AVALS(c,r) * rv;
                        end
                        if (reaction_has_external_inputs(r))
                            CELL_EXTERNAL_INPUT_AVALS(c) = ...
                                CELL_EXTERNAL_INPUT_AVALS(c) + AVALS(c,r);
                        else
                            CELL_INTERNAL_AVALS(c) = ...
                                CELL_INTERNAL_AVALS(c) + AVALS(c,r);
                        end
                    end
                    sum_avals = sum_avals + CELL_INTERNAL_AVALS(c) + ...
                        CELL_EXTERNAL_INPUT_AVALS(c);
                else
                    % other times, only need to recompute external input
                    % reactions for this cell
                    CELL_EXTERNAL_INPUT_AVALS(c) = 0;
                    for r=external_input_reactions
                        % compute AVAL for cell c, reaction r
                        AVALS(c,r) = reaction_rates{r};
                        for s=1:reaction_reactants_num(r)
                            s_index = reaction_reactants_index(r,s);
                            if (s_index > 0)
                                rv = CUR_SPECIES_VALUES(c,s_index);
                            else
                                rv = EXTERNAL_CUR_SPECIES_VALUES(-s_index);
                            end
                            AVALS(c,r)= AVALS(c,r) * rv;
                        end
                        CELL_EXTERNAL_INPUT_AVALS(c) = ...
                            CELL_EXTERNAL_INPUT_AVALS(c) + AVALS(c,r);
                    end
                    sum_avals = sum_avals + CELL_INTERNAL_AVALS(c) + ...
                        CELL_EXTERNAL_INPUT_AVALS(c);
                end     
            end
            
            curr_time = curr_time + log(1/r1) / sum_avals;  % next time
            r2 = r2 * sum_avals;
            record_new_vals = 0;
            
            if( (REC_time(index) < curr_time) && (curr_time < END_TIME) )
                index = index + 1;
                record_new_vals = 1;
                disp(['time ' num2str(curr_time)]);
            end
            
            % disp(['Sum of AVALS = ' num2str(sum_avals) ', r2 = ' num2str(r2)]);
            
            % determine which reaction is next
            tmp_sum = 0.0;
            choose_external_reaction = false;
            % first, consider purely external reactions
            if (num_external_reactions > 0)
                external_reaction_index = 1;
                tmp_sum = tmp_sum + EXTERNAL_AVALS(external_reaction_index);
                while ((tmp_sum < r2) && ...
                        (external_reaction_index < num_external_reactions))
                    external_reaction_index = external_reaction_index + 1;
                    tmp_sum = tmp_sum + EXTERNAL_AVALS(external_reaction_index);
                end
                
                % if chosen, execute external reaction
                if (tmp_sum >= r2)
                    choose_external_reaction = true;
                    for r=1:size(external_reaction_reactants{external_reaction_index,2},2)
                        % each of the reactant species is reduced by one
                        s = external_reaction_reactants{external_reaction_index,2}(r);
                        if (s>0)
                            error('Internal error: external reaction reactant index > 0');
                        else
                            EXTERNAL_CUR_SPECIES_VALUES(-s) = ...
                                EXTERNAL_CUR_SPECIES_VALUES(-s) - 1;
                        end
                    end
                    
                    for r=1:size(external_reaction_products{external_reaction_index,2},2)
                        % each of the product species is increased by one
                        s = external_reaction_products{external_reaction_index,2}(r);
                        if (s>0)
                            error('Internal error: external reaction product index > 0');
                        else
                            EXTERNAL_CUR_SPECIES_VALUES(-s) = ...
                                EXTERNAL_CUR_SPECIES_VALUES(-s) + 1;
                        end
                    end
                end
            end
            
            %
            % if needed, choose and execute internal cell reaction
            %
            %   outer loop: find out which cell had the reaction
            %   inner loop: find out which reaction within the cell
            %
            if (choose_external_reaction == false)
                c = 0;  % cell
                % first, iterate over cells
                while (tmp_sum < r2)
                    c = c + 1;
                    cell_sum_avals = CELL_INTERNAL_AVALS(c) + ...
                        CELL_EXTERNAL_INPUT_AVALS(c);
                    tmp_sum = tmp_sum + cell_sum_avals;
                    if (tmp_sum >= r2)
                        % found cell, now determine reaction
                        tmp_sum = tmp_sum - cell_sum_avals;
                        reaction_index = 1;
                        tmp_sum = tmp_sum + AVALS(c,reaction_index);
                        while (tmp_sum < r2)
                            reaction_index = reaction_index + 1;
                            tmp_sum = tmp_sum + AVALS(c,reaction_index);
                        end
                    end
                end
                
                %tmp_sum = tmp_sum + AVALS(c,reaction_index);
                %while (tmp_sum < r2)
                %    reaction_index = reaction_index + 1;
                %    if (reaction_index > num_reactions)
                %        reaction_index = 1;
                %        c = c + 1;
                %    end
                %    tmp_sum = tmp_sum + AVALS(c,reaction_index);
                %end
                
                % disp(['choose cell ' num2str(c) ' reaction ' num2str(reaction_index)]);
                
                %
                % execute reaction (i.e. update species values)
                %
                for r=1:size(reaction_reactants{reaction_index,2},2)
                    % each of the reactant species is reduced by one in chosen
                    % cell
                    s = reaction_reactants{reaction_index,2}(r);
                    if (s>0)
                        CUR_SPECIES_VALUES(c,s) = CUR_SPECIES_VALUES(c,s) - 1;
                    else
                        EXTERNAL_CUR_SPECIES_VALUES(-s) = ...
                            EXTERNAL_CUR_SPECIES_VALUES(-s) - 1;
                    end
                end
                
                for r=1:size(reaction_products{reaction_index,2},2)
                    % each of the product species is increased by one in chosen
                    % cell
                    s = reaction_products{reaction_index,2}(r);
                    if (s>0)
                        CUR_SPECIES_VALUES(c,s) = CUR_SPECIES_VALUES(c,s) + 1;
                    else
                        EXTERNAL_CUR_SPECIES_VALUES(-s) = ...
                            EXTERNAL_CUR_SPECIES_VALUES(-s) + 1;
                    end
                end
                
                CELL_INTERNAL_AVALS(c) = -1;  % will need to recompute AVALS next time
            end
            
            % record information about species values at this time point
            if (record_new_vals == 1)
                cell_population = 0;
                for cr=1:max_alive_cell_index
                    if (CUR_SPECIES_VALUES(cr,Alive_index) == 1)
                        cell_population = cell_population + 1;
                    end
                    for s=1:num_species
                        REC_SPECIES(cr,s,index) = CUR_SPECIES_VALUES(cr,s);
                    end
                end
                
                if (cell_population ~= num_alive_cells)
                    error('Internal error: cell_population does not equal num_alive_cells');
                end
                
                for e=1:num_external_species
                    EXTERNAL_REC_SPECIES(e,index) = EXTERNAL_CUR_SPECIES_VALUES(e);
                end
                
                REC_POPULATION(index,1) = cell_population;
            end
            
            %
            % if Killer protein greater than threshold, kill cell
            %
            if (CUR_SPECIES_VALUES(c,Killer_index) > KILLER_THRESHOLD)
                for s=1:num_species
                    CUR_SPECIES_VALUES(c,s) = 0;
                end
                CELL_INTERNAL_AVALS(c) = -1;
                num_alive_cells = num_alive_cells - 1;
                disp(['time = ' num2str(curr_time) ', cell ' num2str(c) ' dies']);
                if (num_alive_cells == 0)
                    disp('Population crashed.  Simulation finished.');
                    plot_cell_data(3);
                    return;
                end
            end
            
            %
            % if Volume now exceeds threshold, divide cell
            %
            if (CUR_SPECIES_VALUES(c,Volume_index) > DIVIDE_THRESHOLD)
                next_available_index = 1;
                while (CUR_SPECIES_VALUES(next_available_index,Alive_index) == 1)
                    next_available_index = next_available_index + 1;
                    if (next_available_index > max_cells)
                        disp('Too many cells.  Stopping simulation');
                        return;
                    end
                end
                % distribute species about half/half
                for s=1:num_species
                    CUR_SPECIES_VALUES(next_available_index,s) = ...
                        floor(CUR_SPECIES_VALUES(c,s) / 2);
                    CUR_SPECIES_VALUES(c,s) = CUR_SPECIES_VALUES(c,s) - ...
                        CUR_SPECIES_VALUES(next_available_index,s);
                end
                CUR_SPECIES_VALUES(next_available_index,Alive_index) = 1;
                % add 'magic' species (e.g. for DNA replication)
                for j=1:num_species_to_add_after_cell_division
                    s = add_species_after_cell_division_index(j);
                    num_to_add = add_species_after_cell_division_amount(j);
                    CUR_SPECIES_VALUES(c,s) = CUR_SPECIES_VALUES(c,s) + num_to_add;
                    CUR_SPECIES_VALUES(next_available_index,s) = ...
                        CUR_SPECIES_VALUES(next_available_index,s) + num_to_add;
                end
                CELL_INTERNAL_AVALS(c)                    = -1;
                CELL_INTERNAL_AVALS(next_available_index) = -1;
                max_alive_cell_index = max(max_alive_cell_index, next_available_index);
                num_alive_cells = num_alive_cells + 1;
                disp(['time = ' num2str(curr_time) ', cell ' num2str(c) ...
                    ' divides to occupy position ' num2str(next_available_index)]);
            end
            
            % display population so far and data from two cells
            if (curr_time > next_time_to_display)
                next_time_to_display = next_time_to_display + display_period;
                figure(99);
                hold on;
                plot(REC_POPULATION(1:index));
                title('Cell population');
                plot_cell_data(2);
            end
            
        end
    end
    
    % Summarize the record field at each simulation
    index = 1;
    for j = 0 : RECORD_STEP : END_TIME,
        %     sum_REC_mRNAa(index) = sum_REC_mRNAa(index) + REC_mRNAa(index);
        index = index + 1;
    end
    
    %plot_cell_data(9999);
    
end

end