function[Dopt,no_iter,xi]=Henry(d,L,N,del,mult_ini,add_rnd,order_type,optA,optB,xi_ini,maxt)
% Optimization for "Bridge" constraints on the grid {1,...,L}^d and 
% D-optimality criterion det(F'F/N)^(1/m) with the constraint on time
% Note that "Incomplete LHDs" are a special case of BDs for del=1.
%
% INPUT:
% d ... the number of factors
% L ... the number of levels (the same for each factor)
% N ... required size of the design
% del ... span of each privacy set (in each dimension, the distance of
%         the points must be at least del)
%
% mult_ini ... do we perform a multiple random selection of the best
%              permissible point to be added to the current design
%              (which is then used to initiate the local search)?
% add_rnd ... do we double the size of the set A of "candidate addition
%             points" by random mutation points (one random point is added
%             after every star point)?
% order_type ... 0 = star points in increasing distance from
%                the design points (with possibly interspersed random pts)
%            ... 1 = completely randomized order of all pts in A
%            ... 2 = all pts in A sorted by the variance function
% optA ... do we optimize in the initial greedy phase?
% optB ... do we optimize during the "exchanges" phase?
%
% xi_ini ... initial design (if contains 0 rows, the intial design is
%            constructed using the greedy forward method).
% maxt ... maximum allowed time of computation
%
% IMPLEMENTATION DETAILS:
% Each design point is represented by a row-vector (of the type 1xd).
% Each design xi is represented by an Nxd matrix of design points.
% That is, xi(i,j) is the index of the j-th level of the i-th design point.
% With xi, we always modify its Lxd matrix P_xi of boolean values
% representing permissible levels. That is, P_xi(i,j)=0 if the j-th level
% of dimension i is already "occupied".
% Search for a best design is performed using the variance function, rather
% than determinants of full designs.
%
% OBSERVED BEHAVIOUR (requires a more careful analysis):
% mult_ini ... setting mult_ini=1 seems to be very slightly better.
% add_rnd ... setting add_rnd=0 leads to significantly better results.
% order_type ... it seems that rand_order=0 is very slightly better than 1
%                which is significantly better than 2. 
% optA ... the value optA=1 seems to very slightly improve average results.
%          However, optA=0 gives a greater variance, i.e., for  multiple
%          restarts setting optA=0 could be better.
% optB ... setting optB=1 leads to marginally longer computation time but
%          greatly improves the quality of the results. It is clear that
%          optB=1 is the best.
%
% POSSIBLE FURTHER IMPROVEMENTS:
% 1) Update the criterion (perhaps also M or its inverse) together with the 
% design, using the fact that we always change only one design point. (The 
% program still contains some redundant computations.)
% 2) It might be advantageous to mutate only two or three coordinates, not
% just "one coordinate" (the star pts) or "all coordinates" (if add_rnd=1).

%display(datestr(now,'dd-mm-yyyy HH:MM:SS FFF'));
eps=0.000001; % Against small numerical errors.
tic;

%Full quadratic model
  m=1+2*d+nchoosek(d,2);
  function[reg]=f(x)
     reg=[1;x;x.^2;prod(nchoosek(x,2),2)];
  end

