print("A demonstration of the cocktail algorithm for D-optimal designs (approximate theory)");
print("model: y ~ 1+ x1 + x1^2 + x2 + x1*x2");
print("design variable: -1<x1<1, 0<x2<1 discretized");

n=20;
m=5;
t=2*(1:n)/n-1;
s=(1:n)/n; 
X=matrix(nrow=n*n, ncol=m);
for(i in 1:n)
  for(j in 1:n)
    X[(i-1)*n+j, ]=c(1, t[i], t[i]^2, s[j], t[i]*s[j]); 

source('cocktail.r');
w=cocktail(X)[[1]];

print("output:");
print("optimal weights w(i)");
print(w[w>0]);