iterCntThreshold = 5000;
minErrThreshold = 200;
msgL = 2000;
snrdb = 1:5;

errs  = zeros(length(snrdb),1);
iters = zeros(length(snrdb),1);

berplot = cell(1,5);
numframes = 500;    %GPU version runs 500 frames in parallel.

viterbiTransceiverCPU(-10,1,1);
viterbiTransceiverGPU(-10,1,1,numframes);

N=1; %N tracks which simulation variant is run

fprintf(1,'Bit Error Rate Acceleration Analysis Example\n\n');

fprintf(1,'***Baseline - Standard System object simulation***\n');

% create random stream for each snrdb simulation
s = RandStream.create('mrg32k3a','NumStreams',1,...
    'CellOutput',true,'NormalTransform', 'Inversion');
RandStream.setGlobalStream(s{1});

ts = tic;
for ii=1:numel(snrdb)
    fprintf(1,'Iteration number %d, SNR (dB) = %d\n',ii, snrdb(ii));
    [errs(ii),iters(ii)] =viterbiTransceiverCPU(snrdb(ii), minErrThreshold, iterCntThreshold);
end
ber = errs./ (msgL* iters);
baseTime=toc(ts);
berplot{N} = ber;
desc{N} = 'baseline';
reportResultsCommSysGPU(N, baseTime,baseTime, 'Baseline');

fprintf(1,'\n***Baseline + codegen***\n');
N=N+1; %Increase simulation counter

% Create the coder object and turn off checks which will cause low
% performance.
fprintf(1,'Generating Code ...');
config_obj = coder.config('MEX');
config_obj.EnableDebugging = false;
config_obj.IntegrityChecks = false;
config_obj.ResponsivenessChecks = false;
config_obj.EchoExpressions = false;

% Generate a MEX file
codegen('viterbiTransceiverCPU.m', '-config', 'config_obj', '-args', {snrdb(1), minErrThreshold, iterCntThreshold} )
fprintf(1,'  Done.\n');

%Run once to eliminate startup overhead.
viterbiTransceiverCPU_mex(-10,1,1);

s = RandStream.getGlobalStream;
reset(s);

% Use the generated MEX function viterbiTransceiverCPU_mex in the
% simulation loop.
ts = tic;
for ii=1:numel(snrdb)
    fprintf(1,'Iteration number %d, SNR (dB) = %d\n',ii, snrdb(ii));
    [errs(ii),iters(ii)] = viterbiTransceiverCPU_mex(snrdb(ii), minErrThreshold, iterCntThreshold);
end
ber = errs./ (msgL* iters);
trialtime=toc(ts);
berplot{N} = ber;
desc{N} = 'codegen';
reportResultsCommSysGPU(N, trialtime,baseTime, 'Baseline + codegen');

fprintf(1,'\n***Baseline + parfor***\n');
fprintf(1,'Accessing multiple CPU cores ...\n');
if isempty(gcp('nocreate'))
    pool = parpool;
    poolWasOpen = false;
else
    pool = gcp;
    poolWasOpen = true;
end
nW=pool.NumWorkers;
N=N+1; %Increase simulation counter

snrN = numel(snrdb);

mT = minErrThreshold / nW;
iT = iterCntThreshold / nW;

errN =  zeros(nW, snrN);
itrN =  zeros(nW, snrN);

% replicate snrdb
snrdb_rep=repmat(snrdb,nW,1);

% create an independent stream for each worker
s = RandStream.create('mrg32k3a','NumStreams',nW,...
    'CellOutput',true,'NormalTransform', 'Inversion');

% pre-run
parfor jj=1:nW
    RandStream.setGlobalStream(s{jj});
    viterbiTransceiverCPU(-10, 1, 1);
end

fprintf(1,'Start parfor job ... ');
ts = tic;
parfor jj=1:nW
    for ii=1:snrN
        [err, itr] = viterbiTransceiverCPU(snrdb_rep(jj,ii), mT, iT);
        errN(jj,ii) = err;
        itrN(jj,ii) = itr;
    end
end
ber = sum(errN)./ (msgL*sum(itrN));
trialtime=toc(ts);
fprintf(1,'Done.\n');
berplot{N} = ber;
desc{N} = 'parfor';
reportResultsCommSysGPU(N, trialtime,baseTime, 'Baseline + parfor');

fprintf(1,'\n***Baseline + codegen + parfor***\n');
N=N+1; %Increase simulation counter

% pre-run
parfor jj=1:nW
    RandStream.setGlobalStream(s{jj});
    viterbiTransceiverCPU_mex(1, 1, 1); % use the same mex file
end

fprintf(1,'Start parfor job ... ');
ts = tic;
parfor jj=1:nW
    for ii=1:snrN
        [err, itr] = viterbiTransceiverCPU_mex(snrdb_rep(jj,ii), mT, iT);
        errN(jj,ii) = err;
        itrN(jj,ii) = itr;
    end
end
ber = sum(errN)./ (msgL*sum(itrN));
trialtime=toc(ts);
fprintf(1,'Done.\n');
berplot{N} = ber;
desc{N} = 'codegen + parfor';
reportResultsCommSysGPU(N, trialtime,baseTime, 'Baseline + codegen + parfor');

fprintf(1,'\n***GPU***\n');
N=N+1; %Increase simulation counter

try
    dev = parallel.gpu.GPUDevice.current;
    fprintf(...
        'GPU detected (%s, %d multiprocessors, Compute Capability %s)\n',...
        dev.Name, dev.MultiprocessorCount, dev.ComputeCapability);

    sg = parallel.gpu.RandStream.create('mrg32k3a','NumStreams',1,'NormalTransform','Inversion');
    parallel.gpu.RandStream.setGlobalStream(sg);

    ts = tic;
    for ii=1:numel(snrdb)
        fprintf(1,'Iteration number %d, SNR (dB) = %d\n',ii, snrdb(ii));
        [errs(ii),iters(ii)] =viterbiTransceiverGPU(snrdb(ii), minErrThreshold, iterCntThreshold, numframes);
    end
    ber = errs./ (msgL* iters);
    trialtime=toc(ts);
    berplot{N} = ber;
    desc{N} = 'GPU';
    reportResultsCommSysGPU(N, trialtime,baseTime, 'Baseline + GPU');

    fprintf(1,'  Done.\n');

catch %#ok<CTCH>

    % Report that the appropriate GPU was not found.
    fprintf(1, ['Could not find an appropriate GPU or could not ', ...
        'execute GPU code.\n']);

end

lines = {'kx-.', 'ro-', 'cs--', 'm^:', 'g*-'};
for ii=1:numel(desc)
    semilogy(snrdb, berplot{ii}, lines{ii});
    hold on;
end
hold off;
title('Bit Error Rate for Various Acceleration Strategies');
xlabel('Signal to Noise Ratio (dB)');
ylabel('BER');
legend(desc{:});

if ~poolWasOpen
    delete(gcp);
end