"add"<-
    function(x,...,scope,method=method)
{
    if (missing(scope)) stop("missing scope arg - possible bqtl syntax error")
    new.scope <- scope[grep("^add",scope)]
    if (length(new.scope)==0) stop("no 'add' terms found")
    locus(x,...,
          scope=new.scope,
          method="fooey")    # this forces one term per locus
}
"adjust.linear.bayes" <-
    function(lbo,ana.obj=lbo$call$ana.obj,...)
{
    this.call <- match.call(expand.dots=TRUE)
    specs <- lbo$specs

### get twohk and swap args from lbo
    cl.1 <- lbo$hk$call
    cl.1[[1]]<-cl.1$varcov<-cl.1$ana.obj<-NULL
    cl.2 <- lbo$varcov.call
    cl.2[[1]]<-cl.2$x<-cl.2$ana.obj<-NULL

### this builds up the call to bqtl using args from
### the original call to linear.bayes. It figures out which to bother
### with by inspecting the calls to varcov and twohk as invoked by
### linear.bayes      
    one.gene.expr <- call("bqtl",reg.formula=lbo$varcov.call$x,
                          ana.obj=ana.obj)
### are there extra args to be added ??
    extra.args <- unique(c(names(cl.1),names(cl.2)))
    if (length(extra.args)!=0){
        dummy.call <- expression("$<-"(one.gene.expr,arg.2,arg.3))[[1]]
        for ( i in extra.args ){
            dummy.call[[3]] <- dummy.call[[4]] <- i
            eval(dummy.call)
            one.gene.expr[[i]]<-lbo$call[[i]]
        }
    }
    
### build call to summary.adj
    one.gene.smry.expr <-
        call("summary.adj",adj.obj=as.name("one.gene.models"),
             n.loc=1,coef.znames="",mode.names=NULL,imp.denom="")
    one.gene.smry.expr$coef.znames <- call("$",ana.obj,"reg.names")
    mf.expr <- call("$",ana.obj,"map.frame")
    prior.expr <- call("$",mf.expr,"prior")
    one.gene.smry.expr$imp.denom <- call("/",1,prior.expr)
    
    if ( any(specs$gene.number == 1) ) {
        one.gene.models <- eval(one.gene.expr)
        one.gene.smry <- eval(one.gene.smry.expr)
    }
    else {
        one.gene.smry <- NULL
    }
### now do multi-gene models
    swap.gene.expr <- one.gene.expr
    swaps.expr <-
        expression( configs( unique.config(lbo$swaps[[1]])$uniq ) )[[1]]
    swaps.expr[[2]][[2]][[2]][[2]][[2]] <- as.name(this.call$lbo)
    n.gene.smry.expr <- one.gene.smry.expr
    n.gene.smry.expr$adj.obj <- as.name("swaps.adj")
    n.gene.smry.expr$imp.denom <- NULL
    n.gene <- vector("list",length(lbo$swaps))
    n.swap <- specs$gene.number[specs$n.cycles != 0]

    for( i in n.swap ){
        ## y ~ configs( unique.config(lbo$swaps[[  i  ]])$uniq )
        swaps.expr[[2]][[2]][[2]][[3]] <- i
        swap.gene.expr$reg.formula[[3]] <- swaps.expr
        swaps.adj <- eval(swap.gene.expr)  # do it
        ## now summarize
        n.gene.smry.expr$n.loc <- i
        n.gene.smry.expr$swap.obj <- swaps.expr[[2]][[2]][[2]]
        n.gene[[i]] <- c( eval(n.gene.smry.expr), list(swaps=swaps.adj))
    }
    if ( is.null(lbo$odds) ){
        odds <- loc.posterior <- coefs <- NULL
    }
    else
        if ( any(specs$gene.number == 1) ){
            odds <-
                lbo$odds*cumprod(c(one.gene.smry$adj,
                                   sapply(n.gene[n.swap],"[[","adj")))
            post.pr <- odds/sum(odds)
            loc.posterior <-
                cbind(one.gene.smry$loc,sapply(n.gene[n.swap],"[[","loc")) %*%
                    ( post.pr* c(1,n.swap))
            coefs <-
                cbind(one.gene.smry$coef,sapply(n.gene[n.swap],"[[","coef")) %*%
                    post.pr
        }
        else {
            odds <-
                lbo$odds*cumprod( sapply(n.gene[n.swap],"[[","adj") )/
                    n.gene[n.swap][[1]]$adj  #take first as baseline for odds
            post.pr <- odds/sum(odds)
            loc.posterior <-
                sapply(n.gene[n.swap],"[[","loc") %*%
                    ( post.pr * n.swap )
            coefs <-
                sapply(n.gene[n.swap],"[[","coef") %*%
                    post.pr
        } 
    
    res <- list(odds=odds,loc.posterior=loc.posterior,coefficients=coefs,
                one.gene.adj=one.gene.smry,n.gene.adj=n.gene,
                call=this.call)
    class(res) <- "adjust.linear.bayes"
    
    res
}
"bqtl"<-
    function(reg.formula,ana.obj=analysis.object, 
             scope = ana.obj$reg.names, expand.specials= NULL, ...)
{
    local.covar <-
        function (x, ..., scope = scope, method = method,
                  prefix = NULL, bq.spec=bqtl.specials,chromo=NULL,cM=NULL) 
        {
            if (missing(scope)) 
                stop("missing scope arg - possible bqtl syntax error")
            dots <- list(...)
            if (is.null(chromo)&& length(names(dots))!=0 &&
                any(cloc <- which(1==pmatch(names(dots),"chromo",0)))){
                chromo <- dots[[cloc]]
                if (length(dots)>1)
                    stop("... not allowed")
                else
                    dots<-NULL
            }
            x.arg <- match.call()$x
            if  ( length(x.arg) ==1 && deparse(x.arg) == "all" )  # allow 'all' shorthand 
                x <- seq(along=scope)
            if (is.call(x.arg)){
                if (deparse(x.arg[[1]])%in%bq.spec){
                    if (is.null(x.arg$scope)) x.arg$scope <- as.name("scope")
                    if (is.null(x.arg$method)) x.arg$method <- method
                    if (is.null(x.arg$chromo)) x.arg$chromo <- chromo
                    if (is.null(x.arg$cM)) x.arg$cM <- cM
                    
                    if (!missing(...)) stop("cannot use ... in this context")
                    paste("covar(",eval(x.arg),")",sep="")
                }
                else{ ## x.arg is c(1,2) or 7:8 or whatever
                    x <- eval(x.arg)
                    new.scope <- paste("covar(", scope, ")", sep = "")
                    locus(x, ..., scope = new.scope, method = method,
                          chromo=chromo,cM=cM)
                }
            }
            else {
                deparse.x <- deparse(x.arg)
                if (deparse.x%in%scope){
                    if (!missing(...))
                        stop("cannot use ... args with named locus")
                    else
                        paste("covar(",deparse.x,")",sep="")
                }
                else {
                    new.scope <- paste("covar(", scope, ")", sep = "")
                    locus(x, ..., scope = new.scope, method = method,
                          chromo=chromo,cM=cM)
                }
            }
        }
    
    this.call <- match.call()
    if (any(is.na(pmatch(c("loc.right","map.frame","state.matrix","method"),names(ana.obj)))))
        stop("ana.obj doesn't have required components")

    using.R <- exists("is.R")&&is.R()

### get terms.object, extract 'configs' terms,
    
    bqtl.specials <- c("configs","locus","add","dom","covar","acovar","dcovar")
    reg.terms <- terms(reg.formula,specials=bqtl.specials)
    reg.labels <- labels(reg.terms)
    reg.specials <- attr(reg.terms,"specials")  
    if (length(unlist(reg.specials))==0) { # no specials - pass thru to lapadj
        res <- 
            lapadj(reg.formula,ana.obj,...)
        res$call <- attr(res,"call")<- this.call
        class(res)<- "bqtl"
        return(res)
    }

    method <- ana.obj$method
### work on attr(,"factors")
    reg.factor <- attr(reg.terms,"factors")
    term.rownames <- dimnames(reg.factor)[[1]]
    term.names <- term.rownames[row(reg.factor)[reg.factor>0]]
    terms.conjuncs <-
        if ( sum(reg.factor>0) <2 )
            NULL
        else
            ifelse( diff(col(reg.factor)[reg.factor>0])==0,"colon", "plus")
    n.join <- length(terms.conjuncs)
### flag rows with specials and those without
    reg.specials <- unlist(reg.specials[ !sapply(reg.specials,is.null) ])
    reg.plain <- term.rownames[ - reg.specials ]
    names(reg.plain) <- reg.plain
    response  <- term.rownames[1]
    reg.plain[response] <- paste(response, "~")
    
### let the specials expand themselves
    if (using.R)
        formals(local.covar)$bq.spec <- bqtl.specials # bind bqtl.specials explicitly
    else
        local.covar$bq.spec <- bqtl.specials
    pt.vars <- reg.specials + if (using.R) 1 else 0
    rspec <-
        lapply(attr(reg.terms,"variables")[pt.vars],
               function(x,scope,method,covar) {
                   if ( !is.element("scope",names(x)) ) #typically use default
                       x$scope <- as.name("scope")
                   if ( !is.element("method",names(x)) ) #typically use default
                       x$method <- method
                  if ( !is.element("ana.obj",names(x)) ) #typically use default
                       x$ana.obj <- as.name("ana.obj")
                   eval(x)
               },
               scope=scope,method=method,
               covar=local.covar)
    names(rspec) <- term.rownames[ reg.specials ]
### used a common 'scope' and 'method' for all specials
    rspec.check <-
        sapply(rspec, function(x) any(x %in% c("(NA)","()")))
    if (any(rspec.check)){
        bad.terms <-
            paste(c("invalid term(s) in formula:",names(rspec)[rspec.check]),
                  collapse=" ")
        stop(bad.terms)
    }

### use all combinations of the expanded variables ?
    if (is.null(expand.specials))
        expand.specials <- 
            length(rspec)>1 && any(diff(range(sapply(rspec,length)))!=0)

    if (expand.specials)
        rspec <- do.call("expand.grid",rspec)

    if (using.R ){
        term.list <- c(rspec,as.list(c(reg.plain,plus="+",colon=":")))
    }
    else { #argh! Splus 3.4
        rspec <- lapply(rspec,as.character)
        term.list <- c(rspec,as.list(c(reg.plain,plus="+",colon=":")))
        names(term.list) <- c(names(rspec),names(reg.plain),"plus","colon")
    }
    
### order is <var,conj,var,conj,...,conj,var>
    spec.col.order <-
        if (length(terms.conjuncs)==0)
            term.names
        else
            c(term.names,terms.conjuncs)[c(rep(c(0,n.join+1),n.join)+
                                           rep(1:n.join,rep(2,n.join)),n.join+1)]
####prepend response
    spec.col.order <- c(response, spec.col.order) 
    text.formula <- do.call("paste",term.list[spec.col.order])

### now process it all
    res <- list()
    for (i in seq(along=text.formula)) {
        this.formula <- eval(parse(text=text.formula[i]))
        res[[i]] <-
            lapadj(this.formula, ana.obj,...)
        class(res[[i]])<-"bqtl"
        res[[i]]$call <- substitute(bqtl(this.formula,ana.obj,...))
    }
    if (length(res)==1) {
        res <- res[[1]]
    }
    else
        class(res) <- "bqtl.list"

    attr(res,"call") <- substitute(this.call)
    res
}
configs<-
    function(x,...,scope,method=NULL)
{
###
### configs will return a single variable name when x has length one
### and there are no additional arguments,
###  i.e. configs(x[1],scope=names(vars)) == list(names(vars)[x[1]])
###  additional args are concat'ed w "+" separators
### when dim(x)==NULL or dim(x)[[1]]==1 a vector is returned
###  i.e. configs(x,scope=names(vars)) == names(vars)[x]
### finally when dim(x)[[1]]>1 or dim(x$configs)[[1]]>1
### a vector is returned with elements paste(scope[x[,i]],collapse="+"
###
###  I parenthesize everything. S+3.4 needs this workaround 
###    ~(a):(u+v):(w) parses correctly, but not  ~a:(u+v):w

    if (missing(scope)) stop("missing scope arg - possible bqtl syntax error")
    if ( is.atomic(x) && (length(x) < 2) ) {
        res <- paste(scope[c(x,...)],collapse="+")
        return(paste("(",res,")",sep=""))
    }
    else {
        if (length(list(...))!=0)
            stop("only one arg allowed with vector or matrix")
    }
    
    if (is.element("configs",names(x)))
    {
        x <- x$configs
    }
    if (length( dm <- dim(x))>0)
    {
        dim(x) <- c(dm[1],prod(dm[-1]))
    }
    else
    {
        dim(x) <- c(1,length(x))
    }
    res <-    apply(x,2,function(x,dat) paste(dat[x],collapse="+"),dat=scope)
    paste("(",res,")",sep="")
}

"covar"<-
    function(x)
    x

"dom"<-
    function(x,...,scope=scope,method=method)
{
    if (missing(scope)) stop("missing scope arg - possible bqtl syntax error")
    new.scope <- scope[grep("^dom",scope)]
    if (length(new.scope)==0) stop("no 'dom' terms found")
    locus(x,...,scope=new.scope,method="fooey")
}
"%equiv%" <-
    function(x,y)
{
    if (!is.matrix(x)) dim(x) <- c(length(x),1)
    if (!is.matrix(y)) dim(y) <- c(length(y),1)
    dm <- c(ncol(x),ncol(y));
    nr <- nrow(x);
    if (nr!=nrow(y)) stop("nrow(x) != nrow(y)")
    t.x <- x[,rep(seq(dm[1]),dm[2]),drop=F]
    t.y <- y[,rep(seq(dm[2]),rep(dm[1],dm[2])),drop=F]
    matrix(apply(t.x==t.y,2,sum) == nr,nc=dm[2],dimnames=list(dimnames(x)[[2]],dimnames(y)[[2]]))
}
"formula.bqtl" <-
    function(object) {
        object$call[[2]]
    }
