function cumSol = accumulateSol(f, dfdx, events, deltaIC, IC, inp, abs_tol, rel_tol, tStop, intDir, nonNeg, tGet)
% function cumSol = accumulateSol(f, dfdx, events, deltaIC, IC, inp, abs_tol, rel_tol, tStop, intDir, nonNeg, tGet)
% 
% This function returns a simulated ODE structure as an output argument 
% cumSol.  cumSol is computed using the following input arguments.
%
% Required input arguments are marked with an asterisk.  If an input
% argument is not required, it can be skipped with an empty bracket.
% 
% *f        -   ODE's right-hand side
%               f takes input arguments (t, x, u) and returns a vector of
%               size nX x 1
% 
% *dfdx     -   Jacobian of ODE
%               dfdx takes input arguments (t, x, u) and returns a matrix
%               of size nX x nX
% 
% events    -   terminating events for the ODE
%               events takes input arguments (t, x, u) and returns the 
%               standard events outputs [val, isTerm, direction].
% 
% deltaIC   -   changes in initial conditions at termination time
%           -   deltaIC takes input arguments (t, x, u) and returns a
%               vector of size nX x 1.
% 
% *IC       -   the initial condition vector with which to initialize
%               simulation.  IC is a vector of size nX x 1
% 
% *inp      -   The input function 
% 
% abs_tol   -   the absolute tolerance for integration.  Defaults to 1x10^-6
% 
% rel_tol   -   the relative tolerance for integration.  Defaults to 1x10^-4
% 
% tStop     -   a vector of times at which to stop integration.
%               Integration will be restarted immediately thereafter.  This
%               can be necessary if inputs change discontinuously, or
%               deltaIC is expected to be evaluated at some discrete time.
% 
% intDir    -   the integration direction.  1 == forward in time
%                                          -1 == backwards in time
%               The default value for intDir is 1.
%
% nonNeg    -   Vector indicating the indexes  that are not allowed to go 
%               negative. Fed directly to the ode solver.
%
% tGet      -   Vector, the timepoints for which the solution is wanted.
%               This does not need to be any particular order
% 
% This work is licensed under the Creative Commons Attribution-Noncommercial-No Derivative 
% Works 3.0 United States License. To view a copy of this license, visit 
% http://creativecommons.org/licenses/by-nc-nd/3.0/us/ or send a letter to Creative Commons,
% 171 Second Street, Suite 300, San Francisco, California, 94105, USA.

% All discontinuities must be equal to their right limit
% All events must be terminal
% Events must only exist in forward integration

%% Work-up
% Clean up inputs
if nargin < 12
    tGet = [];
    if nargin < 11
        nonNeg = [];
        if nargin < 10
            intDir = [];
            if nargin < 9
                tStop = [];
                if nargin < 8
                    rel_tol = [];
                    if nargin < 7
                        abs_tol = [];
                    end
                end
            end
        end
    end
end

% Default inputs
if isempty(abs_tol)
    abs_tol = 1e-6;
end
if isempty(rel_tol)
    rel_tol = 1e-4;
end
if isempty(tStop)
    tStop = [0 inp.tF];
end
if isempty(intDir)
    intDir = 1;
end

% Set flag to return a differentiable solution or just particular points
tGet = vec(tGet)';
if isempty(tGet)
    fullSol = true;
else
    fullSol = false;
end

% tStop will be a column vector sorted ascending
tStop = vec(unique(tStop))';

% Error checking inputs
assert(abs(intDir) == 1, 'KroneckerBio:accumulateSol', 'intDir must be 1 (forward) or -1 (reverse)');

% Constants
forwardBoost            = 10;   % Number of epsilons to jump at a discontinuity
maxFalseEvents          = 1e3;  % Number at which integrator should give up
minTimeBetweenEvents    = 1e-3; % Events are invalid if they occur within this timespan
minSpeciesChange        = 1e-6; % Events are invalid if all species change less than this
N = length(tStop);

%% Initialize variables
% If this input doesn't use events or doesn't care about terminating
% events, allow an infinite number of them.  This means integration will
% continue past terminating events until tF is reached.
if ~isfield(inp, 'nEvents') || isempty(inp.nEvents)
    inp.nEvents = inf;
end

if ~isfield(inp, 'firstEvent') || isempty(inp.firstEvent)
    inp.firstEvent = 0;
