%%% Bernoulli Factory Simulation
%%% Let p = P(Z>2), Z~N(0,1). Want to simulate Bernoulli(Cp), C>1 
%%% (c) Copyright (2014) Staci White


C=2;
epsilonbf = 0.2;
gammabf=0.465;
m = 2.324;

H=50000;
dev_Bpvec=gpuArray(zeros(H,1));

K=1000; %%% Guess of max number of Ber(p) needed for each draw from Ber(Cp)

Zrv = randn(H,K);
I=find(Zrv>2);
Brv = zeros(H,K);
Brv(I)=1;
dev_Brv=gpuArray(Brv);  %%% Brv contains the draws from Bernoulli(p)

%%% Need uniform(0,1) random variables to transform to Geometric in first loop
%%% Need uniform(0,1) random variables to get Bernoulli draws in second part of loop

dev_unif1=gpuArray.rand(H,K);
dev_unif2=gpuArray.rand(H,100);

%Load the Kernel
kernel=parallel.gpu.CUDAKernel('BF.ptx', 'BF.cu','BF');

%Set up the threads
kernel.ThreadBlockSize=[1,1,1];
kernel.GridSize=[H,1];

tic

[dev_Bpvec] = feval(kernel,C,epsilonbf,gammabf, m, H,dev_Brv, dev_unif1, dev_unif2, dev_Bpvec);

Bpvec=gather(dev_Bpvec); %%% Bpvec contains the Bernoulli(Cp) draws

toc


