######################################################################
#
# argmax.geno.R
#
# copyright (c) 2001, Karl W Broman, Johns Hopkins University
# Sept, 2001; July, 2001; May, 2001; Apr, 2001; Feb, 2001
# Licensed under the GNU General Public License version 2 (June, 1991)
# 
# Part of the R/qtl package
# Contains: argmax.geno
#
######################################################################

######################################################################
#
# argmax.geno: Use Viterbi algorithm to find most likely sequence of
#              underlying genotypes, given observed marker data
#
######################################################################

argmax.geno <-
function(cross, step=0, off.end=0, error.prob=0,
         map.function=c("haldane","kosambi","c-f"))
{

  # map function
  map.function <- match.arg(map.function)
  if(map.function=="kosambi") mf <- mf.k
  else if(map.function=="c-f") mf <- mf.cf
  else mf <- mf.h

  # don't let error.prob be exactly zero, just in case
  if(error.prob < 1e-14) error.prob <- 1e-14

  n.ind <- nind(cross)
  n.chr <- nchr(cross)
  n.mar <- nmar(cross)

  # loop over chromosomes
  for(i in 1:n.chr) {

    # which type of cross is this?
    if(class(cross)[1] == "f2") {
      one.map <- TRUE
      if(class(cross$geno[[i]]) == "A") { # autosomal
        cfunc <- "argmax_geno_f2"
      }
      else {                             # X chromsome 
        cfunc <- "argmax_geno_bc"
      }
    }
    else if(class(cross)[1] == "bc") {
      cfunc <- "argmax_geno_bc"
      one.map <- TRUE
    }
    else if(class(cross)[1] == "4way") {
      cfunc <- "argmax_geno_4way"
      one.map <- FALSE
    }
    else if(class(cross)[1] == "f2ss") {
      cfunc <- "argmax_geno_f2ss"
      one.map <- FALSE
    }
    else stop(paste("argmax.geno not available for cross type",
                    class(cross)[1], "."))

    # genotype data
    gen <- cross$geno[[i]]$data
    gen[is.na(gen)] <- 0
    
    # recombination fractions
    if(one.map) {
      # recombination fractions
      map <- create.map(cross$geno[[i]]$map,step,off.end)
      rf <- mf(diff(map))
      rf[rf < 1e-14] <- 1e-14

      # new genotype matrix with pseudomarkers filled in
      newgen <- matrix(ncol=length(map),nrow=nrow(gen))
      dimnames(newgen) <- list(NULL,names(map))
      newgen[,colnames(gen)] <- gen
      newgen[is.na(newgen)] <- 0
      n.pos <- ncol(newgen)
    }
    else {
      map <- create.map(cross$geno[[i]]$map,step,off.end)
      rf <- mf(diff(map[1,]))
      rf[rf < 1e-14] <- 1e-14
      rf2 <- mf(diff(map[2,]))
      rf2[rf2 < 1e-14] <- 1e-14

      # new genotype matrix with pseudomarkers filled in
      newgen <- matrix(ncol=ncol(map),nrow=nrow(gen))
      dimnames(newgen) <- list(NULL,dimnames(map)[[2]])
      newgen[,colnames(gen)] <- gen
      newgen[is.na(newgen)] <- 0
      n.pos <- ncol(newgen)
    }
    if(any(is.na(rf))) { # this occurs when there is only one marker
      rf <- rf2 <- 0
      warning(paste("Only one marker on chr ", names(cross$geno)[i],
                    ": argmax results tenuous.", sep=""))
    }


    # call the C function
    if(one.map) {
      z <- .C(cfunc,
              as.integer(n.ind),         # number of individuals
              as.integer(n.pos),         # number of markers
              as.integer(newgen),        # genotype data
              as.double(rf),             # recombination fractions
              as.double(error.prob),     
              argmax=as.integer(newgen), # the output
              PACKAGE="qtl")

      cross$geno[[i]]$argmax <- matrix(z$argmax,ncol=n.pos)
      dimnames(cross$geno[[i]]$argmax) <- list(NULL, names(map))
    }
    else {
      z <- .C(cfunc,
              as.integer(n.ind),         # number of individuals
              as.integer(n.pos),         # number of markers
              as.integer(newgen),        # genotype data
              as.double(rf),             # recombination fractions
              as.double(rf2),            # recombination fractions
              as.double(error.prob),      
              argmax=as.integer(newgen), # the output
              PACKAGE="qtl")

      cross$geno[[i]]$argmax <- matrix(z$argmax,ncol=n.pos)
      dimnames(cross$geno[[i]]$argmax) <- list(NULL, colnames(map))
    }

    # attribute set to the error.prob value used, for later
    #     reference
    attr(cross$geno[[i]]$argmax,"error.prob") <- error.prob
    attr(cross$geno[[i]]$argmax,"step") <- step
    attr(cross$geno[[i]]$argmax,"off.end") <- off.end
  }

  cross
}

  

# end of argmax.geno.R
######################################################################
#
# calc.genoprob.R
#
# copyright (c) 2001, Karl W Broman, Johns Hopkins University
# Sept, 2001; July, 2001; May, 2001; Apr, 2001; Feb, 2001
# Licensed under the GNU General Public License version 2 (June, 1991)
# 
# Part of the R/qtl package
# Contains: calc.genoprob
#
######################################################################

######################################################################
#
# calc.genoprob: calculate genotype probabilities conditional on 
#                observed marker genotypes
#
######################################################################

calc.genoprob <-
function(cross, step=0, off.end=0, error.prob=0,
         map.function=c("haldane","kosambi","c-f"))
{

  # map function
  map.function <- match.arg(map.function)
  if(map.function=="kosambi") mf <- mf.k
  else if(map.function=="c-f") mf <- mf.cf
  else mf <- mf.h
 
  # don't let error.prob be exactly zero, just in case
  if(error.prob < 1e-14) error.prob <- 1e-14

  n.ind <- nind(cross)
  n.chr <- nchr(cross)

  # calculate genotype probabilities one chromosome at a time
  for(i in 1:n.chr) {

    # which type of cross is this?
    if(class(cross)[1] == "f2") {
      one.map <- TRUE
      if(class(cross$geno[[i]]) == "A") { # autosomal
        cfunc <- "calc_genoprob_f2"
        n.gen <- 3
        gen.names <- c("A","H","B")
      }
      else {                             # X chromsome 
        cfunc <- "calc_genoprob_bc"
        n.gen <- 2
        gen.names <- c("A","H")
      }
    }
    else if(class(cross)[1] == "bc") {
      cfunc <- "calc_genoprob_bc"
      n.gen <- 2
      gen.names <- c("A","H")
      one.map <- TRUE
    }
    else if(class(cross)[1] == "4way") {
      cfunc <- "calc_genoprob_4way"
      n.gen <- 4
      one.map <- FALSE
      gen.names <- c("AC","AD","BC","BD")
    }
    else stop(paste("calc.genoprob not available for cross type",
                    class(cross)[1], "."))

    # genotype data
    gen <- cross$geno[[i]]$data
    gen[is.na(gen)] <- 0
    
    # recombination fractions
    if(one.map) {
      # recombination fractions
      map <- create.map(cross$geno[[i]]$map,step,off.end)
      rf <- mf(diff(map))
      rf[rf < 1e-14] <- 1e-14

      # new genotype matrix with pseudomarkers filled in
      newgen <- matrix(ncol=length(map),nrow=nrow(gen))
      dimnames(newgen) <- list(NULL,names(map))
      newgen[,colnames(gen)] <- gen
      newgen[is.na(newgen)] <- 0
      n.pos <- ncol(newgen)
    }
    else {
      map <- create.map(cross$geno[[i]]$map,step,off.end)
      rf <- mf(diff(map[1,]))
      rf[rf < 1e-14] <- 1e-14
      rf2 <- mf(diff(map[2,]))
      rf2[rf2 < 1e-14] <- 1e-14

      # new genotype matrix with pseudomarkers filled in
      newgen <- matrix(ncol=ncol(map),nrow=nrow(gen))
      dimnames(newgen) <- list(NULL,dimnames(map)[[2]])
      newgen[,colnames(gen)] <- gen
      newgen[is.na(newgen)] <- 0
      n.pos <- ncol(newgen)
    }


    # call the C function
    if(one.map) {
      z <- .C(cfunc,
              as.integer(n.ind),         # number of individuals
              as.integer(n.pos),         # number of markers
              as.integer(newgen),        # genotype data
              as.double(rf),             # recombination fractions
              as.double(error.prob),     # 
              genoprob=as.double(rep(0,n.gen*n.ind*n.pos)),
              PACKAGE="qtl")

      cross$geno[[i]]$prob <- array(z$genoprob,dim=c(n.ind,n.pos,n.gen))
      dimnames(cross$geno[[i]]$prob) <- list(NULL, names(map), gen.names)
    }
    else {
      z <- .C(cfunc,
              as.integer(n.ind),         # number of individuals
              as.integer(n.pos),         # number of markers
              as.integer(newgen),        # genotype data
              as.double(rf),             # recombination fractions
              as.double(rf2),            # recombination fractions
              as.double(error.prob),     # 
              genoprob=as.double(rep(0,n.gen*n.ind*n.pos)),
              PACKAGE="qtl")

      cross$geno[[i]]$prob <- array(z$genoprob,dim=c(n.ind,n.pos,n.gen))
      dimnames(cross$geno[[i]]$prob) <-
        list(NULL, colnames(map), gen.names)
    }

    # attribute set to the error.prob value used, for later
    #     reference, especially by calc.errorlod()
    attr(cross$geno[[i]]$prob,"error.prob") <- error.prob
    attr(cross$geno[[i]]$prob,"step") <- step
    attr(cross$geno[[i]]$prob,"off.end") <- off.end
  }

  cross
}

# end of calc.genoprob.R
######################################################################
#
# errorlod.R
#
# copyright (c) 2001, Karl W Broman, Johns Hopkins University
# Oct, 2001; Sept, 2001; May, 2001; Apr, 2001
# Licensed under the GNU General Public License version 2 (June, 1991)
# 
# Part of the R/qtl package
# Contains: calc.errorlod, plot.errorlod, top.errorlod
#
######################################################################

######################################################################
#
# calc.errorlod: Calculate LOD scores indicating likely genotyping
#                errors.
#
######################################################################

calc.errorlod <-
function(cross, error.prob=0.01,
         map.function=c("haldane","kosambi","c-f"))
{

  # map function
  map.function <- match.arg(map.function)

  n.ind <- nind(cross)
  n.chr <- nchr(cross)
  n.mar <- nmar(cross)

  if(class(cross)[1]=="bc") cfunc <- "calc_errorlod_bc"
  else if(class(cross)[1]=="f2") cfunc <- "calc_errorlod_f2"
  else if(class(cross)[1]=="4way") cfunc <- "calc_errorlod_4way"
  else stop(paste("calc.errorlod not available for cross type",
                  class(cross)[1],"."))
  
  # calculate genotype probabilities one chromosome at a time
  for(i in 1:n.chr) {

    # skip chromosomes with only 1 marker
    if(n.mar[i] < 2) next

    if(is.na(match("prob",names(cross$geno[[i]])))) {
      # need to run calc.genoprob
      warning("First running calc.genoprob.")
      cross <- calc.genoprob(cross,error.prob=error.prob,
                             map.function=map.function)
      Pr <- cross$geno[[i]]$prob
    }
    else {
      # if error.prob doesn't correspond to what was used when
      #     running calc.genoprob(), re-run calc.genoprob()
      if(abs(attr(cross$geno[[i]]$prob,"error.prob")
             - error.prob) > 1e-9) {
        warning("Re-running calc.genoprob()")
        cross <-
          calc.genoprob(cross,error.prob=error.prob,
                        step=attr(cross$geno[[i]]$prob,"step"),
                        off.end=attr(cross$geno[[i]]$prob,"off.end"),
                        map.function=map.function)
      }
         
      Pr <- cross$geno[[i]]$prob
      u <- grep("^loc\-*[0-9]+",colnames(Pr))
      if(length(u) > 0) Pr <- Pr[,-u,]
    }
    
    nm <- dim(Pr)[2]
    dat <- cross$geno[[i]]$data
    dat[is.na(dat)] <- 0
    
    z <- .C(cfunc,
            as.integer(n.ind),
            as.integer(nm),
            as.integer(dat),
            as.double(error.prob),
            as.double(Pr),
            errlod=as.double(rep(0,n.ind*nm)),
            PACKAGE="qtl")

    errlod <- array(z$errlod, dim=dim(Pr)[1:2])

    dimnames(errlod) <- list(NULL,colnames(cross$geno[[i]]$data))
    cross$geno[[i]]$errorlod <- errlod

    # attribute set to the error.prob value used, for later
    #     reference.
    attr(cross$geno[[i]]$errorlod,"error.prob") <- error.prob
  }

  cross
}

  

  

######################################################################
#
# plot.errorlod
#
######################################################################

plot.errorlod <-
function(x, chr, ...)
{
  cross <- x
  if(!missing(chr)) cross <- pull.chr(cross,chr)

  # remove chromosomes with < 2 markers
  n.mar <- nmar(cross)
  cross <- pull.chr(cross,names(n.mar)[n.mar >= 2])
  n.chr <- nchr(cross)

  errlod <- NULL
  for(i in 1:n.chr) {
    if(is.na(match("errorlod",names(cross$geno[[i]])))) { # need to run calc.errorlod
      warning("First running calc.errorlod.")
      cross <- calc.errorlod(cross,error.prob=0.01,map.function="haldane")
    }
    errlod <- cbind(errlod,cross$geno[[i]]$errorlod)
  }

  errlod <- t(errlod)

  old.xpd <- par("xpd")
  old.las <- par("las")
  par(xpd=TRUE,las=1)
  on.exit(par(xpd=old.xpd,las=old.las))

  colors <- c("white", "gray85", "hotpink", "firebrick")

  # plot grid 
  image(1:nrow(errlod),1:ncol(errlod),errlod,
        ylab="Individuals",xlab="Markers",col=colors,
        breaks=c(-1,1,2,3,max(errlod)))

  # plot lines at the chromosome boundaries
  n.mar <- nmar(cross)
  chr.names <- names(cross$geno)
  a <- c(0.5,cumsum(n.mar)+0.5)

  # the following makes the lines go slightly above the plotting region
  b <- par("usr")
  segments(a,b[3],a,b[4]+diff(b[3:4])*0.02)

  # this line adds a line above and below the image
  #     (the image function seems to leave these out)
  abline(h=0.5+c(0,ncol(errlod)),xpd=FALSE)

  # add chromosome numbers
  a <- par("usr")
  wh <- cumsum(c(0.5,n.mar))
  for(i in 1:length(n.mar)) 
    text(mean(wh[i+c(0,1)]),a[4]+(a[4]-a[3])*0.025,chr.names[i])

  title(main="Genotyping error LOD scores")

}


######################################################################
#
# top.errorlod
#
# Picks out the genotypes having errorlod values above some cutoff
#
######################################################################

top.errorlod <-
function(cross, chr, cutoff=2, msg=TRUE)  
{
  if(!missing(chr)) cross <- pull.chr(cross,chr)

  mar <- ind <- lod <- chr <- NULL

  # remove chromosomes with < 2 markers
  n.mar <- nmar(cross)
  cross <- pull.chr(cross,names(n.mar)[n.mar >= 2])

  flag <- 0
  for(i in 1:nchr(cross)) {
    
    if(is.na(match("errorlod",names(cross$geno[[i]])))) 
      stop("You first need to run calc.errorlod.")

    el <- cross$geno[[i]]$errorlod

    if(any(el > cutoff)) {
      o <- (el > cutoff)
      mar <- c(mar,colnames(el)[col(el)[o]])
      ind <- c(ind,row(el)[o])
      lod <- c(lod,el[o])
      chr <- c(chr,rep(names(cross$geno)[i],sum(o)))
      flag <- 1
    }
  }
  if(!flag) {
    if(msg) cat("\tNo errorlods above cutoff.\n")
    return(invisible(NULL))
  }
  o <- data.frame(chr=chr,ind=ind,marker=mar,errorlod=lod)[order(-lod,ind),]
  rownames(o) <- 1:nrow(o)
  o
}

# end of errorlod.R
######################################################################
#
# est.map.R
#
# copyright (c) 2001, Karl W Broman, Johns Hopkins University
# Sept, 2001; July, 2001; May, 2001; Apr, 2001
# Licensed under the GNU General Public License version 2 (June, 1991)
# 
# Part of the R/qtl package
# Contains: est.map
#
######################################################################

######################################################################
#
# est.map: re-estimate the genetic map for an experimental cross
#
######################################################################

est.map <- 
function(cross, error.prob=0, map.function=c("haldane","kosambi","c-f"),
         maxit=1000, tol=1e-5, sex.sp=TRUE, print.rf=FALSE)
{

  # map function
  map.function <- match.arg(map.function)
  if(map.function=="kosambi") {
    mf <- mf.k; imf <- imf.k
  }
  else if(map.function=="c-f") {
    mf <- mf.cf; imf <- imf.cf
  }
  else {
    mf <- mf.h; imf <- imf.h
  }

  # don't let error.prob be exactly zero, just in case
  if(error.prob < 1e-14) error.prob <- 1e-14

  n.ind <- nind(cross)
  n.mar <- nmar(cross)
  n.chr <- nchr(cross)

  newmap <- vector("list",n.chr)
  names(newmap) <- names(cross$geno)

  # calculate genotype probabilities one chromosome at a time
  for(i in 1:n.chr) {

    if(n.mar[i] < 2) {
      newmap[[i]] <- cross$geno[[i]]$map
      next
    }

    # which type of cross is this?
    if(class(cross)[1] == "f2") {
      one.map <- TRUE
      if(class(cross$geno[[i]]) == "A") # autosomal
        cfunc <- "est_map_f2"
      else                              # X chromsome 
        cfunc <- "est_map_bc"
    }
    else if(class(cross)[1] == "bc") {
      one.map <- TRUE
      cfunc <- "est_map_bc"
    }
    else if(class(cross)[1] == "4way") {
      one.map <- FALSE
      cfunc <- "est_map_4way"
    }
    else if(class(cross)[1] == "f2ss") {
      one.map <- FALSE
      cfunc <- "est_map_f2ss"
    }
    else stop(paste("est.map not available for cross type",
                    class(cross)[1], "."))

    # genotype data
    gen <- cross$geno[[i]]$data
    gen[is.na(gen)] <- 0
    
    # recombination fractions
    if(one.map) {
      # recombination fractions
      rf <- mf(diff(cross$geno[[i]]$map))
      rf[rf < 1e-14] <- 1e-14
    }
    else {
      rf <- mf(diff(cross$geno[[i]]$map[1,]))
      rf[rf < 1e-14] <- 1e-14
      rf2 <- mf(diff(cross$geno[[i]]$map[2,]))
      rf2[rf2 < 1e-14] <- 1e-14
      if(!sex.sp && class(cross$geno[[i]])=="X")
        temp.sex.sp <- TRUE
      else temp.sex.sp <- sex.sp
    }


    # call the C function
    if(one.map) {
      z <- .C(cfunc,
              as.integer(n.ind),         # number of individuals
              as.integer(n.mar[i]),      # number of markers
              as.integer(gen),           # genotype data
              rf=as.double(rf),          # recombination fractions
              as.double(error.prob),     
              loglik=as.double(0),       # log likelihood
              as.integer(maxit),
              as.double(tol),
              as.integer(print.rf),
              PACKAGE="qtl")

      newmap[[i]] <- cumsum(c(min(cross$geno[[i]]$map),imf(z$rf)))
      names(newmap[[i]]) <- names(cross$geno[[i]]$map)
      attr(newmap[[i]],"loglik") <- z$loglik
    }
    else {
      z <- .C(cfunc,
              as.integer(n.ind),         # number of individuals
              as.integer(n.mar[i]),      # number of markers
              as.integer(gen),           # genotype data
              rf=as.double(rf),          # recombination fractions
              rf2=as.double(rf2),        # recombination fractions
              as.double(error.prob),
              loglik=as.double(0),       # log likelihood
              as.integer(maxit),
              as.double(tol),
              as.integer(temp.sex.sp),
              as.integer(print.rf),
              PACKAGE="qtl")
              
      if(!temp.sex.sp) z$rf2 <- z$rf

      newmap[[i]] <- rbind(cumsum(c(min(cross$geno[[i]]$map[1,]),imf(z$rf))),
                           cumsum(c(min(cross$geno[[i]]$map[2,]),imf(z$rf2))))
      dimnames(newmap[[i]]) <- dimnames(cross$geno[[i]]$map)
      attr(newmap[[i]],"loglik") <- z$loglik
    }

  } # end loop over chromosomes

  class(newmap) <- "map"
  newmap
}

