R/MOLS.R

In blocksdesign: Nested and Crossed Block Designs for Factorial and Unstructured Treatment Sets

#' @title Prime power MOLS from finite fields
#' 
#' @description
#' Constructs r sets of mutually orthogonal Latin squares (MOLS) of dimension p**q for prime p and
#' integer power q where r < p**q. Memory issues mean that the maximum size of the exponent q for 
#' specific p is restricted to the values shown in the table below:
#' 
#'    \tabular{rrrrrrrr}{
#'    \bold{prime p} \tab                                       \bold{maximum q}\cr
#'    2 \tab                                                                13\cr
#'    3 \tab                                                                 8\cr
#'    5 \tab                                                                 6\cr                                                      
#'    7 \tab                                                                 5\cr
#'    11 \tab                                                                4\cr
#'    13 17 19 23 \tab                                                       3\cr
#'    29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 \tab                   2\cr
#'    Any prime >97 \tab                                                     1
#'    }
#' 
#'    
#' @details Generates MOLS by cyclic permutation of a basic Latin square constructed from a vector
#' of ordered elements of a prime-power finite field of size p**q (see Chapter 1 of Raghavarao 1971).
#' 
#' The primitive polynomials for the MOLS generated by this package were
#' extracted from the Table of Primitive Polynomials given in the Supplement to Hansen and Mullen (1992).
#' 
#' The output is a single data frame of size \eqn{p**q x (r+2)} for the required set of MOLS
#' with a column for the rows classification, a column for the columns classification and 
#' separate columns for each treatment set of the required set of squares.
#' 
#' Also see the function \code{GraecoLatin} which will generate pairs of MOLS for a range of 
#' non-prime power design sizes v**2 including
#' all odd values of v and any even valued v <= 30 except for 6 or 2.
#' 
#' @seealso \code{\link{GraecoLatin}}
#' 
#' @references
#' Hansen, T. & Mullen, G. L. (1992) Primitive polynomials over finite fields,
#' Mathematics of Computation, 59, 639-643 and Supplement.
#' 
#' Raghavarao D. (1971) Constructions and Combinatorial Problems in Design of Experiments,
#' Dover Publications, Inc. Section 1.3
#'  
#' @param p is any suitable integer base
#' 
#' @param q is any suitable integer exponent
#' 
#' @param r is any suitable number of squares
#' 
#' @return
#' Data frame of factor levels for rows, columns and treatment sets
#'  
#' @examples
#' MOLS(2,3,7) # Seven MOLS of size 8 x 8
#' MOLS(3,2,4) # Four MOLS of size 9 x 9
#' \donttest{MOLS(3,3,4)} # Four MOLS of size 27 x 27
#' \donttest{MOLS(23,2,2)} # Two MOLS of size 529 x 529
#'  
#' @export
  MOLS=function(p,q,r=1) {
    v=p**q
  # Latin squares
  if (r==1 )
    return(data.frame(Row=rep(1:v,each=v),Col=rep(1:v,v),T1=unlist(lapply(0:(v-1),function(j){(rep(0:(v-1))+j)%%v}))+1))
   
    #MOLS
  Z=isPrimePower(p)
  
  ## tests if p is a prime-power such that N = p**(i*q) for i>1
  if (!isPrime(p) & !isFALSE(Z)) # (!isFALSE(Z)) is TRUE if p = is a prime power  
    return(paste("The parameters (",p,",",q,",",r,") are not valid. Did you mean (",Z$base,",",q*Z$power,",",r,") ?" ))
  
  ## tests if p is a prime given that p is NOT a prime-power 
  if (isFALSE(Z)) ## not a prime and not a prime-power
    return(paste("p must be a prime if r>1" ))

  if ( r>(p**q-1) ) return(paste("The maximum number of orthogonal squares for MOLS of size", p**q,"is", p**q-1))
  if (p > 97 & q>1) return(("Not a valid size - if p is a prime  > 97 then q must equal 1"))  

  # cyclic MOLS based on primes
  if (q==1 ) {
  mols=sapply(0:r,function(z){ sapply(0:(p-1), function(j){ (rep(0:(p-1))*z +j)%%p}) })
  mols=data.frame(rep(1:p,p),mols+1)
  mols[c(1,2)]=mols[c(2,1)]
  colnames(mols)=c("Row","Col",paste("T", 1:r, sep = ""))
  return(mols)
  }
    
  ##  available prime powers
  primes=c(2, 3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97) 
  powers=c(13,8,6,5, 4, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2)
  
  ## primitive polynomials for prime powers less than 100 from the Table of Primitive Polynomials given in the 
  ## Supplement to Hansen and Mullen (1992).

  pp2=list(
  c(1,1,1),
  c(1,0,1,1),
  c(1,0,0,1,1),
  c(1,0,0,1,0,1),
  c(1,0,0,0,0,1,1),
  c(1,0,0,0,0,0,1,1),
  c(1,0,0,0,1,1,1,0,1),
  c(1,0,0,0,0,1,0,0,0,1),
  c(1,0,0,0,0,0,0,1,0,0,1),
  c(1,0,0,0,0,0,0,0,0,1,0,1),
  c(1,0,0,0,0,0,1,0,1,0,0,1,1),
  c(1,0,0,0,0,0,0,0,0,1,1,0,1,1)) 
  pp3=list(
  c(1,1,2),
  c(1,0,2,1),
  c(1,0,0,1,2),
  c(1,0,0,0,2,1),
  c(1,0,0,0,0,1,2),
  c(1,0,0,0,0,2,0,1),
  c(1,0,0,0,0,1,0,0,2))
  pp5=list(
  c(1,1,2),
  c(1,0,3,2),
  c(1,0,1,2,2),
  c(1,0,0,0,4,2),
  c(1,0,0,0,0,1,2))
  pp7=list(
  c(1,1,3),
  c(1,0,3,2),
  c(1,0,1,3,5),
  c(1,0,0,0,1,4))
  pp11=list(
  c(1,1,7),
  c(1,0,1,4),
  c(1,0,0,1,2))
  pp13=list(
  c(1,1,2),
  c(1,0,1,6))
  pp17=list(
  c(1,1,3),
  c(1,0,1,3))
  pp19=list(
  c(1,1,2),
  c(1,0,1,4))
  pp23=list(
    c(1,1,7),
    c(1,0,1,3))
  pp29=list(c(1,1,3))
  pp31=list(c(1,1,12))
  pp37=list(c(1,1,5))
  pp41=list(c(1,1,12))
  pp43=list(c(1,1,3))
  pp47=list(c(1,1,13))
  pp53=list(c(1,1,5))
  pp59=list(c(1,1,2))
  pp61=list(c(1,1,2))
  pp67=list(c(1,1,12))
  pp71=list(c(1,1,11))
  pp73=list(c(1,1,11))
  pp79=list(c(1,1,3))
  pp83=list(c(1,1,2))
  pp89=list(c(1,1,6))
  pp97=list(c(1,1,5))

primpol=list(pp2,pp3,pp5,pp7,pp11,pp13,pp17,pp19,pp23,pp29,pp31,pp37,pp41,pp43,pp47,pp53,pp59,pp61,pp67,pp71,pp73,pp79,pp83,pp89,pp97)
for (i in 1:25)
  primpol[[i]]=lapply(primpol[[i]],rev)

index=which(p==primes)
if ( q>powers[index] ) 
  return(paste("Not a valid power - q must be not more than ", powers[index ]))

checkff = function(p,q,X) {
  rowMax=q
  C=lapply(X,coef)
  M = do.call(rbind, lapply(C, function(x){ 
    length(x) <- rowMax 
    x 
  })) 
  if (any(duplicated(data.frame(M))==TRUE)) stop(paste("Could not find a set of primitive elenents for p = ",p," and q = ",q))
}

## functions ued to calculate Galois Fields (gf)
reduce = function(X,primpol,p,q) {
  c=coef(X)
  if(length(c)>q) 
    c=coef(PolynomF::polynomial(c[1:q])+primpol*c[q+1])
  X=PolynomF::polynomial(c%%p)
  X
}

gf = function(p,q,index) {
  primpol=primpol[[index]][[q-1]]
  primpol=PolynomF::polynomial(  (-primpol[1:q])%%p  )
  X = vector(mode = "list", length = p**q)
  X[[1]]=PolynomF::polynomial(c(0,0))
  X[[2]]=PolynomF::polynomial(c(1,0)) 
  X[[3]]=PolynomF::polynomial(c(0,1))
  for( i in 4:v) { 
    X[[i]]=X[[i-1]]*X[[3]]
    X[[i]]=reduce(X[[i]],primpol,p,q)
  }
  checkff(p,q,X)
  X
}

modcoef=function(z) {
  coef=lapply(z,  function(x) coef(x)%%p)
  rowMax= max(sapply(coef, length))
  M= do.call(rbind, lapply(coef,function(x){ 
    length(x) = rowMax
    x 
  }))  
  M[is.na(M)] = 0
  M
}

## mutually orthogonal Latin squares by method of Raghavarao
molsq=function(p,q,r,index) {
  shift=c(1:v , (1+2*v):(v*v) , (1+v):(2*v) )
  gf=gf(p,q,index)
  M=modcoef(gf)
  LS=rep(0,v*v)
  for (i in 1:ncol(M))
    LS= LS +   (rep(M[,i],v) + rep(M[,i],each=v))%%p * p**(i-1)
  D=data.frame(LS+1)
  if (r>1) {
    for (i in 2:r) {
      LS=LS[shift]
      D=data.frame(D,LS+1)
    }
  }
  colnames(D)=paste("T", 1:ncol(D), sep = "")
  return(data.frame(Row=rep(1:v,each=v),Col=rep(1:v,v),D))
}

return(molsq(p,q,r,index))
}