function [ AlgorithmParams, ModelParams, Diagnostics ] = evolveparticlesdiscretetime( AlgorithmParams, ModelParams, PriorParams, ModelSpec, DataSpec, AlgorithmSpec, SpikeTrainData, PositionData )
    % Perform the particle algorithm forward recursions, accumulating
    % information from observations and updating weights accordingly, then
    % resampling to correct for degeneracy.
    %
    % This version has been adapted back for the discrete time algorithm.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % ModelSpec is a struct with the following fields:
    % - maxNStates - upper bound on state dimension represented by
    % any particle.
    % - useAugStateModel - logical scalar.
    % - augStateComponents - 2-D matrix with 2 columns used for
    % converting augmented state indices to the 'true' states and state
    % space dimension corresponding to each.
    % - stateSpaceTransformation - column vector of the augmented state
    % space dimensions corresponding to the row index 'true' state space
    % dimensions.
    % - spikeTrainModelIndicator
    % - positionModelIndicator
    % - transMatNonzeroRowInds
    % - transMatNonzeroColInds
    % - nNonzerosAugGeneratorMat - these last three needed if augmented state
    % space used. See setupaugmentedstateindices().
    %
    % AlgorithmSpec is a struct with the following fields:
    % - nParticles - total number of particles being used by particle
    % filter.
    % - nTimeBins - number of discrete time bins for particle filtering
    % algorithm.
    % - updateStepTimes - vector of times (in seconds) correspondings to
    % the SMC steps, at which we update weights and (possibly) resample
    % particles.
    % - essThreshold - numeric scalar, for triggering resampling.
    % - usePositiveDiscrimination - logical scalar.
    % - nGibbsSweepsPerResampleMove - nonnegetive integer.
    % - maxBlockSize - positive integer. Should reflect memory constraints.
    % - allowSmallSubpopulationsToGoExtinct - logical scalar.
    % - smallSubpopulationBuffer
    % - positiveDiscriminationThreshold - these last three needed if
    % positive discrimination is used.
    % - timeDiscretisation - 1 for full fixed-width discretisation, 2 for
    % continuous time variable interval widths.
    % - updateStepInterval - time duration, in seconds, of the fixed
    % width bins for the discrete time algorithm.
    % - startTime - time in seconds at the start of this analysis
    % interval.
    % - useResampling - logical scalar.
    % - alwaysResample - logical scalar.
    % - isDatasetNoncontiguous - logical scalar.
    % - datasetStartTimes - vector of positive reals: times (in seconds) at
    % which the data should not be considered to be a contiguous step of
    % the Markov chain model. Can be empty and only required if
    % isDatasetNoncontiguous is true.
    %
    % dominatingProcessJumpTimes - cell array indexed by particles,
    % with cells down the first column (so this is a 'column' array). Each
    % cell (corresponding to a block) is itself a column cell array,
    % containing column vectors, one for each particle in each block. These
    % column vectors are of time points in seconds indicating when the
    % discrete time intervals start for that particle.
    %
    % Note on the form of the parameter cell arrays: these should be
    % 'vector' arrays, i.e. having only one non-singleton dimension,
    % containing the cells that correspond to different blocks. This
    % non-singleton dimension should be the same dimension as the
    % corresponding particles dimension in the parameter structures
    % themselves (the cells). This is to ensure cell2mat and mat2cell
    % maintain the same structure consistently. For example, the
    % AlgorithmParams.logParticleWeights cell array should be a column cell array, and
    % posCovs should have its 4th dimension be the non-singleton.
    %
    % Note: rather than looping over all state dimensions, within the time
    % loop, to update the weights for different subpopulations separately,
    % we could convert all parameter structures to be the same size, padded
    % with blank values (zeros, say) where necessary. This would save an
    % extra loop, but I think would necessitate more computations overall
    % because of all the computations on blank values that don't contribute
    % anything.
    %
    % Note: weights are normalised at the end, after setting new values,
    % although this is not necessary for the algorithm. We do it so we can
    % compute the ESS immediately after resampling, for diagnostic
    % purposes.
    %
    % Note: for the single step forward algorithm functions, when the
    % augmented state space is used, there are some necessary
    % (preallocated) intermediate structures. These depend on the number of
    % particles and states, hence are different for different
    % subpopulations. They need to be updated after every resampling step
    % because the number of particles in each subpopulation could change.
    % Hence these structures are stored in cell arrays, one cell each
    % subpopulation, and are updated using a call to a nested function.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 10.10.2014 - copied from old version.
    % 21.10.2014 - got function working for basic setup.
    % 23.10.2014 - added multiple Gibbs sweeps for sampling MJP
    % trajectories.
    % 05.02.2015 - added the ability to pause the algorithm at some time in the future, and then to resume it.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargin < 7,
        error('Insufficient arguments!')
    end
    if ~isfield(ModelSpec, 'modelIndicator'),
        error('modelIndicator must be a field of ModelSpec!')
    end
    if ModelSpec.modelIndicator ~= 1,
        error('modelIndicator suggests we are to use the continuous time model!')
    end
    if ~isfield(AlgorithmParams, 'nTrueStatesEachBlock'),
        error('nTrueStatesEachBlock must be a field of AlgorithmParams!')
    end
    if ~isfield(AlgorithmParams, 'nParticlesEachBlock'),
        error('nParticlesEachBlock must be a field of AlgorithmParams!')
    end
    if any((rem(AlgorithmParams.nTrueStatesEachBlock, 1) ~= 0) | AlgorithmParams.nTrueStatesEachBlock < 1),
        error('AlgorithmParams.nTrueStatesEachBlock must be a vector of positive integers (state dimensions carried by blocks of particles)!')
    end
    if any(AlgorithmParams.nTrueStatesEachBlock > ModelSpec.maxNStates),
        error('State dimensions in AlgorithmParams.nTrueStatesEachBlock must be bounded above by ModelSpec.maxNStates!')
    end
    if ~isfield(AlgorithmParams, 'nBlocks'),
        error('nBlocks must be a field of AlgorithmParams!')
    end
    if ~isfield(AlgorithmParams, 'nBlocksEachSubpop'),
        error('nBlocksEachSubpop must be a field of AlgorithmParams!')
    end
    if sum(AlgorithmParams.nBlocksEachSubpop) ~= AlgorithmParams.nBlocks,
        error('nBlocksEachSubpop must sum to nBlocks!')
    end
    if ~isfield(AlgorithmParams, 'logParticleWeights'),
        error('logParticleWeights must be a field of AlgorithmParams!')
    end
    if ~isfield(AlgorithmParams, 'logNormalisedParticleWeights'),
        error('logNormalisedParticleWeights must be a field of AlgorithmParams!')
    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 ~isfield(ModelSpec, 'maxNStates'),
        error('maxNStates must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'spikeTrainModelIndicator'),
        error('spikeTrainModelIndicator must be a field of ModelSpec!')
    end
    if ~isfield(ModelSpec, 'positionModelIndicator'),
        error('positionModelIndicator must be a field of ModelSpec!')
    end
    if length(ModelParams) ~= AlgorithmParams.nBlocks,
        error('Parameter struct array ModelParams must be indexed by block and agree in length with AlgorithmParams.nBlocks!')
    end
    if ModelSpec.spikeTrainModelIndicator > 0,
        if ~isfield(ModelParams, 'logSpikeRate'),
            error('logSpikeRate must be a field of ModelParams!')
        end
    end
    if ModelSpec.spikeTrainModelIndicator == 1,
        if ~isfield(ModelParams, 'logNonspikeRate'),
            error('logNonspikeRate must be a field of ModelParams!')
        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, 'logPosDistribution'),