# end of est.map.R
######################################################################
#
# est.rf.R
#
# copyright (c) 2001, Karl W Broman, Johns Hopkins University
# Oct, 2001; May, 2001; Apr, 2001
# Licensed under the GNU General Public License version 2 (June, 1991)
# 
# Part of the R/qtl package
# Contains: est.rf, plot.rf
#
######################################################################

######################################################################
#
# est.rf: Estimate sex-averaged recombination fractions between
#         all pairs of markers
#
######################################################################

est.rf <-
function(cross, maxit=1000, tol=1e-6) 
{
  n.chr <- nchr(cross)
  n.mar <- totmar(cross)
  n.ind <- nind(cross)
  mar.names <- unlist(lapply(cross$geno,function(a) colnames(a$data)))
  
  Geno <- NULL
  # create full genotype matrix
  for(i in 1:n.chr) 
    Geno <- cbind(Geno,cross$geno[[i]]$data)

  # which type of cross is this?
  type <- class(cross)[1]
  if(type == "f2") 
    cfunc <- "est_rf_f2"
  else if(type == "bc") 
    cfunc <- "est_rf_bc"
  else if(type == "4way") 
    cfunc <- "est_rf_4way"
  else stop(paste("est.rf not available for cross type",
                  type, "."))

  Geno[is.na(Geno)] <- 0
  
  if(type=="bc")
    z <- .C(cfunc,
            as.integer(n.ind),         # number of individuals
            as.integer(n.mar),         # number of markers
            as.integer(Geno),
            rf = as.double(rep(0,n.mar*n.mar)),
            PACKAGE="qtl")
  else
    z <- .C(cfunc,
            as.integer(n.ind),         # number of individuals
            as.integer(n.mar),         # number of markers
            as.integer(Geno),
            rf = as.double(rep(0,n.mar*n.mar)),
            as.integer(maxit),
            as.double(tol),
            PACKAGE="qtl")

  cross$rf <- matrix(z$rf,ncol=n.mar)
  dimnames(cross$rf) <- list(mar.names,mar.names)
  cross
}

  

plot.rf <-
function(x, chr, which=c("both","lod","rf"), ...)
{
  cross <- x
  which <- match.arg(which)
  
  if(!missing(chr)) 
    cross <- pull.chr(cross,chr)
  
  g <- cross$rf
  if(is.null(g)) 
    stop("You must run est.rf first.")
  
  old.xpd <- par("xpd")
  old.las <- par("las")
  par(xpd=TRUE,las=1)
  on.exit(par(xpd=old.xpd,las=old.las))

  # if any of the rf's are NA (ie no data), put NAs in corresponding LODs
  some.nas <- FALSE
  if(any(is.na(g))) {
    some.nas <- TRUE
    g[is.na(t(g))] <- NA
  }

  if(which=="lod") { # plot LOD scores 
    # copy upper triangle (LODs) to lower triangle (rec fracs)
    g[row(g) > col(g)] <- t(g)[row(g) > col(g)]
    diag(g) <- 12
    g[g>12] <- 12
    br <- c(seq(-1.01,11.01),12)
  }
  else if(which=="rf") { # plot recombination fractions
    # copy lower triangle (rec fracs) to upper triangle (LODs)
    g[row(g) < col(g)] <- t(g)[row(g) < col(g)]
    diag(g) <- 0
    g[g < 0] <- 0
    g[g > 0.5] <- 0.5
    g <- 0.5 - g
    br <- c(seq(0,0.5,length=13)-0.5/12-0.001,0.5)
  }
  else { # plot both LOD scores and recombination fractions
    # rescale lower triangle (rec fracs) to match range of upper triangle
    diag(g) <- 12
    g[row(g) > col(g) & g > 0.5] <- 0.5
    g[row(g) > col(g)] <- (0.5 - g[row(g) > col(g)])*24
    g[g>12] <- 12
    br <- c(seq(-1.01,11.01),12)
  }
  g[is.na(g)] <- min(br)
  diag(g) <- min(br)

  image(1:ncol(g),1:nrow(g),t(g),ylab="Markers",xlab="Markers",breaks=br,
        col=c("gray50",heat.colors(12)))
  
  # plot lines at the chromosome boundaries
  n.mar <- nmar(cross)
  a <- c(0.5,cumsum(n.mar)+0.5)
  abline(v=a,xpd=FALSE)
  abline(h=a,xpd=FALSE)

  # this line adds a line above the image
  #     (the image function leaves it out)
  abline(h=0.5+nrow(g),xpd=FALSE)
  abline(v=0.5+nrow(g),xpd=FALSE)

  # add chromosome numbers
  a <- par("usr")
  wh <- cumsum(c(0.5,n.mar))
  for(i in 1:length(n.mar)) 
    text(mean(wh[i+c(0,1)]),a[4]+(a[4]-a[3])*0.025,names(cross$geno)[i])
  for(i in 1:length(n.mar)) 
    text(a[2]+(a[2]-a[1])*0.025,mean(wh[i+c(0,1)]),names(cross$geno)[i])

  if(which=="lod") title(main="Pairwise LOD scores")
  else if(which=="rf") title(main="Recombination fractions")
  else title("Pairwise recombination fractions and LOD scores")
  
}

# end of est.rf.R
######################################################################
#
# find.errors.R
#
# copyright (c) 2001, Karl W Broman, Johns Hopkins University
# Oct, 2001; July, 2001; Apr, 2001
# Licensed under the GNU General Public License version 2 (June, 1991)
# 
# Part of the R/qtl package
# Contains: find.errors, plot.errors
#
######################################################################

######################################################################
#
# find.errors: uses the results of argmax.geno to identify possible
#              genotyping errors
#
######################################################################

find.errors <-
function(cross, chr, error.prob=0.01,
         map.function=c("haldane","kosambi","c-f"),msg=TRUE)
{
  map.function <- match.arg(map.function)
  if(!missing(chr)) cross <- pull.chr(cross,chr)

  # remove chromosomes with < 2 markers
  n.mar <- nmar(cross)
  cross <- pull.chr(cross,names(n.mar)[n.mar >= 2])

  chr <- ind <- mar <- obs <- am <- NULL

  for(i in 1:nchr(cross)) {
    if(is.na(match("argmax",names(cross$geno[[i]])))) { # haven't yet run argmax
      warning("First running argmax.geno.")
      cross <- argmax.geno(cross,error.prob=error.prob,
                           map.function=map.function)
    }
    else { # pull out only marker loci from argmax information
      u <- grep("^loc\-*[0-9]+",colnames(cross$geno[[i]]$argmax))
      if(any(u))
        cross$geno[[i]]$argmax <- cross$geno[[i]]$argmax[,-u]
    }

    X <- cross$geno[[i]]$data
    Y <- cross$geno[[i]]$argmax

    if(class(cross)[1]=="bc") {
      wh <- (!is.na(X) & X != Y)
      geno <- c("AA","AB")
    }
    else if(class(cross)[1]=="f2") {
      wh <- (!is.na(X) & X != Y & !(X == 4 & (Y==1 | Y==2)) & !(X==5 & (Y==2 | Y==3)))
      geno <- c("AA","AB","BB","AA/AB","AB/BB")
    }
    else if(class(cross)[1]=="4way") {
      wh <- (!is.na(X) & X != Y & !(X == 5 & (Y==1 | Y==3)) & !(X==6 & (Y==2 | Y==4)) &
             !(X==7 & (Y==1 | Y==2)) & !(X==8 & (Y==3 | Y==4)) &
             !(X==9 & (Y==1 | Y==4)) & !(X==10 & (Y==2 | Y==3)))
      geno <- c("AC","BC","AD","BD","AC/AD","BC/BD","AC/BC","AD/BD","AC/BD","AD/BC")
    }
    else
      stop(paste("find.errors not available for cross type", class(cross)[1],"."))
      
    if(any(wh)) {
      chr <- c(chr,rep(names(cross$geno)[i],sum(wh)))
      ind <- c(ind,row(X)[wh])
      mar <- c(mar,colnames(X)[col(X)[wh]]) 
      obs <- c(obs,geno[X[wh]])
      am <- c(am,geno[Y[wh]])
    }
                    
  }
  if(is.null(am)) {
    if(msg) cat("\tNo errors found.\n");
    return(invisible(NULL))
  }

  errs <- data.frame(chr=chr,ind=ind,marker=mar,
                     geno=obs,argmax.geno=am)[order(ind),]
  rownames(errs) <- 1:nrow(errs)
  errs
}

  

plot.errors <-
function(x, chr, ...)
{
  cross <- x
  if(!missing(chr)) cross <- pull.chr(cross,chr)

  # remove chromosomes with < 2 markers
  n.mar <- nmar(cross)
  cross <- pull.chr(cross,names(n.mar)[n.mar >= 2])

  errs <- NULL

  for(i in 1:nchr(cross)) {

    if(is.na(match("argmax",names(cross$geno[[i]])))) { # haven't yet run argmax 
      warning("First running argmax.geno.")
      cross <- argmax.geno(cross)
    }
    else { # pull out only marker loci from argmax information
      u <- grep("^loc\-*[0-9]+",colnames(cross$geno[[i]]$argmax))
      if(any(u))
        cross$geno[[i]]$argmax <- cross$geno[[i]]$argmax[,-u]
    }

    X <- cross$geno[[i]]$data
    Y <- cross$geno[[i]]$argmax

    if(class(cross)[1]=="bc") 
      wh <- (!is.na(X) & X != Y)
    else if(class(cross)[1]=="f2") 
      wh <- (!is.na(X) & X != Y & !(X == 4 & (Y==1 | Y==2)) & !(X==5 & (Y==2 | Y==3)))
    else if(class(cross)[1]=="4way")
      wh <- (!is.na(X) & X != Y & !(X == 5 & (Y==1 | Y==3)) & !(X==6 & (Y==2 | Y==4)) &
             !(X==7 & (Y==1 | Y==2)) & !(X==8 & (Y==3 | Y==4)) &
             !(X==9 & (Y==1 | Y==4)) & !(X==10 & (Y==2 | Y==3)))
    else
      stop(paste("find.errors not available for cross type", class(cross)[1],"."))

    nc <- ncol(X); nr <- nrow(X)
    X <- matrix(1,ncol=nc,nrow=nr)
    if(any(wh)) X[wh] <- NA
    cross$geno[[i]]$data <- X
  }

  plot.missing(cross,main="Likely genotyping errors")
}
  
         
# end of find.errors.R
######################################################################
#
# plot.R
#
# copyright (c) 2000-2001, Karl W Broman, Johns Hopkins University
# Oct, 2001; Sept, 2001; July, 2001; Apr, 2001; Feb, 2001; Mar, 2000
# Licensed under the GNU General Public License version 2 (June, 1991)
# 
# Part of the R/qtl package
# Contains: plot.missing, plot.map, plot.cross, plot.geno, plot.info
#
######################################################################

plot.missing <-
function(x,chr,reorder=FALSE,main="Missing genotypes",...) 
{
  cross <- x
  if(!missing(chr)) cross <- pull.chr(cross,chr)
  
  # get full genotype data into one matrix
  Geno <- cross$geno[[1]]$data
  if(length(cross$geno) > 1) 
    for(i in 2:length(cross$geno))
      Geno <- cbind(Geno,cross$geno[[i]]$data)

  # reorder the individuals according to their phenotype
  o <- 1:nrow(Geno)
  if(reorder) {
    # if reorder is a number, use the corresponding phenotype
    if(is.numeric(reorder)) o <- order(cross$pheno[,reorder])

    # otherwise, order according to the sum of the phenotypes
    else o <- order(apply(cross$pheno,1,sum))
  }

  # make matrix with  0 where genotype data is missing
  #                   1 where data is not missing
  #                 0.5 where data is partially missing
  type <- class(cross)[1]
  g <- t(Geno[o,])
  g[is.na(g)] <- 0
  if(type == "bc") {
    g[g > 0] <- 1
  }
  else if(type=="f2") {
    g[g > 0 & g < 4] <- 1
    g[g > 3] <- 0.5
  }
  else if(type=="4way") {
    g[g > 0 & g < 5] <- 1
    g[g > 4] <- 0.5
  }
  else {
    g[g > 0] <- 1
  }

  old.xpd <- par("xpd")
  old.las <- par("las")
  par(xpd=TRUE,las=1)
  on.exit(par(xpd=old.xpd,las=old.las))

  colors <- c("#000000", "gray80", "#FFFFFF")

  # plot grid with black pixels where there is missing data
  image(1:nrow(g),1:ncol(g),g,ylab="Individuals",xlab="Markers",col=colors,zlim=c(0,1))

  # plot lines at the chromosome boundaries
  n.mar <- nmar(cross)
  a <- c(0.5,cumsum(n.mar)+0.5)
#  abline(v=a)
  # the following makes the lines go slightly above the plotting region
  b <- par("usr")
  segments(a,b[3],a,b[4]+diff(b[3:4])*0.02)

  # this line adds a line above the image
  #     (the image function seems to leave it out)
  abline(h=0.5+c(0,ncol(g)),xpd=FALSE)

  # add chromosome numbers
  a <- par("usr")
  wh <- cumsum(c(0.5,n.mar))
  for(i in 1:length(n.mar)) 
    text(mean(wh[i+c(0,1)]),a[4]+(a[4]-a[3])*0.025,names(cross$geno)[i])

  title(main=main)
}

plot.map <-
function(x,map2,horizontal=FALSE,...) 
{
  map <- x
  # figure out if the input is a cross (containing a map)
  #    or is the map itself
  if(!is.na(match("geno",names(map)))) 
    map <- pull.map(map)

  sex.sp <- FALSE

  if(is.matrix(map[[1]])) { # sex-specific maps
    one.map <- FALSE
    sex.sp <- TRUE
    if(!missing(map2)) {
      if(is.logical(map2)) {
        horizontal <- map2
        map2 <- lapply(map,function(a) a[2,])
        map <- lapply(map,function(a) a[1,])
      }
      else {
        if(!is.na(match("geno",names(map2))))
          map2 <- pull.map(map2)
        Map1 <- lapply(map,function(a) a[1,])
        Map2 <- lapply(map,function(a) a[2,])
        Map3 <- lapply(map2,function(a) a[1,])
        Map4 <- lapply(map2,function(a) a[2,])
        old.mfrow <- par("mfrow")
        on.exit(par(mfrow=old.mfrow))
        par(mfrow=c(2,1))
        plot.map(Map1,Map3,horizontal)
        plot.map(Map2,Map4,horizontal)
        return(invisible())
      }
    }
    else {
      map2 <- lapply(map,function(a) a[2,])
      map <- lapply(map,function(a) a[1,])
    }
  }
  else {
    # determine whether a second map was given
    if(!missing(map2)) {
      if(is.logical(map2)) { # assume "map2" should be "horizontal"
        horizontal <- map2
        map2 <- NULL
        one.map <- TRUE
      }
      else { # determine if it is a cross object
        if(!is.na(match("geno",names(map2))))
          map2 <- pull.map(map2)
        one.map <- FALSE
      }
    }
    else one.map <- TRUE
  }
       
  if(one.map) {
    n.chr <- length(map)
    map <- lapply(map, function(a) a-min(a))
    maxlen <- max(unlist(lapply(map,max)))

    if(horizontal) {
      old.xpd <- par("xpd")
      old.yaxt <- par("yaxt")
      par(xpd=TRUE,yaxt="n")
      on.exit(par(xpd=old.xpd,yaxt=old.yaxt))
      
      plot(0,0,type="n",xlim=c(0,maxlen),ylim=c(0.5,n.chr+0.5),
	   xlab="Location (cM)", ylab="Chromosome")
      a <- par("usr")
      
      for(i in 1:n.chr) {
	lines(c(min(map[[i]]),max(map[[i]])),n.chr+1-c(i,i))
	nmar <- length(map[[i]])
	for(j in 1:nmar)
	  lines(rep(map[[i]][j],2),n.chr+1-i+c(-1/4,1/4))

	# add chromosome label
	text(a[1]-(a[2]-a[1])*0.02,n.chr+1-i,names(map)[i],adj=1)
	lines(c(a[1],a[1]-(a[2]-a[1])*0.01),rep(n.chr+1-i,2))
      }
    }
    else {
      old.xpd <- par("xpd")
      old.xaxt <- par("xaxt")
      old.las <- par("las")
      par(xpd=TRUE,xaxt="n",las=1)
      on.exit(par(xpd=old.xpd,xaxt=old.xaxt,las=old.las))
      
      plot(0,0,type="n",ylim=c(0,maxlen),xlim=c(0.5,n.chr+0.5),
	   ylab="Location (cM)", xlab="Chromosome")
      
      a <- par("usr")
      
      for(i in 1:n.chr) {
	lines(c(i,i), c(min(map[[i]]),max(map[[i]])))
	nmar <- length(map[[i]])
	for(j in 1:nmar)
	  lines(i+c(-1/4,1/4),rep(map[[i]][j],2))

        # add chromosome label
	text(i,a[3]-(a[4]-a[3])*0.04,names(map)[i])
	lines(rep(i,2),c(a[3],a[3]-(a[4]-a[3])*0.02))
      }
    }
    title(main="Genetic map")
  }
  else {
    # check that maps conform
    if(is.matrix(map2[[1]]))
      stop("Second map appears to be a sex-specific map.")
    if(length(map) != length(map2))
      stop("Maps have different numbers of chromosomes.")
    if(any(sapply(map,length) != sapply(map2,length)))
      stop("Maps have different numbers of markers.")

    map1 <- lapply(map,function(a) a-a[1])
    map2 <- lapply(map2,function(a) a-a[1])

    n.chr <- length(map1)
    maxloc <- max(c(unlist(lapply(map1,max)),unlist(lapply(map2,max))))

    if(!horizontal) {
      old.xpd <- par("xpd")
      old.xaxt <- par("xaxt")
      old.las <- par("las")
      par(xpd=TRUE,xaxt="n",las=1)
      on.exit(par(xpd=old.xpd,xaxt=old.xaxt,las=old.las))

      plot(0,0,type="n",ylim=c(0,maxloc),xlim=c(0.5,n.chr+0.5),
           ylab="Location (cM)", xlab="Chromosome")

      a <- par("usr")
    
      for(i in 1:n.chr) {
      
        if(max(map2[[i]]) < max(map1[[i]])) 
          map2[[i]] <- map2[[i]] + (max(map1[[i]])-max(map2[[i]]))/2
        else 
          map1[[i]] <- map1[[i]] + (max(map2[[i]])-max(map1[[i]]))/2
        
        lines(c(i-0.3,i-0.3), c(min(map1[[i]]),max(map1[[i]])))
        lines(c(i+0.3,i+0.3), c(min(map2[[i]]),max(map2[[i]])))
        
        nmar <- length(map1[[i]])
        for(j in 1:nmar)
          lines(c(i-0.3,i+0.3),c(map1[[i]][j],map2[[i]][j]))

        # add chromosome label
        text(i,a[3]-(a[4]-a[3])*0.04,names(map1)[i])
        lines(rep(i,2),c(a[3],a[3]-(a[4]-a[3])*0.02))
      }
    }
    else {
      old.xpd <- par("xpd")
      old.yaxt <- par("yaxt")
      old.las <- par("las")
      par(xpd=TRUE,yaxt="n",las=1)
      on.exit(par(xpd=old.xpd,yaxt=old.yaxt,las=old.las))

      plot(0,0,type="n",xlim=c(0,maxloc),ylim=c(0.5,n.chr+0.5),
           xlab="Location (cM)", ylab="Chromosome")

      a <- par("usr")
    
      for(i in 1:n.chr) {
      
        if(max(map2[[i]]) < max(map1[[i]])) 
          map2[[i]] <- map2[[i]] + (max(map1[[i]])-max(map2[[i]]))/2
        else 
          map1[[i]] <- map1[[i]] + (max(map2[[i]])-max(map1[[i]]))/2
        
        lines(c(min(map2[[i]]),max(map2[[i]])), c(n.chr-i-0.3+1,n.chr-i+1-0.3))
        lines(c(min(map1[[i]]),max(map1[[i]])), c(n.chr-i+1+0.3,n.chr+1-i+0.3))
        
        nmar <- length(map1[[i]])
        for(j in 1:nmar)
          lines(c(map2[[i]][j],map1[[i]][j]), c(n.chr+1-i-0.3,n.chr+1-i+0.3))

        # add chromosome label
        text(a[1]-diff(a[1:2])*0.04,n.chr+1-i, names(map1)[i])
        lines(c(a[1],a[1]-diff(a[1:2])*0.02), rep(n.chr+1-i,2))
      }

    }
    if(!sex.sp) title(main="Comparison of genetic maps")
    else title(main="Genetic map")
  }    

}


