function [DQ,MDQ,Dvalue]=DQ_gurobi(F,N,arg1,arg2,arg3)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Quadratic integer program for computing exact DQ-optimal designs   %
% (version that uses Gurobi 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 integer quadratic problem for Gurobi

c=zeros(n,1);
Q=zeros(n,n);

for i=1:n
    c(i)=-2*(F(:,i))'*invM*(F(:,i))/N;
end;
	
model.obj=c;

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

l=1;
for i=1:n
    for j=1:n
        if(Q(i,j)~=0)
            opts.QP.qrow(l)=i;
            opts.QP.qcol(l)=j;
            opts.QP.qval(l)=Q(i,j);
            l=l+1;
        end;
    end;
end;    

model.Q=sparse(opts.QP.qrow,opts.QP.qcol,opts.QP.qval);

A=[A; ones(1,n)];
b=[N*b; N];
	
subi=zeros(0,1);
subj=zeros(0,1);
valij=zeros(0,1);
l=1;
for i=1:(k+1)
    for j=1:n
        if(A(i,j)~=0)
            subi(l)=i;
            subj(l)=j;
            valij(l)=A(i,j);
            l=l+1;
        end;
    end;
end;

model.A=sparse(subi,subj,valij);
model.rhs=b;
   
for i=1:k
    model.sense(i)='<';
end;
model.sense(k+1)='=';

model.lb=zeros(n,1);
model.ub=N*ones(n,1);
vartypes=char(n,1);
	
for i=1:n
   	vartypes(i)='I';
end;    
model.vtype=vartypes;

clear params;
params.Presolve = 2;

% Run the Gurobi solver.
result=gurobi(model,params);
DQ=(result.x)./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