%             error('logPosDistribution must be a field of ModelParams!')
%         end
        if ~isfield(PositionData, 'posTimeVec'),
            error('posTimeVec must be a field of PositionData!')
        end
        if ~isfield(PositionData, 'linearisedPositionTrajectory'),
            error('linearisedPositionTrajectory must be a field of PositionData!')
        end
        if ~isfield(PositionData, 'gridCentreCoordsTrajectory'),
            error('gridCentreCoordsTrajectory must be a field of PositionData!')
        end
        if ~isfield(PositionData, 'constrainedDistanceMat'),
            error('constrainedDistanceMat must be a field of PositionData!')
        end
        if ~isfield(PositionData, 'positionIndToGridCentreCoordsMap'),
            error('positionIndToGridCentreCoordsMap must be a field of PositionData!')
        end
        if ~isfield(PositionData, 'nValidDiscretePositions'),
            error('nValidDiscretePositions must be a field of PositionData!')
        end
        if ~isfield(PositionData, 'validDiscretePositionsVector'),
            error('validDiscretePositionsVector must be a field of PositionData!')
        end
    end
    if ~isfield(ModelSpec, 'useAugStateModel'),
        error('useAugStateModel must be a field of ModelSpec!')
    end
    if ModelSpec.useAugStateModel,
        if ~isfield(ModelSpec, 'generatorMatNonzeroRowInds'),
            error('transMatNonzeroRowInds must be a field of ModelSpec!')
        end
        if ~isfield(ModelSpec, 'generatorMatNonzeroColInds'),
            error('transMatNonzeroColInds must be a field of ModelSpec!')
        end
        if ~isfield(ModelSpec, 'nNonzerosAugGeneratorMat'),
            error('nNonzerosAugGeneratorMat must be a field of ModelSpec!')
        end
        if ~isfield(ModelSpec, 'augStateComponents'),
            error('augStateComponents must be a field of ModelSpec!')
        end
        if ~isfield(ModelSpec, 'stateSpaceTransformation'),
            error('stateSpaceTransformation must be a field of ModelSpec!')
        end
    end
%     if ~isfield(AlgorithmSpec, 'nResampleLeadTimeBins'),
%         error('nResampleLeadTimeBins must be a field of AlgorithmSpec!')
%     end
    if ~isfield(AlgorithmSpec, 'debugMode'),
        error('debugMode must be a field of AlgorithmSpec!')
    end
    if ~isfield(AlgorithmSpec, 'nParticles'),
        error('nParticles must be a field of AlgorithmSpec!')
    end
    if ~isfield(AlgorithmSpec, 'nUpdateSteps'),
        error('nUpdateSteps must be a field of AlgorithmSpec!')
    end
    if ~isfield(AlgorithmSpec, 'updateStepTimes'),
        error('updateStepTimes must be a field of AlgorithmSpec!')
    end
    if ~isfield(AlgorithmSpec, 'essThreshold'),
        error('essThreshold must be a field of AlgorithmSpec!')
    end
    if ~isfield(AlgorithmSpec, 'usePositiveDiscrimination'),
        error('usePositiveDiscrimination must be a field of AlgorithmSpec!')
    end
    if ~isfield(AlgorithmSpec, 'nGibbsSweepsPerResampleMove'),
        error('nGibbsSweepsPerResampleMove must be a field of AlgorithmSpec!')
    end
    if AlgorithmSpec.usePositiveDiscrimination,
        if ~isfield(AlgorithmSpec, 'allowSmallSubpopulationsToGoExtinct'),
            error('allowSmallSubpopulationsToGoExtinct must be a field of AlgorithmSpec!')
        end
        if ~isfield(AlgorithmSpec, 'smallSubpopulationBuffer'),
            error('smallSubpopulationBuffer must be a field of AlgorithmSpec!')
        end
        if ~isfield(AlgorithmSpec, 'positiveDiscriminationThreshold'),
            error('positiveDiscriminationThreshold must be a field of AlgorithmSpec!')
        end
    end
    if (ModelSpec.spikeTrainModelIndicator > 0) && (size(SpikeTrainData.discreteSpikeTrainsMat, 1) ~= AlgorithmSpec.nUpdateSteps),
        error('We do not have enough spike observations on this discretisation level to perform nTimeBins time bins of particle filtering!')
    end
