function [lastPos, nPasses, lastDir, visits] = portPassingSimulation(nSims, N, p)
% PORTPASSINGSIMULATION Random walk simulation of port around a table
%
%   lastPos = portPassingSimulation(nSims, N, p) runs nSims simulations of
%   passing port around a table with N people with probability p of passing
%   to the left and (1-p) to the right. We always start at person 0 so
%   there are N-1 possible ending spots (1:N-1). The return argument
%   lastPos is of length nSims and holds the value of the last position to
%   receive the port.
%
%   [lastPos, nPass, lastDir, visits] = portPassingSimulation(nSims, N, p)
%   returns not just the last position of the port, but also the number of
%   passes needed to reach lastPos (nPass) and the direction in which last
%   position was reached (lastDir). nPass and lastDir are the same size as
%   lastPos. In addition visits (of size 1xN) indicates how many times each
%   of the positions around the table was visited during all the
%   simulations. Note that visits is costly to compute so not requesting it
%   as an output will speed up the computation.

% Copyright 2017 The MathWorks, Inc.

assert(N>1, 'There must be more than one person at the table');
% Where is the port at the end of the last block of simulations (it always
% starts at position zero)
blockEndPosition = zeros(nSims, 1);
% For all simulations what is the next left-most position (column 1) and
% next right-most (column 2) position the port has reached? A simulation
% ends when the current position of the port is the same as both columns.
nextLimit = [ones(nSims, 1) zeros(nSims, 1)+(N-1)];
% How many of the simulations are still running?
nRunning = nSims;
% Outputs from the simulation that will be filled in as it progresses.
% 'lastPos' is the location that the port arrived at when all possible
% posistions had been visited, and 'nPasses' is the number of passes needed
% to reach that location, 'lastDir' the last direction the port move on
% arrival at lastPos.
lastPos = NaN(nSims, 1);
nPasses = NaN(nSims, 1);
lastDir = NaN(nSims, 1);
% visits is the total number of times each position at the table is
% visited. This is computationally intensive to add to the simulation so is
% only computed if requested
visits = zeros(1, N);
computeVisits = nargout > 3;

% Which pass of the simulation are we on? Updated at the end of the loop so
% starts at 1.
iteration = 1;
% Blocking value - how many passes do we simulate in a block to get
% vectorization of the simulation. I haven't spent much time optimizing
% this parameter - it might be better to have it depend on N in some way?
BLOCK_SIZE = 100;

% A brief note on performance - most of the time in the simulation should
% be in the calls to rand and cumsum below. All the rest should be trivial
% in comparison. If you profile this code you should see that (which is a
% result of using the BLOCK_SIZE vectorization approach).

% While there are simulation still not finished, keep running
while nRunning > 0
    % Generate a matrix that randomly chooses to pass left (+1) or right
    % (-1) with probabilities p and 1-p. This matrix simulates BLOCK_SIZE
    % passes across all the nNotFinished simulations. We don't want too
    % many simulations to finish each time but we do want to take advantage
    % of MATLAB's vectorization of rand & cumsum.
    d = (2*(rand(nRunning, BLOCK_SIZE) < p) - 1);
    % From the pass directions we need to deduce where the port actually
    % goes. We do this by cumulatively summing the pass matrix, starting at
    % the actual location of the port from the previous block. 
    cR = mod(cumsum(d, 2) + blockEndPosition, N);
    % Which simulations have not finished as we iterate over the loop? This
    % allows us to NOT reallocate memory for the cumulative cR matrix which
    % is a slow operation. As we loop over the actual passes we will update
    % this vector so that we never double count a finished simulation.
    notFinished = true(nRunning, 1);
    % Loop over each pass seeing which simulations finish on that pass
    for ii = 1:BLOCK_SIZE
        % If we are tracking the number of visits to each location around
        % the table then iterate over each of the simulations checking
        % which are finished, where not finished we simply increment 
        if computeVisits 
            for jj = 1:numel(notFinished)
                if notFinished(jj)
                    thisPos = cR(jj,ii) + 1;
                    visits(thisPos) = visits(thisPos) + 1;
                end
            end
        end
        % Deduce which simulations are first visiting one of their limits
        % as we will need to update those limits. In addition, where the
        % position updates both limits simultaneously we know this is the
        % last position to be visited
        nextLimitVisited = nextLimit == cR(:, ii);
        justFinished = all(nextLimitVisited, 2) & notFinished;
        % Note the singleton expansion multiplication by [1 -1] allows us
        % to increment column 1 and decrement column 2 of nextLimit in one
        % operation
        nextLimit = nextLimit + nextLimitVisited.*[1 -1];
        
        % For finished simulations we need to set output values and
        % indicate that those simulations are now finished
        nJustFinished = sum(justFinished);
        if nJustFinished > 0            
            nNowRunning = nRunning - nJustFinished;
            % Fill in the output variables and update notFinished
            lastPos(nNowRunning+1:nRunning) = cR(justFinished, ii);
            nPasses(nNowRunning+1:nRunning) = iteration;
            lastDir(nNowRunning+1:nRunning) = d(justFinished, ii);
            notFinished(justFinished) = false;
            % Update nRunning to indicate how many simulations are still
            % running (used also to fill output variables correctly)
            nRunning = nNowRunning; 
        end
        iteration = iteration + 1;
    end
    % After looping over all the passes in the simulation remove all the
    % finished simulations and set the starting location for the next
    % batch.
    blockEndPosition = cR(notFinished, end);
    nextLimit = nextLimit(notFinished, :);
end
% Ignore number of visits to location 0 (the head of the table where we
% started)
visits = visits(2:end);

