function firstArrivals = formfirstarrivalsmatrix( stateTrajectories, nStates, maxNTransitions, nPages )
    % Given a matrix of multiple hidden state trajectories, build a matrix
    % of logical indicators of from which source state did each trajectory
    % first arrive at each destination state.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % firstArrivals - a logical 3-D matrix. First two dimensions hold
    % matrices corresponding to transition matrices, for which the (i, j)th
    % element (of page p) is true if and only if the first appearance of
    % state j follows an appearance of state i (in the trajectory stored in
    % column p of stateTrajectories).
    %
    % stateTrajectories - when using the augmented state space, this is a
    % matrix of the 'true' state trajectories.
    %
    % nStates - when using the augmented state space, this is the number of
    % true states.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 06.10.2014 - copied from old version. Changed dimensions of output to
    % (nStates, nStates - 1) to conform to requirements for sampling
    % generator matrices (never needed the first column anyway).
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if nargin < 2,
        error('Insufficient arguments!')
    end
    if nargin < 3,
        maxNTransitions = size(stateTrajectories, 1);
    else
        if size(stateTrajectories, 1) ~= maxNTransitions,
            error('maxNTransitions does not match the number of time bins implied by stateTrajectories!')
        end
    end
    if nargin < 4,
        nPages = size(stateTrajectories, 2);
    else
        if size(stateTrajectories, 2) ~= nPages,
            error('nPages does not match the number of pages implied by stateTrajectories!')
        end
    end
    
    firstArrivals = false([nStates, nStates - 1, nPages]);
    
    if (nStates == 1) || (maxNTransitions == 0),
        return;
    end
    
    trajectoriesColumnStartInds = 0:maxNTransitions:(maxNTransitions * (nPages - 1));
%     matrixStartInds = 0:(nStates ^ 2):((nStates ^ 2) * (nPages - 1));
    matrixStartInds = 0:(nStates * (nStates - 1)):((nStates * (nStates - 1)) * (nPages - 1));

    for iDestinationState = 2:nStates,
        % Get row vector of time bins for which this state was first
        % exhibited.
        [~, firstArrivalTimeBins] = max(stateTrajectories == iDestinationState, [], 1);
        % States that are never exhibited will have given us time bins of 1
        % in firstArrivalTimeBins. Since we don't care about destination
        % states at time bin 1 (they don't exist), we can simply ignore
        % these using these logical indices.
        validFirstArrivalLog = firstArrivalTimeBins > 1;
        % Get previous states by linear indexing of the appropriate columns
        % and rows of the trajectories matrix.
        prevStates = stateTrajectories(firstArrivalTimeBins(validFirstArrivalLog) + trajectoriesColumnStartInds(validFirstArrivalLog) - 1);
        % Linearly index the appropriate matrix elements; set them to true.
        firstArrivals(prevStates + matrixStartInds(validFirstArrivalLog) + (iDestinationState - 2) * nStates) = true;
    end
end