myrect <- function(hor, col = 1, step = 1, vert, steps)
{
  offset <- step / 10
  if(length(vert) == 1) {
    vert <- vert[1]
    if(missing(steps))
      steps <- seq(hor[1], hor[2], by = step)
    rect(steps + offset - step / 2, vert + offset,
         steps - offset + step / 2, vert - offset + step,
         col = gray(col))
  }
  else {
    hor <- hor[1]
    if(missing(steps))
      steps <- seq(vert[1], vert[2], by = step)
    rect(hor + offset - step / 2, steps + offset - step / 2,
         hor - offset + step / 2, steps - offset + step / 2,
         col = gray(col))
  }
}

tmpfn <- function(pos, geno, step, vert = 1,
  fill = c("shade","marker","sample"))
{
  steps <- seq(pos[1], pos[2], by = step)
  nstep <- length(steps)
  fill <- match.arg(fill)
  switch(fill,
    shade = {
      myrect(pos, seq(geno[1], geno[2], length = nstep),
        step, vert, steps)
    },
    marker = {
      myrect(pos, geno, step, vert, pos)
    },
    sample = {
      if(geno[1] == geno[2])
        col <- rep(geno[1], nstep)
      else {
        crossover <- sample(seq(nstep - 1), 1)
        col <- c(rep(geno[1], crossover), rep(geno[2], nstep - crossover))
      }
      myrect(pos, col, step, vert, steps)
    })
}
mygo <- function(vert,fill)
{
  tmpfn(c(0,20),c(0,1),1,vert,fill)
  tmpfn(c(20,30),c(1,1),1,vert,fill)
  tmpfn(c(30,45),c(1,0),1,vert,fill)
  tmpfn(c(45,50),c(0,0),1,vert,fill)
}
########################################################
png("imputegeno.png", width = 960, height = 720)
plot(c(0,50),c(0,12.5), type = "n", bty = "n", yaxt = "n",
  xlab = "position (cM)", ylab = "filled in genotypes")
axis(2, 0.5 + 3:12, 1:10)
mtext("markers", 2, 0, at = 0.5, las = 1, adj = 1)
mtext("genoprob", 2, 0, at = 2, las = 1, adj = 1)
mtext("multiple imputations of genotypes", 3, 1, cex = 1.5)
mygo(0,"marker")
mygo(1.5,"shade")
for(i in 3:12)
  mygo(i,"sample")
dev.off()

png("genarch.png", width = 960, height = 960)
par(mar = c(10,0,4,10))
plot(c(-2,55),c(-5,50), type = "n", bty = "n", yaxt = "n",
  xlab = "position (cM)", ylab = "position (cM)")
axis(4)
mtext("main QTL", 4, 0, at = -2, las = 1, adj = 0, cex = 2)
mtext("main QTL", 1, 0, at = 53, las = 3, adj = 1, cex = 2)
mtext("epistatic pairs of QTL", 3, 0, at = 25, las = 1, adj = 0.6, cex = 3)

## Margins.
myrect(c(0,50), 1, 2, -4)
myrect(53, 1, 2, c(0,50))
## Main QTL.
myrect(, 0, 2, -4, c(4,16,36,50))
myrect(53, 0, 2, c(0,50), c(4,16,36,50))

## Grid.
for(i in seq(0,50, by = 2))
  myrect(c(i+2,50), 1, 2, i - 1)
## Epistatic pairs of QTL.
myrect(, 0, 2, 3, 36)
myrect(, 0, 2, 15, 50)

segments(c(4,4), c(0,4), c(4,34), c(4,4), col = "black", lwd = 4, lty = 2)
segments(c(16,16), c(0,16), c(16,48), c(16,16), col = "black", lwd = 4, lty = 2)

dev.off()