function postMean = computeposteriormeanasympmean( prevPostMean, res, PositionData, driftRate, currentPosCountEachState, previousPosCountEachState, nExtendedStates, weightedPosCountEachStateByTarget, stateInds, targetPosInds )
    % This function is designed to be called recursively. We compute the
    % posterior mean asymptotic mean parameter, in a geometrically
    % constrained environment, for multiple particles.
    %
    % We use the 'sandpile' algorithm at a very low resolution. In every
    % iteration, there are at most 4 positions to choose between. The
    % chosen position is divided into 4, those subregions assessed for
    % being valid or not ('inside' or 'outside' the maze), and those
    % positions passed to the next iteration. This is repeated until we
    % reach the desired (or maximum) resolution.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % prevPostMean - 2D matrix (states index rows, particles columns) of
    % the mean positions identified for each sample at the previous
    % iteration (one level coarser resolution).
    %
    % res - the current resolution: there are 2^res discretisation grid
    % cells covering the longest spatial dimension; equivalently, the
    % longest spatial dimension has been divided into parts of size
    % 2^(-res) times the whole length.
    %
    % currentPosCountEachState, previousPosCountEachState - 3-D matrices
    % listing the count of visits to each position whilst in each state,
    % and for each particle. Former includes all observations AFTER start
    % time. Latter includes position at start time (initial position) but
    % does NOT include the last position observation.
    %
    % PositionData.sandMats - cell array of the 'sand' matrices for each resolution.
    % Resolution indexes cells. Each cell is a matrix with D rows and d
    % columns, where D is the number of valid positions at the highest
    % resolution, and d is the number of valid positions at the current
    % resolution.
