function viterbiPath = computeviterbipath( EstModelParams, ModelSpec, SpikeTrainData, PositionData, DataSpec, logLikelihood, timePointsVec, nTimePoints, logInitialStateDist, logTransitionMat, nStates, startTime, endTime )
    % Given a fitted model, compute the Viterbi path of hidden states on a
    % discrete time grid (i.e. we turn the model into a HMM).
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % logInitialStateDist - to be supplied separately from ModelParams, and
    % appropriate to the interval of data (e.g. the invariant distribution
    % of the process).
    %
    % logTransitionMat - to be supplied separately from ModelParams.
    % Required only for the exact approach to computing the smoothed
    % posteriors by imposing a time grid on the data and considering the
    % model as an HMM. Must be computed (from the generator matrix) for the
    % time discretisation, which must also be specified (time interval must
    % be known).
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 27.01.2015 - copied from old version.
    % 28.01.2015 - coded. Added functionality for 'sub-Viterbi' paths.
    % 05.02.2015 - removed functionality for 'sub-Viterbi' paths - our approach was in error. See computelistviterbipaths.m for correct approach.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargin < 9,
        error('Insufficient arguments!')
    end
    if nargin < 11,
        nStates = length(logInitialStateDist);
    end
    if nargin < 12,
        startTime = DataSpec.startTime;
    end
    if nargin < 13,
        endTime = DataSpec.endTime;
    end
    if isempty(timePointsVec),
        % Setup time discretisation.
        if isempty(nTimePoints),
            error('At least one of timePointsVec or nTimePoints must be nonempty!')
        end
        timeStep = (endTime - startTime) / nTimePoints;
        timePointsVec = (startTime:timeStep:endTime)';
        % We can't have a discretisation point actually AT the end time,
        % because we exclude data at that point by convention.
        timePointsVec = timePointsVec(1:(end - 1));

        % Must now also compute the transition matrix, knowing the time
        % step.
        logTransitionMat = log(matrixexponentialmulti(timeStep .* EstModelParams.generatorMat));
    else
        if size(timePointsVec, 2) ~= 1,
            error('timePointsVec must be a column vector!');
        end
    end
    nTimePoints = length(timePointsVec);
    if isempty(logLikelihood),
        haveLikelihood = false;
    else
        if ~isequal(size(logLikelihood), [nTimePoints, nStates]),
            error('Size of logLikelihood does not match nTimePoints and nStates!')
        end
        haveLikelihood = true;
    end

    % We will not use augmented states.
    ModelSpec.useAugStateModel = false;

    % Compute likelihood at each time, if not supplied.
    if ~haveLikelihood,
        timeBinInds = true([nTimePoints, 1]);
        interjumpIntervalsArray = {[timePointsVec(2:end) - timePointsVec(1:(end - 1)); endTime - timePointsVec(end)]};
        % Compute likelihood.
        logLikelihood = computelikelihood(ModelSpec, SpikeTrainData, PositionData, EstModelParams, false, {timePointsVec}, timeBinInds, startTime, endTime, ModelParams.nTrueStatesThisBlock, interjumpIntervalsArray, nStateTimePoints, 1);
    end
    
    % Compute the optimal path probabilities up to each time t, for each
    % state at time t (forward recursions).
    logOptimalPathProb = -inf([nTimePoints, nStates]);
    % The back pointer at time step 1 is left undefined.
    backpointers = nan([nTimePoints, nStates]);
    % Initialisation.
    logOptimalPathProb(1, :) = logInitialStateDist' + logLikelihood(1, :);
    for t = 2:nTimePoints,
        % Store backpointers at each time.
        [logOptimalPathProb(t, :), backpointers(t, :)] = max(bsxfun(@plus, logOptimalPathProb(t - 1, :)', logTransitionMat), [], 1);
        logOptimalPathProb(t, :) = logOptimalPathProb(t, :) + logLikelihood(t, :);
    end

    % Reconstruct optimal path.
    viterbiPath = nan([nTimePoints, 1]);
    [~, rankedEndStateInds] = sort(logOptimalPathProb(end, :), 2, 'descend');
    viterbiPath(end) = rankedEndStateInds(1);
    for t = (nTimePoints - 1):(-1):1,
        viterbiPath(t) = backpointers(t + 1, viterbiPath(t + 1));
    end
end
