clear all
close all
settings = createSettings();

% ------------------------------------------------------------------------
% Parameters
a = .9;
b = .1;
sigma = .25;
nStep = 1000;

% ------------------------------------------------------------------------
% Reference run
[x,y] = createData(nStep,a,b,sigma);

settings.param = [a,b,sigma];
settings.nStep = nStep;
settings.x = x;
settings.y = y;
settings.method = 'bootstrap';
settings.essThreshold = .5;
settings.nStep = nStep;
settings.N = 2^17;
settings.param = [a,b,sigma];
settings.nLag = 8;
settings.silent = false;
settings.store = false;

fprintf('Settings:\n');
disp(settings);
    
[~,refMU,refMUW] = alphaSMC(settings);
 
% ------------------------------------------------------------------------
% Test runs
NMC = 1000;
settings.N = 2^9;
settings.silent = true;
settings.resetFigs = false;
essThresholds = .05:.05:1;
nTH = length(essThresholds);
settings.store = false;

methodLabels = {'adaptiveBootstrap','blockstrategy1','greedy','SimpleRearrange'};
progressLabels = {'Running adaptive bootstrap',...
    'Running alpha SMC with block strategy',...
    'Running alpha SMC with greedy strategy',...
    'Running alpha SMC with simple strategy with rearrangement'};
nMethod = length(methodLabels);

fncLabels = {'$\phi(x) = x$','$\phi(x) = x^4$','\phi(x) = 1_{x>0}$'};
fncYRanges = {[0 .04],[0 .3],[0, .3]};
nFnc = length(refMU);

burnIn = 30;

results = zeros(nMethod,nFnc,settings.nLag,nTH,NMC);
wresults = zeros(nMethod,nFnc,settings.nLag,nTH,NMC);

for meth = 1:nMethod
    
    fprintf('%s\n',methodLabels{meth});
    settings.method = methodLabels{meth};
    
    fprintf('Settings:\n');
    disp(settings);
    
    for th = 1:nTH
        
        fprintf('\tESS threshold %6.2f ',essThresholds(th));
        settings.essThreshold = essThresholds(th);
        divisor = 0;
        
        fprintf('*****');
        for k = 1:NMC
            
            fprintf('\b\b\b\b\b%5i',k);
            
            [X,MU,MUW,ESS,NSTRATA,W] = alphaSMC(settings);
            
            divisor = divisor + settings.nStep - burnIn + 1;            
            
            for lag = 1:settings.nLag
            
                for fnc = 1:nFnc
                    
                    results(meth, fnc, lag, th, k) = ...
                        results(meth, fnc, lag, th, k) + ...
                        sum((MU{fnc}(settings.nLag-lag+1,burnIn:end)' - ...
                        refMU{fnc}(settings.nLag-lag+1,burnIn:end)').^2);
                    
                    wresults(meth, fnc, lag, th, k) = ...
                        wresults(meth, fnc, lag, th, k) + ...
                        sum((MUW{fnc}(settings.nLag-lag+1,burnIn:end)' - ...
                        refMUW{fnc}(settings.nLag-lag+1,burnIn:end)').^2);
                                        
                end
            end
        end
        
        fprintf('Done\n');
    end
end

save(['largerun' datestr(now,'yyyy-mm-dd-HH-MM-SS') '.mat']);
