function varargout = decodestate( ModelParams, ModelSpec, AlgorithmSpec, SpikeTrainData, PositionData, DataSpec, method, startTime, endTime, timePointsVec, nTimePoints, logInitialStateDist, logTransitionMat, nSamples )
    % This function controls the decoding (estimation) of state via the
    % smoothed posteriors (computed using one of two different methods) or
    % the Viterbi algorithm.
    %
    % We can decode state on any (sub)interval of data. For this we
    % typically assume the state process is stationary, and use the
    % invariant distribution over state to initialise the process at the
    % start of the interval (i.e. for logInitialStateDist).
    %
    % We assume that all samples have the same number of states (state
    % space dimension). Thus ModelParams is a scalar struct, not a struct
    % array.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % startTime, endTime - of the data interval.
    %
    % 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).
    %
    % method - 1 for exact smoothed posteriors on a discrete time grid (HMM
    % approach); 2 for estimated smoothed posteriors on a similarly
    % specified time grid; 3 for Viterbi algorithm on a time grid. For 1
    % and 2 MAP estimates are available as output as well as the posterior
    % distributions.
    %
    % ModelParams - is the struct of parameter point estimates, but in the
    % case of method 2 this can be a sample of parameter values, i.e. one
    % value for each of many particles.
    %
    % decodedStateTrajectory - an output matrix with two columns; first is
    % the sequence of time points, second is the estimate (MAP or Viterbi)
    % at that time.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 28.01.2015 - created.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargout > 3,
        error('Incorrect number of outputs!')
    end
    if (method == 3) && (nargout > 1),
        error('Incorrect number of outputs!')
    end
    if nargin < 7,
        error('Insufficient arguments!')
    end
    if nargin < 8,
        startTime = DataSpec.startTime;
    end
    if nargin < 9,
        endTime = DataSpec.endTime;
    end
    
    if isempty(timePointsVec),
        % Setup time discretisation.
        if (nargin < 11) || isempty(nTimePoints),
            error('At least one of timePointsVec or nTimePoints must be supplied and 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));
        nTimePoints = length(timePointsVec);
    else
        if size(timePointsVec, 2) ~= 1,
            error('timePointsVec must be a column vector!');
        end
        nTimePoints = length(timePointsVec);
        timeStep = (endTime - startTime) / nTimePoints;
    end

    % Must now also compute the transition matrix, knowing the time
    % step.
    if (nargin < 13) || isempty(logTransitionMat),
        logTransitionMat = log(matrixexponentialmulti(timeStep .* ModelParams.generatorMat));
    end
    
    if size(logTransitionMat, 3) == 1,
        havePointEstimates = true;
    else
        havePointEstimates = false;
    end
    if ~havePointEstimates && (method ~= 2),
        error('For methods 1 and 3 we require ModelParams to be point estimates!')
    end
    % For now, when using point estimates we do not use augmented states.
    if havePointEstimates && ModelSpec.useAugStateModel,
        ModelSpec.useAugStateModel = false;
    end
    if ModelSpec.useAugStateModel,
        if size(ModelParams.generatorMat, 1) ~= ModelSpec.maxNAugStates,
            error('Size of generatorMat does not match maxNAugStates! Did you mean to NOT use the augmented state model?')
        end
    end

    if (nargin < 12) || isempty(logInitialStateDist),
        if havePointEstimates && (startTime > DataSpec.startTime),
            logInitialStateDist = computemjpstationarydistribution( ModelParams.generatorMat );
        elseif startTime == DataSpec.startTime,
            logInitialStateDist = ModelParams.logInitialStateDist;
        else
            error('Must supply an initial state distribution!')
        end
    end
    if havePointEstimates,
        ModelParams.nStatesThisBlock = size(logInitialStateDist, 1);
        ModelParams.nTrueStatesThisBlock = ModelParams.nStatesThisBlock;
    else
        if ~isfield(ModelParams, 'nTrueStatesThisBlock'),
            error('nTrueStatesThisBlock must be a field of ModelParams!')
        end
        if ~isfield(ModelParams, 'nStatesThisBlock'),
            error('nStatesThisBlock must be a field of ModelParams!')
        end
    end
    if ~isfield(ModelParams, 'nStatesThisBlock'),
        error('nStatesThisBlock must be a field of ModelParams!')
    end
%     if ~isfield(ModelParams, 'generatorMat'),
%         error('generatorMat must be a field of ModelParams!')
%     end
%     if ~isfield(ModelParams, 'jumpRate'),
%         error('jumpRate must be a field of ModelParams!')
%     end
%     if ~isfield(ModelParams, 'dominatingRate'),
%         error('dominatingRate must be a field of ModelParams!')
%     end
    if ~isfield(ModelParams, 'logTransitionMat'),
        error('logTransitionMat must be a field of ModelParams!')
    end
    if ~isfield(ModelParams, 'logInitialStateDist'),
        error('logInitialStateDist must be a field of ModelParams!')
    end
    if ModelSpec.spikeTrainModelIndicator > 0,
        if ~isfield(ModelParams, 'logSpikeRate'),
            error('logSpikeRate must be a field of ModelParams!')
        end
        if ~isfield(SpikeTrainData, 'spikeTimesArray'),
            error('spikeTimesArray must be a field of SpikeTrainData!')
        end
        if ~isfield(SpikeTrainData, 'nSpikeTrains'),
            error('nSpikeTrains must be a field of SpikeTrainData!')
        end
    end
    if ModelSpec.positionModelIndicator > 0,
        if ~isfield(ModelParams, 'asympMeanPosInd'),
            error('asympMeanPosInd must be a field of ModelParams!')
        end
        if ~isfield(ModelParams, 'posCov'),
            error('posCov must be a field of ModelParams!')
        end
        if ~isfield(ModelParams, 'driftRate'),
            error('driftRate must be a field of ModelParams!')
        end
    end
    
    % Exact smoothed posteriors (HMM approach).
    if method == 1,
        timeBinInds = true([nTimePoints, 1]);
        interjumpIntervalsArray = {[timePointsVec(2:end) - timePointsVec(1:(end - 1)); endTime - timePointsVec(end)]};
        % Compute likelihood.
        logLikelihood = computelikelihood(ModelSpec, SpikeTrainData, PositionData, ModelParams, false, {timePointsVec}, timeBinInds, startTime, endTime, ModelParams.nTrueStatesThisBlock, nTimePoints, 1, interjumpIntervalsArray);

        logForwardProbMat = forwardfiltering( ModelSpec, [], logLikelihood, logTransitionMat, logInitialStateDist, timeBinInds, nTimePoints, ModelParams.nTrueStatesThisBlock );
        logBackwardProbMat = backwardsmoothing( ModelSpec, [], logLikelihood, logTransitionMat, timeBinInds, nTimePoints, ModelParams.nTrueStatesThisBlock );

        logSmoothedPost = logForwardProbMat + logBackwardProbMat;
        logSmoothedPost = normaliselogdistributionsmatrix( logSmoothedPost, 2 );
        
        % MAP estimates.
        [~, decodedStateTrajectory] = max(logSmoothedPost, [], 2);
        decodedStateTrajectory = [timePointsVec, decodedStateTrajectory];
        
        % Marginal (data) likelihood, if we want it.
        if nargout == 3,
            logMarginalLikelihood = sumlogprobmat(logForwardProbMat(end, :), 2);
        end
        
    % Estimated smoothed posteriors.
    elseif method == 2,
        if nargin < 14,
            nSamples = 1000;
        end
        if nargout == 1,
            decodedStateTrajectory = estimatesmoothedposteriors( ModelParams, ModelSpec, AlgorithmSpec, SpikeTrainData, PositionData, DataSpec, timePointsVec, nTimePoints, nSamples, startTime, endTime );
        elseif nargout == 2,
            [decodedStateTrajectory, logSmoothedPost] = estimatesmoothedposteriors( ModelParams, ModelSpec, AlgorithmSpec, SpikeTrainData, PositionData, DataSpec, timePointsVec, nTimePoints, nSamples, startTime, endTime );
        elseif nargout == 3,
            % We optionally compute an estimate of the state transition
            % probabilities conditional on observations, used for decoding
            % position in another function.
            [logSmoothedPost, logObsCondStateTransitionProbs, decodedStateTrajectory] = estimatesmoothedposteriors( ModelParams, ModelSpec, AlgorithmSpec, SpikeTrainData, PositionData, DataSpec, timePointsVec, nTimePoints, nSamples, startTime, endTime );
        end
        
    % Viterbi path.
    elseif method == 3,
        timeBinInds = true([nTimePoints, 1]);
        interjumpIntervalsArray = {[timePointsVec(2:end) - timePointsVec(1:(end - 1)); endTime - timePointsVec(end)]};
        % Compute likelihood.
        logLikelihood = computelikelihood(ModelSpec, SpikeTrainData, PositionData, ModelParams, false, {timePointsVec}, timeBinInds, startTime, endTime, ModelParams.nTrueStatesThisBlock, nTimePoints, 1, interjumpIntervalsArray);
        
        decodedStateTrajectory = [timePointsVec, computeviterbipath( ModelParams, ModelSpec, SpikeTrainData, PositionData, DataSpec, logLikelihood, timePointsVec, nTimePoints, logInitialStateDist, logTransitionMat, ModelParams.nTrueStatesThisBlock, startTime, endTime )];
    end

    if nargout == 1,
        varargout = {decodedStateTrajectory};
    elseif nargout == 2,
        varargout = {decodedStateTrajectory, logSmoothedPost};
    elseif nargout == 3,
        if method == 1,
            varargout = {decodedStateTrajectory, logSmoothedPost, logMarginalLikelihood};
        elseif method == 2,
            varargout = {decodedStateTrajectory, logSmoothedPost, logObsCondStateTransitionProbs};
        end
    end
end