end

% ODE options
simOpts = odeset('Jacobian', dfdx, 'AbsTol', abs_tol, 'RelTol', rel_tol, 'NonNegative', nonNeg);
if ~isempty(events)
    simOpts.Events = events;
end

% Initialize event counters
ne = 0;             % number of events past inp.firstEvent
neAll = 0;          % total number of notable events
nFalseEvents = 0;   % total number of false events

% Initialize empty solution
if fullSol
    cumSol = [];
else
    cumSol.solver = 'ode15s';
    cumSol.extdata = [];
    cumSol.x = tGet;
    cumSol.y = zeros(length(IC),length(tGet));
    cumSol.stats = [];
    cumSol.idata = [];
end

% Apply the appropriate deltaIC for the first tStop
if ~isempty(deltaIC)
    if intDir == 1
        IC = IC + intDir*deltaIC(tStop(1), inp.u);
    else
        IC = IC + intDir*deltaIC(tStop(N), inp.u);
    end
end

%% Integrate through each time interval
for k = 1:(N-1)
    % Update the indexes of the time points
    if intDir == 1   	% forward integration
        i0 = k;
        i1 = k+1;
    else                % backward integration
        i0 = N-(k-1);
        i1 = N-k;
    end
    
    % Integration time interval
    if i0 == N || i1 == N
        % At final interval evaluate both points
        tInt = [tStop(i0), tStop(i1)];
    else
        % At intermediate intervals, do not evaluate the forwardmost point
        % Instead, stop a little before going forward or start a little
        % going backward.
        if intDir == 1  % forward
            tInt = [tStop(i0), tStop(i1) - forwardBoost*eps(tStop(i1))];
        else            % backward
            tInt = [tStop(i0) - forwardBoost*eps(tStop(i0)), tStop(i1)];
        end
    end
    
    % Initialize variables for tracking events
    tEventInt = tInt;
    t0 = tInt(1);
    
    % Integrate until the total number of events is equal to inp.nEvents or
    % this time interval is completed
    while ne < inp.nEvents && t0 ~= tEventInt(2) && nFalseEvents < maxFalseEvents
        if fullSol
            % We want the full solution in the interval tEventInt
            tOde = tEventInt;
        else
            % Find tGet points for this interval
            if tEventInt(1) == tStop(end) || tEventInt(2) == tStop(end)
                % Find points matching the end
                itGet = tGet(tGet >= min(tEventInt) & tGet <= max(tEventInt));
            else
                % The end point is covered by another interval
                itGet = tGet(tGet >= min(tEventInt) & tGet < max(tEventInt));
            end
            
            % Construct time vector for ode solver
            tOde = unique([tEventInt, itGet]);
            if numel(tOde) < 3
                % Include a third element to prevent the solver from saving
                % all timepoints
                tOde = [tOde(1), (tOde(1)+tOde(2))/2, tOde(2)];
            end
            if intDir == -1
                % Unique will have sorted tOde the wrong way. Reverse order
                tOde = tOde(end:1);
            end
        end
        
        % Start or restart integration after an event
        if fullSol
            % Ask for single output
            simSol = ode15sf(f, tOde, IC, simOpts, inp.u);
        else
            % Multiple outputs are required
            if isempty(events)
                [simx simy] = ode15sf(f, tOde, IC, simOpts, inp.u);
                simSol.x = simx';
                simSol.y = simy';
            else
                [simx simy simte simye simie] = ode15sf(f, tOde, IC, simOpts, inp.u);
                simSol.x = simx';
                simSol.y = simy';
                simSol.te = simte';
                simSol.ye = simye';
                simSol.ie = simie';
            end
        end
        
        % Check if integration failed when it shouldn't have
        if isempty(events) && simSol.x(end) ~= tInt(2)
            error('KroneckerBio:accumulateSol:IntegrationFailure', 'Did not integrate through entire interval!');
        end
        
        % Discard all but the last event to prevent repeats
        if isfield(simSol, 'ie') && length(simSol.ie) > 1
            simSol.ie = simSol.ie(end);
            simSol.xe = simSol.xe(end);
            simSol.ye = simSol.ye(:,end);
        end
        
        % Combine new solution with cumulative
        if fullSol
            cumSol = mergeSols(cumSol, simSol);
        else
            cumSol = pasteSols(cumSol, simSol);
        end
        
        % Handle any discrete times or event dependent times
        if ~isempty(deltaIC)
            IC = simSol.y(:, end).' + intDir*deltaIC(simSol.x(end), inp.u);
        else
            IC  = simSol.y(:, end).';
        end
        
        % Update time interval that remains to be integrated
        t0         = simSol.x(end);
        tEventInt  = [t0 + forwardBoost*intDir*eps(t0), tInt(2)];
        
        % If we have not reached tF, then a terminating event was
        % detected. Check to ensure that this is not due to numerical
        % problems, but rather that it constitutes a 'true' event.
        %
        % There are at least two criteria that can invalidate an event:
        % 1 -   It is next to the previous event (i.e. a true event that
        %       causes termination twice in succession at very nearby
        %       timepoints)
        % 2 -   It is due to very small change in a state state variable,
        %       a 'wiggle' with low amplitude.
        
        if simSol.x(end) ~= tInt(2)
            % **** IMPORTANT SECTION ****
            % This line makes sure events are not double-counted. This is done
            % by checking if two of the SAME event are within 1e-5 time units
            % of one another.
            
            % Concatenate the initial conditions to the front
            augmentedEventTimes = [0 cumSol.xe];
            augmentedEventAmp   = [cumSol.y(:,1) cumSol.ye];
            
            % Check for a duplicated event
            isDuplicate = length(cumSol.ie) > 1 && ...              % there has been more than one event
                ~diff(cumSol.ie(end-1:end)) && ...                  % AND it's the same type as the previous one
                diff(cumSol.xe(end-1:end) < minTimeBetweenEvents);  % AND it's really near the previous one in time
            
            % Check for a wiggle event
            isWiggle = diff(augmentedEventTimes(end-1:end)) < minTimeBetweenEvents || ...     % this event happened really near the previous one in time
                all(abs(diff(augmentedEventAmp(:,end-1:end), [], 2)) < minSpeciesChange);     % AND no species changed appreciably
            
            if isDuplicate || isWiggle
                % Discard the duplicate event
                cumSol.xe = cumSol.xe(1:neAll);
                cumSol.ie = cumSol.ie(1:neAll);
                cumSol.ye = cumSol.ye(:,1:neAll);
                
                % Record its duplicity
                nFalseEvents = nFalseEvents + 1;
            else
                % Note the event
                neAll = neAll + 1;
                if simSol.x(end) > inp.firstEvent
                    % Note the real event
                    ne = ne + 1;
                end
            end
        end
    end
