function stateTrajectoryMat = backwardsampling( ModelSpec, logForwardProbMat, stateTrajectoryMat, logTransitionMat, trajectoryMatBinInds, maxNTransitions, nStates, nParticles )
    % Perform the stochastic backward recursions given the already-computed
    % forward probabilities and jump times. Multiple trajectories can be
    % simulated simultaneously, i.e. multiple particles. The forward
    % probabilities can be computed on varying numbers of jumps for each
    % particle, and the sampled states will occur on the same jumps.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % stateTrajectoryMat - 2D matrix with number of rows equal to the
    % greatest number of timesteps required of any particle, and number of
    % columns equal to the total number of samples required. This should be
    % preallocated, hence it is a parameter, but not essential. If not
    % preallocated it should be an empty matrix.
    %
    % logForwardProbMat - the precomputed forward probabilities matrix (log
    % domain).
    %
    % ModelSpec.useAugStateModel - logical scalar indicating whether or not
    % we must use the augmented indexing of the transition matrix (that it
    % is on an augmented state space and therefore requires special
    % treatment).
    %
    % trajectoryMatBinInds is optional; if supplied it is a 2D
    % matrix of logical indices, one for each timestep (rows) between 1 and
    % the upper limit for any particle. In each column there is a binary
    % indicator of inclusion or exclusion for that particle at that time
    % step.
    %
    % [particleIndsByCopyNumber - a column vector of particle indices
    % repeated the number of times each particle is required to be sampled
    % from.] OBSOLETE:
    % Do this 'expansion' of structures so that we sample the correct
    % number of times for each particle outside of this function.
    % logForwardProbMatByCopyNumber = logForwardProbMat(:, :, particleIndsByCopyNumber);
    % logTransitionMatByCopyNumber = logTransitionMat(:, particleIndsByCopyNumber);
    % trajectoryMatBinIndsByCopyNumber = trajectoryMatBinInds(:, particleIndsByCopyNumber);
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 14.10.2014 - copied from old version, backwardsamplingshortinline.m.
    % Changed names and format of parameters slightly.
    % 16.10.2014 - changed stateTrajectoryMat to have NaN in unused
    % positions i.e. where a particle does not make as many jumps as some
    % other particle and therefore does not use some rows.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargin < 4,
        error('Incorrect number of arguments!')
    end
    if nargin < 6,
        maxNTransitions = size(logForwardProbMat, 1);
    else
        if size(logForwardProbMat, 1) ~= maxNTransitions,
            error('Size of logForwardProbMat does not match maxNTransitions!')
        end
    end
    if nargin < 7,
        nStates = size(logForwardProbMat, 2);
    else
        if size(logForwardProbMat, 2) ~= nStates,
            error('Size of logForwardProbMat does not match nStates!')
        end
    end
    if nargin < 8,
        nParticles = size(logForwardProbMat, 3);
    else
        if size(logForwardProbMat, 3) ~= nParticles,
            error('Size of logForwardProbMat does not match nParticles!')
        end
    end
    % Better for memory management if this is preallocated.
    if isempty(stateTrajectoryMat),
        stateTrajectoryMat = nan([maxNTransitions, nParticles]);
    elseif ~isequal(size(stateTrajectoryMat), [maxNTransitions, nParticles]),
        error('stateTrajectoryMat must be a 2-D matrix with the time dimension of logForwardProbMat indexing rows and particles indexing columns!')
    end

    if ~islogical(ModelSpec.useAugStateModel),
        error('ModelSpec.useAugStateModel must be a logical!')
    end
    if ModelSpec.useAugStateModel,
        if ~isfield(ModelSpec, 'generatorMatNonzeroBinInds'),
            error('generatorMatNonzeroBinInds must be a field of ModelSpec!')
        end
        if ~isfield(ModelSpec, 'augStateComponents'),
            error('augStateComponents must be a field of ModelSpec!')
        end
        stateDim = ModelSpec.augStateComponents(nStates, 2);
    end
    if nargin < 5,
        trajectoryMatBinInds = true([maxNTransitions, nParticles]);
    else
        if size(trajectoryMatBinInds, 1) ~= maxNTransitions,
            error('Time dimension of trajectoryMatBinInds and logForwardProbMat must match!')
        end
        if size(trajectoryMatBinInds, 2) ~= nParticles,
            error('Particle dimension (2) of trajectoryMatBinInds and logForwardProbMat (3) must match!')
        end
        if any(any(~islogical(trajectoryMatBinInds))),
            error('trajectoryMatBinInds must be a logical matrix!')
        end
        if any(any(trajectoryMatBinInds(2:end, :) > trajectoryMatBinInds(1:(end - 1), :))),
            error('Any particles excluded from a time point (row of trajectoryMatBinInds) must not be included in later time points!')
        end
    end
       
    % It might be better to preallocate this outside, but not a big deal.
    logSamplingDistributionMat = -inf([nStates, nParticles]);
    
    % This makes some later computations more straightforward.
    logForwardProbMat = permute(logForwardProbMat, [2, 3, 1]);

    if ~ModelSpec.useAugStateModel,
        % A logical indexing vector that says, of those particles
        % involved at THIS time, t, were they also involved at t + 1? (At
        % initialisation, none of them were.)
        nParticlesEachTime = sum(trajectoryMatBinInds, 2);
        nParticlesAtEnd = nParticlesEachTime(end);
        particlesAlreadyInitialisedBinInds = false([1, nParticlesAtEnd]);
        stateTrajectoryMat(end, trajectoryMatBinInds(end, :)) = backwardsamplingsinglestepregular([], logForwardProbMat(:, trajectoryMatBinInds(end, :), end), logSamplingDistributionMat(:, trajectoryMatBinInds(end, :)), logTransitionMat(:, :, trajectoryMatBinInds(end, :)), particlesAlreadyInitialisedBinInds, nStates, nParticlesAtEnd);
        for t = (maxNTransitions - 1):-1:2,
            % A logical indexing vector that says, of those particles
            % involved at THIS time, t, were they also involved at t + 1?
            particlesAlreadyInitialisedBinInds = trajectoryMatBinInds(t, trajectoryMatBinInds(t, :)) & trajectoryMatBinInds(t + 1, trajectoryMatBinInds(t, :));

            logSamplingDistributionMatThisTime = logSamplingDistributionMat(:, trajectoryMatBinInds(t, :));
            logForwardProbVec = logForwardProbMat(:, trajectoryMatBinInds(t, :), t);
            % Condition the transition probabilities on the sampled
            % destination.
            % First, those trajectories that require initialisation.
            if any(~particlesAlreadyInitialisedBinInds),
                logSamplingDistributionMatThisTime(:, ~particlesAlreadyInitialisedBinInds) = logForwardProbVec(:, ~particlesAlreadyInitialisedBinInds);
            end
            % Second, those trajectories that have already been initialised.
            if any(particlesAlreadyInitialisedBinInds),
                stateSampleVecNext = stateTrajectoryMat(t + 1, trajectoryMatBinInds(t, :));
                % 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.
