function stateSampleVec = backwardsamplingsinglestepaugmented( stateSampleVecNext, logForwardProbVec, logSamplingDistributionMat, logAugTransitionMat, generatorMatNonzeroBinInds, trajectoryIndexingMat, particlesAlreadyInitialisedBinInds, nStates, nParticles )
    % A single step of the backward iterative sampling of state
    % trajectories when the transition matrix is for augmented state space
    % (hence special sparse structure handling is required). Can handle
    % multiple trajectories being sampled simultaneously. Some of these
    % trajectories will be starting at this time, some will have started at
    % later times (previous iterations). For the latter case, the sample
    % from one step in the future must be supplied.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % stateSampleVec and stateSampleVecNext are row vectors, with columns
    % indexed by particle (copies) to sample from. The latter are the state
    % samples from the next timestep (previous iteration).
    %
    % logForwardProbVec is a matrix with one rows index by state (the
    % forward probabilities over states at this timestep), and columns
    % indexed by particle (copies) to be sampled from.
    %
    % logSamplingDistributionMat is used to store an intermediate. Should
    % be initialised with all -Inf (log probability of 0) and be nStates by
    % total number of particles' copies.
    %
    % particlesAlreadyInitialisedBinInds says, of those particles involved
    % at THIS time, t, were they also involved at t + 1?
    %
    % generatorMatNonzeroBinInds is a matrix of the logical indices of the
    % nonzero locations in a single augmented transition matrix.
    %
    % trajectoryIndexingMat is a matrix the same size as the
    % logSamplingDistributionMat (states index rows, particles (copies)
    % index columns), but with each column filled with the label of that
    % particle/trajectory.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 14.10.2014 - copied from old version.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargin < 6,
        error('Incorrect number of arguments!')
    end
    if nargin < 8,
        nStates = size(logForwardProbVec, 1);
    else
        if size(logForwardProbVec, 1) ~= nStates,
            error('Size of logForwardProbVec does not match nStates!')
        end
    end
    if nargin < 9,
        nParticles = size(logForwardProbVec, 2);
    else
        if size(logForwardProbVec, 2) ~= nParticles,
            error('Size of logForwardProbVec does not match nParticles!')
        end
    end
    if size(logAugTransitionMat, 3) ~= nParticles,
        error('n pages of logAugTransitionMat must match that of other data structures!')
    end
    if size(logAugTransitionMat, 1) ~= nStates,
        error('n states of logAugTransitionMat must match that of other data structures, and source states must be in dimension 1!')
    end
    if size(logAugTransitionMat, 2) ~= nStates,
        error('n states of logAugTransitionMat must match that of other data structures, and destination states must be in dimension 2!')
    end
    if any(any(any(logAugTransitionMat > 0))),
        error('logAugTransitionMat contains log probabilities greater than 0!')
    end
%     if any(any(logForwardProbVec > 0)),
%         error('logForwardProbVec contains log probabilities greater than 0!')
%     end

    % Better for memory management if this is preallocated.
    if isempty(logSamplingDistributionMat),
        logSamplingDistributionMat = -inf([nStates, nParticles]);
    elseif ~isequal(size(logSamplingDistributionMat), [nStates, nParticles]),
        error('logSamplingDistributionMat must be a 2D matrix with states indexing rows and particle copies (trajectories requires) indexing columns!')
    end

    if any(any(~islogical(generatorMatNonzeroBinInds))),
        error('generatorMatNonzeroBinInds must be a logical matrix!')
    end
    if ~isequal(size(generatorMatNonzeroBinInds), [nStates, nStates]),
        error('generatorMatNonzeroBinInds must be a 2D matrix with the same dimensions as a transition matrix!')
    end
    if ~isequal(size(trajectoryIndexingMat), [nStates, nParticles]),
        error('trajectoryIndexingMat must be a 2D matrix with states indexing rows and particle copies (trajectories requires) indexing columns!')
    end

    if nargin < 7,
        particlesAlreadyInitialisedBinInds = true([1, nParticles]);
    else
        if size(particlesAlreadyInitialisedBinInds, 1) ~= 1,
            error('Time dimension of particlesAlreadyInitialisedBinInds must be 1!')
        end
        if size(particlesAlreadyInitialisedBinInds, 2) ~= nParticles,
            error('Particle dimension (2) of particlesAlreadyInitialisedBinInds and logForwardProbMat (3) must match!')
        end
        if any(any(~islogical(particlesAlreadyInitialisedBinInds))),
            error('particlesAlreadyInitialisedBinInds must be a logical matrix!')
        end
        if any(any(particlesAlreadyInitialisedBinInds(2:end, :) > particlesAlreadyInitialisedBinInds(1:(end - 1), :))),
            error('Any particles excluded from a time point (row of particlesAlreadyInitialisedBinInds) must not be included in later time points!')
        end
    end

    % Condition the transition probabilities on the sampled destination.
    % First, those trajectories that require initialisation.
    if any(~particlesAlreadyInitialisedBinInds),
        logSamplingDistributionMat(:, ~particlesAlreadyInitialisedBinInds) = logForwardProbVec(:, ~particlesAlreadyInitialisedBinInds);
        % Normalise these.