if (!exists("is.element")) {
    "is.element" <- function (el, set) match(el, set, 0) > 0
    "%in%"<-is.element
}
"lapadj"<-
function(reg.formula, ana.obj,  
         rparm = NULL,  tol = 9.900000000000001e-09,
	return.hess = F, mode.names = NULL, mode.mat = NULL,
          maxit = 100, nem = 1,setup.only=FALSE,subset=NULL,casewt=NULL,
         ...)
{
  this.call <- match.call(expand.dots=TRUE)
  n.data<-nrow(ana.obj$data)
  if (class(ana.obj)!="analysis.object") stop("ana.obj must be an analysis.object") 

  if (is.character(ana.obj$method))
    crstype <- if (is.element(ana.obj$method,c("BC1","RI.self","RI.sib"))) 1 else 2

  if(is.null(mode.mat))
      mode.mat <- ana.obj$mode.mat
  if(is.null(mode.names))
      mode.names <- dimnames(mode.mat)[[2]]
  
  iter <- maxit
  regr.names <- ana.obj$reg.names
  form.terms <- terms(reg.formula[-2],specials="covar")
  form.factors <- attr(form.terms,"factors")
  if (length(form.factors) == 0) { # intercept only, assume rparm[1]==0
      
      mf.expr <- expression(model.frame(reg.formula,ana.obj$data,
                        na.action=na.omit,subset=subset))[[1]]
      mf.expr$subset<- subset
      mf <- eval(mf.expr)
      y <-
        model.extract(   mf,                     "response")
      n  <- length(y)
      mean.y <- mean(y)
      var.y <- var(y)*(n-1)/n
      logpost <- sum(log(dnorm(y, mean.y , sqrt(var.y))))
      post <-
          exp((-(n - 1)/2*log(2*pi)-log(n)/2 -
               (n-1)/2*log(var.y*(n - 1)/2) +
               lgamma((n - 1)/2))-log(2))
      parm <- c(mean.y,log(var.y)/2)
      return(list(logpost=logpost,posterior=post,hk.exact=post,hk.approx=post,
                  parm=parm,reg.vec=character(0),call=this.call))
      ## return the null posterior
  }
  e.marker <- seq(nrow(form.factors)) %in% attr(form.terms,"specials")$covar
  actual.covar <- ifelse(e.marker, F,
                         !dimnames(form.factors)[[1]] %in% ana.obj$reg.names)
  covar.ind <- ( e.marker | actual.covar )
  covar.terms <- covar.ind%*%form.factors>0
  plain.terms <- (1-covar.ind)%*%form.factors>0
  both.terms <- covar.terms & plain.terms
  if (any(both.terms)) {
      covar.terms <- ( !both.terms ) & covar.terms
      plain.terms <- ( !both.terms ) & plain.terms
      both.index <- seq(along=both.terms)[both.terms]
      both.covar <- form.factors[covar.ind,both.terms] %equiv% form.factors[covar.ind,covar.terms]
      if (!apply(both.covar,1,any)) stop("can't find covariate main effect for interaction")
      both.plain <- form.factors[!covar.ind,both.terms] %equiv% form.factors[!covar.ind,plain.terms]
      if (!apply(both.plain,1,any)) stop("can't find locus main effect for interaction")
  }
  
  if (any(covar.ind)) {
      lhs <- paste(deparse(reg.formula[[2]]),"~")
      covar.formula <-
          eval(parse(text=paste(c(lhs,dimnames(form.factors)[[2]][covar.terms]),
                     collapse="+")))
      covar.terms <- terms(covar.formula)
      cfe <- expression(model.frame(covar.terms,ana.obj$data,
                                 na.action=na.omit,subset=subset))[[1]]
      cfe$subset <- subset
      covar.frame <- eval(cfe)
      subset <- row.names(ana.obj$data)%in%row.names(covar.frame)
      covar.matrix <- model.matrix(covar.terms,covar.frame)
      y <- model.extract(covar.frame,"response")
      covar.assign <-
          if (exists("is.R") && is.R() )
              attr(covar.matrix,"assign")[-1]
          else
              rep(seq(along=attr(covar.matrix,"assign"))-1,sapply(attr(covar.matrix,"assign"),length))[-1] 

      covar.reps <- table(factor(covar.assign,sort(unique(covar.assign))))
      if (any(both.terms)) {
          pttz <- apply(both.covar,1,function(x) seq(x)[x])-1
          pttx <- apply(both.plain,1,function(x) seq(x)[x])-1
          nt2uz <- length(pttx)
      }
      else {
          pttz <- pttx <- nt2uz <- 0
      }
      ncz <- ncol(covar.matrix)-1
      nz2uz <- ncz
      ptz <- seq(along=covar.assign)-1 
  }
  else {
      y <- eval(reg.formula[[2]],ana.obj$data)
      
      if (is.null(subset))
        subset <- !is.na(y)
      if (length(y)!=nrow(ana.obj$loc.right) )
        stop("length of dependent variable does not match other args")
      y <- y[subset]
      covar.matrix <- NULL
      ptz <- pttx <- pttz <- nt2uz <- nz2uz <- ncz <- 0
  }
  
  if (any(plain.terms)) {
      plain.formula <-
          eval(parse(text=paste(c("~",dimnames(form.factors)[[2]][plain.terms]),collapse="+")))
      reg.vec <- if (exists("is.R") && is.R() )
          as.character(attr(form.terms,"variables")[-1])[plain.terms]
      else
          as.character(attr(form.terms,"variables"))[plain.terms]
      
      which.vars <- sort(match(reg.vec, regr.names))
      
      if (length(which.vars)==0)
          stop("regressors not recognized - check spelling")
      ncx <- length(which.vars)
      cols.used <- col(regr.names)[which.vars]
      z.used <- unique(cols.used)
      nloc <- length(z.used)
      nrx <- if (crstype==1)
          2^nloc
      else
          3^nloc
      
      if (crstype==1) {
          if (nloc==1) reg.frame <- as.data.frame(mode.mat)
          else
              reg.frame <- as.data.frame(do.call("expand.grid", rep(list(mode.mat), nloc)))
          marker.distances <- switch(ana.obj$method,
                                     "RI.self"={ana.obj$map.frame$lambda/(2-ana.obj$map.frame$lambda)},
                                     "RI.sib"={ana.obj$map.frame$lambda/(4-3*ana.obj$map.frame$lambda)},
                                     ana.obj$map.frame$lambda)
      }
      else
      {
          marker.distances <- ana.obj$map.frame$lambda
          modes.used <- row(regr.names)[which.vars]
          col.seq <- match(cols.used,z.used)
          loc.indx <- as.matrix(do.call("expand.grid",rep(list(1:3),nloc)))
          reg.frame <-
              as.data.frame(matrix(mode.mat[cbind(c(loc.indx[,col.seq]),
                                                  rep(modes.used,rep(nrx,ncx)))],
                                   nc=ncx))
      }
      names(reg.frame) <- regr.names[which.vars]
      reg.proto <- model.matrix(plain.formula,reg.frame)[,-1,drop=F]
      ncx <- ncol(reg.proto)
      ptx <- seq(from=0,by=1,length=ncx)
  } # if (any(plain.terms)) {
  else
  {
      nrx <- 1
      ncx <- nloc <- reg.proto <- ptx <-  z.used <- modes.used <- 0
  }

### set up nparm
  nx2uz <- ncx
  nreg <- 1 + nx2uz + nz2uz + nt2uz
  np <- nreg+1
  
  n <- length(y)
  nparm <- c(n,nrx,ncx,nloc,ncz,
             nx2uz,nz2uz,nt2uz,
             nreg,np,ptx,ptz,pttx,pttz)
  reg.vec <- c(dimnames(reg.proto)[[2]],dimnames(covar.matrix)[[2]][-1],dimnames(form.factors)[[2]][both.terms])
  
### rparm for intercept is assumed to be in the first element, named "intercept",
###  and names(rparm) must be matched, except when length matches number of variables
###  in regr.names
### interactions can ONLY be handled by scalar rparm unless
### names(rparm)%in%dimnames(reg.proto)[[2]]
###
  rparm.named <- !is.null(names(rparm))
  
  if(is.null(rparm)) {
      rparm <- rep(0, nreg)
  }
  else {
      if (rparm.named)
          rparm <- c(rparm["intercept"],rparm[reg.vec])
      else
      {
          len.choices <- c(1,length(regr.names)+1,nreg,2)
          
          rparm <-
              switch(match(length(rparm),len.choices,5),
                     c(0,rep(rparm,nreg-1)),
                 {
                     names(rparm) <- c("intercept",regr.names)
                     c(rparm["intercept"],rparm[dimnames(reg.proto)[[2]]])
                 },
                     rparm,
                     if (method=="F2") c(0,rparm[modes.used]) else NULL,
                     stop("unamed 'rparm' has wrong length")
                       )
      }
  }
  if (rparm.named&&is.na(rparm["intercept"]))
      rparm["intercept"] <- 0.0 # provide unstated intercept
  if (any(is.na(rparm)))
      stop(paste("rparm undefined for term:",reg.vec[is.na(rparm)][1]))
  if(nloc > 1) {
      need.prod <- as.numeric(ana.obj$loc.right[subset, z.used[ - nloc], drop = F] > 
                              rep(z.used[-1], rep(n, nloc - 1)))
      lambda <- rep(0, nloc - 1) 
      for(i in 2:nloc)
          lambda[i - 1] <-
              prod(marker.distances[z.used[i - 1]:(z.used[i] - 1)])
  }
  else need.prod <- lambda <- 0
  if (is.null(ana.obj$casewt)) {
      casewt <- -1.0 #flag for no weights
  }
  else
  {
      if (length(ana.obj$casewt)!=length(subset))
          stop("length(casewt) must  == 1 or nrow(ana.obj$data)")
      else
          casewt <- ana.obj$casewt[subset]
  }
  if (setup.only) {
      res <-
          list("lapadj",
               csrtype=as.integer(crstype),
               nparm = as.integer(nparm),
               y = as.double(y),
               x = as.double(reg.proto),
               z = as.double(if (ncz>0) covar.matrix[,-1] else 0),
               rpar = as.double(rparm),
               tab =
               if (nrx>1) {
                   as.double(ana.obj$state.matrix[subset, z.used,  ])
               }
               else {
                   as.double(rep(1,n))
               } ,
               needProd = as.integer(need.prod),
               lambda = as.double(lambda),
               coefs = double(nreg),
               casewt=as.double(casewt),
               ln.sigma = double(1),
               hkExact = double(1),
               hkApprox = double(1),
               llk = double(1),
               hess = double((np)^2),
               postApprox = double(1),
               iter = as.integer(iter),
               tol = as.double(tol),
               nem = as.integer(nem) )
      attr(res,"subset") <- subset
      attr(res,"reg.names") <- dimnames(reg.proto)[[2]]
      attr(res,"N")<-c(N=n.data,N.omit=n.data-n,N.used=n)
      res
  }
  else
  {
      z <- .C("lapadj",
              csrtype=as.integer(crstype),
              nparm = as.integer(nparm),
              y = as.double(y),
              x = as.double(reg.proto),
              z = as.double(if (ncz>0) covar.matrix[,-1] else 0),
              rpar = as.double(rparm),
              tab =
              if (nrx>1) {
                  as.double(ana.obj$state.matrix[subset, z.used,  ])
              }
              else {
                  as.double(rep(1,n))
              } ,
              needProd = as.integer(need.prod),
              lambda = as.double(lambda),
              coefs = double(nreg),
              casewt=as.double(casewt),
              ln.sigma = double(1),
              hkExact = double(1),
              hkApprox = double(1),
              llk = double(1),
              hess = double((np)^2),
              postApprox = double(1),
              iter = as.integer(iter),
              tol = as.double(tol),
              nem = as.integer(nem))
      if ( (z$postApprox==0) && (z$hkApprox==0))
        adj <- 1.0
      else
        adj <- z$postApprox/z$hkApprox
      list(adj = adj, logpost = z$llk, parm =
           c(z$coefs[1: (nreg)], z$ln.sigma), posterior = z$postApprox,
           hk.approx = z$ hkApprox, hk.exact = z$hkExact, reg.vec = reg.vec,
           rparm = rparm, hess = if(return.hess) matrix(z$hess, nc = np) else NULL,
           iter = z$iter,N=c(N=n.data,N.omit=n.data-n,N.used=n),
           call=this.call)
  }
} 
"linear.bayes" <-
    function(x, ana.obj, partial=NULL,rparm,specs,
             scope,subset,casewt,...){
        x.call <- match.call(expand.dots=TRUE)
        if (class(x)!="varcov"){
            vc.call <- x.call
            vc.call[[1]] <- as.name("varcov")
            OK.args <-
                is.element(names(vc.call),
                           c("subset","casewt","ana.obj","partial",
                             "scope","x"))
            OK.args[1] <- TRUE
            vc.call[!OK.args] <- NULL
            x <- eval(vc.call)
        }
        else {
            vc.call <- x$call
        }
        default.specs <-
            list(
                 gene.number=c(1,2,3,4,5),
                 burn.in=1,
                 n.cycles=c(0,200,200,100,100))
        
        if (missing(specs)){              
            specs <- default.specs
        }
        else {
            named <- pmatch( names(specs), names(default.specs), 0)
            if (any(named == 0)) stop("illegal name in specs arg")
            names(specs) <- names(default.specs)[named]
            specs <- c(specs, default.specs[-named] )
        }
        do.hk <- any(specs$gene.number<3)
        if (missing(rparm)) {
            warning("default: rparm  <- 0  may be unstable")
            rparm <- 0
        }
        if (do.hk)
            hk <- twohk(x,ana.obj,rparm=rparm,...)
        else
            hk <- NULL
        nums.2.sample <- sort(specs$gene.number[specs$n.cycles!=0])
        contig.nums <- all(diff(nums.2.sample)==1)
        if (!contig.nums) stop("specs$gene.number must be contiguous")
        var.indx <-
            if (ana.obj$method=="F2")
                seq(from=1,by=2,length=ncol(x$var.x)/2)
            else
                seq(ncol(x$var.x))
        if (specs$burn.in > 0){
            starts <- vector("list",max(nums.2.sample))
            for (i in nums.2.sample) {
                tmp <-  swap(x,sample(var.indx,i),
                             rparm,1,ana.obj=ana.obj,...)$configs[,i,1]
                starts[[i]] <- tmp[tmp != 0]
            }
            
        }
        swaps <- vector("list",max(nums.2.sample))
        for (i in nums.2.sample){
            j <- match(i,specs$gene.number)
            cycles <- specs$n.cycles[[j]]
            if (cycles != 0)
                swaps[[i]] <-
                    swap(x,starts[[i]], rparm,cycles, ana.obj=ana.obj,...)
        }
        swaps.OK <- sapply(swaps,function(x) 
                           if (length(unlist(x))==0)
                           TRUE          # NULL is OK
                           else {
                               call.indx <- match("call",names(x),0)
                               !any(is.na(unlist( x[-call.indx] )))})
        smry <- vector("list",length(swaps))
        
        for (i in seq(along=swaps))
            if (swaps.OK[[i]]) smry[[i]] <- summary(swaps[[i]])

        if (any(!swaps.OK)){
            warning("attempt to summarize samples failed. rparm = 0 ??")
            odds <- loc.posterior <- coefs <- NULL          
        }
        else{
            swap.odds <- sapply(smry,function(x) 1/x$ratio$mean)
            if (do.hk){
                odds <-
                    c(sum(hk$loc.1),
                      cumprod(c(sum(hk$loc.2),
                                swap.odds[nums.2.sample[nums.2.sample>2]])))
                post.pr <- odds/sum(odds)
                if (ana.obj$method=="F2") {
                    loc.posterior <-
                        c(hk$loc.1%*%rep( 1/sum(hk$loc.1), 3)) * post.pr[1] +
                            2*c(hk$loc.2%*%rep( 1/sum(hk$loc.2), 3)) *
                                post.pr[2] +
                                    sapply(smry[-(1:2)],
                                           "[[", "loc.posterior" ) %*%
                                               post.pr[-(1:2)]
                }
                else {
                    loc.posterior <-
                        hk$loc.1 / sum(hk$loc.1) * post.pr[1] +
                            2 * hk$loc.2 / sum(hk$loc.2) * post.pr[2] +
                                sapply(smry[-(1:2)],
                                       "[[", "loc.posterior" ) %*%
                                           post.pr[-(1:2)]
                }
                coefs <-
                    cbind(c(hk$coefs.1),c(hk$coefs.2),
                          sapply(smry[-(1:2)],"[[","coefs")) %*% post.pr
            }
            else {
                swap.nums <- nums.2.sample[nums.2.sample>2]
                odds <-
                    cumprod(c(1,swap.odds[ swap.nums[-1] ]))
                post.pr <- odds/sum(odds)
                loc.posterior <-
                    sapply(smry[swap.nums],"[[", "loc.posterior" ) %*% post.pr
                coefs <- sapply(smry[swap.nums],"[[","coefs") %*% post.pr
            }
            
        }
        res <- list(hk=hk, swaps=swaps, smry=smry, odds=odds,
                    coefficients=coefs, loc.posterior=loc.posterior,
                    specs=specs,varcov.call=vc.call,call=x.call)
        class(res) <- "linear.bayes"
        res
    }

