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);