** Find the Phi_p optimal design for polynomial regression model with order (d-1): eta(x,theta)=theta_1+theta_2*x+...+theta_d*x**(d-1);
** lambda(x)=1-x^2 is the efficiency function;
** 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;
** the order of the polynomial is d-1;
d=6;
** WB is the matrix of interested parameters;
WB=i(d);  
** space is the design space;
space={-1 1};
** there is no fixed design points because the efficiency function lambda(x)=1-x^2;

** find z=(x1,...,xd,w2,...,wd) to minimize the variance-covariance matrix;
** z0 is the starting point for the search;
z0=((2*(1:d)`-1)/d-1)//j(d-1,1,1)/d;  
print z0;

** lb, ub are the lower bound, upper bound for z, respectively;
lb=j(d,1,space[1])//j(d-1,1,0);
ub=j(d,1,space[2])//j(d-1,1,1);

** z is under the constraint of constr_A*z<constr_b, here it translates into w2+...+wd<1;
constr_A=j(1,d,0)||j(1,d-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) global (d);
    ft=j(d,1,1);
	ft[1,1]=1;
	do i=2 to d;
		ft[i,1]=x**(i-1);
	end;
    ** recall that (1-x^2) is the efficiency function;
    infor=ft*ft`*(1-x**2);
    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,wb) global(d);

    wd=1-sum(z[d+1:2*d-1]);
	temp=z//wd;
	infor=j(d,d,0);
	do i=1 to d;
    	infor=infor+temp[d+i]*infor1(temp[i]);
	end;

	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,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,wb)-phi_value(p,zminus,wb))/2/eps;
		
	end;
	return(grad);
finish;

** the main part, find the critical point using quasi-Newton algorithm;
start QNewton_BFGS(p,z0,wb,lb, ub,constr_A,constr_b,max_iter);
	nvar=nrow(z0);
	grad0=gradient(p,z0,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,wb)-phi_value(p,z,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,wb)-phi_value(p,z,wb);
		end;

        ** update;
        old_grad=grad;
		z=newz ;
  
     	grad=gradient(p,z,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,wb,lb, ub,constr_A,constr_b,max_iter);
value=phi_value(p,crit_pt,wb);
print crit_pt value;

quit;