plot.cross <-
function(x,...)
{
  cross <- x
  
  old.yaxt <- par("yaxt")
  if(ncol(cross$pheno) > 2) {
    old.ask <- par("ask")
    par(ask=TRUE)
    on.exit(par(ask=old.ask,yaxt=old.yaxt))
  }
  else {
    old.mfrow <- par("mfrow")
    par(mfrow=c(2,2))
    on.exit(par(mfrow=old.mfrow,yaxt=old.yaxt))
  }

  plot.missing(cross)
  plot.map(cross)
  par(yaxt = "n")
  for(i in 1:ncol(cross$pheno)) 
    hist(cross$pheno[,i],nclass=round(sqrt(nrow(cross$pheno))+5),
	 xlab=colnames(cross$pheno)[i],prob=TRUE, ylab="",
	 main=paste("Histogram of", colnames(cross$pheno)[i]))
}


##################################################r####################
#
# plot.geno: Plot genotypes for a specified chromosome, with likely
#           genotyping errors indicated. 
#
######################################################################

plot.geno <-
function(x, chr, ind, horizontal=FALSE, cutoff=2,
         method=c("lod","argmax"),min.sep=1,...)
{
  cross <- x  
  method <- match.arg(method)
  cross <- pull.chr(cross,chr)
  type <- class(cross)[1]
  
  if(type != "bc" && type != "f2")
    stop("This function has only been coded for bc and f2 crosses.")

  if(method=="lod") {
    if(is.na(match("errorlod",names(cross$geno[[1]])))) {
      warning("First running calc.errorlod.")
      cross <- calc.errorlod(cross,error.prob=0.01)
    }
  }
  else {
    if(is.na(match("argmax",names(cross$geno[[1]])))) {
      warning("First running argmax.geno.")
      cross <- argmax.geno(cross,error.prob=0.01)
    }
  }
  
  # if necessary, discard parts of argmax that are not at markers
  if(method=="argmax") {
    wh <- grep("^loc\-*[0-9]+",colnames(cross$geno[[1]]$argmax))
    if(length(wh) > 0) 
      cross$geno[[1]]$argmax <- cross$geno[[1]]$argmax[,-wh]
  }

  # indicators for apparent errors
  errors <- matrix(0,ncol=ncol(cross$geno[[1]]$data),
                   nrow=nrow(cross$geno[[1]]$data))
  dimnames(errors) <- dimnames(cross$geno[[1]]$data)

  if(method=="lod") top <- top.errorlod(cross,1,cutoff,FALSE)
  else top <- find.errors(cross,1,msg=FALSE)
  if(length(top) > 0)
    for(i in 1:nrow(top))
      errors[top[i,2],as.character(top[i,3])] <- 1

  # map, data, errors
  map <- cross$geno[[1]]$map
  if(is.matrix(map)) map <- map[1,] # if sex-specific map
  L <- diff(range(map))
  min.d <- L*min.sep/100
  d <- diff(map)
  d[d < min.d] <- min.d
  map <- cumsum(c(0,d))

  data <- cross$geno[[1]]$data
  if(!missing(ind)) {
    data <- data[ind,]
    errors <- errors[ind,]
  }
  n.ind <- nrow(errors)

  color <- c("white","gray60","black","green","orange","red")

  if(horizontal==TRUE) {
    plot(0,0,type="n",xlab="Position (cM)",ylab="Individual",
         main=paste("Chromosome",names(cross$geno)[1]),
         ylim=c(0.5,n.ind+0.5),xlim=c(0,max(map)))
    segments(0,1:n.ind,max(map),1:n.ind)

    # AA genotypes
    tind <- rep(1:n.ind,length(map));tind[is.na(data)] <- NA
    ind <- tind; ind[!is.na(data) & data!=1] <- NA
    x <- rep(map,rep(n.ind,length(map)))
    points(x,ind,pch=16,col=color[1])
    points(x,ind,pch=1)

    # AB genotypes
    ind <- tind; ind[!is.na(data) & data!=2] <- NA
    if(type=="bc")
      points(x,ind,pch=16,col=color[3]) 
    else {
      points(x,ind,pch=16,col=color[2])
      points(x,ind,pch=1)
    }

    if(type=="f2") {
      # BB genotypes
      ind <- tind; ind[!is.na(data) & data!=3] <- NA
      points(x,ind,pch=16,col=color[3])

      # not BB (D in mapmaker/qtl) genotypes
      ind <- tind; ind[!is.na(data) & data!=4] <- NA
      points(x,ind,pch=16,col=color[4])
      points(x,ind,pch=1)

      # not AA (C in mapmaker/qtl) genotypes
      ind <- tind; ind[!is.na(data) & data!=5] <- NA
      points(x,ind,pch=16,col=color[5])
      points(x,ind,pch=1)
    }

    # plot map
    u <- par("usr")
    segments(map,u[3],map,(u[3]+1)/2)
    segments(map,u[4],map,(n.ind+u[4])/2)

    if(any(errors)) {
      ind <- rep(1:n.ind,length(map));ind[errors!=1]<-NA
      points(x,ind,pch=0,col=color[6],cex=1.5)
    }

  }
  else {
    plot(0,0,type="n",ylab="Position (cM)",xlab="Individual",
         main=paste("Chromosome",names(cross$geno)[1]),
         xlim=c(0.5,n.ind+0.5),ylim=c(0,max(map)))
    segments(1:n.ind,0,1:n.ind,max(map))
    
    # AA genotypes
    tind <- rep(1:n.ind,length(map));tind[is.na(data)] <- NA
    ind <- tind; ind[!is.na(data) & data!=1] <- NA
    y <- rep(map,rep(n.ind,length(map)))
    points(ind,y,pch=16,col="white")
    points(ind,y,pch=1)

    # AB genotypes
    ind <- tind; ind[!is.na(data) & data!=2] <- NA
    if(type=="bc")
      points(ind,y,pch=16,col=color[3])
    else {
      points(ind,y,pch=16,col=color[2])
      points(ind,y,pch=1)
    }

    if(type=="f2") {
      # BB genotypes
      ind <- tind; ind[!is.na(data) & data!=3] <- NA
      points(ind,y,pch=16,col=color[3])

      # not BB genotypes
      ind <- tind; ind[!is.na(data) & data!=4] <- NA
      points(ind,y,pch=16,col=color[4])
      points(ind,y,pch=1)

      # not AA genotypes
      ind <- tind; ind[!is.na(data) & data!=5] <- NA
      points(ind,y,pch=16,col=color[5])
      points(ind,y,pch=1)
    }

    # plot map
    u <- par("usr")
    segments(u[1],map,(u[1]+1)/2,map)
    segments(u[2],map,(n.ind+u[2])/2,map)

    if(any(errors)) {
      ind <- rep(1:n.ind,length(map));ind[errors!=1]<-NA
      points(ind,y,pch=0,col=color[6],cex=1.5)
    }
  }
}
    
######################################################################
#
# plot.info: Plot the proportion of missing information in the
#            genotype data.
#
######################################################################

plot.info <-
function(x,chr,which=c("both","entropy","variance"),return.result=FALSE,...)
{
  cross <- x
  which <- match(match.arg(which),c("entropy","variance","both"))-1

  if(!missing(chr)) cross <- pull.chr(cross,chr)
  results <- NULL

  n.chr <- nchr(cross)
  if(is.na(match("prob",names(cross$geno[[1]])))) { # need to run calc.genoprob
    warning("First running calc.genoprob.")
    cross <- calc.genoprob(cross)
  }
  gap <- attr(cross$geno[[1]]$prob,"off.end")*2+10 # gap between chr in plot
  n.ind <- nind(cross)
  for(i in 1:n.chr) {
    n.gen <- dim(cross$geno[[i]]$prob)[3]
    n.pos <- ncol(cross$geno[[i]]$prob)

    # calculate information (between 0 and 1)
    info <- .C("R_info",
               as.integer(n.ind),
               as.integer(n.pos),
               as.integer(n.gen),
               as.double(cross$geno[[i]]$prob),
               info1=as.double(rep(0,n.pos)),
               info2=as.double(rep(0,n.pos)),
               as.integer(which))

    if(which != 1) { # rescale entropy version
      if(n.gen==3) maxent <- 1.5*log(2)
      else maxent <- log(n.gen)
      info$info1 <- -info$info1/maxent
    }
    if(which != 0) { # rescale variance version
      maxvar <- c(0.25,0.5,1.25)[n.gen-1]
      info$info2 <- info$info2/maxvar
    }

    # reconstruct map
    map <- create.map(cross$geno[[i]]$map,
                      attr(cross$geno[[i]]$prob,"step"),
                      attr(cross$geno[[i]]$prob,"off.end"))
    if(is.matrix(map)) map <- map[1,]

    z <- data.frame(chr=rep(names(cross$geno)[i],length(map)),pos=map,
                    "Missing information"=info$info1,
                    "Missing information"=info$info2)
    w <- names(map)
    o <- grep("^loc\-*[0-9]+",w)
    if(length(o) > 0) 
      w[o] <- paste(w[o],names(cross$geno)[i],sep=".c")
    rownames(z) <- w
    results <- rbind(results,z)
  }

  if(which==0) plot.scanone(results,ylim=c(0,1),gap=gap)
  else if(which==1) plot.scanone(results[,-3],ylim=c(0,1),gap=gap)
  else if(which==2) plot.scanone(results,results[,-3],ylim=c(0,1),gap=gap)
  
  if(return.result) {
    colnames(results)[3:4] <- c("misinfo.entropy","misinfo.variance")
    class(results) <- c("scanone","data.frame")
    return(results)
  }
  else invisible()
}

# end of plot.R
######################################################################
#
# read.cross.R
#
# copyright (c) 2000-2001, Karl W Broman, Johns Hopkins University
# Oct, 2001; Sept, 2001; Feb, 2001; Sept, 2000; Aug, 2000 
# Licensed under the GNU General Public License version 2 (June, 1991)
#
# Part of the R/qtl package
# Contains: read.cross, read.cross.karl, read.cross.mm,
#           read.cross.gary, read.cross.csv
#
######################################################################


######################################################################
#
# read.cross: read data from an experimental cross
#
######################################################################

read.cross <-
function(format=c("csv","mm","gary","karl"),...)
{
  format <- match.arg(format)
  
  if(format=="karl") # karl's format
    cross <- read.cross.karl(...)
  if(format=="mm") # mapmaker format
    cross <- read.cross.mm(...)
  if(format=="gary") # gary's format
    cross <- read.cross.gary(...)
  if(format=="csv") # comma-delimited format
    cross <- read.cross.csv(...)
  
  cross
}




######################################################################
#
# read.cross.karl
#
# read data in Karl's format
#
######################################################################

read.cross.karl <-
function(dir,genfile,mapfile,phefile)  
{
  # create file names 
  if(missing(genfile)) genfile <- "gen.txt"
  if(missing(mapfile)) mapfile <- "map.txt"
  if(missing(phefile)) phefile <- "phe.txt"

  if(!missing(dir)) {
    # remove ending "/" if it exists
    n <- nchar(dir)
    if(substr(dir,n,n) == "/")
      dir <- substr(dir,0,n-1)
    genfile <- paste(dir,genfile, sep="/")
    mapfile <- paste(dir,mapfile, sep="/")
    phefile <- paste(dir,phefile, sep="/")
  }

  # read data
  geno <- as.matrix(read.table(genfile,na.strings="0"))
  pheno <- as.matrix(read.table(phefile,na.strings="-",header=TRUE))
  tempmap <- scan(mapfile, what=character(),quiet=TRUE)

  # fix up map information
  # number of chromosomes
  n.chr <- as.numeric(tempmap[1])
  n.mar <- 1:n.chr
  g <- map <- geno.data <- vector("list", n.chr)
  cur <- 2
  min.mar <- 1
  names(g) <- as.character(1:n.chr)
  for(i in 1:n.chr) { # loop over chromosomes
    # number of markers
    n.mar[i] <- as.numeric(tempmap[cur])
    cur <- cur+1

    # pull out appropriate portion of genotype data
    geno.data[[i]] <- geno[,min.mar:(min.mar+n.mar[i]-1)]
    min.mar <- min.mar + n.mar[i]

    # recombination fractions
    r <- as.numeric(tempmap[cur:(cur+n.mar[i]-2)])

    # convert to cM distances (w/ Kosambi map function)
    d <- 0.25*log((1+2*r)/(1-2*r))*100

    # convert to locations
    map[[i]] <- round(c(0,cumsum(d)),2)
    cur <- cur+n.mar[i]-1

    # marker names
    names(map[[i]]) <- tempmap[cur:(cur+n.mar[i]-1)]
    dimnames(geno.data[[i]]) <- list(NULL, names(map[[i]]))
    cur <- cur+n.mar[i]

    g[[i]] <- list(data=geno.data[[i]],map=map[[i]])

    # attempt to pull out chromosome number
    mar.names <- names(map[[i]])
    twodig <- grep("[Dd][1-9][0-9][Mm]", mar.names)
    onedig <- grep("[Dd][1-9][Mm]", mar.names)
    xchr <- grep("[Dd][Xx][Mm]", mar.names)

    chr.num <- NULL
    if(length(twodig) > 0)
      chr.num <- c(chr.num,substr(mar.names[twodig],2,3))
    if(length(onedig) > 0)
      chr.num <- c(chr.num,substr(mar.names[onedig],2,2))
    if(length(xchr) > 0)
      chr.num <- c(chr.num,rep("X",length(xchr)))

    # no marker names of the form above
    if(is.null(chr.num)) {
      chr.num <- length(mar.names)
      names(chr.num) <- "1"
    }
    else {
      chr.num <- table(chr.num)
    }
  
    m <- max(chr.num)
    if(m > sum(chr.num)/2 && m > 1) 
      names(g)[i] <- names(chr.num)[chr.num==m][1]

    if(names(g)[i] == "X" || names(g)[i] == "x") class(g[[i]]) <- "X"
    else class(g[[i]]) <- "A"
  }

  # check that data dimensions match 
  n.mar1 <- sapply(g,function(a) ncol(a$data))
  n.mar2 <- sapply(g,function(a) length(a$map))
  n.phe <- ncol(pheno)
  n.ind1 <- nrow(pheno)
  n.ind2 <- sapply(g,function(a) nrow(a$data))
  if(any(n.ind1 != n.ind2)) {
    print(c(n.ind1,n.ind2))
    stop("Number of individuals in genotypes and phenotypes do not match.");
  }
  if(any(n.mar1 != n.mar2)) {
    print(c(n.mar,n.mar2))
    stop("Numbers of markers in genotypes and marker names files do not match.");
  }

  # print some information about the amount of data read
  cat(" --Read the following data:\n");
  cat("\t", n.ind1, " individuals\n");
  cat("\t", sum(n.mar1), " markers\n");
  cat("\t", n.phe, " phenotypes\n");

  # add phenotype names, if missing
  if(is.null(colnames(pheno))) 
    dimnames(pheno) <- list(NULL, paste("phenotype", 1:n.phe,sep=""))

  # determine map type: f2 or bc or 4way?
  if(max(geno[!is.na(geno)])<=2) type <- "bc"
  else if(max(geno[!is.na(geno)])<=5) type <- "f2"
  else type <- "4way"
  cross <- list(geno=g,pheno=pheno)
  class(cross) <- c(type,"cross")

  # check that nothing is strange in the genotype data
  cross.type <- class(cross)[1]
  if(cross.type=="f2") max.gen <- 5
  else if(cross.type=="bc") max.gen <- 2
  else max.gen <- 10

  u <- unique(geno)
  if(any(!is.na(u) & (u > max.gen | u < 1))) 
    stop(paste("There are stange values in the genotype data :",
               paste(u,collapse=":"), "."))

  cross
}
  

######################################################################
#
# read.cross.mm: read data from an experimental cross in mapmaker
#                format.
#
# We need two files: a "raw" file containing the genotype and
# phenotype data and a "map" file containing the chromosomes
# assignments and (optionally) map positions.
#
# The map file contains two or three columns, separated by white
# space, with the chromosome number, marker name (with markers in
# order along the chromosomes) and (optionally) the map position.
#
######################################################################

read.cross.mm <-
function(dir,rawfile,mapfile,estimate.map=TRUE)
{
  # create file names
  if(missing(mapfile))
    stop("You must specify the mapfile containing the chromosome assignments of markers.")
  if(missing(rawfile))
    stop("You must specify the rawfile containing the genotype and phenotype data.")
  if(!missing(dir)) {
    # remove ending "/" if it exists
    n <- nchar(dir)
    if(substring(dir,n,n) == "/")
      dir <- substring(dir,0,n-1)
    mapfile <- paste(dir, mapfile, sep="/")
    rawfile <- paste(dir, rawfile, sep="/")
  }

  # count lines in rawfile
  n.lines <- length(scan(rawfile, what=character(), skip=0, nlines=0,
                         blank.lines.skip=FALSE,quiet=TRUE))
  map <- read.table(mapfile,header=FALSE,as.is=TRUE,blank=FALSE)
  o <- (1:nrow(map))[map[,1]==""]
  if(length(o) > 0) map <- map[-o,]

  # remove any leading *'s from the marker names
  g <- grep("^*",map[,2],extended=FALSE)
  if(length(g) > 0) 
    map[g,2] <- substr(map[g,2],2,nchar(map[g,2]))

  cur.mar <- 0
  cur.phe <- 0

  NEW.symb <- c("1","2","3","4","5","0")
  OLD.symb <- c("A","H","B","D","C","-")

  flag <- 0
  for(i in 1:n.lines) {
    a <- scan(rawfile,what=character(),skip=i-1,nlines=1,
              blank.lines.skip=TRUE,quiet=TRUE)

    if(length(a) == 0) next

    if(length(grep("#", a[1])) != 0) next

    if(flag == 0) {
      flag <- 1
      type <- a[length(a)]
      if(type == "intercross") type <- "f2"
      else if(type == "backcross") type <- "bc"
      else
        stop(paste("File indicates invalid cross type: ", type,
                   ".", sep=""))
    }
    else if(flag == 1) {
      flag <- 2
      n.ind <- as.numeric(a[1])
      n.mar <- as.numeric(a[2])
      n.phe <- as.numeric(a[3])
      cat("\tType of cross:         ", type, "\n")
      cat("\tNumber of individuals: ", n.ind, "\n")
      cat("\tNumber of markers:     ", n.mar, "\n")
      cat("\tNumber of phenotypes:  ", n.phe, "\n")
      convert.symb <- FALSE
      if(length(a) > 3 && a[4] == "symbols") {
        b <- a[-(1:4)]
        convert.symb <- TRUE
        new.symb <- substring(b,1,1)
        old.symb <- substring(b,3,3)
      }

      marnames <- rep("", n.mar)
      geno <- matrix(0,ncol=n.mar,nrow=n.ind)
      if(n.phe == 0) {
        pheno <- matrix(1:n.ind,ncol=1)
        phenames <- c("number")
      }
      else {
        pheno <- matrix(0,ncol=n.phe,nrow=n.ind)
        phenames <- rep("", n.phe)
      }

    }
    else {

      if(substring(a[1],1,1) == "*") {
        cur.mar <- cur.mar+1
        cur.row <- 1

        if(cur.mar > n.mar) { # now reading phenotypes
          cur.phe <- cur.phe+1
          if(cur.phe > n.phe) next 
          phenames[cur.phe] <- substring(a[1],2)
          p <- a[-1]
          p[p=="-"] <- NA
          n <- length(p)
          pheno[cur.row+(0:(n-1)),cur.phe] <- as.numeric(p)
          cur.row <- cur.row + n
        }

        else { # reading genotypes
          marnames[cur.mar] <- substring(a[1],2)
          g <- paste(a[-1],collapse="")
          g <- unlist(strsplit(g,""))

          if(convert.symb) {
            h <- g
            for(j in 1:length(new.symb)) {
              if(any(h==new.symb[j]))
                g[h==new.symb[j]] <- old.symb[j]
            }
          }

          h <- g
          for(j in 1:length(NEW.symb)) {
            if(any(h==OLD.symb[j]))
              g[h==OLD.symb[j]] <- NEW.symb[j]
          }

          n <- length(g)

          geno[cur.row+(0:(n-1)),cur.mar] <- as.numeric(g)
          cur.row <- cur.row + n
        }

      }
      else { # continuation lines
        if(cur.mar > n.mar) { # now reading phenotypes
          a[a=="-"] <- NA
          n <- length(a)
          pheno[cur.row+(0:(n-1)),cur.phe] <- as.numeric(a)
          cur.row <- cur.row + n
        }
        else {
          g <- paste(a,collapse="")
          g <- unlist(strsplit(g,""))
          if(convert.symb) {
            h <- g
            for(j in 1:length(new.symb)) {
              if(any(h==new.symb[j]))
                g[h==new.symb[j]] <- old.symb[j]
            }
          }
          h <- g
          for(j in 1:length(NEW.symb)) {
            if(any(h==OLD.symb[j]))
              g[h==OLD.symb[j]] <- NEW.symb[j]
          }
          n <- length(g)
          geno[cur.row+(0:(n-1)),cur.mar] <- as.numeric(g)
          cur.row <- cur.row + n
        }

      }          

    }
  }
  dimnames(pheno) <- list(NULL, phenames)

  # parse map file
  includes.pos <- FALSE
  if(ncol(map) == 3) includes.pos <- TRUE
  
  chr <- as.character(map[,1])
  markers <- map[,2]
  if(includes.pos) pos <- map[,3]

  Geno <- vector("list",length(unique(chr)))
  names(Geno) <- unique(chr)

  for(i in unique(chr)) {
    mar <- markers[chr == i]

    # create map
    if(includes.pos) map <- pos[chr == i]
    else map <- seq(0,by=5,length=length(mar))
    names(map) <- mar
    
    # pull out genotype data
    o <- match(mar,marnames)
    if(any(is.na(o))) {
      stop(paste("Cannot find markers in genotype data: ",
           paste(mar[is.na(o)],collapse=" "), ".",sep=""))
    }

    if(length(o)==1) data <- matrix(geno[,o],ncol=1)
    else data <- geno[,o]
    # add marker names to data
    colnames(data) <- mar
    # changes 0's to NA's
    data[!is.na(data) & data==0] <- NA

    Geno[[i]] <- list(data=data,map=map)
    if(i=="X" || i=="x") class(Geno[[i]]) <- "X"
    else class(Geno[[i]]) <- "A"
  }

  cross <- list(geno=Geno,pheno=pheno)
  class(cross) <- c(type,"cross")

  if(!includes.pos && estimate.map) {
    cat(" --Estimating genetic map\n")
    newmap <- est.map(cross)
    cross <- replace.map(cross,newmap)
  }

  cross
}


