function stateSampleVec = backwardsamplingsinglestepregular( stateSampleVecNext, logForwardProbVec, logSamplingDistributionMat, logTransitionMat, particlesAlreadyInitialisedBinInds, nStates, nParticles )
    % A single step of the backward iterative sampling of state
    % trajectories when the transition matrix is of 'regular' form (no
    % special sparse structure handling 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?
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 14.10.2014 - copied from old version.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
    if nargin < 4,
        error('Incorrect number of arguments!')
    end
    if nargin < 6,
        nStates = size(logForwardProbVec, 1);
    else
        if size(logForwardProbVec, 1) ~= nStates,
            error('Size of logForwardProbVec does not match nStates!')
        end
    end
    if nargin < 7,
        nParticles = size(logForwardProbVec, 2);
    else
        if size(logForwardProbVec, 2) ~= nParticles,
            error('Size of logForwardProbVec does not match nParticles!')
        end
    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, and source states must be in dimension 1!')
    end
    if size(logTransitionMat, 2) ~= nStates,
        error('n states of logTransitionMat must match that of other data structures, and destination states must be in dimension 2!')
    end
    if any(any(any(logTransitionMat > 0))),
        error('logTransitionMat 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 nargin < 5,
        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);
    end
    % Second, those trajectories that have already been initialised.
    if any(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.
        logTransitionMat = logTransitionMat(:, :, particlesAlreadyInitialisedBinInds);
        nThisTimeAndNext = sum(particlesAlreadyInitialisedBinInds);
        % Form linear indices corresponding to each column we are
        % interested in.
        transitionMatIndVec = bsxfun(@plus, (1:nStates)', (stateSampleVecNext(particlesAlreadyInitialisedBinInds) - 1) * nStates + (0:(nThisTimeAndNext - 1)) * nStates^2);
        % Get back in the necessary shape for matrix addition.
        logTransitionMatConditionalOnDestination = reshape(logTransitionMat(transitionMatIndVec(:)), [nStates, nThisTimeAndNext]);
        logSamplingDistributionMat(:, particlesAlreadyInitialisedBinInds) = bsxfun(@plus, logForwardProbVec(:, particlesAlreadyInitialisedBinInds), logTransitionMatConditionalOnDestination);
    end
    % Normalise and sample.
    logSamplingDistributionMat = exp(normaliselogdistributionsmatrix( logSamplingDistributionMat, 1 ));
    stateSampleVec = samplefromcategoricaldistributions(logSamplingDistributionMat, 1);
end