%         try
        logSamplingDistributionMat(:, ~particlesAlreadyInitialisedBinInds) = exp(normaliselogdistributionsmatrix( logSamplingDistributionMat(:, ~particlesAlreadyInitialisedBinInds), 1 ));
%         catch err
%         logSamplingDistributionMat
%         error(err.message)
%         error(err.stack)
%         end
    end
    % Second, those trajectories that have already been initialised.
    if any(particlesAlreadyInitialisedBinInds),
        samplingDistributionMatInter = logSamplingDistributionMat(:, particlesAlreadyInitialisedBinInds);
        % Some clever indexing is required to get the right columns of the
        % transition matrix from the appropriate pages (particles).
        % First restrict to the particles we are interested in here.
        logAugTransitionMat = logAugTransitionMat(:, :, particlesAlreadyInitialisedBinInds);
        % Form a vector of column indices corresponding to the destinations
        % we are conditioning on.
        nThisTimeAndNext = sum(particlesAlreadyInitialisedBinInds);
        transitionMatIndVec = stateSampleVecNext(particlesAlreadyInitialisedBinInds) + (0:(nThisTimeAndNext - 1)) * nStates;
        % Pick out those columns.
        logAugTransitionMat = reshape(logAugTransitionMat, [nStates, nStates * nThisTimeAndNext]);
        logAugTransitionMat = logAugTransitionMat(:, transitionMatIndVec);
        % This forms a matrix of the logical indices of the nonzeros in
        % columns of an augmented transition matrix of the correct size,
        % corresponding to those destination states previously sampled
        % (next time step).
        nonzerosConditionedOnDestinationLogInds = generatorMatNonzeroBinInds(:, stateSampleVecNext(particlesAlreadyInitialisedBinInds));
        logForwardProbVec = logForwardProbVec(:, particlesAlreadyInitialisedBinInds);
        samplingDistributionMatInter(nonzerosConditionedOnDestinationLogInds) = logForwardProbVec(nonzerosConditionedOnDestinationLogInds) + logAugTransitionMat(nonzerosConditionedOnDestinationLogInds);
        % Normalise these.
        maxLogProbEachCol = max(samplingDistributionMatInter, [], 1);
        maxLogProbEachCol(maxLogProbEachCol == -Inf) = 0;
        maxLogProbEachColNonzerosOnly = (maxLogProbEachCol(:, trajectoryIndexingMat(nonzerosConditionedOnDestinationLogInds)))';
        if nStates > 1,
            % Normalising constant for each column.
            logNormalisingConstantsVec = log(accumarray(trajectoryIndexingMat(nonzerosConditionedOnDestinationLogInds), exp(samplingDistributionMatInter(nonzerosConditionedOnDestinationLogInds) - maxLogProbEachColNonzerosOnly))) + permute(maxLogProbEachCol, [2, 1]);
            samplingDistributionMatInter(nonzerosConditionedOnDestinationLogInds) = samplingDistributionMatInter(nonzerosConditionedOnDestinationLogInds) - logNormalisingConstantsVec(trajectoryIndexingMat(nonzerosConditionedOnDestinationLogInds));
        else
            % Normalising constant for each column.
            logNormalisingConstantsVec = log(accumarray((trajectoryIndexingMat(nonzerosConditionedOnDestinationLogInds))', exp((samplingDistributionMatInter(nonzerosConditionedOnDestinationLogInds))' - maxLogProbEachColNonzerosOnly))) + permute(maxLogProbEachCol, [2, 1]);
            samplingDistributionMatInter(nonzerosConditionedOnDestinationLogInds) = (samplingDistributionMatInter(nonzerosConditionedOnDestinationLogInds))' - logNormalisingConstantsVec(trajectoryIndexingMat(nonzerosConditionedOnDestinationLogInds));
        end
        if any(any(isnan(samplingDistributionMatInter))),
            error('NaNs in sampling distribution!')
        end
        logSamplingDistributionMat(:, particlesAlreadyInitialisedBinInds) = exp(samplingDistributionMatInter);
    end

    stateSampleVec = samplefromcategoricaldistributions(logSamplingDistributionMat, 1, [nStates, nParticles]);
%     stateSampleVec = samplefromcategoricaldistributions(logSamplingDistributionMat, 1);
end