######################################################################
#
# read.cross.gary
#
# read data in Gary's format
#
######################################################################

read.cross.gary <-
function(dir,genfile,mnamesfile,chridfile,phefile,pnamesfile,mapfile)
{
  # create file names 
  if(missing(genfile)) genfile <- "geno.dat"
  if(missing(mnamesfile)) mnamesfile <- "mnames.txt"
  if(missing(chridfile)) chridfile <- "chrid.dat"
  if(missing(phefile)) phefile <- "pheno.dat"
  if(missing(pnamesfile)) pnamesfile <- "pnames.txt"
  if(missing(mapfile)) mapfile <- "markerpos.txt"

  if(!missing(dir)) {
    # remove ending "/" if it exists
    n <- nchar(dir)
    if(substr(dir,n,n) == "/")
      dir <- substr(dir,0,n-1)
    genfile <- paste(dir,genfile, sep="/")
    mnamesfile <- paste(dir,mnamesfile, sep="/")
    chridfile <- paste(dir,chridfile, sep="/")
    phefile <- paste(dir,phefile, sep="/")
    pnamesfile <- paste(dir,pnamesfile, sep="/")
    mapfile <- paste(dir,mapfile, sep="/")
  }

  # read data
  allgeno <- as.matrix(read.table(genfile,na.strings="9"))+1
  pheno <- as.matrix(read.table(phefile,na.strings="-",header=FALSE))
  chr <- scan(chridfile,what=character(),quiet=TRUE)
  mnames <- scan(mnamesfile,what=character(),quiet=TRUE)
  map <- read.table(mapfile,row.names=1)
  pnames <- scan(pnamesfile,what=character(),quiet=TRUE)

  map <- map[mnames,1]

  # fix up map information
  # number of chromosomes
  uchr <- unique(chr)
  n.chr <- length(uchr)
  geno <- vector("list",n.chr)
  names(geno) <- uchr
  min.mar <- 1
  for(i in 1:n.chr) { # loop over chromosomes
    # create map
    temp.map <- map[chr==uchr[i]]

    # deal with any markers that didn't appear in the marker pos file
    if(any(is.na(temp.map))) {
      o <- (1:length(temp.map))[is.na(temp.map)]
      for(j in o) {
        if(j==1 || all(is.na(temp.map[1:(j-1)]))) {
          z <- min((1:length(temp.map))[-o])
          temp.map[j] <- min(temp.map,na.rm=TRUE)-(z-j+1)
        }
        else if(j==length(temp.map) || all(is.na(temp.map[-(1:j)]))) {
          z <- max((1:length(temp.map))[-o])
          temp.map[j] <- max(temp.map,na.rm=TRUE)+(j-z+1)
        }
        else {
          temp.map[j] <- (min(temp.map[-(1:j)],na.rm=TRUE)+
                          max(temp.map[1:(j-1)],na.rm=TRUE))/2
        }
      }
    }
    
    names(temp.map) <- mnames[chr==uchr[i]]

    # pull out appropriate portion of genotype data
    data <- allgeno[,min.mar:(length(temp.map)+min.mar-1)]
    min.mar <- min.mar + length(temp.map)
    colnames(data) <- names(temp.map)

    geno[[i]] <- list(data=data,map=temp.map)
    if(uchr[i] == "X" || uchr[i] == "x") 
      class(geno[[i]]) <- "X"
    else class(geno[[i]]) <- "A"
  }
  colnames(pheno) <- pnames

  # check that data dimensions match 
  n.mar1 <- sapply(geno,function(a) ncol(a$data))
  n.mar2 <- sapply(geno,function(a) length(a$map))
  n.phe <- ncol(pheno)
  n.ind1 <- nrow(pheno)
  n.ind2 <- sapply(geno,function(a) nrow(a$data))
  if(any(n.ind1 != n.ind2)) {
    print(c(n.ind1,n.ind2))
    stop("Number of individuals in genotypes and phenotypes do not match.");
  }
  if(any(n.mar1 != n.mar2)) {
    print(c(n.mar,n.mar2))
    stop("Numbers of markers in genotypes and marker names files do not match.");
  }

  # print some information about the amount of data read
  cat(" --Read the following data:\n");
  cat("\t", n.ind1, " individuals\n");
  cat("\t", sum(n.mar1), " markers\n");
  cat("\t", n.phe, " phenotypes\n");

  # determine map type: f2 or bc or 4way?
  if(max(allgeno[!is.na(allgeno)])<=2) type <- "bc"
  else type <- "f2"
  cross <- list(geno=geno,pheno=pheno)
  class(cross) <- c(type,"cross")

  # check that nothing is strange in the genotype data
  cross.type <- class(cross)[1]
  if(cross.type=="f2") max.gen <- 5
  else max.gen <- 2

  u <- unique(allgeno)
  if(any(!is.na(u) & (u > max.gen | u < 1))) 
    stop(paste("There are stange values in the genotype data :",
               paste(sort(u),collapse=":"), "."))

  cross
}
  

######################################################################
#
# read.cross.csv
#
# read data in comma-delimited format
#
######################################################################

read.cross.csv <-
function(dir,file,sep=",",na.strings="-",genotypes=c("A","H","B","C","D"))
{
  # create file names 
  if(missing(file)) file <- "data.csv"

  if(!missing(dir)) {
    # remove ending "/" if it exists
    n <- nchar(dir)
    if(substr(dir,n,n) == "/")
      dir <- substr(dir,0,n-1)
    file <- paste(dir,file, sep="/")
  }

  # read data
  data <- as.matrix(read.table(file,sep=sep,na.strings=na.strings,
                               as.is=TRUE))
  data[!is.na(data) & data == ""] <- NA

  # pull apart phenotypes, genotypes and map
  n.phe <- max((1:ncol(data))[is.na(data[2,])])
  pheno <- matrix(as.numeric(data[-(1:3),1:n.phe]),nrow=nrow(data)-3)
  colnames(pheno) <- data[1,1:n.phe]
  mnames <- data[1,-(1:n.phe)]
  chr <- data[2,-(1:n.phe)]
  map <- as.numeric(data[3,-(1:n.phe)])
  allgeno <- matrix(match(data[-(1:3),-(1:n.phe)],genotypes),
                    ncol=ncol(data)-n.phe,nrow=nrow(data)-3)

  # fix up map information
  # number of chromosomes
  uchr <- unique(chr)
  n.chr <- length(uchr)
  geno <- vector("list",n.chr)
  names(geno) <- uchr
  min.mar <- 1
  for(i in 1:n.chr) { # loop over chromosomes
    # create map
    temp.map <- map[chr==uchr[i]]
    names(temp.map) <- mnames[chr==uchr[i]]

    # pull out appropriate portion of genotype data
    data <- allgeno[,min.mar:(length(temp.map)+min.mar-1)]
    min.mar <- min.mar + length(temp.map)
    colnames(data) <- names(temp.map)

    geno[[i]] <- list(data=data,map=temp.map)
    if(uchr[i] == "X" || uchr[i] == "x") 
      class(geno[[i]]) <- "X"
    else class(geno[[i]]) <- "A"
  }

  # check that data dimensions match 
  n.mar1 <- sapply(geno,function(a) ncol(a$data))
  n.mar2 <- sapply(geno,function(a) length(a$map))
  n.phe <- ncol(pheno)
  n.ind1 <- nrow(pheno)
  n.ind2 <- sapply(geno,function(a) nrow(a$data))
  if(any(n.ind1 != n.ind2)) {
    print(c(n.ind1,n.ind2))
    stop("Number of individuals in genotypes and phenotypes do not match.");
  }
  if(any(n.mar1 != n.mar2)) {
    print(c(n.mar,n.mar2))
    stop("Numbers of markers in genotypes and marker names files do not match.");
  }

  # print some information about the amount of data read
  cat(" --Read the following data:\n");
  cat("\t", n.ind1, " individuals\n");
  cat("\t", sum(n.mar1), " markers\n");
  cat("\t", n.phe, " phenotypes\n");

  # determine map type: f2 or bc or 4way?
  if(max(allgeno[!is.na(allgeno)])<=2) type <- "bc"
  else type <- "f2"
  cross <- list(geno=geno,pheno=pheno)
  class(cross) <- c(type,"cross")

  # check that nothing is strange in the genotype data
  cross.type <- class(cross)[1]
  if(cross.type=="f2") max.gen <- 5
  else max.gen <- 2

  u <- unique(allgeno)
  if(any(!is.na(u) & (u > max.gen | u < 1))) 
    stop(paste("There are stange values in the genotype data :",
               paste(sort(u),collapse=":"), "."))

  cross
}
  

# end of read.cross.R
######################################################################
#
# ripple.R
#
# copyright (c) 2001, Karl W Broman, Johns Hopkins University
# Oct, 2001
# Licensed under the GNU General Public License version 2 (June, 1991)
# 
# Part of the R/qtl package
# Contains: ripple, summary.ripple
#           ripple.perm1, ripple.perm2, ripple.perm.sub
#
######################################################################

######################################################################
#
# ripple: Check marker orders for a given chromosome, comparing all
#         possible permutations of a sliding window of markers
#
######################################################################

ripple <-
function(cross, chr, window=4, error.prob=0,
         map.function=c("haldane","kosambi","c-f"),
         maxit=1000,tol=1e-5,sex.sp=TRUE)
{
  # pull out relevant chromosome
  if(length(chr) > 1)
    stop("ripple only works for one chromosome at a time.")
  cross <- pull.chr(cross,chr)
  chr.name <- names(cross$geno)[1]

  map.function <- match.arg(map.function)

  # get marker orders to test
  n.mar <- totmar(cross)
  if(n.mar <= window) # look at all possible orders
    orders <- ripple.perm2(n.mar)
  else { 
    temp <- ripple.perm1(window)
    n <- nrow(temp)
    orders <- cbind(temp,matrix(rep((window+1):n.mar,n),
                                  byrow=TRUE,ncol=n.mar-window))
    for(i in 2:(n.mar-window+1)) {
      left <- matrix(rep(1:(i-1),n),byrow=TRUE,ncol=i-1)
      if(i < n.mar-window+1)
        right <- matrix(rep((i+window):n.mar,n),byrow=TRUE,ncol=n.mar-window-i+1)
      else
        right <- NULL
      orders <- rbind(orders,cbind(left,temp+i-1,right))
    }
    # keep only distinct orders
    orders <- as.numeric(unlist(strsplit(unique(apply(orders,1,paste,collapse=":")),":")))
    orders <- matrix(orders,ncol=n.mar,byrow=TRUE)
  }

  m <- seq(0,by=5,length=n.mar)
  temcross <- cross
  if(is.matrix(cross$geno[[1]]$map)) 
    temcross$geno[[1]]$map <- rbind(m,m)
  else temcross$geno[[1]]$map <- m

  # calculate log likelihoods (and est'd chr length) for each marker order
  n.orders <- nrow(orders)
  loglik <- 1:n.orders
  chrlen <- 1:n.orders

  # how often to print information about current order being considered
  if(n.orders > 49) print.by <- 10
  else if(n.orders > 14) print.by <- 5
  else print.by <- 2

  for(i in 1:n.orders) {
    if(i==1) cat("  ", n.orders,"total orders\n")
    if((i %/% print.by)*print.by == i) cat("    --Order", i, "\n")
    temcross$geno[[1]]$data <- cross$geno[[1]]$data[,orders[i,]]
    newmap <- est.map(temcross,error.prob,map.function,maxit,tol,sex.sp)
    loglik[i] <- attr(newmap[[1]],"loglik")
    chrlen[i] <- newmap[[1]][n.mar]
  }

  # re-scale log likelihoods and convert to lods
  loglik <- (loglik - loglik[1])/log(10)
  
  # sort orders by lod
  o <- rev(order(loglik[-1])+1)

  orders <- cbind(orders,lod=loglik,chrlen)[c(1,o),]
  class(orders) <- c("ripple","matrix")
  attr(orders,"chr") <- chr.name
  attr(orders,"window") <- window
  attr(orders,"error.prob") <- error.prob
  orders
}

######################################################################
#
# summary.ripple: print top results from ripple().  We do this so
#                 that we can return *all* results but allow easy
#                 view of only the important ones
#
######################################################################

summary.ripple <-
function(object,lod.cutoff=2,...)
{
  n <- ncol(object)
  if(!any(object[-1,n-1] > -lod.cutoff)) object <- object[1:2,]
  else # make sure first row is included
    object <- rbind(object[1,],object[-1,][object[-1,n-1] > -lod.cutoff,])

  class(object) <- c("summary.ripple","matrix")
  object
}

######################################################################
#
# ripple.perm1: Utility function for ripple().  Returns all possible
#               permutations of {1, 2, ..., n}
#
######################################################################

ripple.perm1 <-  
function(n)
{
  if(n == 1) return(rbind(1))
  o <- rbind(c(n-1,n),c(n,n-1))
  if(n > 2)
    for(i in (n-2):1)
      o <- ripple.perm.sub(i,o)
  dimnames(o) <- NULL
  o
}

######################################################################
#
# ripple.perm2: Utility function for ripple().  Returns all possible
#               permutations of {1, 2, ..., n}, up to orientation of
#               the entire group
#
######################################################################

ripple.perm2 <- 
function(n)
{
  if(n < 3) return(rbind(1:n))
  o <- rbind(c(n-2,n-1,n),c(n-1,n-2,n),c(n-1,n,n-2))
  if(n > 3)
    for(i in (n-3):1)
      o <- ripple.perm.sub(i,o)
  dimnames(o) <- NULL
  o
}

######################################################################
#
# ripple.perm.sub: Subroutine used for ripple().  I'm too tired to
#                  explain.
#
######################################################################

ripple.perm.sub <-
function(x,mat)
{
  res <- cbind(x,mat)
  if(ncol(mat) > 1) {
    for(i in 1:ncol(mat))
        res <- rbind(res,cbind(mat[,1:i],x,mat[,-(1:i)]))
  }
  res
}

# end of ripple.R
######################################################################
#
# scanone.R
#
# copyright (c) 2001, Karl W Broman, Johns Hopkins University
# Oct, 2001; Sept, 2001; May, 2001, Apr, 2001; Feb, 2001
# Licensed under the GNU General Public License version 2 (June, 1991)
# 
# Part of the R/qtl package
# Contains: scanone, plot.scanone, scanone.perm
#           summary.scanone, print.summary.scanone
#
######################################################################

######################################################################
#
# scanone: scan genome, calculating LOD scores with single QTL model
#          (currently covariates not allowed; will be added later)
#
######################################################################

scanone <-
function(cross, chr, pheno.col=1, method=c("im","anova","hk","imp"),
         start=NULL, maxit=1000, tol=1e-8)
{
  method <- match.arg(method)
  if(method=="imp")
    warning("We don't have the imputation method working yet.")

  if(!missing(chr)) cross <- pull.chr(cross,chr)

  # check phenotypes
  if(length(pheno.col) > 1) pheno.col <- pheno.col[1]
  if(pheno.col < 1 || pheno.col > nphe(cross))
    stop("Specified phenotype column is invalid.")

  pheno <- cross$pheno[,pheno.col]
  keep.ind <- (1:length(pheno))[!is.na(pheno)]
  pheno <- pheno[keep.ind]
  n.ind <- length(keep.ind)
  n.chr <- nchr(cross)
  type <- class(cross)[1]

  if(method == "im") {
    m <- mean(pheno)
    s <- sd(pheno)*sqrt((n.ind-1)/n.ind)
    llik0 <- sum(log10(dnorm(pheno,m,s)))
  }
  else if(method=="hk")
    lrss0 <- log10(sum((pheno-mean(pheno))^2))

  results <- NULL

  if(method=="im") {
    if(is.null(start)) std.start <- 1
    else if(length(start)==1) std.start <- -1
    else std.start <- 0
  }

  # calculate genotype probabilities one chromosome at a time
  for(i in 1:n.chr) {

    # which type of cross is this?
    if(type == "f2") {
      if(class(cross$geno[[i]]) == "A") { # autosomal
        n.gen <- 3
        gen.names <- c("A","H","B")
      }
      else {                             # X chromsome 
        n.gen <- 2
        gen.names <- c("A","H","B") 
      }
    }
    else if(type == "bc") {
      n.gen <- 2
      gen.names <- c("A","H")
    }
    else if(type == "4way") {
      n.gen <- 4
      gen.names <- c("AC","AD","BC","BD")
    }
    else stop(paste("scanone not available for cross type",
                    type, "."))

    # starting values for interval mapping
    if(method=="im") {
      this.start <- rep(0,n.gen+1)
      if(std.start == 0) {
        if(length(start) < n.gen+1)
          stop(paste("Length of start argument should be 0, 1 or", n.gen+1))
        this.start <- c(start[1:n.gen],start[length(start)])
      }
    }

    # pull out reconstructed genotypes (anova)
    # or genotype probabilities (im or hk)

    if(method == "anova") {
      cfunc <- "R_scanone_anova"
      newgeno <- cross$geno[[i]]$data
      newgeno <- newgeno[keep.ind,]
      newgeno[is.na(newgeno)] <- 0 

      # discard partially informative genotypes
      if(type=="f2" || type=="f2ss") newgeno[newgeno>3] <- 0
      if(type=="4way") newgeno[newgeno>4] <- 0

      n.pos <- ncol(newgeno)
      map <- cross$geno[[i]]$map
      if(is.matrix(map)) map <- map[1,]
    }
    else if(method == "imp") {
      if(is.na(match("draws",names(cross$geno[[i]])))) { # need to run sim.geno
        warning("First running sim.geno.")
        cross <- sim.geno(cross)
      }

      draws <- cross$geno[[i]]$draws
      n.pos <- ncol(draws)
      n.draws <- dim(draws)[3]
      draws <- draws[keep.ind,,]

      map <- create.map(cross$geno[[i]]$map,
                        attr(cross$geno[[i]]$draws,"step"),
                        attr(cross$geno[[i]]$draws,"off.end"))
      if(is.matrix(map)) map <- map[1,]

      cfunc <- "R_scanone_imp"
    }
    else {
      if(is.na(match("prob",names(cross$geno[[i]])))) { # need to run calc.genoprob
        warning("First running calc.genoprob.")
        cross <- calc.genoprob(cross)
      }
      genoprob <- cross$geno[[i]]$prob
      n.pos <- ncol(genoprob)
      genoprob <- genoprob[keep.ind,,]

      map <- create.map(cross$geno[[i]]$map,
                        attr(cross$geno[[i]]$prob,"step"),
                        attr(cross$geno[[i]]$prob,"off.end"))
      if(is.matrix(map)) map <- map[1,]

      if(method == "im") cfunc <- "R_scanone_im"
      else cfunc <- "R_scanone_hk"
    }

    # call the C function
    if(method == "anova") {
      z <- .C(cfunc,
              as.integer(n.ind),         # number of individuals
              as.integer(n.pos),         # number of markers
              as.integer(n.gen),         # number of possible genotypes
              as.integer(newgeno),       # genotype data
              as.double(pheno),          # phenotype data
              result=as.double(rep(0,n.pos*(n.gen+2))),
              PACKAGE="qtl")
    }
    else if(method=="imp") {
      z <- .C(cfunc,
              as.integer(n.ind),
              as.integer(n.pos),
              as.integer(n.gen),
              as.integer(n.draws),
              as.integer(draws),
              as.double(pheno),
              result=as.double(rep(0,n.pos)),
              PACKAGE="qtl")
    }
    else if(method=="hk") {
      z <- .C(cfunc,
              as.integer(n.ind),         # number of individuals
              as.integer(n.pos),         # number of markers
              as.integer(n.gen),         # number of possible genotypes
              as.double(genoprob),       # genotype probabilities
              as.double(pheno),          # phenotype data
              result=as.double(rep(0,n.pos*(n.gen+2))),
              PACKAGE="qtl")
    }
    else { # interval mapping
      z <- .C(cfunc,
              as.integer(n.ind),         # number of individuals
              as.integer(n.pos),         # number of markers
              as.integer(n.gen),         # number of possible genotypes
              as.double(genoprob),       # genotype probabilities
              as.double(pheno),          # phenotype data
              result=as.double(rep(0,n.pos*(n.gen+2))),
              as.integer(std.start),
              as.double(this.start),
              as.integer(maxit),
              as.double(tol),
              PACKAGE="qtl")
    }

    z <- matrix(z$result,nrow=n.pos)

    if(method == "im")
      z[,1] <- z[,1] - llik0
    else if(method == "hk")
      z[,1] <- (n.ind/2)*(lrss0-z[,1])
#    else if(method == "imp")
#      z[,1] <- z[,1] - (n.gen-1)/2*log10(n.ind)

    if(method != "imp") {
      if(type=="f2" && class(cross$geno[[i]])=="X") # add BB column
        z <- cbind(z[,1:3],rep(NA,n.pos),z[,4])

      colnames(z) <- c("lod",gen.names,"sigma")
    }
    else colnames(z) <- c("lod")
      
    w <- names(map)
    o <- grep("^loc\-*[0-9]+",w)
    if(length(o) > 0) 
      w[o] <- paste(w[o],names(cross$geno)[i],sep=".c")
    rownames(z) <- w
    
    z <- as.data.frame(z)
    z <- cbind(chr=rep(names(cross$geno)[i],length(map)), pos=map, z)
    rownames(z) <- w
    results <- rbind(results,z)
  }

  # replace any lod = NaN with 0
  results[is.na(results[,3]),3] <- 0

  class(results) <- c("scanone",class(results))
  results
}

  