locus<-
    function(x,...,scope=scope,method=NULL,chromo=NULL,cM=NULL,ana.obj=NULL)
{
###
###   
### When method != "F2", locus will return a single variable name when x
### has length one, there are no additional arguments. When method ==
### "F2", then there will be two values returned if x has length one,
### hopefully the additive and dominance terms of that locus.
### 
### additional args are concat-ed w "+" separators when dim(x)==NULL or
### dim(x)[[1]]==1 a vector is returned.
### 
### I parenthesize everything. S+3.4 needs this workaround 
###  ~(a):(u+v):(w) parses correctly, but not  ~a:(u+v):w
###
  if (missing(scope)) stop("missing scope arg - possible bqtl syntax error")
  dots <- list(...)
  if (is.null(chromo)&& length(names(dots))!=0 &&
      any(cloc <- which(1==pmatch(names(dots),"chromo",0)))){
    chromo <- dots[[cloc]]
    if (length(dots)>1)
      stop("... not allowed")
    else
      dots<-NULL
  }
  if (!is.null(chromo)){
    if (!missing(x)) stop("using both x and chromo args not allowed")
    x <- if (is.null(cM))
        map.index(ana.obj,chromo=chromo)
    else
        map.index(ana.obj,chromo=chromo,cM=cM)
  }

  x.call <- match.call()$x
  if  ( length(x.call) ==1 && deparse(x.call) == "all" ){  # all loci ? 
    x <- seq(along=scope)
    if (method=="F2") dim(x) <- c(2,length(x)%/%2)
    if (length(dots) != 0) stop("... not allowed")
    return(configs(x,scope=scope))
  } #else
  
  if (length(x)>1 && length(dots)!=0)
    stop("only one arg allowed with vector or matrix")
  if (method == "F2") {
    if (length(x)==1) {
      x.1 <- 2*x-1
      y <- 2*x
      if ( length(dots)==0 )
        dots <- list(y)
      else
        dots <- c(list(y),lapply(dots,function(x) c(2*x-1,2*x)))
      configs.args <-  c(list(x=x.1),dots,scope=list(scope))
    }
        else {
            x.1 <- 2*x-1
            y <- 2*x
            if ( length( dm <- dim(x) )>0 && dm[1]>1 ) {
                x.1 <- aperm(array(c(x.1,y),c(dm[1],prod(dm[-1]),2)),c(3,1,2))
                dim(x.1) <- c(dm[1]*2,prod(dm[-1]))
            }
            else {
                x.1 <- rbind(x.1,y)
            }
            configs.args <-  c(list(x.1),scope=list(scope))
        }
    do.call("configs",configs.args)
  }
  else {
    configs.args <- c(list(x),dots,scope=list(scope))
    do.call("configs",configs.args)
  }
  
}

"make.analysis.obj"<-
    function(data, map.frame, marker.frame, marker.levels=NULL,method="F2",
             casewt=NULL,varcov=FALSE,mode.mat=NULL)
{
    
### map.frame - a map.frame.object OR a
### lambda - vector of pseudo-marker distance parms - has names attribute
### marker.frame has dimnames attribute to match 
###
### methods are: F2, BC1, RI.self, and RI.sib
    "check.markers" <-
        function(x,y)
        {
            tab <- table(unlist(sapply(x,as.character)))
            bad.x <- !(names(tab) %in% c(y,"NA"))
            bad.y<- !(y %in%names(tab))
            list(bad.markers=tab[bad.x],unused.levels=y[bad.y])
        }
    
    this.call <- match.call()
    
    if( inherits(data,"map.frame") )
        warning("first arg should not be a map.frame! ?? ")
    
    if( inherits(marker.frame,"map.frame") )
        warning("third arg should not be a map.frame! ?? ")
    
    method.types <- c("F2","BC1","RI.self","RI.sib")
    method.code <- pmatch(method,method.types)
    if (is.na(method.code))
        stop(paste("method =",method, "not found"))
    method <- method.types[method.code]
    
    if (length(dm.data <- dim(data))==0) dm.data<- c(length(data),1)
    if (dm.data[1]!=nrow(marker.frame)) stop(
               "nrow(data)!=nrow(marker.frame)")
    
    
### RI strains use coding of 1,2,6
    
    if (length(marker.levels)==0) marker.levels <-
        switch(method,
               F2=f2.levels(),
               BC1=bc1.levels(),
               RI.self=,
               RI.sib=ri.levels() 
               )
### set mode.mat
    if (is.null(mode.mat)){
        mode.mat <-
            if (method=="F2")
                matrix(c(-1,0,1,-1,1,-1),nc=2,
                       dimnames=list(marker.levels[1:3],c("add","dom")))
            else
                matrix(c(-1,1),nc=1,dimnames=list(marker.levels[1:2],NULL))
    }
    else {   # check the supplied version
        dmmt <- dim(mode.mat)
        if (method=="F2") {
            if ( length(dmmt)==0 || dmmt[1]!=3 || dmmt[2] > 2)
                stop("F2 mode.mat must be 3 x 2 matrix")
            if (dmmt[2]==1)
                warning("F2 mode.mat has just one column")
            else
                if (length(dimnames(mode.mat)[[2]])==0)
                    stop("need column names for F2 mode.mat")
        }
        else {
            if ( length(dmmt)==0 || dmmt[1]!=2 || dmmt[2] > 1)
                stop("non-F2 mode.mat must be 2 x 1 matrix")
        }
    }
    
    if(inherits(map.frame,"map.frame")) {
        lambda <- map.frame$lambda
        marker.names <-
            as.character(map.frame$marker.name[is.element(map.frame$is.marker,
                                                          c("TRUE","T"))])
        names.lambda <- as.character(map.frame$marker.name)
        m.chromo <- map.frame$chr.num
    }
    else {
        ## not actually a map.frame
        names.lambda <- names(map.frame)
        if (length(map.frame) != length(names.lambda))
            stop("problem with map.frame arg - must be data frame or vector with names()")
        lambda <- map.frame
        marker.names <- dimnames(marker.frame)[[2]]
        m.chromo <- cumsum(c(1,lambda==0))[-length(lambda)-1]
    }
    
### check the setup
### markers found to match marker.names in map.frame
    if (any(lost.markers <-
            !is.element(marker.names,dimnames(marker.frame)[[2]])
            )){
        cat(deparse(this.call$map.frame),"defines markers NOT found in",
            deparse(this.call$marker.frame),":\n")
        print(marker.names[lost.markers])
        stop("bailing out")
    }
    
### marker values match marker.levels??

    ck.mark <- check.markers(marker.frame[,marker.names],marker.levels)
    if (length(ck.mark$bad) != 0)
    {
        cat("Marker values found, but not defined in marker.levels:\n")
        print(ck.mark$bad)
        stop("bailing out")
    }
    if (any(marker.levels[1:2]%in%ck.mark$unused.levels))
        warning(paste("unused marker.levels:",
                      paste(ck.mark$unused.levels,collapse=" ")))
    if (method=="F2" && all(marker.levels[3:5]%in%ck.mark$unused.levels) )
        warning(paste("method=F2, but did not find",
                      paste(marker.levels[3:5],collapse=", "),
                      "in",deparse(this.call$marker.frame),"\n"))
    m.state.matrix <-
        array(0,c(dm.data[1],length(lambda),if (method=="F2") 3 else 2),dimnames=list(NULL,names.lambda,NULL))
    m.loc.right <- array(0,c(dm.data[1],length(lambda)),dimnames=list(NULL,names.lambda))
    
    
### the principal reason for the explicit loop is to avoid large objects in R
###  --- the space required sometimes kills the whole function
    for (i in unique(m.chromo)) {
        which.chromo <- m.chromo==i
        m.frame <- as.data.frame(matrix(NA, nr = dm.data[1], nc = sum(which.chromo)))
        names(m.frame) <- names.lambda[m.chromo==i]
        names.to.use <-
            names(m.frame)[is.element(names(m.frame),marker.names)]
        
        if (length(names.to.use) == 0) stop(
                  "names(lambda) did not match names(marker.frame)")
        for (i.name in names.to.use) m.frame[,i.name] <- marker.frame[,i.name]
        
        lambda.i <- lambda[which.chromo][ - sum(which.chromo)]
        m.state.matrix[,which.chromo,] <-
            make.state.matrix( make.marker.numeric(m.frame,marker.levels),
                              lambda.i , method)
        left.loc <- min(seq(which.chromo)[which.chromo])-1
        m.loc.right[,which.chromo] <- left.loc + make.loc.right(m.frame, lambda.i)
    }
    
    m.regressor.matrix <- make.regressor.matrix(m.state.matrix,mode.mat)
    
    m.varcov <-
        if (varcov!=FALSE)
            make.varcov(m.regressor.matrix, data[,varcov],casewt=casewt)
        else
            NULL
    reg.names <- dimnames(m.regressor.matrix)[[2]]
### sanity check on m.regressor.matrix
    if (any(null.var <-
            apply(apply(m.regressor.matrix,2,range),2,diff) <
            0.5/nrow(m.regressor.matrix)))
        warning(paste("null or tiny variance in",
                      paste(reg.names[null.var],collapse=" "),"\n"))
    
    dim(reg.names) <- c(dim(m.state.matrix)[[3]]-1,dim(m.state.matrix)[[2]])
    dimnames(reg.names) <-
        list(
             if (method!="F2") NULL else dimnames(mode.mat)[[2]],
             names.lambda
             )
    version.stamp <- if (exists("is.R")&&is.R()) version.bqtl else "splus"
    structure(
              list(data=data.frame(data,m.regressor.matrix), varcov = m.varcov, 
                   reg.names=reg.names,method=method,
                   state.matrix = m.state.matrix, loc.right = m.loc.right,
                   map.frame=map.frame,casewt=casewt,mode.mat=mode.mat,
                   call = this.call,version=version.stamp),
              class="analysis.object"
              )
}
"make.loc.right"<-
    function(marker.frame, marker.distances)
{
    ncm <- ncol(marker.frame)
    nrm <- nrow(marker.frame)
    res <- matrix(0, nr = nrm, nc = ncm)
    dimnames(res) <- list(NULL, names(marker.frame))
    init <- rep(ncm, nrm)
    for(i in rev(1:ncm)) {
        is.terminal <-
            if(i < ncm && marker.distances[i] == 0)
                rep(TRUE, nrm)
            else
                !is.na(marker.frame[, i])
        init[is.terminal] <- i
        res[, i] <- init
    }
    res
}
"make.location.prior"<-
    function(x, add.2.end = 0, normalize = T)
{
###x is a lambda parametter vector
    chromo.break <- x == 0
    map.dist <- map.dx(x)
    len.map <- length(map.dist)
    res <-
        ifelse(chromo.break,
               add.2.end/2 + c(0, map.dist[ - len.map]), 
               map.dist/2) + ifelse(c(TRUE, chromo.break[ - len.map]), 
                                    map.dist,
                                    add.2.end/2 + c(0, map.dist[ - len.map])/2)
    if(normalize)
        res/sum(res)
    else res
}
"make.map.frame"<-
    function(dx,chr.num=NULL,
             prior=make.location.prior(lambda),morgan=100,nint=NULL,reso=NULL)
{
    if (is.null(nint)&&!is.null(reso))
        nint <- marker.fill(make.map.frame(dx),reso, TRUE )
    int.suffix <- function(nint)
    {
        res <- sapply(nint, function(x)
                  {
                      strng <- if(x > 1) {
                          c("", format((1:(x - 1))/10^(trunc(log10(x - 1) +
                                                             1))))
                      }
                      else {
                          ""
                      }
                      substring(strng, 2)
                  }
                      )
        unlist(res)
    }
    if (is.data.frame(dx)) {
        ## allow a data.frame with the first column containing the distances
        if (ncol(dx)>1) {
            dx.synonyms <- c("cM","cm","dx","M","m")
            dx.name <- dx.synonyms[is.element(dx.synonyms,names(dx))][1]
            if ( length(dx.name)!=0 && dx.name!="NA") {# found explicitly labelled distance
                marker.names <-
                    if (any(pmatch(names(dx),"marker.names",0)!=0))
                        dx$marker.name
                    else {
                        warning("could not find marker.names using row.names(dx) instead")
                        row.names(dx)
                    }
                if (is.element(dx.name,c("m","M"))){
                    morgan <- 100
                    mgd <- dx[,dx.name]/100
                }
                else
                    mgd <- dx[,dx.name]
                if ( missing(chr.num) && ("chr.num"%in%names(dx)) )
                    chr.num <- dx$chr.num
                
            }
            else {
                warning(paste("using first column",names(dx)[[1]], "as map distance"))
                marker.names <- row.names(dx)
                mgd <- dx[[1]]
            }
        }
        else { # 1 column data.frame
            marker.names <- row.names(dx)
            mgd <- dx[[1]]
        }
    }
    else { #not a  data.frame
        mgd <- dx 
        if (is.null(names(dx)))
            names(dx) <- paste("marker",seq(along=chr.num),sep=".")
        marker.names <-  names(dx)
    }
    chr.num <-
        if (is.null(chr.num))
            cumsum( diff(c(Inf,mgd))<0 )
        else
            eval(chr.num)
### check ordering within chromos
    mk.order <- order(chr.num,mgd)
    if (any(out.order<-
            tapply(mgd,chr.num,function(x) length(x)>1 && any(diff(x)<0))))
        stop(paste("markers not in increasing order on chr:",
                   paste(names(out.order)[out.order],collapse=",")))
    
    ##  check arg lengths
    if (!all.equal( length(nint), length(mgd), length(chr.num) ) )
        stop("lengths of nint, chr.num, dx not all equal")
    ## default is 'right' if only one marker on a chromo
    pos.type <- as.numeric(c(0,chr.num[-length(chr.num)])==chr.num) +
        as.numeric(c(chr.num[-1],0)==chr.num)*2
    pos.type <- factor(c(1,1,2,3)[pos.type+1],1:3,c("right","left","center"))
    is.marker <- rep(TRUE,length(chr.num))
    
    increments <- diff(mgd/morgan)

    not.r <- ifelse(pos.type=="right",0.0,
                    c(exp(-increments),0.0))
    if (!is.null(nint)) {
        orig <- cumsum(c(1,nint[-length(nint)]))
        not.r <- rep(not.r^(1/nint),nint)
        marker.names <- paste(rep(marker.names,nint),int.suffix(nint),sep="")
        is.marker <- rep(is.marker,nint)
        is.marker[-orig] <- FALSE
        chr.num <- rep(chr.num,nint)
        mgd <-
            rep(mgd,nint) +
                morgan*rep(c(increments,0)/nint,nint)*
                    (unlist(sapply(nint,seq))-1)
        pos.type <- rep(pos.type,nint)
        pos.type[-orig] <- "center"
    }
    lambda <- ifelse(not.r==0,0,2*not.r-1)
    prior <- eval(prior)
    pos.plot <- mgd
    locus<-paste("C",chr.num,as.numeric(format(100*mgd/morgan,digits=4)),
                 sep=".")
    res <- data.frame(
                      marker.name=
                      I(as.character(make.names(marker.names,TRUE))),
                      cM=as.numeric(mgd),
                      prior=as.numeric(prior),
                      pos.type=I(pos.type),
                      is.marker=I(is.marker),
                      pos.plot=as.numeric(pos.plot),
                      lambda=as.numeric(lambda),
                      locus=I(locus),
                      chr.num=I(chr.num))
    attr(res,"morgan") <- morgan
    class(res)<- c("map.frame",class(res))
    res
}
"make.marker.numeric" <- function(marker.frame,
                                  level.names=NULL )
{
### convert factor or character frames into numeric
### for use in make.state.matrix
  if (is.null(level.names)) {
    level.names <- c("AA","Aa","aa","A-","a-","--")}
  else
    if ( length(level.names) != 6 )
      stop("length(level.names) must = 6 if none NULL value is used")

  char.frame <- as.matrix(marker.frame)
  char.frame[is.na(char.frame)] <- level.names[6]
  char.frame[char.frame=="NA"] <- level.names[6]
  level.numeric <- seq(length(level.names))
  names(level.numeric) <- level.names
  dm <- dim(marker.frame)

  res <- matrix(level.numeric[char.frame],nc=dm[2],
                dimnames=list(NULL,names(marker.frame)))
  if (any(is.na(res))) {
    stop(paste("unrecognized genotype(s):",unique(res[is.na(res)])))
  }
  res
  
}
"make.regressor.matrix"<-
    function(state.matrix, mode.mat = NULL)
{
    dsm <- dim(state.matrix)
    dsn <- dimnames(state.matrix)[[2]]
    dimnames(state.matrix) <- NULL
    if(is.null(mode.mat))
        if(dsm[3] == 3) {
            ## assume F2 with both additive and dominance effects
            mode.mat <- cbind(add = c(-1, 0, 1), dom = c(-1, 1, -1))
        }
        else {
            ## assume BC1
            mode.mat <- matrix(c(-1, 1), nr = 2)
        }
    if(nrow(mode.mat) == 1) {
        dim(state.matrix) <- c(dsm[1] * dsm[2], dsm[3])
        regressor.matrix <- state.matrix %*% mode.mat
        dim(regressor.matrix) <- dsm[1:2]
        dimnames(regressor.matrix) <- list(NULL, dsn)
    }
    else {
        dim(state.matrix) <- c(dsm[1] * dsm[2], dsm[3])
        regressor.matrix <- state.matrix %*% mode.mat
        dimnames(regressor.matrix) <- NULL
        dim(regressor.matrix) <- c(dsm[1:2], ncol(mode.mat))
        regressor.matrix <- aperm(regressor.matrix, c(1, 3, 2))
        mode.names <-  dimnames(mode.mat)[[2]]
        dim(regressor.matrix) <- c(dsm[1], dsm[2] * ncol(mode.mat))
        dimnames(regressor.matrix) <-
            list(NULL,
                 if (is.null(mode.names)) dsn else
                 paste(rep(mode.names, 
                           dsm[2]), rep(dsn, rep(ncol(mode.mat), dsm[2])),
                       sep = "."))
    }
    regressor.matrix
}