%                 logTransitionMatThisTime = logTransitionMat(:, :, particlesAlreadyInitialisedBinInds);
                logTransitionMatThisTime = logTransitionMat(:, :, trajectoryMatBinInds(t, :));
                logTransitionMatThisTime = logTransitionMatThisTime(:, :, particlesAlreadyInitialisedBinInds);
                nParticlesAlreadyInitialised = sum(particlesAlreadyInitialisedBinInds);
                % Form linear indices corresponding to each column we are
                % interested in.
                transitionMatIndVec = bsxfun(@plus, (1:nStates)', (stateSampleVecNext(particlesAlreadyInitialisedBinInds) - 1) * nStates + (0:(nParticlesAlreadyInitialised - 1)) * nStates ^ 2);
                % Get back in the necessary shape for matrix addition.
                logTransitionMatConditionalOnDestination = reshape(logTransitionMatThisTime(transitionMatIndVec(:)), [nStates, nParticlesAlreadyInitialised]);
                logSamplingDistributionMatThisTime(:, particlesAlreadyInitialisedBinInds) = bsxfun(@plus, logForwardProbVec(:, particlesAlreadyInitialisedBinInds), logTransitionMatConditionalOnDestination);
            end
            % Normalise and sample.
            logSamplingDistributionMatThisTime = exp(normaliselogdistributionsmatrix( logSamplingDistributionMatThisTime, 1 ));
            stateTrajectoryMat(t, trajectoryMatBinInds(t, :)) = samplefromcategoricaldistributions(logSamplingDistributionMatThisTime, 1, [nStates, nParticlesEachTime(t)]);
        end
        if maxNTransitions > 1,
            % Do the last step for ALL particles, assuming that they will
            % all be involved at time step 1 at least.
            particlesAlreadyInitialisedBinInds = trajectoryMatBinInds(2, :);
            stateTrajectoryMat(1, :) = backwardsamplingsinglestepregular(stateTrajectoryMat(2, :), logForwardProbMat(:, :, 1), logSamplingDistributionMat, logTransitionMat, particlesAlreadyInitialisedBinInds, nStates, nParticles);
        end
    else
        % A matrix with columns filled with the labels of the corresponding
        % particle/trajectory. We will only use as many columns as
        % necessary each time step, as dictated by the number of particles
        % involved indicator.
        trajectoryIndexingMat = repmat(1:nParticles, [nStates, 1]);
        % A logical indexing vector that says, of those particles
        % involved at THIS time, t, were they also involved at t + 1? (At
        % initialisation, none of them were.)
        nParticlesInvolvedEachTime = sum(trajectoryMatBinInds, 2);
%         nParticlesInvolved = sum(trajectoryMatBinInds(end, :));
        particlesAlreadyInitialisedBinInds = false([1, nParticlesInvolvedEachTime(end)]);
        stateTrajectoryMat(end, trajectoryMatBinInds(end, :)) = backwardsamplingsinglestepaugmented([], logForwardProbMat(:, trajectoryMatBinInds(end, :), end), logSamplingDistributionMat(:, trajectoryMatBinInds(end, :)), logTransitionMat(:, :, trajectoryMatBinInds(end, :)), ModelSpec.generatorMatNonzeroBinInds{stateDim}, trajectoryIndexingMat(:, 1:nParticlesInvolvedEachTime(end)), particlesAlreadyInitialisedBinInds, nStates, nParticlesInvolvedEachTime(end));
        for t = (maxNTransitions - 1):-1:2,
%             nParticlesInvolved = sum(trajectoryMatBinInds(t, :));
            trajectoryIndexingMatThisTime = trajectoryIndexingMat(:, 1:nParticlesInvolvedEachTime(t));
            % A logical indexing vector that says, of those particles
            % involved at THIS time, t, were they also involved at t + 1?
            particlesAlreadyInitialisedBinInds = trajectoryMatBinInds(t, trajectoryMatBinInds(t, :)) & trajectoryMatBinInds(t + 1, trajectoryMatBinInds(t, :));
            
            logSamplingDistributionMatThisTime = logSamplingDistributionMat(:, trajectoryMatBinInds(t, :));
            logForwardProbVec = logForwardProbMat(:, trajectoryMatBinInds(t, :), t);
            % Condition the transition probabilities on the sampled destination.
            % First, those trajectories that require initialisation.
            if any(~particlesAlreadyInitialisedBinInds),
                logSamplingDistributionMatThisTime(:, ~particlesAlreadyInitialisedBinInds) = logForwardProbVec(:, ~particlesAlreadyInitialisedBinInds);
                % Normalise these.
                logSamplingDistributionMatThisTime(:, ~particlesAlreadyInitialisedBinInds) = exp(normaliselogdistributionsmatrix( logSamplingDistributionMatThisTime(:, ~particlesAlreadyInitialisedBinInds), 1 ));
            end
            % Second, those trajectories that have already been
            % initialised.
            if any(particlesAlreadyInitialisedBinInds),
                stateSampleVecNext = stateTrajectoryMat(t + 1, trajectoryMatBinInds(t, :));
                logAugTransitionMatThisTime = logTransitionMat(:, :, trajectoryMatBinInds(t, :));
                samplingDistributionMatInter = logSamplingDistributionMatThisTime(:, 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.
                logAugTransitionMatThisTime = logAugTransitionMatThisTime(:, :, particlesAlreadyInitialisedBinInds);
                % Form a vector of column indices corresponding to the
                % destinations we are conditioning on.
                nParticlesAlreadyInitialised = sum(particlesAlreadyInitialisedBinInds);
                transitionMatIndVec = stateSampleVecNext(particlesAlreadyInitialisedBinInds) + (0:(nParticlesAlreadyInitialised - 1)) * nStates;
                % Pick out those columns.
                logAugTransitionMatThisTime = reshape(logAugTransitionMatThisTime, [nStates, nStates * nParticlesAlreadyInitialised]);
                logAugTransitionMatThisTime = logAugTransitionMatThisTime(:, 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).
                nonzerosConditionedOnDestinationBinInds = ModelSpec.generatorMatNonzeroBinInds{stateDim}(:, stateSampleVecNext(particlesAlreadyInitialisedBinInds));
                logForwardProbVec = logForwardProbVec(:, particlesAlreadyInitialisedBinInds);
                samplingDistributionMatInter(nonzerosConditionedOnDestinationBinInds) = logForwardProbVec(nonzerosConditionedOnDestinationBinInds) + logAugTransitionMatThisTime(nonzerosConditionedOnDestinationBinInds);
                % Normalise these.
                maxLogProbEachCol = max(samplingDistributionMatInter, [], 1);
                maxLogProbEachCol(maxLogProbEachCol == -Inf) = 0;
                maxLogProbEachColNonzerosOnly = (maxLogProbEachCol(:, trajectoryIndexingMatThisTime(nonzerosConditionedOnDestinationBinInds)))';
                if nStates > 1,
                    % Normalising constant for each column.
                    logNormalisingConstantsVec = log(accumarray(trajectoryIndexingMatThisTime(nonzerosConditionedOnDestinationBinInds), exp(samplingDistributionMatInter(nonzerosConditionedOnDestinationBinInds) - maxLogProbEachColNonzerosOnly))) + permute(maxLogProbEachCol, [2, 1]);
                    samplingDistributionMatInter(nonzerosConditionedOnDestinationBinInds) = samplingDistributionMatInter(nonzerosConditionedOnDestinationBinInds) - logNormalisingConstantsVec(trajectoryIndexingMatThisTime(nonzerosConditionedOnDestinationBinInds));
                else
                    % Normalising constant for each column.
                    logNormalisingConstantsVec = log(accumarray((trajectoryIndexingMatThisTime(nonzerosConditionedOnDestinationBinInds))', exp((samplingDistributionMatInter(nonzerosConditionedOnDestinationBinInds))' - maxLogProbEachColNonzerosOnly))) + permute(maxLogProbEachCol, [2, 1]);
                    samplingDistributionMatInter(nonzerosConditionedOnDestinationBinInds) = (samplingDistributionMatInter(nonzerosConditionedOnDestinationBinInds))' - logNormalisingConstantsVec(trajectoryIndexingMatThisTime(nonzerosConditionedOnDestinationBinInds));
                end
                if any(any(isnan(samplingDistributionMatInter))),
                    error('NaNs in sampling distribution!')
                end
                logSamplingDistributionMatThisTime(:, particlesAlreadyInitialisedBinInds) = exp(samplingDistributionMatInter);
            end

            stateTrajectoryMat(t, trajectoryMatBinInds(t, :)) = samplefromcategoricaldistributions(logSamplingDistributionMatThisTime, 1, [nStates, nParticlesInvolvedEachTime(t)]);
%             stateTrajectoryMat(t, trajectoryMatBinInds(t, :)) = samplefromcategoricaldistributions(logSamplingDistributionMatThisTime, 1);
        end
        if maxNTransitions > 1,
            % Do the last step for ALL particles, assuming that they will
            % all be involved at time step 1 at least.
            particlesAlreadyInitialisedBinInds = trajectoryMatBinInds(2, :);
            stateTrajectoryMat(1, :) = backwardsamplingsinglestepaugmented(stateTrajectoryMat(2, :), logForwardProbMat(:, :, 1), logSamplingDistributionMat, logTransitionMat, ModelSpec.generatorMatNonzeroBinInds{stateDim}, trajectoryIndexingMat, particlesAlreadyInitialisedBinInds, nStates, nParticles);
        end
    end
end