function [ varargout ] = resizeblocks( ModelParams, AlgorithmParams, ModelSpec, AlgorithmSpec, stateTrajectory, trajectoryMatBinInds, jumpTimes, interjumpIntervals, maxNTransitions, iUpdate, logPartialMarginalLikelihoodsCell, firstCopyBinInds )
    % With this function, we optimise the size of the 'blocks' (matrices in
    % cells of various parameter cell arrays) following a resampling step
    % of the particle algorithm. Blocks that are too big are split up, and
    % if some blocks have become sufficiently small that they can be
    % combined, we do so. We also clean up any empty blocks (removing
    % them).
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % logParticleWeights - cell array of column vectors, containing the
    % (log) weights of particles in each block, which correspond to
    % different state dimensions.
    %
    % AlgorithmParams.nTrueStatesEachBlock - column vector specifying the
    % state dimension carried by the corresponding blocks (cells of the
    % cell arrays).
    %
    % logCondObsLikelihoodPrev - cell array indexed by particles, like
    % logParticleWeights and the parameter structures, containing
    % conditional observations likelihood from the current step of the
    % particle filter weight updates algorithm. These are blocked and need
    % to be adapted in the same way as all other particle-indexed
    % structures. Contents are scalar numerics (log probabilities) for each
    % particle. This parameter is not necessary. If we do not wish to
    % 'adapt' and such structure, supply a [] in its place.
    %
    % logForwardProbPrev - similar cell array to logCondObsLikelihoodPrev,
    % containing forward probabilities computed at the current step of the
    % particle filter algorithm. This parameter is not necessary. If we do
    % not wish to 'adapt' and such structure, supply a [] in its place.
    %
    % particleDiscretisationTimePoints - 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.
    %
    % AlgorithmParams - struct with the following fields:
    % - maxBlockSize - maximum number of particles that we allow to
    % constitute any individual block.
    % - nParticles - total number of particles across all blocks.
    % - useAugmentedStateSpace - logical scalar.
    %
    % ModelSpec - struct with the following fields:
    % - maxStateDimension - upper bound on the various state dimensions we
    % consider.
    % - stateSpaceTransformation - column vector of the augmented state
    % space dimensions corresponding to the row index 'true' state space
    % dimensions.
    %
    % 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
    % logParticleWeights cell array should be a column cell array, and
    % posCovs should have its 4th dimension be the non-singleton.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 17.10.2014 - copied from old version, resizeblocks.m.
    % 27.10.2014 - added position model parameters.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if nargin < 10,
        error('Insufficient arguments!')
    end
    if nargin > 10,
        if (nargin == 11) && (nargout ~= 8),
            error('Incorrect number of outputs!')
        end
        havePartialMarginalLikelihoods = true;
    else
        havePartialMarginalLikelihoods = false;
        if nargout ~= 7,
            error('Incorrect number of outputs!')
        end
    end
    if nargin > 11,
        if nargout ~= 9,
            error('Incorrect number of outputs!')
        end
        haveFirstCopyBinInds = true;
    else
        haveFirstCopyBinInds = false;
    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 ~isfield(ModelSpec, 'useAugStateModel'),
        error('useAugStateModel must be a field of ModelSpec!')
    end
    if ~isfield(AlgorithmSpec, 'updateStepTimes'),
        error('updateStepTimes must be a field of AlgorithmSpec!')
    end
    if ~isfield(AlgorithmSpec, 'nParticles'),
        error('nParticles must be a field of AlgorithmSpec!')
    end
    if ~isfield(AlgorithmSpec, 'availableMemory'),
        error('availableMemory must be a field of AlgorithmSpec!')
    end
    if ~isequal(AlgorithmParams.nParticlesEachBlock, cellfun(@length, AlgorithmParams.logParticleWeights)),
        error('nParticlesEachBlock does not match size of particle weights vectors!')
    end
    if ~isfield(AlgorithmParams, 'nBlocks'),
        error('nBlocks must be a field of AlgorithmParams!')
    end
    if ~isfield(AlgorithmParams, 'nTrueStatesEachBlock'),
        error('nTrueStatesThisBlock must be a field of AlgorithmParams!')
    end
    if ~isfield(AlgorithmParams, 'nParticlesEachBlock'),
        error('nParticlesThisBlock must be a field of AlgorithmParams!')
    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, 'nStatesThisBlock'),
        error('nStatesThisBlock must be a field of ModelParams!')
    end
    if ~isfield(ModelParams, 'nTrueStatesThisBlock'),
        error('nTrueStatesThisBlock must be a field of ModelParams!')
    end
    if ~isfield(ModelParams, 'nParticlesThisBlock'),
        error('nParticlesThisBlock must be a field of ModelParams!')
    end
    if ModelSpec.modelIndicator == 2,
        if ~isfield(AlgorithmSpec, 'maxNJumpsLimit'),
            error('maxNJumpsLimit must be a field of AlgorithmSpec!')
        end
        if ~isfield(ModelParams, 'generatorMat'),
            error('logTransitionMat must be a field of ModelParams!')
        end
        if ~isfield(ModelParams, 'dominatingRate'),
            error('dominatingRate must be a field of ModelParams!')
        end
        if ~isfield(ModelParams, 'jumpRate'),
            error('jumpRate must be a field of ModelParams!')
        end
    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 sum(AlgorithmParams.nParticlesEachBlock) ~= AlgorithmSpec.nParticles,
        error('nParticlesThisBlock does not sum to nParticles!')
    end
    if length(ModelParams) ~= AlgorithmParams.nBlocks,
        error('Model parameter structure must be indexed by number of blocks!')
    end
    if length(stateTrajectory) ~= AlgorithmParams.nBlocks,
        error('stateTrajectory must be indexed by number of blocks!')
    end
    if length(trajectoryMatBinInds) ~= AlgorithmParams.nBlocks,
        error('trajectoryMatBinInds must be indexed by number of blocks!')
    end
    if length(jumpTimes) ~= AlgorithmParams.nBlocks,
        error('jumpTimes must be indexed by number of blocks!')
    end
    if length(interjumpIntervals) ~= AlgorithmParams.nBlocks,
        error('interjumpIntervals must be indexed by number of blocks!')
    end
    if length(maxNTransitions) ~= AlgorithmParams.nBlocks,
        error('maxNTransitions must be indexed by number of blocks!')
    end
    if havePartialMarginalLikelihoods,
        if length(logPartialMarginalLikelihoodsCell) ~= AlgorithmParams.nBlocks,
            error('logPartialMarginalLikelihoodsCell must be indexed by number of blocks!')
        end
    end
    if haveFirstCopyBinInds,
        if length(firstCopyBinInds) ~= AlgorithmParams.nBlocks,
            error('firstCopyBinInds must be indexed by number of blocks!')
        end
    end
    if size(AlgorithmParams.nTrueStatesEachBlock) ~= [AlgorithmParams.nBlocks, 1],
        error('AlgorithmParams.nTrueStatesEachBlock must be a column vector the same length as all parameter cell arrays (number of blocks)!')
    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 (rem(iUpdate, 1) ~= 0) || (iUpdate < 1),
        error('iUpdate must be a positive integer!')
    end
    if isfield(ModelParams, 'logNonspikeRate');
        haveNonspikeRate = true;
    else
        haveNonspikeRate = false;
    end
    
    if ModelSpec.useAugStateModel,
        nStates = ModelSpec.maxNAugStates;
    else
        nStates = ModelSpec.maxNStates;
    end
    
    % Initialise new structures with one block for each subpopulation.
    % How many subpops are there? Some may have gone extinct already.
    nSubpops = ModelSpec.maxNStates - AlgorithmParams.nTrueStatesEachBlock(1) + 1;

    TempAlgorithmParams.nParticlesEachBlock = zeros([ModelSpec.maxNStates, 1]);
    TempAlgorithmParams.nBlocksEachSubpop = ones([ModelSpec.maxNStates, 1]);
    TempAlgorithmParams.nTrueStatesEachBlock = zeros([ModelSpec.maxNStates, 1]);
    TempAlgorithmParams.logParticleWeights = cell([ModelSpec.maxNStates, 1]);
    TempAlgorithmParams.logNormalisedParticleWeights = cell([ModelSpec.maxNStates, 1]);
    
    tempStateTrajectory = cell([1, ModelSpec.maxNStates]);
    tempTrajectoryMatBinInds = cell([1, ModelSpec.maxNStates]);
    tempJumpTimes = cell([ModelSpec.maxNStates, 1]);
    tempInterjumpIntervals = cell([ModelSpec.maxNStates, 1]);
    tempMaxNTransitions = zeros([ModelSpec.maxNStates, 1]);
    if havePartialMarginalLikelihoods,
        tempLogPartialMarginalLikelihoodsCell = cell([1, ModelSpec.maxNStates]);
    end
    
    if ModelSpec.modelIndicator == 2,
        [TempModelParams(1:ModelSpec.maxNStates).generatorMat] = deal([]);
        [TempModelParams(1:ModelSpec.maxNStates).jumpRate] = deal([]);
        [TempModelParams(1:ModelSpec.maxNStates).dominatingRate] = deal([]);
    end
    [TempModelParams(1:ModelSpec.maxNStates).logInitialStateDist] = deal([]);
    [TempModelParams(1:ModelSpec.maxNStates).logTransitionMat] = deal([]);
    if ModelSpec.spikeTrainModelIndicator > 0,
        nSpikeTrains = size(ModelParams(1).logSpikeRate, 2);
        [TempModelParams(1:ModelSpec.maxNStates).logSpikeRate] = deal([]);
        if haveNonspikeRate,
            [TempModelParams(1:ModelSpec.maxNStates).logNonspikeRate] = deal([]);
        end
    end
    if ModelSpec.positionModelIndicator > 0,
        [TempModelParams(1:ModelSpec.maxNStates).driftRate] = deal([]);
