function ModelSpec = setupaugmentedstateindices( ModelSpec )
    % Works out the indices necessary for accessing the nonzero entries in
    % an augmented state transition matrix. Included are the row indices
    % and column indices of these entries ordered by appearance when
    % linearly indexed, the source state associated with each nonzero entry
    % linearly indexed, also and similarly, the destination state, and the
    % number of nonzeros.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % We create the following fields of ModelSpec:
    % nNonzerosAugGeneratorMat, generatorMatNonzeroRowInds,
    % generatorMatNonzeroColInds, generatorMatLinearInds,
    % generatorMatNonzeroBinInds, generatorMatLowerZeroBinInds,
    % generatorMatNonzeroOrUpperRowInds, generatorMatNonzeroOrUpperColInds,
    % nNonzerosAugGeneratorMatUpper and generatorMatUpperLinearInds.
    %
    % nNonzerosAugGeneratorMat, a column vector, number of nonzeros for
    % each state dimension from 1 to the maximum allowed
    % (maxStateDimension).
    %
    % generatorMatNonzeroRowInds, generatorMatNonzeroRowInds, similar cell
    % arrays to above with column vectors of the row indices and column
    % indices of nonzero entries in order of access by linear indexing.
    % Used for deriving the linear indices.
    %
    % generatorMatLinearInds, similar cell array of column vectors, this
    % time of the linear indices of the non-zero entries of an augmented
    % transition matrix.
    %
    % generatorMatNonzeroBinInds is a logical matrix of the same dimensions
    % as an augmented transition matrix, and is true everywhere we expect a
    % non-zero element of the augmented transition matrix, and false
    % everywhere that is always zero in the augmented transition matrix.
    %
    % generatorMatLowerZeroBinInds is the negation of
    % generatorMatNonzeroBinInds in logical conjunction with row indices
    % being greater than the corresponding column indices. I.e., this
    % matrix identifies the zeros of the transition matrix in the lower
    % half.
    %
    % generatorMatNonzeroOrUpperRowInds, generatorMatNonzeroOrUpperColInds,
    % and generatorMatUpperLinearInds are similar to their similarly named
    % counterparts above, but additionally index those locations that would
    % be zero in the discrete-time Markov transition matrix and appear
    % above the diagonal. In the cts-time Markov process transition
    % matrices, the transitions corresponding to these locations are not
    % necessarily impossible.
    %
    % ModelSpec.maxNStates must be defined.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 01.10.2014 - copied from old version.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if ~isfield(ModelSpec, 'maxNStates'),
        error('maxNStates must be a field of ModelSpec!')
    end

    stateVec = (1:ModelSpec.maxNStates)';
    ModelSpec.nNonzerosAugGeneratorMat = stateVec .^ 2 + [0; cumsum(stateVec(1:(end - 1)))] + [0; cumsum(stateVec(1:(end - 1)) .^ 2)];
    
    ModelSpec.generatorMatNonzeroRowInds = cell(1, ModelSpec.maxNStates);
    ModelSpec.generatorMatNonzeroColInds = cell(1, ModelSpec.maxNStates);
    ModelSpec.generatorMatLinearInds = cell(1, ModelSpec.maxNStates);
    ModelSpec.generatorMatNonzeroBinInds = cell(1, ModelSpec.maxNStates);
    ModelSpec.generatorMatLowerZeroBinInds = cell(1, ModelSpec.maxNStates);
    ModelSpec.generatorMatNonzeroOrUpperRowInds = cell(1, ModelSpec.maxNStates);
    ModelSpec.generatorMatNonzeroOrUpperColInds = cell(1, ModelSpec.maxNStates);
    ModelSpec.generatorMatUpperLinearInds = cell(1, ModelSpec.maxNStates);
    for iStateDim = 1:ModelSpec.maxNStates,
        augStateDim = iStateDim * (iStateDim + 1) / 2;
        dummyGeneratorMat = 1 / iStateDim * (1 - eye(iStateDim)) - 1 / iStateDim * eye(iStateDim);
        augmentedDummyGeneratorMat = formaugmentedgeneratormatrix3d( dummyGeneratorMat, 1, iStateDim, augStateDim );

        [ModelSpec.generatorMatNonzeroRowInds{iStateDim}, ModelSpec.generatorMatNonzeroColInds{iStateDim}] = find(augmentedDummyGeneratorMat);

        ModelSpec.generatorMatLinearInds{iStateDim} = ModelSpec.generatorMatNonzeroRowInds{iStateDim} + (ModelSpec.generatorMatNonzeroColInds{iStateDim} - 1) * augStateDim;
        
        ModelSpec.generatorMatNonzeroBinInds{iStateDim} = augmentedDummyGeneratorMat ~= 0;
        
        lowerTransMatLogInds = bsxfun(@gt, (1:augStateDim)', 1:augStateDim);
        ModelSpec.generatorMatLowerZeroBinInds{iStateDim} = ~ModelSpec.generatorMatNonzeroBinInds{iStateDim} & lowerTransMatLogInds;
        
        % Additional indexing structures for the cts-time Markov process
        % model (some transitions not possible in the discrete time model
        % are now possible).
        upperGeneratorMatLogInds = bsxfun(@lt, (1:augStateDim)', 1:augStateDim);
        augmentedDummyGeneratorMat(upperGeneratorMatLogInds) = 1 / iStateDim;
        
        [ModelSpec.generatorMatNonzeroOrUpperRowInds{iStateDim}, ModelSpec.generatorMatNonzeroOrUpperColInds{iStateDim}] = find(augmentedDummyGeneratorMat);
        
        ModelSpec.generatorMatUpperLinearInds{iStateDim} = ModelSpec.generatorMatNonzeroOrUpperRowInds{iStateDim} + (ModelSpec.generatorMatNonzeroOrUpperColInds{iStateDim} - 1) * augStateDim;
    end
    
    ModelSpec.nNonzerosAugGeneratorMatUpper = cellfun(@(iMat) sum(sum(~iMat)), ModelSpec.generatorMatLowerZeroBinInds)';
end