%     % The last cell is the full, highest resolution sand
%     % matrix, which is sparse and upper triangular.
    %
    % PositionData.resShifterMats - a cell array, with cells indexed by resolution of 2D
    % matrices that convert discrete positions at the previous resolution
    % (coarser, indexing rows) to discrete positions at the current
    % resolution (i.e., res, from 1 to maximum). Thus there are 4 columns
    % for each such matrix, and the first cell has only one row (zeroth
    % resolution is just one discretisation cell). Only valid positions are
    % listed in the columns. If there are less than 4, the remainder is
    % padded with NaNs.
    %
    % PositionData.validDiscretePositionsVector - column vector of the
    % valid position indices.
    %
    % Note: When multiple positions are tied for the sample sum of 'sand',
    % the tie is broken by selecting as mean position the position of
    % lowest linear index. (But small numerical errors may upset this.)
    %
    % Note: a discrete position - at any resolution - is considered valid
    % if it contains a discrete position at the highest resolution that is
    % known to be valid via validPositionsLogMat.
    %
    % See also computemeanpositionlowres(). That function is less general;
    % we can reproduce its behaviour here by properly converting the
    % discrete-time position data into durations spent at each position,
    % given knowledge of the discretisation bin width.
    % computemeanpositionlowres is also slower and bad for memory when
    % using large data sets.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % CHANGELOG
    % 30.10.2014 - copied from old version,
    % computemeanpositionlowrescontinuoustime.m.
    % 31.10.2014 - changed position count matrices to be indexed by state
    % and particle separately (was state * particle before). Coded new
    % model: the mean is of a weighting of observations with the
    % observations immediately preceding ('current' with 'previous'
    % observations), rather than a straight mean of observed positions. The
    % weighting is by drift rate.
    % 04.11.2014 - changed back to old approach using 'extended states':
    % only those state/particle combinations that actually have
    % observations and thus require a posterior mean.
    % 27.11.2014 - madeweightedPosCountEachStateByTarget, stateInds and
    % targetPosInds optional inputs so they aren't computed (unnecessarily)
    % more than once.
    %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if nargin < 6,
        error('Insufficient arguments!')
    end
    if ~isfield(PositionData, 'sandMats'),
        error('sandMats must be a field of PositionData!')
    end
    if ~isfield(PositionData, 'resShifterMats'),
        error('resShifterMats must be a field of PositionData!')
    end
    maxRes = length(PositionData.sandMats);
    if size(currentPosCountEachState, 1) ~= PositionData.nValidDiscretePositions,
        error('currentPosCountEachState must be a matrix with rows corresponding to valid positions (at the highest resolution)!')
    end
    if nargin < 7,
        nExtendedStates = size(currentPosCountEachState, 2);
    elseif nExtendedStates ~= size(currentPosCountEachState, 2),
        error('currentPosCountEachState must be a matrix with columns corresponding to ''extended'' states!')
    end
    if nargin < 10,
        % We multiply by number of times each position is visited as a
        % weighted sum of 'current' observation and 'previous' observation.
        weightedPosCountEachStateByTarget = repmat(permute(currentPosCountEachState - bsxfun(@times, driftRate', previousPosCountEachState), [1, 3, 2]), [1, 4, 1]);
%         weightedPosCountEachStateByTarget = repmat(permute(abs(currentPosCountEachState - bsxfun(@times, driftRate', previousPosCountEachState)), [1, 3, 2]), [1, 4, 1]);
%         weightedPosCountEachStateByTarget = repmat(permute(currentPosCountEachState + bsxfun(@times, driftRate', previousPosCountEachState), [1, 3, 2]), [1, 4, 1]);
        % Used for an accumarray.
        stateInds = repmat(permute(1:nExtendedStates, [1, 3, 2]), [PositionData.nValidDiscretePositions, 4, 1]);
        targetPosInds = repmat(1:4, [PositionData.nValidDiscretePositions, 1, nExtendedStates]);
    end
    if ~isequal(size(currentPosCountEachState), size(previousPosCountEachState)),
        error('currentPosCountEachState and previousPosCountEachState must be the same size!')
    end
    if ~isequal(size(driftRate), [nExtendedStates, 1]),
        error('driftRate should be a column vector with nExtendedStates elements!')
    end
    if (res > 1) && (isempty(prevPostMean)),
        error('Target regions at resolution res - 1 must be supplied!')
    end

    if res == 0,
        postMean = ones([nExtendedStates, 1]);
        return;
    end
    
    % The first thing we must do is take the target region of each sample
    % and identify the new targets at this resolution, rejecting them if
    % invalid.

    % For each sample, the positions we consider valid targets.
    if res == 1,
        targetPosMat = repmat(PositionData.resShifterMats{1}, [nExtendedStates, 1]);
    else
        targetPosMat = PositionData.resShifterMats{res}(prevPostMean, :);
    end
    validTargetsBinInds = targetPosMat > 0;

    % This gives us a 3D matrix, with rows indicating the discrete
    % positions at the highest resolution, 4 columns corresponding to the
    % target positions, and pages corresponding to the extended states.
    sandMatLinearInds = bsxfun(@plus, (permute(targetPosMat, [3, 2, 1]) - 1) * PositionData.nValidDiscretePositions, PositionData.validDiscretePositionsVector);
    validTargetsByPosIndBinInds = sandMatLinearInds > 0;
    % This is a column vector giving the sand contribution to each target
    % position.
    sandContribToTargets = PositionData.sandMats{res}(sandMatLinearInds(validTargetsByPosIndBinInds));
%     weightedPosCountEachStateByTarget = repmat(permute(currentPosCountEachState - bsxfun(@times, driftRate', previousPosCountEachState), [1, 3, 2]), [1, 4, 1]);
%     stateInds = repmat(permute(1:nExtendedStates, [1, 3, 2]), [PositionData.nValidDiscretePositions, 4, 1]);
%     targetPosInds = repmat(1:4, [PositionData.nValidDiscretePositions, 1, nExtendedStates]);
    sampleSandContrib = sandContribToTargets .* weightedPosCountEachStateByTarget(validTargetsByPosIndBinInds);
    
    % Add the 'sand' contributions to each target position.
    sumSandContribMat = accumarray([stateInds(validTargetsByPosIndBinInds), targetPosInds(validTargetsByPosIndBinInds)], sampleSandContrib, [nExtendedStates, 4]);
    clear validTargetsByPosIndBinInds sampleSandContrib;
    
    % We want to find the maximum; some could be negative but zeros may
    % exist in positions of nonvalid targets. Make these -inf so we can
    % find the maximum of the valid targets.
    sumSandContribMat(~validTargetsBinInds) = -inf;
        
    % Find target position of maximum accumulated 'sand' for each sample.
    % This returns the indices of the target positions as a matrix, indexed
    % by state and particle.
    [~, meanPosInds] = max(sumSandContribMat, [], 2);
    
    % Identify the 'mean' position of each sample. This first line returns
    % the linear indices of the targetPosMat corresponding to the mean
    % position indices found above.
    targetPosMatInds = sub2ind([nExtendedStates, 4], (1:nExtendedStates)', meanPosInds);
    % Obtains a column vector.
    postMean = targetPosMat(targetPosMatInds);

    % If this is not the highest resolution, use the 'mean' position as the
    % target region of a higher resolution discretisation.
    if res < maxRes,
        clear sumSandContribMat targetPosMat;
        postMean = computeposteriormeanasympmean( postMean, res + 1, PositionData, driftRate, currentPosCountEachState, previousPosCountEachState, nExtendedStates, weightedPosCountEachStateByTarget, stateInds, targetPosInds );
    end
end