%         [TempModelParams(1:ModelSpec.maxNStates).instDriftRate] = deal([]);
        [TempModelParams(1:ModelSpec.maxNStates).asympMeanPosInd] = deal([]);
        [TempModelParams(1:ModelSpec.maxNStates).posCov] = deal([]);
%         [TempModelParams(1:ModelSpec.maxNStates).stationaryPosCov] = deal([]);
%         if (ModelSpec.positionJitterModelIndicator == 0) || (ModelSpec.positionJitterModelIndicator == 1) || (ModelSpec.positionJitterModelIndicator == 2),
            [TempModelParams(1:ModelSpec.maxNStates).jitterCov] = deal([]);
%             [TempModelParams(1:ModelSpec.maxNStates).instJitterCov] = deal([]);
%         end
%         if ModelSpec.modelIndicator == 2,
%             [TempModelParams(1:ModelSpec.maxNStates).posObsVar] = deal([]);
%         end
    end
    TempModelParams = TempModelParams';
    
    if haveFirstCopyBinInds,
        tempFirstCopyBinInds = cell([ModelSpec.maxNStates, 1]);
    end

    % First we estimate the number of dominating process arrivals in the
    % next iteration of the algorithm. This affects the largest matrix
    % size, which determines the size of blocks.
    % We estimate the mean number of dominating arrivals. Depends on amount
    % of time covered in next update step, and the greatest dominating rate
    % of any particle.
    if ModelSpec.modelIndicator == 2,
        estMeanNJumps = (AlgorithmSpec.updateStepTimes(min([iUpdate + 1, end])) - AlgorithmSpec.updateStepTimes(1)) * mean(cat(1, ModelParams.dominatingRate)) + mean(maxNTransitions);
        % Make sure this is no more than the hard limit.
        estMeanNJumps = min([estMeanNJumps, AlgorithmSpec.maxNJumpsLimit]);
    elseif ModelSpec.modelIndicator == 1,
        estMeanNJumps = max(maxNTransitions);
    end
    % Calculate the mean size of the largest matrix in the next iteration.
    if ModelSpec.spikeTrainModelIndicator > 0,
        meanLargestMatrixSizePerParticle = max([estMeanNJumps * ModelSpec.maxNStates * nSpikeTrains, ...
                                            estMeanNJumps * nStates ^ 2]);
    else
        meanLargestMatrixSizePerParticle = estMeanNJumps * nStates ^ 2;
    end
    % Now we can say how many blocks we will need.
    maxBlockSize = ceil(AlgorithmSpec.availableMemory / meanLargestMatrixSizePerParticle / 8);
    
    % Need this to combine state trajectory matrices.
    maxMaxNJumps = max(maxNTransitions);

    % First we combine all blocks and sort them in ascending order of
    % dominating rate, within each subpopulation. The sorting is so that
    % after we partition into separate blocks, particles of similar
    % dominating rate are grouped together, which makes their largest
    % structures more efficient.
    for iSubpop = 1:ModelSpec.maxNStates,
        blocksThisSubpopBinInds = AlgorithmParams.nTrueStatesEachBlock == iSubpop;
        TempAlgorithmParams.nBlocksEachSubpop(iSubpop) = sum(blocksThisSubpopBinInds);
        if TempAlgorithmParams.nBlocksEachSubpop(iSubpop) == 0,
            continue;
        end

        TempAlgorithmParams.nParticlesEachBlock(iSubpop, 1) = sum(AlgorithmParams.nParticlesEachBlock(blocksThisSubpopBinInds));
        TempAlgorithmParams.nTrueStatesEachBlock(iSubpop, 1) = iSubpop;
        
        if ModelSpec.modelIndicator == 2,
            dominatingRate = cat(1, ModelParams(blocksThisSubpopBinInds).dominatingRate);
            % Ascending order: the biggest dominating rates (numerically)
            % produce on average the most dominating process arrivals.
            [dominatingRate, sortedInds] = sort(dominatingRate, 1, 'ascend');
        elseif ModelSpec.modelIndicator == 1,
            sortedInds = (1:TempAlgorithmParams.nParticlesEachBlock(iSubpop, 1))';
        end
        
        % Algorithm parameters.
        TempAlgorithmParams.logParticleWeights{iSubpop, 1} = cell2mat(AlgorithmParams.logParticleWeights(blocksThisSubpopBinInds));
        TempAlgorithmParams.logParticleWeights{iSubpop, 1} = TempAlgorithmParams.logParticleWeights{iSubpop, 1}(sortedInds);
        TempAlgorithmParams.logNormalisedParticleWeights{iSubpop, 1} = cell2mat(AlgorithmParams.logNormalisedParticleWeights(blocksThisSubpopBinInds));
        TempAlgorithmParams.logNormalisedParticleWeights{iSubpop, 1} = TempAlgorithmParams.logNormalisedParticleWeights{iSubpop, 1}(sortedInds);
        
        % MJP path structures.
        maxNTransitionsCell = num2cell(maxNTransitions(blocksThisSubpopBinInds))';
        nParticlesCell = num2cell(AlgorithmParams.nParticlesEachBlock(blocksThisSubpopBinInds))';
        tempStateTrajectory{1, iSubpop} = cell2mat(cellfun(@(stateTrajectoryMat, maxNJumps, nParticles) [stateTrajectoryMat; nan([maxMaxNJumps - maxNJumps, nParticles])], stateTrajectory(blocksThisSubpopBinInds), maxNTransitionsCell, nParticlesCell, 'UniformOutput', false));
        tempStateTrajectory{1, iSubpop} = tempStateTrajectory{1, iSubpop}(:, sortedInds);
        tempMaxNTransitions(iSubpop, 1) = size(tempStateTrajectory{1, iSubpop}, 1);
        tempTrajectoryMatBinInds{1, iSubpop} = cell2mat(cellfun(@(binIndsMat, maxNJumps, nParticles) [binIndsMat; false([maxMaxNJumps - maxNJumps, nParticles])], trajectoryMatBinInds(blocksThisSubpopBinInds), maxNTransitionsCell, nParticlesCell, 'UniformOutput', false));
        tempTrajectoryMatBinInds{1, iSubpop} = tempTrajectoryMatBinInds{1, iSubpop}(:, sortedInds);
        tempJumpTimes{iSubpop, 1} = cat(1, jumpTimes{blocksThisSubpopBinInds});
        tempJumpTimes{iSubpop} = tempJumpTimes{iSubpop}(sortedInds);
        if size(tempJumpTimes{iSubpop}, 2) > 1,
            tempJumpTimes{iSubpop, 1} = tempJumpTimes{iSubpop}';
        end
        tempInterjumpIntervals{iSubpop, 1} = cat(1, interjumpIntervals{blocksThisSubpopBinInds});
        tempInterjumpIntervals{iSubpop, 1} = tempInterjumpIntervals{iSubpop}(sortedInds);
        if size(tempInterjumpIntervals{iSubpop}, 2) > 1,
            tempInterjumpIntervals{iSubpop, 1} = tempInterjumpIntervals{iSubpop}';
        end
        if havePartialMarginalLikelihoods,
            tempLogPartialMarginalLikelihoodsCell{1, iSubpop} = cell2mat(logPartialMarginalLikelihoodsCell(blocksThisSubpopBinInds));
            tempLogPartialMarginalLikelihoodsCell{1, iSubpop} = tempLogPartialMarginalLikelihoodsCell{1, iSubpop}(sortedInds);
        end
        
        % Model parameter structures.
        if ModelSpec.modelIndicator == 2,
            TempModelParams(iSubpop, 1).generatorMat = cat(3, ModelParams(blocksThisSubpopBinInds).generatorMat);
            TempModelParams(iSubpop, 1).generatorMat = TempModelParams(iSubpop, 1).generatorMat(:, :, sortedInds);
            TempModelParams(iSubpop, 1).jumpRate = cat(2, ModelParams(blocksThisSubpopBinInds).jumpRate);
            TempModelParams(iSubpop, 1).jumpRate = TempModelParams(iSubpop, 1).jumpRate(:, sortedInds);
        end
        TempModelParams(iSubpop, 1).logInitialStateDist = cat(2, ModelParams(blocksThisSubpopBinInds).logInitialStateDist);
        TempModelParams(iSubpop, 1).logInitialStateDist = TempModelParams(iSubpop, 1).logInitialStateDist(:, sortedInds);
        TempModelParams(iSubpop, 1).logTransitionMat = cat(3, ModelParams(blocksThisSubpopBinInds).logTransitionMat);
        TempModelParams(iSubpop, 1).logTransitionMat = TempModelParams(iSubpop, 1).logTransitionMat(:, :, sortedInds);
        if ModelSpec.spikeTrainModelIndicator > 0,
            TempModelParams(iSubpop, 1).logSpikeRate = cat(3, ModelParams(blocksThisSubpopBinInds).logSpikeRate);
            TempModelParams(iSubpop, 1).logSpikeRate = TempModelParams(iSubpop, 1).logSpikeRate(:, :, sortedInds);
            if haveNonspikeRate,
                TempModelParams(iSubpop, 1).logNonspikeRate = cat(3, ModelParams(blocksThisSubpopBinInds).logNonspikeRate);
                TempModelParams(iSubpop, 1).logNonspikeRate = TempModelParams(iSubpop, 1).logNonspikeRate(:, :, sortedInds);
            end
        end
        if ModelSpec.positionModelIndicator > 0,
            TempModelParams(iSubpop, 1).driftRate = cat(2, ModelParams(blocksThisSubpopBinInds).driftRate);
            TempModelParams(iSubpop, 1).driftRate = TempModelParams(iSubpop, 1).driftRate(:, sortedInds);