%     if (ModelSpec.positionModelIndicator == 1) && (size(PositionData.linearisedPositionTrajectoryDiscreteTime, 1) ~= AlgorithmSpec.nUpdateSteps),
%         error('We do not have enough position observations on this discretisation level to perform nTimeBins time bins of particle filtering!')
%     end

    % We seem to get a problem related to real number precision unless we
    % do this: in position likelihood function, reducing position data
    % structures to time interval of interest, small numerical differences
    % between update step times and position observation times cause us to
    % select the wrong positions.
    PositionData.posTimeVec = AlgorithmSpec.updateStepTimes(1:(end - 1));

    % Initialise everything (if we are not resuming from a pause).
    if ~isfield(AlgorithmParams, 'startUpdateStep') || (AlgorithmParams.startUpdateStep <= 2),
        startUpdateStep = 1;

        Diagnostics.effectiveSampleSize = nan([AlgorithmSpec.nUpdateSteps, 2]);
        Diagnostics.resampleEachStepBinInd = false([AlgorithmSpec.nUpdateSteps, 1]);
        Diagnostics.estLogStateDimensionPosterior = -inf([AlgorithmSpec.nUpdateSteps, ModelSpec.maxNStates]);
        if AlgorithmSpec.inferenceAlgorithmIndicator == 1,
            if ModelSpec.spikeTrainModelIndicator > 0,
                blocksThisSubpopBinInds = AlgorithmParams.nTrueStatesEachBlock == ModelSpec.maxNStates;
                Diagnostics.spikeRateEstimatesSmc = -inf([AlgorithmSpec.nUpdateSteps + 1, ModelSpec.maxNStates, SpikeTrainData.nSpikeTrains]);
                sample = permute([ModelParams(blocksThisSubpopBinInds).logSpikeRate], [3, 1, 2]);
                sample = sumlogprobmat(bsxfun(@plus, [AlgorithmParams.logParticleWeights{blocksThisSubpopBinInds}], sample), 1) - sumlogprobmat([AlgorithmParams.logParticleWeights{blocksThisSubpopBinInds}], 1);
                Diagnostics.spikeRateEstimatesSmc(1, :, :) = sample;

                Diagnostics.spikeRateSquaredEstimatesSmc = -inf([AlgorithmSpec.nUpdateSteps + 1, ModelSpec.maxNStates, SpikeTrainData.nSpikeTrains]);
                sample = permute([2 .* ModelParams(blocksThisSubpopBinInds).logSpikeRate], [3, 1, 2]);
                sample = sumlogprobmat(bsxfun(@plus, [AlgorithmParams.logParticleWeights{blocksThisSubpopBinInds}], sample), 1) - sumlogprobmat([AlgorithmParams.logParticleWeights{blocksThisSubpopBinInds}], 1);
                Diagnostics.spikeRateSquaredEstimatesSmc(1, :, :) = sample;
            end
        end

        % Initialise Markov jump process trajectory structures.
        stateTrajectory = cell([1, AlgorithmParams.nBlocks]);
        jumpTimes = cell([AlgorithmParams.nBlocks, 1]);
        interjumpIntervals = cell([AlgorithmParams.nBlocks, 1]);
        trajectoryMatBinInds = cell([1, AlgorithmParams.nBlocks]);
        maxNTransitions = zeros([AlgorithmParams.nBlocks, 1]);
        for iBlock = 1:AlgorithmParams.nBlocks,
            [ stateTrajectory{iBlock}, trajectoryMatBinInds{iBlock}, jumpTimes{iBlock}, interjumpIntervals{iBlock}, maxNTransitions(iBlock) ] = sampleinitialstate( ModelParams(iBlock), AlgorithmSpec, ModelParams(iBlock).nParticlesThisBlock );
        end
        % Marginal likelihood for no data: probability one.
        logPartialMarginalLikelihoodsCell = cellfun(@(nParticles) zeros([1, nParticles]), num2cell(AlgorithmParams.nParticlesEachBlock), 'UniformOutput', false);

        % Initialise particle blocks.
        [ ModelParams, AlgorithmParams, stateTrajectory, trajectoryMatBinInds, jumpTimes, interjumpIntervals, maxNTransitions, logPartialMarginalLikelihoodsCell ] = resizeblocks( ModelParams, AlgorithmParams, ModelSpec, AlgorithmSpec, stateTrajectory, trajectoryMatBinInds, jumpTimes, interjumpIntervals, maxNTransitions, 1, logPartialMarginalLikelihoodsCell );

        logEndForwardProbMat = cell([1, AlgorithmParams.nBlocks]);

%         % Set the first time of a forced resample.
%         % This establishes how many states we need to sample before the
%         % first iteration (see loop below).
%         forcedResampleTimeBin = startUpdateStep + AlgorithmSpec.nResampleLeadTimeBins - 1;

    % Restart from a pause.
    else
        startUpdateStep = AlgorithmParams.startUpdateStep;
        Diagnostics = AlgorithmParams.Diagnostics;
        stateTrajectory = AlgorithmParams.stateTrajectory;
        jumpTimes = AlgorithmParams.jumpTimes;
        interjumpIntervals = AlgorithmParams.interjumpIntervals;
        trajectoryMatBinInds = AlgorithmParams.trajectoryMatBinInds;
        maxNTransitions = AlgorithmParams.maxNTransitions;
        logPartialMarginalLikelihoodsCell = AlgorithmParams.logPartialMarginalLikelihoodsCell;
        logEndForwardProbMat = AlgorithmParams.logEndForwardProbMat;
%         forcedResampleTimeBin = AlgorithmParams.forcedResampleTimeBin;
        AlgorithmParams = rmfield(AlgorithmParams, {'Diagnostics', 'stateTrajectory', 'jumpTimes', 'interjumpIntervals', 'trajectoryMatBinInds', 'maxNTransitions', 'logPartialMarginalLikelihoodsCell'});
    end

    for iUpdate = startUpdateStep:AlgorithmSpec.nUpdateSteps,
        if AlgorithmSpec.inferenceAlgorithmIndicator == 2,
            % For Gibbs sampling only.
