#================================================================================
# FONCTIONS DE CALCUL D'UNE FRACTION REGULIERE DE RESOLUTION DONNEE
# Auteur: H. Monod, INRA Jouy en Josas
# Copyright INRA 2009
# Ce programme est une version preliminaire et simplifiee d'une librairie R
# en preparation par A. Kobilinsky, H. Monod, A. Bouvier 
#================================================================================
regular.fraction <- function(s,p,r,resolution){
  # DESCRIPTION
  #  generates a regular fractional factorial design of given resolution,
  #  for s factors at p levels in p^r units
  # ARGUMENTS
  #  s : number of input factors
  #  p : unique prime number of levels of all input and unit factors
  #  r : number of unit factors (so that there are N=p^r units)
  #  resolution : resolution of the fraction
  #  max.sol : maximum number of solutions
  # DETAILS
  #  This is a simplified version of a more general library in preparation.
  #  In this version, all factors must have the same prime number of levels
  #  and only fractions with a given resolution can be constructed. The first
  #  q factors are used as basic factors. The first solution is kept although
  #  it may not be the most interesting one (no control of aberration). This
  #  function is programmed entirely in R and so it is not efficient with respect
  #  to computer time. There is no explicit check on the arguments and so it
  #  is up to the user to restrict p to a prime number such as 2, 3, 5 or 7.
  # OUTPUT:
  #  a list with two components: plan (the design in base p) and matrice.cle
  #  (the design key). The design has N=p^r rows (units) and s columns (factors).
  #  All its elements are integers modulo p that represent the factor levels.

  # ensemble ineligible
  cat("Determination des termes ineligibles: ")
  ineligible <- diag(s)
  for(reso in 2:(resolution-1)){
    combis <- combn(s,reso)
    ncombi <- ncol(combis)
    select <- cbind( c(combis), rep(seq(ncombi),rep(reso,ncombi)) )
    ineli <- matrix(0,s,ncombi)
    ineli[select] <- 1
    ineligible <- cbind(ineligible,ineli)
  }
  cat(ncol(ineligible)," termes ineligibles.\n")
  if( (p!=2) ){
    ineligible <- representative.basep(ineligible,p)
  }
  # Identification of the last non-zero coefficients in each ineligible trt character
  ineligible.lnz <- apply(ineligible, 2, function(x){max(seq(along=x)[x!=0])})
  # initialisation of PhiStar by using the first q factors as basic factors
  PhiStar <- diag(r)
  #
  f <- ncol(PhiStar)
  if(s == f){
    check <- !any(apply(((PhiStar %*% ineligible)%%p)==0, 2, all))
    if(check) return(list(PhiStar))
  }
  # Calculation of the set of initially admissible elements of U*
  admissible <- t(convertinto.basep(seq((p^r)-1),p)) 
  nb.admissible <- ncol(admissible)
  # Backtrack search - preliminaries
  eeU <- list(length=s-f)
  leeU <- rep(NA,s-f)
  neeU <- rep(0,s-f)
  # Backtrack search
  cat("Recherche d'une solution (algorithme backtrack).\n")
  jprev <- 0  ;  j <- 1
  solved <- FALSE
  while((j > 0)&(!solved)){
    PhiStar <- PhiStar[,seq(f+j-1), drop=FALSE]
    if(jprev < j){
      ineligible.j <- ineligible[ seq(f+j-1), ineligible.lnz==(f+j), drop=FALSE ]
      admissible.keep <- planor.kernelcheck.basep(PhiStar, admissible, ineligible.j, p)
      eeU[[j]] <- seq(nb.admissible)[admissible.keep]
      leeU[j] <- length(eeU[[j]])
      neeU[j] <- 0
    }
    if(neeU[j] < leeU[j]){
      neeU[j] <- neeU[j]+1
      newcolj <- (eeU[[j]])[neeU[j]]
      PhiStar <- cbind(PhiStar,admissible[,newcolj])
      if(j == (s-f)){
        cat("Solution obtenue. ")
        solved <- TRUE
        jprev <- j ; j <- j
      }
      else{
        jprev <- j ; j <- j+1
      }
    }
    else{
      jprev <- j ; j <- j-1
    }
  }
  if(solved){
    # Construction du plan
    plan <- crossing(rep(p,r),start=0) %*% PhiStar %%p
    # Sortie
    out <- list(plan=plan, matrice.cle=PhiStar, p=p)
  }
  else{
    cat("Pas de solution. ")
    out <- NULL
  }
  cat("Recherche terminee.\n")
  return(out)
}
#---------------------------------------------------------------------------
planor.kernelcheck.basep <- function(PhiStar, admissible, IneligibleSet, p){
  ImagesIS <- (- PhiStar %*% IneligibleSet)%%p
  avoid <- convertfrom.basep( t(ImagesIS), p)
  candidate <- convertfrom.basep( t(admissible), p)
  test <- !(candidate %in% avoid)
  return(test)
}
#---------------------------------------------------------------------------
convertinto.basep <- function (x, p) {
  # Conversion of an integer or integer vector x into base p
  # The coefficients are ordered by increasing powers of p
  if (!is.numeric(x)) 
    stop("cannot decompose non-numeric arguments")
  if (length(x) > 1) {
    l <- matrix(0, length(x), length(Recall(max(x),p)))
    for(i in seq(along = x)){
      dec.i <- Recall(x[i],p)
      l[i, seq(along=dec.i) ] <- dec.i
    }
    return(l)
  }
  if (x != round(x) || x < 0) 
    return(x)
  val <- x%%p
  while ( (x <- x%/%p) > 0 ) {
    newval <- x%%p
    val <- c(val,newval)
  }
    return(val)
}
#---------------------------------------------------------------------------
convertfrom.basep <- function (x, p) {
  # Conversion of integers x coded as vectors of coefficients in base p
  # to classical integers in base 10
  
  if (!is.numeric(x)) 
    stop("cannot recompose non-numeric arguments")
  if( (max(x)>p) || (min(x)<0) )
    stop("x must be reduced modulo p")
  if (is.matrix(x)) {
    l <- rep(NA, nrow(x))
    for(i in seq(along = l)){
      l[i] <- Recall(x[i,],p)
    }
    return(l)
  }
  val <- sum( x * p^(seq(along=x)-1) )
  return(val)
}
#---------------------------------------------------------------------------
inverses.basep <- function(p){
  # Raw calculation of the inverses modulo p

  if(p==2) return(1)
  else if(p==3) return(c(1,2))
  products <- outer(seq(2,p-2), seq(2,p-2), "*")%%p
  inverses <- 1 + apply(products, 1, function(x){ seq(along=x)[x==1] })
  return( c(1,inverses,p-1) )
}
#---------------------------------------------------------------------------
representative.basep <- function(mat,p){
  # generates the minimal set of representatives in base p
  # of the columns x of matrix mat 

  mat <- as.matrix(mat)
  #
  if(p==2) return(mat %%2)
  #
  representative <- NULL
  for(j in seq(ncol(mat))){
    x <- mat[,j]
    select <- seq(x)[x != 0]
    nbtocross <- length(select)-1
    if( nbtocross <= 0 ) mat.j <- x
    else{
      select <- select[seq(nbtocross)]
      N <- (p-1)^nbtocross
      mat.j <- matrix(x, nrow(mat), N)
      mat.j[select,] <- t( crossing(rep(p-1,nbtocross),start=1) )
    }
    representative <- cbind(representative, mat.j)
  }
  return(representative %%p)
}
#---------------------------------------------------------------------------
crossing <- function(n,start=1){
  # Generates all n1 x n2 x ... x ns combinations of size s with n1,...,ns integers

  N <- prod(n)
  s <- length(n)
  n <- c(n,1)
  crosses <- matrix(NA, N, s)
  for(i in seq(s))
    {
      motif <- start + seq(n[s+1-i])-1
      repet1 <- rep( prod(n[s+1-i+seq(i)]), n[s+1-i] )
      if(i==s){ repet2 <- 1 }
      else{ repet2 <- prod(n[seq(s-i)]) }
      crosses[,s-i+1] <- rep( rep( motif, repet1 ), repet2 )
    }
  return(crosses)
}  
#---------------------------------------------------------------------------