function [Tn, Eps, N, T] = poisson_process(n, lambda, T, graph, n_instances)
% Simulation of the Poisson process
%
%[TN, N]
%   TN  Times of jumps
%   EPS Sizes of jumps
%   N   Number of jumps
%   T   Time needed for 'n' jumps (if 'T' not given)
%
%POISSON_PROCESS(N,LAMBDA,T,GRAPH,N_INSTANCES)
%   N   Number of jumps to be performed
%       Default value: 10
%       Note: If the 'T' parameter is set, than the jumps are generated
%       until the given time t='T' is reached.
%   LAMBDA  Parameter for exponencial distribution
%           Default value: 1
%
%   T   Time interval [0,T]
%       Default value: 0, that means that the 'n' jumps are performed
%   GRAPH   'true' or 'false' whether to plot the result or not
%           Default value: false
%   N_INSTANCES Number of instances of the process to be plotted
%               Default value: 1
%               Note: If the number of instances to plot is bigger than 2,
%               then the returned vector W is highlighted with RED colored
%               DASH-DOT line

%% Input parameters handling
if nargin < 5; n_instances = 1; end
if nargin < 4; graph = false; end
if nargin < 3; T = 0; end
if nargin < 2; lambda = 1; end
if nargin < 1; n = 10; end

def_T = T;  % Default 'T', needed for the recursion in plotting
%TODO: general parameter's check

%% Computation task
% If 'T==0', than no concrete time interval is given,
% Generate only 'n' exponential variables.
% Else, generate points until the sum of epsilons is less than given 'T'.
if T == 0  
    Eps = exprnd(lambda, [n, 1]);
    T = ceil(sum(Eps));  % If 'T' not set, set it accordingly to the jumps
else
    Eps = exprnd(lambda);
    if Eps > T
        fprintf("No jumps occured to time t=" + T + "\n")
        N = 0; Tn = 0; return
    end
    while 1
        e = exprnd(lambda);
        if sum(Eps) + e < T
            Eps = cat(1, Eps, e);
        else
            break
        end
    end
end

Tn = [0; cumsum(Eps)];  % Times of jumps
N = Nt(Tn, T);


%% Plotting
if graph
    figure
    maxN = N; maxT = T; % Needed for axes
    t = [Tn; T];
    jump_values = (0:N);
    stairs2(t, jump_values, 'left', '-.r'); hold on;  % Plot the first one
    if n_instances > 1  % If the number of instances is bigger, plot others
        for i = 2:n_instances
            [Tn_new, ~, N_new, T_new] = poisson_process(n, lambda, def_T);
            maxN(maxN < N_new) = N_new;
            maxT(maxT < T_new) = T_new;
            t = [Tn_new; T_new];
            jump_values = (0:N_new);
            stairs2(t, jump_values, 'left'); hold on;
        end
    end
    % Plot editing
    a = gca;  % Get the current axis
    a.Position([1 2 4]) = [.11 .3 .6];  % [1, 2] position, [3, 4] w, h
    title("Simulation of the Poisson process")
    xlabel("\it t")
    ylabel("\it N_t")
    xlim([0 maxT]);
    ylim([0 maxN]);
    %yticks(0:1:maxN+1);
    grid on;
    message = "\lambda: " + lambda + newline + ...
              "Number of runs: " + n_instances + newline + ...
              "Max T: " + maxT + newline + ...
              "Max N: " + maxN;
    annotation('textbox',[.1 .2 0 0], ...
               'String',message, ...
               'FitBoxToText','on', ...
               'BackgroundColor', 'white')
end

end

%% Computes the number of jumps to time 't'
function N = Nt(Tn, t)
    N = binary_search(Tn, t) - 1;  % Tn(N) <= t < Tn(N+1)
end
