
function[DQ,MDQ,Dvalue]=DQ_cplex(F,N,arg1,arg2,arg3)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Quadratic integer program for computing exact DQ-optimal designs   %
% (version that uses CPLEX solver of integer quadratic programming)  %
%                                                                    %
% Based on the paper: Harman R, Filova L: Computing efficient exact  %
% designs of experiments using integer quadratic programming,        %
% Computational Statistics & Data Analysis, Volume 71, pp. 1159–1167 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Input:
% F ... an mxn matrix with columns equal to regressors of the model
% (n is the size of the design space and m is the dimension of parameter) 
% N ... an integer number giving the required number of trials (the size)
%
% arg1, arg2, arg3 ... additional parameters with the following meaning:
%
% if only arg1 is supplied, then it is interpreted as M
% if both arg1, arg2 are supplied, they are interpreted as A,b.
% if all three arg1-3 are supplied, they are interpreted as M,A,b, where:
%
% M ... the mxm information matrix of the D-optimal approximate design
% A ... a kxn matrix of real coefficients
% b ... a kx1 vector of real coefficients
% (that is the constraints Aw<=b are added to w(1)+...+w(n)=1)
%
% Notes: The design w above is normalized, i.e., w(i) means the proportion
% of trials in the design point i. The elements of A,b can be arbitrary
% real numbers, but note that for some A,b the system Aw<=b, w(1)+...+w(n)=1 
% has no solution w.

switch nargin
    case 2
       M=0; A=[]; b=[];
    case 3
       M=arg1; A=[]; b=[];
    case 4
       M=0; A=arg1; b=arg2; 
    case 5
       M=arg1; A=arg2; b=arg3;
    otherwise
       error('quip:IncorrectInput','The number of input arguments should be between 2 and 5.');
end;

dF=size(F);
m=dF(1); % number of parameters
n=dF(2); % size of the (discrete) design space
dA=size(A);
k=dA(1);

% If the D-optimal approximate infomation matrix is missing,
% compute it using sedumi (SDP solver).
if M==0
    M=zeros(m,m);
	w=sdpvar(n,1);
	constraints = [w >= 0, sum(w) == 1, A*w<=b];
    for i=1:n
        M=M+w(i)*F(:,i)*F(:,i)';
    end;
	objective = -logdet(M);
	options = sdpsettings('solver','sedumi','savesolveroutput','1');
	approx=solvesdp(constraints, objective, options);
	appweights=approx.solveroutput.y(1:n); 
	for i=1:n
		M=M+appweights(i)*F(:,i)*F(:,i)';
	end;
    M=double(M);
end;

invM=eye(m)/M;

% Set up the quadratic integer program for the CPLEX solver
f=zeros(n,1);
H=zeros(n,n);

for i=1:n
    f(i)=-2*(F(:,i))'*invM*(F(:,i))/N;
	end;

for i=1:n
    for j=i:n
        H(i,j)=((F(:,i))'*invM*(F(:,j))/N)^2;
        H(j,i)=H(i,j);
    end;
end;

sostype=[]; sosind=[]; soswt=[];
Aineq=A; bineq=N*b; Aeq=ones(1,n); beq=N;
lb=zeros(n,1); ub=N*ones(n,1);

ctype=[];
for i=1:n
    ctype=[ctype 'I'];
end;    

%use the cplex solver to compute the DQ-optimal design
DQ=cplexmiqp(H,f,Aineq,bineq,Aeq,beq,sostype,sosind,soswt,lb,ub,ctype);

DQ=DQ./N;
%compute the DQ-optimal information matrix
MDQ=zeros(m,m);
for i=1:n
    MDQ=MDQ+DQ(i)*F(:,i)*F(:,i)';
end;

Dvalue=(det(MDQ))^(1/m);

end

