function logForwardProbMat = forwardfiltering( ModelSpec, logForwardProbMat, logLikelihood, logTransitionMat, logInitialStateDist, trajectoryMatBinInds, maxNTransitions, nStates, nParticles )
    % Perform the forward recursions for accumulating data for filtering.
    % This version permits an arbitrary discretisation of time, i.e. in the
    % cts-time (MJP) model we condition on jump times, which give us the
    % times at which we compute forward probabilities.
    %
    % 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).
    %
    % Information about which particles are involved at each time is
    % supplied in trajectoryMatBinInds.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % logLikelihood is a matrix containing the likelihood vectors: columns
    % are indexed by state, giving the likelihood function on each state at
    % each time, which indexes rows. On 'true' states, even if augmented
    % state model used.
    %
    % logForwardProbMat stores the forward probabilities matrix (log
    % domain). Should be preallocated - initialised with all -Inf (i.e. the
    % log of a probability of 0) - but not essential. If not, then supply
    % an empty matrix.
    %
    % logTransitionMat the log of the transition matrix of the MJP
    % conditioned on a set on jump times. If it is 3D, it must have the
    % same number of pages as the other probability structures.
    %
    % logInitialStateDist is the initial state distribution; if multiple
    % processes/particles, then these must correspond to the same indices
    % as the transition matrices. For a fixed initial state, one entry will
    % have (log) probability 0, the rest -Inf.
    %
    % 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. In each column there is a binary
    % indicator of inclusion or exclusion for that particle at that time
    % step. Where this is false, logLikelihood will be -Inf, and so will
    % logForwardProbMat at the end (rather than NaN, since this is useful
    % for identifying when 0/0 occurs).
    %
    % ModelSpec must contain the generatorMatNonzeroRowInds,
    % generatorMatNonzeroColInds and nNonzerosAugGeneratorMat parameters if
    % the augmented state version is used. See
    % setupaugmentedstateindices().
    %
    % 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.
    % Thus the preallocated forward matrix supplied must be of the correct
    % size for this dimension.
    %
    % The convention for forward probability structures is:
    % dim1 = time
    % dim2 = state
    % dim3 = particle
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 13.10.2014 - copied from old version, forwardfilteringshortinline.
    % Altered parameter names and format to new version.
    % 14.10.2014 - changed logLikelihood to have n 'true' states columns in
    % the augmented state model (uses less memory and should be faster).
    % 17.10.2014 - removed requirement for logLikelihood to be negative -
    % not necessary with contious observations.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargin < 5,
        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(logInitialStateDist, 1);
    else
        if size(logInitialStateDist, 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
    % Better for memory management if this is preallocated.
    if isempty(logForwardProbMat),
        logForwardProbMat = -inf([maxNTransitions, nStates, nParticles]);
    else
        if size(logForwardProbMat, 1) ~= maxNTransitions,
            error('Size of logForwardProbMat must be [maxNTransitions, nStates, nParticles]!')
        end
        if size(logForwardProbMat, 2) ~= nStates,
            error('Size of logForwardProbMat must be [maxNTransitions, nStates, nParticles]!')
        end
        if size(logForwardProbMat, 3) ~= nParticles,
            error('Size of logForwardProbMat must be [maxNTransitions, nStates, 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) ~= size(logTransitionMat, 2),
        error('First two dimensions of logTransitionMat must be a square matrix!')
    end
    if size(logInitialStateDist, 1) ~= size(logTransitionMat, 1),
        error('n states of logInitialStateDist must match that of logTransitionMat!')
    end
    if size(logInitialStateDist, 2) ~= nParticles,
        error('n columns of logInitialStateDist must match the number of pages (particles) of other data structures!')
    end
    if any(any(any(logTransitionMat > 0))),
        error('logTransitionMat contains log probabilities greater than 0!')
    end
    if any(any(any(logInitialStateDist > 0))),
        error('logInitialStateDist contains log probabilities greater than 0!')
    end
%     if any(any(any(logLikelihood > 0))),
%         error('logLikelihood contains log probabilities greater than 0!')
%     end
    if nargin < 6,
        trajectoryMatBinInds = true([maxNTransitions, nParticles]);
    else
        if size(trajectoryMatBinInds, 2) ~= nParticles,
            error('Size of trajectoryMatBinInds must match that of other structures!')
        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
        if ~all(trajectoryMatBinInds(1, :)),
            error('No particles must be excluded from EVERY time point!')
        end
    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 ~ModelSpec.useAugStateModel,
        % Initialse with initial state distribution and first observation.
        % Note we do not restrict particles, assuming that they will all be
        % involved in the first step at least.
        logForwardProbMat(1, :, :) = permute(logInitialStateDist, [3, 1, 2]) + logLikelihood(1, :, :);
        % Best to do this transformation outside of the time loop.
        logTransitionMat = permute(logTransitionMat, [4, 1, 3, 2]);
        for t = 2:maxNTransitions,
            % Multiply source states' probabilities with transition matrix
            % and marginalise.
            forwardInter =  bsxfun(@plus, logTransitionMat(:, :, trajectoryMatBinInds(t, :), :), logForwardProbMat(t - 1, :, trajectoryMatBinInds(t, :)));
            % Take account of observations via the likelihood.
            logForwardProbMat(t, :, trajectoryMatBinInds(t, :)) = logLikelihood(t, :, trajectoryMatBinInds(t, :)) + permute(sumlogprobmat(forwardInter, 2), [1, 4, 3, 2]);
        end
    else
        % Initialse with initial state distribution and first observation.
        % Note we do not restrict particles, assuming that they will all be
        % involved in the first step at least.
        logForwardProbMat(1, :, :) = permute(logInitialStateDist, [3, 1, 2]) + logLikelihood(1, ModelSpec.augStateComponents(1:nStates, 1), :);
        % 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 '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]), floor((1:(1/ModelSpec.nNonzerosAugGeneratorMat(stateDim)):(nParticles + 1 - 1/ModelSpec.nNonzerosAugGeneratorMat(stateDim)))')];
        accumarryIndexingMat = [repmat(ModelSpec.generatorMatNonzeroColInds{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.
        generatorMatNonzeroRowInds = ModelSpec.generatorMatNonzeroRowInds{stateDim};
        % The 'destination state' associated with each linear-indexed
        % non-zero entry in the augmented transition matrix.
        generatorMatNonzeroColInds = 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;
        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]);
        nParticlesIncludedEachTime = sum(trajectoryMatBinInds, 2);
        for t = 2:maxNTransitions,
            repeatedParticleIndsIncluded = 1:(nParticlesIncludedEachTime(t) * ModelSpec.nNonzerosAugGeneratorMat(stateDim));
            generatorMatLinearIndsExtendedThisTime = generatorMatLinearIndsExtended(repeatedParticleIndsIncluded, :);
            
            forwardInterMatThisTime = forwardInterMat(:, :, trajectoryMatBinInds(t, :));
            forwardInter1 = logForwardProbMat(t - 1, generatorMatNonzeroRowInds, trajectoryMatBinInds(t, :));
            % Multiply source states' probabilities.
            logAugTransitionMat = logTransitionMat(:, :, trajectoryMatBinInds(t, :));
            if nStates > 1,
                forwardInterMatThisTime(generatorMatLinearIndsExtendedThisTime) = logAugTransitionMat(generatorMatLinearIndsExtendedThisTime) + forwardInter1(:);
            else
                forwardInterMatThisTime(generatorMatLinearIndsExtendedThisTime) = permute(logAugTransitionMat(generatorMatLinearIndsExtendedThisTime), [3, 2, 1]) + forwardInter1(:);
            end
            % Prepare to sum over source states.
            maxLogProbEachCol = max(forwardInterMatThisTime, [], 1);
            maxLogProbEachCol(maxLogProbEachCol == -Inf) = 0;
            maxLogProbEachColNonzerosOnly = maxLogProbEachCol(:, generatorMatNonzeroColInds, :);
            if nStates > 1,
                forwardInterMatThisTime(generatorMatLinearIndsExtendedThisTime) = forwardInterMatThisTime(generatorMatLinearIndsExtendedThisTime) - maxLogProbEachColNonzerosOnly(:);
            else
                forwardInterMatThisTime(generatorMatLinearIndsExtendedThisTime) = permute(forwardInterMatThisTime(generatorMatLinearIndsExtendedThisTime), [3, 2, 1]) - maxLogProbEachColNonzerosOnly(:);
            end
            % Sum over source states.
            if nStates > 1,
                forwardInter2 = log(accumarray(accumarryIndexingMat(repeatedParticleIndsIncluded, :), exp(forwardInterMatThisTime(generatorMatLinearIndsExtendedThisTime)), [nStates, nParticlesIncludedEachTime(t)])) + permute(maxLogProbEachCol, [2, 3, 1]);
            else
                forwardInter2 = log(accumarray(accumarryIndexingMat(repeatedParticleIndsIncluded, :), permute(exp(forwardInterMatThisTime(generatorMatLinearIndsExtendedThisTime)), [3, 2, 1]), [nStates, nParticlesIncludedEachTime(t)])) + permute(maxLogProbEachCol, [2, 3, 1]);
            end
            % Take account of observations via the likelihood.
            logForwardProbMat(t, :, trajectoryMatBinInds(t, :)) = logLikelihood(t, ModelSpec.augStateComponents(1:nStates, 1), trajectoryMatBinInds(t, :)) + permute(forwardInter2, [3, 1, 2]);
        end
    end
end

