# R/mainfunctions.R In IsoCheck: Isomorphism Check for Multi-Stage Factorial Designs with Randomization Restrictions

```#================================================
# a function that applies a collineation C to a spread spr in order to relabel it.
applyCollineation <- function(C, spr){
ddd <- dim(spr)
out <- array(NA, ddd)
for(i in 1:ddd[3]){
for(j in 1:ddd[2]){
out[,j,i] <- (C %*% spr[,j,i])%%2
}
}
return(out)
}
#================================================

#================================================
# a function that converts a spread spr into its bitstring characterization (for checking equivalence)
getBitstrings <- function(spr){
ddd <- dim(spr)
n <- ddd[1]
mult <- 2^(0:(n-1))
strings <- matrix(0, nrow = ddd[3], ncol = 2^n -1)
smallests <- rep(NA, ddd[3])
for(flat in 1:(ddd[3])){
for(pencil in 1:(ddd[2])){
ind <- sum(mult * spr[,pencil, flat])
strings[flat, ind] <- 1
}
smallests[flat] <- which.max(strings[flat,] == 1)
}
sort_order <- sort(smallests, index.return = TRUE)\$ix
return(strings[sort_order,])
}
#================================================

#================================================
# Returns a list consisting of two components --- a boolean indicating whether
# or not spread1 is isomorphic to spread2, and a list of isomorphism establishing
# collineations. Note that when returnfirstIEC = TRUE, only the first IEC is
# returned (much faster runtime if they are isomorphic)
message("Spreads are not of same dimension.")
return(FALSE)
}
ell <- n/t
mu <- (2^n - 1)/(2^t - 1)
flatsize <- 2^t-1
spacesize <- 2^n -1

# note: for the code to work, the pencils in each flat in the spreads being checked for ismorphism
# should be arranged such that they are isomorphic to the yates order (basis elements in slots 1,2,4,8,etc.)
# furthermore, the union of the first ell flats should span PG(n-1,2)
# the following 2 steps permute within the equivalence class to guarantee these criteria are met.

pencilmappingoptions <- .getpenciloptions(t)

flatchoices <- gtools::combinations(mu, ell) # all of the ways to choose ell out of mu flats
perms <- gtools::permutations(ell,ell) # all of the ways to permute a sequence of ell flats
basischoices <- gtools::permutations(nrow(pencilmappingoptions),ell, repeats.allowed = TRUE) # all of the ways to choose bases from the ell flats

count <- 0 # the number of IECs found so far
validCBys <- list() # a list to collect our CBys corresponding to IECs
totaloptions <- nrow(flatchoices) * nrow(basischoices) * nrow(perms) # this should be equal to the size of the space to be searched

# iterate across flat choices for y1...yn
for(kk in 1:nrow(flatchoices)){
# first check if this choice of flats form a linearly independent set, if not skip.
if(is.na(xx[1,1]) == FALSE){
# iterate through the different options for bases of the flats
for(vv in 1:nrow(basischoices)){
# iterate through the different permutations
for(ww in 1:nrow(perms)){
# create the CBy
CBy <- matrix(nrow = n, ncol = n)
for(i in 1:ell){
for(j in 1:t){
CBy[,t * (i-1) +j] <- spread2[,pencilmappingoptions[basischoices[vv, i],j], flatchoices[kk, perms[ww,i]]]
}
}
bitsnew <- getBitstrings(relabelled) # get the bitstring characterization
if(sum(bitsnew != bits2) == 0){ # check equivalence. If equivalent, add to list
if(returnfirstIEC == TRUE){
return(list(result = TRUE, IECs = list((CBy %*% CxB) %%2)))
}
count <- count + 1
validCBys[[count]] <- CBy
}
}
message("percent done: ", 100 * round(((kk-1) *nrow(basischoices) * nrow(perms) + (vv)*nrow(perms))  /totaloptions ,4)) # to help keep track of how far along we are
}
}
}
if(count == 0){
if(printstatement == TRUE){
message("The two spreads are not isomorphic.")
}
return(list(result = FALSE, IECs = NA))
}

IECs <- validCBys
for(i in 1:length(IECs)){
IECs[[i]] <- (validCBys[[i]] %*% CxB) %%2
}
if(printstatement == TRUE){
message("The two spreads are isomorphic. For example, one isomorphism establishing collineation is")
print(IECs[[1]])
}

return(list(result = TRUE , IECs = IECs))
}
#================================================

#================================================
# a function that converts a star to its corresponding spread
ddd <- dim(star)
n <- ddd[1]
t <- log(ddd[2] + 1)/log(2)

# find nucleus
df <- suppressMessages(plyr::match_df(data.frame(t(star[,,1])), data.frame(t(star[,,2]))))

# find dimension of nucleus
t0 <- log(nrow(df) + 1)/log(2)

# find basis for nucelus
if(nrow(df) == 1){
nucleus <- as.matrix(df)
}else{
nucleus <- .rrefmod2(as.matrix(df))[1:t0,]
}

alls <- matrix(0, nrow = 2^n, ncol = n)
for(i in 1:n){
alls[(floor((0:(2^n-1))/(2^(i-1))) %%2 == 1),i] <- 1
}
alls <- alls[-1,]
newbasis <- nucleus
for(i in t0:(n - 1)){
newbasisspan <- t((t(newbasis) %*% t(alls[1:(2^i - 1), 1:i])) %% 2)
remalls <- suppressMessages(dplyr::anti_join(data.frame(alls), data.frame(newbasisspan)))
newbasis <- rbind(remalls[1,], newbasis)
}

# remove nucleus with collineation

collineation <- solve(t(newbasis))%%2
transformed <- applyCollineation(collineation, star)
spread <- array(NA, c(n - t0, 2^(t - t0) - 1 , ddd[3] ))
for(i in 1:ddd[3]){
spread[,,i] <- transformed[1:(n-t0),which(colSums(matrix(transformed[n-t0 + (1:t0), ,i], nrow = t0)) == 0),i]
}
}

# receives two stars star1 and star2 as arrays(n, 2^t-1, mu) and checks their isomorphism. Returns a list containing FALSE is the two stars
# are not isomorphic. If they are isomorphic, it returns a list consisting of T as well
# as a list of IECs from star1 to star2. If returnfirstIEC = TRUE, then only the
# If returnfirstIEC = TRUE the algorithm terminates as soon as the first IEC is found.
#================================================

#================================================
checkStarIsomorphism <- function(star1, star2, returnfirstIEC = FALSE){

if(sum(dim(star1) == dim(star2)) != 3){
message("Stars are not of same dimension.")
return(FALSE)
}
starcol1 <- reduce1\$collineation
starcol2 <- reduce2\$collineation

if(isocheck[[1]] == FALSE){
message("The two stars are not isomorphic.")
return(list(result = FALSE, IECs = NA))
}
else{
results <- list()
transafter <- starcol1
transbeforeinv <- starcol2
transbefore <- solve(transbeforeinv) %% 2
for(i in 1:length(isocheck[[1]])){
fullmat <- diag(dim(star1)[1])
results[[i]] <- (transbefore %*% fullmat %*% transafter) %%2
}
message("The two stars are isomorphic. For example, one isomorphism establishing collineation is")
print(results[[1]])
return(list(result = TRUE, IECs = results))
}
}
#================================================

#================================================
# checks the equivalence of star1 and star2.
checkStarEquivalence <- function(star1, star2){

if(sum(dim(star1) == dim(star2)) != 3){
message("Stars are not of same dimension.")
return(FALSE)
}
strings1 <- getBitstrings(star1)
strings2 <- getBitstrings(star2)
mat <- matrix(log(1:ncol(strings1)), nrow = nrow(strings1), ncol = ncol(strings1), byrow = TRUE)
if(sum(abs(sort(rowSums(strings1 * mat)) - sort(rowSums(strings2 * mat))) < 0.00001) == nrow(strings1)){
return(T)
}
else{
return(F)
}
}
#================================================

#================================================

message("Spreads are not of same dimension.")
return(FALSE)
}
if(sum(strings1 != strings2) == 0){
return(T)
}
return(F)
}
#================================================

#================================================
# printstatements = FALSE allows you to silence the explanation for the failure
is.spread <- function(spr, printstatements = TRUE){
dims <- dim(spr)
if(length(dims) != 3){
if(printstatements == TRUE){
message("Spread should be given as a 3d binary array")
}
return(FALSE)
}

# weird corner case to address
if(sum(dims == c(1,1,1)) == 3){
if(sum(spr != array(1,c(1,1,1))) == 0){
return(TRUE)
}
}

# check if the representation is binary
if(sum(sort(unique(c(spr))) != c(0,1)) > 0){
if(printstatements == TRUE){
message("Spread should be expressed in binary")
}
return(FALSE)
}

t <- log(dims[2] + 1)/log(2)

# check if the flats are of a valid size
test1 <- floor(t) == t
if(test1 == FALSE){
if(printstatements == TRUE){
message("Invalid size of flats. Should be 2^t-1 for some t.")
}
return(FALSE)
}

# make sure the right number of flats are present
test2 <- ((2^dims[1] - 1)/(dims[2]) == dims[3])
if(test2 == FALSE){
if(printstatements == TRUE){
message("The number of flats is not equal to (2^n-1)/(2^t-1).")
}
return(FALSE)
}

# ensure that no null effects are included
for(i in 1:dims[2]){
for(j in 1:dims[3]){
if(sum(spr[,i,j]) == 0){
if(printstatements == TRUE){
message("Null effects should not be present.")
}
return(FALSE)
}
}
}

# check for repeated entries
bitstrings <- getBitstrings(spr)
if(max(bitstrings) > 1){
if(printstatements == TRUE){
message("Duplicate entries in a flat.")
}
return(FALSE)
}

# make sure the flats partition the space
if(sum(colSums(bitstrings) != rep(1, 2^(dims[1]) - 1)) != 0){
if(printstatements == TRUE){
message("The flats do not form a partition of the space.")
}
return(FALSE)
}

# make sure each flat is actually a flat (is the span of t elements)
for(i in 1:dims[3]){
if(max(which(rowSums(.rrefmod2(t(spr[,,i]))) > 0)) != t){
if(printstatements == TRUE){
message("Flat ", i, " is not a valid flat.")
}
return(FALSE)
}
}
return(TRUE)
}
#================================================

#================================================
# printstatements = FALSE allows you to silence the explanation for the failure
is.star <- function(star, printstatements = TRUE){

dims <- dim(star)
if(length(dims) != 3){
if(printstatements == TRUE){
message("Star should be given as a 3d binary array")
}
return(FALSE)
}

# weird corner case to address
if(sum(dims == c(1,1,1)) == 3){
if(sum(star != array(1,c(1,1,1))) == 0){
return(TRUE)
}
}

# check if the representation is binary
if(sum(sort(unique(c(star))) != c(0,1)) > 0){
if(printstatements == TRUE){
message("Star should be expressed in binary")
}
return(FALSE)
}

t <- log(dims[2] + 1)/log(2)

# check if the flats are of a valid size
test1 <- floor(t) == t
if(test1 == FALSE){
if(printstatements == TRUE){
message("Invalid size of flats. Should be 2^t-1 for some t.")
}
return(FALSE)
}

# ensure that no null effects are included
for(i in 1:dims[2]){
for(j in 1:dims[3]){
if(sum(star[,i,j]) == 0){
if(printstatements == TRUE){
message("Null effects should not be present.")
}
return(FALSE)
}
}
}

# check for repeated entries
bitstrings <- getBitstrings(star)
if(max(bitstrings) > 1){
if(printstatements == TRUE){
message("Duplicate entries in a flat.")
}
return(FALSE)
}

inclusioncounts <- colSums(bitstrings)

# check if all effects are included
if(sum(0 == inclusioncounts) > 0){
if(printstatements == TRUE){
message("The elements do not form a cover.")
}
return(FALSE)
}

# check if the only overlap among flats is the nucleus.
if(length(unique(inclusioncounts)) > 2){
if(printstatements == TRUE){
message("Nucleus is not common to all flats or a.")
}
return(FALSE)
}

nucleus_size <- sum(inclusioncounts > 1)

t0 <- log(nucleus_size + 1)/log(2)

if(t0 == 0){
message("It is a spread (i.e., trivial star).")
}

# check if the nucleus is of a valid size
test1 <- floor(t0) == t0
if(test1 == FALSE){
if(printstatements == TRUE){
message("Invalid size of nucleus. Should be 2^t-1 for some t.")
}
return(FALSE)
}

# make sure each flat is actually a flat (is the span of t elements)
for(i in 1:dims[3]){
if(max(which(rowSums(.rrefmod2(t(star[,,i]))) > 0)) != t){
if(printstatements == TRUE){
message("Flat ", i, " is not a valid flat.")
}
return(FALSE)
}
}
return(TRUE)
}

#=================================================
vectortostring <- function(arry){
dim_ord = dim(arry)
if(sum(dim_ord) == 0){
arry = as.matrix(arry)
dim_ord = dim(arry)
}
n = dim_ord[1]
k = dim_ord[2]
if(length(dim_ord)==2){
dim(arry) <- c(dim_ord,1)
dim_ord = c(dim_ord,1)
}
mu = dim_ord[3]
str_arry <- array(0,c(dim_ord[-1]))

for(i_mu in 1:mu){
for(i_k in 1:k){

vec = arry[,i_k,i_mu]
len_vec = length(vec)
lett_val = LETTERS[vec*seq(1:len_vec)]
final_lett_val=NULL;
for(i_n in 1:length(lett_val)){
final_lett_val = paste(final_lett_val, lett_val[i_n],sep="")
}
str_arry[i_k,i_mu]=final_lett_val
}
}
return(str_arry)
}
#=================================================

stringtovector <- function(string,n){
ret <- rep(0, n)
for(i in 1:n){
ret[i]<- (grepl(LETTERS[i], string)==TRUE)
}
return(ret)
}

#================================================
```

## Try the IsoCheck package in your browser

Any scripts or data that you put into this service are public.

IsoCheck documentation built on March 25, 2020, 5:09 p.m.