function logForwardProbVec = forwardfilteringsinglestepaugmented( ModelSpec, logForwardProbVecPrev, logLikelihoodVec, logTransitionMat, nStates, nParticles )
    % Perform a single step of the forward accumulation recursions, when
    % the transition matrix is of the augmented form (hence, sparse, and
    % requires special handling to minimise computational cost). We
    % potentially have multiple processes/particles to compute the forward
    % probabilities for simultaneously.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % logForwardProbVec will be a 3D matrix, but the first dimension is of
    % size 1 (time).
    %
    % logForwardProbVecPrev is the same, and contains the probabilities
    % from the previous timestep.
    %
    % logLikelihoodVec has the same structure.
    %
    % transitionMatSourceStates is a column vector of the source states as
    % they appear in a linear-indexed augmented transition matrix of the
    % appropriate dimension.
    %
    % generatorMatLinearIndsExtended is a column vector of the linear
    % indices required to select the non-zero entries of the augmented
    % transition matrix, one column at a time. Indices for different
    % particles are stacked.
    %
    % accumarryIndexingMat holds indices for the summation via accumarray:
    % two columns, first lists the 'destination' state, or column index of
    % the transition matrix, corresponding to each nonzero entry, second
    % lists particle indices (repeated for each nonzero entry of the
    % transition matrix).
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 15.10.2014 - copied from old version. Changed likelihood vector so
    % that it is now based on the 'true' state space. Brought computation
    % of intermediate structures in to function and removed them from
    % arguments.
    % 17.10.2014 - removed requirement for logLikelihood to be negative -
    % not necessary with contious observations.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargin < 4,
        error('Incorrect number of arguments!')
    end
    if nargin < 5,
        nStates = size(logForwardProbVecPrev, 2);
    else
        if size(logForwardProbVecPrev, 2) ~= nStates,
            error('Size of logForwardProbVecPrev does not match nStates!')
        end
    end
    if nargin < 6,
        nParticles = size(logForwardProbVecPrev, 3);
    else
        if size(logForwardProbVecPrev, 3) ~= nParticles,
            error('Size of logForwardProbVecPrev does not match nParticles!')
        end
    end
    if size(logLikelihoodVec, 3) ~= nParticles,
        error('Size of logLikelihoodVec does not match nParticles!')
    end
    if size(logTransitionMat, 3) ~= nParticles,
        error('n pages of logTransitionMat must match that of other data structures!')
    end
    if size(logTransitionMat, 1) ~= nStates,
        error('n states of logTransitionMat must match that of other data structures!')
    end
    if any(any(any(logTransitionMat > 0))),
        error('logTransitionMat contains log probabilities greater than 0!')
    end
%     if any(any(logLikelihoodVec > 0)),
%         error('logLikelihoodVec contains log probabilities greater than 0!')
%     end
    if ~isfield(ModelSpec, 'augStateComponents'),
        error('augStateComponents must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'generatorMatNonzeroRowInds'),
        error('generatorMatNonzeroRowInds must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'generatorMatNonzeroColInds'),
        error('generatorMatNonzeroColInds must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'nNonzerosAugGeneratorMat'),
        error('nNonzerosAugGeneratorMat must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'generatorMatLinearInds'),
        error('generatorMatLinearInds must be a field of ModelSpec!')
    end
    % Number of 'true' states even in 'augmented' model.
    stateDim = size(logLikelihoodVec, 2);
    if stateDim ~= ModelSpec.augStateComponents(nStates, 2),
        error('Size of logLikelihoodVec does not match stateDim!')
    end

    % Indices for the summation via accumarray: two columns, first lists
    % the 'destination' state, or column index of the transition matrix,
    % corresponding to each nonzero entry, second lists particle indices
    % (repeated for each nonzero entry of the transition matrix).
    accumarryIndexingMat = [repmat(ModelSpec.generatorMatNonzeroColInds{stateDim}, [nParticles, 1]), reshape(repmat(1:nParticles, [ModelSpec.nNonzerosAugGeneratorMat(stateDim), 1]), [ModelSpec.nNonzerosAugGeneratorMat(stateDim) * nParticles, 1])];
    % The linear indices of the non-zero entries of the augmented
    % transition matrix, stacked for each particle.
    particleIndAugmentationVector = (0:(nParticles - 1)) .* nStates .^ 2;
    generatorMatLinearIndsExtended = bsxfun(@plus, ModelSpec.generatorMatLinearInds{stateDim}, particleIndAugmentationVector);
    generatorMatLinearIndsExtended = reshape(generatorMatLinearIndsExtended, [ModelSpec.nNonzerosAugGeneratorMat(stateDim) * nParticles, 1]);
    % Preallocate this intermediary structure. It will hold the product of
    % the transition probabilities, in augmented format, with 'source'
    % state forward probabilities from the previous timestep.
    forwardInterMat = -inf([nStates, nStates, nParticles]);

    forwardInter1 = logForwardProbVecPrev(:, ModelSpec.generatorMatNonzeroRowInds{stateDim}, :);
    % Multiply source states' probabilities.
    if nStates > 1,
        forwardInterMat(generatorMatLinearIndsExtended) = logTransitionMat(generatorMatLinearIndsExtended) + forwardInter1(:);
    else
        forwardInterMat(generatorMatLinearIndsExtended) = permute(logTransitionMat(generatorMatLinearIndsExtended), [3, 2, 1]) + forwardInter1(:);
    end
    % Prepare to sum over source states.
    maxLogProbEachCol = max(forwardInterMat, [], 1);
    maxLogProbEachCol(maxLogProbEachCol == -Inf) = 0;
    maxLogProbEachColNonzerosOnly = maxLogProbEachCol(:, ModelSpec.generatorMatNonzeroColInds{stateDim}, :);
    if nStates > 1,
        forwardInterMat(generatorMatLinearIndsExtended) = forwardInterMat(generatorMatLinearIndsExtended) - maxLogProbEachColNonzerosOnly(:);
    else
        forwardInterMat(generatorMatLinearIndsExtended) = permute(forwardInterMat(generatorMatLinearIndsExtended), [3, 2, 1]) - maxLogProbEachColNonzerosOnly(:);
    end
    % Sum over source states.
    if nStates > 1,
        logSummationMat = forwardInterMat(generatorMatLinearIndsExtended);
        forwardInter2 = log(accumarray(accumarryIndexingMat, exp(logSummationMat), [nStates, nParticles])) + permute(maxLogProbEachCol, [2, 3, 1]);
    else
        logSummationMat = permute(forwardInterMat(generatorMatLinearIndsExtended), [3, 2, 1]);
        forwardInter2 = log(accumarray(accumarryIndexingMat, exp(logSummationMat), [nStates, nParticles])) + permute(maxLogProbEachCol, [2, 3, 1]);
    end
    % Take account of observations via the likelihood.
    logForwardProbVec = logLikelihoodVec(1, ModelSpec.augStateComponents(1:nStates, 1), :) + permute(forwardInter2, [3, 1, 2]);
end