function [ PR ] = Packing_Radius( Xn, sequence )
% function [ PR ] = Packing_Radius( Xn, sequence )
% Computes the packing radius of the n point design Xn (a d*n matrix) 
% sequence is a sorted vector of m integers in {1,...,n}, m<=n
% PR is m row vector of packing radii for designs Xn(:,(1:i)), i in sequence
%   (if sequence = m, PR is just PR(Xn))
%-----
% 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,n]=size(Xn);

LPR=length(sequence);
PR=NaN(1,LPR);
timebar= waitbar(0,'packing radius...');
if LPR>1 && (sequence(end)-sequence(1)+1==LPR)
    % The sequence contains all LPR>1 elements by increment 1
    % Compute PR(i) for first i>=2 
    %       (there is no packing radius for a one point design)
    iSgt1=find(sequence>1);
    iX1=sequence(iSgt1(1));
    PR(iSgt1(1))=min(pdist(Xn(:,1:iX1)'))/2;
    for i=iSgt1(1)+1:LPR
        waitbar(i/LPR,timebar);
        iX=sequence(i);
        iXm1=sequence(i-1);
        PR(i)=min(PR(i-1),min(pdist2(Xn(:,1:iXm1)',Xn(:,iX)'))/2);
    end    
else % The sequence contains a few elements only --> brute force   
    for i=1:LPR
        waitbar(i/LPR,timebar);
        iCR=sequence(i);
        if iCR>1 % no packing radius for a one point design
            PR(i)=min(pdist(Xn(:,1:iCR)'))/2;
        end 
    end
end
close(timebar)
end