"make.state.matrix"<-
    function(marker.frame, marker.distances, method = "F2")
{
###it is assumed
### F2 method has as.numeric(marker.states)%in%c( 1,2,3,4,5,6 )
### for BC1 method as.numeric(marker.states)%in%c( 1,2,5,6 )
### or NA's
###  if NA's are present, they are recoded to 6
  
                                        # no this is wrong NA cause error

###  the assumed setup is as follows:
###    strains are A and a
###
###  type   F2.code   BC.code RI.code
###  AA      1         1       1
###  Aa      2         2 
###  aa      3                 2 
###  A-      4           
###  a-      5           
###  --      6         6       6
###
### marker. frame is assumed to be all numeric
###
    sum.or.one <- function(x)  {y <- sum(x);if (y>0) y else 1.0}
    md.unique <- unique(marker.distances)
    md.ind <- factor(marker.distances, md.unique)
    n.ind <- length(md.unique)
###
### RI strains use an approximation - they aren't quite a Markov Chain
    
    if (method=="RI.self") md.unique <- md.unique/(2-md.unique)
    if (method=="RI.sib") md.unique <- md.unique/(4-3*md.unique)
    
    
    mf.num <- as.matrix(marker.frame)
    if (any(is.na(mf.num)))
        stop("marker.frame must be all numeric")

    if (n.ind==0){            # only one marker
      if (method=="F2")
        return(
               array(rbind(diag(3),c(1,2,0)/3,c(0,2,1)/3,c(1,2,1)/4)[mf.num,],
                     c(dim(mf.num),3)))
      else
        return(array(rbind(diag(2),NA,NA,NA,c(1,1)/2)[mf.num,],
                     c(dim(mf.num),2)))
    }

    switch(method,
           F2 = {
               tmat <- matrix(c(1, 1, 1, -1, 0, 1, 1, -1, 1), 3, byrow
                              = F)
               tinv <- matrix(c(1, -2, 1, 2, 0, -2, 1, 2, 1)/4, 3, 
                              byrow = F)
               t1 <- tmat[, 1, drop = F] %*% tinv[1,  , drop = F]
               t2 <- tmat[, 2, drop = F] %*% tinv[2,  , drop = F]
               t3 <- tmat[, 3, drop = F] %*% tinv[3,  , drop = F]
               tx.mat <- array(rep(t1, n.ind), c(3, 3, n.ind)) +
                   outer(t2, md.unique) + outer(t3, md.unique^2)
               unc.pr <- t1[1,  ]
               dom.A.pr <- c(1,2,0)/3
               dom.a.pr <- c(0,2,1)/3
               lp <- 3
           }
           ,
           RI.self=,
           RI.sib=,
           BC1 = {
               t1 <- matrix(1/2, nc = 2, nr = 2)
               t2 <- matrix(c(1, -1, -1, 1)/2, nc = 2)
               tx.mat <- array(rep(t1, n.ind), c(2, 2, n.ind)) +
                   outer(t2, md.unique)
               unc.pr <- t1[1,  ]                        
               lp <- 2
           }
           ,
           stop(paste(method, " method unknown")))
    left.marker <-
        matrix(c(diag(lp),rep(NA, if (lp==2) 8 else 9)),nr=lp)[,mf.num]
    dmm <- c(lp, dim(marker.frame))
    dim(left.marker) <- dmm
    center.marker <- right.marker <- left.marker
    dimnames(center.marker) <- list(NULL, NULL, dimnames(marker.frame)[[2]])
    zero.dx <- marker.distances == 0
    chr.number <- sum(zero.dx)
    is.miss <- mf.num[,c(TRUE, zero.dx)]==6
    is.miss <- rep(is.miss,rep(lp,length(is.miss)))
    left.marker[,  , c(TRUE, zero.dx)][is.miss] <- unc.pr
    is.miss <- mf.num[,c(zero.dx, TRUE)]==6
    is.miss <- rep(is.miss,rep(lp,length(is.miss)))
    right.marker[,  , c(zero.dx, TRUE)][is.miss] <- unc.pr
    if (method=="F2") {

        ## initialize first marker if necessary
        is.A.minus <- mf.num[,c(TRUE, zero.dx)]==4
        is.A.minus <- rep(is.A.minus,rep(3,length(is.A.minus)))
        is.a.minus <- mf.num[,c(TRUE, zero.dx)]==5
        is.a.minus <- rep(is.a.minus,rep(3,length(is.a.minus)))
        left.marker[,  , c(TRUE, zero.dx)][is.A.minus] <- dom.A.pr
        left.marker[,  , c(TRUE, zero.dx)][is.a.minus] <- dom.a.pr
        ## initialize last marker if necessary
        is.A.minus <-mf.num[,c(zero.dx, TRUE)]==4
        is.A.minus <- rep(is.A.minus,rep(3,length(is.A.minus)))
        is.a.minus <-mf.num[,c(zero.dx, TRUE)]==5
        is.a.minus <- rep(is.a.minus,rep(3,length(is.a.minus)))
        right.marker[,  , c(zero.dx, TRUE)][is.A.minus] <- dom.A.pr
        right.marker[,  , c(zero.dx, TRUE)][is.a.minus] <- dom.a.pr
    }
    for(i in 2:dmm[3]) {
        i.tx <- as.numeric(md.ind[i - 1])
        is.miss <- mf.num[,i]==6
        if(any(is.miss)) {
            should.be <- t(as.matrix(left.marker[, is.miss, i - 1])
                           ) %*% tx.mat[,  , i.tx]
            left.marker[, is.miss, i] <- t(should.be)
        }
        if (method == "F2") {
            is.A.minus <- mf.num[,i]==4
            is.a.minus <- mf.num[,i]==5
            if(any(is.A.minus)) {
                
                should.be <- t(t(as.matrix(left.marker[, is.A.minus, i - 1])
                                 ) %*% (tx.mat[,  , i.tx]*
                                        rep(dom.A.pr/unc.pr,rep(3,3))))

                left.marker[, is.A.minus, i ] <-
                    sweep(should.be, 2, apply(should.be, 2, sum.or.one), "/")
            }
            if (any(is.a.minus)) {
                should.be <- t(t(as.matrix(left.marker[, is.a.minus, i - 1])
                                 ) %*% (tx.mat[,  , i.tx]*
                                        rep(dom.a.pr/unc.pr,rep(3,3)))/unc.pr)

                left.marker[, is.a.minus, i ] <-
                    sweep(should.be, 2, apply(should.be, 2, sum.or.one), "/")
            }
        }
    }
    for(i in rev(1:(dmm[3] - 1))) {
        i.tx <- as.numeric(md.ind[i])
        is.miss <- mf.num[,i]==6
        if(any(is.miss)) {
            should.be <- t(as.matrix(right.marker[, is.miss, i + 1]
                                     )) %*% tx.mat[,  , i.tx]
            right.marker[, is.miss, i] <- t(should.be)
        }
        if (method == "F2") {
            is.A.minus <- mf.num[,i]==4
            is.a.minus <- mf.num[,i]==5
            if(any(is.A.minus)) {
                
                should.be <- t(t(as.matrix(right.marker[, is.A.minus, i + 1])
                                 ) %*% (tx.mat[,  , i.tx]*
                                        rep(dom.A.pr/unc.pr,rep(3,3))))

                right.marker[, is.A.minus, i ] <-
                    sweep(should.be, 2, apply(should.be, 2, sum.or.one), "/")
            }
            if (any(is.a.minus)) {
                should.be <- t(t(as.matrix(right.marker[, is.a.minus, i + 1])
                                 ) %*% (tx.mat[,  , i.tx]*
                                        rep(dom.a.pr/unc.pr,rep(3,3))))

                right.marker[, is.a.minus, i ] <-
                    sweep(should.be, 2, apply(should.be, 2, sum.or.one), "/")
            }
        }
    }
    for(i in 1:dmm[3]) {
        is.miss <- mf.num[,i]==6
        if(any(is.miss)) {
            should.be <- as.matrix((left.marker[, is.miss, i] * 
                                    right.marker[, is.miss, i])/unc.pr)
            
            center.marker[, is.miss, i] <-
                sweep(should.be, 2, 
                      apply(should.be, 2, sum.or.one), "/")
        }
        if (method == "F2") {
            is.A.minus <- mf.num[,i]==4
            is.a.minus <- mf.num[,i]==5
            if(any(is.A.minus)) {
                mpy <- ifelse(dom.A.pr==0.0,0.0,(1/dom.A.pr/unc.pr^2))
                should.be <- as.matrix((left.marker[, is.A.minus, i] * 
                                        right.marker[, is.A.minus, i])*mpy)
                

                center.marker[, is.A.minus, i ] <-
                    sweep(should.be, 2, apply(should.be, 2, sum.or.one), "/")
            }
            if (any(is.a.minus)) {
                mpy <- ifelse(dom.a.pr==0.0,0.0,(1/dom.a.pr/unc.pr^2))
                should.be <- as.matrix((left.marker[, is.a.minus, i] * 
                                        right.marker[, is.a.minus, i])*mpy)
                

                center.marker[, is.a.minus, i ] <-
                    sweep(should.be, 2, apply(should.be, 2, sum.or.one), "/")
            }
        }
    }

    aperm(center.marker, c(2, 3, 1))
}
"make.varcov"<-
    function(regressor.matrix, y, subset = is.finite(y), casewt = NULL)
{
    df <-
        if (is.logical(subset))
            sum(subset) - 1
        else
            length(y[subset])-1
    if(is.null(casewt)) {
        var.x <- var(regressor.matrix[subset,  ]) * df
        cov.xy <- var(regressor.matrix[subset,  ], y[subset]) * df
        var.y <- var(y[subset]) * df
    }
    else {
        casewt <- as.vector(casewt[subset])
        mean.x <- matrix(casewt/sum(casewt), nr = 1) %*%
            as.matrix(regressor.matrix[subset,  ])
        mean.y <- sum(y[subset] * casewt)/sum(casewt)
        center.x <-
            as.matrix(regressor.matrix[subset,  ]) -
                rep(mean.x, rep(sum(subset), length(mean.x)))
        center.y <- y[subset] - mean.y
        var.x <- t(center.x) %*% (casewt * center.x)
        cov.xy <- t(center.x) %*% (casewt * center.y)
        var.y <- sum(casewt * center.y^2)
    }
    res <- list(var.x = var.x, var.y = var.y, cov.xy = cov.xy, df = df)
    attr(res,"call") <- match.call()
    class(res) <- "varcov"
    res
}
"map.dx"<-
    function(lambda = (1 - 2 * theta), theta = NULL, min.lambda = 0)
{
    if(length(lambda) && any(lambda < min.lambda)) {
        lambda <- pmax(lambda, min.lambda)
        warning("small value(s) for lambda corrected")
    }
    - log(lambda)/2
}
"map.index"<-
    function(x,...)
    UseMethod("map.index")
