function [ X ] = boundary_phobic_GP( Xcand, nmax, beta )
% function [ X ] = boundary_phobic_GP( Xcand, nmax, beta )
% Xcand (d*Q) : a set of Q candidate points well spread in [0,1]^d
% nmax = size of the design X constructed in [0,1]^d (a subset of Xcand)
% 2 <= nmax << Q (Q can be large as not all Q*(Q-1) interdistances are calculated)
% beta>0 forces design points to remain away from the boundary of the hypercube
%          beta=inf --> standard coffee house (efficiency >= 50% for covering and packing)
%          beta=2   --> corresponds to standard spacing
% see the paper "Incremental space-filling design based 
%   on coverings and spacings: improving upon low discrepancy sequences", by 
%   Nogalez-Gómes, Pronzato and Rendas, submitted to the Journal of Statistical Theory 
%   and Practice (2021).           
% The covering and packing efficiency of each designs X(:,1:n), 2 <= n <= nmax  
%   is >= 0.5/(1+sqrt(d)/beta)
% The first point X(:,1) is always the closest one in Xcand to ones(d,1)/2 (i.e., the 
%   closest to the center of the hypercube)
%-----
% 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.

[d,Q]=size(Xcand);
IndX=zeros(1,nmax);

if isinf(beta)==0
    Xboundary=beta*[Xcand' -Xcand'+1]; % precalculation to accelerate a bit
end    

%------------------
% 1) Initialization
%------------------
Xall=NaN(d,nmax);
% 1st point is closest to the center
Center=ones(d,1)/2; 
D_to_center=pdist2(Center',Xcand');
[~,i1]=min(D_to_center); 
X=Xcand(:,i1); IndX(1)=i1;
Xall(:,1)=X;
D_to_X=pdist2(Xcand',X');  % distances to design points

%--------------------
% 2) Start iterations
%--------------------
timebar= waitbar(0,'Coffee-house...'); 
for n=2:nmax
    waitbar(n/nmax,timebar);
    
    [distmax,~]=max(D_to_X);  
    Imax=find(D_to_X==distmax);
    ibest=Imax(1);
    if isinf(beta)==0
        Dum=[D_to_X Xboundary];        
        % impose absence of repetitions of design points (which may happen
        % when beta < inf
        Ind_free=setdiff(1:Q,IndX);        
        Dum=Dum(Ind_free,:);
        [~,ibest]=max(min(Dum,[],2)); ibest=ibest(1);
        ibest=Ind_free(ibest);
    end
    Xnew=Xcand(:,ibest);
    %
    D_to_X=min(D_to_X,pdist2(Xcand',Xnew')); % update all D_to_X     
    x_last=Xnew; 
    Xall(:,n)=x_last;
    X=Xall(:,1:n);
    IndX(n)=ibest;  
end
close(timebar) 
end
        