
function demo_QMC_2D

NV=2.^(1:18);
f=@(x,y) exp(-x.^2-y.^2); I=(sqrt(pi)*erf(1)/2)^2;
% f=@(x,y) ((x-0.5).^2+(y-0.5).^2).^(1/2); I=3.825978534775538e-01;
trials=1000;
AE_mat=[];

% 0. MC, 1. Sobol, 2. Haltonset, 3. Hammersley
QMC_type=3; 

for k=1:trials
    AE=[];
    for n=NV
        P=QMC_pointset(QMC_type,n,2);

        fP=feval(f,P(:,1),P(:,2));
        In=sum(fP)/n;
        AE(end+1,1)=abs(I-In);
    end
    AE_mat(:,k)=AE;
end

AE_mat_aver=abs(sum(AE_mat,2)/trials);

fprintf('\n \n \t \t * AVERAGE ABSOLUTE ERRORS \n \n');
for k=1:length(NV)
    fprintf('\n \t n: %7.0f ae: %1.5e',NV(k),AE_mat_aver(k));
end

fprintf('\n \n \t \t \t * AVERAGE RATIOS \n \n');
rat=AE_mat_aver(1:end-1)./AE_mat_aver(2:end);
for k=1:length(NV)-1
    fprintf('\n \t n1: %7.0f n2: %7.0f rat: %1.5e',NV(k),NV(k+1),rat(k));
end
fprintf('\n \n');



switch QMC_type
    case 0
        fprintf('\n \t * Pointset: random');
    case 1
        fprintf('\n \t * Pointset: Sobol set');
    case 2
        fprintf('\n \t * Pointset: Halton set');
end

fprintf('\n \n');






function P=QMC_pointset(QMC_type,n,dim)

switch QMC_type

    case 0
        P=rand(n,dim);

    case 1

        Pset = sobolset(dim);
        P=Pset(1:n,:);

    case 2
        Pset = haltonset(dim);
        P=Pset(1:n,:);

    case 3
        P = Hammersley(n,dim); P=P';


end