%             if ModelSpec.spikeTrainModelIndicator > 0,
%                 Diagnostics.spikeRateEstimatesGibbs = -inf([AlgorithmSpec.nGibbsSweepsPerResampleMove + 1, ModelSpec.maxNStates, SpikeTrainData.nSpikeTrains]);
%                 Diagnostics.spikeRateEstimatesGibbs(1, :, :) = permute([ModelParams(:).logSpikeRate], [3, 1, 2]);
%             end
            Diagnostics.SampleParams = samplemodelparametersfromprior( PriorParams, ModelSpec, DataSpec, AlgorithmSpec.nGibbsSweepsPerResampleMove );
            Diagnostics.SampleParams.logInitialStateDist(:, 1) = ModelParams.logInitialStateDist;
            Diagnostics.SampleParams.logTransitionMat(:, :, 1) = ModelParams.logTransitionMat;
            if ModelSpec.spikeTrainModelIndicator > 0,
                Diagnostics.SampleParams.logSpikeRate(:, :, 1) = ModelParams.logSpikeRate;
            end
            if ModelSpec.positionModelIndicator > 0,
                Diagnostics.SampleParams.driftRate(:, 1) = ModelParams.driftRate;
                Diagnostics.SampleParams.asympMeanPosInd(:, 1) = ModelParams.asympMeanPosInd;
                if ModelSpec.positionJitterModelIndicator == 3,
                    Diagnostics.SampleParams.jitterCov(:, :, :, 1) = ModelParams.jitterCov;
                    Diagnostics.SampleParams.posCov(:, :, :, 1) = ModelParams.posCov;
                    Diagnostics.SampleParams.stationaryPosCov(:, :, :, 1) = ModelParams.stationaryPosCov;
                else
                    Diagnostics.SampleParams.jitterCov(:, :, 1) = ModelParams.jitterCov;
                    Diagnostics.SampleParams.posCov(:, :, 1) = ModelParams.posCov;
                    Diagnostics.SampleParams.stationaryPosCov(:, :, 1) = ModelParams.stationaryPosCov;
                end
            end
            move();
            return;
        end
        % Check if we are due to pause the algorithm.
        if isfield(AlgorithmSpec, 'pauseTime') && ~isempty(AlgorithmSpec.pauseTime),
            theTime = clock;
            if (theTime(3) >= AlgorithmSpec.pauseTime(1)) && ((theTime(4) * 60 + theTime(5)) >= (AlgorithmSpec.pauseTime(2) * 60 + AlgorithmSpec.pauseTime(3))),
%             if iUpdate >= (AlgorithmSpec.nUpdateSteps / 10),
                % Pause: store all working algorithm parameters and samples as they are.
                AlgorithmParams.startUpdateStep = iUpdate;
                AlgorithmParams.Diagnostics = Diagnostics;
                AlgorithmParams.stateTrajectory = stateTrajectory;
                AlgorithmParams.jumpTimes = jumpTimes;
                AlgorithmParams.interjumpIntervals = interjumpIntervals;
                AlgorithmParams.trajectoryMatBinInds = trajectoryMatBinInds;
                AlgorithmParams.maxNTransitions = maxNTransitions;
                AlgorithmParams.logPartialMarginalLikelihoodsCell = logPartialMarginalLikelihoodsCell;
                AlgorithmParams.logEndForwardProbMat = logEndForwardProbMat;
%                 AlgorithmParams.forcedResampleTimeBin = forcedResampleTimeBin;
                return;
            end
        end
        reportloopprogress( iUpdate, AlgorithmSpec.nUpdateSteps, 1/100 );

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % If we resampled last step, must now obtain marginal likelihood
        % for data up to previous step, the denominator of the weight
        % update factor. This is essential for correct continuation of the
        % weight-updating forward recursions: the previously calculated
        % marginal likelihoods are now based on different (and redundant)
        % parameter values. We do this by inserting into the dominating
        % process an dummy jump time at the last update step time. Then
        % after computing the forward probabilities, retain the one we need
        % and discard the dummy jump time.

        logMarginalLikelihoodsCell = cell([1, AlgorithmParams.nBlocks]);

        if AlgorithmSpec.debugMode,
            sprintf('%d ', AlgorithmParams.nParticlesEachBlock)
        end

        for iBlock = 1:AlgorithmParams.nBlocks,

            % Our objective in this loop is to sample new state process
            % trajectories and compute marginal likelihoods as required by
            % our algorithm specification.

            % If we moved parameter estimates last update step we need to
            % compute the marginal likelihood at that step as well as at
            % this step to calculate the weight update factor.
            if (iUpdate > 1) && Diagnostics.resampleEachStepBinInd(iUpdate - 1),
%             if iUpdate > 1,
                [ logEndForwardProbMat{iBlock}, jumpTimes{iBlock}, interjumpIntervals{iBlock}, trajectoryMatBinInds{iBlock}, maxNTransitions(iBlock), logMarginalLikelihoodsCell{iBlock}, logPartialMarginalLikelihoodsCell{iBlock} ] = samplehmmstatetrajectory( ModelParams(iBlock), ModelSpec, SpikeTrainData, PositionData, [], jumpTimes{iBlock}, trajectoryMatBinInds{iBlock}, AlgorithmSpec.updateStepTimes(iUpdate + 1), AlgorithmSpec.updateStepTimes(iUpdate), maxNTransitions(iBlock), ModelParams(iBlock).nParticlesThisBlock );
            elseif iUpdate == 1,
                [ logEndForwardProbMat{iBlock}, jumpTimes{iBlock}, interjumpIntervals{iBlock}, trajectoryMatBinInds{iBlock}, maxNTransitions(iBlock), logMarginalLikelihoodsCell{iBlock} ] = samplehmmstatetrajectory( ModelParams(iBlock), ModelSpec, SpikeTrainData, PositionData, [], jumpTimes{iBlock}, trajectoryMatBinInds{iBlock}, AlgorithmSpec.updateStepTimes(iUpdate + 1), AlgorithmSpec.updateStepTimes(iUpdate + 1), maxNTransitions(iBlock), ModelParams(iBlock).nParticlesThisBlock );
            % Just need marginal likelihood.
            else
                [ logEndForwardProbMat{iBlock}, jumpTimes{iBlock}, interjumpIntervals{iBlock}, trajectoryMatBinInds{iBlock}, maxNTransitions(iBlock), logMarginalLikelihoodsCell{iBlock} ] = samplehmmstatetrajectory( ModelParams(iBlock), ModelSpec, SpikeTrainData, PositionData, logEndForwardProbMat{iBlock}, jumpTimes{iBlock}, trajectoryMatBinInds{iBlock}, AlgorithmSpec.updateStepTimes(iUpdate + 1), AlgorithmSpec.updateStepTimes(iUpdate), maxNTransitions(iBlock), ModelParams(iBlock).nParticlesThisBlock );
            end

        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Weight updates.
        % This should be done by subpopulation; i.e. all blocks
        % corresponding to maxNStates 1, followed by all blocks
        % corresponding to maxNStates 2, and so on.
        for iSubpop = 1:ModelSpec.maxNStates,
            if AlgorithmParams.nParticlesEachSubpop(iSubpop) == 0,
                continue;
            end
            blocksThisSubpopBinInds = AlgorithmParams.nTrueStatesEachBlock == iSubpop;
            AlgorithmParams.nParticlesEachBlockThisSubpop = [ModelParams(blocksThisSubpopBinInds).nParticlesThisBlock];
            
            logMarginalLikelihoods = cell2mat(logMarginalLikelihoodsCell(blocksThisSubpopBinInds));

            logPartialMarginalLikelihoods = cell2mat(logPartialMarginalLikelihoodsCell(blocksThisSubpopBinInds));

            logWeightUpdateFactor = logMarginalLikelihoods - logPartialMarginalLikelihoods;

            % 'Correct' for 0 probabilities?
            logWeightUpdateFactor((logMarginalLikelihoods == -Inf) & (logPartialMarginalLikelihoods == -Inf)) = -Inf;

            % Update weights.
            AlgorithmParams.logParticleWeights(blocksThisSubpopBinInds) = mat2cell(cell2mat(AlgorithmParams.logNormalisedParticleWeights(blocksThisSubpopBinInds)) + permute(logWeightUpdateFactor, [2, 1]), AlgorithmParams.nParticlesEachBlockThisSubpop, 1);
        end
        
        % Normalise weights. Note this must be done over all
        % subpopulations, hence outside the above loop.
        AlgorithmParams.logNormalisedParticleWeights = mat2cell(normaliselogdistributionsmatrix(cell2mat(AlgorithmParams.logParticleWeights), 1), AlgorithmParams.nParticlesEachBlock, 1);
        
        % Retain marginal likelihoods at this time step for the weight
        % update at the next time step (in case we don't resample and have
        % to compute it again anyway).
        logPartialMarginalLikelihoodsCell = logMarginalLikelihoodsCell;
        
        % Update the estimate for the posterior over number of states. This
        % is used in the resampling routine.
        for iSubpop = 1:ModelSpec.maxNStates,
            if AlgorithmParams.nParticlesEachSubpop(iSubpop) == 0,
                continue;
            end
            blocksThisSubpopBinInds = AlgorithmParams.nTrueStatesEachBlock == iSubpop;
            Diagnostics.estLogStateDimensionPosterior(iUpdate, iSubpop) = sumlogprobmat(cell2mat(AlgorithmParams.logNormalisedParticleWeights(blocksThisSubpopBinInds)), 1);
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Resample?

        if AlgorithmSpec.useResampling && AlgorithmSpec.alwaysResample,
            Diagnostics.resampleEachStepBinInd(iUpdate) = true;
            Diagnostics.effectiveSampleSize(iUpdate, 1) = computeess( cell2mat(AlgorithmParams.logNormalisedParticleWeights(:)), AlgorithmSpec.nParticles );
        elseif AlgorithmSpec.useResampling,
            % The ESS is calculated using the normalised particle weights.
            [ Diagnostics.resampleEachStepBinInd(iUpdate), Diagnostics.effectiveSampleSize(iUpdate, 1) ] = resamplecheck( cell2mat(AlgorithmParams.logNormalisedParticleWeights(:)), AlgorithmSpec.essThreshold, AlgorithmSpec.nParticles );
        end
        
        if AlgorithmSpec.debugMode,
