function [samplingWeights, discriminateSubpopLog] = positivediscrimination( AlgorithmParams, estLogStateDimensionPosterior, AlgorithmSpec, maxNStates )
    % Apply the 'positive discrimination' algorothim of Chopin (2007) to
    % keep subpopulations of particles from going 'extinct' as a result of
    % resampling.
    % We return adjusted particle weights that should be used in the
    % resampling algorithm.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % AlgorithmParams.logNormalisedParticleWeights - cell array of log
    % particle weights. Each cell is a 'block' corresponding to a paticular
    % state dimension. Multiple blocks may represent the same dimension.
    % Should be normalised (i.e. over all blocks).
    %
    % AlgorithmParams.nTrueStatesEachBlock - row vector of the state
    % dimensions represented by each particle block.
    %
    % estLogStateDimensionPosterior - column vector of log posterior
    % probabilities over state dimension recently estimated using the
    % particle population.
    %
    % AlgorithmSpec - a structure containing the following fields: -
    % allowSmallSubpopulationsToGoExtinct - logical. -
    % smallSubpopulationBuffer - nonnegative integer: state dimensions of
    % label at least the estimated MAP minus this buffer are not allowed to
    % go extinct. - positiveDiscriminationThreshold: the threshold
    % posterior probability of a state dimension, below which positive
    % discrimination is triggered for the corresponding subpopulation of
    % particles.
    %
    % smallSubpopulationBuffer - cell array the same size as
    % logNormalisedParticleWeights, with cells of the same size and
    % dimensions. Contains adjusted weights, normalised, that should be
    % used in the resampling algorithm, and put into effect the positive
    % discrimination.
    %
    % samplingWeights - will be the same size and dimension of cell array
    % as logNormalisedParticleWeights, with cells of the same size and
    % dimension, but will NOT be in the log domain.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 16.10.2014 - copied from old version.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    if nargin < 4,
        error('Insufficient parameters!')
    end
    if nargin < 5,
        maxNStates = length(estLogStateDimensionPosterior);
    else
        if length(estLogStateDimensionPosterior) ~= maxNStates,
            error('Number of state dimensions implied by estLogStateDimensionPosterior does not match maxNStates!')
        end
    end
    if any((rem(AlgorithmParams.nTrueStatesEachBlock, 1)) ~= 0) || any(AlgorithmParams.nTrueStatesEachBlock <= 0),
        error('AlgorithmParams.nTrueStatesEachBlock should be a vector of positive integers!')
    end
    if max(AlgorithmParams.nTrueStatesEachBlock) > maxNStates,
        error('Values in AlgorithmParams.nTrueStatesEachBlock must be bounded above by maxNStates!')
    end
    if length(AlgorithmParams.logNormalisedParticleWeights) ~= length(AlgorithmParams.nTrueStatesEachBlock),
        error('AlgorithmParams.logNormalisedParticleWeights and AlgorithmParams.nTrueStatesEachBlock must have the same length (number of blocks)!')
    end
    if any(cellfun(@(iVec) any(iVec > 0), AlgorithmParams.logNormalisedParticleWeights)),
        error('Cells of AlgorithmParams.logNormalisedParticleWeights must be vectors of log probabilities!')
    end
    if ~isfield(AlgorithmSpec, 'allowSmallSubpopulationsToGoExtinct'),
        error('allowSmallSubpopulationsToGoExtinct must be a field of AlgorithmSpec!')
    else
        if ~islogical(AlgorithmSpec.allowSmallSubpopulationsToGoExtinct),
            error('AlgorithmSpec.allowSmallSubpopulationsToGoExtinct should be a logical flag!')
        end
        if ~isfield(AlgorithmSpec, 'smallSubpopulationBuffer'),
            error('smallSubpopulationBuffer must be a field of AlgorithmSpec!')
        else
            if (rem(AlgorithmSpec.allowSmallSubpopulationsToGoExtinct, 1)) ~= 0 || (AlgorithmSpec.allowSmallSubpopulationsToGoExtinct < 0),
                error('AlgorithmSpec.allowSmallSubpopulationsToGoExtinct should be a nonnegative integer!')
            end
        end
    end
    if ~isfield(AlgorithmSpec, 'positiveDiscriminationThreshold'),
        error('positiveDiscriminationThreshold must be a field of AlgorithmSpec!')
    else
        if (AlgorithmSpec.positiveDiscriminationThreshold < 0) || (AlgorithmSpec.positiveDiscriminationThreshold > 1),
            error('AlgorithmSpec.positiveDiscriminationThreshold must be a valid probability!')
        end
    end
    
    % Which subpopulations require positive discrimination because their
    % estimated posterior state dimension probability is low?
    discriminateSubpopLog = (estLogStateDimensionPosterior < log(AlgorithmSpec.positiveDiscriminationThreshold)) & ~isinf(estLogStateDimensionPosterior);
    if AlgorithmSpec.allowSmallSubpopulationsToGoExtinct,
        % Compute the MAP state dimension.
        [~, estMapStateDimension] = max(estLogStateDimensionPosterior);
        subpopDiscriminationOkLog = (1:maxNStates) >= (estMapStateDimension - AlgorithmSpec.smallSubpopulationBuffer);
        discriminateSubpopLog = discriminateSubpopLog & subpopDiscriminationOkLog;
    end
    nSubPopsDiscriminated = sum(discriminateSubpopLog);

    % Normalise weights of subpopulations to be discriminated for.
    normalisedSubpopLogParticleWeights = AlgorithmParams.logNormalisedParticleWeights;
    for iStateDimension = 1:maxNStates,
        % To identify the blocks that constitute this subpopulation.
        stateDimLog = AlgorithmParams.nTrueStatesEachBlock == iStateDimension;
        if sum(stateDimLog) == 0,
            continue;
        end
        nParticlesEachBlock = cellfun(@length, AlgorithmParams.logNormalisedParticleWeights(stateDimLog));
        if discriminateSubpopLog(iStateDimension),
            % Normalisation must be over all blocks that constitute this
            % subpopulation.
            normalisedSubpopLogParticleWeights(stateDimLog) = mat2cell(normaliselogdistributionsmatrix( cell2mat(AlgorithmParams.logNormalisedParticleWeights(stateDimLog)), 1 ), nParticlesEachBlock, 1);
        else
            % If not, set all weights to zero (so they won't be
            % discriminated for).
