function [WG,I,WSMP,WSTR,ess,esss] = adaptation(wg,essTh,method)
%
% WG        NSTR-by-NSMP matrix, where NSTR is the number of 
%           strata and NSMP is the number of samples per stratum, i.e.
%           NSTR * NSMP = N.
%           
%           Each row of WG defines a resampling distribution within single
%           stratum. NOTE: The rows are not normalised to probabilities.
%
% I         NSTR-by-NSMP matrix of sample indices that determine how the
%           samples are allocated in the strata. ith row of I contains the
%           indices of the samples contained in the ith stratum.
%
% WSMP      N-by-1 array where the ith element is the weight of the stratum
%           associated with the ith particle. The particles are assumed
%           sorted such that the first k samples belong to the first
%           stratum , the next k in the second etc. Ordering within stratum
%           is irrelevant because the weights within stratum are equal.
%
% WSTR      NSTR-by-1 array of stratum weights. The elements of this array
%           are the normalisation terms that make the rows of WG
%           probability distributions.
%
% ess       The effective sample size the just exceeds the threshold.
%
% esss      All effective sample sizes throughout the iteration.

TOL = 1e-5; % tolerance for comapring doubles (ESS thresholds)

I = (1:length(wg))';
N = length(wg);

% Force single stratum for non adapting bootstrap
if strcmp(method,'bootstrap')
    wg = wg';
    I = I';    
end

% Iteration counter
it = 1; 
% Preallocate storage for ESS values in different iterations
esss = zeros(log2(N)+1,1); 
found = false;
num = mean(wg(:))^2;
while true
    
    % Number of strata and samples within a stratum
    [nStr,k] = size(wg); 
    
    % Weigths of the strata
    wstr = sum(wg,2);
        
    % Evaluate ESS
%     esss(it) = (sum(wsmp)^2 / N) / (wsmp(:)' * wsmp(:));
    esss(it) = num / ((wstr' * wstr)/N/k);
%     keyboard
    % If ESS is the first ESS value exceeding the threshold store the
    % result to be set as output.
    if ~found && esss(it) > essTh - TOL
        
        ess = esss(it);
        WG = wg;
        WSTR = wstr;
        % Assign the stratum specific weights to the particles
        wsmp = wstr(ceil((1:N)/k))';
        WSMP = wsmp(:);
        indexer = I;
        found = true;
    end
    
    % End iteration if only one stratum is left
    if nStr == 1
        break;
    end
    
    % Choose stratum coupling strategy
    switch method
        case 'greedy'
            [wg,I] = greedyCoupling(wg,I);
        case 'blockstrategy1'
            [wg,I] = simpleCoupling(wg,I);
        case 'blockstrategy2'
            [wg,I] = simpleCoupling2(wg,I);
        case 'adaptiveBootstrap'
            wg = wg';
            I = I';
        case 'SimpleRearrange'
            [wg,I] = simpleCoupling(wg,I);
        otherwise
            error('unknown method');
    end
    
    % Update iteration counter
    it = it + 1;
    
end

I = indexer;