%             TempModelParams(iSubpop, 1).instDriftRate = cat(2, ModelParams(blocksThisSubpopBinInds).instDriftRate);
%             TempModelParams(iSubpop, 1).instDriftRate = TempModelParams(iSubpop, 1).driftRate(:, sortedInds);
            TempModelParams(iSubpop, 1).asympMeanPosInd = cat(2, ModelParams(blocksThisSubpopBinInds).asympMeanPosInd);
            TempModelParams(iSubpop, 1).asympMeanPosInd = TempModelParams(iSubpop, 1).asympMeanPosInd(:, sortedInds);
            if (ModelSpec.positionJitterModelIndicator == 0) || (ModelSpec.positionJitterModelIndicator == 1) || (ModelSpec.positionJitterModelIndicator == 2),
                TempModelParams(iSubpop, 1).jitterCov = cat(3, ModelParams(blocksThisSubpopBinInds).jitterCov);
                TempModelParams(iSubpop, 1).jitterCov = TempModelParams(iSubpop, 1).jitterCov(:, :, sortedInds);
                TempModelParams(iSubpop, 1).posCov = cat(3, ModelParams(blocksThisSubpopBinInds).posCov);
                TempModelParams(iSubpop, 1).posCov = TempModelParams(iSubpop, 1).posCov(:, :, sortedInds);