%             [iUpdate, Diagnostics.effectiveSampleSize(iUpdate, 1)]
            sprintf('ESS: %0.1f', Diagnostics.effectiveSampleSize(iUpdate, 1))
            
            Diagnostics.estLogStateDimensionPosterior(iUpdate, :)
%             [~, kappaMax] = max(Diagnostics.estLogStateDimensionPosterior(iUpdate, :));
%             sprintf('kappa hat: %d', kappaMax)

            %[ModelParams.nParticlesThisBlock]
        end
        
        % Yes. Resample and 'move' parameters.
        if Diagnostics.resampleEachStepBinInd(iUpdate),
            if AlgorithmSpec.debugMode,
                sprintf('Resampling. . .')
            end
                        
            resample;
            move;
            
            % ESS after resampling.
            Diagnostics.effectiveSampleSize(iUpdate, 2) = computeess( cell2mat(AlgorithmParams.logNormalisedParticleWeights(:)), AlgorithmSpec.nParticles );

        else
            
            Diagnostics.effectiveSampleSize(iUpdate, 2) = Diagnostics.effectiveSampleSize(iUpdate, 1);
        end
        
%         % Whether or not we have sampled new parameter values, we should
%         % look at resizing blocks, because in the next update step there
%         % will be more dominating process arrivals, incurring a greater
%         % memory demand on each block.
%         [ ModelParams, AlgorithmParams, stateTrajectory, trajectoryMatBinInds, jumpTimes, interjumpIntervals, maxNTransitions, logPartialMarginalLikelihoodsCell ] = resizeblocks( ModelParams, AlgorithmParams, ModelSpec, AlgorithmSpec, stateTrajectory, trajectoryMatBinInds, jumpTimes, interjumpIntervals, maxNTransitions, iUpdate, logPartialMarginalLikelihoodsCell );

        if AlgorithmSpec.inferenceAlgorithmIndicator == 1,
            try
            if ModelSpec.spikeTrainModelIndicator > 0,
                blocksThisSubpopBinInds = AlgorithmParams.nTrueStatesEachBlock == ModelSpec.maxNStates;
%                 nParticlesEachBlockThisSubpop = [ModelParams(blocksThisSubpopBinInds).nParticlesThisBlock];
                nParticlesEachBlock = [ModelParams.nParticlesThisBlock]';
                sample = permute([ModelParams(blocksThisSubpopBinInds & (nParticlesEachBlock > 0)).logSpikeRate], [3, 1, 2]);
                sample = sumlogprobmat(bsxfun(@plus, [AlgorithmParams.logParticleWeights{blocksThisSubpopBinInds & (nParticlesEachBlock > 0)}], sample), 1) - sumlogprobmat([AlgorithmParams.logParticleWeights{blocksThisSubpopBinInds & (nParticlesEachBlock > 0)}], 1);
                Diagnostics.spikeRateEstimatesSmc(iUpdate + 1, :, :) = sample;

                sample = permute([2 .* ModelParams(blocksThisSubpopBinInds & (nParticlesEachBlock > 0)).logSpikeRate], [3, 1, 2]);
                sample = sumlogprobmat(bsxfun(@plus, [AlgorithmParams.logParticleWeights{blocksThisSubpopBinInds & (nParticlesEachBlock > 0)}], sample), 1) - sumlogprobmat([AlgorithmParams.logParticleWeights{blocksThisSubpopBinInds & (nParticlesEachBlock > 0)}], 1);
                Diagnostics.spikeRateSquaredEstimatesSmc(iUpdate + 1, :, :) = sample;
            end
            end
        end
        
    end
    
    % Finalise output structures.
    AlgorithmParams.logNormalisedParticleWeights = mat2cell(normaliselogdistributionsmatrix(cell2mat(AlgorithmParams.logParticleWeights), 1), AlgorithmParams.nParticlesEachBlock, 1);
   
    
    
    
    
    function resample()
    
        % Adjust the weights for the resampling algorithm according to
        % positive discrimination?
        % Also returns discriminatedSubpopBinInds: a row vctor that
        % indicates which subpopulations are being positively
        % discriminated.

        % Normalise weights.
        % (Should already be done.)
