% Author: Elsa W. Birch, ewbirch@stanford.edu, Covert Lab - Stanford Bioengineering
% Date: 4/24/2011
% Version: associated with production of the figures for Birch et al. 2012
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Summary: This function initializes a dynamic rFBA model by running it with inexhastable
% media for longer than the timescale of regulatory changes to approximate pre-culture of
% the organism in log phase growth in a dilute culture.  
% 
%
% Inputs:   organism - this is the organism struture as produced by
%                      FBA_EB_organismFileParse.m, and with fields .model, .simulation, etc,
%                      and may also have the field .modelInfected (see accompanying
%                      documentation.
%           media - this is a column vector of media associated with the simulation as a
%                   whole, *not* already translated into particular organism format.
%           tStep - length of time step used in model, currently this is not used, but
%                   rather a timestep is determined within the function (because the
%                   timestep required by the overall simulation may be much smaller - and
%                   therefore time consuming - as required by ODE portions of model that
%                   are not used in initalization.
%           stimuliCurrent - cell array of stimuli valid for initialization state, assumed
%                            to be in order/format of *this* organism, not whole
%                            simulation.
%
%
% Outputs:  organismSimulation - outputs the simulation structure of organism with first
%                                column/element filled with the relevent value from the
%                                pre-culture simulation.  in the case that the organism
%                                being intialized is infected the uninfected model is used
%                                for all initialization procedures and then the output
%                                simulation is of the infected demensions.
%           initializationWarnings - cell array of warning message strings.
%
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
%
% Notes: 
%
% Currently it tests two initialization cases of regulation, either all rules on or all
% off, and outputs (as argument) a warning if these simulations do not result in idential
% end timepoints.  Ideally this function would be modified both to intake a 'defualt'
% state of regulation for each regulatory element, and/or test formally for convergance to
% a steady state solution.
%
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

% Main function
function [organismSimulation, initializationWarnings] = ...
                    FBA_EB_intitializerPreculture(organism, media, tStep, stimuliCurrent)

% time parameters for testing steady state convergence of fluxes (assume 5x(max delay for gene activity)
maxTimeDelay = max( organism.model.activityDelayList );
tStep = maxTimeDelay / 5; % this is a reset of the timestep
nTestSteps = 7 * maxTimeDelay / tStep;
organismInput = organism; % save unmodified organism

% call subfunction that sets organism.simualtion structure equal to that of uninfected model run for
% required initialization time
organism = setUpInitializationSimStruct( organism, nTestSteps );

% initialize warnings
initializationWarnings = {};

% First calculate growth optimal fluxes without regulation
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

% steady State Right Hand Side
rhsSpecifications = {[], [], [] ,[] ,[]}; % no set (not zero), free, upper, lower, or double bounds set

initialBiomass = organism.simulation.initialBiomass;

mediaFluxBounds = FBA_EB_setCurrentExchangeFluxBounds(tStep, initialBiomass, media(organism.mediaIdxTranslation),...
                        organism.model.fluxBoundsDefault, organism.model.mediaExchangeRxnIdxList);

% biomass max objective function
objectiveCoefficientVector = zeros( organism.model.nRxns, 1 );
objectiveCoefficientVector(  organism.model.cycleIdxsUser ) = ones(size(organism.model.cycleIdxsUser) )*(-1e-5);
objectiveCoefficientVector(organism.model.rxns.aBIOMASS.idxRxn) = 1;

% simulate using linear programming optimiation (currently glpk, wtih 0 RHS)
[mediaOnlyFlux, status] = ...
    FBA_EB_glpkSolve( organism.model.stoichMatrix,...
                mediaFluxBounds(:,1),...
                mediaFluxBounds(:,2),...
                objectiveCoefficientVector,...
                rhsSpecifications);

% test for solver failure
if ~( (status == 2) || (status == 5) || ~any(mediaOnlyFlux))
    initializationWarnings = [initializationWarnings; 'Initialization using only media bounds failure/questionable/zero.']
end % if solver not optimal

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 


% Calculate growth optimal fluxes with regulation all on/off
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

% store separate organism structures for alternate initialization state tests
organismTestOn = organism;
organismTestOff = organism;

organismTestOn = firstTimepointSet( organismTestOn, 1, mediaOnlyFlux, mediaFluxBounds );
organismTestOff = firstTimepointSet( organismTestOff, 0, mediaOnlyFlux, mediaFluxBounds );

iSim = 1;
% loop through initialization simulation iterations
while iSim < nTestSteps
    
    % update counter
    iSim = iSim + 1;
    
    % biomass update (no growth ~ inf media)
    organismTestOn.simulation.biomass(iSim) = organismTestOn.simulation.biomass(iSim - 1);
    organismTestOff.simulation.biomass(iSim) = organismTestOff.simulation.biomass(iSim - 1);
    
    % set reg/ media constraints
    fluxConstraintFunction = organism.methods.fluxConstraintHandle;
    [organismTestOn.simulation] = fluxConstraintFunction(tStep, iSim, organismTestOn,...
                                                            media(organism.mediaIdxTranslation),stimuliCurrent);
    [organismTestOff.simulation] = fluxConstraintFunction(tStep, iSim, organismTestOff,...
                                                            media(organism.mediaIdxTranslation),stimuliCurrent);
    
    % simulate using linear programming optimiation (currently glpk, wtih 0 RHS)
    [organismTestOn.simulation.fluxes(:,iSim), organismTestOn.simulation.status(iSim)] = ...
    FBA_EB_glpkSolve( organismTestOn.model.stoichMatrix,...
                organismTestOn.simulation.fluxBoundsLower(:,iSim),...
                organismTestOn.simulation.fluxBoundsUpper(:,iSim),...
                objectiveCoefficientVector,...
                rhsSpecifications);
    [organismTestOff.simulation.fluxes(:,iSim), organismTestOf.simulation.status(iSim)] = ...
    FBA_EB_glpkSolve( organismTestOff.model.stoichMatrix,...
                organismTestOff.simulation.fluxBoundsLower(:,iSim),...
                organismTestOff.simulation.fluxBoundsUpper(:,iSim),...
                objectiveCoefficientVector,...
                rhsSpecifications);
    
	% clean up cycles
	organismTestOff.simulation.fluxes(:,iSim) = FBA_EB_cycleRemove(organismTestOff, iSim);
	organismTestOn.simulation.fluxes(:,iSim) = FBA_EB_cycleRemove(organismTestOn, iSim);
    
end % for while iSim

% compare two final states, within threshold
fluxDistEuclid = sum((organismTestOff.simulation.fluxes(:,iSim) - organismTestOn.simulation.fluxes(:,iSim)).^2);
if fluxDistEuclid > 1e-4 % arbitrary thresh
    initializationWarnings = [ initializationWarnings; 'Check initilization'];
end % if difference between fulx distribution

% 
% figure
% plot(organismTestOff.simulation.fluxes(organism.model.rxns.aBIOMASS.idxRxn,:), 'x')
% hold on
% plot(organismTestOn.simulation.fluxes(organism.model.rxns.aBIOMASS.idxRxn,:), 'o')

% set initialized organism
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

% restore simulation structure to that input to function
organism = organismInput;

% check if initializing infected or uninfected model

if organism.infectionFlag
    
    % case when organism being intitialized  is not infeted
    %   flux bounds    
    organism.simulation.fluxBoundsLower(1:organism.modelInfected.nHostRxns,1) = organismTestOn.simulation.fluxBoundsLower(:,iSim);
    organism.simulation.fluxBoundsUpper(1:organism.modelInfected.nHostRxns,1) = organismTestOn.simulation.fluxBoundsUpper(:,iSim);

    %   fluxes
    organism.simulation.fluxes(1:organism.modelInfected.nHostRxns,1) = organismTestOn.simulation.fluxes(:,iSim);
    %   biomass, growth, status
    organism.simulation.biomass(1) = organism.simulation.initialBiomass; 
    organism.simulation.growth(1) = organismTestOn.simulation.fluxes(organism.model.rxns.aBIOMASS.idxRxn, iSim);
    organism.simulation.status(1) = organismTestOn.simulation.status(iSim);
    %   genes, geneProducts, TFs,
    organism.simulation.genes(:,1) = organismTestOn.simulation.genes(:,iSim);
    organism.simulation.geneProducts(:,1) = organismTestOn.simulation.geneProducts(:,iSim);
    organism.simulation.tfs(:,1) = organismTestOn.simulation.tfs(:,iSim);
    organism.simulation.rxnRegulation(1:organism.modelInfected.nHostRxns,1) = organismTestOn.simulation.rxnRegulation(:,iSim);
else
    
    % case when organism being intitialized  is not infeted
    %   flux bounds    
    organism.simulation.fluxBoundsLower(:,1) = organismTestOn.simulation.fluxBoundsLower(:,iSim);
    organism.simulation.fluxBoundsUpper(:,1) = organismTestOn.simulation.fluxBoundsUpper(:,iSim);

    %   fluxes
    organism.simulation.fluxes(:,1) = organismTestOn.simulation.fluxes(:,iSim);
    %   biomass, growth, status
    organism.simulation.biomass(1) = organism.simulation.initialBiomass; 
    organism.simulation.growth(1) = organismTestOn.simulation.fluxes(organism.model.rxns.aBIOMASS.idxRxn, iSim);
    organism.simulation.status(1) = organismTestOn.simulation.status(iSim);
    %   genes, geneProducts, TFs,
    organism.simulation.genes(:,1) = organismTestOn.simulation.genes(:,iSim);
    organism.simulation.geneProducts(:,1) = organismTestOn.simulation.geneProducts(:,iSim);
    organism.simulation.tfs(:,1) = organismTestOn.simulation.tfs(:,iSim);
    organism.simulation.rxnRegulation(:,1) = organismTestOn.simulation.rxnRegulation(:,iSim);
    
end % if infectionFlag



organismSimulation = organism.simulation;

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

end % function FBA_EB_rFBAintitializerPreculture

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

% Subfunctions

function initializedOrganism = firstTimepointSet( organism, regulationState, fluxState, fluxBounds )

% set up simuation storage for this organism
    %   flux bounds    
    organism.simulation.fluxBoundsLower(:,1) = fluxBounds(:,1);
    organism.simulation.fluxBoundsUpper(:,1) = fluxBounds(:,1);
    
    %   fluxes
    organism.simulation.fluxes(:,1) = fluxState;
    %   biomass, growth, status
    organism.simulation.biomass(1) = organism.simulation.initialBiomass; 
    organism.simulation.growth(1) = fluxState(organism.model.rxns.aBIOMASS.idxRxn);
    organism.simulation.status(1) = 0;
    %   genes, geneProducts, TFs, stimuli
    if regulationState
        organism.simulation.genes(:,1) = true(organism.model.nGene, 1);
        organism.simulation.geneProducts = true(organism.model.nGene, 1);
        organism.simulation.tfs = true(organism.model.nTFs, 1);
%         organism.simulation.stimuli = true(organism.model.nStimuli, 1);
    else
        organism.simulation.genes(:,1) = false(organism.model.nGene, 1);
        organism.simulation.geneProducts(:,1) = false(organism.model.nGene, 1);
        organism.simulation.tfs(:,1) = false(organism.model.nTFs, 1);
%         organism.simulation.stimuli(:,1) = false(organism.model.nStimuli, 1);
    end % if regulationState

    initializedOrganism = organism;
    
end

function organism = setUpInitializationSimStruct( organism, nTimeStepsInit )
    % set up simuation storage for this organism
        %   flux bounds
        organism.simulation.fluxBoundsLower = ...
            zeros(organism.model.nRxns, nTimeStepsInit );
        organism.simulation.fluxBoundsUpper = ...
            zeros(organism.model.nRxns, nTimeStepsInit );
        %   fluxes
        organism.simulation.fluxes = ...
            zeros(organism.model.nRxns, nTimeStepsInit );
        %   biomass, growth, status
        organism.simulation.biomass = zeros(nTimeStepsInit, 1); 
        organism.simulation.growth = zeros(nTimeStepsInit, 1);
        organism.simulation.status = zeros(nTimeStepsInit, 1);
        %   genes, geneProducts, TFs, stimuli
        organism.simulation.genes = ...
            false(organism.model.nGene, nTimeStepsInit);
        organism.simulation.geneProducts = ...
            false(organism.model.nGene, nTimeStepsInit);
        organism.simulation.tfs = ...
            false(organism.model.nTFs, nTimeStepsInit);
        organism.simulation.rxnRegulation = ...
            false(organism.model.nRxns, nTimeStepsInit);
        
        % set stimuli at all timepoints equal to default
        organism.simulation.stimuli = ...
             repmat(organism.model.stimuliValuesDefault, 1, nTimeStepsInit);
        
        % initial biomass
        organism.simulation.biomass(1) = organism.simulation.initialBiomass;

end