######################################################################
#
# plot.scanone: plot output from scanone
#
######################################################################

plot.scanone <- 
function(x,output2,output3,chr,incl.markers=TRUE,ylim,
         lty=c(1,2,3),col="black",lwd=2,add=FALSE,gap=25,...)
{
  output <- x
  second <- third <- TRUE
  if(missing(output2) && missing(output3)) 
     second <- third <- FALSE
  if(missing(output3))
    third <- FALSE
  if(missing(output2))
    second <- FALSE

  if(length(lty)==1) lty <- rep(lty,3)
  if(length(lwd)==1) lwd <- rep(lwd,3)
  if(length(col)==1) col <- rep(col,3)

  # pull out desired chromosomes
  if(missing(chr))
    chr <- unique(as.character(output[,1]))

  if(length(chr) == 0) chr <- sort(unique(output[,1]))
  else if(all(chr < 0)) { 
    a <- sort(unique(output[,1]))
    chr <- a[-match(-chr,a)]
  }
  output <- output[!is.na(match(output[,1],chr)),]
  if(second) output2 <- output2[!is.na(match(output2[,1],chr)),]
  if(third) output3 <- output3[!is.na(match(output3[,1],chr)),]
  
  # beginning and end of chromosomes
  temp <- grep("^loc\-*[0-9]+",rownames(output))
  if(length(temp)==0) temp <- output
  else temp <- output[-temp,]
  begend <- matrix(unlist(tapply(temp[,2],temp[,1],range)),ncol=2,byrow=TRUE)
  len <- begend[,2]-begend[,1]

  # locations to plot start of each chromosome
  start <- gap/2+c(0,cumsum(len+gap))-c(begend[,1],0)

  maxx <- sum(len+gap)
  maxy <- max(output[,3])
  if(second) maxy <- max(c(maxy,output2[,3]))
  if(third) maxy <- max(c(maxy,output3[,3]))

  # graphics parameters
  old.xpd <- par("xpd")
  old.las <- par("las")
  par(xpd=TRUE,las=1)
  on.exit(par(xpd=old.xpd,las=old.las))

  # make frame of plot
  if(missing(ylim)) ylim <- c(0,maxy)

  if(!add)
    plot(0,0,ylim=ylim,xlim=c(0,maxx),type="n",
         xlab="Map position (cM)",ylab=dimnames(output)[[2]][3])

  for(i in 1:length(chr)) {
    # plot first output
    x <- output[output[,1]==chr[i],2]+start[i]
    y <- output[output[,1]==chr[i],3]
    lines(x,y,lwd=lwd[1],lty=lty[1],col=col[1])

    # plot chromosome number
    a <- par("usr")
    if(!add) {
      tloc <- mean(c(start[i],start[i+1]-gap))
      text(tloc,a[4]+(a[4]-a[3])*0.03,as.character(chr[i]))
      lines(rep(tloc,2),c(a[4],a[4]+(a[4]-a[3])*0.015))
    }

    # plot second output
    if(second) {
      x <- output2[output2[,1]==chr[i],2]+start[i]
      y <- output2[output2[,1]==chr[i],3]
      lines(x,y,lty=lty[2],col=col[2],lwd=lwd[2])
    }

    if(third) {
      x <- output3[output3[,1]==chr[i],2]+start[i]
      y <- output3[output3[,1]==chr[i],3]
      lines(x,y,lty=lty[3],col=col[3],lwd=lwd[3])
    }

    # plot lines at marker positions
    if(incl.markers && !add) {
      nam <- dimnames(output)[[1]][output[,1]==chr[i]]
      wh.genoprob <- (1:length(nam))[grep("^loc\-*[0-9]+",nam)]
      if(length(wh.genoprob)==0) wh.genoprob <- 1:length(nam)
      else wh.genoprob <- (1:length(nam))[-wh.genoprob]
      pos <- output[output[,1]==chr[i],2][wh.genoprob]+start[i]
      for(j in pos) 
	lines(c(j,j),c(a[3],a[3]+(a[4]-a[3])*0.02))
    }

  }

}



######################################################################
#
# scanone.perm: Permutation test of scanone
#
######################################################################

scanone.perm <-
function(cross, chr, pheno.col=1, method=c("im","anova","hk"),
         start=NULL,n.perm=1000, maxit=1000, tol=1e-8)
{
  method <- match.arg(method)

  if(method=="hk") {
    warning("We don't have Haley-Knott working yet; running IM instead.")
    method <- "im"
  }

  if(!missing(chr))
    cross <- pull.chr(cross,chr)

  # check phenotypes
  if(length(pheno.col) > 1) pheno.col <- pheno.col[1]
  if(pheno.col < 1 || pheno.col > nphe(cross))
    stop("Specified phenotype column is invalid.")

  if(method=="im") {
    if(is.null(start)) std.start <- 1
    else if(length(start)==1) std.start <- -1
    else std.start <- 0
  }

  pheno <- cross$pheno[,pheno.col]
  keep.ind <- (1:length(pheno))[!is.na(pheno)]
  pheno <- pheno[keep.ind]
  n.ind <- length(keep.ind)
  n.chr <- nchr(cross)
  type <- class(cross)[1]

  if(method == "im") {
    m <- mean(pheno)
    s <- sd(pheno)*sqrt((n.ind-1)/n.ind)
    llik0 <- sum(log10(dnorm(pheno,m,s)))
  }
  else lrss0 <- log10(sum((pheno-mean(pheno))^2))

  results <- NULL

  # calculate genotype probabilities one chromosome at a time
  for(i in 1:n.chr) {

    # which type of cross is this?
    if(type == "f2") {
      if(class(cross$geno[[i]]) == "A")  # autosomal
        n.gen <- 3
      else                              # X chromsome 
        n.gen <- 2
    }
    else if(type == "bc") 
      n.gen <- 2
    else if(type == "4way") 
      n.gen <- 4
    else stop(paste("scanone not available for cross type",
                    type, "."))

    if(method=="im") {
      # starting values
      this.start <- rep(0,n.gen+1)
      if(std.start == 0) {
        if(length(start) < n.gen+1)
          stop(paste("Length of start argument should be 0, 1 or", n.gen+1))
        this.start <- c(start[1:n.gen],start[length(start)])
      }
    }

    # pull out genotypes (anova) or genotype probabilities (im or hk)

    if(method == "anova") {
      cfunc <- "scanone_anova_perm"
      newgeno <- cross$geno[[i]]$data
      newgeno <- newgeno[keep.ind,]
      newgeno[is.na(newgeno)] <- 0 

      # discard partially informative genotypes
      if(type=="f2" || type=="f2ss") newgeno[newgeno>3] <- 0
      if(type=="4way") newgeno[newgeno>4] <- 0

      n.pos <- ncol(newgeno)
      map <- cross$geno[[i]]$map
      if(is.matrix(map)) map <- map[1,]
    }
    else {
      if(is.na(match("prob",names(cross$geno[[i]])))) { # need to run calc.genoprob
        warning("First running calc.genoprob.")
        cross <- calc.genoprob(cross)
      }
      genoprob <- cross$geno[[i]]$prob
      n.pos <- ncol(genoprob)
      genoprob <- genoprob[keep.ind,,]

      map <- create.map(cross$geno[[i]]$map,
                        attr(cross$geno[[i]]$prob,"step"),
                        attr(cross$geno[[i]]$prob,"off.end"))
      if(is.matrix(map)) map <- map[1,]

      if(method == "im") cfunc <- "scanone_im_perm"
      else cfunc <- "scanone_hk_perm"
    }

    # call the C function
    if(method=="anova") 
      z <- .C(cfunc,
              as.integer(n.ind),         # number of individuals
              as.integer(n.pos),         # number of markers
              as.integer(n.gen),         # number of possible genotypes
              as.integer(newgeno),       # genotype data
              as.double(pheno),          # phenotype data
              as.integer(n.perm),
              result=as.double(rep(0,n.perm)),
              PACKAGE="qtl")
    else if(method=="hk") 
      z <- .C(cfunc,
              as.integer(n.ind),         # number of individuals
              as.integer(n.pos),         # number of markers
              as.integer(n.gen),         # number of possible genotypes
              as.double(genoprob),       # genotype probabilities
              as.double(pheno),          # phenotype data
              as.integer(n.perm),
              result=as.double(rep(0,n.perm)),
              PACKAGE="qtl")
    else # interval mapping
      z <- .C(cfunc,
              as.integer(n.ind),         # number of individuals
              as.integer(n.pos),         # number of markers
              as.integer(n.gen),         # number of possible genotypes
              as.double(genoprob),       # genotype probabilities
              as.double(pheno),          # phenotype data
              as.integer(n.perm),
              result=as.double(rep(0,n.perm)),
              as.integer(std.start),
              as.double(this.start),
              as.integer(maxit),
              as.double(tol),
              PACKAGE="qtl")

    if(method == "im")
      z <- z$result - llik0
    else 
      z <- (n.ind/2)*(lrss0-z$result)


    results <- cbind(results,z)
  }

  results <- cbind(apply(results,1,max),results)
  colnames(results) <- c("max",names(cross$geno))

  results
}


# give, for each chromosome, the output line at the maximum LOD
summary.scanone <-
function(object,threshold=0,...)
{
  output <- lapply(split(object,object[,1]),
                   function(b) b[b[,3]==max(b[,3]),])
  results <- output[[1]]
  if(length(output) > 1)
    for(i in 2:length(output))
      results <- rbind(results,output[[i]])
  class(results) <- c("summary.scanone","data.frame")
  if(!any(results[,3] >= threshold)) {
    cat("    There were no LOD peaks above the threshold.\n")
    invisible()
  }
  else {
    return(results[results[,3] >= threshold,])
  }
}

# print output of summary.scanone
print.summary.scanone <-
function(x,...)
{
  x[,-(1:2)] <- round(data.frame(x[,-(1:2)]),6)
  print.data.frame(x,digits=2)
}


# end of scanone.R
######################################################################
#
# sim.geno.R
#
# copyright (c) 2001, Karl W Broman, Johns Hopkins University
# Sept, 2001; July, 2001; May, 2001; Apr, 2001; Feb, 2001
# Licensed under the GNU General Public License version 2 (June, 1991)
# 
# Part of the R/qtl package
# Contains: sim.geno
#
######################################################################

######################################################################
#
# sim.geno: simulate from the joint distribution Pr(g | O)
#
######################################################################

sim.geno <-
function(cross, n.draws=1, step=0, off.end=0, error.prob=0,
         map.function=c("haldane","kosambi","c-f"))
{

  # map function
  map.function <- match.arg(map.function)
  if(map.function=="kosambi") mf <- mf.k
  else if(map.function=="c-f") mf <- mf.cf
  else mf <- mf.h

  # don't let error.prob be exactly zero, just in case
  if(error.prob < 1e-14) error.prob <- 1e-14

  n.ind <- nind(cross)
  n.chr <- nchr(cross)

  # calculate genotype probabilities one chromosome at a time
  for(i in 1:n.chr) {

    # which type of cross is this?
    if(class(cross)[1] == "f2") {
      n.gen <- 3
      gen.names <- c("A","H","B")
      one.map <- TRUE
      if(class(cross$geno[[i]]) == "A") # autosomal
        cfunc <- "sim_geno_f2"
      else                              # X chromsome
        cfunc <- "sim_geno_bc"
    }
    else if(class(cross)[1] == "bc") {
      cfunc <- "sim_geno_bc"
      n.gen <- 2
      gen.names <- c("A","H")
      one.map <- TRUE
    }
    else if(class(cross)[1] == "4way") {
      n.gen <- 4
      gen.names <- c("AC","AD","BC","BD")
      cfunc <- "sim_geno_4way"
      one.map <- FALSE
    }
    else stop(paste("sim_geno not available for cross type",
                    class(cross)[1], "."))

    # genotype data
    gen <- cross$geno[[i]]$data
    gen[is.na(gen)] <- 0
    
    # recombination fractions
    if(one.map) {
      # recombination fractions
      map <- create.map(cross$geno[[i]]$map,step,off.end)
      rf <- mf(diff(map))
      rf[rf < 1e-14] <- 1e-14

      # new genotype matrix with pseudomarkers filled in
      newgen <- matrix(ncol=length(map),nrow=nrow(gen))
      dimnames(newgen) <- list(NULL,names(map))
      newgen[,colnames(gen)] <- gen
      newgen[is.na(newgen)] <- 0
      n.pos <- ncol(newgen)
    }
    else {
      map <- create.map(cross$geno[[i]]$map,step,off.end)
      rf <- mf(diff(map[1,]))
      rf[rf < 1e-14] <- 1e-14
      rf2 <- mf(diff(map[1,]))
      rf2[rf2 < 1e-14] <- 1e-14

      # new genotype matrix with pseudomarkers filled in
      newgen <- matrix(ncol=ncol(map),nrow=nrow(gen))
      dimnames(newgen) <- list(NULL,colnames(map))
      newgen[,colnames(gen)] <- gen
      newgen[is.na(newgen)] <- 0
      n.pos <- ncol(newgen)
    }

    
    # call C function
    if(one.map) {
      z <- .C(cfunc,
              as.integer(n.ind),         # number of individuals
              as.integer(n.pos),         # number of markers
              as.integer(n.draws),       # number of simulation replicates
              as.integer(newgen),        # genotype data
              as.double(rf),             # recombination fractions
              as.double(error.prob),     # 
              draws=as.integer(rep(0,n.draws*n.ind*n.pos)), 
              PACKAGE="qtl")

      cross$geno[[i]]$draws <- array(z$draws,dim=c(n.ind,n.pos,n.draws))
      dimnames(cross$geno[[i]]$draws) <- list(NULL, names(map), NULL)
    }
    else {
      z <- .C(cfunc,
              as.integer(n.ind),         # number of individuals
              as.integer(n.pos),         # number of markers
              as.integer(n.draws),       # number of simulation replicates
              as.integer(newgen),        # genotype data
              as.double(rf),             # recombination fractions
              as.double(rf2),            # recombination fractions
              as.double(error.prob),     # 
              draws=as.integer(rep(0,n.draws*n.ind*n.pos)),
              PACKAGE="qtl")

      cross$geno[[i]]$draws <- array(z$draws,dim=c(n.ind,n.pos,n.draws))
      dimnames(cross$geno[[i]]$draws) <- list(NULL, colnames(map), NULL)

    }

    # attribute set to the error.prob value used, for later
    #     reference
    attr(cross$geno[[i]]$draws,"error.prob") <- error.prob
    attr(cross$geno[[i]]$draws,"step") <- step
    attr(cross$geno[[i]]$draws,"off.end") <- off.end
  }

  cross
}

  

# end of sim.geno.R
######################################################################
#
# simulate.R
#
# copyright (c) 2001, Karl W Broman, Johns Hopkins University
# Oct, 2001; July, 2001; May, 2001; April, 2001
# Licensed under the GNU General Public License version 2 (June, 1991)
# 
# Part of the R/qtl package
# Contains: sim.map, sim.cross, sim.cross.bc, sim.cross.f2,
#           sim.cross.4way
#
######################################################################

######################################################################
#
# sim.map: simulate a genetic map
#
######################################################################

sim.map <-
function(len=rep(100,20), n.mar=10, anchor.tel=TRUE, include.x=TRUE,
         sex.sp=FALSE)
{
  if(length(len)!=length(n.mar) && length(len)!=1 && length(n.mar)!=1)
    stop("Lengths of vectors len and n.mar do not conform.")

  # make vectors the same length
  if(length(len) == 1) len <- rep(len,length(n.mar))
  else if(length(n.mar) == 1) n.mar <- rep(n.mar,length(len))

  map <- vector("list",length(len))
  names(map) <- as.character(1:length(map))
  if(include.x) names(map)[length(map)] <- "X"

  for(i in 1:length(len)) {
    if(anchor.tel) {
      if(n.mar[i] < 2) n.mar[i] <- 2
      map[[i]] <- c(0,len[i])
      if(n.mar[i] > 2) map[[i]] <- sort(c(map[[i]],runif(n.mar[i]-2,0,len[i])))
    }
    else {
      map[[i]] <- sort(runif(n.mar[i],0,len[i]))
      map[[i]] <- map[[i]] - min(map[[i]])
    }
    names(map[[i]]) <- paste("D", names(map)[i], "M", 1:n.mar[i], sep="")
    class(map[[i]]) <- "A"
  }

  if(sex.sp) {
    for(i in 1:length(len)) {
      if(anchor.tel) {
        if(n.mar[i] < 2) n.mar[i] <- 2
        tempmap <- c(0,len[i])
        if(n.mar[i] > 2) tempmap <- sort(c(tempmap,runif(n.mar[i]-2,0,len[i])))
      }
      else {
        tempmap <- sort(runif(n.mar[i],0,len[i]))
        tempmap <- tempmap - min(tempmap)
      }
      map[[i]] <- rbind(map[[i]],tempmap)
      dimnames(map[[i]]) <- list(NULL,paste("D", names(map)[i], "M", 1:n.mar[i], sep=""))
      class(map[[i]]) <- "A"

      if(include.x && i==length(len))  # if X chromosome, force no recombination in male
        map[[i]][2,] <- rep(0,ncol(map[[i]]))

    }
  }

  if(include.x) class(map[[length(map)]]) <- "X"

  class(map) <- "map"
  map
}

    
    
  
######################################################################
#
# sim.cross: Simulate an experimental cross
#
# Note: These functions are a bit of a mess.  I was in the "get it to
#       work without worrying about efficiency" mode while writing it.
#       Sorry!
#
######################################################################

