R/orchard.r

In mlstiburek/orchards: Minimum-Inbreeding seed orchard layout

orchMI <- function(ID, GR, GM = NULL, xy = 1, ...) {

set.seed(5)

if (sum(ID[,2]) != max(GR))
  stop("Number of positions on the grid does not correspond to the number of trees.")

MI <- GR #matrix to store the final MI layout
NC <- max(ID[,1]) #total number of clones
NP <- max(GR) #total number of positions
if(is.null(GM)) { #create empty IBD if not provided by the user
  GM <- matrix(0:0, nrow=NC, ncol=NC)
  diag(GM) <- 0.5
}
A <- GM

#PART 1: FLOW MATRIX
#expand genetic relationship by clonal replicates
count <- 1
for (i in 1:NC){
  NR <- ID[i,2]
    if(NR > 1){
      rm <- matrix(0:0, nrow = NR-1, ncol = ncol(A))
      for (j in 1:(NR-1)){
        rm[j,] <- A[count,]
      }
      A <- Append(A, values=rm, after=count, rows = TRUE)
      cm <- matrix(0:0, nrow = nrow(A), ncol = NR-1)
      for (j in 1:(NR-1)){
        cm[,j] <- A[,count]
      }
      A <- Append(A, values=cm, after=count, rows = FALSE)
    }
  count <- count + NR
}
A <- trunc(100*(1-A))
diag(A) <- 0

#PART 2: DISTANCE MATRIX
x <- 1:NP
y <- 1:NP
count <- 1
for(i in seq_len(nrow(GR))){
  for(j in seq_len(ncol(GR))){
    if(GR[i,j] > 0) { #record x and y coordinates of all non-empty positions
      x[count] <- i * xy
      y[count] <- j
      count <- count + 1
    }
  }
}

my_data_points <- data.frame(x,y)
my_distances <- as.matrix(distances(my_data_points)) # Euclidean distances
B <- trunc(100*my_distances)

#find solution of the quadratic assignment problem
a <- qap(A, B, verbose=TRUE, ...)

#format qap solution to the Minimum-Inbreeding layout
count <- 1
for (i in 1:NC){
  NR <- ID[i,2]
  for (j in 1:NR){
    position <- a[count]
    loc <- which(GR == position, TRUE)
    if(!is.null(rownames(ID))){
      MI[loc[1],loc[2]] <- rownames(ID)[i] #use real names
    } else {
      MI[loc[1],loc[2]] <- i #use numbers as clonal codes
    }
    count <- count + 1
  }
}
  MI <- noquote(MI)
  return(MI)
}