function augmentedStateTransitionMat = transformregularmatrixtoaugmented( trueStateTransitionMat, nTrueStates, nPages, nAugmentedStates )
    % Given a transition matrix in 'regular' form, return the equivalent
    % transition matrix in 'augmented' form, assuming the original states
    % are the 'true' ones.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % trueStateTransitionMat - 3-D matrix; first two dimensions hold square
    % transition matrices on 'true' state space, one for each page. Must
    % not be in log domain.
    %
    % augmentedStateTransitionMat - similar 3-D matrix, but pages are
    % 'augmented' state transition matrices. Is not in log domain.
    %
    % nAugmentedStates - number of augmented states. If not supplied, we
    % can work it out.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 06.10.2014 - copied from old version. Changed so we now permit
    % generator matrices with rowsums of zero.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if nargin < 2,
        nTrueStates = size(trueStateTransitionMat, 1);
    else
        if size(trueStateTransitionMat, 1) ~= nTrueStates,
            error('Number of rows of trueStateTransitionMat does not match nTrueStates!')
        end
    end
    if nargin < 3,
        nPages = size(trueStateTransitionMat, 3);
    else
        if size(trueStateTransitionMat, 3) ~= nPages,
            error('Number of pages of trueStateTransitionMat does not match nPages!')
        end
    end
    if nargin < 4,
        nAugmentedStates = nTrueStates * (nTrueStates + 1) / 2;
    end
    if size(trueStateTransitionMat, 2) ~= nTrueStates,
        error('Each page of trueStateTransitionMat must be a square matrix!')
    end
    
    augmentedStateTransitionMat = zeros([nAugmentedStates, nAugmentedStates, nPages]);

    % Form the triangle numbers plus one: the kth one is the first
    % augmented state corresponding to true state k.
    triNums = cumsum(0:nTrueStates) + 1;

    for iTrueState = 1:nTrueStates,
        % Augmented state indices associated with this true state.
        augStates = triNums(iTrueState):(triNums(iTrueState + 1) - 1);
        % Take them straight out of the true state transition matrix.
        augmentedStateTransitionMat(augStates, augStates, :) = trueStateTransitionMat(augStates - triNums(iTrueState) + 1, augStates - triNums(iTrueState) + 1, :);
        % This is the only other (non-zero probability) transition possible
        % from this source state: that which increases the augmentation
        % variable (number of states exhibited) by one and transitions the
        % true state to the same value.
        if iTrueState < nTrueStates,
            augmentedStateTransitionMat(augStates, triNums(iTrueState + 2) - 1, :) = sum(trueStateTransitionMat(augStates - triNums(iTrueState) + 1, (iTrueState + 1):nTrueStates, :), 2);
        end
    end
end