sim.cross <-
function(map, model=NULL, n.ind=100, type=c("f2","bc","4way"),
         error.prob=0, missing.prob=0, partial.missing.prob=0,
         keep.qtlgeno=FALSE, error.ind=TRUE,
         map.function=c("haldane","kosambi","c-f"))
{
  type <- match.arg(type)
  map.function <- match.arg(map.function)

  if(type=="bc")
    cross <- sim.cross.bc(map,model,n.ind,error.prob,missing.prob,
                          error.ind,map.function)
  else if(type=="f2")
    cross <- sim.cross.f2(map,model,n.ind,error.prob,missing.prob,
                          partial.missing.prob,error.ind,map.function)
  else
    cross <- sim.cross.4way(map,model,n.ind,error.prob,missing.prob,
                            partial.missing.prob,error.ind,map.function)
  if(!keep.qtlgeno)
    cross <- drop.qtlgeno(cross)

  cross

}


######################################################################
#
# sim.cross.bc
#
######################################################################

sim.cross.bc <-
function(map,model,n.ind,error.prob,missing.prob,
         error.ind,map.function)
{
  if(map.function=="kosambi") mf <- mf.k
  else if(map.function=="c-f") mf <- mf.cf
  else mf <- mf.h

  if(any(sapply(map,is.matrix)))
    stop("Map must not be sex-specific.")

  n.chr <- length(map)

  if(is.null(model)) n.qtl <- 0
  else {
    if(!((!is.matrix(model) && length(model) == 3) ||
         (is.matrix(model) && ncol(model) == 3))) 
      stop("Model must be a matrix with 3 columns (chr, pos and effect).")
    if(!is.matrix(model)) model <- rbind(model)
    n.qtl <- nrow(model)
    if(any(model[,1] < 0 | model[,1] > length(map)))
      stop("Chromosome indicators in model matrix out of range.")
    model[,2] <- model[,2]+1e-14 # so QTL not on top of marker
  }

  # if any QTLs, place qtls on map
  if(n.qtl > 0) {
    for(i in 1:n.qtl) {
      temp <- map[[model[i,1]]]
      if(model[i,2] < min(temp)) {
        temp <- c(model[i,2],temp)
        names(temp)[1] <- paste("QTL",i,sep="")
      }
      else if(model[i,2] > max(temp)) {
        temp <- c(temp,model[i,2])
        names(temp)[length(temp)] <- paste("QTL",i,sep="")
      }
      else {
        j <- max((1:length(temp))[temp < model[i,2]])
        temp <- c(temp[1:j],model[i,2],temp[(j+1):length(temp)])
        names(temp)[j+1] <- paste("QTL",i,sep="")
      }
      map[[model[i,1]]] <- temp
    }
  }
  
  geno <- vector("list", n.chr)
  names(geno) <- names(map)
  n.mar <- sapply(map,length)
  mar.names <- lapply(map,names)
  chr.type <- sapply(map,function(a)
                     if(is.null(class(a))) return("A")
                     else return(class(a)))
  
  for(i in 1:n.chr) {
    data <- matrix(nrow=n.ind,ncol=n.mar[i])
    dimnames(data) <- list(NULL,mar.names[[i]])

    # simulate genotype data
    d <- diff(map[[i]]) # inter-marker distances (cM)
    r <- mf(d) # recombination fractions (Kosambi map function)
    rbar <- 1-r

    # first locus on chromosome
    data[,1] <- sample(1:2,n.ind,repl=TRUE)

    # rest of markers
    if(n.mar[i] > 1) {
      for(j in 1:(n.mar[i]-1)) {
        rec <- sample(0:1,n.ind,repl=TRUE,prob=c(1-r[j],r[j]))
        data[rec==0,j+1] <- data[rec==0,j]
        data[rec==1,j+1] <- 3-data[rec==1,j]
      }
    } # if n.mar[i] > 1

    geno[[i]] <- list(data = data, map = map[[i]])
    class(geno[[i]]) <- chr.type[i]
    class(geno[[i]]$map) <- NULL
    
  } # end loop over chromosomes

  # simulate phenotypes
  pheno <- rnorm(n.ind,0,1)

  if(n.qtl > 0) {
    # find QTL positions in genotype data
    QTL.chr <- QTL.loc <- NULL
    for(i in 1:n.chr) {
      o <- grep("^QTL[0-9]+",mar.names[[i]])
      if(length(o)>0) {
        QTL.chr <- c(QTL.chr,rep(i,length(o)))
        QTL.loc <- c(QTL.loc,o)
      }
    }

    # incorporate QTL effects
    for(i in 1:n.qtl) {
      QTL.geno <- geno[[QTL.chr[i]]]$data[,QTL.loc[i]]
      pheno[QTL.geno==1] <- pheno[QTL.geno==1] - model[i,3]
      pheno[QTL.geno==2] <- pheno[QTL.geno==2] + model[i,3]
    }

  } # end simulate phenotype
      
  n.mar <- sapply(geno, function(a) length(a$map))

  # add errors
  if(error.prob > 0) {
    for(i in 1:n.chr) {
      a <- sample(0:1,n.mar[i]*n.ind,repl=TRUE,
                  prob=c(1-error.prob,error.prob))
      geno[[i]]$data[a == 1] <- 3 - geno[[i]]$data[a == 1]
      if(error.ind) {
        errors <- matrix(0,n.ind,n.mar[i])
        errors[a==1] <- 1
        colnames(errors) <- colnames(geno[[i]]$data)
        geno[[i]]$errors <- errors
      }
    } 
  } 

  # add missing
  if(missing.prob > 0) {
    for(i in 1:n.chr) {
      o <- grep("^QTL[0-9]+",mar.names[[i]])
      if(length(o)>0)
        x <- geno[[i]]$data[,o]
      geno[[i]]$data[sample(c(TRUE,FALSE),n.mar[i]*n.ind,repl=TRUE,
                            prob=c(missing.prob,1-missing.prob))] <- NA
      if(length(o)>0)
        geno[[i]]$data[,o] <- x
    }
  }

  pheno <- as.matrix(pheno)
  dimnames(pheno) <- list(NULL, "phenotype")

  cross <- list(geno=geno,pheno=pheno)
  class(cross) <- c("bc","cross")

  cross
}  
       
######################################################################
#
# sim.cross.f2
#
######################################################################

sim.cross.f2 <-              
function(map,model,n.ind,error.prob,missing.prob,partial.missing.prob,
         error.ind,map.function)
{
  if(map.function=="kosambi") mf <- mf.k
  else if(map.function=="c-f") mf <- mf.cf
  else mf <- mf.h

  if(any(sapply(map,is.matrix)))
    stop("Map must not be sex-specific.")

  n.chr <- length(map)
  if(is.null(model)) n.qtl <- 0
  else {
    if(!((!is.matrix(model) && length(model) == 4) ||
         (is.matrix(model) && ncol(model) == 4))) {
      stop("Model must be a matrix with 4 columns (chr, pos and effects).")
    }
    if(!is.matrix(model)) model <- rbind(model)
    n.qtl <- nrow(model)
    if(any(model[,1] < 0 | model[,1] > length(map)))
      stop("Chromosome indicators in model matrix out of range.")
    model[,2] <- model[,2]+1e-14 # so QTL not on top of marker
  }

  # if any QTLs, place qtls on map
  if(n.qtl > 0) {
    for(i in 1:n.qtl) {
      temp <- map[[model[i,1]]]
      if(model[i,2] < min(temp)) {
        temp <- c(model[i,2],temp)
        names(temp)[1] <- paste("QTL",i,sep="")
      }
      else if(model[i,2] > max(temp)) {
        temp <- c(temp,model[i,2])
        names(temp)[length(temp)] <- paste("QTL",i,sep="")
      }
      else {
        j <- max((1:length(temp))[temp < model[i,2]])
        temp <- c(temp[1:j],model[i,2],temp[(j+1):length(temp)])
        names(temp)[j+1] <- paste("QTL",i,sep="")
      }
      map[[model[i,1]]] <- temp
    }
  }
  
  geno <- vector("list", n.chr)
  names(geno) <- names(map)
  n.mar <- sapply(map,length)
  mar.names <- lapply(map,names)
  chr.type <- sapply(map,function(a)
                     if(is.null(class(a))) return("A")
                     else return(class(a)))
  
  for(i in 1:n.chr) {

    data <- matrix(nrow=n.ind,ncol=n.mar[i])
    dimnames(data) <- list(NULL,mar.names[[i]])

    # simulate genotype data
    d <- diff(map[[i]]) # inter-marker distances (cM)
    r <- mf(d) # recombination fractions (Kosambi map function)
    rbar <- 1-r

    # first locus on chromosome
    if(chr.type[i]=="X") data[,1] <- sample(1:2,n.ind,repl=TRUE)
    else data[,1] <- sample(1:3,n.ind,repl=TRUE,prob=c(1,2,1))
    
    # rest of markers
    if(n.mar[i] > 1) {
      for(j in 1:(n.mar[i]-1)) {
        if(chr.type[i]=="X") { # X chromosome (like a backcross)
          rec <- sample(0:1,n.ind,repl=TRUE,prob=c(1-r[j],r[j]))
          data[rec==0,j+1] <- data[rec==0,j]
          data[rec==1,j+1] <- 3-data[rec==1,j]
        }
        else { # F2 autosome
          data[data[,j]==1,j+1] <- sample(1:3,sum(data[,j]==1),repl=TRUE,
                     prob=c(rbar[j]*rbar[j],2*r[j]*rbar[j],r[j]*r[j]))
          data[data[,j]==2,j+1] <- sample(1:3,sum(data[,j]==2),repl=TRUE,
                     prob=c(r[j]*rbar[j],rbar[j]*rbar[j]+r[j]*r[j],r[j]*rbar[j]))
          data[data[,j]==3,j+1] <- sample(1:3,sum(data[,j]==3),repl=TRUE,
                     prob=c(r[j]*r[j],2*r[j]*rbar[j],rbar[j]*rbar[j]))
        }
      } # end loop over intervals
    } # if n.mar[i] > 1

    geno[[i]] <- list(data = data, map = map[[i]])
    class(geno[[i]]) <- chr.type[i]
    class(geno[[i]]$map) <- NULL
    
  } # end loop over chromosomes

  # simulate phenotypes
  pheno <- rnorm(n.ind,0,1)

  if(n.qtl > 0) {
    # find QTL positions in genotype data
    QTL.chr <- QTL.loc <- NULL
    for(i in 1:n.chr) {
      o <- grep("^QTL[0-9]+",mar.names[[i]])
      if(length(o)>0) {
        QTL.chr <- c(QTL.chr,rep(i,length(o)))
        QTL.loc <- c(QTL.loc,o)
      }
    }

    # incorporate QTL effects
    for(i in 1:n.qtl) {
      QTL.geno <- geno[[QTL.chr[i]]]$data[,QTL.loc[i]]
      pheno[QTL.geno==1] <- pheno[QTL.geno==1] - model[i,3]
      pheno[QTL.geno==2] <- pheno[QTL.geno==2] + model[i,4]
      pheno[QTL.geno==3] <- pheno[QTL.geno==3] + model[i,3]
    }

  } # end simulate phenotype
      
  n.mar <- sapply(geno, function(a) length(a$map))

  # add errors
  if(error.prob > 0) {
    for(i in 1:n.chr) {
      if(chr.type[i]=="X") {
        a <- sample(0:1,n.mar[i]*n.ind,repl=TRUE,
                    prob=c(1-error.prob,error.prob))
        geno[[i]]$data[a == 1] <- 3 - geno[[i]]$data[a == 1]
      }
      else {
        a <- sample(0:2,n.mar[i]*n.ind,repl=TRUE,
                    prob=c(1-error.prob,error.prob/2,error.prob/2))
        if(any(a>0 & geno[[i]]$data==1))
          geno[[i]]$data[a>0 & geno[[i]]$data==1] <-
            (geno[[i]]$data+a)[a>0 & geno[[i]]$data==1]
        if(any(a>0 & geno[[i]]$data==2)) {
          geno[[i]]$data[a>0 & geno[[i]]$data==2] <-
            (geno[[i]]$data+a)[a>0 & geno[[i]]$data==2]
          geno[[i]]$data[geno[[i]]$data>3] <- 1
        }
        if(any(a>0 & geno[[i]]$data==3))
          geno[[i]]$data[a>0 & geno[[i]]$data==3] <-
            (geno[[i]]$data-a)[a>0 & geno[[i]]$data==3]
      }

      if(error.ind) {
        errors <- matrix(0,n.ind,n.mar[i])
        errors[a>0] <- 1
        colnames(errors) <- colnames(geno[[i]]$data)
        geno[[i]]$errors <- errors
      }

    } # end loop over chromosomes
  } # end simulate genotyping errors

  # add partial missing
  if(partial.missing.prob > 0) {
    for(i in 1:n.chr) {
      if(chr.type[i] != "X") {
        o <- sample(c(TRUE,FALSE),n.mar[i],repl=TRUE,
                    prob=c(partial.missing.prob,1-partial.missing.prob))
        if(any(o)) {
          o2 <- grep("^QTL[0-9]+",mar.names[[i]])
          if(length(o2)>0)
            x <- geno[[i]]$data[,o2]
          m <- (1:n.mar[i])[o]
          for(j in m) {
            if(runif(1) < 0.5) 
              geno[[i]]$data[geno[[i]]$data[,j]==1 | geno[[i]]$data[,j]==2,j] <- 4
            else 
              geno[[i]]$data[geno[[i]]$data[,j]==3 | geno[[i]]$data[,j]==2,j] <- 5
          }
          if(length(o2)>0)
            geno[[i]]$data[,o2] <- x
        }
      }

    } # end loop over chromosomes
  } # end simulate partially missing data
            
  # add missing
  if(missing.prob > 0) {
    for(i in 1:n.chr) {
      o <- grep("^QTL[0-9]+",mar.names[[i]])
      if(length(o)>0)
        x <- geno[[i]]$data[,o]
      geno[[i]]$data[sample(c(TRUE,FALSE),n.mar[i]*n.ind,repl=TRUE,
                            prob=c(missing.prob,1-missing.prob))] <- NA
      if(length(o)>0)
        geno[[i]]$data[,o] <- x
    }
  }

  pheno <- as.matrix(pheno)
  dimnames(pheno) <- list(NULL, "phenotype")

  cross <- list(geno=geno,pheno=pheno)
  class(cross) <- c("f2","cross")

  cross
}

######################################################################
#
# sim.cross.4way
#
######################################################################

sim.cross.4way <-              
function(map,model,n.ind,error.prob,missing.prob,partial.missing.prob,
         error.ind,map.function)
{
  if(map.function=="kosambi") mf <- mf.k
  else if(map.function=="c-f") mf <- mf.cf
  else mf <- mf.h

  if(!all(sapply(map,is.matrix)))
    stop("Map must be sex-specific.")

  n.chr <- length(map)
  if(is.null(model)) n.qtl <- 0
  else {
    if(!((!is.matrix(model) && length(model) == 5) ||
         (is.matrix(model) && ncol(model) == 5))) {
      stop("Model must be a matrix with 5 columns (chr, pos and effects).")
    }
    if(!is.matrix(model)) model <- rbind(model)
    n.qtl <- nrow(model)
    if(any(model[,1] < 0 | model[,1] > length(map)))
      stop("Chromosome indicators in model matrix out of range.")
    model[,2] <- model[,2]+1e-14 # so QTL not on top of marker
  }

  # if any QTLs, place qtls on map
  if(n.qtl > 0) {
    for(i in 1:n.qtl) {
      temp <- map[[model[i,1]]]
      temp1 <- temp[1,]
      temp2 <- temp[2,]
      qtlloc <- model[i,2]

      if(qtlloc < min(temp1)) {
        temp1 <- c(qtlloc,temp1)
        temp2 <- min(temp2) - (min(temp1)-qtlloc)/diff(range(temp1))*diff(range(temp2))
        temp1 <- temp1-min(temp1)
        temp2 <- temp2-min(temp2)
        n <- c(paste("QTL",i,sep=""),colnames(temp))
      }
      else if(qtlloc > max(temp1)) {
        temp1 <- c(temp1,qtlloc)
        temp2 <- (qtlloc-max(temp1))/diff(range(temp1))*diff(range(temp2))+max(temp2)
        n <- c(colnames(temp),paste("QTL",i,sep=""))
      }
      else {
        temp1 <- c(temp1,qtlloc)
        o <- order(temp1)
        wh <- (1:length(temp1))[order(temp1)==length(temp1)]
        temp2 <- c(temp2[1:(wh-1)],NA,temp2[-(1:(wh-1))])
        temp2[wh] <- temp2[wh-1] + (temp1[wh]-temp1[wh-1])/(temp1[wh+1]-temp1[wh-1]) *
          (temp2[wh+1]-temp2[wh-1])
        temp1 <- sort(temp1)
        n <- c(colnames(temp),paste("QTL",i,sep=""))[o]
      }
      map[[model[i,1]]] <- rbind(temp1,temp2)
      dimnames(map[[model[i,1]]]) <- list(NULL, n)
    }
  }
  
  geno <- vector("list", n.chr)
  names(geno) <- names(map)
  n.mar <- sapply(map,ncol)
  mar.names <- lapply(map,function(a) colnames(a))
  chr.type <- sapply(map,function(a)
                     if(is.null(class(a))) return("A")
                     else return(class(a)))
  
  for(i in 1:n.chr) {

    data <- matrix(nrow=n.ind,ncol=n.mar[i])
    dimnames(data) <- list(NULL,mar.names[[i]])

    # simulate genotype data
    d <- diff(map[[i]][1,])
    r <- mf(d)
    rbar <- 1-r
    d2 <- diff(map[[i]][2,])
    r2 <- mf(d2)
    rbar2 <- 1-r2
    
    data2 <- data

    # first locus on chromosome
    data[,1] <- sample(1:2,n.ind,repl=TRUE)  # mother's chromosome
    data2[,1] <- sample(1:2,n.ind,repl=TRUE) # father's chromosome

    sex <- NULL
    if(chr.type[i]=="X") {
      sex <- rep(0,n.ind)
      sex[data2[,1]==2] <- 1
    }

    # rest of markers
    if(n.mar[i] > 1) {
      for(j in 1:(n.mar[i]-1)) {
        rec <- sample(0:1,n.ind,repl=TRUE,prob=c(1-r[j],r[j]))
        data[rec==0,j+1] <- data[rec==0,j]
        data[rec==1,j+1] <- 3-data[rec==1,j]
          
        rec <- sample(0:1,n.ind,repl=TRUE,prob=c(1-r2[j],r2[j]))
        data2[rec==0,j+1] <- data2[rec==0,j]
        data2[rec==1,j+1] <- 3-data2[rec==1,j]
      }
    } 
    data <- data + (data2-1)*2 

    geno[[i]] <- list(data = data, map = map[[i]])
    class(geno[[i]]) <- chr.type[i]
    class(geno[[i]]$map) <- NULL
    
  } # end loop over chromosomes

  # simulate phenotypes
  pheno <- rnorm(n.ind,0,1)

  if(n.qtl > 0) {
    # find QTL positions
    QTL.chr <- QTL.loc <- NULL
    for(i in 1:n.chr) {
      o <- grep("^QTL[0-9]+",mar.names[[i]])
      if(length(o)>0) {
        QTL.chr <- c(QTL.chr,rep(i,length(o)))
        QTL.loc <- c(QTL.loc,o)
      }
    }

    # incorporate QTL effects
    for(i in 1:n.qtl) {
      QTL.geno <- geno[[QTL.chr[i]]]$data[,QTL.loc[i]]
      pheno[QTL.geno==1] <- pheno[QTL.geno==1] + model[i,3]
      pheno[QTL.geno==2] <- pheno[QTL.geno==2] + model[i,4]
      pheno[QTL.geno==3] <- pheno[QTL.geno==3] + model[i,5]
    }

  } # end simulate phenotype
      
  n.mar <- sapply(geno, function(a) ncol(a$map))

  # add errors
  if(error.prob > 0) {
    for(i in 1:n.chr) {
      if(chr.type[i] != "X") { # 4-way cross; autosomal
        a <- sample(0:3,n.mar[i]*n.ind,repl=TRUE,
                    prob=c(1-error.prob,rep(error.prob/3,3)))
        if(any(a>0 & geno[[i]]$data==1))
          geno[[i]]$data[a>0 & geno[[i]]$data==1] <-
            geno[[i]]$data[a>0 & geno[[i]]$data==1] + a[a>0 & geno[[i]]$data==1]
        if(any(a>0 & geno[[i]]$data==2))
          geno[[i]]$data[a>0 & geno[[i]]$data==2] <-
            geno[[i]]$data[a>0 & geno[[i]]$data==2] + c(-1,1,2)[a[a>0 & geno[[i]]$data==2]]
        if(any(a>0 & geno[[i]]$data==3))
          geno[[i]]$data[a>0 & geno[[i]]$data==3] <-
            geno[[i]]$data[a>0 & geno[[i]]$data==3] + c(-2,-1,1)[a[a>0 & geno[[i]]$data==3]]
        if(any(a>0 & geno[[i]]$data==4))
          geno[[i]]$data[a>0 & geno[[i]]$data==4] <-
            geno[[i]]$data[a>0 & geno[[i]]$data==4] - a[a>0 & geno[[i]]$data==4]
      }
      else {
        a <- sample(0:1,n.mar[i]*n.ind,repl=TRUE,
                    prob=c(1-error.prob,error.prob))
        if(any(a>0 & geno[[i]]$data==1))
          geno[[i]]$data[a>0 & geno[[i]]$data==1] <-
            geno[[i]]$data[a>0 & geno[[i]]$data==1] + 1
        if(any(a>0 & geno[[i]]$data==2))
          geno[[i]]$data[a>0 & geno[[i]]$data==2] <-
            geno[[i]]$data[a>0 & geno[[i]]$data==2] - 1
        if(any(a>0 & geno[[i]]$data==3))
          geno[[i]]$data[a>0 & geno[[i]]$data==3] <-
            geno[[i]]$data[a>0 & geno[[i]]$data==3] + 1
        if(any(a>0 & geno[[i]]$data==4))
          geno[[i]]$data[a>0 & geno[[i]]$data==4] <-
            geno[[i]]$data[a>0 & geno[[i]]$data==4] - 1
      }

      if(error.ind) {
        errors <- matrix(0,n.ind,n.mar[i])
        errors[a>0] <- 1
        colnames(errors) <- colnames(geno[[i]]$data)
        geno[[i]]$errors <- errors
      }

    } # end loop over chromosomes
  } # end simulate genotyping errors

  # add partial missing
  if(partial.missing.prob > 0) {
    for(i in 1:n.chr) {
      if(chr.type[i] != "X") {
        o <- sample(c(TRUE,FALSE),n.mar[i],repl=TRUE,
                    prob=c(partial.missing.prob,1-partial.missing.prob))

        if(any(o)) {
          o2 <- grep("^QTL[0-9]+",mar.names[[i]])
          if(length(o2)>0)
            x <- geno[[i]]$data[,o2]
          m <- (1:n.mar[i])[o]
          for(j in m) {
            a <- sample(1:4,1)
            if(a==1) { # AB:AA marker
              geno[[i]]$data[geno[[i]]$data[,j]==1 | geno[[i]]$data[,j]==3,j] <- 5
              geno[[i]]$data[geno[[i]]$data[,j]==2 | geno[[i]]$data[,j]==4,j] <- 6
            }
            else if(a==2) { # AA:AB marker
              geno[[i]]$data[geno[[i]]$data[,j]==1 | geno[[i]]$data[,j]==2,j] <- 7
              geno[[i]]$data[geno[[i]]$data[,j]==3 | geno[[i]]$data[,j]==4,j] <- 8
            }
            else if(a==3)  # AB:AB marker
              geno[[i]]$data[geno[[i]]$data[,j]==2 | geno[[i]]$data[,j]==3,j] <- 10
            else  # AB:BA marker
              geno[[i]]$data[geno[[i]]$data[,j]==1 | geno[[i]]$data[,j]==4,j] <- 9
          }
          if(length(o2) > 0)
            geno[[i]]$data[,o2] <- x
        }
      }

    } # end loop over chromosomes
  } # end simulate partially missing data
            
  # add missing
  if(missing.prob > 0) {
    for(i in 1:n.chr) {
      o <- grep("^QTL[0-9]+",mar.names[[i]])
      if(length(o)>0)
        x <- geno[[i]]$data[,o]
      geno[[i]]$data[sample(c(TRUE,FALSE),n.mar[i]*n.ind,repl=TRUE,
                            prob=c(missing.prob,1-missing.prob))] <- NA
      if(length(o)>0)
        geno[[i]]$data[,o] <- x
    }
  }

  if(!is.null(sex)) {
    pheno <- cbind(pheno,sex)
    dimnames(pheno) <- list(NULL, c("phenotype", "sex"))
  }
  else {
    pheno <- cbind(pheno)
    dimnames(pheno) <- list(NULL, "phenotype")
  }

  cross <- list(geno=geno,pheno=pheno)
  class(cross) <- c("4way","cross")

  cross
}



