function logLikelihood = computediscretetimespiketrainslikelihood( discreteSpikeTrainsMat, logPoissonLikelihoodFactor, logSpikeRate, logNonspikeRate, model, updateStepInterval, trajectoryMatBinInds, maxNTransitions, nParticles )
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % logLikelihood is a matrix containing the likelihood vectors: each row
    % is indexed by state, giving the likelihood function on each state at
    % each time, which indexes rows.
    %
    % trajectoryMatBinInds - matrix of logicals indicating how
    % matrices of likelihoods are to be indexed by streams with potentially
    % different numbers of intervals.
    %
    % Note: With discrete time data, we insist that all particles sample on
    % the same discrete time scale. That is, if we discretised time so that
    % time t is the start of time bin i, then time bin i of the likelihood
    % corresponds to the discrete interval from t for every stream.
    % Continuous time data could be different: in that case the likelihood
    % may be evaluated over different time intervals for different streams.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 14.10.2014 - copied from old version,
    % computepoissoneventtimelikelihood.m. Changed names and formats of
    % parameters. We now return the likelihood without 'augmenting' the
    % state space - leave number of columns as number of 'true' states.
    % This is faster and uses less memory.
    % 15.10.2014 - removed requirement for initial jump times to all be the
    % same - we may compute likelihoods for multiple particles from
    % different lower time bound.
    % 23.10.2014 - added the possibility of empty spike trains.
    % 27.11.2014 - moved the removal of the last element of histc output to
    % outside the loop over spike trains, for greater efficiency.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargin < 7,
        error('Insufficient arguments!')
    end
    if nargin < 8,
        maxNTransitions = size(trajectoryMatBinInds, 1);
    else
        if size(trajectoryMatBinInds, 1) ~= maxNTransitions,
            error('trajectoryMatBinInds must have as many rows as the biggest column vector of time points in jumpTimes!')
        end
    end
%     if size(discreteSpikeTrainsMat, 1) ~= maxNTransitions,
%         error('discreteSpikeTrainsMat must have as many rows as the biggest column vector of time points in jumpTimes!')
%     end
%     if size(logPoissonLikelihoodFactor, 1) ~= maxNTransitions,
%         error('logPoissonLikelihoodFactor must have as many rows as the biggest column vector of time points in jumpTimes!')
%     end
    if nargin < 9,
        nParticles = size(logSpikeRate, 3);
    else
        if size(logSpikeRate, 3) ~= nParticles,
            error('Size of logSpikeRate does not match nParticles!')
        end
    end
    if size(trajectoryMatBinInds, 2) ~= nParticles,
        error('trajectoryMatBinInds must have the same number of column as the cell array jumpTimes has cells!')
    end
    
    timeBinsVec = (1:maxNTransitions)';
    timeBinsVec = timeBinsVec(trajectoryMatBinInds(:, 1));
    
    % We assume all columns of trajectoryMatBinInds are identical.
    spikingPart = bsxfun(@times, permute(discreteSpikeTrainsMat(timeBinsVec, :), [1, 4, 3, 2]), permute(logSpikeRate, [4, 1, 3, 2]));
    if model == 1,
        nonSpikesMat = ~discreteSpikeTrainsMat(timeBinsVec, :);
        nonSpikingPart = bsxfun(@times, permute(nonSpikesMat, [1, 4, 3, 2]), permute(logNonspikeRate, [4, 1, 3, 2]));
        logLikelihood = sum(spikingPart + nonSpikingPart, 4);
    elseif model == 2,
        commonFactor = -updateStepInterval * permute(exp(sumlogprobmat(logSpikeRate, 2)), [2, 1, 3]);
        logLikelihood = bsxfun(@plus, bsxfun(@plus, sum(spikingPart, 4), commonFactor), logPoissonLikelihoodFactor(timeBinsVec));
    end

end