"map.index.default"<-
    function(x,chromo,cM=NULL)
{
    if (length(chromo)==1 && missing(cM) || length(cM)!=2){
        subset <-    x$chr.num==chromo
        if (!missing(cM)){
            this.cM <- x[ subset ,"cM" ]
            loc <-
                seq(along=this.cM)[ abs(this.cM-cM) == min(abs(this.cM-cM)) ][1]
            seq(along=subset)[subset][loc]
        }
        else
        {
            seq(along=subset)[subset]
        }
    }
    else
    {# two locations find a loci in span
      if (length(chromo)==1) chromo[2] <- chromo
        if (length(chromo)>2) stop("> 2 values for chromo not allowed")
        if (chromo[1]>chromo[2]) stop("chromo[1]>chromo[2]")
        subset <-
            ( (x$chr.num <= chromo[2]) & (x$chr.num >= chromo[1] ) )
        if (missing(cM)){
            return( seq(along=subset)[subset] )
        }
        else {
            if (length(cM) != 2) stop("cM must have 2 elements")
            if (chromo[2]==chromo[1] && cM[1]>cM[2]) stop("cM[1]>cM[2]")
            seq(from=map.index(x,chromo[1],cM[1]),
                to=map.index(x,chromo[2],cM[2]))
        }
    }
}
"map.index.analysis.object"<-
    function(x,...)
    map.index(x$map.frame,...)
"map.names" <-
    function(x,...)
{
    UseMethod("map.names")
}
"map.names.default"<-
    function(x,...)
{
    names(x,...)
}
"map.names.map.frame" <-
    function(x,...)
{
    c(x$marker.name[map.index(x,...)])
}
"map.names.analysis.object"<-
    function(x,...)
    map.names(x$map.frame,...)
"map.names.bqtl"<-
    function(x,...,ana.obj=NULL)
{
    if (exists("is.R") && is.R())
        tokens <- unlist(strsplit(x$reg.vec,":"))
    else
        tokens <-  unlist(lapply(x$reg.vec,unpaste,":"))
    covar.used <- tokens[grep("covar\\(",tokens)]
    
    if (is.null(ana.obj))
        ana.obj <- eval(match.call(bqtl,x$call)$ana.obj)
    
    
    reg.names.indx <- ana.obj$reg.names%in%tokens
    if (length(covar.used) != 0)
        reg.names.indx <-
            reg.names.indx | ( paste("covar(",ana.obj$reg.names,")",sep="")%in%
                              covar.used)
    res <-
        if (ana.obj$method == "F2")
            unique(
               {x<-dimnames(ana.obj$reg.names)[[2]];
                rep(x,rep(2,length(x)))[reg.names.indx]})
        else
            unique(dimnames(ana.obj$reg.names)[[2]][reg.names.indx])
    res
}
"map.names.bqtl.list"<-
    function(x,...)
{
    res <- list()
    for (i in seq(along=x))
        res[[i]]<-map.names.bqtl(x[[i]],...)
    res
}
"map.location"<-
    function(x,...)
    UseMethod("map.location")
"map.loc"<-map.location
"map.location.default"<-
    function(x,y,chromo=NULL,cM=NULL,map.names=NULL)
{
    if (missing(y)){
        sw.arg <-    1+2*is.null(chromo)+is.null(map.names)
        y <- switch(sw.arg,
                    ##both non-null
                    stop("cannot use both chromo= and map.name ="),
                    ## chromo
                    if (all(chromo%in%x$chr.num)){
                        if (is.null(cM))
                            x$chr.num %in% chromo
                        else
                            map.index(x,chromo,cM)
                    }
                    else
                    stop("bad chromo number specified"),
                    ## map.name
                    if (all(map.names%in%x$marker.name)){
                        x$marker.name%in%map.names}
                    else
                    stop(paste("unrecognized map.name:",
                               map.names[!map.names%in%x$marker.name])),
                    ## both null use all names
                    seq(nrow(x))
                    )
    }
    res <- if (mode(y)=="character")
        map.location(x,map.names=y)
    else
        x[y,c("chr.num","cM","marker.name")]
    class(res) <- c("map.location",class(x))
    res
}
"map.location.analysis.object"<-
    function(x,...)
    map.location(x$map.frame,...)
"map.location.bqtl"<-
    function(x,ana.obj=NULL)
{
    if (is.null(ana.obj))
        ana.object <- eval(match.call(bqtl,x$call)$ana.obj)
    mnames <-    map.names(x,ana.obj=ana.object)
    map.location(ana.object,mnames)
    
}
"map.location.bqtl.list"<-
    function(x,...)
{
    res <- list()
    for ( i in seq(along=x))
        res[[i]] <- map.location.bqtl(x[[i]],...)
    res
}
"marker.fill"<-
    function(map.frame, reso , return.nint=FALSE)
{
### given some markers and a desired resolution
### add pseudo-markers at 'good' distances
### and create names that show their relative
### positions
### map.frame$mgd is the distance in M or cM from the left telomere
### reso is in same metric
    int.suffix <- function(nint)
    {
        res <- sapply(nint, function(x)
                  {
                      strng <- if(x > 1) {
                          c("", format((1:(x - 1))/10^(trunc(log10(x - 1) +
                                                             1))))
                      }
                      else {
                          ""
                      }
                      substring(strng, 2)
                  }
                      )
        unlist(res)
    }
    
### assume that map.frame is a map.frame.object
### start of each chromosome
    if (any(is.na(pmatch(c("cM","marker.name","chr.num"),names(map.frame)))))
        stop("need $cM, $marker.name, and $chr.num components in first arg")
    dx <- map.frame$cM
    names(dx) <-   map.frame$marker.name
    mgd <- tapply(dx, map.frame$chr.num,c)
    names(mgd) <- rep("",length(mgd))
    
    res.main <- lapply(mgd, function(x, reso, int.suffix, return.nint)
                   {
                       genes <- grep("gene", names(x))
                       if(length(genes) > 0)
                           mark.loc <- x[ - genes]
                       else mark.loc <- x
###     nint <- pmax(1, (diff(mark.loc)+reso/2) %/% reso)
                       nint <- pmax(1, ceiling(diff(mark.loc) / reso))
                       if(return.nint) {
                           nint <- c(nint, 1)
                           names(nint) <- names(mark.loc)
                           return(nint)
                       }
                       else
                       {
                           dxes <- c(rep(diff(mark.loc)/nint, nint), Inf)
                           names(dxes) <-
                               paste(rep(names(mark.loc), c(nint, 1)), 
                                     int.suffix(c(nint, 1)), sep = "")
                           return(dxes)
                       }
                   }
                       , reso, int.suffix,return.nint)
    unlist(res.main)
}
### convenience functions for defining marker levels
###
if (!( exists("is.R") && is.R() ))          # S 3.4 needs this
    formals <- function(x) x[-length(x)]
###
"f2.levels" <-
    function(AA="AA", Aa="Aa", aa="aa", not.aa="A-", not.AA="a-",
             miss.val="--")
{
    this.call<-match.call()
    formal.names <- names(formals(eval(this.call[[1]])))
    res <- c(AA,Aa,aa,not.aa,not.AA,miss.val)
    names(res) <- formal.names
    res
}

"bc1.levels" <-
    function(AA="AA", Aa="Aa", miss.val="--")
{
    this.call<-match.call()
    formal.names <- names(formals(eval(this.call[[1]])))
    res <- c(AA, Aa, miss.val )
    names(res) <- formal.names
    c(res[1:2],rep("nil",3),res[3])
}

"ri.levels" <-
    function(AA="AA", aa="aa", miss.val="--")
{
    this.call<-match.call()
    formal.names <- names(formals(eval(this.call[[1]])))
    res <- c(AA, aa, miss.val )
    names(res) <- formal.names
    c(res[1:2],rep("nil",3),res[3])
}

"logpost"<-
    function(x,...)
    UseMethod("logpost")
"loglik"<-
    function(x,...)
    UseMethod("loglik")


"logpost.bqtl" <-
    function(bqtl.obj,...)
{
    bqtl.obj$logpost
}

"logpost.default"<-logpost.bqtl
"loglik.bqtl"<-logpost.bqtl
"loglik.default"<-logpost.bqtl

"logpost.bqtl.list"<-
    function(bqtl.list,...)
{
    res <- sapply(bqtl.list,logpost.bqtl)
    res
}

"loglik.bqtl.list"<-logpost.bqtl.list

"coef.bqtl"<-
    function(bqtl.obj,...)
{
    ncoefs <- length(bqtl.obj$parm)-1
    coefs <- bqtl.obj$parm[1:ncoefs]
    names(coefs) <- c("Intercept",bqtl.obj$reg.vec)
    coefs
}

### method dispatch for coef takes care of this: "coefficients.bqtl"<-coef.bqtl

"coef.bqtl.list"<-
    function(bqtl.list,...)
{
    res <- sapply(bqtl.list,coef.bqtl)
    res
}

### method dispatch for coef takes care of this: "coefficients.bqtl.list"<-coef.bqtl.list

"summary.bqtl"<-
    function(bqtl.obj,...)
{
    nparm <-    length(bqtl.obj$parm)
    coefs <-
        if (is.null(bqtl.obj$hess))
            coef(bqtl.obj)
        else
        {
            vals <- coef(bqtl.obj)
            hess.qr <- qr(-bqtl.obj$hess)
            if (hess.qr$rank < nparm){
                use.coef <-  vals != 0
                hrc <- c(use.coef,TRUE)
                std.err <- ifelse( use.coef, 0, NA)
                std.err[use.coef]<-
                    sqrt(diag(solve(-bqtl.obj$hess[hrc,hrc]))[-hess.qr$rank])
            }
            else{
                std.err <- sqrt(diag(solve(hess.qr))[-nparm])
            }
            t.stats <- vals/std.err
            p.vals <- pnorm(-abs(t.stats))*2
            cbind(Value=vals,Std.Err=std.err,t.statistic=t.stats,p.value=p.vals)
        }
    
    loglik <- loglik(bqtl.obj)
    sigma <- exp(bqtl.obj$parm[nparm])
    list(coefficients=coefs,loglik=loglik,std.res=sigma,N=bqtl.obj$N)
}