%                 TempModelParams(iSubpop, 1).stationaryPosCov = cat(3, ModelParams(blocksThisSubpopBinInds).stationaryPosCov);
%                 TempModelParams(iSubpop, 1).stationaryPosCov = TempModelParams(iSubpop, 1).stationaryPosCov(:, :, sortedInds);
%                 TempModelParams(iSubpop, 1).instJitterCov = cat(3, ModelParams(blocksThisSubpopBinInds).instJitterCov);
%                 TempModelParams(iSubpop, 1).instJitterCov = TempModelParams(iSubpop, 1).instJitterCov(:, :, sortedInds);
            elseif ModelSpec.positionJitterModelIndicator == 3,
                TempModelParams(iSubpop, 1).jitterCov = cat(4, ModelParams(blocksThisSubpopBinInds).jitterCov);
                TempModelParams(iSubpop, 1).jitterCov = TempModelParams(iSubpop, 1).jitterCov(:, :, :, sortedInds);
                TempModelParams(iSubpop, 1).posCov = cat(4, ModelParams(blocksThisSubpopBinInds).posCov);
                TempModelParams(iSubpop, 1).posCov = TempModelParams(iSubpop, 1).posCov(:, :, :, sortedInds);
%                 TempModelParams(iSubpop, 1).stationaryPosCov = cat(4, ModelParams(blocksThisSubpopBinInds).stationaryPosCov);
%                 TempModelParams(iSubpop, 1).stationaryPosCov = TempModelParams(iSubpop, 1).stationaryPosCov(:, :, :, sortedInds);
            end
