function [W, t] = wiener_process(n, T, graph, n_instances, tplot)
% Simulation of the Wiener process
% 
%[W, t]
%   W   Vector with length of 'n' representing one scenario of the Wiener
%       process
%   t   Time series
%
%WIENER_PROCESS(N,T,GRAPH,N_INSTANCES)
%   N   Number of subintervals of [0,T] (granularity)
%       Default value: 1000
%   T   Time interval [0,T]
%       Default value: 1
%   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
%   TPLOT   Number from (0,T] interval specific for plotting.
%
%Note:
%   When printing, only 1e3 of realizations are plotted due to the
%   execution time. But also the runs with max and min values are printed!
%
%Example usage: 
%   wiener_process(100, 1, true, 1000)

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

if tplot > T; error('Ratio of the data to be plotted has to be <= T'); end

%% Computation task
dt = T/n;
dW = sqrt(dt)*randn(n,n_instances);
dW(1,:) = 0;
W = cumsum(dW,1);
t = linspace(0, tplot, tplot/T*n);  % print only the ratio given by 'tplot'
st = length(t);
%% Plotting
if graph
    figure
    % Statistics
    maxW = max(W, [], 1); minW = min(W, [], 1);
    % Keep the indices for plotting
    % To plot the run with maximum and minimum value
    [maxW, imax] = max(maxW); [minW, imin] = min(minW);
    lastW = W(end,:);
    % Plot itself
    % Plot only reasonable number of instances (1000.)
    % If the max and min runs are in those first 1000, then don't plot
    % them again.
    n_to_print = min(1e3, n_instances);
    if imax <= n_to_print; imax = []; end
    if imin <= n_to_print; imin = []; end
    plot(t,W(1:st,[1:n_to_print imax imin])); hold on;
    % 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 Wiener process")
    xlabel("\it t")
    ylabel("\it W_t")
    grid on;
    message = "Minimal value of W: " + minW + newline + ...
              "Maximal value of W: " + maxW + newline + ...
              "Mean of values of W at time T: " + mean(lastW) + newline + ...
              "St. dev. of values of W at time T: " + std(lastW);
    annotation('textbox',[.1 .19 0 0],'String',message,'FitBoxToText','on')
end

