function nCopiesToRetainEachParticle = residualresampling( particleWeights, nParticlesEachBlock )
    % Will return column vectors of counts of how many times each particle
    % is resampled, in cells of a cell array corresponding to blocks of
    % particles, or as a single cell array, depending on the form of the
    % input.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % particleWeights - in general, a cell array of column vectors, but can
    % be supplied as a single column vector (matrix). Weights should be
    % normalised. Should not be in log domain.
    %
    % nCopiesToRetainEachParticle - a column cell array of column vectors.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 16.10.2014 - copied from old version.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if iscell(particleWeights),
        isCell = true;
        if nargin < 2.
            nParticlesEachBlock = cellfun(@length, particleWeights);
        end
        if any(nParticlesEachBlock == 0),
            error('There are empty particle blocks!')
        end
        % We convert to a column vector.
        particleWeights = cell2mat(particleWeights(:));
    else
        isCell = false;
    end
    
    if any(particleWeights < 0) || any(particleWeights > 1),
        error('particleWeights must be valid probabilities!')
    end
    
    nParticles = length(particleWeights);
    
    proportionalCounts = particleWeights * nParticles;

    % How many copies of each particle should definitely be retained?
    nCopiesToRetainEachParticle = floor(proportionalCounts);

    % This is how many extra times we need to resample particles - dealing
    % with the 'residual'.
    sizeOfResidual = nParticles - sum(nCopiesToRetainEachParticle);

    if sizeOfResidual > 0,
        % This leaves the residual for each particle in sampling probs.
        samplingProbs = proportionalCounts - nCopiesToRetainEachParticle;
        % Normalise, use as sampling distribution.
        samplingProbs = samplingProbs / sum(samplingProbs);

        % Choose which particles to sample with random routine based on
        % normalised weights.
        chosenParticleInds = samplefromcategoricaldistribution( samplingProbs, sizeOfResidual, 1 );

        % Add in the residual resampling counts to the total resampling count
        % for each particle.
        nCopiesToRetainEachParticle = nCopiesToRetainEachParticle + histc(chosenParticleInds, 1:nParticles, 1);
    end
    
    % Split counts into relevant blocks, if required.
    if isCell,
        nCopiesToRetainEachParticle = mat2cell(nCopiesToRetainEachParticle, nParticlesEachBlock, 1);
    end
end