%         AlgorithmParams.logNormalisedParticleWeights = mat2cell(normaliselogdistributionsmatrix(cell2mat(AlgorithmParams.logParticleWeights), 1), AlgorithmParams.nParticlesEachBlock, 1);

        if AlgorithmSpec.usePositiveDiscrimination,
            [samplingWeights, discriminatedSubpopBinInds] = positivediscrimination( AlgorithmParams, Diagnostics.estLogStateDimensionPosterior(iUpdate, :), AlgorithmSpec, ModelSpec.maxNStates );
        else
            samplingWeights = exp(cell2mat(AlgorithmParams.logNormalisedParticleWeights(:)));
            discriminatedSubpopBinInds = false([1, ModelSpec.maxNStates]);
        end
        % Expand out to 'blocks' rather than just subpopulations.
        discriminatedBlockBinInds = discriminatedSubpopBinInds(AlgorithmParams.nTrueStatesEachBlock);
        nTrueStatesEachDiscriminatedBlock = AlgorithmParams.nTrueStatesEachBlock(discriminatedBlockBinInds);
        
        % Returns a cell array of the count of each particle (corresponding
        % to those in logParticleWeights) to reappear in the population.
        % Particles are either culled (nCopiesToRetainEachParticle(i) ==
        % 0) or retained (nCopiesToRetainEachParticle(i) > 0).
        nCopiesToRetainEachParticle = residualresampling(samplingWeights, AlgorithmParams.nParticlesEachBlock);

        % Reset weights.
        if any(~discriminatedBlockBinInds),
            AlgorithmParams.logParticleWeights(~discriminatedBlockBinInds) = cellfun(@(weightsBlock) zeros(size(weightsBlock)), AlgorithmParams.logParticleWeights(~discriminatedBlockBinInds), 'UniformOutput', false);
        end
