function [X,MU,MUW,ESS,NSTRATA,W,ESSCurves] = alphaSMC(settings)

if settings.resetFigs
    close all
end

% Parameters
N = settings.N;
model = settings.model;
nStep = settings.nStep;
a = settings.param(1);
b = settings.param(2);
sigma = settings.param(3);
nLag = settings.nLag;

% Preallocate data arrays
xs = [zeros(N,nLag-1) randn(N,1)]; % Draw the first sample
xw = xs;
if settings.store
    X = zeros(N,nStep,nLag);
    ESSCurves = zeros(log2(settings.N)+1,nStep);
    NSTRATA = zeros(1,nStep);
    ESS = zeros(1,nStep);
    W = ones(N,nStep); % Stratum weights
else
    W = [];
    X = [];
    ESS = [];
    ESSCurves = [];
    NSTRATA = [];
end

MU = {zeros(nLag,nStep),zeros(nLag,nStep)};

if ~settings.silent
    % Show progress bar for the impatient user
    hwb = waitbar(0,'Running alphaSMC...');
end

Wprev = ones(N,1);
% Main loop
for i = 1:nStep
    
    % Evaluate the likelihoods
    g = likelihood(model,xs(:,end),settings.y(i),b);
    
    if strcmp(settings.method,'SimpleRearrange')
        % Rearrange
        [~,rndI] = sort(rand(N,1));
        xs = xs(rndI,:);
        Wprev = Wprev(rndI);
        g = g(rndI);
    end
    
    % Adaptation
    [wg,I,wsmp,wstr,ess,esscurve] = adaptation(Wprev .* g,...
        settings.essThreshold,settings.method);
    
    % Reorder samples so that the first k samples belong to the first
    % stratum, the next k samples to the second stratum etc.
    xs = xs(I',:);
    [nStrata,k] = size(wg);

    % Sample
    if nStrata == N % i.e. SIS step
        
        % In priciple the case nStrata == N need not be considered
        % separately. This is just a speedup.
        
        % This is for weighted estimates
        c2 = wg ./ sum(wg);
        xw(:,1:(end-1)) = xs(:,2:end);
        xw(:,end) = transitionSample(model,xs(:,end),[a,0]);
        
        % Shift particle histories
        xs(:,1:(nLag-1)) = xs(:,2:end);        
        % Predict
        xs(:,end) = transitionSample(model,xs(:,end),[a,sigma]);
        
    else
        
        u = rand(k,nStrata) + ones(k,1)*(0:(nStrata-1));
        c = wg' ./ (ones(k,1)*wstr(:)');
        c2 = c .* (ones(k,1)*(wstr(:)/sum(wstr))');
        
        xw(:,1:(end-1)) = xs(:,2:end);
        xw(:,end) = transitionSample(model,xs(:,end),[a,0]);
        
        cdf = cumsum(c(:));
        [~,inds] = histc(u(:),cdf);
        
        % Shift particle histories
        xs(:,1:(end-1)) = xs(inds+1,2:end);
        % Predict
        xs(:,end) = transitionSample(model,xs(inds+1,end),[a,sigma]);
    end
    Wprev = wsmp / sum(wsmp);
    
    % --------------------------------------------------------------------
    % Algorithmically irrelevant output manipulations
    
    % Store results for output    
    if settings.store
        ESS(i) = ess;
        NSTRATA(i) = nStrata;
        W(:,i) = Wprev;
        ESSCurves(:,i) = esscurve;
        NSTRATA(i) = nStrata;
        X(:,i,:) = permute(xs,[1 3 2]); % X is N x nStep x nLag
    end
    
    % Here MU{n}(end,k) is the sample particle approximation of function n
    % at time i + 2 - k. That is, MU{n}(1,i) is the prediction estimate,
    % MU{n}(2,i) is the filtering estimate, MU{n}(3,i) is single step
    % smooted estimate.
    xsfl = fliplr(xs);
    MU{1}(:,i) = (Wprev'*xsfl)'; 
    MU{2}(:,i) = (Wprev'*(xsfl.*xsfl.*xsfl.*xsfl))'; 
    MU{3}(:,i) = (Wprev'*(xsfl>0))'; 
    % Strangely x.*x.*x.*x appears to be faster than x.^4
    
    xwfl = fliplr(xw);
    MUW{1}(:,i) = (c2(:)'*xwfl)'; 
    MUW{2}(:,i) = (c2(:)'*(xwfl.*xwfl.*xwfl.*xwfl))'; 
    MUW{3}(:,i) = (c2(:)'*(xwfl>0))'; 
    
    if ~settings.silent
        waitbar(i/nStep,hwb);
    end
end

if ~settings.silent
    close(hwb);
    if settings.store
        visualiseSingleRun(X,MU,ESS,NSTRATA,settings);
    end
end