%             if ModelSpec.modelIndicator == 2,
%                 TempModelParams(iSubpop, 1).posObsVar = cat(1, ModelParams(blocksThisSubpopBinInds).posObsVar);
%                 TempModelParams(iSubpop, 1).posObsVar = TempModelParams(iSubpop, 1).posObsVar(sortedInds);
%             end
        end
        % Sorted already.
        if ModelSpec.modelIndicator == 2,
            TempModelParams(iSubpop, 1).dominatingRate = dominatingRate;
        end
        TempModelParams(iSubpop, 1).nParticlesThisBlock = TempAlgorithmParams.nParticlesEachBlock(iSubpop, 1);
        TempModelParams(iSubpop, 1).nStatesThisBlock = size(TempModelParams(iSubpop, 1).logTransitionMat, 1);
        TempModelParams(iSubpop, 1).nTrueStatesThisBlock = iSubpop;

        if haveFirstCopyBinInds,
            tempFirstCopyBinInds{iSubpop, 1} = cell2mat(firstCopyBinInds(blocksThisSubpopBinInds));
            tempFirstCopyBinInds{iSubpop, 1} = tempFirstCopyBinInds{iSubpop, 1}(sortedInds);
        end
    end

    % Clear originals.
    clear AlgorithmParams ModelParams stateTrajectory trajectoryMatBinInds jumpTimes interjumpIntervals logPartialMarginalLikelihoodsCell;
    
    % Upper limit on how many blocks we will need.
    atMostBlocksNeeded = ceil(AlgorithmSpec.nParticles / maxBlockSize) + nSubpops;

    NewAlgorithmParams.nParticlesEachBlock = zeros([atMostBlocksNeeded, 1]);
    NewAlgorithmParams.nTrueStatesEachBlock = zeros([atMostBlocksNeeded, 1]);
    NewAlgorithmParams.logParticleWeights = cell([atMostBlocksNeeded, 1]);
    NewAlgorithmParams.logNormalisedParticleWeights = cell([atMostBlocksNeeded, 1]);

    newStateTrajectory = cell([1, atMostBlocksNeeded]);
    newTrajectoryMatBinInds = cell([1, atMostBlocksNeeded]);
    newJumpTimes = cell([atMostBlocksNeeded, 1]);
    newInterjumpIntervals = cell([atMostBlocksNeeded, 1]);
    newMaxNTransitions = zeros([atMostBlocksNeeded, 1]);
    if havePartialMarginalLikelihoods,
        newLogPartialMarginalLikelihoodsCell = cell([1, atMostBlocksNeeded]);
    end
    
    if ModelSpec.modelIndicator == 2,
        [NewModelParams(1:atMostBlocksNeeded).generatorMat] = deal([]);
        [NewModelParams(1:atMostBlocksNeeded).jumpRate] = deal([]);
        [NewModelParams(1:atMostBlocksNeeded).dominatingRate] = deal([]);
    end
    [NewModelParams(1:atMostBlocksNeeded).logInitialStateDist] = deal([]);
    [NewModelParams(1:atMostBlocksNeeded).logTransitionMat] = deal([]);
    if ModelSpec.spikeTrainModelIndicator > 0,
        [NewModelParams(1:atMostBlocksNeeded).logSpikeRate] = deal([]);
        if haveNonspikeRate,
            [NewModelParams(1:atMostBlocksNeeded).logNonspikeRate] = deal([]);
        end
    end
    if ModelSpec.positionModelIndicator > 0,
        [NewModelParams(1:atMostBlocksNeeded).driftRate] = deal([]);
%         [NewModelParams(1:atMostBlocksNeeded).instDriftRate] = deal([]);
        [NewModelParams(1:atMostBlocksNeeded).asympMeanPosInd] = deal([]);
        [NewModelParams(1:atMostBlocksNeeded).posCov] = deal([]);
%         [NewModelParams(1:atMostBlocksNeeded).stationaryPosCov] = deal([]);
%         if (ModelSpec.positionJitterModelIndicator == 0) || (ModelSpec.positionJitterModelIndicator == 1) || (ModelSpec.positionJitterModelIndicator == 2),
            [NewModelParams(1:atMostBlocksNeeded).jitterCov] = deal([]);
