######################################################
# This function takes data, pos (marker positions), and
# nlocus (the grid size for test suppressor gene loci)
#
# A sequence, lsupp, between .001 and .999 of size nlocus is first generated.
# After fixing the position of x_s in lsupp, log-likelihood function is
# maximized over delta, lambda, and omega such that omega>=delta.
# The value of maximized log-likelihood is stored as a vector (loglik).
# Other parameter values where log-likelihood is maximized
# are also temporarily stored.
#
# Max(loglik) is found and the corresponding parameter values, delta
# lambda, omega were returned with lsupp, and loglik.

mle.alt<-function(data, pos, nlocus) {
  lsupp <- seq(0.001, 0.999,length = nlocus)
  loglik <- rep(NA, nlocus)
  messages <- rep(NA, nlocus)
  theta <- matrix(NA, nlocus, 3)
  theta.start <- c(0.5, 10, 0)
  
  for(i in 1:nlocus) {
    print(i)
    profile <- nlminb(start = theta.start, mloglik.theta.alt, locus=lsupp[i],
		      data = data, pos = pos, lower = c(0.001, 1,  -Inf), 
		      upper = c(0.999, Inf, Inf))
                        #see below for the definition of mloglik.theta.alt
                        #leta=log(1/lambda)

    
    loglik[i] <-  - profile$objective	#log-likelihood
    theta[i,  ] <- profile$parameters	#optimized theta value
    messages[i] <- profile$message	
  }

  #compute (delta, omega, lambda) from theta
  #(see below for the definition of theta)
  
  deltahat <- theta[, 1]
  omegahat <- theta[, 1] + (1 - theta[, 1])/theta[, 2]
  lambdahat <- 1/exp(theta[, 3])

  locus <- lsupp
  loglik.max <- max(loglik)
  index <- (loglik == loglik.max)

  deltahat.opt <- deltahat[index]
  omegahat.opt <- omegahat[index]
  lambdahat.opt <- lambdahat[index]
  locushat.opt <- locus[index]

  print(list(loglik.max = loglik.max, deltahat.opt = deltahat.opt, 
	     omegahat.opt = omegahat.opt, lambdahat.opt = lambdahat.opt, 
	     locushat.opt = locushat.opt))
  return(list(loglik.max = loglik.max, deltahat.opt = deltahat.opt, 
	      omegahat.opt = omegahat.opt, lambdahat.opt = lambdahat.opt, 
	      locushat.opt = locushat.opt, loglik = loglik, lsupp = lsupp, 
	      messages = messages))
}

######################################################
# This function computes (-(log-likelihood)) under the alternative hypo,
# which is exactly same as what mloglik.alt does.
# It is just modified to be used as an object function of nlminb.

mloglik.theta.alt<-function(theta, locus, data, pos) {
        #re-parametrization   
        #theta[1]=delta, theta[2]= (1-delta)/(omega-delta), theta[3]=leta
	delta <- theta[1]
	omega <- theta[1] + (1 - theta[1])/theta[2]
	leta <- theta[3]
	return(mloglik.alt(delta, omega, leta, locus, data, pos))
}

######### Examples ##################################
pro<-dget("pro.data")
thy<-dget("thy.data")
fjp<-dget("fjp.data")

pro.mle.alt<-mle.alt(data=pro$data, pos=pro$pos, nlocus=25)
thy.mle.alt<-mle.alt(data=thy$data, pos=thy$pos, nlocus=25)
fjp.mle.alt<-mle.alt(data=fjp$data, pos=fjp$pos, nlocus=25)

######################################################