# end of simulate.R
######################################################################
#
# summary.cross.R
#
# copyright (c) 2001, Karl W Broman, Johns Hopkins University
# Oct, 2001; Sept, 2001; Feb, 2001
# Licensed under the GNU General Public License version 2 (June, 1991)
# 
# Part of the R/qtl package
# Contains: summary.cross, print.summary.cross, nind, nchr, nmar,
#           totmar, nphe
#
######################################################################

summary.cross <-
function(object,...)
{
  if(is.na(match("cross",class(object))))
    stop("This is not an object of class cross.")
    
  n.ind <- nind(object)
  tot.mar <- totmar(object)
  n.phe <- nphe(object)
  n.chr <- nchr(object)
  n.mar <- nmar(object)
  type <- class(object)[1]

  Geno <- object$geno[[1]]$data
  if(n.chr > 1)
    for(i in 2:n.chr)
      Geno <- cbind(Geno,object$geno[[i]]$data)

  missing.gen <- mean(is.na(Geno))
  
  if(type=="f2") {
    typings <- table(factor(Geno[!is.na(Geno)], levels=1:5))
    names(typings) <- c("AA","AB","BB","not BB","not AA")
  }
  else if(type=="bc") {
    typings <- table(factor(Geno[!is.na(Geno)], levels=1:2))
    names(typings) <- c("AA","AB")
  }
  else 
    typings <- table(factor(Geno[!is.na(Geno)]))

  typings <- typings/sum(typings)

  missing.phe <- as.numeric(cbind(apply(object$pheno,2,function(a) mean(is.na(a)))))

  # check that object$geno[[i]]$data has colnames and that they match
  #     the names in object$geno[[i]]$map
  for(i in 1:n.chr) {
    nam1 <- colnames(object$geno[[i]]$data)
    map <- object$geno[[i]]$map
    if(is.matrix(map)) nam2 <- colnames(map)
    else nam2 <- names(map)
    chr <- names(object$geno)[[i]]
    if(is.null(nam1)) {
      warn <- paste("The data matrix for chr", chr,
                    "lacks column names")
      warning(warn)
    }
    if(is.null(nam2)) {
      warn <- paste("The genetic map for chr", chr,
                    "lacks column names")
      warning(warn)
    }
    if(any(nam1 != nam2)) {
      warn <- paste("Marker names in the data matrix and genetic map\n",
                    "for chr ", chr, " do not match.",sep="")
      stop(warn)
    }
      
  }
    

  cross.summary <- list(type=type, n.ind = n.ind, n.phe=n.phe, 
			n.chr=n.chr, n.mar=n.mar,
			missing.gen=missing.gen,typing.freq=typings,
			missing.phe=missing.phe)
  class(cross.summary) <- "summary.cross"
  cross.summary
  
}


print.summary.cross <-
function(x,...)
{
  cat("\n")
  if(x$type=="f2") cat("    F2 intercross\n\n")
  else if(x$type=="bc") cat("    Backcross\n\n")
  else if(x$type=="4way") cat("    4-way cross\n\n")
  else cat(paste("    cross", x$type, "\n\n",sep=" "))

  cat("    No. individuals: ", x$n.ind,"\n\n")
  cat("    No. phenotypes:  ", x$n.phe,"\n\n")
  cat("    No. chromosomes: ", x$n.chr,"\n")
  cat("    Total markers:   ", sum(x$n.mar), "\n")
  cat("    No. markers:     ", x$n.mar, "\n")
  cat("\n")
  cat("    Percent genotyped: ", round((1-x$missing.gen)*100,1), "\n")
  cat("    Genotypes (%):     ", 
      paste(names(x$typing.freq),round(x$typing.freq*100,1),sep=":", collapse="  "),
      "\n")
  cat("\n")
  cat("    Percent phenotyped: ", round((1-x$missing.phe)*100,1), "\n")
  cat("\n")
}



nind <-
function(object)
{
  if(any(is.na(match(c("pheno","geno"),names(object)))))
    stop("This is not an object of class cross.")

  n.ind1 <- nrow(object$pheno)
  n.ind2 <- sapply(object$geno,function(x) nrow(x$data))
  if(any(n.ind2 != n.ind1))
    stop("Different numbers of individuals in genotypes and phenotypes.")
  n.ind1
}

nchr <-
function(object)
{
  if(any(is.na(match(c("pheno","geno"),names(object)))))
    stop("This is not an object of class cross.")

  length(object$geno)
}

nmar <- 
function(object)
{
  if(any(is.na(match(c("pheno","geno"),names(object)))))
    stop("This is not an object of class cross.")

  if(!is.matrix(object$geno[[1]]$map))
    n.mar1 <- sapply(object$geno, function(x) length(x$map))
  else # sex-specific maps
    n.mar1 <- sapply(object$geno, function(x) ncol(x$map))

  n.mar2 <- sapply(object$geno, function(x) ncol(x$data))
  if(any(n.mar1 != n.mar2))
    stop("Different numbers of markers in genotypes and maps.")
  n.mar1
}

totmar <-
function(object)
{
  if(any(is.na(match(c("pheno","geno"),names(object)))))
    stop("This is not an object of class cross.")

  if(!is.matrix(object$geno[[1]]$map))
    totmar1 <- sum(sapply(object$geno, function(x) length(x$map)))
  else # sex-specific maps
    totmar1 <- sum(sapply(object$geno, function(x) ncol(x$map)))
  totmar2 <- sum(sapply(object$geno, function(x) ncol(x$data)))
  if(totmar1 != totmar2)
    stop("Different numbers of markers in genotypes and maps.")
  totmar1
}

nphe <-
function(object) {
  if(any(is.na(match(c("pheno","geno"),names(object)))))
    stop("This is not an object of class cross.")

  ncol(object$pheno)
}

# end of summary.cross.R
######################################################################
#
# util.R
#
# copyright (c) 2001, Karl W Broman, Johns Hopkins University
# Oct, 2001; Sept, 2001; July, 2001; Apr, 2001; Feb, 2001
# Licensed under the GNU General Public License version 2 (June, 1991)
# 
# Part of the R/qtl package
# Contains: pull.map, replace.map, pull.chr, create.map,
#           convert.cross, clean, drop.qtlgeno, drop.nullmarkers
#           drop.markers, geno.table, mf.k, mf.h, imf.k, imf.h
#           mf.cf, imf.cf, convert2ss, switch.order
#
######################################################################

######################################################################
#
# pull.map
#
# pull out the map portion of a cross object, as a list
#
######################################################################

pull.map <-
function(cross)
{
  a <- lapply(cross$geno,function(a) { b <- a$map; class(b) <- class(a); b })
  class(a) <- "map"
  a
  
}

######################################################################
#
# replace.map
#
# replace the map portion of a cross object with a list defining a map
#
######################################################################

replace.map <-
function(cross, map)
{
  n.chr <- nchr(cross) 
  n.mar <- nmar(cross)

  n.chr2 <- length(map)
  n.mar2 <- sapply(map,length)

  type <- class(cross)[1]
  if(type=="4way" || type=="f2ss") {
    mnames <- unlist(lapply(cross$geno, function(a) colnames(a$map)))
    mnames2 <- unlist(lapply(map, function(a) colnames(a)))
    n.mar2 <- n.mar2/2
  }
  else if(type == "bc" || type == "f2") {
    mnames <- unlist(lapply(cross$geno, function(a) names(a$map)))
    mnames2 <- unlist(lapply(map, function(a) names(a)))
  }
  else {
    stop(paste("Cross type", type, "not yet supported."))
  }

  # check that things line up properly
  if(n.chr != n.chr2)
    stop("Numbers of chromosomes don't match.")
  if(any(names(cross$geno) != names(map)))
    stop("Chromosome names don't match.")
  if(any(n.mar != n.mar2))
    stop("Number of markers don't match.")
  if(any(mnames != mnames2))
    stop("Marker names don't match.")

  # proceed if no errors
  for(i in 1:length(cross$geno))
    cross$geno[[i]]$map <- map[[i]]

  cross
}



######################################################################
#
# pull.chr
#
# Pull out a portion of the chromosomes from a cross object
#
######################################################################

pull.chr <-
function(cross, chr)
{
  if(is.numeric(chr)) {
    if(any(chr < 1 | chr > nchr(cross)))
      stop("Chromosome numbers out of range.")
  }
  else {
    if(any(is.na(match(chr,names(cross$geno)))))
      stop("Not all chromosome names found.")
  }

  if(!is.na(match("rf",names(cross)))) { # pull out part of rec fracs
    n.mar <- nmar(cross)
    n.chr <- nchr(cross)
    wh <- rbind(c(0,cumsum(n.mar)[-n.chr])+1,cumsum(n.mar))
    dimnames(wh) <- list(NULL, names(n.mar))
    wh <- as.matrix(wh[,chr])
    wh <- unlist(apply(wh,2,function(a) a[1]:a[2]))
    cross$rf <- cross$rf[wh,wh]
  }

  cross$geno <- cross$geno[chr]
  cross
}



######################################################################
#
# create.map
#
# create a new map with inserted inter-marker locations
#
# Note: map is a vector or a matrix with 2 rows
# 
######################################################################

create.map <-
function(map, step, off.end)
{
  if(!is.matrix(map)) { # sex-ave map
    if(step==0 && off.end==0) return(map)
    else if(step==0 && off.end > 0) {
      a <- c(floor(min(map)-off.end),ceiling(max(map)+off.end))
      names(a) <- paste("loc", a, sep="")
      return(sort(c(a,map)))
    }
    else if(step>0 && off.end == 0) {
      a <- seq(floor(min(map)),max(map),
               by = step)
      if(any(is.na(match(a,map)))) {
        a <- a[is.na(match(a,map))]
        names(a) <- paste("loc",a,sep="")
        return(sort(c(a,map)))
      }
      else return(map)
    }
    else {
      a <- seq(floor(min(map)-off.end),ceiling(max(map)+off.end+step),
               by = step)
      a <- a[is.na(match(a,map))]
      
      # no more than one point above max(map)+off.end
      z <- (1:length(a))[a >= max(map)+off.end]
      if(length(z) > 1) a <- a[-z[-1]]
      
      names(a) <- paste("loc",a,sep="")
      return(sort(c(a,map)))
    }
  } # end sex-ave map
  else { # sex-specific map
    if(step==0 && off.end==0) return(map)
    else if(step==0 && off.end > 0) {
      L1 <- diff(range(map[1,]))
      L2 <- diff(range(map[2,]))
      a <- c(floor(min(map[1,])-off.end),ceiling(max(map[1,])+off.end))
      names(a) <- paste("loc", a, sep="")
      b <- c(floor(min(map[2,])-off.end)*L2/L1,
             ceiling(max(map[1,])+off.end)*L2/L1)
      n <- c(names(a)[1],colnames(map),names(a)[2])
      map <- cbind(c(a[1],b[1]),map,c(a[2],b[2]))
      dimnames(map) <- list(NULL,n)
      return(map)
    }
    else if(step>0 && off.end == 0) {
      a <- seq(floor(min(map[1,])),max(map[1,]),
               by = step)
      a <- a[is.na(match(a,map[1,]))]
      names(a) <- paste("loc",a,sep="")
      b <- sapply(a,function(x,y,z) {
          I <- min((1:length(y))[y > x])
          (x-y[I-1])/(y[I]-y[I-1])*(z[I]-z[I-1])+z[I-1] }, map[1,],map[2,])
      names(b) <- names(a)
      return(rbind(sort(c(a,map[1,])),sort(c(b,map[2,]))))
    }
    else {
      a <- seq(floor(min(map[1,])-off.end),ceiling(max(map[1,])+off.end+step),
               by = step)
      a <- a[is.na(match(a,map[1,]))]
      # no more than one point above max(map)+off.end
      z <- (1:length(a))[a >= max(map[1,])+off.end]
      if(length(z) > 1) a <- a[-z[-1]]
      names(a) <- paste("loc",a,sep="")

      b <- sapply(a,function(x,y,z) {
        if(x < min(y))
          return( min(z) - (min(y)-x)/diff(range(y))*diff(range(z)) )
        else if(x > max(y))
          return( max(z) + (x - max(y))/diff(range(y))*diff(range(z)) )
        else {
          I <- min((1:length(y))[y > x])
          (x-y[I-1])/(y[I]-y[I-1])*(z[I]-z[I-1])+z[I-1]
        }
        }, map[1,],map[2,])
      names(b) <- names(a)

      return(rbind(sort(c(a,map[1,])), sort(c(b,map[2,]))))
    }
    

  }
}

  
  
######################################################################
#
# convert.cross: convert a "qtl.cross" data set from the format
#                used in old versions (<= 0.65) of R/qtl to the
#                updated data structure (versions >= 0.70).
#
######################################################################

convert.cross <-
function(cross)   
{
  require(qtl)
  nchr <- length(cross$map)
  geno <- vector("list",nchr)
  nmar <- c(0,cumsum(sapply(cross$map,length)))
  for(i in 1:nchr) {
    whichpos <- (nmar[i]+1):nmar[i+1]
    geno[[i]] <- list(data=cross$geno[,whichpos],map=cross$map[[i]])
    dimnames(geno[[i]]$data) <- list(NULL, names(cross$map[[i]]))
    class(geno[[i]]) <- "A"
    chr.name <- names(cross$map)[i]
    if(chr.name == "X" || chr.name == "x")
      class(geno[[i]]) <- "X"
  }
  names(geno) <- names(cross$map)
  cross$geno <- geno
  type <- cross$type
  cross <- cross[1:2]
  class(cross) <- c(type,"cross")
  cross
}

    

######################################################################
#
# clean
# 
# remove all of the extraneous stuff from a cross object, to get back
# to just the data
#
######################################################################

clean <-
function(cross)
{
  cross2 <- list(geno=cross$geno,pheno=cross$pheno)

  for(i in 1:length(cross$geno)) {
    cross2$geno[[i]] <- list(data=cross$geno[[i]]$data,
                             map=cross$geno[[i]]$map)
    class(cross2$geno[[i]]) <- class(cross$geno[[i]])
  }
    
  class(cross2) <- class(cross)
  cross2
}


######################################################################
#
# drop.qtlgeno
#
# remove any QTLs from the genotype data and the genetic maps
# from data simulated via sim.cross. (They all have names "QTL*")
#
######################################################################

drop.qtlgeno <-
function(cross)  
{
  n.chr <- nchr(cross)
  mar.names <- lapply(cross$geno, function(a) {
    m <- a$map
    if(is.matrix(m)) return(colnames(m))
    else return(names(m)) } )
    
  for(i in 1:n.chr) {
    o <- grep("^QTL[0-9]+",mar.names[[i]])
    if(length(o) != 0) {
      cross$geno[[i]]$data <- cross$geno[[i]]$data[,-o]
      if(is.matrix(cross$geno[[i]]$map)) 
        cross$geno[[i]]$map <- cross$geno[[i]]$map[,-o]
      else
        cross$geno[[i]]$map <- cross$geno[[i]]$map[-o]
    }
  }
  cross
}

######################################################################
#
# drop.nullmarkers
#
# remove markers that have no genotype data from the data matrix and
# genetic maps
#
######################################################################