%             [NewModelParams(1:atMostBlocksNeeded).instJitterCov] = deal([]);
%         end
%         if ModelSpec.modelIndicator == 2,
%             [NewModelParams(1:atMostBlocksNeeded).posObsVar] = deal([]);
%         end
    end
    NewModelParams = NewModelParams';
    
    if haveFirstCopyBinInds,
        newFirstCopyBinInds = cell([atMostBlocksNeeded, 1]);
    end
    
    % Now we partition particles into new blocks with maxBlockSize
    % particles. Any surplus in each subpopulation goes into a smaller
    % block.
    % This tells us which position in the new blocked structures is where
    % the subpopulation's new blocks should start from.
    startBlockInd = 1;
    for iSubpop = 1:ModelSpec.maxNStates,
        if (TempAlgorithmParams.nBlocksEachSubpop(iSubpop) == 0) || (TempModelParams(iSubpop).nParticlesThisBlock == 0),
            continue;
        end
        % Doing this case separately suppresses warnings about having too
        % many indices.
        if TempModelParams(iSubpop).nParticlesThisBlock == 1,
            NewAlgorithmParams.nParticlesEachBlock(startBlockInd) = TempAlgorithmParams.nParticlesEachBlock(iSubpop);
            NewAlgorithmParams.nTrueStatesEachBlock(startBlockInd) = TempAlgorithmParams.nTrueStatesEachBlock(iSubpop);
            NewAlgorithmParams.logParticleWeights{startBlockInd} = TempAlgorithmParams.logParticleWeights{iSubpop};
            NewAlgorithmParams.logNormalisedParticleWeights{startBlockInd} = TempAlgorithmParams.logNormalisedParticleWeights{iSubpop};

            newStateTrajectory{1, startBlockInd} = tempStateTrajectory{iSubpop};
            newTrajectoryMatBinInds{1, startBlockInd} = tempTrajectoryMatBinInds{iSubpop};
            newJumpTimes{startBlockInd, 1} = tempJumpTimes{iSubpop};
            newInterjumpIntervals{startBlockInd, 1} = tempInterjumpIntervals{iSubpop};
            newMaxNTransitions(startBlockInd, 1) = size(newStateTrajectory{1, startBlockInd}, 1);
            if havePartialMarginalLikelihoods,
                newLogPartialMarginalLikelihoodsCell{1, startBlockInd} = tempLogPartialMarginalLikelihoodsCell{iSubpop};
            end

            if haveFirstCopyBinInds,
                newFirstCopyBinInds{startBlockInd, 1} = tempFirstCopyBinInds{iSubpop};
            end

            NewModelParams(startBlockInd, 1).nParticlesThisBlock = NewAlgorithmParams.nParticlesEachBlock(startBlockInd, 1);
            NewModelParams(startBlockInd, 1).nStatesThisBlock = TempModelParams(iSubpop).nStatesThisBlock;
            NewModelParams(startBlockInd, 1).nTrueStatesThisBlock = TempModelParams(iSubpop).nTrueStatesThisBlock;
            if ModelSpec.modelIndicator == 2,
                % Generator.
                [NewModelParams(startBlockInd, 1).generatorMat] = TempModelParams(iSubpop).generatorMat;
                % Jump rate.
                [NewModelParams(startBlockInd, 1).jumpRate] = TempModelParams(iSubpop).jumpRate;
                % Dominating rate.
                [NewModelParams(startBlockInd, 1).dominatingRate] = TempModelParams(iSubpop).dominatingRate;
            end
            % Transition matrix.
            [NewModelParams(startBlockInd, 1).logTransitionMat] = TempModelParams(iSubpop).logTransitionMat;
            % Initial state distribution.
            [NewModelParams(startBlockInd, 1).logInitialStateDist] = TempModelParams(iSubpop).logInitialStateDist;
            if ModelSpec.spikeTrainModelIndicator > 0,
                % Spike rate.                
                [NewModelParams(startBlockInd, 1).logSpikeRate] = TempModelParams(iSubpop).logSpikeRate;
                if haveNonspikeRate,
                    [NewModelParams(startBlockInd, 1).logNonspikeRate] = TempModelParams(iSubpop).logNonspikeRate;
                end
            end
            if ModelSpec.positionModelIndicator > 0,
                [NewModelParams(startBlockInd, 1).driftRate] = TempModelParams(iSubpop).driftRate;
%                 [NewModelParams(startBlockInd, 1).instDriftRate] = TempModelParams(iSubpop).instDriftRate;
                [NewModelParams(startBlockInd, 1).asympMeanPosInd] = TempModelParams(iSubpop).asympMeanPosInd;
                [NewModelParams(startBlockInd, 1).posCov] = TempModelParams(iSubpop).posCov;
%                 [NewModelParams(startBlockInd, 1).stationaryPosCov] = TempModelParams(iSubpop).stationaryPosCov;
%                 if (ModelSpec.positionJitterModelIndicator == 0) || (ModelSpec.positionJitterModelIndicator == 1) || (ModelSpec.positionJitterModelIndicator == 2),
                    [NewModelParams(startBlockInd, 1).jitterCov] = TempModelParams(iSubpop).jitterCov;
%                     [NewModelParams(startBlockInd, 1).instJitterCov] = TempModelParams(iSubpop).instJitterCov;
%                 end
%                 if ModelSpec.modelIndicator == 2,
%                     [NewModelParams(startBlockInd, 1).posObsVar] = TempModelParams(iSubpop).posObsVar;
%                 end
            end

            startBlockInd = startBlockInd + 1;
            continue;
        end

        % Only need to split this subpopulation into new blocks if it is
        % larger than the max block size. Otherwise nBlocksNeeded will be 1
        % and nSurplusParticles will be n particles in this subpopulation.

        % How many blocks will partition this subpopulation?
        nBlocksNeeded = ceil(TempAlgorithmParams.nParticlesEachBlock(iSubpop) / maxBlockSize);
        % Vector of block sizes that partition this subpopulation.
        blockPartitionVec = maxBlockSize .* ones([nBlocksNeeded, 1]);
        % The last block may be smaller - how many particles should go
        % into a smaller block at the end?