"posterior"<-
    function(x,...)
    UseMethod("posterior")

"posterior.bqtl" <-
    function(bqtl.obj,...)
{
    bqtl.obj$posterior
}

"posterior.default" <-posterior.bqtl

"posterior.bqtl.list"<-
    function(bqtl.list)
{
    res <- sapply(bqtl.list,posterior.bqtl)
    res
}

"plot.analysis.object"<-
  function(x,...)
  plot(x$map.frame,...)

  

"plot.map.frame"<-
    function(mf.obj,y, fun = if (y.type=="matrix") matlines else lines,
             type="l", include.rug=TRUE, rug.lwd=0.1,title.string=NULL,
             y.range=NULL,
             ylab=deparse(substitute(y)),xlab="Location", ...)
{
    if (any(is.na(pmatch(c("chr.num","pos.plot"),names(mf.obj)))))
        stop("need $chr.num and $pos.plot components in first arg")
### try to be smart: guess if this is a vector or matrix arg.
    if (missing(y)){
        y.range <- range(mf.obj$chr.num)
        x.range <- range(unlist(mf.obj$cM))
        plot(x.range,y.range,type="n",xlab=xlab,
             ylab= if (missing(ylab)) "Chromosome" else ylab)
        for (i in unique(mf.obj$chr.num)) {
            i.index <- i==mf.obj$chr.num
            x <- mf.obj$pos.plot[i.index]
            x.range <- range(x)
            lines(x.range,c(i,i))
            seg.col <- c(1,2)
            seg.lty <- c(1,3)
            if (any( which.seg <- mf.obj$is.marker[i.index] )) {
                n.seg <- sum(which.seg)
                segments(x[which.seg],rep(i,n.seg)-0.1,
                         x[which.seg],rep(i,n.seg)+0.1,lty=seg.lty[1],
                         col=seg.col[1])}
            if (any( which.seg <- !mf.obj$is.marker[i.index] )) {
                n.seg <- sum(which.seg)
                segments(x[which.seg],rep(i,n.seg)-0.1,
                         x[which.seg],rep(i,n.seg)+0.1,lty=seg.lty[2],
                         col=seg.col[2])}
        }
        title(paste(title.string,deparse(substitute(mf.obj))))
        return(invisible())
    }
    
    y.type <-
        if ( length(dmm <- dim(y)) && max(dmm)<length(y))
            "matrix"
        else
            "vector"
    if (is.null(y)) stop("NULL value for y not accepted in plot.map.frame")
    if (is.null(y.range))
        y.range <- range(y,na.rm=T)

    if (exists("is.R") && is.R() && y.type=="matrix") { # S compat workaround
        matlines <-
            function(x, y, type = 'l', lty=1:5, lwd = 1, pch=NULL, col=1:6, ...)
                matplot(x=x, y=y, type=type, lty=lty, lwd=lwd, pch=pch, col=col,
                        add=TRUE, ...)
        matpoints <-
            function(x, y, type = 'p', lty=1:5, lwd = 1, pch=NULL, col=1:6, ...)
                matplot(x=x, y=y, type=type, lty=lty, lwd=lwd, pch=pch, col=col,
                        add=TRUE, ...)
    }

    for (i in unique(mf.obj$chr.num)) {
        i.index <- i==mf.obj$chr.num
        yp <- if (y.type=="matrix") y[i.index,,drop=F] else y[i.index]
        if (length(yp)==0) 
            print(paste("no data for chr.num = ",i))
        else { 
            x <- mf.obj$pos.plot[i.index]
            x.range <- range(x)
            plot(x.range,y.range,type="n",ylab=ylab,xlab=xlab )
            fun( x, yp, type=type, ...)
            if (include.rug) rug(x[mf.obj$is.marker[i.index]],lwd=rug.lwd)
            title(paste(title.string,"Chromosome", i))
        }
    }
    NULL
    invisible()
}


"predict.bqtl"<-
    function(object,newdata)
{
    coefs <- coef(object)
    coef.names <- names(coefs)[-1]        # drop intercept
    bqtl.call <- match.call(bqtl,object$call)
    reg.form <- eval(bqtl.call$reg.formula)
###    bqtl.terms <- terms(bqtl.call$reg.formula)
    if (missing(newdata)) {
        data.obj <- bqtl.call[["ana.obj"]]
        newdata <- substitute(data.obj$data)
        bqtl.terms <- terms(reg.form)
        model.fr <-
            model.frame(bqtl.terms,data=eval(newdata),
                        na.action=na.omit)
    }
    else {
        bqtl.terms <- delete.response(terms(reg.form))
        if (class(newdata) == "analysis.object")
            model.fr <- model.frame(bqtl.terms,data=eval(newdata$data),
                                    na.action=na.omit)
        else
            model.fr <- model.frame(bqtl.terms,data=eval(newdata),
                                    na.action=na.omit)
    }
    
    model.mat <- model.matrix(delete.response(bqtl.terms),model.fr)
    
    coef.order <-  1 +                   #assume intercept is first
        match(dimnames(model.mat)[[2]][-1],coef.names)

    fit <-  model.mat %*% as.matrix( coefs[c(1,coef.order)] )
    dimnames(fit)[[2]] <- "fitted.value"
    fit
}
"fitted.bqtl"<-
    predict.bqtl
"residuals.linear.bayes" <-
    function(object,...)
    predict.linear.bayes(object,...,return.resids=TRUE)

### since most of the code is the same in residuals and predict, why
### duplicate it?

"fitted.linear.bayes" <-
    function(...)
    predict.linear.bayes(...)

"predict.linear.bayes"<-
    function(object,newdata=lb.call$ana.obj,return.resids=FALSE)
{
    coefs <- coef(object)
    lb.call <- match.call(linear.bayes,object$call)
    x <- eval(lb.call$x)
    partial <- eval(lb.call$partial)
    if (x[[1]]!="~"){          #not a formula:try to go on...
        if (return.resids)
            stop("cannot create residuals without formula in linear.bayes()")
        newdata <- eval(substitute(newdata))
        preds <- as.matrix(newdata)%*%coefs
        preds <- preds - mean(preds)    #center the results
        return(preds)
    }
    ana.obj <- eval(newdata)
    scope <- ana.obj$reg.names
    
### take this from varcov:
    bqtl.specials <- c("configs","locus") # dont use: add or dom so
                                        # swap, hk wont jam up    
    reg.terms <- terms(x,specials=bqtl.specials)
    if (any(attr(reg.terms,"order")>1))
        stop("cannot use interactions in this formula")
    covs <-
        rhs.bqtl(reg.terms,ana.obj,bqtl.specials,NULL,scope,
                 expand.specials=FALSE)
    response <-
        dimnames(attr(reg.terms,"factors"))[[1]][attr(reg.terms,"response")]
    if (!is.null(partial)) {
        p.terms <- terms(partial,specials=bqtl.specials)
        if (any(attr(p.terms,"order")>1))
            stop("cannot use interactions in this formula")
        cntl <-
            rhs.bqtl(p.terms,ana.obj,bqtl.specials,NULL,scope,
                     expand.specials=FALSE)
        dummy.formula <- eval(parse(text=paste(response,"~",covs,"+",cntl)))
    }
    else {
        cntl <- NULL
        dummy.formula <- eval(parse(text=paste(response,"~",covs)))
    }
    dummy.terms <- terms(dummy.formula)
    mf <- model.frame(dummy.terms,ana.obj$data, na.action=na.omit)
    if (is.null(partial)) {
        x <- model.matrix(dummy.terms,mf)[,-1]
        y <- model.extract(mf,"response")
        preds  <- x%*%matrix(coefs,nc=1)
        ## can't always get an intercept from linear.bayes, hence:
        preds<-preds + mean(y) - mean(preds)
    }
    else
    {
        
        cntl.out <- eval(parse(text=paste("~ . -", cntl)))
        x.formula <- update(dummy.formula,
                            eval(parse(text=paste("~ . -", cntl)))  )
        x <- model.matrix(terms(x.formula),mf)[,-1]
        y <- model.extract(mf,"response")
        z <- model.matrix(terms(eval(parse(text=paste("~", cntl)))),
                          mf)[,-1,drop=F]
        pred.z <- y - lsfit(z,y)$resid
        preds <- x%*%matrix(coefs,nc=1) + pred.z
        preds <- preds + mean(y) - mean(preds)
    }
    if (return.resids){
        y - preds
    }
    else
        preds
}

"residuals.bqtl"<-
    function(object)
{
    coefs <- coef(object)
    coef.names <- names(coefs)[-1]        # drop intercept
    bqtl.call <- match.call(bqtl,object$call)  
    bqtl.terms <- terms(eval(bqtl.call$reg.formula))
    data.obj <- bqtl.call[["ana.obj"]]
    newdata <- substitute(data.obj$data)
    model.fr <-
        model.frame(bqtl.terms,data=eval(newdata),
                    na.action=na.omit)
    model.mat <- model.matrix(delete.response(bqtl.terms),model.fr)
    response <- model.extract(model.fr,"response")
    
    coef.order <-  1 +                   #assume intercept is first
        match(dimnames(model.mat)[[2]][-1],coef.names)
    
    resids <-   response -model.mat %*% as.matrix( coefs[c(1,coef.order)] )
    resids
    
}
"rhs.bqtl"<-
    function(reg.terms,ana.obj,bqtl.specials,local.covar,
             scope, expand.specials=NULL,method,...)
{
    using.R <- exists("is.R")&&is.R()
    reg.labels <- labels(reg.terms)
    reg.specials <- attr(reg.terms,"specials")  
    if (missing(method))
        method <- ana.obj$method
### work on attr(,"factors")
    reg.factor <- attr(reg.terms,"factors")
    term.rownames <- dimnames(reg.factor)[[1]]
    term.names <- term.rownames[row(reg.factor)[reg.factor>0]]
    terms.conjuncs <-
        if ( sum(reg.factor>0) <2 )
            NULL
        else
            ifelse( diff(col(reg.factor)[reg.factor>0])==0,"colon", "plus")
    n.join <- length(terms.conjuncs)
### flag rows with specials and those without
    reg.specials <- unlist(reg.specials[ !sapply(reg.specials,is.null) ])
    reg.plain <-
        if (is.null(reg.specials))
            term.rownames
        else
            term.rownames[ - reg.specials ]
    names(reg.plain) <- reg.plain
    
    
### let the specials expand themselves
    if (using.R)
        formals(local.covar)$bq.spec <- bqtl.specials # bind bqtl.specials explicitly
    else
        local.covar$bq.spec <- bqtl.specials
    pt.vars <- reg.specials + if (using.R) 1 else 0
    if (length(reg.specials) != 0) {
        rspec <-
            lapply(attr(reg.terms,"variables")[pt.vars],
                   function(x,scope,method,covar) {
                       if ( !is.element("scope",names(x)) ) #typically use default
                           x$scope <- as.name("scope")
                       if ( !is.element("method",names(x)) ) #typically use default
                           x$method <- method
                       eval(x)
                   },
                   scope=scope,method=method,
                   covar=local.covar)
        names(rspec) <- term.rownames[ reg.specials ]
### used a common 'scope' and 'method' for all specials
        rspec.check <-
            sapply(rspec, function(x) any(x %in% c("(NA)","()")))
        if (any(rspec.check)){
            bad.terms <-
                paste(c("invalid term(s) in formula:",names(rspec)[rspec.check]),
                      collapse=" ")
            stop(bad.terms)
        }
        
### use all combinations of the expanded variables ?
        if (is.null(expand.specials))
            expand.specials <- 
                length(rspec)>1 && any(diff(range(sapply(rspec,length)))!=0)

        if (expand.specials)
            rspec <- do.call("expand.grid",rspec)
        else
            rspec <- lapply(rspec,paste,collapse="+")
    }
    else { # no specials
        rspec <- NULL
    }
    if (using.R ){
        term.list <- c(rspec,as.list(c(reg.plain,plus="+",colon=":")))
    }
    else { #argh! Splus 3.4
        rspec <- lapply(rspec,as.character)
        term.list <- c(rspec,as.list(c(reg.plain,plus="+",colon=":")))
        names(term.list) <- c(names(rspec),names(reg.plain),"plus","colon")
  }
  
### order is <var,conj,var,conj,...,conj,var>
    spec.col.order <-
        if (length(terms.conjuncs)==0)
            term.names
        else
            c(term.names,terms.conjuncs)[c(rep(c(0,n.join+1),n.join)+
                                           rep(1:n.join,rep(2,n.join)),n.join+1)]
    do.call("paste",term.list[spec.col.order])
}
"summary.adj"<-
function(adj.obj, n.loc, coef.znames, mode.names = c("add", "dom"), imp.denom
	 = NULL, swap.obj = NULL)
{

### imp.denom is any non-null to switch on use of swap.obj$hk.exact
### this amounts to importance sampling from a swap.obj that differs
### in params from adj.obj to use enumerations (like one gene case)
### use imp.denom=1/elementwise.prior
###

  adj.unlist <-
    function(x,component="adj")
      sapply(x,"[[",component)
  adj.type <- if (is.null(imp.denom)) {
    if (is.null(swap.obj))
      "uniform"
    else
      "swap"
  }
  else {
    if (is.null(swap.obj))
      {
        imp.denom <- 1
        "uniform"
      }
    else
      "imp.swap"
  } 
  if(is.null(mode.names)) {
    coef.unames <- coef.znames
    n.mode <- 1
  }
  else {
    coef.unames <- outer(mode.names, coef.znames, paste, sep = ".")
    n.mode <- length(mode.names)
  }

  
  if(is.element(adj.type,c("swap","imp.swap"))) {
    config <- unique.config(swap.obj)
    if(max(config$match) != length(adj.obj))
      stop("incompatible argument lengths:adj.obj and swap.obj"
           )
  }
  else
    config <- list(match=seq(along = adj.obj),
                   uniq=rbind(seq(along = adj.obj)))
  
  if (is.null(adj.obj[[1]]$reg.vec)) { # take locs from swap obj or default
    locs <- (config$uniq[,config$match] - 1)%/%n.mode +1
    coef.ind <- config$uniq[,config$match]
    coef.ind <- coef.ind[coef.ind>0]
  }
  else { # figure out locs from adj.obj[[i]]$reg.vec
    locs <-
      factor(unlist(lapply(adj.obj[config$match], function(x, nm, n.mode)
                           unique((match(x$reg.vec, nm) - 1) %/% n.mode) + 1,
                           nm = coef.unames, n.mode = n.mode)),
             1:length(coef.znames))
    coef.names <- unlist(lapply(adj.obj, "[[", "reg.vec")[config$match])
    coef.ind <- factor(coef.names, coef.unames)
  }
  coef.multiple <- unlist(lapply(adj.obj, function(x) length(x$parm) - 2))[config$match]
  
  coefs <- unlist(lapply(adj.obj, function(x) x$parm[ - c(1, length(x$parm))])[config$match])
  
  
  adj <- adj.unlist(adj.obj)[config$match]
  post <- adj*unlist(adj.unlist(adj.obj,component="hk.exact"))[config$match]
  switch(adj.type,
       uniform={
         post <- post/imp.denom
         adj.mean <- mean(adj/ (imp.denom/sum(imp.denom)) )
         adj.var <- 0
         hk.ratio.mean <- 1
         loc.post <- tapply(rep(post, rep(n.loc, length(post))),
                            locs, sum)
         loc.post[is.na(loc.post)] <- 0
         loc.post <- loc.post/sum(loc.post)
         coef.avg <- tapply(coefs*rep(post,coef.multiple), coef.ind, sum)/sum(post)
         coef.avg[is.na(coef.avg)] <- 0
         },
       swap={
         adj.mean <- mean(adj)
         adj.var <- var(c(apply(matrix(adj, nr = n.loc), 2, mean)))/(
                                              length(adj)/n.loc)
         hk.ratio.mean <- 1
         loc.post <- tapply(rep(adj, rep(n.loc, length(adj))), locs, sum)
         loc.post[is.na(loc.post)] <- 0
         loc.post <- loc.post/sum(loc.post)
         
         coef.adj <- rep(adj, coef.multiple)
         coef.avg <- tapply(coef.adj * coefs, coef.ind, sum)/sum(adj)
         coef.avg[is.na(coef.avg)] <- 0
       },
       imp.swap={
         post <- post/swap.obj$hk.exact
         adj.mean <- mean(post)
         hk.ratio <- adj.unlist(adj.obj, comp="hk.exact")[config$match]/swap.obj$hk.exact
         hk.ratio.mean <- mean(hk.ratio)	
         adj.var <- var(c(apply(matrix(adj * hk.ratio, nr = n.loc), 2, 
                                mean)))/(length(adj)/n.loc)
         loc.post <- tapply(rep(post, rep(n.loc, length(post))), locs, mean)
         loc.post[is.na(loc.post)] <- 1
         loc.post <- loc.post * swap.obj$alt.marg
         loc.post <- loc.post/sum(loc.post)
         coef.adj <- rep(post, coef.multiple)
         coef.avg <- tapply(coef.adj * coefs, coef.ind, sum)/sum(post)
         coef.avg[is.na(coef.avg)] <- 0
       }
       )
  list(adj = adj.mean, var = adj.var, coef = coef.avg, loc = loc.post, 
       hk.ratio.mean = hk.ratio.mean)
}
"summary.analysis.object" <-
    function(x)
{
    mf <- summary(x$map.frame)
    norg <- nrow(x$data)
    method <- x$method
    origin<-x$call
    pheno.cov <- setdiff(names(x$data),x$reg.names)
    list(map=mf,N=norg,method=method,variable.names=pheno.cov,
         mode.mat=x$mode.mat,call=origin)
}

