function logBackwardProbMat = backwardsmoothing( ModelSpec, logBackwardProbMat, logLikelihood, logTransitionMat, trajectoryMatBinInds, maxNTransitions, nStates, nParticles )
    % Perform the backward recursions of the forward-backward algorithm,
    % used for 'correcting' the fowrward probabilities and obtaining
    % smoothed posteriors. The likelihood vector at each timepoint must be
    % supplied. The transition matrix can either be "regular" - fully
    % parameterised - or "augmented", which is sparse and requires special
    % handling for faster computations. The default is "regular".
    % Recursions are performed simultaneously for any number of processes
    % (i.e. particles in the algorithm of Chopin).
    %
    % We are able to perform the recursions for particles on different time
    % points, and with different numbers of time points. Information about
    % which particles are involved at each time is supplied in
    % trajectoryMatBinInds.
    %
    % With this version we do not call a subfunction at each step of the
    % time loop. This improves performance slightly.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % logLikelihood is a matrix containing the likelihood vectors.
    %
    % logBackwardProbMat stores the preallocated backward probabilities
    % matrix (log domain). Should be preallocated - in which case MUST be
    % initialised with all zeros (i.e. the log of a probability of 1) - but
    % not essential. If not, then supply an empty matrix.
    %
    % logTransitionMat the log of the transition matrix. If it is 3D, it
    % must have the same number of pages as the other probability
    % structures.
    %
    % 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 2-D
    % matrix of logical indices, one for each timestep (rows) between 1 and
    % the upper limit for any particle. We assume all particles take part
    % in the initial state distribution, so this has one fewer rows than
    % the forward probabilities matrix. In each column there is a binary
    % indicator of inclusion or exclusion for that particle at that time
    % step.
    %
    % ModelSpec must contain the transMatNonzeroRowInds,
    % transMatNonzeroColInds and nNonzerosAugTransMat parameters if the
    % augmented version is used. See setupaugmentedstateindices().
    %
    % Note that in previous versions we had particleIndices and
    % nParticles as arguments for specifying which transition matrices
    % and how many pages of the forward probabilities matrix should be
    % used, according to the 'block size' (memory management idea). Now
    % that management should be performed BEFORE passing data to this
    % function, for added generality.
    %
    % Note: if multiple particles are being processed together, they must
    % all have the same dimension. The first way to 'block' particles is
    % according to dimension, and submit them to this function separately.
    %
    % The above two notes together mean that the preallocated forward
    % matrix supplied must be of the correct size for this dimension.
    %
    % ModelSpec must also contain augStateTransformationToMaxState.
    %
    % Note that logBackwardProbMat is supplied pre-initialised; the last
    % time step (for any particle stream) has probabilities of 1 in all
    % states.
    %
    % The convention for all our data is to have:
    % dim1 = time
    % dim2 = state
    % dim3 = particle
    %
    % or (for parameter structures)
    % dim1 = state/position
    % dim2 = state/neuron/position/1
    % dim3 = particle
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 27.01.2015 - copied from old version, backwardsmoothinginline.m, made
    % necessary changes.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if nargin < 4,
        error('Incorrect number of arguments!')
    end
    if nargin < 7,
        maxNTransitions = size(logLikelihood, 1);
    else
        if size(logLikelihood, 1) ~= maxNTransitions,
            error('Size of logLikelihood does not match maxNTransitions!')
        end
    end
    if nargin < 8,
        nStates = size(logTransitionMat, 1);
    else
        if size(logTransitionMat, 1) ~= nStates,
            error('n states of logInitialStateDist must be nStates!')
        end
    end
    if nargin < 9,
        nParticles = size(logLikelihood, 3);
    else
        if size(logLikelihood, 3) ~= nParticles,
            error('Size of logLikelihood does not match nParticles!')
        end
    end
    if nargin < 5,
        trajectoryMatBinInds = true([maxNTransitions, nParticles]);
    else
        if size(trajectoryMatBinInds, 1) ~= maxNTransitions,
            error('Time dimension of trajectoryMatBinInds and logBackwardProbMat must match!')
        end
        if size(trajectoryMatBinInds, 2) ~= nParticles,
            error('Particle dimension (2) of trajectoryMatBinInds and logBackwardProbMat (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
    % Better for memory management if this is preallocated.
    if isempty(logBackwardProbMat),
        logBackwardProbMat = zeros([maxNTransitions, nStates, nParticles]);
    end
    if size(logBackwardProbMat, 1) ~= maxNTransitions,
        error('Size of logLikelihood and logBackwardProbMat must match except for rows, logForwardProbMat should have one more!')
    end
    if size(logBackwardProbMat, 2) ~= nStates,
        error('Size of logLikelihood and logBackwardProbMat must match except for rows, logForwardProbMat should have one more!')
    end
    if size(logBackwardProbMat, 3) ~= nParticles,
        error('Size of logLikelihood and logBackwardProbMat must match except for rows, logForwardProbMat should have one more!')
    end
    if size(logTransitionMat, 3) ~= size(logLikelihood, 3),
        error('n pages of logTransitionMat must match that of other data structures!')
    end
    if size(logTransitionMat, 1) ~= size(logTransitionMat, 2),
        error('First two dimensions of logTransitionMat must be a square matrix!')
    end
%     if any(any(any(logBackwardProbMat ~= 0))),
%         error('logBackwardProbMat must be initialised with all 0!')
%     end
    if any(any(any(logTransitionMat > 0))),
        error('logTransitionMat contains log probabilities greater than 0!')
    end
%     if any(any(any(logLikelihood > 0))),
%         error('logLikelihood contains log probabilities greater than 0!')
%     end
    if ~islogical(ModelSpec.useAugStateModel),
        error('ModelSpec.useAugStateModel must be a logical!')
    end
    if ModelSpec.useAugStateModel,       
        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(logLikelihood, 2);
        if stateDim ~= ModelSpec.augStateComponents(nStates, 2),
            error('Size of logLikelihood does not match stateDim!')
        end
    end
    
    if maxNTransitions < 1,
        return;
    end
    
    if ~ModelSpec.useAugStateModel,
        for t = (maxNTransitions - 1):-1:1,
            % Multiply destination states' probabilities, and take account
            % of observation at destination time.
            backwardInter =  bsxfun(@plus, logTransitionMat(:, :, trajectoryMatBinInds(t + 1, :)), logLikelihood(t + 1, :, trajectoryMatBinInds(t + 1, :)) + logBackwardProbMat(t + 1, :, trajectoryMatBinInds(t + 1, :)));
            % Marginalise out destination state, leaving backward
            % probability at current time.
            logBackwardProbMat(t, :, trajectoryMatBinInds(t + 1, :)) = permute(sumlogprobmat(backwardInter, 2), [2, 1, 3]);
        end
%         % Do the last step for ALL particles, assuming that they will all
%         % be involved in first time step at least.
%         % Multiply destination states' probabilities, and take account of
%         % observation at destination time.
%         backwardInter = bsxfun(@plus, logTransitionMat, logLikelihood(2, :, :) + logBackwardProbMat(2, :, :));
%         % Marginalise out destination state, leaving backward probability
%         % at current time.
%         logBackwardProbMat(1, :, :) = permute(sumlogprobmat(backwardInter, 2), [2, 1, 3]);
    else
        % There is some setup necessary before we begin these recursions,
        % which depends on the number of particles, and the state dimension
        % they represent.
        % Indices for the summation via accumarray: two columns, first
        % lists the 'source' state, or row 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.generatorMatNonzeroRowInds{stateDim}, [nParticles, 1]), floor((1:(1/ModelSpec.nNonzerosAugGeneratorMat(stateDim)):(nParticles + 1 - 1/ModelSpec.nNonzerosAugGeneratorMat(stateDim)))')];
        accumarryIndexingMat = [repmat(ModelSpec.generatorMatNonzeroRowInds{stateDim}, [nParticles, 1]), reshape(repmat(1:nParticles, [ModelSpec.nNonzerosAugGeneratorMat(stateDim), 1]), [ModelSpec.nNonzerosAugGeneratorMat(stateDim) * nParticles, 1])];
        % The 'source state' associated with each linear-indexed non-zero
        % entry in the augmented transition matrix.
        transMatNonzeroRowInds = ModelSpec.generatorMatNonzeroRowInds{stateDim};
        % The 'destination state' associated with each linear-indexed
        % non-zero entry in the augmented transition matrix.
        transMatNonzeroColInds = ModelSpec.generatorMatNonzeroColInds{stateDim};
        % The linear indices of the non-zero entries of the augmented
        % transition matrix, stacked for each particle.
        particleIndAugmentationVector = (0:(nParticles - 1)) .* nStates .^ 2;
        transMatLinearIndsExtended = bsxfun(@plus, ModelSpec.generatorMatLinearInds{stateDim}, particleIndAugmentationVector);
        transMatLinearIndsExtended = reshape(transMatLinearIndsExtended, [ModelSpec.nNonzerosAugGeneratorMat(stateDim) * nParticles, 1]);
        % Preallocate this intermediary structure. It will hold the product
        % of the transition probabilities, in augmented format, with
        % 'destination' state backward probabilities from the next
        % timestep, and with the likelihood at the next time step.
        backwardInterMat = -inf([nStates, nStates, nParticles]);
        nParticlesIncludedEachTime = sum(trajectoryMatBinInds, 2);
        for t = (maxNTransitions - 1):-1:1,
            repeatedParticleIndsIncluded = 1:(nParticlesIncludedEachTime(t + 1) * ModelSpec.nNonzerosAugGeneratorMat(stateDim));
            transMatLinearIndsExtendedThisTime = transMatLinearIndsExtended(repeatedParticleIndsIncluded, :);

            backwardInterMatThisTime = backwardInterMat(:, :, trajectoryMatBinInds(t + 1, :));
            backwardInter1 = logBackwardProbMat(t + 1, transMatNonzeroColInds, trajectoryMatBinInds(t + 1, :));
            likelihoodInter = logLikelihood(t + 1, ModelSpec.augStateComponents(transMatNonzeroColInds, 1), trajectoryMatBinInds(t + 1, :));
            % Multiply destination states' probabilities and take account of the
            % next observation.
            logAugTransitionMat = logTransitionMat(:, :, trajectoryMatBinInds(t + 1, :));
            if nStates > 1,
                backwardInterMatThisTime(transMatLinearIndsExtendedThisTime) = logAugTransitionMat(transMatLinearIndsExtendedThisTime) + backwardInter1(:) + likelihoodInter(:);
            else
                backwardInterMatThisTime(transMatLinearIndsExtendedThisTime) = permute(logAugTransitionMat(transMatLinearIndsExtendedThisTime), [3, 2, 1]) + backwardInter1(:) + likelihoodInter(:);
            end
            % Prepare to sum over destination states.
            maxLogProbEachRow = max(backwardInterMatThisTime, [], 2);
            maxLogProbEachRow(maxLogProbEachRow == -Inf) = 0;
            maxLogProbEachRowNonzerosOnly = maxLogProbEachRow(transMatNonzeroRowInds, :, :);
            if nStates > 1,
                backwardInterMatThisTime(transMatLinearIndsExtendedThisTime) = backwardInterMatThisTime(transMatLinearIndsExtendedThisTime) - maxLogProbEachRowNonzerosOnly(:);
            else
                backwardInterMatThisTime(transMatLinearIndsExtendedThisTime) = permute(backwardInterMatThisTime(transMatLinearIndsExtendedThisTime), [3, 2, 1]) - maxLogProbEachRowNonzerosOnly(:);
            end
            % Sum over destination states.
            if nStates > 1,
                backwardInter2 = log(accumarray(accumarryIndexingMat(repeatedParticleIndsIncluded, :), exp(backwardInterMatThisTime(transMatLinearIndsExtendedThisTime)), [nStates, nParticlesIncludedEachTime(t + 1)])) + permute(maxLogProbEachRow, [1, 3, 2]);
            else
                backwardInter2 = log(accumarray(accumarryIndexingMat(repeatedParticleIndsIncluded, :), exp(permute(backwardInterMatThisTime(transMatLinearIndsExtendedThisTime), [3, 2, 1])), [nStates, nParticlesIncludedEachTime(t + 1)])) + permute(maxLogProbEachRow, [1, 3, 2]);
            end
            logBackwardProbMat(t, :, trajectoryMatBinInds(t + 1, :)) = permute(backwardInter2, [3, 1, 2]);
        end
%         % Do the last step for ALL particles, assuming that they will all
%         % be involved in first time step at least.
%         logBackwardProbMat(1, :, :) = backwardsmoothingsinglestepaugmented(logBackwardProbMat(2, :, :), logLikelihood(1, :, :), logTransitionMat, backwardInterMat, accumarryIndexingMat, transMatNonzeroRowInds, transMatNonzeroColInds, transMatLinearIndsExtended);
    end
end