%         nSurplusParticles = rem(TempAlgorithmParams.nParticlesEachBlock(iSubpop), nBlocksNeeded - 1);
        nSurplusParticles = TempAlgorithmParams.nParticlesEachBlock(iSubpop) - (nBlocksNeeded - 1) * maxBlockSize;
        if nSurplusParticles > 0,
            blockPartitionVec(end) = nSurplusParticles;
        end

        % Algorithm params.
        NewAlgorithmParams.nParticlesEachBlock(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1) = blockPartitionVec;
        NewAlgorithmParams.nTrueStatesEachBlock(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1) = iSubpop;
        NewAlgorithmParams.logParticleWeights(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1) = mat2cell(TempAlgorithmParams.logParticleWeights{iSubpop}, blockPartitionVec, 1);
        NewAlgorithmParams.logNormalisedParticleWeights(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1) = mat2cell(TempAlgorithmParams.logNormalisedParticleWeights{iSubpop}, blockPartitionVec, 1);

        % MJP path structures.
        newStateTrajectory(1, startBlockInd:(startBlockInd + nBlocksNeeded - 1)) = mat2cell(tempStateTrajectory{iSubpop}, tempMaxNTransitions(iSubpop), blockPartitionVec);
        newTrajectoryMatBinInds(1, startBlockInd:(startBlockInd + nBlocksNeeded - 1)) = mat2cell(tempTrajectoryMatBinInds{iSubpop}, tempMaxNTransitions(iSubpop), blockPartitionVec);
        newJumpTimes(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1) = mat2cell(tempJumpTimes{iSubpop}, blockPartitionVec, 1);
        newInterjumpIntervals(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1) = mat2cell(tempInterjumpIntervals{iSubpop}, blockPartitionVec, 1);
        newMaxNTransitions(startBlockInd:(startBlockInd + nBlocksNeeded - 1), :) = cellfun(@(trajectoryMat) size(trajectoryMat, 1), newStateTrajectory(startBlockInd:(startBlockInd + nBlocksNeeded - 1)))';
        if havePartialMarginalLikelihoods,
%             tempCell = mat2cell(tempLogPartialMarginalLikelihoodsCell{iSubpop}, 1, 1, blockPartitionVec);
%             [newLogPartialMarginalLikelihoodsCell(1, 1, startBlockInd:(startBlockInd + nBlocksNeeded - 1))] = tempCell{:};
            newLogPartialMarginalLikelihoodsCell(1, startBlockInd:(startBlockInd + nBlocksNeeded - 1)) = mat2cell(tempLogPartialMarginalLikelihoodsCell{iSubpop}, 1, blockPartitionVec);
        end

        % Model parameter structures.
        tempCell = num2cell(NewAlgorithmParams.nParticlesEachBlock(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1));
        [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).nParticlesThisBlock] = tempCell{:};
        [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).nStatesThisBlock] = deal(TempModelParams(iSubpop).nStatesThisBlock);
        [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).nTrueStatesThisBlock] = deal(TempModelParams(iSubpop).nTrueStatesThisBlock);
        if ModelSpec.modelIndicator == 2,
            % Generator.
            tempCell = mat2cell(TempModelParams(iSubpop).generatorMat, TempModelParams(iSubpop).nStatesThisBlock, TempModelParams(iSubpop).nStatesThisBlock, blockPartitionVec);
            [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).generatorMat] = tempCell{:};
            % Jump rate.
            tempCell = mat2cell(TempModelParams(iSubpop).jumpRate, TempModelParams(iSubpop).nStatesThisBlock, blockPartitionVec);
            [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).jumpRate] = tempCell{:};
            % Dominating rate.
            tempCell = mat2cell(TempModelParams(iSubpop).dominatingRate, blockPartitionVec, 1);
            [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).dominatingRate] = tempCell{:};
        end
        % Transition matrix.
        tempCell = mat2cell(TempModelParams(iSubpop).logTransitionMat, TempModelParams(iSubpop).nStatesThisBlock, TempModelParams(iSubpop).nStatesThisBlock, blockPartitionVec);
        [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).logTransitionMat] = tempCell{:};
        % Initial state distribution.
        tempCell = mat2cell(TempModelParams(iSubpop).logInitialStateDist, TempModelParams(iSubpop).nStatesThisBlock, blockPartitionVec);
        [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).logInitialStateDist] = tempCell{:};
        if ModelSpec.spikeTrainModelIndicator > 0,
            % Spike rate.
            tempCell = mat2cell(TempModelParams(iSubpop).logSpikeRate, TempModelParams(iSubpop).nTrueStatesThisBlock, nSpikeTrains, blockPartitionVec);
            [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).logSpikeRate] = tempCell{:};
            if haveNonspikeRate,
                tempCell = mat2cell(TempModelParams(iSubpop).logNonspikeRate, TempModelParams(iSubpop).nTrueStatesThisBlock, nSpikeTrains, blockPartitionVec);
                [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).logNonspikeRate] = tempCell{:};
            end
        end
        if ModelSpec.positionModelIndicator > 0,
            tempCell = mat2cell(TempModelParams(iSubpop).driftRate, TempModelParams(iSubpop).nTrueStatesThisBlock, blockPartitionVec);
            [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).driftRate] = tempCell{:};
%             tempCell = mat2cell(TempModelParams(iSubpop).instDriftRate, TempModelParams(iSubpop).nTrueStatesThisBlock, blockPartitionVec);
%             [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).instDriftRate] = tempCell{:};
            tempCell = mat2cell(TempModelParams(iSubpop).asympMeanPosInd, TempModelParams(iSubpop).nTrueStatesThisBlock, blockPartitionVec);
            [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).asympMeanPosInd] = tempCell{:};
            if (ModelSpec.positionJitterModelIndicator == 0) || (ModelSpec.positionJitterModelIndicator == 1) || (ModelSpec.positionJitterModelIndicator == 2),
                tempCell = mat2cell(TempModelParams(iSubpop).jitterCov, 2, TempModelParams(iSubpop).nTrueStatesThisBlock, blockPartitionVec);
                [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).jitterCov] = tempCell{:};
                tempCell = mat2cell(TempModelParams(iSubpop).posCov, 2, TempModelParams(iSubpop).nTrueStatesThisBlock, blockPartitionVec);
                [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).posCov] = tempCell{:};