%    Linear model
%     m=1+d;
%     function[reg]=f(x)
%         reg=[1;x];
%     end

    function[val]=Phi_norm(xi)
        % Value of the normalized xi on [-1,1].
        N_xi=size(xi,1);
        F=zeros(N_xi,m);
        lvl=(0:(L-1))/((L-1)/2)-1;
        for i=1:N_xi
            F(i,:)=f(lvl(xi(i,:))')';
        end
        val=(det(F'*F/N))^(1/m);
    end

    % Check whether points x and y collide.
    function[col]=Collides(x,y)
        if min(abs(x-y))<del
            col=1;
        else
            col=0;
        end
    end

   % Remove the k-th design point from a permissible xi.
    function[zeta,P_zeta]=Remove_point(xi,P_xi,k)
        zeta=xi;
        P_zeta=P_xi;
        zetak=zeta(k,:);
        zeta(k,:)=[];
        for j=1:d
            zeta_j=zeta(:,j);
            smaller=(zeta_j<zetak(j));
            if sum(smaller)==0
                i_left=1;
            else
                i_left=max(zeta_j(smaller))+del;
            end
            greater=(zeta_j>zetak(j));
            if sum(greater)==0
                i_right=L;
            else
                i_right=min(zeta_j(greater))-del;
            end
            P_zeta(i_left:i_right,j)=1;
        end
    end
       
    % Add a permissible point x to a permissible design xi.
    function[zeta,P_zeta]=Add_point(xi,P_xi,x)
        zeta=xi;
        P_zeta=P_xi;
        N_zeta=size(zeta,1);
        zeta(N_zeta+1,:)=x;
        for j=1:d
            min0=max(x(j)-del+1,1);
            max0=min(x(j)+del-1,L);
            P_zeta(min0:max0,j)=0;
        end
    end
 
    % Generate a random permissible x to be included in permissible xi.
    function[x]=Generate_point(P_xi)
        x=zeros(0,d);
        for j=1:d
            levels=1:L;
            levels=levels(logical(P_xi(:,j)));
            if size(levels,2)>1
                x(j)=randsample(levels,1);
            else
                x(j)=levels(1);
            end
        end
    end

    % Optimize the position of the k-th design point in a permissible xi.
    function[zeta,P_zeta]=Optimize(xi,P_xi,k)
   
        % Compute xi and P obtained after removing the k-th point from xi.
        % (P_xi will be updated at the end of the local optimization.)
        [xi0,P]=Remove_point(xi,P_xi,k);
        
        % Compute the inverse of M(xi0) to be used in the var-function.
        N_xi0=size(xi0,1);
        F=zeros(N_xi0,m);
        for i=1:N_xi0
            F(i,:)=f(xi0(i,:)')';
        end
        M1=eye(m)/(F'*F+eps*eye(m));    
        
        % Pos is an Lxd matrix of permissible numbers of increases.
        Pos=zeros(L,d);
        for j=1:d
            s=0;
            for i=L:(-1):1
                s=s+1;
                if P(i,j)
                    Pos(i,j)=s;
                    s=0;
                end
            end
        end

        % Neg is an Lxd matrix of permissible numbers of decreases.
        Neg=zeros(L,d);
        for j=1:d
            s=0;
            for i=1:L
                s=s+1;
                if P(i,j)
                    Neg(i,j)=s;
                    s=0;
                end
            end
        end

        % During the optimization, we will move the k-th point of zeta.
        zeta=xi;
        var_best=f(zeta(k,:)')'*M1*f(zeta(k,:)');
        
        % Main body of the local optimization.
        l_opt=0;
        while ~l_opt
            l_opt=1;
            for j=1:d
                original=zeta(k,j);
                % Attempt to increase the level of the j-th coordinate.
                candidate=zeta(k,j)+Pos(zeta(k,j),j);
                if candidate<=L
                    zeta(k,j)=candidate;
                    var_act=f(zeta(k,:)')'*M1*f(zeta(k,:)');
                    if var_act>var_best
                        var_best=var_act;
                        l_opt=0;
                        break;
                    else
                        zeta(k,j)=original;
                    end
                end
                % Attempt to decrease the level of the j-th coordinate.
                candidate=zeta(k,j)-Neg(zeta(k,j),j);
                if candidate>=1
                    zeta(k,j)=candidate;
                    var_act=f(zeta(k,:)')'*M1*f(zeta(k,:)');
                    if var_act>var_best
                        var_best=var_act;
                        l_opt=0;
                        break;
                    else
                        zeta(k,j)=original;
                    end
                end
            end
        end
   
        % Update the matrix P_zeta.
        P_zeta=P;
        for j=1:d
            min0=max(zeta(k,j)-del+1,1);
            max0=min(zeta(k,j)+del-1,L);
            P_zeta(min0:max0,j)=0;
        end
    end

    % Maximize permissible xi in a greedy way (xi can be empty).
    function[zeta,P_zeta]=Maximize(xi,P_xi,do_optimize)

        zeta=xi;
        P_zeta=P_xi;
        N_zeta=size(zeta,1);
        
        for t=(N_zeta+1):N
            x_bst=Generate_point(P_zeta);
            if mult_ini==1
                if N_zeta>0
                    F_tmp=zeros(N_zeta,m);
                    for i_tmp=1:N_zeta
                        F_tmp(i_tmp,:)=f(zeta(i_tmp,:)')';
                    end
                    M1=eye(m)/(F_tmp'*F_tmp+eps*eye(m));
                else
                    M1=zeros(m,m);
                end
                var_bst=f(x_bst')'*M1*f(x_bst');
                for ix=2:(mult_ini*round(2.5^d))
                    x=Generate_point(P_zeta);
                    var_xact=f(x')'*M1*f(x');
                    if var_xact>var_bst
                        x_bst=x; var_bst=var_xact;
                    end
                end
            end
            [zeta,P_zeta]=Add_point(zeta,P_zeta,x_bst);
            N_zeta=N_zeta+1;
            if do_optimize
                [zeta,P_zeta]=Optimize(zeta,P_zeta,N_zeta);
            end
        end
   
    end

    % Generate the zeta-mutation of a maximal permissible xi by adding x.
    function[zeta,P_zeta]=Mutate(xi,x,opt)
  
        % Create zeta as xi without the point colliding with x.
        [zeta,P_zeta]=Add_point(zeros(0,d),ones(L,d),x);
        N_zeta=1; % This could be slightly improved by a top-down method.
        for i=1:N
            if ~Collides(xi(i,:),x)
                [zeta,P_zeta]=Add_point(zeta,P_zeta,xi(i,:));
                N_zeta=N_zeta+1;
            end
        end
  
        if N_zeta==N+1
            % If no point had been removed, delete the "worst" point.
            zeta_bt=zeta; zeta_bt(1,:)=[];
            Phi_bt=Phi_norm(zeta_bt);
            for i=2:N
                zeta_ac=zeta; zeta_ac(i,:)=[];
                Phi_ac=Phi_norm(zeta_ac);
                if Phi_ac>Phi_bt
                    Phi_bt=Phi_ac;
                    zeta_bt=zeta_ac;
                end
            end
            zeta=zeta_bt;
        elseif N_zeta<=N-1
            % If more than 1 point had been removed, maximize the design.
            [zeta,P_zeta]=Maximize(zeta,P_zeta,opt);
        end
    end

% Create an initial design.
%tic; 
xi=zeros(0,d); P_xi=ones(L,d);
if size(xi_ini,1)==0
    [xi,~]=Maximize(xi,P_xi,optA);
else
    for id=1:size(xi_ini,1)
        [xi,P_xi]=Add_point(xi,P_xi,xi_ini(id,:));
    end
end
scatter(xi(:,1),xi(:,2),'filled'); drawnow;

% Main cycle. 
finish=0;
Phi_norm_xi=Phi_norm(xi);

no_iter=0;
while (~finish && toc<maxt)
    finish=1;
    no_iter=no_iter+1;

    % Compute the "star" set A(xi) of mutation points.
     if add_rnd
         A=zeros(2*d*(L-1)*N,d);
     else
         A=zeros(d*(L-1)*N,d);
     end
     iA=0; I=eye(d);
     pperm=randsample(N,N);
     for dist=1:(L-1)
         for ip=1:N
             xip=xi(pperm(ip),:);
             for dir=1:d
                 cand_plus=xip+dist*I(dir,:);
                 if max(cand_plus)<=L
                     iA=iA+1;
                     A(iA,:)=cand_plus;
                     if add_rnd
                         iA=iA+1;
                         A(iA,:)=randsample(L,d,1)';
                     end
                 end
                 cand_minus=xip-dist*I(dir,:);
                 if min(cand_minus)>=1
                     iA=iA+1;
                     A(iA,:)=cand_minus;
                     if add_rnd
                         iA=iA+1;
                         A(iA,:)=randsample(L,d,1)';
                     end
                 end
             end
         end
     end
     
     N_mut=size(A,1);
     if order_type==1
         A=A(randsample(N_mut,N_mut),:);
     elseif order_type==2
        F_xi=zeros(N,m);
        for i_xi=1:N
            F_xi(i_xi,:)=f(xi(i_xi,:)')';
        end
        M1=eye(m)/(F_xi'*F_xi+eps*eye(m));
        varf=zeros(N_mut,1);
        for ivar=1:N_mut
            varf(ivar)=f(A(ivar,:)')'*M1*f(A(ivar,:)');
        end
        [~,index]=sort(varf,'descend');
        A=A(index,:);
     end
     
    % A natural idea is to clean A to contain only unique rows, but this
    % seems to be not worth doing, even using the fast function "unique".
  
    % Search the set of mutations; transfer to the first mutation better
    % than the current best.
    %next_mut=1;
    for k=1:N_mut
 
        % Stop if time limit is exceeded.
        if(toc>maxt)
            break
        end
        % Perform mutation and check its value.
        zeta_act=Mutate(xi,A(k,:),optB);
        Phi_act=Phi_norm(zeta_act);
        if Phi_act>Phi_norm_xi
            xi=zeta_act;
            Phi_norm_xi=Phi_act;
            finish=0;  %display(xi);
            scatter(xi(:,1),xi(:,2),'filled'); drawnow;
            break;
        end
      
        % Display information about the current state of computation.
%         if k>=next_mut
%           % display('current time / iteration no. / mutation no. / required mutations');
%           % display([round(toc),no_iter,k,N_mut]);
%            next_mut=next_mut+20;
%         end
 
    end
end

Dopt=Phi_norm_xi;
end

