%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%    This file replicate the figure 1 of the paper.
%    You need to set the number of iterations  
%    N.B: The program is not optimized for GPU we where not iterested in    
%    timinig in this exercise.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gav = zeros(n, 3, iIter);
gavvar = NaN(n, 3, iIter);
gavGPU = gpuArray.zeros(n, 3, iIter);
gavvarGPU = gpuArray(NaN(n, 3, iIter));
begin = 10;
if fixgenrator==1
    sc = RandStream('CombRecursive', 'Seed', 1);
    RandStream.setGlobalStream(sc);
    sg = parallel.gpu.RandStream('CombRecursive', 'Seed', 1);
    parallel.gpu.RandStream.setGlobalStream(sg);
end
for j = 1:iIter
    x = rand(n,1);
    gav(1, 1, j) = x(1, 1);
    gav(1, 2, j) = x(1, 1)^2;
    gav(1, 3, j) = cos(pi * x(1, 1));
    for i = 2:n
        gav(i, 1, j) = sum(x(1:i))/i;
        gav(i, 2, j) = sum(x(1:i).^2)/i;
        gav(i, 3, j) = sum(cos(pi * x(1:i)))/i;
        gavvar(i, 1, j) = var(x(1:i));
        gavvar(i, 2, j) = var(x(1:i).^2);
        gavvar(i, 3, j) = var(cos(pi * x(1:i)));
    end
    %%%%%%%%% GPU implementation %%%%%%%%%%%%%%%%%%%
    xGPU = gpuArray.rand(n, 1);
    gavGPU(1, 1, j) = xGPU(1, 1);
    gavGPU(1, 2, j) = xGPU(1, 1)^2;
    gavGPU(1, 3, j) = cos(pi * xGPU(1, 1));
    for i = 2:n
        gavGPU(i, 1, j) = sum(xGPU(1:i))/i;
        gavGPU(i, 2, j) = sum(xGPU(1:i).^2)/i;
        gavGPU(i, 3, j) = sum(cos(pi * xGPU(1:i)))/i;
        gavvarGPU(i, 1, j) = var(xGPU(1:i));
        gavvarGPU(i, 2, j) = var(xGPU(1:i).^2);
        gavvarGPU(i, 3, j) = var(cos(pi * xGPU(1:i)));
    end
