% Exchange Algorithm to find the exact optimal design (no replications)
% using A-optimal; for mixed 2- and 3-level control and noise factors

%INPUT
%kC: kC(1) number of 2-level control factors; kC(2) number of 
%3-level qualitative control factors; kC(3) number of 3-level
%quantitative control factors
%kn: kn is the number of 2-level noise factors, noise factors only have
%2-levels.
%n: number of runs
%rho is the correlation parameters defined in Joseph and Delaney (2007)

%OUTPUT
%D: index of the points selected in the design
%Design: returned design
%Design=(noise factors|2-level control| 3-level qualitative control| 3-level
%quantitative control)
%Aopt: the objective value

function [D, Design, Aopt]=EAmixed(kC,kn,n,rho)
p2=kC(1)+kn; % number of 2-level factors
p3=kC(2)+kC(3); % number of 3-level factors
N=2^p2*3^p3;

%U is the full model matrix 
%R is the full correlation matrix
%points is all the candidates
U2=[1,-1;1,1];
% the model matrix for qualitative 3-level factor 
U3=[1,-(1.5)^0.5,0.5^0.5;
    1, 0, -2^0.5;
    1, 1.5^0.5, 0.5^0.5];
% the model matrix for quantitative 3-level factor
R2=[1,0; 
    0,(1-rho)/(1+rho)];
% correlation matrix for 3-level qualitative factor
R31=[1, 0, 0;
    0, (1-rho)/(1+2*rho), 0;
    0, 0, (1-rho)/(1+2*rho)];
% correlation matrix for 3-level quantitative factor
tempvalue=3+4*rho+2*rho^4;
R32=[1, 0, -2^0.5*(rho-rho^4)/tempvalue;
     0, 3*(1-rho^4)/tempvalue, 0;
     -2^0.5*(rho-rho^4)/tempvalue, 0, (rho^4-4*rho+3)/tempvalue];
D2=[-1,1]';
D3=[-1,0,1]';
tempU=1;
tempR=1;
points=[-1,1]';

for i=1:p2
    tempU=kron(U2,tempU);
    tempR=kron(R2,tempR);
end
for i=2:p2
    points=[kron([1,1]',points),kron(D2,ones(size(points,1),1))];
end
for i=1:kC(2)
    tempU=kron(U3,tempU);
    tempR=kron(R31,tempR);
    points=[kron([1,1,1]',points),kron(D3,ones(size(points,1),1))];
end
for i=1:kC(3)
    tempU=kron(U3,tempU);
    tempR=kron(R32,tempR);
    points=[kron([1,1,1]',points),kron(D3,ones(size(points,1),1))];
end
U=tempU;
R=tempR;

%flag: indicates the number of factors in a term (=log2(flag))
flag=1;
wt=zeros(1,N);
for i=1:kn
    flag=kron([1,2],flag);
end
for i=1:kC(1)
    flag=kron([1,1],flag);
end
for i=1:p3
    flag=kron([1,1,1],flag);
end
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,:);