% Example of use of boundary_phobic_GP.m
% produces the plots in the paper "Quasi-uniform designs with optimal and 
% near-optimal uniformity constant", by L. Pronzato and A. Zhigljavsky

% Copyright 2021 CNRS [L. Pronzato]
%
% Redistribution and use in source and binary forms, with or without modification, 
% are permitted provided that the following conditions are met:
% 
% 1. Redistributions of source code must retain the above copyright notice, 
% this list of conditions and the following disclaimer.
% 
% 2. Redistributions in binary form must reproduce the above copyright notice, 
% this list of conditions and the following disclaimer in the documentation and/or 
% other materials provided with the distribution.
% 
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY 
% EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 
% OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT 
% SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, 
% INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED 
% TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR 
% BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN 
% ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH 
% DAMAGE.

rng('shuffle');

d=2;

% Candidate set = Q Sobol' points in [0,1]^d
pS = sobolset(d);
Q=2^22;
Xgrid=(pS(1:Q,:))';
% Complement by the vertices and the center of the cube 
k=2; Complement=((fullfact(k*ones(1,d))-1)/(k-1))';
Xgrid=[Xgrid Complement]; 
Xgrid=[Xgrid 1/2*ones(d,1)];
Xgrid=unique(Xgrid','rows'); Xgrid=Xgrid';
%--------------------------------------------------------------------------
% 1) Pure greedy-packing
beta=inf;
nmax=85;
X_inf = boundary_phobic_GP( Xgrid, nmax, beta );
% compute covering and packing radius and mesh ratio of X_inf(:,1:n), n=1,...,nmax
CR_inf=Covering_Radius(X_inf, 1:nmax, 0, Xgrid, 0  );
PR_inf=Packing_Radius(X_inf, 1:nmax );
MR_inf=CR_inf./PR_inf;

figure(1)
drop=find(MR_inf<1.95);

Ldrop=length(drop);
plot([2 nmax],[2 2],'r-','LineWidth',2)
axis([1 nmax 0.9 2.1])
hold on
for i=1:Ldrop
    plot([drop(i) drop(i)],[2 MR_inf(drop(i))],'r--',...
         'LineWidth',2) 
end
xlabel('$n$','Interpreter','latex','FontSize',20)
ylabel('MR$(X_n)$','Interpreter','latex','FontSize',20)
hold off

figure(2)
plot_design( X_inf(:,1:80), CR_inf(80) )

%--------------------------------------------------------------------------
% 2) Greedy-packing with boundary avoidance
beta=4;
nmax=80;
X_4 = boundary_phobic_GP( Xgrid, nmax, beta );
% compute covering and packing radius and mesh ratio of X_inf(:,1:n), n=1,...,nmax
CR_4=Covering_Radius(X_4, 1:nmax, 0, Xgrid, 0  );
PR_4=Packing_Radius(X_4, 1:nmax );
MR_4=CR_4./PR_4;

figure(3)
MR_max=2*(1+sqrt(d)/beta);
plot(1:nmax,MR_4,'r-',...
    [1 nmax],[MR_max MR_max],'k--',...
     'LineWidth',2)
figure(4)
plot_design( X_4(:,1:nmax), CR_4(nmax) )