drop.nullmarkers <-
function(cross)
{
  n.chr <- nchr(cross)

  keep.chr <- rep(TRUE,n.chr)
  for(i in 1:n.chr) {
    o <- !apply(cross$geno[[i]]$data,2,function(a) sum(!is.na(a)))
    if(any(o)) { # remove from genotype data and map
      mn.drop <- colnames(cross$geno[[i]]$data)[o]
      if(length(mn.drop) == ncol(cross$geno[[i]]$data)) 
        keep.chr[i] <- FALSE # removing all markers from this chromosome

      if(sum(!o) == 1) mn <- colnames(cross$geno[[i]]$data)[!o]

      cross$geno[[i]]$data <- cross$geno[[i]]$data[,!o]

      if(is.matrix(cross$geno[[i]]$map)) {
        if(sum(!o) == 1) {
          x <- as.matrix(cross$geno[[i]]$map[,!o])
          colnames(x) <- mn
        }
        else 
          cross$geno[[i]]$map <- cross$geno[[i]]$map[,!o]
      }
      else 
        cross$geno[[i]]$map <- cross$geno[[i]]$map[!o]

      if(sum(!o) == 1) {
        cross$geno[[i]]$data <- as.matrix(cross$geno[[i]]$data)
        colnames(cross$geno[[i]]$data) <- mn
      }

      # results of calc.genoprob
      if(!is.na(match("prob",names(cross$geno[[i]])))) {
        o <- match(mn.drop,colnames(cross$geno[[i]]$prob))
        cross$geno[[i]]$prob <- cross$geno[[i]]$prob[,-o,]
      }

      # results of argmax.geno
      if(!is.na(match("argmax",names(cross$geno[[i]])))) {
        o <- match(mn.drop,colnames(cross$geno[[i]]$argmax))
        cross$geno[[i]]$argmax <- cross$geno[[i]]$argmax[,-o]
      }

      # results of sim.geno
      if(!is.na(match("draws",names(cross$geno[[i]])))) {
        o <- match(mn.drop,colnames(cross$geno[[i]]$draws))
        cross$geno[[i]]$draws <- cross$geno[[i]]$draws[,-o,]
      }

      # results of est.rf
      if(!is.na(match("rf",names(cross)))) {
        o <- match(mn.drop,colnames(cross$rf))
        cross$rf <- cross$rf[-o,-o]
      }
    }
  }

  cross$geno <- cross$geno[keep.chr]

  cross
}

    
######################################################################
#
# drop.markers
#
# remove a vector of markers from the data matrix and genetic maps
#
######################################################################

drop.markers <-
function(cross,markers)
{
  n.chr <- nchr(cross)

  keep.chr <- rep(TRUE,n.chr)
  for(i in 1:n.chr) {
    # find markers on this chromosome
    o <- match(markers,colnames(cross$geno[[i]]$data))
    o <- o[!is.na(o)]
    a <- rep(FALSE,ncol(cross$geno[[i]]$data))
    a[o] <- TRUE
    o <- a
    
    if(any(o)) { # remove from genotype data and map
      mn.drop <- colnames(cross$geno[[i]]$data)[o]
      if(length(mn.drop) == ncol(cross$geno[[i]]$data)) 
        keep.chr[i] <- FALSE # removing all markers from this chromosome

      if(sum(!o) == 1) mn <- colnames(cross$geno[[i]]$data)[!o]

      cross$geno[[i]]$data <- cross$geno[[i]]$data[,!o]

      if(is.matrix(cross$geno[[i]]$map)) {
        if(sum(!o) == 1) {
          x <- as.matrix(cross$geno[[i]]$map[,!o])
          colnames(x) <- mn
        }
        else 
          cross$geno[[i]]$map <- cross$geno[[i]]$map[,!o]
      }
      else 
        cross$geno[[i]]$map <- cross$geno[[i]]$map[!o]

      if(sum(!o) == 1) {
        cross$geno[[i]]$data <- as.matrix(cross$geno[[i]]$data)
        colnames(cross$geno[[i]]$data) <- mn
      }

      # results of calc.genoprob
      if(!is.na(match("prob",names(cross$geno[[i]])))) {
        o <- match(mn.drop,colnames(cross$geno[[i]]$prob))
        cross$geno[[i]]$prob <- cross$geno[[i]]$prob[,-o,]
      }

      # results of argmax.geno
      if(!is.na(match("argmax",names(cross$geno[[i]])))) {
        o <- match(mn.drop,colnames(cross$geno[[i]]$argmax))
        cross$geno[[i]]$argmax <- cross$geno[[i]]$argmax[,-o]
      }

      # results of sim.geno
      if(!is.na(match("draws",names(cross$geno[[i]])))) {
        o <- match(mn.drop,colnames(cross$geno[[i]]$draws))
        cross$geno[[i]]$draws <- cross$geno[[i]]$draws[,-o,]
      }

      # results of est.rf
      if(!is.na(match("rf",names(cross)))) {
        o <- match(mn.drop,colnames(cross$rf))
        cross$rf <- cross$rf[-o,-o]

      }
    }
  }

  cross$geno <- cross$geno[keep.chr]

  cross
}

    
######################################################################
#
# geno.table
#
# create table showing observed numbers of individuals with each
# genotype at each marker
#
######################################################################

geno.table <- 
function(cross)
{
  n.chr <- nchr(cross)

  type <- class(cross)[1]
  if(type == "f2") {
    n.gen <- 5
    gen.names <- c("AA","AB","BB","AA/AB","AB/BB")
  }
  else if(type == "bc") {
    n.gen <- 2
    gen.names <- c("AA","AB")
  }
  else if(type == "4way") {
    n.gen <- 10
    gen.names <- c("AC","BC","AD","BD","AC/AD","BC/BD",
                   "AC/BC","AD/BD","AC/BD","AD/BC")
  }
  else
    stop(paste("Unknown cross type:",type))
    
  res <- lapply(cross$geno, function(a,ngen) {
                a <- a$data; a[is.na(a)] <- 0
                apply(a,2,function(b,ngen) table(factor(b,levels=0:ngen)),ngen)
                },n.gen)
  results <- NULL
  for(i in 1:length(res)) 
    results <- rbind(results,t(res[[i]]))
  colnames(results) <- c("NA",gen.names)
  rownames(results) <- unlist(lapply(cross$geno,function(a) colnames(a$data)))
  results
}
    
# map functions
mf.k <- function(d) 0.5*tanh(d/50)
mf.h <- function(d) 0.5*(1-exp(-d/50))
imf.k <- function(r) 50*atanh(2*r)
imf.h <- function(r) -50*log(1-2*r)

# carter-falconer: mf.cf, imf.cf
imf.cf <- function(r) 12.5*(log(1+2*r)-log(1-2*r))+25*atan(2*r)

mf.cf <-
function(d)
{
  icf <- function(r,d)
    imf.cf(r)-d

  sapply(d,function(a) {
    if(a==0) return(0)
    uniroot(icf, c(0,0.5-1e-10),d=a)$root })
}


# convert F2 intercross to sex-specific
convert2ss <-
function(cross)  
{
  if(class(cross)[1] != "f2")
    stop("This function applies only to f2 crosses.")

  class(cross)[1] <- "f2ss"

  for(i in 1:nchr(cross)) 
    cross$geno[[i]]$map <- rbind(cross$geno[[i]]$map,
                                 cross$geno[[i]]$map)

  cross
}

######################################################################
#
# switch.order: change the marker order on a given chromosome to some
#               specified order
#
######################################################################

switch.order <-
function(cross, chr, order)
{
  # check chromosome argument
  if(is.character(chr)) {
    old.chr <- chr
    chr <- match(chr, names(cross$geno))
    if(length(chr) > 1) chr <- chr[1]
    if(is.na(chr))
      stop("There is no chromosome named", chr)
  }

  # check order argument
  n.mar <- nmar(cross)
  if(n.mar[chr] == length(order)-2) # useful for output from ripple()
    order <- order[1:n.mar[chr]]
  if(n.mar[chr] != length(order))
    stop("Incorrect number of markers.")

  # remove any intermediate calculations, as they
  #   will no longer be meaningful
  cross <- clean(cross)

  # re-order markers
  cross$geno[[chr]]$data <- cross$geno[[chr]]$data[,order]
  m <- seq(0,by=5,length=ncol(cross$geno[[chr]]$data))
  names(m) <- colnames(cross$geno[[chr]]$data)
  if(is.matrix(cross$geno[[chr]]$map)) 
    cross$geno[[chr]]$map <- rbind(m,m)
  else
    cross$geno[[chr]]$map <- m

  # re-estimate map
  newmap <- est.map(pull.chr(cross,chr))
  cross$geno[[chr]]$map <- newmap[[1]]

  cross
}

# end of util.R
######################################################################
#
# vbscan.R
#
# copyright (c) 2001, Karl W Broman, Johns Hopkins University
# July, 2001; May, 2001
# Licensed under the GNU General Public License version 2 (June, 1991)
# 
# Part of the R/qtl package
# Contains: vbscan, vbscan.perm
#
######################################################################

######################################################################
#
# vbscan: scan genome for a quantitative phenotype for which some
# individuals' phenotype is undefined (for example, the size of a
# lesion, where some individuals have no lesion).
#
######################################################################

vbscan <-
function(cross, chr, pheno.col=1, upper=FALSE,
	 maxit=1000, tol=1e-8)
{
  if(!missing(chr)) cross <- pull.chr(cross, chr)

  # check arguments are okay
  if(length(pheno.col) > 1) pheno.col <- pheno.col[1]
  if(pheno.col > nphe(cross))
    stop("Specified phenotype column exceeds the number of phenotypes")
  y <- cross$pheno[,pheno.col]

  # remove individuals with missing phenotypes
  if(any(is.na(y))) {
    drop <- (1:length(y))[is.na(y)]
    y <- y[-drop]
    for(i in 1:length(cross$geno)) 
      cross$geno[[i]]$prob <- cross$geno[[i]]$prob[-drop,,]
  }

  # modify phenotypes
  if(upper) {
    if(!any(y == Inf)) y[y==max(y)] <- Inf
  }
  else {
    if(!any(y == -Inf)) y[y==min(y)] <- -Inf
  }
  survived <- rep(0,length(y))
  survived[y == -Inf | y == Inf] <- 1

  # The following line is included since .C() doesn't accept Infs
  y[y == -Inf | y == Inf] <- 99999

  results <- NULL
  for(i in 1:length(cross$geno)) {
    # make sure inferred genotypes or genotype probabilities are available
    if(is.na(match("prob",names(cross$geno[[i]])))) {
      cat(" -Calculating genotype probabilities\n")
      cross <- calc.genoprob(cross)
    }

    n.pos <- dim(cross$geno[[i]]$prob)[2]
    n.ind <- length(y)
    n.gen <- dim(cross$geno[[i]]$prob)[3]

    z <- .C("R_vbscan",
            as.integer(n.pos),
            as.integer(n.ind),
            as.integer(n.gen),
            as.double(cross$geno[[i]]$prob),
            as.double(y),
            as.integer(survived),
            lod=as.double(rep(0,(4+2*n.gen)*n.pos)),
            as.integer(maxit),
            as.double(tol))

    map <- create.map(cross$geno[[i]]$map,
                      attr(cross$geno[[i]]$prob,"step"),
                      attr(cross$geno[[i]]$prob,"off.end"))
    if(is.matrix(map)) map <- map[1,]

    res <- data.frame(chr=rep(names(cross$geno)[i],length(map)),
                      pos = map,
                      matrix(z$lod,nrow=n.pos,byrow=TRUE))

    w <- names(map)
    o <- grep("^loc\-*[0-9]+",w)
    if(length(o) > 0)
      w[o] <- paste(w[o],names(cross$geno)[i],sep=".c")
    rownames(res) <- w
    
    colnames(res) <- c("chr","pos","lod","lod.p","lod.mu",
                       paste("pi",c("AA","AB","BB")[1:n.gen]),
                       paste("mu",c("AA","AB","BB")[1:n.gen]),
                       "sigma")

    # the following is for the case where there is a sex chromosome
    if(!is.null(results) && ncol(results) != ncol(res)) {
      if(ncol(results) > ncol(res)) {
        o <- match(colnames(res),colnames(results))
        missing <- (1:ncol(results))[-o]
        temp <- as.data.frame(matrix(ncol=ncol(results),nrow=nrow(res)))
        colnames(temp) <- colnames(results)
        rownames(temp) <- rownames(res)
        temp[,o] <- res
        res <- temp
      }
      else {
        o <- match(colnames(results),colnames(res))
        missing <- (1:ncol(res))[-o]
        temp <- as.data.frame(matrix(ncol=ncol(res),nrow=nrow(results)))
        colnames(temp) <- colnames(res)
        rownames(temp) <- rownames(results)
        temp[,o] <- results
        results <- temp
      }
    }
    
    results <- rbind(results,res)
  }
  
  class(results) <- c("scanone","data.frame")
  results
}


# permutation test for vbscan
vbscan.perm <-
function(cross, chr, pheno.col=1, upper=FALSE,
	 n.perm=1000, maxit=1000, tol=1e-8)
{
  if(!missing(chr)) cross <- pull.chr(cross, chr)

  n.ind <- nind(cross)

  res <- matrix(ncol=3,nrow=n.perm)
  for(i in 1:n.perm) {
    cross$pheno <- as.matrix(cross$pheno[sample(1:n.ind),])
    tem <- vbscan(cross,pheno.col=pheno.col,upper=upper,
                      maxit=maxit,tol=tol)
#    print(c(i,dim(tem)))
    res[i,] <- apply(tem[,3:5],2,max)
  }
  colnames(res) <- c("lod","lod.p","lod.mu")
  res
}
######################################################################
#
# write.cross.R
#
# copyright (c) 2000-2001, Karl W Broman, Johns Hopkins University
# Oct, 2001; Sept, 2001; Feb, 2001
# Licensed under the GNU General Public License version 2 (June, 1991)
#
# Part of the R/qtl package
# Contains: write.cross, write.cross.mm, write.cross.csv
#
######################################################################


######################################################################
#
# write.cross: Wrapper for the other write.cross functions
#
######################################################################

write.cross <-
function(cross, format=c("csv","mm"), filestem="data", chr, digits=5)
{
  match.arg <- match.arg(format)
  if(missing(chr)) chr <- 1:nchr(cross)

  if(format=="csv") write.cross.csv(cross,filestem,chr,digits)
  else write.cross.mm(cross,filestem,chr,digits)
}


######################################################################
#
# write.cross.mm: Write data for an experimental cross in Mapmaker
#                 format
#
#           creates two files: "raw" file with geno & pheno data
#                              "prep" file with map information
#
######################################################################

write.cross.mm <-
function(cross, filestem="data", chr, digits=5)
{
  if(!missing(chr)) 
    cross <- pull.chr(cross,chr)

  n.ind <- nind(cross)
  tot.mar <- totmar(cross)
  n.phe <- nphe(cross)
  n.chr <- nchr(cross)
  n.mar <- nmar(cross)
  
  type <- class(cross)[1]
  if(type != "f2" && type != "bc")
    stop("write.cross.mm only works for intercross and backcross data.")

  # write genotype and phenotype data
  file <- paste(filestem, ".raw", sep="")
  
  # write experiment type
  if(type == "f2") 
    write("data type f2 intercross", file, append=FALSE)
  else 
    write("data type f2 backcross", file, append=FALSE)

  # write numbers of progeny, markers and phenotypes
  write(paste(n.ind, tot.mar, n.phe), file, append=TRUE)

  # max length of marker name
  mlmn <- max(nchar(unlist(lapply(cross$geno,function(a) colnames(a$data)))))+1

  # write marker data
  for(i in 1:n.chr) {
    for(j in 1:n.mar[i]) {
      mn <- paste("*", colnames(cross$geno[[i]]$data)[j], sep="")
      if(nchar(mn) < mlmn)
        mn <- paste(mn,paste(rep(" ", mlmn-nchar(mn)),collapse=""),sep="")
      g <- cross$geno[[i]]$data[,j]
      
      x <- rep("", n.ind)
      x[is.na(g)] <- "-"
      x[!is.na(g) & g==1] <- "A"
      x[!is.na(g) & g==2] <- "H"
      if(type == "f2") {
        x[!is.na(g) & g==3] <- "B"
        x[!is.na(g) & g==4] <- "D"
        x[!is.na(g) & g==5] <- "C"
      }

      if(n.ind < 60)
        write(paste(mn, paste(x,collapse="")), file, append=TRUE)
      else {
        lo <- seq(1,n.ind-1,by=60)
        hi <- c(lo[-1]-1,n.ind)
        for(k in 1:length(lo)) {
          if(k==1) write(paste(mn,paste(x[lo[k]:hi[k]],collapse="")),file,append=TRUE)
          else write(paste(paste(rep(" ", mlmn),collapse=""),
                           paste(x[lo[k]:hi[k]],collapse="")),file,append=TRUE)
        }
      }

    }
  } # end writing marker data

  # max length of phenotype name
  mlpn <- max(nchar(colnames(cross$pheno)))+1

  # write phenotypes
  for(i in 1:n.phe) {
    pn <- paste("*",colnames(cross$pheno)[i],sep="")
    if(nchar(pn) < mlpn)
      pn <- paste(pn, paste(rep(" ", mlpn-nchar(pn)),collapse=""),sep="")

    x <- as.character(round(cross$pheno[,i],digits))
    x[x=="NA"] <- "-"

    if(n.ind < 10)
      write(paste(pn, paste(x,collapse="")), file, append=TRUE)
    else {
      lo <- seq(1,n.ind-1,by=10)
      hi <- c(lo[-1]+1,n.ind)
      for(k in 1:length(lo)) {
        if(k==1) write(paste(pn,paste(x[lo[k]:hi[k]],collapse=" ")),file,append=TRUE)
        else write(paste(paste(rep(" ", mlpn),collapse=""),
                         paste(x[lo[k]:hi[k]],collapse=" ")),file,append=TRUE)
      }
    }
  }
    

  # make "prep" file with map information
  file <- paste(filestem, ".prep", sep="")

  for(i in 1:n.chr) {
    cname <- paste("chr", names(cross$geno)[i], sep="")
    line <- paste("make chromosome", cname)
    if(i==1) write(line, file, append=FALSE)
    else write(line, file, append=TRUE)

    mn <- names(cross$geno[[i]]$map)
    dis <- round(diff(cross$geno[[i]]$map),2)
    dis <- paste("=", dis, sep="")
    write(paste(paste("sequence", mn[1]), paste(dis,mn[-1],collapse=" ")),
          file, append=TRUE)

    write(paste("anchor", cname), file, append=TRUE)
    

  }

} 

######################################################################
#
# write.cross.csv: Write data for an experimental cross in
#                  comma-delimited format (the same format as is read
#                  by read.cross.csv)
#
######################################################################

write.cross.csv <-
function(cross, filestem="data", chr, digits=5)
{
  if(!missing(chr)) cross <- pull.chr(cross,chr)

  n.ind <- nind(cross)
  tot.mar <- totmar(cross)
  n.phe <- nphe(cross)
  n.chr <- nchr(cross)
  n.mar <- nmar(cross)
  
  type <- class(cross)[1]
  if(type != "f2" && type != "bc")
    stop("write.cross.csv only works for intercross and backcross data.")

  file <- paste(filestem, ".csv", sep="")
  
  geno <- matrix(ncol=tot.mar,nrow=n.ind)
  alleles <- c("A","H","B","D","C")
  firstmar <- 1
  for(i in 1:n.chr) {
    # replace allele numbers with 
    geno[,firstmar:(firstmar+n.mar[i]-1)] <-
      alleles[match(cross$geno[[i]]$data,1:5)]
    firstmar <- firstmar + n.mar[i]
  }
  geno[geno=="NA"] <- "-"
  data <- cbind(matrix(as.character(round(cross$pheno,digits)),ncol=n.phe),geno)
  colnames(data) <- c(colnames(cross$pheno),unlist(lapply(cross$geno,function(a) colnames(a$data))))
  chr <- rep(names(cross$geno),n.mar)
  pos <- unlist(lapply(cross$geno,function(a) a$map))
  chr <- c(rep("",n.phe),chr)
  pos <- c(rep("",n.phe),as.character(round(pos,digits)))

  # write names
  write.table(matrix(colnames(data),nrow=1),file,append=FALSE,quote=FALSE,
              sep=",",row.names=FALSE,col.names=FALSE)
  # write chr IDs
  write.table(matrix(chr,nrow=1),file,append=TRUE,quote=FALSE,sep=",",
              row.names=FALSE,col.names=FALSE)
  # write marker positions
  write.table(matrix(pos,nrow=1),file,append=TRUE,quote=FALSE,sep=",",
              row.names=FALSE,col.names=FALSE)
  # write phenotype and genotype data
  write.table(data,file,append=TRUE,quote=FALSE,sep=",",row.names=FALSE,
              col.names=FALSE)

}

# end of write.cross.R 
######################################################################
#
# zzz.R
#
# copyright (c) 2001, Karl W Broman, Johns Hopkins University
# Feb, 2001
# Licensed under the GNU General Public License version 2 (June, 1991)
#
# Part of the R/qtl package
#
# .First.lib is run when the package is loaded with library(qtl)
#
######################################################################

.First.lib <- function(lib, pkg) library.dynam("qtl", pkg, lib)

# end of zzz.R
