% Exchange Algorithm to find the exact optimal design (no replications)
% For 2-level control and 2-level noise factors only
% INPUT
% kC: number of control factors;
% kn: number of noise variables 
% n: number of runs
% r: parameter of the prior correlation matrix of beta (usually set to be
% 1/3)

% OUTPUT
% D: index of the points selected in the design 
% Design: design matrix=(control factors | noise factors)
% Aopt: the objective value

function [D,Design, Aopt]=EA2level(kC,kn,n,r)
p=kC+kn; %p is the number of factors
N=2^p;   %N is the number of candidates

%U is the full model matrix for all the candidate points
%points is the design matrix (all the candidates points)
%R is the correlation matrix
%points is the design settings for all the factors. 
U=[1,-1;1,1];
R=diag([1,r]);
points=[-1,1]';
for i=2:p
    U=[U,-U;U,U];
    points=[kron([1,1]',points),kron([-1,1]',ones(size(points,1),1))];
    temp=diag(R)';
    R=diag([temp,r*temp]);
end

%flag is used to indicate how many noise factors in the term =log2(flag)
wt=zeros(1,N);
flag=1;
for i=1:kC
    flag=kron([1,1],flag);
end
for i=1:kn
    flag=kron([1,2],flag);
end
% wt is the diagonal matrix A, wt is user specified
wt(flag==2)=1;

% pick up the initial points
n0=floor(n/3);
pick=0;
D=myrandint(1,n0,1:N,'noreplace');   % introduce some randomness to avoid trapping in the local optimal
temp=U(D,:)*R;
V=temp*U(D,:)';  % as long as D has no replicates, V is not singular
M=temp'*inv(V)*temp;
opt=wt*diag(M);
left=N-n0;
candidate=setdiff(1:N,D); % all the left candidate points
    
%choose the next n-n0 points
l=n0+1;
while l<=n
   F=R-M;
   delta0=0;
   i=1;
   while i<=left
       f=U(candidate(i),:);
       a=f*F;
       d=a*f';
       delta1=wt*(a.*a)'/d;
       if delta1>delta0
           a0=a;
           d0=d;
           pick=candidate(i);
           delta0=delta1;
       end
       i=i+1;
   end
  M=M+(a0'*a0)/d0;
  D=[D,pick];
  candidate=candidate(candidate~=pick);
  l=l+1;
  opt=opt+delta0;
  left=left-1;
end

while 1
    compopt=opt;
    for i=1:n
        pick=0;
        D1=setdiff(D,D(i));
        temp=U(D1,:)*R;
        V=inv(temp*U(D1,:)');
        M1=temp'*V*temp;
        opt1=wt*diag(M1);
        F=R-M1;
        delta0=opt-opt1;
        for j=1:left
            f=U(candidate(j),:);
            a=f*F;
            d=a*f';
            delta1=wt*(a.*a)'/d;
            if delta1>delta0
                pick=candidate(j);
                delta0=delta1;
            end
        end  
        if pick~=0
           opt=opt1+delta0;
           candidate=candidate(candidate~=pick);
           candidate=[candidate,D(i)];
           D(i)=pick;
        end
    end
    if abs(compopt-opt)<1e-8
        break
    end
end
D=sort(D);
Aopt=opt/(wt*diag(R));
Design=points(D,:);