Nothing
## Densest packing-based maximum projection designs
## Author: Xu He, hexu@amss.ac.cn
ProjSepD <- function(design) {
.C(ProjSep, PACKAGE="LatticeDesign",
as.double(design), as.integer(dim(design)[2]), as.integer(dim(design)[1]), as.double(rep(0,dim(design)[2])) )[[4]]
}
DPMPD <- function(p=2, n, rotation="magic", w=100){
if(!p>=2|!p<=8|!p==round(p)) stop("p must be an integer greater than one and no greater than eight.")
if(!n>=2|!n==round(n)) stop("n must be an integer greater than one.")
if(!w>=1|!w==round(w)|!w<=10000) stop("w must be a positive integer greater no greater than 10000.")
#### Generate the unrotated candidate sets.
detG <- 2
if(p==2) detG <- 2*sqrt(3)
if(p==6) detG <- sqrt(3)
if(p==7) detG <- sqrt(2)
if(p==8) detG <- 1
r <- (n*detG)^(1/p)/2*sqrt(p)+1
ss <- floor( r )
FE <- matrix(0,ss+1,2)
FE[,2] = 0:ss
FE[,1] = FE[,2]^2
for(j in 2:p) { # Generate all nonnegative integer points with radius <= r
nFE <- matrix(0,0,j+1)
for(k in 0:ss) { oFE = cbind(FE,rep(k,dim(FE)[1])); oFE[,1] = oFE[,1]+k^2; nFE <- rbind(nFE,oFE); }
FE <- nFE[nFE[,1]<= floor( r^2 ), ]
}
FE[,1]=0; for(j in 2:(p+1)) FE[,1]=FE[,1]+FE[,j]; FE = FE[ floor(FE[,1]/2)*2==FE[,1] ,]; ## Delete half
if(p>=6) { ## Add half while reflecting to negative first dimension
FE1 = FE[ FE[,2]>0, ]
FE1[,2] = -FE1[,2]
FE2 = FE[ FE[,2]==1, ]
FE2[,2] = -FE2[,2]
FE3 = rbind(FE,FE2)
FE4 = FE3 +0.5
FE5 = FE +0.5
FE5[,2] = -FE5[,2]-1
FE <- rbind(FE,FE1,FE4,FE5)
}
if(p<6) { ## Reflect to negative first dimension
nFE = FE[FE[,2]>0,]
nFE[,2] <- -nFE[,2]
FE <- rbind(FE,nFE)
}
if(p==2) { FE[,3] <- FE[,3]*sqrt(3); } ## Rescale the last dimension
if(p==6) { FE[,7] <- FE[,7]*sqrt(3); }
if(p==7) { FE[,8] <- FE[,8]*sqrt(2); }
FE[,1] <- 0;
for(j in 2:(p+1)) FE[,1] <- FE[,1] +FE[,j]^2
FE <- FE[FE[,1]<=r^2,]
FE <- FE[,-1]
for(j in 2:p) { ## Reflect to negative second-to-pth dimensions
nFE = FE[FE[,j]>0,]
nFE[,j] <- -nFE[,j]
nFE[,1] <- -nFE[,1]
FE <- rbind(FE,nFE)
}
## Prepare rotation for p=4 and 8
mv = 25
if(w>100) mv = 55
vlist <- matrix(0,0,3)
for(v1 in 1:mv) for(v2 in v1:mv) {
vlist = rbind(vlist,c(v1,v2,v1^2+v2^2))
}
vvlist <- matrix(0,0,6)
for(i in 1:dim(vlist)[1]) for(j in 1:dim(vlist)[1]) {
ratio = vlist[j,3]/vlist[i,3]
if(ratio>1) if(floor(sqrt(ratio))!=sqrt(ratio)) if(floor(ratio)*vlist[i,3]==vlist[j,3]) vvlist = rbind(vvlist,c(vlist[i,1:2],vlist[j,1:2],ratio,0))
}
vvlist[,6] = abs( (vvlist[,3]-vvlist[,1]*sqrt(vvlist[,5])) / sqrt( (vvlist[,3]-vvlist[,1]*sqrt(vvlist[,5]))^2 + (vvlist[,4]-vvlist[,2]*sqrt(vvlist[,5]))^2 ) )
j=2
while(j<=dim(vvlist)[1]) {
for(k in 1:j) {
if(k<j) if(abs(vvlist[j,6]-vvlist[k,6])<10^-10) { vvlist = vvlist[-j,]; break; }
}
if(k==j) j=j+1
}
for(i in 1:dim(vvlist)[1]) vvlist[i,6] = max(vvlist[i,1:4])
vvlist1 = vvlist[vvlist[,5]==2,]
vvlist2 = vvlist[vvlist[,5]==5,]
vvlist3 = vvlist[vvlist[,5]==13,]
if(w<=100) {
vvlist1 = vvlist1[ order(vvlist1[,6])[1:10], ]
vvlist2 = vvlist2[ order(vvlist2[,6])[1:10], ]
vvlist3 = vvlist3[ order(vvlist3[,6])[1:10], ]
}
# set.seed(1)
#### Try w rotation matrices
ress = matrix(0,w,p+1) # Minimum of separation distance of 1,...,p dimensional projections, overall score.
l <- (n*detG)^(1/p)
maxscore <- -10^10
for(ll in 1:w) {
# set.seed(ll)
#### Generate rotation matrix
if(rotation!="magic"|p==5|p==7) { ## Random rotation
CPair <- matrix(0,p*(p-1)/2,2)
row <- 1
for(i in 1:(p-1)) for(j in (i+1):p) { CPair[row,1] <- i; CPair[row,2] <- j; row <- row+1; }
R <- diag(p)
for(a in 1:((p*(p-1)/2))) {
alpha <- runif(1,0,2*pi)
thepair <- CPair[a,]
W <- diag(p)
W[thepair[1],thepair[1]] <- cos(alpha); W[thepair[2],thepair[2]] <- cos(alpha);
W[thepair[1],thepair[2]] <- sin(alpha); W[thepair[2],thepair[1]] <- -sin(alpha);
R <- R%*%W
}
}
if(rotation=="magic"&p==2) { ## Magic rotation matrices
q1 <- 3
bb1 <- 20; bb2 <- 4;
if(w>100) { bb1 <- floor(w^(2/3)); bb2 <- floor(w^(1/3)); }
u1 = ceiling(ll/bb1)
u2 = ceiling( ( ll - (u1-1)*bb1 ) / bb2 )
u3 = ll - (u1-1)*bb1 - (u2-1)*bb2 -floor(bb2/2)
if(u3^2==3*u1^2+u2^2) u3 = bb2 -floor(bb2/2)+1
R = rbind( c(u1*sqrt(q1)+u3,-u2), c(u2,u1*sqrt(q1)+u3) )
for(j in 1:2) R[,j] = R[,j]/sqrt(sum(R[,j]^2))
R1 = rbind( c(sqrt(3)+1,sqrt(3)-1), c(-sqrt(3)+1,sqrt(3)+1) ) /2/sqrt(2)
R = R1 %*% R
}
if(rotation=="magic"&p==3) { ## Magic rotation matrices
qlist = 2:(w+10)
for(k in 2:20) qlist = qlist[ floor(qlist/k^3)*k^3!=qlist ]
bb1 <- 10; bb2 <- 55;
if(w>100) bb2 <- floor((w+10)/2)
q1 = qlist[ll];
if(ll>bb1) q1 = qlist[ll-bb1]/8
if(ll>bb2) q1 = qlist[ll-bb2]/27
R3=matrix(0,3,3)
R3[1,1]=R3[3,2]=R3[2,3] = 1-q1
R3[2,1]=R3[1,2]=R3[3,3] = q1-q1^(1/3)
R3[3,1]=R3[2,2]=R3[1,3] = (1-q1)*(q1^(1/3)+q1^(2/3))
R <- R3 / sqrt(sum(R3[,3]*R3[,3]))
}
if(rotation=="magic"&p==4) { ## Magic rotation matrices
bb1 <- 10;
if(w>100) bb1 <- floor(sqrt(w))
if(w>bb1*76) bb1 <- ceiling(w/76)
rot1 = vvlist1[ceiling(ll/bb1),]
GRR1 = rbind( rot1[c(3,1)], rot1[c(4,2)] )
q1 = rot1[5]
R1 <- GRR1 %*% rbind( c(1,1), c(-sqrt(q1),sqrt(q1)) )
for(j in 1:2) R1[,j] = R1[,j]/sqrt(sum(R1[,j]^2))
rot2 = vvlist2[ll-ceiling(ll/bb1)*bb1+bb1,]
GRR2 = rbind( rot2[c(3,1)], rot2[c(4,2)] )
q2 = rot2[5]
R2 <- GRR2 %*% rbind( c(1,1), c(-sqrt(q2),sqrt(q2)) )
for(j in 1:2) R2[,j] = R2[,j]/sqrt(sum(R2[,j]^2))
R <- rbind( cbind(R2[1,1]*R1,R2[1,2]*R1), cbind(R2[2,1]*R1,R2[2,2]*R1) )
}
if(rotation=="magic"&p==6) { ## Magic rotation matrices
bb1 <- 50; bb2 <- 25;
if(w>100) { bb2 <- floor((w/5)^(1/3)*5); bb1 <- floor((w/5)^(2/3)*5); }
q1 <- 3
u1 = ceiling(ll/bb1)
u2 = ceiling( ( ll - (ceiling(ll/bb1)-1)*bb1 ) / bb2 )
u3 = ceiling( ( ll - (ceiling(ll/bb2)-1)*bb2 ) / 5 ) - ceiling(bb2/10)
if(u3==0) u3 = ceiling(bb2/5) - ceiling(bb2/10) + 1
if(w<=100) if(u1==2&u2==2&abs(u3)==2) u1 = 3
R2 = rbind( c(u1*sqrt(q1)+u3,-u2), c(u2,u1*sqrt(q1)+u3) )
for(j in 1:2) R2[,j] = R2[,j]/sqrt(sum(R2[,j]^2))
qlist = 2:100
for(k in 2:4) qlist = qlist[ floor(qlist/k^3)*k^3!=qlist ]
q1 = qlist[ ll-(ceiling(ll/5)-1)*5 ];
R3=matrix(0,3,3)
R3[1,1]=R3[3,2]=R3[2,3] = 1-q1
R3[2,1]=R3[1,2]=R3[3,3] = q1-q1^(1/3)
R3[3,1]=R3[2,2]=R3[1,3] = (1-q1)*(q1^(1/3)+q1^(2/3))
R3 <- R3 / sqrt(sum(R3[,3]*R3[,3]))
R <- rbind( cbind(R3[1,1]*R2,R3[1,2]*R2,R3[1,3]*R2), cbind(R3[2,1]*R2,R3[2,2]*R2,R3[2,3]*R2), cbind(R3[3,1]*R2,R3[3,2]*R2,R3[3,3]*R2) )
}
if(rotation=="magic"&p==8) { ## Magic rotation matrices
bb2 <- floor(w^(1/3)); bb1 <- floor(w^(2/3));
rot1 = vvlist1[ ceiling(ll/bb1) ,]
GRR1 = rbind( rot1[c(3,1)], rot1[c(4,2)] )
q1 = rot1[5]
R1 <- GRR1 %*% rbind( c(1,1), c(-sqrt(q1),sqrt(q1)) )
for(j in 1:2) R1[,j] = R1[,j]/sqrt(sum(R1[,j]^2))
rot2 = vvlist2[ ceiling( (ll-ceiling(ll/bb1)*bb1+bb1) / bb2 ) ,]
GRR2 = rbind( rot2[c(3,1)], rot2[c(4,2)] )
q2 = rot2[5]
R2 <- GRR2 %*% rbind( c(1,1), c(-sqrt(q2),sqrt(q2)) )
for(j in 1:2) R2[,j] = R2[,j]/sqrt(sum(R2[,j]^2))
rot3 = vvlist3[ ll-ceiling(ll/bb2)*bb2+bb2 ,]
GRR3 = rbind( rot3[c(3,1)], rot3[c(4,2)] )
q3 = rot3[5]
R3 <- GRR3 %*% rbind( c(1,1), c(-sqrt(q3),sqrt(q3)) )
for(j in 1:2) R3[,j] = R3[,j]/sqrt(sum(R3[,j]^2))
R <- rbind( cbind(R2[1,1]*R1,R2[1,2]*R1), cbind(R2[2,1]*R1,R2[2,2]*R1) )
R <- rbind( cbind(R3[1,1]*R,R3[1,2]*R), cbind(R3[2,1]*R,R3[2,2]*R) )
}
#### Find delta and thus epsilon and design
E <- FE %*% R
isE <- FALSE
isB <- FALSE
isL <- FALSE
epsilonV <- matrix(0,3,p)
for(i in 1:100) {
while(1) {
delta <- runif(p,-1,1)
if( sum( sort(abs(delta))[(p-1):p] ) > 1) next
if(p==2) delta[p] <- delta[p]*sqrt(3)
if(p==6) delta[p] <- delta[p]*sqrt(3)
if(p==7) delta[p] <- delta[p]*sqrt(2)
if( sum( (abs(delta)+l/2)^2 ) <= r^2) break
}
epsilon <- delta %*% R
design <- E
for(j in 1:p) { if(length(design)<=p) break; design <- design[ abs(design[,j]-epsilon[j])<l/2 ,]; }
if(length(design)<=p) next;
if(dim(design)[1]==n) { isE <- TRUE; epsilonV[1,] <- epsilon; break; }
if(isL==FALSE & dim(design)[1]<n) { isL <- TRUE; epsilonV[2,] <- epsilon; }
if(isB==FALSE & dim(design)[1]>n) { isB <- TRUE; epsilonV[3,] <- epsilon; }
if(isL==TRUE & isB==TRUE) break;
}
if(isE==FALSE) for(j in 1:(p-1)) {
epsilon <- c(epsilonV[2,1:j],epsilonV[3,(j+1):p])
design <- E
for(j in 1:p) { if(length(design)<=p) break; design <- design[ abs(design[,j]-epsilon[j])<l/2 ,]; }
if(length(design)<=p) next;
if(dim(design)[1]==n) { isE <- TRUE; epsilonV[1,] <- epsilon; break; }
if(dim(design)[1]<n) { epsilonV[2,] <- epsilon; break; }
if(dim(design)[1]>n) { epsilonV[3,] <- epsilon; }
}
if(isE==FALSE) for(k in 1:100) {
epsilon <- (epsilonV[2,]+epsilonV[3,])/2
design <- E
for(j in 1:p) { if(length(design)<=p) break; design <- design[ abs(design[,j]-epsilon[j])<l/2 ,]; }
if(length(design)<=p) next;
if(dim(design)[1]==n) { isE <- TRUE; epsilonV[1,] <- epsilon; break; }
if(dim(design)[1]<n) { epsilonV[2,] <- epsilon; }
if(dim(design)[1]>n) { epsilonV[3,] <- epsilon; }
}
for(j in 1:p) design[,j] <- (design[,j]-epsilon[j]) / l + .5
#### Evaluate design
ress[ll,1:p] = ProjSepD(design)
ress[ll,p+1] = sum( log(ress[ll,1:p]*n^(1/(1:p))) * c(1/2,2:p) )
if(ll==1|ress[ll,p+1]>maxscore) {
maxscore <- ress[ll,p+1]
designBest <- design
}
}
for(j in 1:p) designBest[,j] = designBest[,j]-(min(designBest[,j])+max(designBest[,j])-1)/2
return(list(Design=designBest,ProjectiveSeparationDistance=sqrt(ProjSepD(designBest))))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.