"summary.map.frame" <-
    function(x)
{
    chrs<-table(x$chr.num,
                factor(x$is.marker,c(TRUE,FALSE),c("marker","pseudo-marker")))
    lens <- do.call("rbind",tapply(x$cM,x$chr.num,range))
    res <- cbind(chrs,lens)
    dimnames(res)[[2]][3:4]<-c("cM.low","cM.high")
    res
}

"summary.swap"<-
function(sc.obj, method=NULL, ncoef=length(sc.obj$alt.coef),nloc=sc.obj$nloc)
{
### summarize a swap.cycle object
  this.call <- sys.call()
  if (missing(method))
    method <-
      if(diff(dim(sc.obj$config))[1] < 0 )
        "F2"
      else
        "BC1" #includes RIL cases
  if(method=="F2") {
    if (is.null(ncoef)) ncoef <- ((max(sc.obj$config) - 1) %/% 2 + 1) * 2
    if (is.null(nloc)) nloc <- ncoef/ 2
  }
  else
    {
      if (is.null(ncoef)) {
        if (is.null(nloc))
          ncoef <- nloc <- max(sc.obj$config)
        else
          ncoef <- nloc
      }
      else
        nloc <- ncoef
    }
  coefs <- rep(0, ncoef)
  rep.post <- rep(sc.obj$post, rep(nrow(sc.obj$coefs), length(sc.obj$post
                                                              )))
  cp <- sc.obj$coefs * rep.post
  tap.coef <- tapply(cp, sc.obj$conf, sum)
  coefs[sort(unique(sc.obj$config))] <- if (method=="F2") tap.coef[-1] else tap.coef
  coefs <- coefs/sum(sc.obj$post)
  if (method=="F2") {
    locs <- apply((sc.obj$config - 1) %/% 2 + 1, 2:3, function(x) {unique(x[x > 0])})
  }
  else
    locs <- sc.obj$config
  diml <- dim(locs)
  if(length(diml) == 2)
    diml <- c(1, diml)
  loc.post <- rep(0, nloc)
  rep.marg <- rep(sc.obj$marg, rep(diml[1], diml[2] * diml[3]))
  loc.post[sort(unique(locs))] <-
    tapply(rep.marg, locs, sum)/sum(rep.marg)*diml[2]
  ## note assumption that second subscript tells size of model
  blk.ratio <- apply(sc.obj$cond/sc.obj$marg, 2, mean)
  res <- list(loc.posterior = loc.post, coefs = coefs,
              ratio = list(mean = mean(blk.ratio),
                var = var(blk.ratio)/length(blk.ratio)),
              setup=c(genes=diml[2],nreps=diml[3],nloc=nloc)
              ,call=this.call)
  class(res) <- "summary.swap"
  res
}
"swap" <-
    function(varcov, invars, rparm, nreps, ana.obj, ...)
{
    this.call <- match.call()
    if (class(ana.obj) != "analysis.object")
        stop("arg 5 must be analysis.object")
    method <- ana.obj$method
    res <-
        if (method == "F2") {
            swapf2(varcov, invars, rparm, nreps, ana.obj, ...)
        }
        else {
            swapbc1(varcov, invars, rparm, nreps, ana.obj, ...)
        }
    res$call <- this.call
    res
}
"swapbc1" <-
    function(varcov,invars,rparm,nreps,ana.obj,
             locs=NULL,
             locs.prior = NULL,
             tol=1e-10)
{
    this.call <- match.call()
    if (any(duplicated(invars)))
        stop("duplicate values of invars not allowed")
    nin <- length(invars)

    if (!missing(ana.obj)) {
        if (class(ana.obj) != "analysis.object")
            stop("arg \'ana.obj\' must be analysis object")
        if (missing(locs.prior)) {
            match.names <-
                match(dimnames(varcov$var.x)[[1]],ana.obj$reg.names,0)
            if (any(match.names==0))
                stop("variable names in varcov not found in ana.obj")
            locs.prior <-
                ana.obj$map.frame[match.names,"prior"]

            locs.prior <- locs.prior/sum(locs.prior)
        }
    }
    else {
        if (missing(locs.prior))
            locs.prior <- rep(1, ncol(varcov$var.x))
    }
    if (missing(locs)) 
        locs <- seq(ncol(varcov$var.x))

    in.locs <- unique(locs[invars])
    nlocs <- length(in.locs)
    nstep <- nlocs
    
    len.locs <- length(unique(locs))
    
    nopt <- len.locs
    optmax <- 1
    nrx <- nrow(varcov$var.x)
    if (nrx!=length(rparm)) {
        if (length(rparm)==1)
            rparm <- rep(rparm,nrx)
        else
            stop(paste("rparm must have length 1 or ",nrx))
    }
    if (length(locs.prior)!=len.locs)
        stop("length(locs.prior)!=length(unique(locs))")
    if (length(nreps)>1) stop("nreps should be an integer")
    optpri <- c(locs.prior)
    optvars <- locs - 1
    optcur <- in.locs
    noptuse <- length(optcur)

    optblock <- (1:nopt)-1
    useopt <- rep(1,nopt)
    useopt[c(optcur,optblock[optcur]+1)] <- 0
    inpri <- 1
    xmax <- (nlocs-1)
    zmax <- 1
    nmax <- xmax+zmax
    invars <- c(invars-1)
    z <- .C("swapbc1",
            nreps=as.integer(nreps),
            nstep=as.integer(nstep),
            varx=as.double(varcov$var.x),
            covxy=as.double(varcov$cov.xy),
            vary=as.double(varcov$var.y),
            df=as.double(varcov$df),
            rparm=as.double(rparm),
            optpri=as.double(optpri),
            inpri=as.double(inpri),
            nrx=as.integer(nrx),
            noptuse=as.integer(noptuse),
            optblock=as.integer(optblock),
            optcur=as.integer(optcur-1),
            invars=as.integer(invars),
            nin=as.integer(nin),
            optvars=as.integer(optvars),
            nopt=as.integer(nopt),
            optmax=as.integer(optmax),
            useopt=as.integer(useopt),
            locs=as.integer(locs-1),
            pvt=as.integer(1:nmax),
            rank=integer(1),
            wrksp=double(2*nmax),
            gama=double(zmax*xmax),
            bee=double(xmax),
            xx=double(xmax^2),
            xy=double(xmax),
            zz=double(zmax^2),
            zy=double(zmax),
            xz=double(xmax*zmax),
            beta=double(xmax),
            posteriors=double(nreps*nstep),
            marg=double(nreps*nstep),
            cond=double(nreps*nstep),
            coefs=double(nreps*nstep*nmax),
            configs=integer(nreps*nstep*nmax),
            qraux=double(xmax),
            zraux=double(zmax),
            zrank=integer(1),
            tol=as.double(tol),
            postwk=double(nopt),
            coefwk=double(nopt*nmax),
            alt.marginal=double(nopt),
            alt.coef=double(nrx))
    res <- z[c("configs","posteriors","coefs","cond","marg",
               "alt.marginal","alt.coef")]
    res$configs <- array(res$configs+1,c(nmax,nlocs,nreps))
    dim(res$coefs) <- c(nmax,nlocs,nreps)
    dim(res$marg) <- dim(res$cond) <- dim(res$posteriors) <- c(nlocs,nreps)
    res$nloc <- nopt
    res[["call"]] <- this.call
    class(res)<-"swap"
    res
}
"swapf2" <-
    function(varcov,invars,rparm,nreps,ana.obj,
             locs=NULL,
             locs.prior = NULL,
             combo.prior = rep(1, 3)/3,tol=1e-10)
{
    this.call <- match.call()
    if (any(duplicated(invars)))
        stop("duplicate values of invars not allowed")
    nin<- length(invars)
    if (!missing(ana.obj)) {
        if (class(ana.obj) != "analysis.object")
            stop("arg \'ana.obj\' must be analysis object")
        if (missing(locs.prior)) {
            match.names <-
                match(dimnames(varcov$var.x)[[1]],ana.obj$reg.names,0)
            match.cols <- unique(col(ana.obj$reg.names)[match.names])
            if (any(match.names==0))
                stop("variable names in varcov not found in ana.obj")
            locs.prior <-
                ana.obj$map.frame[match.cols,"prior"]
            locs.prior <- locs.prior/sum(locs.prior)
        }
    }
    else {
        if (missing(locs.prior))
            locs.prior <- rep(1, ncol(varcov$var.x)/2)
    }

    if (missing(locs))
        locs <- rep(seq(ncol(varcov$var.x)/2), rep(2, ncol(varcov$var.x)/2))
    
    in.locs <- locs[invars]
    uniq.locs <- unique(in.locs)
    loc.tab <- table(in.locs)
    mode.offset <- ifelse(loc.tab[as.character(in.locs)]==2,2,(invars-1)%%2)
    mode.offset <- mode.offset[ match( uniq.locs, in.locs ) ]
    nlocs <- length(uniq.locs)
    nstep <- nlocs
    len.locs <- length(unique(locs))
    
    nopt <- len.locs*3
    optmax <- 2
    nrx <- nrow(varcov$var.x)
    if (length(combo.prior)!=3) stop("combo.prior must = 3")
    if (nrx!=length(rparm)) {
        if (length(rparm)<3)
            rparm <- rep(rparm,nrx/length(rparm))
        else
            stop(paste("rparm must have length 1, 2, or ",nrx))
    }
    if (length(locs.prior)!=len.locs)
        stop("length(locs.prior)!=length(unique(locs))")
    if (length(nreps)>1) stop("nreps should be an integer")
    optpri <- c(locs.prior)%o%c(combo.prior)
    optvars <- matrix(c(rbind(unique(locs)*2-2,-1),rbind(unique(locs)*2-1,-1),
                        locs*2-rep(2:1,nopt/3)),nrow=2)

    optcur <- uniq.locs + mode.offset*len.locs
    noptuse <- length(optcur)
    optblock <-
        rbind(rep(1:len.locs,3)+rep(c(len.locs,0,0),rep(len.locs,3)),
              rep(1:len.locs,3)+
              rep(c(2*len.locs,2*len.locs,len.locs),rep(len.locs,3)))-1
    useopt <- rep(1,nopt)
    useopt[c(optcur,optblock[,optcur]+1)] <- 0
    inpri <- 1
    xmax <- 2*(nlocs-1)
    zmax <- 2
    nmax <- xmax+zmax
    invars <- c(invars-1,rep(-1,nmax-nin))
    
    z <- .C("swapf2",
            nreps=as.integer(nreps),
            nstep=as.integer(nstep),
            varx=as.double(varcov$var.x),
            covxy=as.double(varcov$cov.xy),
            vary=as.double(varcov$var.y),
            df=as.double(varcov$df),
            rparm=as.double(rparm),
            optpri=as.double(optpri),
            inpri=as.double(inpri),
 	    nrx=as.integer(nrx),
            noptuse=as.integer(noptuse),
            optblock=as.integer(optblock),
            optcur=as.integer(optcur-1),
            invars=as.integer(invars),
            nin=as.integer(nin),
            optvars=as.integer(optvars),
            nopt=as.integer(nopt),
            optmax=as.integer(optmax),
            useopt=as.integer(useopt),
            locs=as.integer(locs-1),
            pvt=as.integer(1:nmax),
            rank=integer(1),
            wrksp=double(2*nmax),
            gama=double(zmax*xmax),
            bee=double(xmax),
            xx=double(xmax^2),
            xy=double(xmax),
            zz=double(zmax^2),
            zy=double(zmax),
            xz=double(xmax*zmax),
            beta=double(xmax),
            posteriors=double(nreps*nstep),
            marg=double(nreps*nstep),
            cond=double(nreps*nstep),
            coefs=double(nreps*nstep*nmax),
            configs=integer(nreps*nstep*nmax),
            qraux=double(xmax),
            zraux=double(zmax),
            zrank=integer(1),
            tol=as.double(tol),
            postwk=double(nopt),
            coefwk=double(nopt*nmax),
            alt.marginal=double(nopt),
            alt.coef=double(nrx))

    res <- z[c("configs","posteriors","coefs","cond","marg","alt.marginal","alt.coef")]
    res$configs <- array(res$configs+1,c(nmax,nlocs,nreps))
    dim(res$coefs) <- c(nmax,nlocs,nreps)
    dim(res$marg) <- dim(res$cond) <- dim(res$posteriors) <- c(nlocs,nreps)
    res$nloc <- len.locs
    res[["call"]] <- this.call
    class(res) <- "swap"
    res
}
"twohk" <-
    function(varcov, ana.obj, ...)
{
    this.call <- match.call()
    if (class(ana.obj) != "analysis.object")
        stop("second arg must be analysis.object")
    method <- ana.obj$method
    res <-
        if (method == "F2") {
            twohkf2(varcov, ana.obj, ...)
        }
        else {
            twohkbc1(varcov, ana.obj, ...)
        }
    res$call <- this.call
    res
}
"twohkbc1"<-
    function(varcov, ana.obj, rparm=0, locs = NULL ,
             locs.prior = NULL )
{
### step thru all locations and compute the marginal and posterior
### probability for each location
    this.call <- sys.call()
    if (!missing(ana.obj)) {
        if (class(ana.obj) != "analysis.object")
            stop("arg \'ana.obj\' must be analysis object")
        if (missing(locs.prior)) {
            match.names <-
                match(dimnames(varcov$var.x)[[1]],ana.obj$reg.names,0)
            if (any(match.names==0))
                stop("variable names in varcov not found in ana.obj")
            locs.prior <-
                ana.obj$map.frame[match.names,"prior"]
        }
    }
    else {
        if (missing(locs.prior))
            locs.prior <- rep(1, ncol(varcov$var.x))
    }
    if (missing(locs)) 
        locs <- seq(ncol(varcov$var.x))
    nlocs <- length(locs)
    nin <- 1
    nopt <-  nlocs
    optmax <- 1
    nmax <- 2
    nrx <- nrow(varcov$var.x)
    optpri <- c(locs.prior)
    optvars <- locs - 1
    marg <- post <- cond <- rep(0, nlocs)
    coefs <- rep(0, nlocs)
    useopt <- rep(1,nlocs)
    if (nrx!=length(rparm)) {
        if (length(rparm)==1)
            rparm <- rep(rparm,nrx)
        else
            stop(paste("rparm must have length 1 or ",nrx))
    }
    if (nlocs != length(locs.prior))
        stop("lengths of locs.prior and locs not equal")
    
    res <- .C("twohkbc1",
              varx = as.double(varcov$var.x),
              covxy = as.double(varcov$cov.xy),
              vary = as.double(varcov$var.y),
              df = as.double(varcov$df),
              rparm = as.double(rparm),
              optpri = as.double(optpri),
              nrx=as.integer(nrx),
              optvars=as.integer(optvars) ,
              nopt = as.integer(nopt),
              optmax = as.integer(optmax),
              useopt=as.integer(useopt),
              pvt = as.integer(1:nmax),
              rank = integer(1),
              wrksp = double(2*nmax),
              gama = double(nin*optmax),
              bee = double(nin),
              xx = double(nin*nin),
              xy = double(nin),
              zz = double(optmax*optmax),
              zy = double(optmax),
              xz = double(nin*optmax),
              beta =double(optmax),
              posterior = double(nopt),
              loc.2 =as.double(marg),
              loc.1 = as.double(cond),
              coefwk = double(nmax*nopt),
              coefs.2 = as.double(coefs),
              coefs.1=as.double(coefs),
              qraux =double(nmax),
              zraux =double(nmax),
              zrank=integer(1),
              tol=as.double(1e-10))[c("loc.2","loc.1","coefs.2","coefs.1")]
    dim(res$loc.2) <- c(nlocs,1)
    dim(res$loc.1) <- c(nlocs,1)
    dim(res$coefs.2) <- c(1,nlocs)
    dim(res$coefs.1) <- c(1,nlocs)
    res[["call"]] <- this.call
    res
}
"twohkf2"<-
    function(varcov, ana.obj, rparm=0, locs = NULL , locs.prior = NULL ,
             combo.prior = rep(1, 3)/3)
{
### step thru all locations and compute the marginal and posterior
### probability for each location
    this.call <- sys.call()
    if (!missing(ana.obj)) {
        if (class(ana.obj) != "analysis.object")
            stop("arg \'ana.obj\' must be analysis object")
        if (missing(locs.prior)) {
            match.names <-
                match(dimnames(varcov$var.x)[[1]],ana.obj$reg.names,0)
            match.cols <- unique(col(ana.obj$reg.names)[match.names])
            if (any(match.names==0))
                stop("variable names in varcov not found in ana.obj")
            locs.prior <-
                ana.obj$map.frame[match.cols,"prior"]
        }
    }
    else {
        if (missing(locs.prior))
            locs.prior <- rep(1, ncol(varcov$var.x)/2)
    }

    if (missing(locs))
        locs <- rep(seq(ncol(varcov$var.x)/2), rep(2, ncol(varcov$var.x)/2))
    
    uniq.locs <- unique(locs)
    nlocs <- length(uniq.locs)
    nin <- 2
    nopt <-  nlocs*3
    optmax <- 2
    nmax <- 4
    nrx <- nrow(varcov$var.x)
    if (length(combo.prior)!=3) stop("combo.prior must = 3")
    if (nrx!=length(rparm)) {
        if (length(rparm)<3)
            rparm <- rep(rparm,nrx/length(rparm))
        else
            stop(paste("rparm must have length 1, 2, or ",nrx))
    }
    if (length(locs.prior)!=nlocs)
        stop("length(locs.prior)!=length(unique(locs))")
    
    optpri <- c(locs.prior)%o%c(combo.prior)
    optvars <- matrix(c(rbind(seq(from=1,by=2,length=nlocs),0),
                        rbind(seq(from=2,by=2,length=nlocs),0),
                        seq(from=1,by=1,length=2*nlocs))-1,nrow=2)
    marg <- post <- cond <- rep(0, nlocs*3)
    coefs <- rep(0, 2*nlocs)
    useopt <- rep(1,nlocs*3)
    
    res <- .C("twohkf2",
              varx = as.double(varcov$var.x),
              covxy = as.double(varcov$cov.xy),
              vary = as.double(varcov$var.y),
              df = as.double(varcov$df),
              rparm = as.double(rparm),
              optpri = as.double(optpri),
              nrx=as.integer(nrx),
              optvars=as.integer(optvars) ,
              nopt = as.integer(nopt),
              optmax = as.integer(optmax),
              useopt=as.integer(useopt),
              pvt = as.integer(1:nmax),
              rank = integer(1),
              wrksp = double(2*nmax),
              gama = double(nin*optmax),
              bee = double(nin),
              xx = double(nin*nin),
              xy = double(nin),
              zz = double(optmax*optmax),
              zy = double(optmax),
              xz = double(nin*optmax),
              beta =double(optmax),
              posterior = double(nopt),
              loc.2 =as.double(marg),
              loc.1 = as.double(cond),
              coefwk = double(nmax*nopt),
              coefs.2 = as.double(coefs),
              coefs.1=as.double(coefs),
              qraux =double(nmax),
              zraux =double(nmax),
              zrank=integer(1),
              tol=as.double(1e-10))[c("loc.2","loc.1","coefs.2","coefs.1")]
    dim(res$loc.2) <- c(nlocs,3)
    dim(res$loc.1) <- c(nlocs,3)
    dim(res$coefs.2) <- c(2,nlocs)
    dim(res$coefs.1) <- c(2,nlocs)
    res[["call"]] <- this.call
    res
}
"unique.config"<-
    function(swap.obj)
{
    configs <- apply(swap.obj$config, 2:3, function(x)
                     paste(sort(x), collapse = "."))
    unique.con <- unique(configs)
    conf.match <- match(configs, unique.con)
    uniq.mat <-
        matrix(swap.obj$config, nr =
               dim(swap.obj$config)[1])[, ! duplicated(conf.match), drop = F]
    list(uniq = uniq.mat, match = conf.match)
}

