.packageName <- "qtlDesign"
"deflate.bc" <-
function(theta)
  {
    theta1 <- recomb(0.5*genetic.dist(theta))
    q <- ((1-theta1)^2)/(theta1^2 + (1-theta1)^2)
    A <- (1-4*q*(1-q))*(1-theta)
    A
  }

"deflate.f2" <-
function(theta)
  {
    t <- recomb(0.5*genetic.dist(theta))
    num <- 6*t^4 - 12*t^3 +10*t^2 - 4*t + 1
    den <- 8*t^4 - 16*t^3 +12*t^2 - 4*t + 1
    add.factor <- deflate.bc(theta)
    dom.factor <- (add.factor^2)*num/den
    list(add=add.factor,dom=dom.factor)
  }

"fracmiss" <-
function(y,delta,m=NULL,theta=NULL)
  {
    if(is.null(m)&is.null(theta))
      {
        ans <- fracmiss0(y,delta)
      }
    else
      {
        if( (length(m)==1) & (length(theta)==1) )
          {
            ans <- fracmiss1(y,m,delta,theta)
          }
        else
          {
            if( (length(m)==2) & (length(theta)==2) )
              {
                ans <- fracmiss2(y,m[1],m[2],delta,theta[1],theta[2])
              }
            else
              {
                stop("Incompatible m and theta lengths.")
              }
          }
      }
    ans
  }

"fracmiss0" <-
function(y,delta)
{
  density <- 0.5 * ( dnorm(y,mean=delta) + dnorm(y,mean=-delta) )
  A <- exp(delta*y)
  B <- exp(-delta*y)
  genovar <- (A*B)/((A+B)^2)
  # print(c(density,genovar))
  4*y*y*density*genovar
}

"fracmiss1" <-
function(y,m,delta,theta)
{
  if(m==0)
    {
      density <- 0.5 * ( theta * dnorm(y,mean=delta)
                        + (1-theta) * dnorm(y,mean=-delta) )
      A <- exp(delta*y) * theta
      B <- exp(-delta*y) * (1-theta)
    }
  if(m==1)
    {
      density <- 0.5 * ( (1-theta) * dnorm(y,mean=delta)
                        + theta * dnorm(y,mean=-delta) )
      A <- exp(delta*y) * (1-theta)
      B <- exp(-delta*y) * theta
    }

  genovar <- (A*B)/((A+B)^2)
  # print(c(density,genovar))
  4*y*y*density*genovar
}

"fracmiss2" <-
function(y,m1,m2,delta,theta1,theta2)
{
  # y = phenotype
  # m1 = left marker genotype
  # m2 = right marker genotype
  # delta = QTL effect; means are -delta and +delta
  # theta1 = recombination fraction to left marker
  # theta2 = recombination fraction to right marker

  lik0.pheno <- dnorm(y,mean=-delta)
  lik0.m1 <- theta1*m1 + (1-theta1)*(1-m1)
  lik0.m2 <- theta2*m2 + (1-theta2)*(1-m2)      
  lik0 <- lik0.pheno * lik0.m1 * lik0.m2 * 0.5

  lik1.pheno <- dnorm(y,mean=delta)
  lik1.m1 <- (1-theta1)*m1 + theta1*(1-m1)
  lik1.m2 <- (1-theta2)*m2 + theta2*(1-m2)      
  lik1 <- lik1.pheno * lik1.m1 * lik1.m2 * 0.5

  density <- lik0 + lik1
  genovar <- (lik0*lik1)/((lik0+lik1)^2)
  4*y*y*density*genovar
}

"genetic.dist" <-
function(theta)
  {
    -0.5*log(1-2*theta)
  }

"info.null" <-
function(alpha)
  {
    z <- -qnorm(alpha/2)
    2*( z*dnorm(z) + alpha/2 )
  }

"info.bc.null" <-
function(alpha,theta=0)
  {
    info.null(alpha)*deflate.bc(theta)
  }

"info.f2.null" <-
function(alpha,theta=0)
  {
    defl <- deflate.f2(theta)
    list( add=info.null(alpha)*defl$add, dom=info.null(alpha)*defl$dom )
  }