%                 tempCell = mat2cell(TempModelParams(iSubpop).stationaryPosCov, 2, TempModelParams(iSubpop).nTrueStatesThisBlock, blockPartitionVec);
%                 [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).stationaryPosCov] = tempCell{:};
%                 tempCell = mat2cell(TempModelParams(iSubpop).instJitterCov, 2, TempModelParams(iSubpop).nTrueStatesThisBlock, blockPartitionVec);
%                 [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).instJitterCov] = tempCell{:};
            elseif ModelSpec.positionJitterModelIndicator == 3,
                tempCell = mat2cell(TempModelParams(iSubpop).jitterCov, 2, 2, TempModelParams(iSubpop).nTrueStatesThisBlock, blockPartitionVec);
                [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).jitterCov] = tempCell{:};
                tempCell = mat2cell(TempModelParams(iSubpop).posCov, 2, 2, TempModelParams(iSubpop).nTrueStatesThisBlock, blockPartitionVec);
                [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).posCov] = tempCell{:};
%                 tempCell = mat2cell(TempModelParams(iSubpop).stationaryPosCov, 2, 2, TempModelParams(iSubpop).nTrueStatesThisBlock, blockPartitionVec);
%                 [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).stationaryPosCov] = tempCell{:};
            end
%             if ModelSpec.modelIndicator == 2,
%                 tempCell = mat2cell(TempModelParams(iSubpop).posObsVar, blockPartitionVec, 1);
%                 [NewModelParams(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1).posObsVar] = tempCell{:};
%             end
        end

        if haveFirstCopyBinInds,
            try
            newFirstCopyBinInds(startBlockInd:(startBlockInd + nBlocksNeeded - 1), 1) = mat2cell(tempFirstCopyBinInds{iSubpop}, blockPartitionVec, 1);
            catch
                blockPartitionVec
                size(tempFirstCopyBinInds{iSubpop})
                firstCopyBinInds
                error('Here.')
            end
        end

        startBlockInd = startBlockInd + nBlocksNeeded;

    end

    % Remove empty blocks.
    usedBlockBinInds = NewAlgorithmParams.nParticlesEachBlock > 0;
    
    NewAlgorithmParams.nParticlesEachBlock = NewAlgorithmParams.nParticlesEachBlock(usedBlockBinInds);
    NewAlgorithmParams.nTrueStatesEachBlock = NewAlgorithmParams.nTrueStatesEachBlock(usedBlockBinInds);
    NewAlgorithmParams.logParticleWeights = NewAlgorithmParams.logParticleWeights(usedBlockBinInds);
    NewAlgorithmParams.logNormalisedParticleWeights = NewAlgorithmParams.logNormalisedParticleWeights(usedBlockBinInds);

    newStateTrajectory = newStateTrajectory(usedBlockBinInds);
    newTrajectoryMatBinInds = newTrajectoryMatBinInds(usedBlockBinInds);
    newJumpTimes = newJumpTimes(usedBlockBinInds);
    newInterjumpIntervals = newInterjumpIntervals(usedBlockBinInds);
    newMaxNTransitions = newMaxNTransitions(usedBlockBinInds);
    if havePartialMarginalLikelihoods,
        newLogPartialMarginalLikelihoodsCell = newLogPartialMarginalLikelihoodsCell(usedBlockBinInds);
    end
    
    NewModelParams = NewModelParams(usedBlockBinInds);
    
    if haveFirstCopyBinInds,
        newFirstCopyBinInds = newFirstCopyBinInds(usedBlockBinInds);
    end

    % Finish NewAlgorithmParams.
    NewAlgorithmParams.nBlocks = length(NewModelParams);
    NewAlgorithmParams.nBlocksEachSubpop = zeros([ModelSpec.maxNStates, 1]);
    NewAlgorithmParams.nParticlesEachSubpop = zeros([ModelSpec.maxNStates, 1]);
    for iSubpop = 1:ModelSpec.maxNStates,
        NewAlgorithmParams.nBlocksEachSubpop(iSubpop) = sum(NewAlgorithmParams.nTrueStatesEachBlock == iSubpop);
        NewAlgorithmParams.nParticlesEachSubpop(iSubpop) = sum(NewAlgorithmParams.nParticlesEachBlock(NewAlgorithmParams.nTrueStatesEachBlock == iSubpop));
    end
    
    if nargout == 7,
        varargout = {NewModelParams, NewAlgorithmParams, newStateTrajectory, newTrajectoryMatBinInds, newJumpTimes, newInterjumpIntervals, newMaxNTransitions};
    elseif nargout == 8,
        varargout = {NewModelParams, NewAlgorithmParams, newStateTrajectory, newTrajectoryMatBinInds, newJumpTimes, newInterjumpIntervals, newMaxNTransitions, newLogPartialMarginalLikelihoodsCell};
    elseif nargout == 9,
        varargout = {NewModelParams, NewAlgorithmParams, newStateTrajectory, newTrajectoryMatBinInds, newJumpTimes, newInterjumpIntervals, newMaxNTransitions, newLogPartialMarginalLikelihoodsCell, newFirstCopyBinInds};
    end
    
end