"update.bqtl"<-
    function(setup,loc.mat,ana.obj)
{
    dmloc <- dim(loc.mat)
    if (length(dmloc) == 0){
        nloc <-   length(loc.mat)
        n.alt <- 1
    }
    else {
        nloc <- dmloc[1]
        n.alt <- dmloc[2]
    }
    if ( nloc != setup$nparm[4] )
        stop("nrow(loc.mat) not = number of loci in setup")
    if (any( loc.mat > ncol(ana.obj$loc.right))||any(loc.mat<1) )
        stop("invalid locus number")
    subset <- attr(setup,"subset")
    lambda <- switch(ana.obj$method,
                   RI.self = {ana.obj$map.frame$lambda/(2 - ana.obj$map.frame$lambda)},
                   RI.sib = {ana.obj$map.frame$lambda/(4 - 3 * ana.obj$map.frame$lambda)},
                   ana.obj$map.frame$lambda)
    setup <-
        c(list("upbqtl"),
          setup[-c(1,22)],
          list( loc.mat = as.integer(loc.mat-1)) ,
          list( n.alt = as.integer(n.alt)) ,
          list( loc.order = as.integer(apply(as.matrix(loc.mat),2,order)-1)) ,
          list( res = double(n.alt)),
          list( orig.x = as.double(setup$x)) ,
          list( loc.right = as.integer(ana.obj$loc.right[subset,]-1)) ,
          list( map.lambda = as.double(lambda)) ,
          list( state.matrix = as.double(ana.obj$state.matrix[subset,,])) ,
          list( n.state.loc = as.integer(dim(ana.obj$state.matrix)[2])))
    
    res <- do.call(".C",setup)
    
    res$res
}

"varcov"<-
    function(x,ana.obj,partial=NULL,scope=ana.obj$reg.names,...)
{
    this.call <- match.call(expand=TRUE)
    using.R <- exists("is.R")&&is.R()
    bqtl.specials <- c("configs","locus") # dont use: add or dom so
                                        # swap, hk wont jam up    
    reg.terms <- terms(x,specials=bqtl.specials)
    if (any(attr(reg.terms,"order")>1))
        stop("cannot use interactions in this formula")

    covs <-
        rhs.bqtl(reg.terms,ana.obj,bqtl.specials,NULL,scope,
                 expand.specials=FALSE)
    response <-
        dimnames(attr(reg.terms,"factors"))[[1]][attr(reg.terms,"response")]
    if (!is.null(partial)) {
        p.terms <- terms(partial,specials=bqtl.specials)
        if (any(attr(p.terms,"order")>1))
            stop("cannot use interactions in this formula")
        cntl <-
            rhs.bqtl(p.terms,ana.obj,bqtl.specials,NULL,scope,
                     expand.specials=FALSE)
        dummy.formula <- eval(parse(text=paste(response,"~",covs,"+",cntl)))
    }
    else {
        cntl <- NULL
        dummy.formula <- eval(parse(text=paste(response,"~",covs)))
    }
    dummy.terms <- terms(dummy.formula)
    mf <- model.frame(dummy.terms,ana.obj$data, na.action=na.omit)
    if (is.null(partial)) {
        x <- model.matrix(dummy.terms,mf)[,-1]
        y <- model.extract(mf,"response")
        make.varcov(x,y,...)
    }
    else
    {
        
        cntl.out <- eval(parse(text=paste("~ . -", cntl)))
        x.formula <- update(dummy.formula,
                            eval(parse(text=paste("~ . -", cntl)))  )
        x <- model.matrix(terms(x.formula),mf)[,-1]
        y <- model.extract(mf,"response")
        z <- model.matrix(terms(eval(parse(text=paste("~", cntl)))),mf)[,-1,drop=F]
        y[] <- lsfit(z,y)$resid
        x[] <- lsfit(z,x)$resid
        res <- make.varcov(x,y,...)
        r <- ncol(z)
        res$df <- res$df - r
        res$call <- this.call
        ## further adjustments unneeded since posterior is scaled by variance
        res
    }
    
}
version.bqtl<-"Version:1.0-3"
"zero.dup"<-
function(x, dig = 6)
{
	ifelse(duplicated(signif(x, dig)), 0, x)
}
###
### The obligatory zzz file
###
.First.lib <- function(lib, pkg) {
  library.dynam( "bqtl", pkg, lib )
}