"info2cost.bc.null" <-
function(alpha,cost,d=0,G=NULL)
  {
    if((d==0) & is.null(G))
      {
        ans <- info.bc.null(alpha,theta=0)/(1+cost*alpha)
      }
    else
      {
        if((d==0)|is.null(G)|(G<=0))
          {
            stop("Cannot compute with given d and G.")
          }
        else
          {
            theta <- recomb(d)
            ans <- info.bc.null(alpha,theta)/(1+cost*alpha*G/d)
          }
      }
    ans
  }

"missinfo" <-
function(delta,alpha,theta=NULL)
  {
    if(is.null(theta))
      {
        ans <- missinfo0(delta,alpha)
      }
    else
      if(length(theta)==1)
        {
          ans <- missinfo1(delta,alpha,theta)
        }
      else
        {
          ans <- missinfo2(delta,alpha,theta[1],theta[2])
        }
    mean(ans$mi)
  }

"missinfo.sim" <-
function(delta,n,alpha,theta=NULL)
  {
    if(is.null(theta))
      {
        ans <- missinfo0.sim(delta,n,alpha)
      }
    else
      if(length(theta)==1)
        {
          ans <- missinfo1.sim(delta,n,alpha,theta)
        }
      else
        {
          ans <- missinfo2.sim(delta,n,alpha,theta[1],theta[2])
        }
    mean(ans$mi)
  }

"missinfo0" <-
function(delta,alpha)
{
  limit <- uniroot(pmixnorm,interval=c(0,delta+5),mean=c(-delta,delta),
               level=1-alpha/2)$root
  mi <- 2*integrate(fracmiss0,rel.tol=1e-7,lower=0,upper=limit,
                    delta=delta)$value
  list( mi=mi, lim=limit )
}

"missinfo0.sim" <-
function(delta,n,alpha)
{
y <- rnorm(n=n)
g <- rbinom(n=n,size=1,prob=0.5)
y <- y+(2*g-1)*delta

srt <- order(y)

y <- y[srt]
g <- g[srt]

m <- round(n*alpha/2)

if(m>=1)
	qq <- c(g[1:m], rep(0.5,n-2*m), g[(n-m+1):n])
else
	qq <- rep(0.5,n)

a <- qq*dnorm(y,mean=delta)
b <- (1-qq)*dnorm(y,mean=-delta)
qstar <- a/(a+b)

ans <- 4*y*y*qstar*(1-qstar)
list(mi=ans,y=y,g=g,delta=delta,n=n,alpha=alpha,qstar=qstar,
     ulim=y[n-m+1],llim=y[m])
}

"missinfo1" <-
function(delta,alpha,theta)
{
  limit <- uniroot(pmixnorm,interval=c(0,delta+5),mean=c(-delta,delta),
               level=1-alpha/2)$root
  mi.mid <- 2*integrate(fracmiss0,sub=10000,rel.tol=1e-7,lower=0,upper=limit,
                         delta=delta)$value
  mi0.upp <- 2*integrate(fracmiss1,sub=10000,rel.tol=1e-7,lower=limit,
                         upper=delta+10,
                         m=0,delta=delta,theta=theta)$value
  mi1.upp <- 2*integrate(fracmiss1,sub=10000,rel.tol=1e-7,lower=limit,
                         upper=delta+10,
                         m=1,delta=delta,theta=theta)$value

  list( mi=mi.mid+mi0.upp+mi1.upp, lim=limit )
}

"missinfo1.sim" <-
function(delta,n,alpha,theta)
{
y <- rnorm(n=n)
# simulate the QTL
g <- rbinom(n=n,size=1,prob=0.5)
# simulate the phenotype
y <- y+(2*g-1)*delta
# simulate a marker theta recombination fraction away
g1 <- rbinom(n=n,size=1,prob=theta)
g1 <- (g+g1) %% 2

srt <- order(y)

y <- y[srt]
g <- g[srt]
g1 <- g1[srt]

m <- round(n*alpha/2)

if(m>=1)
	qq <- c(g1[1:m]*(1-theta)+theta*(1-g1[1:m]),
                rep(0.5,n-2*m),
                g1[(n-m+1):n]*(1-theta)+theta*(1-g1[(n-m+1):n]))
else
	qq <- rep(0.5,n)

a <- qq*dnorm(y,mean=delta)
b <- (1-qq)*dnorm(y,mean=-delta)
qstar <- a/(a+b)

ans <- 4*y*y*qstar*(1-qstar)
list(mi=ans,y=y,g=g,delta=delta,n=n,alpha=alpha,theta=theta,qstar=qstar)
}