end

%% Work-down
% Add on any changes to initial conditions at the final timepoint
if ~isempty(deltaIC)
    if fullSol
        cumSol.y(:, end) = cumSol.y(:, end) + intDir*deltaIC(cumSol.x(end), inp.u).';
    else
        % Find the final indexes and add to the solution at those points
        finalInds = find(simSol.x(end) == cumSol.x);
        cumSol.y(:,finalInds) = cumSol.y(:,finalInds) + repmat(intDir*deltaIC(simSol.x(end), inp.u).', 1, length(finalInds));
    end
end

% if no integration was performed return an appropriate structure
if N == 1
    if fullSol
        cumSol.solver = 'ode15s';
        cumSol.extdata = [];
        cumSol.x = tStop;
        cumSol.y = IC;
        cumSol.stats = [];
        cumSol.idata = [];
    else
        simSol.x = tStop;
        cumSol.y(:,tStop == tGet) = repmat(IC, 1,nnz(tStop == tGet));
    end
end

% Make non-events solutions compatible with those that expect them
if isempty(events)
    cumSol.xe = [];
    cumSol.ye = [];
    cumSol.ie = [];
end

% Put NaN in columns that were never reached during integration
if ~fullSol
    if intDir == 1
        cumSol.y(:,cumSol.x > simSol.x(end)) = NaN;
    elseif intDir == -1
        cumSol.y(:,cumSol.x < simSol.x(end)) = NaN;
    end
end

% Order solution structure fields (in case cumSol.xe, .ye and .ie are out of order)
cumSol = orderfields(cumSol, {'solver', 'extdata', 'x', 'y',  'xe', 'ye', 'ie', 'stats', 'idata'});
end