%             normalisedSubpopLogParticleWeights(stateDimLog) = cellfun(@(iBlock) -Inf * iBlock, normalisedSubpopLogParticleWeights(stateDimLog), 'UniformOutput', false);
            nParticlesTotal = sum(nParticlesEachBlock);
            normalisedSubpopLogParticleWeights(stateDimLog) = mat2cell(-inf([nParticlesTotal, 1]), nParticlesEachBlock, 1);
        end
    end

    % Matrix version.
%     samplingWeights = exp(AlgorithmParams.logNormalisedParticleWeights) * (1 - nSubPopsDiscriminated) * AlgorithmSpec.positiveDiscriminationThreshold ...
%                             + exp(normalisedSubpopLogParticleWeights) * AlgorithmSpec.positiveDiscriminationThreshold;
    % Cell array version.
    samplingWeights = cellfun(@(iBlockWeights, iBlockDiscriminatedWeights) ...
                                        exp(iBlockWeights) * (1 - nSubPopsDiscriminated * AlgorithmSpec.positiveDiscriminationThreshold) ...
                                        + exp(iBlockDiscriminatedWeights) * AlgorithmSpec.positiveDiscriminationThreshold, ...
                                    AlgorithmParams.logNormalisedParticleWeights, normalisedSubpopLogParticleWeights, 'UniformOutput', false);
end