"missinfo2" <-
function(delta,alpha,theta1,theta2)
{
  limit <- uniroot(pmixnorm,interval=c(0,delta+5),mean=c(-delta,delta),
               level=1-alpha/2)$root
  mi.mid <- 2*integrate(fracmiss0,sub=10000,rel.tol=1e-7,lower=0,upper=limit,
                         delta=delta)$value

  mi00.upp <- 2*integrate(fracmiss2,sub=10000,rel.tol=1e-7,lower=limit,
                          upper=delta+10,m1=0,m2=0,delta=delta,
                          theta1=theta1,theta2=theta2)$value

  mi01.upp <- 2*integrate(fracmiss2,sub=10000,rel.tol=1e-7,lower=limit,
                          upper=delta+10,m1=0,m2=1,delta=delta,
                          theta1=theta1,theta2=theta2)$value

  mi10.upp <- 2*integrate(fracmiss2,sub=10000,rel.tol=1e-7,lower=limit,
                          upper=delta+10,m1=1,m2=0,delta=delta,
                          theta1=theta1,theta2=theta2)$value

  mi11.upp <- 2*integrate(fracmiss2,sub=10000,rel.tol=1e-7,lower=limit,
                          upper=delta+10,m1=1,m2=1,delta=delta,
                          theta1=theta1,theta2=theta2)$value


  list( mi=mi.mid+mi00.upp+mi01.upp+mi10.upp+mi11.upp, lim=limit )
}

"missinfo2.sim" <-
function(delta,n,alpha,theta1,theta2)
{
y <- rnorm(n=n)
# simulate the QTL
g <- rbinom(n=n,size=1,prob=0.5)
# simulate the phenotype
y <- y+(2*g-1)*delta
# simulate a marker theta1 recombination fraction away
g1 <- rbinom(n=n,size=1,prob=theta1)
g1 <- (g+g1) %% 2
# simulate a marker theta2 recombination fraction away
g2 <- rbinom(n=n,size=1,prob=theta2)
g2 <- (g+g2) %% 2

# find order of phenotypes
srt <- order(y)
# sort the phenotypes, QTL genotypes, and flanking marker genotypes
y <- y[srt]
g <- g[srt]
g1 <- g1[srt]
g2 <- g2[srt]

# calculate prior probabilities of QTL=1 given the flanking markers
# prob QTL=1
qq.1 <- g1*(1-theta1) + (1-g1)*theta1 + g2*(1-theta2) + (1-g2)*theta2
# prob QTL=0
qq.0 <- (1-g1)*(1-theta1) + g1*theta1 + (1-g2)*(1-theta2) + g2*theta2
# prob QTL=1 given flanking markers
qq <- qq.1/(qq.0+qq.1)

# number genotyped at each flank
m <- round(n*alpha/2)

# if number genotyped is greater than 2, get the prior probs
# else replace it by the equilibrium distribution
if(m>=1)
	qq <- c(qq[1:m], rep(0.5,n-2*m), qq[(n-m+1):n])
else
	qq <- rep(0.5,n)

# likelihood for y given QTL=1
a <- qq*dnorm(y,mean=delta)
# likelihood for y given QTL=0
b <- (1-qq)*dnorm(y,mean=-delta)
# posterior prob of QTL=1
qstar <- a/(a+b)

ans <- 4*y*y*qstar*(1-qstar)
list(mi=ans,y=y,g=g,delta=delta,n=n,alpha=alpha,theta=c(theta1,theta2),
     qstar=qstar)
}

"optalpha.bc" <-
function(cost,d=0,G=NULL)
  {
    optimize(f=info2cost.bc.null,interval=c(0.0001,0.9999),max=TRUE,
             G=G,d=d,cost=cost)$maximum
  }

"pmixnorm" <-
function(x,mean=c(0,0),sd=c(1,1),mix.prop=0.5,level=0)
{
  ans <- mix.prop * pnorm(x,mean[1],sd[1]) +
    (1-mix.prop) * pnorm(x,mean[2],sd[2]) - level
  ans
}

