** Find the Phi_p optimal design for LINEXP model: eta(x,theta)=theta1+theta2*exp(theta3*x)+theta4*x;
** A quasi-Newton algorithm using BFGS update is implemented to find the critical point, which gives the optimal design;
** You need to change the fisher information function and settings about the parameters, design boundary, 
   etc, to accommodate other models in the paper;

proc iml;
** p=0 meand D-optimality, p=1 means A-optimality;
p=1;
** WB is the matrix of interested parameters;
WB=i(4);  
** theta is the value of parameter;
theta={1 0.5 -1 1}; 
** space is the design space;
space={0 1};
** x1 and x4 are fixed design points, by complete class result;
x1=space[1]; x4=space[2];

** find z=(x2,x3,w2,w3,w4) to minimize the variance-covariance matrix, z0 is the starting point for the search;
z0={0.1,0.2,0.3,0.3,0.2};  
** lb, ub are the lower bound, upper bound for z, respectively;
lb={0,0,0,0,0};
ub={1,1,1,1,1};
** z is under the constraint of constr_A*z<constr_b, here it translates into w2+w3+w4<1;
constr_A={0 0 1 1 1};
constr_b=1;

max_iter=100;

** information matrix for a single design point;
** you need to change this function for a different model;
start infor1(x,theta);
    ft=j(4,1,1);
    ft[1,1]=1;
    ft[2,1]=exp(theta[3]*x);
    ft[3,1]=theta[2]*x*exp(theta[3]*x);
    ft[4,1]=x;
    infor=ft*ft`;
    return(infor);
finish;


** the objective function based on Phi_p optimality, to be minimized;
** you need to change this for a different model;
start phi_value(p,z,theta,wb) global (x1,x4);
    x2=z[1];
	x3=z[2];
	w2=z[3];
	w3=z[4];
	w4=z[5];
    w1=1-w2-w3-w4;
    infor=w1*infor1(x1,theta)+w2*infor1(x2,theta)+w3*infor1(x3,theta)+w4*infor1(x4,theta);
	
	if p=0 then value=log(det(wb*ginv(infor)*wb`));
	else value=log(trace((wb*ginv(infor)*wb`)**p)/nrow(wb))*(1/p);
    return(value);
finish;

** the function that computes gradient nummerically;
start gradient(p,z,theta,wb);
	nvar=nrow(z);
	grad=j(nvar,1,0);
	eps=0.000001;
	do i=1 to nvar;
		zplus=z;
		zplus[i]=z[i]+eps;
		zminus=z;
		zminus[i]=z[i]-eps;
		grad[i]=(phi_value(p,zplus,theta,wb)-phi_value(p,zminus,theta,wb))/2/eps;
		
	end;
	return(grad);
finish;

** the main part, find the critical point using quasi-Newton algorithm;
start QNewton_BFGS(p,z0,theta,wb,lb, ub,constr_A,constr_b,max_iter);
	nvar=nrow(z0);
	grad0=gradient(p,z0,theta,wb);
	 
	grad=grad0;
	z=z0;
	invhess=i(nvar);
    iter=1;
    ** the converging criteria;
	do while ( (iter<max_iter)&(grad`*grad>0.000000001*grad0`*grad0) );
	    
    	delta=-invhess*grad;
		
		newz=z+delta;
        ** adjust the step so that newz is within the bounds and constraint is satisfied; 
		do while ( (any(newz-lb<0) | any(newz-ub>0) | constr_A*newz>1 ) );
			delta=delta/2;
			newz=z+delta;
			
		end;

		delta=delta/2;
	    newz=z+delta;
		diff=phi_value(p,newz,theta,wb)-phi_value(p,z,theta,wb); 
 
        ** further adjust the step so that we don't overstep;
		do while ( diff>0.001*grad`*delta );
			delta=delta/2;
			newz=z+delta;
			diff=phi_value(p,newz,theta,wb)-phi_value(p,z,theta,wb);
		end;

        ** update;
        old_grad=grad;
		z=newz ;
  
     	grad=gradient(p,z,theta,wb);
  		gamma=grad-old_grad;

        ** BFGS estimate of the inverse hessian;
  		if (delta`*gamma>0) then do;
    		c1=delta`*gamma;
    		c2=(1+gamma`*invhess*gamma/c1)/c1;
    		dd=delta*delta`;
    		dg=delta*gamma`;
    		gd=gamma*delta`;  
    		M=c2*dd-(dg*invhess+invhess*gd)/c1;
    		invhess=invhess+M;
		end;
  
  		else invhess=i(nvar);
		iter=iter+1;
		
	end;
	print iter;
	return (z);
finish;

** Now everything is defined, and we can call QNewton_BFGS function to find the critical point;
crit_pt=QNewton_BFGS(p,z0,theta,wb,lb, ub,constr_A,constr_b,max_iter);
value=phi_value(p,crit_pt,theta,wb);
print crit_pt value;

quit;