function [ CR ] = Covering_Radius( Xn, sequence, d_threshold, Xtest, k_fracfact )
% function [ CR ] = Covering_Radius( Xn, sequence, d_threshold, Xtest, k_fracfact )
% Computes the covering radius of the designs Xm in [0,1]^d 
% for m in sequence = a vector of m consecutive integers in {1,...,n}
% CR is m row vector of covering radii for designs Xn(:,(1:i)), i in sequence
%   (if sequence = n, CR is just CR(Xn))
% for d<=d_threshold, the exact value is calculated (via Voronoi diagram)
% for d>d_threshold, CR is underestimated by evaluation on the Q grid points in Xtest 
% Xtest is a d*Q matrix, for instance of scrambled Sobol' points:
%   pSmm = sobolset(d); pSmm = scramble(pSmm,'MatousekAffineOwen');
%   Xtest=(net(pS,2^19))';
% To keep the computational cost reasonable, d_threshold should be <= 5 (this is enforced)
% When k_fracfact=k>0, Xtest is completed by a k^d fractional factorial design 
%      (==> k should be kept small form large d, k=2 is reasonable)
%-----
% 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,~]=size(Xn);
d_threshold=min(d_threshold,5);
%
nmax=sequence(end);
[~,n]=size(Xn);
if nmax>n
    display('sequence is too long for the nb. of design points Xn') 
    CR=NaN(size(sequence));
    return
end

if isempty(Xtest)==0 && k_fracfact>0
    % complete by a k^d fractional factorial design
    k=k_fracfact;
    Complement=((fullfact(k*ones(1,d))-1)/(k-1))';
    Xtest=[Xtest Complement];
end
LCR=length(sequence);
CR=NaN(1,LCR);
timebar= waitbar(0,'covering radius...'); 
if d>d_threshold
    % Compute all what we need once for all
    Dist_all=pdist2(Xn(:,1:nmax)',Xtest');
end    
for i=1:LCR
    waitbar(i/LCR,timebar);
    iCR=sequence(i);
    if d==1
        CR(i) = Minimax_d1( Xn(:,1:iCR), 0,1 );
    elseif d<=d_threshold
        [ CR(i), ~, ~, ~ ] = voronoi_truncated01( Xn(:,1:iCR), [], [], 0 );
    else
        %CR(i)=max(min(Dist_all(1:iCR,:),[],1));    % --> faster to compute
                                                    %     recursively as below
        if i==1
            D2Xn=min(Dist_all(1:iCR,:),[],1);
%             [CR(i),jmax]=max(D2Xn);
        else
            D2Xn=min(D2Xn,Dist_all(iCR,:));
%             jbig=find(D2Xn>=D2Xn(jmax));
%             D2Xnbig=D2Xn(D2Xn>=D2Xn(jmax));
%             [CR(i),jbigmax]=max(D2Xnbig);
%             jmax=jbig(jbigmax);
        end    
        CR(i)=max(D2Xn);
    end
end
close(timebar)
end
