#######################################################
#
# Profile LOD plotting code: 
#
#   Make sure a graphics device is active.
#
########################################################
# Data
pro<-dget("pro.data")
thy<-dget("thy.data")
fjp<-dget("fjp.data")

########### 1. Based on mles ###########################
# MLE calculations

# See 'mle.null.s'
pro.mle.null<-mle.null(data=pro$data, pos=pro$pos)
thy.mle.null<-mle.null(data=thy$data, pos=thy$pos)
fjp.mle.null<-mle.null(data=fjp$data, pos=fjp$pos)

# See 'mle.alt.s'
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)


plot(pro.mle.alt$lsupp,
     (pro.mle.alt$loglik-pro.mle.null$loglik)/log(10),
     type="n",xlab="",ylab="LOD",xlim=c(0,1),
     xaxs="i",yaxs="i",axes=F)
axis(side=1) 
axis(side=2) 
lines(pro.mle.alt$lsupp,
     (pro.mle.alt$loglik-pro.mle.null$loglik)/log(10),type='l')


########### 2. Based on simple estimators  ##############

# See 'simple.est.null.s'

pro.simple.est.null<-simple.est.null(data=pro$data, pos=pro$pos)
thy.simple.est.null<-simple.est.null(data=thy$data, pos=thy$pos)
fjp.simple.est.null<-simple.est.null(data=fjp$data, pos=fjp$pos)

# See 'simple.est.alt.s'

pro.simple.est.alt<-simple.est.alt(lambda= pro.simple.est.null$lambdahat,
                                   data=pro$data, pos=pro$pos, nlocus=25)

thy.simple.est.alt<-simple.est.alt(lambda= thy.simple.est.null$lambdahat,
                                   data=thy$data, pos=thy$pos, nlocus=25)

fjp.simple.est.alt<-simple.est.alt(lambda= fjp.simple.est.null$lambdahat,
                                   data=fjp$data, pos=fjp$pos, nlocus=25)


plot(pro.simple.est.alt$lsupp,
     (pro.simple.est.alt$loglik-pro.simple.est.null$loglik)/log(10),
     type="n",xlab="",ylab="LOD",xlim=c(0,1),
     xaxs="i",yaxs="i",axes=F)
axis(side=1) 
axis(side=2) 
lines(pro.simple.est.alt$lsupp,
     (pro.simple.est.alt$loglik-pro.simple.est.null$loglik)/log(10),type='l')

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