#####################################################################
##
## $Id: qb.sim.tour.R,v 1.3.2.3 2006/09/06 18:43:28 byandell Exp $
##
##     Copyright (C) 2006 Brian S. Yandell
##
## This program is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by the
## Free Software Foundation; either version 2, or (at your option) any
## later version.
##
## These functions are distributed in the hope that they will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## The text of the GNU General Public License, version 2, is available
## as http://www.gnu.org/copyleft or by writing to the Free Software
## Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
##
##############################################################################
library(qtlbim)

sim.compare <- function(n.ind = 100, effect = 1, n.qtl = 2,
                        n.mark = 200, bymark = 1, pos = n.mark / 2,
                        remove = TRUE,
                        ... )
{
  ## Simulate cross data with 1 QTL.
  markers <- seq(0, n.mark, by = bymark)
  names(markers) <- paste("M",markers,sep="")
  sim <- sim.cross(list(ch1 = markers),
                   matrix(c(1, pos, -effect), 1, 3),
                   n.ind,
                   "bc" )

  ## Classical R/qtl.
  sim <- calc.genoprob(sim, 1, error = 0)
  sim.scan <- scanone(sim)

  ## Bayesian R/qtlbim.
  sim <- qb.genoprob(sim, step=1)
  cat("\n*** WARNING: object sim created\n\n")
  assign("sim", sim, 1)
  qbSim <- qb.mcmc(sim,
                       chr.nqtl=n.qtl,
                       main.nqtl=1, mean.nqtl=1, max.nqtl=n.qtl,
                       epistasis=FALSE)
  cat("\n*** WARNING: object qbSim created\n\n")
  assign("qbSim", qbSim, 1)
  bone.lod <- qb.scanone(qbSim, type = "LOD", verbose = FALSE)
  bone.add <- qb.scanone(qbSim, type="estimate", scan="add", verbose = FALSE)
  ## Kludge to match effectscan.
  bone.add[seq(bone.add)] <- -bone.add
  
  cat("Classical R/qtl summary\n")
  print(summary(sim.scan))
  cat("Bayesian R/qtlbim summary\n")
  print(summary(bone.lod))
  print(summary(bone.add))

  ## Idealized confounded effect.
  poss <- seq(0, n.mark, by = 1)
  expect <- 2 * effect * exp(-2 * abs(poss - pos) / 100)
  rss <- (n.ind - 2) * (1 + ((2 * effect)^2 - expect^2) / 4)
  lod <- n.ind * log10(((n.ind - 1) * (1 + (2 * effect)^2 / 4)) / rss) / 2
 
  ## LOD scans.
  tmpar <- par( mfrow = c(2,1), mar = c(5.1,4.1,2.1,0.1) )
  on.exit(par(tmpar))
  ylim <- range(bone.lod, sim.scan$lod, lod)
  plot(sim.scan, ylim = ylim, lty = 2, ...)
  mtext("black dash = qtl; red solid = qtlbim; blue dot = ideal 1 QTL",
        1, 2)
  plot(bone.lod, add = TRUE, col = "red",
       verbose = FALSE)
  lines(poss, lod, lty = 3, lwd = 2, col = "blue")

  ## Substitution effect scans.
  ## Classical one QTL.
  ylim <- range( bone.add, 0, 2 * effect )
  effectscan(sim,
       ylim = ylim, lty = 2, lwd = 2, ... )
  axis( 1, at = markers, labels = FALSE )
  mtext("black dash = qtl; red solid = qtlbim; blue dot = ideal 1 QTL",
        1, 2)

  ## Bayesian 1-2 QTL.
  plot(bone.add, add = TRUE, col = "red", verbose = FALSE)

  lines(poss, expect, lty = 3, lwd = 2, col = "blue")

  if(remove)
    qb.remove(qbSim)

  invisible()
}

## ask before plot routines
tmpar <- par(ask = dev.interactive())
on.exit({
  par(tmpar)
})

## Large effect simulation.
sim.compare(n.ind = 100, effect = 1, n.qtl = 2)
## Smaller effect simulation.
sim.compare(n.ind = 200, effect = 0.5, n.qtl = 2)