end
MeanGPU = gather(gavGPU);
VarianceGPU = gather(gavvarGPU);
MeanCPU = gav;
VarianceCPU = gavvar;
%%%%%%%%%%%%%%%%%%%%% Take the know quantities %%%%%%%%%%%%%%%%%%%%%%%%%%%%
mTheoricalMean1 =  kron(ones(n - begin + 1, 1)/2, ones(1, iIter));
mTheoricalMean2 =  kron(ones(n - begin + 1, 1)/3, ones(1, iIter));
mTheoricalMean3 =  kron(ones(n - begin + 1, 1) * 0, ones(1, iIter));
mTheoricalVar1 =  kron(ones(n - begin + 1, 1)/12, ones(1, iIter));
mTheoricalVar2 =  kron(ones(n - begin + 1, 1)*4/45, ones(1, iIter));
mTheoricalVar3 =  kron(ones(n - begin + 1, 1)*1/2, ones(1, iIter));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       GPU
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% Distance Mean GPU %%%%%%%%%%%%%%%%%%%
mDistanceMeanGPU1 = mean((squeeze(MeanGPU(begin:end, 1, :)) - mTheoricalMean1).^2);
mDistanceMeanGPU2 = mean((squeeze(MeanGPU(begin:end, 2, :)) - mTheoricalMean2).^2);
mDistanceMeanGPU3 = mean((squeeze(MeanGPU(begin:end, 3, :)) - mTheoricalMean3).^2);
%%%%%%%%%%%%%%%%% Distance Variance GPU %%%%%%%%%%%%%%%%%%%
mDistanceVarGPU1 = mean((squeeze(VarianceGPU(begin:end, 1, :)) - mTheoricalVar1).^2);
mDistanceVarGPU2 = mean((squeeze(VarianceGPU(begin:end, 2, :)) - mTheoricalVar2).^2);
mDistanceVarGPU3 = mean((squeeze(VarianceGPU(begin:end, 3, :)) - mTheoricalVar3).^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       CPU
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%% Distance Mean CPU %%%%%%%%%%%%%%%%%%%
mDistanceMeanCPU1 = mean((squeeze(MeanCPU(begin:end, 1, :)) - mTheoricalMean1).^2);
mDistanceMeanCPU2 = mean((squeeze(MeanCPU(begin:end, 2, :)) - mTheoricalMean2).^2);
mDistanceMeanCPU3 = mean((squeeze(MeanCPU(begin:end, 3, :)) - mTheoricalMean3).^2);
%%%%%%%%%%%%%%%%% Distance Variance CPU %%%%%%%%%%%%%%%%%%%
mDistanceVarCPU1 = mean((squeeze(VarianceCPU(begin:end, 1, :)) - mTheoricalVar1).^2);
mDistanceVarCPU2 = mean((squeeze(VarianceCPU(begin:end, 2, :)) - mTheoricalVar2).^2);
mDistanceVarCPU3 = mean((squeeze(VarianceCPU(begin:end, 3, :)) - mTheoricalVar3).^2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   Make the graph for figure 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fs=2;
figure(1);
subplot(3, 1, 1);
[fCPU,xCPU]=ecdf(mDistanceMeanCPU1');
[fGPU,xGPU]=ecdf(mDistanceMeanGPU1');
plot(xGPU,fGPU,'r','Linewidth',fs);hold on
plot(xCPU,fCPU,'b','Linewidth',fs);hold off 
legend('GPU','CPU');
subplot(3, 1, 2);
[fCPU,xCPU]=ecdf(mDistanceMeanCPU2');
[fGPU,xGPU]=ecdf(mDistanceMeanGPU2');
plot(xGPU,fGPU,'r','Linewidth',fs);hold on
plot(xCPU,fCPU,'b','Linewidth',fs);hold off 
legend('GPU','CPU');
subplot(3, 1, 3);
[fCPU,xCPU]=ecdf(mDistanceMeanCPU3');
[fGPU,xGPU]=ecdf(mDistanceMeanGPU3');
plot(xGPU,fGPU,'r','Linewidth',fs);hold on
plot(xCPU,fCPU,'b','Linewidth',fs);hold off 
legend('GPU','CPU');
if fixgenrator==1
    figname='output/Figure1meanFixed.pdf';
else
    figname='output/Figure1mean.pdf';
end
saveas(gcf, figname);


fs=2;
figure(1);
subplot(3, 1, 1);
[fCPU,xCPU]=ecdf(mDistanceVarCPU1');
[fGPU,xGPU]=ecdf(mDistanceVarGPU1');
plot(xGPU,fGPU,'r','Linewidth',fs);hold on
plot(xCPU,fCPU,'b','Linewidth',fs);hold off 
legend('GPU','CPU');
subplot(3, 1, 2);
[fCPU,xCPU]=ecdf(mDistanceVarCPU2');
[fGPU,xGPU]=ecdf(mDistanceVarGPU2');
plot(xGPU,fGPU,'r','Linewidth',fs);hold on
plot(xCPU,fCPU,'b','Linewidth',fs);hold off 
legend('GPU','CPU');
subplot(3, 1, 3);
[fCPU,xCPU]=ecdf(mDistanceVarCPU3');
[fGPU,xGPU]=ecdf(mDistanceVarGPU3');
plot(xGPU,fGPU,'r','Linewidth',fs);hold on
plot(xCPU,fCPU,'b','Linewidth',fs);hold off 
legend('GPU','CPU');
if fixgenrator==1
    figname='output/Figure1varFixed.pdf';
else
    figname='output/Figure1var.pdf';
end
saveas(gcf, figname);

% figure(1);
% subplot(3, 1, 1);
% hist(mDistanceMeanCPU1 - mDistanceMeanGPU1, 75); 
% ecdf(mDistanceMeanCPU1 - mDistanceMeanGPU1); 
% ecdf(mDistanceMeanCPU1');hold on 
% ecdf(mDistanceMeanGPU1');
% 
% h = findobj(gca,'Type','patch');
% set(h,'FaceColor','r','EdgeColor','w','facealpha',0.75);
% subplot(3, 1, 2);
% hist(mDistanceMeanCPU2 - mDistanceMeanGPU2, 75);
% ecdf(mDistanceMeanCPU2');hold on 
% ecdf(mDistanceMeanGPU2');
% 
% h = findobj(gca,'Type','patch');
% set(h,'FaceColor','r','EdgeColor','w','facealpha',0.75);
% subplot(3, 1, 3);
% hist(mDistanceMeanCPU3 - mDistanceMeanGPU3, 75); 
% h = findobj(gca,'Type','patch');
% set(h,'FaceColor','r','EdgeColor','w','facealpha',0.75);
% figname='output/Figure1mean.pdf';
% saveas(gcf, figname);
% 
% 
% figure(2);
% subplot(3, 1, 1);
% hist(mDistanceVarCPU1 - mDistanceVarGPU1, 75); 
% h = findobj(gca,'Type','patch');
% set(h,'FaceColor','r','EdgeColor','w','facealpha',0.75);
% subplot(3, 1, 2);
% hist(mDistanceVarCPU2 - mDistanceVarGPU2, 75); 
% h = findobj(gca,'Type','patch');
% set(h,'FaceColor','r','EdgeColor','w','facealpha',0.75);
% subplot(3, 1, 3);
% hist(mDistanceVarCPU3 - mDistanceVarGPU3, 75); 
% h = findobj(gca,'Type','patch');
% set(h,'FaceColor','r','EdgeColor','w','facealpha',0.75);
% figname='output/Figure1var.pdf';
% saveas(gcf, figname);