"power.bc" <- function(n,prop,thresh=3,alpha=1,theta=0,effective.n=FALSE)
  {
    # convert to deviation from overall mean
    delta <- prop2delta.bc(prop)
    # effective sample size
    m <- n * info.bc.null(alpha,theta)
    # non-centrality parameter
    ncp <- m*delta^2
    if( m<=30 )
      {
        stop("Approximation not reliable as effective sample size < 30.")
      }
    # threshold in 2*loglikelihood units
    T <- 2*log(10)*thresh
    # power using non-central chi-square
    pow <- 1-pchisq(T,df=1,ncp=ncp)
    if(!effective.n)
     {
       return(pow)
     }
    else
      {
        return(list(power=pow,effective.n=m))
      }
  }

"detectable.bc" <- function (n, power = 0.8, thresh = 3,
                           alpha = 1, theta = 0, delta = FALSE) 
{
  # proportion of variance explained for given sample size,
  # power, threshold, selection fraction, and size of marker interval
    prop <- uniroot(function(x) {
        power.bc(n, x, thresh, alpha, theta, effective.n = FALSE) - 
            power
    }, interval = c(0, (1000/n)/(1+1000/n)))$root
    
  # decide what to return depending on delta flag
    if (!delta) {
        return(prop)
    }
    else {
        return(prop2delta.bc(prop))
    }
}


"prop2delta.bc" <- function(prop)
  {
    sqrt(1/(1-prop)-1)
  }

"delta2prop.bc" <- function(delta)
  {
    delta^2/(1+delta^2)
  }
"power.f2" <-  function (n, model, thresh = 3, alpha = 1,
                       theta = 0, effective.n = FALSE) 
{
  # get info per individual
  iii <- info.f2.null(alpha, theta)
  # additive and dominance components
  a <- model[1]
  d <- model[2]
  # calculate non-centrality parameter
  ncp <- n * ( iii$add*a^2/2 + iii$dom*d^2/4 )
  m <- n*min(iii$add,iii$dom)
  # if effective sample size not big enough, stop
  if (m <= 30) {
    stop("Approximation not reliable as effective sample size < 30.")
  }
  # calculate threshold in chi-square scale
  T <- 2 * log(10) * thresh
  # calculate power
  pow <- 1 - pchisq(T, df = 2, ncp = ncp)
  # decide what to return depending on effective.n flag
  if (!effective.n) {
    return(pow)
  }
  else {
    return(list(power = pow, effective.n = m))
  }
}


"detectable.f2" <- function (n, model="add", power = 0.8, thresh = 3,
                           alpha = 1, theta = 0, delta=FALSE)
{
  # model
  if(model=="add")
    {
      a <- 1
      d <- 0
      model <- c(a,d)
    }
  else if(model=="dom")
    {
      a <- 1
      d <- 1
      model <- c(a,d)
    }
  else if( is.numeric(model) && (length(model)==2) )
    {
      model <- model
    }
  else
    {
      stop("Cannot understand model argument.")
    }

  # proportion of variance explained for given sample size,
  # power, threshold, selection fraction, and size of marker interval
  del <- uniroot(function(x) {
    power.f2(n, x*model, thresh, alpha, theta, effective.n = FALSE) - 
      power  }, interval = c(0, (1000/n)/(1+1000/n)))$root

    # decide what to return depending on delta flag
    if (!delta)
      {
        return(delta2prop.f2(del,model))
      }
    else
      {
        return(del*model)
      }
  
}

# from proportion of variance explained to a delta relative to a model
"prop2delta.f2" <- function(prop,model)
  {
    # additive and dominance components
    a <- model[1]
    d <- model[2]
    # get genetic variance
    gv <- prop/(1-prop)
    # convert to delta
    delta <- sqrt(gv/(d^2/4+a^2/2))
    delta <- delta*model
    delta
  }

"delta2prop.f2" <- function(delta,model)
  {
    # additive and dominance components
    a <- delta*model[1]
    d <- delta*model[2]
    gv <- (d^2/4+a^2/2)
    gv/(1+gv)
  }

# proportion of variance explained to genetic variance 
"prop2gv" <- function (prop) 
{
    prop/(1 - prop)
}
"gv2prop" <- function (gv) 
{
    gv/(1 + gv)
}

"recomb" <-
function(d)
  {
    0.5*(1-exp(-2*d))
  }