%         try
        if any(discriminatedBlockBinInds),
            AlgorithmParams.logParticleWeights(discriminatedBlockBinInds) = cellfun(@(weightsBlock, iLogStateDimPosterior) (iLogStateDimPosterior - log(AlgorithmSpec.positiveDiscriminationThreshold)) * ones(size(weightsBlock)), AlgorithmParams.logParticleWeights(discriminatedBlockBinInds), num2cell((Diagnostics.estLogStateDimensionPosterior(iUpdate, nTrueStatesEachDiscriminatedBlock))'), 'UniformOutput', false);
        end
%         catch err
%             nTrueStatesEachDiscriminatedBlock
%             error(err.message)
%             error(err.stack)
%         end

        % Normalise weights. Note this must be done over all
        % subpopulations.
        AlgorithmParams.logNormalisedParticleWeights = mat2cell(normaliselogdistributionsmatrix(cell2mat(AlgorithmParams.logParticleWeights), 1), AlgorithmParams.nParticlesEachBlock, 1);

        % Copy out structures indexed by particle the requisite number of
        % times.
        firstCopyBinInds = cell([AlgorithmParams.nBlocks, 1]);
        for jBlock = 1:AlgorithmParams.nBlocks,
            nParticlesThisBlock = sum(nCopiesToRetainEachParticle{jBlock});
            repeatedCopiedParticleInds = repeatindicesbycounts( nCopiesToRetainEachParticle{jBlock} );
            stateTrajectory{jBlock} = stateTrajectory{jBlock}(:, repeatedCopiedParticleInds);
            trajectoryMatBinInds{jBlock} = trajectoryMatBinInds{jBlock}(:, repeatedCopiedParticleInds);
            jumpTimes{jBlock} = jumpTimes{jBlock}(repeatedCopiedParticleInds);
            interjumpIntervals{jBlock} = interjumpIntervals{jBlock}(repeatedCopiedParticleInds);
            logPartialMarginalLikelihoodsCell{jBlock} = logPartialMarginalLikelihoodsCell{jBlock}(repeatedCopiedParticleInds);
            if size(logPartialMarginalLikelihoodsCell{jBlock}, 1) > 1,
                logPartialMarginalLikelihoodsCell{jBlock} = logPartialMarginalLikelihoodsCell{jBlock}';
            end
            AlgorithmParams.logParticleWeights{jBlock} = AlgorithmParams.logParticleWeights{jBlock}(repeatedCopiedParticleInds);
            AlgorithmParams.logNormalisedParticleWeights{jBlock} = AlgorithmParams.logNormalisedParticleWeights{jBlock}(repeatedCopiedParticleInds);
            AlgorithmParams.nParticlesEachBlock(jBlock) = nParticlesThisBlock;
            ModelParams(jBlock).nParticlesThisBlock = nParticlesThisBlock;
            ModelParams(jBlock).logInitialStateDist = ModelParams(jBlock).logInitialStateDist(:, repeatedCopiedParticleInds);
            ModelParams(jBlock).logTransitionMat = ModelParams(jBlock).logTransitionMat(:, :, repeatedCopiedParticleInds);
            if ModelSpec.spikeTrainModelIndicator > 0,
                ModelParams(jBlock).logSpikeRate = ModelParams(jBlock).logSpikeRate(:, :, repeatedCopiedParticleInds);
                if ModelSpec.spikeTrainModelIndicator == 1,
                    ModelParams(jBlock).logNonspikeRate = ModelParams(jBlock).logNonspikeRate(:, :, repeatedCopiedParticleInds);
                end
            end
            if ModelSpec.positionModelIndicator > 0,
                ModelParams(jBlock).asympMeanPosInd = ModelParams(jBlock).asympMeanPosInd(:, repeatedCopiedParticleInds);
                if (ModelSpec.positionJitterModelIndicator == 0) || (ModelSpec.positionJitterModelIndicator == 1) || (ModelSpec.positionJitterModelIndicator == 2),
                    ModelParams(jBlock).jitterCov = ModelParams(jBlock).jitterCov(:, :, repeatedCopiedParticleInds);
                    ModelParams(jBlock).posCov = ModelParams(jBlock).posCov(:, :, repeatedCopiedParticleInds);
%                     ModelParams(jBlock).stationaryPosCov = ModelParams(jBlock).stationaryPosCov(:, :, repeatedCopiedParticleInds);
                elseif ModelSpec.positionJitterModelIndicator == 3,
                    ModelParams(jBlock).jitterCov = ModelParams(jBlock).jitterCov(:, :, :, repeatedCopiedParticleInds);
                    ModelParams(jBlock).posCov = ModelParams(jBlock).posCov(:, :, :, repeatedCopiedParticleInds);
%                     ModelParams(jBlock).stationaryPosCov = ModelParams(jBlock).stationaryPosCov(:, :, :, repeatedCopiedParticleInds);
                end
                if ModelSpec.positionModelIndicator == 1,
                    ModelParams(jBlock).posObsVar = ModelParams(jBlock).posObsVar(repeatedCopiedParticleInds);
                end
                ModelParams(jBlock).driftRate = ModelParams(jBlock).driftRate(:, repeatedCopiedParticleInds);
%                 ModelParams(jBlock).logPosDistribution = ModelParams(jBlock).logPosDistribution(:, :, repeatedCopiedParticleInds);
            end

            % We will also need logical indicators of which particles need
            % to have MJP paths sampled (i.e. only 'new copies').
            % Indicator of which particles are being retained in any
            % quantity.
            particlesResampledBinInds = nCopiesToRetainEachParticle{jBlock} > 0;
            % Number of particles to populate this block.
%             nParticlesResampled = sum(nCopiesToRetainEachParticle{jBlock});
            nParticlesResampled = nParticlesThisBlock;
            % This tells us the number of indices to count along to get to
            % the first in each set of copies from one 'ancestor'.
            nIndsToFirstOfSubsequentCopies = cumsum(nCopiesToRetainEachParticle{jBlock}(particlesResampledBinInds)) + 1;
            % And these are the indices of the first in each set of copies.
            firstCopyInds = [1; nIndsToFirstOfSubsequentCopies(1:(end - 1))];
            % In logical indexing form.
            firstCopyBinInds{jBlock} = false([nParticlesResampled, 1]);
            firstCopyBinInds{jBlock}(firstCopyInds) = true;

        end
        
        % Resize blocks based on the new subpopulation sizes and memory
        % constraints.
        % We double up trajectoryMatBinInds here to use it in place of
        % stateTrajectory - which has not been sampled yet. This doesn't
        % affect anything.
        [ ModelParams, AlgorithmParams, ~, trajectoryMatBinInds, jumpTimes, interjumpIntervals, maxNTransitions, logPartialMarginalLikelihoodsCell, ~ ] = resizeblocks( ModelParams, AlgorithmParams, ModelSpec, AlgorithmSpec, trajectoryMatBinInds, trajectoryMatBinInds, jumpTimes, interjumpIntervals, maxNTransitions, iUpdate, logPartialMarginalLikelihoodsCell, firstCopyBinInds );
        
    end

    function move()

        % Parameters 'move' kernel: sample from full conditional given MJP
        % state trajectory.
        for jSweep = 1:AlgorithmSpec.nGibbsSweepsPerResampleMove,
            if AlgorithmSpec.inferenceAlgorithmIndicator == 2,
                reportloopprogress( jSweep, AlgorithmSpec.nGibbsSweepsPerResampleMove, 1/1000 );
            end
            for jBlock = 1:AlgorithmParams.nBlocks,
               
                % Sample new paths.
                % Potential for multiple sweeps of MJP sampler.
                for kSweep = 1:AlgorithmSpec.nSweepsForMjpPost,
                    [ stateTrajectory{jBlock}, jumpTimes{jBlock}, interjumpIntervals{jBlock}, trajectoryMatBinInds{jBlock}, maxNTransitions(jBlock) ] = samplehmmstatetrajectory( ModelParams(jBlock), ModelSpec, SpikeTrainData, PositionData, [], jumpTimes{jBlock}, trajectoryMatBinInds{jBlock}, AlgorithmSpec.updateStepTimes(iUpdate + 1), AlgorithmSpec.updateStepTimes(iUpdate + 1), maxNTransitions(jBlock), ModelParams(jBlock).nParticlesThisBlock );
                end

                % Now we have a full set of sampled state paths. The next
                % step of the 'move' kernel is, conditioned on those
                % sampled paths, sample parameter values.

                % Compute statistics needed for full conditionals and
                % sample.
                stateEachObsTime = [];
                currentTimeParticleInds = [];
                stateParticleIndsEachObsTime = [];
                currentPosCountEachState = [];
                nObsEachState = [];
                spikeCountsMat = [];
                firstArrivals = [];
                if ModelSpec.useAugStateModel,
                    if ModelSpec.spikeTrainModelIndicator > 0,
                        if ModelSpec.positionModelIndicator > 0,
                            [ stateOccupancyDurations, transitionCountsMat, firstArrivals, spikeCountsMat, stateEachObsTime, currentTimeParticleInds, stateParticleIndsEachObsTime, currentPosCountEachState, previousPosCountEachState, nObsEachState ] = computesamplingstatistics( ModelSpec, SpikeTrainData, PositionData, stateTrajectory{jBlock}, jumpTimes{jBlock}, trajectoryMatBinInds{jBlock}, interjumpIntervals{jBlock}, ModelParams(jBlock).nTrueStatesThisBlock, AlgorithmSpec.updateStepTimes(iUpdate + 1), maxNTransitions(jBlock), ModelParams(jBlock).nParticlesThisBlock );
                        else
                            [ stateOccupancyDurations, transitionCountsMat, firstArrivals, spikeCountsMat ] = computesamplingstatistics( ModelSpec, SpikeTrainData, PositionData, stateTrajectory{jBlock}, jumpTimes{jBlock}, trajectoryMatBinInds{jBlock}, interjumpIntervals{jBlock}, ModelParams(jBlock).nTrueStatesThisBlock, AlgorithmSpec.updateStepTimes(iUpdate + 1), maxNTransitions(jBlock), ModelParams(jBlock).nParticlesThisBlock );
                        end
                    else
                        if ModelSpec.positionModelIndicator > 0,
                            [ stateOccupancyDurations, transitionCountsMat, firstArrivals, stateEachObsTime, currentTimeParticleInds, stateParticleIndsEachObsTime, currentPosCountEachState, previousPosCountEachState, nObsEachState ] = computesamplingstatistics( ModelSpec, SpikeTrainData, PositionData, stateTrajectory{jBlock}, jumpTimes{jBlock}, trajectoryMatBinInds{jBlock}, interjumpIntervals{jBlock}, ModelParams(jBlock).nTrueStatesThisBlock, AlgorithmSpec.updateStepTimes(iUpdate + 1), maxNTransitions(jBlock), ModelParams(jBlock).nParticlesThisBlock );
                        else
                            [ stateOccupancyDurations, transitionCountsMat, firstArrivals ] = computesamplingstatistics( ModelSpec, SpikeTrainData, PositionData, stateTrajectory{jBlock}, jumpTimes{jBlock}, trajectoryMatBinInds{jBlock}, interjumpIntervals{jBlock}, ModelParams(jBlock).nTrueStatesThisBlock, AlgorithmSpec.updateStepTimes(iUpdate + 1), maxNTransitions(jBlock), ModelParams(jBlock).nParticlesThisBlock );
                        end
                    end
                else
                    if ModelSpec.spikeTrainModelIndicator > 0,
                        if ModelSpec.positionModelIndicator > 0,
                            [ stateOccupancyDurations, transitionCountsMat, spikeCountsMat, stateEachObsTime, currentTimeParticleInds, stateParticleIndsEachObsTime, currentPosCountEachState, previousPosCountEachState, nObsEachState ] = computesamplingstatistics( ModelSpec, SpikeTrainData, PositionData, stateTrajectory{jBlock}, jumpTimes{jBlock}, trajectoryMatBinInds{jBlock}, interjumpIntervals{jBlock}, ModelParams(jBlock).nTrueStatesThisBlock, AlgorithmSpec.updateStepTimes(iUpdate + 1), maxNTransitions(jBlock), ModelParams(jBlock).nParticlesThisBlock );
                        else
                            [ stateOccupancyDurations, transitionCountsMat, spikeCountsMat ] = computesamplingstatistics( ModelSpec, SpikeTrainData, PositionData, stateTrajectory{jBlock}, jumpTimes{jBlock}, trajectoryMatBinInds{jBlock}, interjumpIntervals{jBlock}, ModelParams(jBlock).nTrueStatesThisBlock, AlgorithmSpec.updateStepTimes(iUpdate + 1), maxNTransitions(jBlock), ModelParams(jBlock).nParticlesThisBlock );
                        end
                    else
                        if ModelSpec.positionModelIndicator > 0,
                            [ stateOccupancyDurations, transitionCountsMat, stateEachObsTime, currentTimeParticleInds, stateParticleIndsEachObsTime, currentPosCountEachState, previousPosCountEachState, nObsEachState ] = computesamplingstatistics( ModelSpec, SpikeTrainData, PositionData, stateTrajectory{jBlock}, jumpTimes{jBlock}, trajectoryMatBinInds{jBlock}, interjumpIntervals{jBlock}, ModelParams(jBlock).nTrueStatesThisBlock, AlgorithmSpec.updateStepTimes(iUpdate + 1), maxNTransitions(jBlock), ModelParams(jBlock).nParticlesThisBlock );
%                         else
%                             [ stateOccupancyDurations, transitionCountsMat ] = computesamplingstatistics( ModelSpec, SpikeTrainData, PositionData, stateTrajectory{jBlock}, jumpTimes{jBlock}, trajectoryMatBinInds{jBlock}, interjumpIntervals{jBlock}, ModelParams(jBlock).nTrueStatesThisBlock, AlgorithmSpec.updateStepTimes(iUpdate), maxNTransitions(jBlock), ModelParams(jBlock).nParticlesThisBlock );
                        end
                    end
                end
                
                % Sample.
                ModelParams(jBlock) = samplemodelparametersfromfullconditional( ModelParams(jBlock), PriorParams, ModelSpec, DataSpec, PositionData, AlgorithmSpec.updateStepTimes(iUpdate + 1), stateTrajectory{jBlock}, maxNTransitions(jBlock), stateOccupancyDurations, transitionCountsMat, firstArrivals, spikeCountsMat, stateEachObsTime, currentPosCountEachState, previousPosCountEachState, nObsEachState, currentTimeParticleInds, stateParticleIndsEachObsTime );
                
                if any(cellfun(@(x) any(x(2:end) <= 0), interjumpIntervals{jBlock})),
                    error('Bad interjumpIntervals!')
                end
                                
            end

            if AlgorithmSpec.inferenceAlgorithmIndicator == 2,
                % For Gibbs sampling only.
%                 if ModelSpec.spikeTrainModelIndicator > 0,
%                     Diagnostics.spikeRateEstimatesGibbs(jSweep + 1, :, :) = permute([ModelParams(:).logSpikeRate], [3, 1, 2]);
%                 end
                Diagnostics.SampleParams.logInitialStateDist(:, jSweep) = ModelParams.logInitialStateDist;
                Diagnostics.SampleParams.logTransitionMat(:, :, jSweep) = ModelParams.logTransitionMat;
                if ModelSpec.spikeTrainModelIndicator > 0,
                    Diagnostics.SampleParams.logSpikeRate(:, :, jSweep) = ModelParams.logSpikeRate;
                end
                if ModelSpec.positionModelIndicator > 0,
                    Diagnostics.SampleParams.driftRate(:, jSweep) = ModelParams.driftRate;
                    Diagnostics.SampleParams.asympMeanPosInd(:, jSweep) = ModelParams.asympMeanPosInd;
                    if ModelSpec.positionJitterModelIndicator == 3,
                        Diagnostics.SampleParams.jitterCov(:, :, :, jSweep) = ModelParams.jitterCov;
                        Diagnostics.SampleParams.posCov(:, :, :, jSweep) = ModelParams.posCov;
                        Diagnostics.SampleParams.stationaryPosCov(:, :, :, jSweep) = ModelParams.stationaryPosCov;
                    else
                        Diagnostics.SampleParams.jitterCov(:, :, jSweep) = ModelParams.jitterCov;
                        Diagnostics.SampleParams.posCov(:, :, jSweep) = ModelParams.posCov;
                        Diagnostics.SampleParams.stationaryPosCov(:, :, jSweep) = ModelParams.stationaryPosCov;
                    end
                end
            end
        end
        
    end
end
