.packageName <- "phyclust"
#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/b_seq.r"
### This file contains functions for PARAMETERIC BOOTSTRAP that generate
### star trees based on an pcobj fitted by phyclust(), and generate sequences
### according to the star trees.

bootstrap.seq <- function(pcobj, star.trees){
  if(pcobj$code.type == "NUCLEOTIDE"){
    seq.boot <- bootstrap.seq.nucleotide(pcobj, star.trees)
  } else if(pcobj$code.type == "SNP"){
    seq.boot <- bootstrap.seq.snp(pcobj, star.trees)
  } else{
    stop("code.type is not found.")
  }
  seq.boot
} # End of bootstrap.seq()

bootstrap.seq.nucleotide <- function(pcobj, star.trees){
  if(pcobj$code.type != "NUCLEOTIDE"){
    stop("Only for NUCLEOTIDE data.")
  }

  if(is.null(pcobj$QA$kappa)){
    kappa <- 1
  } else{
    kappa <- pcobj$QA$kappa
  }
  if(length(kappa) == 1){		# EE or EV
    kappa <- rep(kappa, pcobj$K)
  }

  if(is.null(pcobj$QA$pi)){
    pi <- rep(0.25, 4)
  } else{
    pi <- pcobj$QA$pi
  }
  if(length(pi) == 4){			# EE or EV
    pi <- matrix(rep(pi, pcobj$K), nrow = pcobj$K, byrow = TRUE)
  }

  seq.nucleo <- NULL
  for(k in 1:pcobj$K){
    seq.nucleo[[k]] <- gen.seq.HKY(star.trees[[k]], pi[k,], kappa[k],
                                   pcobj$L, anc.seq = pcobj$Mu[k,])
    if(star.trees[[k]]$n.tip == 1){
      seq.nucleo[[k]] <- seq.nucleo[[k]][1:2]
      tmp <- unstrsplit(seq.nucleo[[k]][1], " ")
      tmp <- tmp[tmp != ""]
      seq.nucleo[[k]][1] <- paste(" ", as.numeric(tmp[1]) - 1,
                                  " ", tmp[-1], sep = "")
      class(seq.nucleo[[k]]) <- "seqgen"
    }
  }

  seq.nucleo
} # End of bootstrap.seq.nucleotide().

bootstrap.seq.snp <- function(pcobj, star.trees){
  if(pcobj$code.type != "SNP"){
    stop("Only for SNP data.")
  }

  if(is.null(pcobj$QA$pi)){
    pi <- rep(0.5, 2)
  } else{
    pi <- pcobj$QA$pi
  }
  if(length(pi) == 2){			# EE or EV
    pi <- matrix(rep(pi, pcobj$K), nrow = pcobj$K, byrow = TRUE)
  }

  seq.snp <- NULL
  for(k in 1:pcobj$K){
    seq.snp[[k]] <- gen.seq.SNP(star.trees[[k]], pi[k,],
                                pcobj$L, anc.seq = pcobj$Mu[k,])
    if(star.trees[[k]]$n.tip == 1){
      seq.snp[[k]] <- seq.snp[[k]][1:2]
      tmp <- unstrsplit(seq.snp[[k]][1], " ")
      tmp <- tmp[tmp != ""]
      seq.snp[[k]][1] <- paste(" ", as.numeric(tmp[1]) - 1,
                               " ", tmp[-1], sep = "")
      class(seq.snp[[k]]) <- "seqgen"
    }
  }

  seq.snp
} # End of bootstrap.seq.snp().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/b_star_tree.r"
### This file contains functions for PARAMETERIC BOOTSTRAP that generate
### star trees based on an pcobj fitted by phyclust(), and generate sequences
### according to the star trees.

bootstrap.star.tree <- function(n.tip, height = NULL){
  if(n.tip > 1){
    ms.tree <- ms(n.tip, opts = "-T")
    rooted.tree <- read.tree(text = ms.tree[3])
    star.tree <- as.star.tree(rooted.tree)
  } else if(n.tip == 1){
    star.tree <- list(edge = matrix(c(3, 3, 1, 2), nrow = 2),
                      Nnode = 1,
                      tip.label = c("1", "2"),
                      edge.length = c(1, 0))
    class(star.tree) <- "phylo"
  } else{
    stop("n.tip > 0")
  }

  if(! is.null(height)){
    star.tree$edge.length[star.tree$edge[, 2] <= n.tip] <- height
  }
  star.tree$n.tip <- n.tip

  star.tree
} # End of bootstrap.star.tree().

bootstrap.star.trees <- function(pcobj, min.n.class = 1){
  if(min.n.class * pcobj$K >= pcobj$N.X.org){
    stop("min.n.class * pcobj$k >= pcobj$N.X.org")
  }

  repeat{
    n.class <- rmultinom(1, pcobj$N.X.org, pcobj$Eta)
    if(all(n.class >= min.n.class)){
      n.class <- as.vector(n.class)
      break
    }
  }

  Tt <- pcobj$QA$Tt
  if(length(Tt) == 1){		# EE or VE
    Tt <- rep(Tt, pcobj$K)
  }

  star.trees <- NULL
  for(k in 1:pcobj$K){
    star.trees[[k]] <- bootstrap.star.tree(n.class[k], height = Tt[k])
  }

  star.trees
} # End of bootstrap.star.trees().

bootstrap.merge.seq <- function(seqs, code.type = .code.type[1]){
  K <- length(seqs)
  da <- read.seqgen(seqs[[1]], code.type = code.type[1])
  if(K > 1){
    for(k in 2:K){
      tmp.da <- read.seqgen(seqs[[k]], code.type = code.type[1])
      da$nseq <- da$nseq + tmp.da$nseq
      da$seqname <- c(da$seqname, tmp.da$seqname)
      da$org.code <- rbind(da$org.code, tmp.da$org.code)
      da$org <- rbind(da$org, tmp.da$org)
    }
    da$info <- paste(" ", da$nseq, " ", da$seqlen, sep = "")
  }
  da
} # End of bootstrap.merge.seq().

bootstrap.star.trees.seq <- function(pcobj, min.n.class = 1){
  star.trees <- bootstrap.star.trees(pcobj, min.n.class = min.n.class)
  seq.boot <- bootstrap.seq(pcobj, star.trees)
  list(trees = star.trees, seq = seq.boot)
} # End of bootstrap.star.trees.seq().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_base.r"
### Object Name (Type): Definaiton (Examples)
###   codeseq (string): A nucleotide sequence. (ACGT...)
###   nidseq (vector[L]): A nucleotide id sequence. (1234....)
###                       The length should be L = 2 + 3 * n + 3 for n >= 1.


### Transfer an nucleotide sequence to a code id sequence.
### Eg. ACGT-... => 01234...
code2nid <- function(codeseq){
  nidseq <- codeseq
  for(i in 1:nrow(.nucleotide)){
    nidseq[codeseq == as.character(.nucleotide$code[i])] <- .nucleotide$nid[i]
    nidseq[codeseq == as.character(.nucleotide$code.l[i])] <- .nucleotide$nid[i]
  }
  nidseq[! codeseq %in%
         c(as.character(.nucleotide$code), as.character(.nucleotide$code.l))] <- 4
  nidseq <- as.numeric(nidseq)

  if(is.matrix(codeseq)){
    dim(nidseq) <- dim(codeseq)
  }
  nidseq
} # End of code2nid().

### Transfer a nucleotide id sequence to an nucleotide sequence.
### Eg. 01234... => ACGT-...
nid2code <- function(nidseq, lower.case = TRUE){
  codeseq <- nidseq
  for(i in 1:nrow(.nucleotide)){
    if(lower.case){
      codeseq[nidseq == .nucleotide$nid[i]] <- as.character(.nucleotide$code.l[i])
    } else{
      codeseq[nidseq == .nucleotide$nid[i]] <- as.character(.nucleotide$code[i])
    }
  }
  codeseq[!nidseq %in% .nucleotide$nid] <- "-"

  if(is.matrix(nidseq)){
    dim(codeseq) <- dim(nidseq)
  }
  codeseq
} # End of nid2code().

### Transfer an SNP sequence to a SNP id sequence.
### Eg. 1212-... => 01012...
snp2sid <- function(snpseq){
  sidseq <- snpseq
  for(i in 1:nrow(.snp)){
    sidseq[snpseq == as.character(.snp$code[i])] <- .snp$sid[i]
  }
  sidseq[! snpseq %in% as.character(.snp$code)] <- 2
  sidseq <- as.numeric(sidseq)

  if(is.matrix(snpseq)){
    dim(sidseq) <- dim(snpseq)
  }
  sidseq
} # End of snp2sid().

### Transfer a SNP id sequence to an SNP sequence.
### Eg. 01012... => 1212-..
sid2snp <- function(sidseq){
  snpseq <- sidseq
  for(i in 1:nrow(.snp)){
    snpseq[sidseq == .snp$sid[i]] <- as.character(.snp$code[i])
  }
  snpseq[!sidseq %in% .snp$sid] <- "-"

  if(is.matrix(sidseq)){
    dim(snpseq) <- dim(sidseq)
  }
  snpseq
} # End of sid2snp().


### Transfer a nucleotide seqeunce to a SNP sequence.
### Eg. AGCT-... => 1122-...
code2snp <- function(codeseq){
  snpseq <- codeseq
  snpseq[codeseq %in% c("A", "a", "G", "g")] <- "1"
  snpseq[codeseq %in% c("C", "c", "T", "t")] <- "2"
  snpseq[!(snpseq %in% c("1", "2"))] <- "-"

  if(is.matrix(codeseq)){
    dim(snpseq) <- dim(codeseq)
  }
  snpseq
} # End of code2snp().

### Transfer a SNP seqeunce to a nucleotide sequence.
### Eg. 1122-... => AACC-...
snp2code <- function(snpseq, half = TRUE){
  codeseq <- snpseq

  if(half){
    id <- snpseq == "1"
    tl <- sum(id)
    codeseq[id] <- rep(c("A", "G"), ceiling(tl / 2))[1:tl]
    id <- snpseq == "2"
    tl <- sum(id)
    codeseq[id] <- rep(c("C", "T"), ceiling(tl / 2))[1:tl]
  } else{
    codeseq[snpseq == "1"] <- "A"
    codeseq[snpseq == "2"] <- "C"
  }

  codeseq[!(snpseq %in% c("1", "2"))] <- "-"

  if(is.matrix(snpseq)){
    dim(codeseq) <- dim(snpseq)
  }
  codeseq
} # End of code2snp().


### Transfer nid and sid.
nid2sid <- function(nidseq){
  snp2sid(code2snp(nid2code(nidseq)))
} # End of nid2sid()

sid2nid <- function(sidseq, half = TRUE){
  code2nid(snp2code(sid2snp(sidseq), half = half))
} # End of nid2sid()
#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_control.r"
### This file contains all control functions for EM.

### Control for EM.
.EMControl <- function(
    exhaust.iter = 10,
    fixed.iter = 5,
    short.iter = 100,
    EM.iter = 1000,
    short.eps = 1e-2,
    EM.eps = 1e-10,

    cm.reltol = 1e-16,
    cm.maxit = 5000,

    nm.abstol.Mu.given.QA = 1e-16,
    nm.reltol.Mu.given.QA = 1e-16,
    nm.maxit.Mu.given.QA = 500,
    nm.abstol.QA.given.Mu = 1e-16,
    nm.reltol.QA.given.Mu = 1e-16,
    nm.maxit.QA.given.Mu = 5000,
    est.non.seg.site = FALSE,

    max.init.iter = 100,
    init.procedure = .init.procedure[1],
    init.method = .init.method[1],
    substitution.model = .substitution$model[1],
    edist.model = .edist.model[1],
    identifier = .identifier[1],
    code.type = .code.type[1],
    em.method = .em.method[1],
    boundary.method = .boundary.method[1],

    min.n.class = 1
){
  list(exhaust.iter = as.integer(exhaust.iter),
       fixed.iter = as.integer(fixed.iter),
       short.iter = as.integer(short.iter),
       EM.iter = as.integer(EM.iter),
       short.eps = as.double(short.eps),
       EM.eps = as.double(EM.eps),

       cm.reltol = as.double(cm.reltol),
       cm.maxit = as.integer(cm.maxit),

       nm.abstol.Mu.given.QA = as.double(nm.abstol.Mu.given.QA),
       nm.reltol.Mu.given.QA = as.double(nm.reltol.Mu.given.QA),
       nm.maxit.Mu.given.QA = as.integer(nm.maxit.Mu.given.QA),
       nm.abstol.QA.given.Mu = as.double(nm.abstol.QA.given.Mu),
       nm.reltol.QA.given.Mu = as.double(nm.reltol.QA.given.Mu),
       nm.maxit.QA.given.Mu = as.integer(nm.maxit.QA.given.Mu),
       est.non.seg.site = as.logical(est.non.seg.site),

       max.init.iter = as.integer(max.init.iter),
       init.procedure = init.procedure[1],
       init.method = init.method[1],
       substitution.model = as.character(substitution.model[1]),
       edist.model = edist.model[1],
       identifier = identifier[1],
       code.type = code.type[1],
       em.method = em.method[1],
       boundary.method = boundary.method[1],

       min.n.class = as.integer(min.n.class)
      )
} # End of .EMControl().
#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_file.r"
### Functions to read and write files.

read.fasta <- function(filename, byrow = TRUE, code.type = .code.type[1]){
  if(code.type[1] == "NUCLEOTIDE"){
    ret <- read.fasta.nucleotide(filename, byrow = byrow)
  } else if(code.type[1] == "SNP"){
    stop("SNP is not implemented yet.")
  } else{
    stop("code.type is not found.")
  }

  ret
} # End of read.fasta().

read.phylip <- function(filename, byrow = TRUE, code.type = .code.type[1]){
  if(code.type[1] == "NUCLEOTIDE"){
    ret <- read.phylip.nucleotide(filename, byrow = byrow)
  } else if(code.type[1] == "SNP"){
    ret <- read.phylip.snp(filename, byrow = byrow)
  } else{
    stop("code.type is not found.")
  }

  ret
} # End of read.phylip().

write.fasta <- function(seqdata, filename, classid = NULL,
    seqname = NULL, width.line = 60, lower.case = FALSE,
    code.type = .code.type[1]){
  if(code.type[1] == "NUCLEOTIDE"){
    ret <- write.fasta.nucleotide(seqdata, filename, classid = classid,
             seqname = seqname, width.line = width.line,
             lower.case = lower.case)
  } else if(code.type[1] == "SNP"){
    stop("SNP is not implemented yet.")
  } else{
    stop("code.type is not found.")
  }
} # End of write.fasta().

write.phylip <- function(seqdata, filename, classid = NULL,
    seqname = NULL, width.seqname = 10, width.line = 60, lower.case = FALSE,
    code.type = .code.type[1]){
  if(code.type[1] == "NUCLEOTIDE"){
    ret <- write.phylip.nucleotide(seqdata, filename, classid = classid,
             seqname = seqname, width.seqname = width.seqname,
             width.line = width.line, lower.case = lower.case)
  } else if(code.type[1] == "SNP"){
    ret <- write.phylip.snp(seqdata, filename, classid = classid,
             seqname = seqname, width.seqname = width.seqname,
             width.line = width.line)
  } else{
    stop("code.type is not found.")
  }
} # End of write.phylip().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_file_fasta.r"
### Read FASTA file.
### Return a list contains
###   nseq: Number of sequence.
###   seqlen: Sequence length.
###   seqname: Sequence's name.
###   org: Original sequence in nid. array[nsite, nseq]

read.fasta.nucleotide <- function(filename, byrow = TRUE){
  tmp <- readLines(filename)
  org <- NULL
  tmp.org <- NULL
  seqname <- NULL

  i <- 1
  nseq <- 0
  repeat{
    if(regexpr("^>", tmp[i]) == 1){
      if(!is.null(tmp.org)){
        org <- cbind(org, unlist(strsplit(tmp.org, "")))
        tmp.org <- NULL
      }

      if(regexpr("\\|", tmp[i]) == 1){
        tmp.seqname <- unlist(strsplit(tmp[i], "\\|"))[1]
      } else{
        tmp.seqname <- unlist(strsplit(tmp[i], " "))[1]
      }

      seqname <- c(seqname, gsub(">(.*)", "\\1", tmp.seqname))
      nseq <- nseq + 1
    } else if(tmp[i] %in% c("", " ")){
    } else{
      tmp.org <- paste(tmp.org, tmp[i], sep = "")
    }

    i <- i + 1
    if(i > length(tmp)){
      org <- cbind(org, unlist(strsplit(tmp.org, "")))
      break
    }
  }
  
  fasta <- list(code.type = "NUCLEOTIDE",
              nseq = nseq, seqlen = nrow(org), seqname = seqname,
              org.code = org, org = NULL)
  fasta$org <- matrix(code2nid(org), nrow = fasta$seqlen, ncol = fasta$nseq)

  if(byrow){
    fasta$org.code <- t(fasta$org.code)
    fasta$org <- t(fasta$org)
  }

  class(fasta) <- "seq.data"
  fasta
} # End of read.fasta.nucleotide().


write.fasta.nucleotide <- function(seqdata, filename, classid = NULL,
    seqname = NULL, width.line = 60, lower.case = FALSE){
  seqdata <- apply(seqdata, 2, nid2code, lower.case = lower.case)
  
  n.seq <- nrow(seqdata)
  tl.seq <- ncol(seqdata)

  if(is.null(seqname)){
    seqname <- as.character(1:n.seq)
  }

  if(!is.null(classid)){
    seqname <- cbind(seqname, rep("-", n.seq) , classid)
    seqname <- apply(seqname, 1, paste, collapse = "")
  }

  tl.show <- ceiling(tl.seq / width.line)
  show.range <- list()
  for(i in 1:tl.show){
    show.range[[i]] <- (1:width.line) + (i - 1) * width.line
    if(show.range[[i]][width.line] > tl.seq){
      show.range[[i]] <- show.range[[i]][1:which(show.range[[i]] == tl.seq)]
    }
  }

  ret <- NULL
  for(i in 1:n.seq){
    ret <- c(ret, paste(">", seqname[i], sep = ""))
    for(j in 1:tl.show){
      ret <- c(ret, paste(seqdata[i, show.range[[j]]], collapse = ""))
    }
  }

  write.table(ret, file = filename, quote = FALSE, sep = "", row.names = FALSE,
              col.names = FALSE)
} # End of write.fasta.nucleotide().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_file_phylip.r"
### Read PHYLIP file.
### Return a list contains
###   nseq: Number of sequence.
###   seqlen: Sequence length.
###   seqname: Sequence's name.
###   org: Original sequence in nid. array[nseq, nsite]

read.phylip.nucleotide <- function(filename, byrow = TRUE){
  phylip <- list()
  phylip$code.type <- "NUCLEOTIDE"

### Read header.
  phylip$info <- readLines(filename, n = 1)
  tmp <- unstrsplit(phylip$info, " ")
  tmp <- tmp[tmp != ""]
  phylip$nseq <- as.numeric(tmp[1])
  phylip$seqlen <- as.numeric(tmp[2])

### Read data and transfer to nid.
  op <- options("stringsAsFactors")
  options(stringsAsFactors = FALSE)
  tmp <- read.table(filename, sep = "", quote = "", skip = 1, fill = TRUE)
  options(op)

### Split the data by reading blocks and rejoin them by sequences.
  org <- split(tmp, gl(nrow(tmp) / phylip$nseq, phylip$nseq))
  org <- do.call("cbind", org)
  org <- as.matrix(org, nrow = phylip$nseq)
  phylip$seqname <- org[, 1]
  org <- apply(as.matrix(org[, -1], nrow = phylip$nseq), 1,
               function(y){ unstrsplit(y, "") })
  phylip$org.code <- matrix(org, nrow = phylip$seqlen, ncol = phylip$nseq)
  phylip$org <- matrix(code2nid(org), nrow = phylip$seqlen, ncol = phylip$nseq)

  if(byrow){
    phylip$org.code <- t(phylip$org.code)
    phylip$org <- t(phylip$org)
  }

  class(phylip) <- "seq.data"
  phylip
} # End of read.phylip.nucleotide().


write.phylip.nucleotide <- function(seqdata, filename, classid = NULL,
    seqname = NULL, width.seqname = 10, width.line = 60, lower.case = FALSE){
  seqdata <- apply(seqdata, 2, nid2code, lower.case = lower.case)
  
  n.seq <- nrow(seqdata)
  tl.seq <- ncol(seqdata)

  if(is.null(seqname)){
    seqname <- as.character(1:n.seq)
  }

  if(!is.null(classid)){
    seqname <- cbind(seqname, rep("-", n.seq) , classid)
    seqname <- apply(seqname, 1, paste, collapse = "")
  }

  tmp.seqname <- strsplit(seqname, "")
  for(i in 1:n.seq){
    tl <- length(tmp.seqname[[i]])
    if(tl > width.seqname){
      tmp.seqname[[i]] <- paste(tmp.seqname[[i]][1:width.seqname],
                                collapse = "")
    } else{
      tmp.seqname[[i]] <- paste(c(tmp.seqname[[i]],
                                  rep(" ", width.seqname - tl)), collapse = "")
    }
  }
  tmp.seqname <- do.call("c", tmp.seqname)

  head <- paste(n.seq, tl.seq, collapse = " ")
  write(head, file = filename)

  tl.show <- ceiling(tl.seq / width.line)
  for(i in 1:tl.show){
    show.range <- (1:width.line) + (i - 1) * width.line
    if(show.range[width.line] > tl.seq){
      show.range <- show.range[1:which(show.range == tl.seq)]
    }

    tmp.seqdata <- apply(seqdata[, show.range], 1, paste, collapse = "")

    if(i == 1){
      tmp.seqdata <- cbind(tmp.seqname, tmp.seqdata)
      tmp.seqdata <- apply(tmp.seqdata, 1, paste, collapse = "")
    } 

    write.table(tmp.seqdata, file = filename, append = TRUE, col.names = FALSE,
                row.names = FALSE, sep = "", quote = FALSE)
    write("", file = filename, append = TRUE)
  }
} # End of write.phylip.nucleotide().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_file_seqgen.r"
### Read seq-gen output.
### Return a list contains
###   nseq: Number of sequence.
###   seqlen: Sequence length.
###   seqname: Sequence's name.
###   org: Original sequence in nid. array[nseq, nsite]
read.seqgen <- function(text, byrow = TRUE, code.type = .code.type[1]){
  phylip <- list()
  phylip$code.type <- code.type[1]

### Read header.
  phylip$info <- text[1]
  tmp <- unstrsplit(phylip$info, " ")
  tmp <- tmp[tmp != ""]
  phylip$nseq <- as.numeric(tmp[1])
  phylip$seqlen <- as.numeric(tmp[2])

### Read data and transfer to nid.
  org <- do.call("rbind", lapply(text[2:length(text)], unstrsplit, ""))
  org <- as.matrix(org, nrow = phylip$nseq)
  org.names <- apply(as.matrix(org[, 1:10], nrow = phylip$nseq), 1,
                     paste, collapse = "")
  org.names <- gsub("  *", "", org.names)
  org <- org[, -(1:10)]

### Split the data by reading blocks and rejoin them by sequences.
  phylip$seqname <- org.names
  if(code.type[1] == "NUCLEOTIDE"){
    phylip$org.code <- matrix(org, ncol = phylip$seqlen, nrow = phylip$nseq)
    phylip$org <- matrix(code2nid(org), ncol = phylip$seqlen, nrow = phylip$nseq)
  } else if(code.type[1] == "SNP"){
    org <- code2snp(org)
    phylip$org.code <- matrix(org, ncol = phylip$seqlen, nrow = phylip$nseq)
    phylip$org <- matrix(snp2sid(org), ncol = phylip$seqlen, nrow = phylip$nseq)
  } else{
    stop("code.type is not found.")
  }

  if(!byrow){
    phylip$org.code <- t(phylip$org.code)
    phylip$org <- t(phylip$org)
  }

  class(phylip) <- "seq.data"
  phylip
} # End of read.seqgen().


print.seq.data <- function(x, ...){
  seq.data <- x
  cat("code.type: ", seq.data$code.type,
      ", n.seq: ", seq.data$nseq,
      ", seq.len: ", seq.data$seqlen,
      ".\n", sep = "")
}

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_file_snp.r"
### Read SNP file. (PHYLIP format)
### Return a list contains
###   nseq: Number of sequence.
###   seqlen: Sequence length.
###   seqname: Sequence's name.
###   org: Original sequence in nid. array[nseq, nsite]
read.phylip.snp <- function(filename, byrow = TRUE){
  snp <- list()
  snp$code.type <- "SNP"

### Read header.
  snp$info <- readLines(filename, n = 1)
  tmp <- unstrsplit(snp$info, " ")
  tmp <- tmp[tmp != ""]
  snp$nseq <- as.numeric(tmp[1])
  snp$seqlen <- as.numeric(tmp[2])

### Read data and transfer to nid.
  op <- options("stringsAsFactors")
  options(stringsAsFactors = FALSE)
  tmp <- read.table(filename, sep = "", quote = "", skip = 1, fill = TRUE)
  options(op)

### Split the data by reading blocks and rejoin them by sequences.
  org <- split(tmp, gl(nrow(tmp) / snp$nseq, snp$nseq))
  org <- do.call("cbind", org)
  org <- as.matrix(org, nrow = snp$nseq)
  snp$seqname <- org[, 1]
  org <- apply(as.matrix(org[, -1], nrow = nrow(org)), 1,
               function(y){ unstrsplit(as.character(y), "") })
  snp$org.code <- matrix(org, nrow = snp$seqlen, ncol = snp$nseq)
  snp$org <- matrix(snp2sid(org), nrow = snp$seqlen, ncol = snp$nseq)

  if(byrow){
    snp$org.code <- t(snp$org.code)
    snp$org <- t(snp$org)
  }

  class(snp) <- "seq.data"
  snp
} # End of read.phylip.snp().


write.phylip.snp <- function(seqdata, filename, classid = NULL,
    seqname = NULL, width.seqname = 10, width.line = 60){
  seqdata <- apply(seqdata, 2, sid2snp)

  n.seq <- nrow(seqdata)
  tl.seq <- ncol(seqdata)

  if(is.null(seqname)){
    seqname <- as.character(1:n.seq)
  }

  if(!is.null(classid)){
    seqname <- cbind(seqname, rep("-", n.seq) , classid)
    seqname <- apply(seqname, 1, paste, collapse = "")
  }

  tmp.seqname <- strsplit(seqname, "")
  for(i in 1:n.seq){
    tl <- length(tmp.seqname[[i]])
    if(tl > width.seqname){
      tmp.seqname[[i]] <- paste(tmp.seqname[[i]][1:width.seqname],
                                collapse = "")
    } else{
      tmp.seqname[[i]] <- paste(c(tmp.seqname[[i]],
                                  rep(" ", width.seqname - tl)), collapse = "")
    }
  }
  tmp.seqname <- do.call("c", tmp.seqname)

  head <- paste(n.seq, tl.seq, collapse = " ")
  write(head, file = filename)

  tl.show <- ceiling(tl.seq / width.line)
  for(i in 1:tl.show){
    show.range <- (1:width.line) + (i - 1) * width.line
    if(show.range[width.line] > tl.seq){
      show.range <- show.range[1:which(show.range == tl.seq)]
    }

    tmp.seqdata <- apply(seqdata[, show.range], 1, paste, collapse = "")

    if(i == 1){
      tmp.seqdata <- cbind(tmp.seqname, tmp.seqdata)
      tmp.seqdata <- apply(tmp.seqdata, 1, paste, collapse = "")
    } 

    write.table(tmp.seqdata, file = filename, append = TRUE, col.names = FALSE,
                row.names = FALSE, sep = "", quote = FALSE)
    write("", file = filename, append = TRUE)
  }
} # End of write.phylip.snp().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_find_best.r"
# This file contains functions to run all init method and init procedure.

#.init.procedure <- c("exhaustEM", "emEM", "RndEM", "RndpEM")
#.init.method <- c("randomMu", "NJ", "randomNJ", "PAM", "K-Medoids", "manualMu")
find.best <- function(X, K, EMC = .EMC, manual.id = NULL, byrow = TRUE,
    init.procedure = .init.procedure, init.method = .init.method,
    file.tmp = NULL, visible = FALSE, save.all = FALSE){
  if(K <= 0){
    stop("K > 0.")
  }
  if(byrow){
    X.bycol <- t(X)
  } else{
    X.bycol <- X
  }
  org.logL <- -Inf
  ret <- NULL
  new.ret <- NULL

  if(! exists(".Random.seed")){
    set.seed(12345)
  }

  if(save.all){
    ret.all <- NULL
    count.all <- 1
  }

  for(init.proc in init.procedure){
    for(init.meth in init.method){
      if(((init.proc != "exhaustEM")) &&
          (init.meth %in% c("NJ", "PAM", "manualMu")) ||
         ((K == 1) && (init.proc %in% c("NJ", "randomNJ"))) ||
         ((init.meth == "manualMu") && is.null(manual.id))
        ){
        next
      }

      EMC$init.procedure <- init.proc
      EMC$init.method <- init.meth

      seed.start <- .Random.seed
      if(!is.null(file.tmp)){
        save(list = c("EMC", "seed.start",
                      "new.ret", "ret", "org.logL"), file = file.tmp)
      }

      if(visible){
        cat("Run: ", init.proc, " and ", init.meth, "\n", sep = "")
      }

      new.ret <- try(phyclust(X.bycol, K, EMC, manual.id = manual.id,
                              byrow = FALSE), silent = TRUE)

      if(visible){
        print(new.ret)
        cat("\n")
      }

      seed.end <- .Random.seed
      if(!is.null(file.tmp)){
        save(list = c("EMC", "seed.start", "seed.end",
                      "new.ret", "ret", "org.logL"), file = file.tmp)
      }

      if(save.all){
        if(class(new.ret) == "try-error"){
          new.ret$init.procedure <- init.proc
          new.ret$init.method <- init.meth
        }

        new.ret$seed.start <- seed.start
        new.ret$seed.end <- seed.end

        ret.all[[count.all]] <- new.ret
        count.all <- count.all + 1
      }

      if(class(new.ret) != "try-error"){
        if(is.finite(new.ret$logL) && new.ret$logL > org.logL &&
           new.ret$conv$flag < 2){
          org.logL <- new.ret$logL
          ret <- new.ret
        }
      }
    }
  }

  if(save.all){
    ret$save.all <- save.all
  }
  ret
} # End of find.best().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_mytool.r"
### To unlist a strsplit object.
unstrsplit <- function(x, y, ...){
  unlist(strsplit(x, y, ...))
} # End of unstrsplit().


### Functions for summary() and print.summary().
my.format <- function(x, digits = max(4, getOption("digits") - 3)){
  paste(formatC(x, format = "f", width = -1, digits = digits), collapse = " ")
} # End of my.format().

my.cat <- function(...){
  cat(..., sep = "")
} # End of my.cat().

my.print <- function(x, digits = max(4, getOption("digits") - 3)){
  print.default(x, na.print = "", quote = FALSE, right = TRUE, digits = digits)
} # End of my.print().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_phyclust.r"
# This file contains functions to call phyclust in C.

phyclust <- function(X, K, EMC = .EMC, manual.id = NULL, byrow = TRUE){
  if(K <= 0){
    stop(K > 0)
  }
  if(byrow){
    X <- t(X)
  }

  N.X.org <- ncol(X)
  L <- nrow(X)

  if(! is.null(manual.id)){
    if(max(manual.id) != K){
      stop("manual.id is not correct.")
    } else{
      manual.id <- manual.id - 1
    }
  } else{
    if(EMC$init.method == "manualMu"){
      stop("manual.id is missing.")
    }
  }

  EMC <- check.EMC(EMC)
  EMC <- translate.EMC(EMC)

  ret <- .Call("R_phyclust",
               as.integer(N.X.org),
               as.integer(L),
               as.integer(K),
               as.integer(X),
               EMC,
               as.integer(manual.id),
               PACKAGE = "phyclust")

  if(!is.finite(ret$logL)){
    stop("logL is not finite.\n")
  }

  ret$class.id <- ret$class.id + 1
  ret <- translate.ret(ret, EMC)
  class(ret) <- "phyclust"
  ret
} # End of phyclust().


### For internal used.
check.EMC <- function(EMC){
  if(EMC$init.method %in% c("NJ", "PAM", "manualMu") &&
     EMC$init.procedure != "exhaustEM"){
    my.cat("init procedure: ", EMC$init.procedure, " -> exhaustEM",
           " (", EMC$init.method, ")\n")
    EMC$init.procedure <- "exhaustEM"
    EMC$exhaust.iter <- 1
  }

  if((EMC$code.type == "SNP") &&
     (!EMC$substitution.model %in% c("SNP_JC69", "SNP_F81", "E_SNP_F81"))){
    my.cat("substitution model: ", EMC$substitution.model, " -> SNP_JC69",
           " (", EMC$code.type, ")\n")
    EMC$substitution.model <- "SNP_JC69"
  }

  if(EMC$substitution.model %in% c("SNP_JC69", "SNP_F81", "E_SNP_F81") &&
     EMC$edist.model != "D_HAMMING"){
    my.cat("edist model: ", EMC$edist.model, " -> D_HAMMING",
           " (", EMC$substitution.model, ")\n")
    EMC$edist.model <- "D_HAMMING"
  }

  EMC
} # End of check.EMC().

translate.ret <- function(ret, EMC = NULL){
  ret$Mu <- matrix(ret$Mu, nrow = ret$K, byrow = TRUE)
  if(! is.null(ret$Z.normalized)){
    ret$Z.normalized <- matrix(ret$Z.normalized, nrow = ret$N.X.org, byrow = TRUE)
  }
  ret$QA$pi <- matrix(ret$QA$pi, nrow = ret$K, byrow = TRUE)
  ret$QA$kappa <- matrix(ret$QA$kappa, nrow = ret$K, byrow = TRUE)
  ret$QA$Tt <- matrix(ret$QA$Tt, nrow = ret$K, byrow = TRUE)

  if(!is.null(EMC)){
    ret$init.procedure <- .init.procedure[EMC$init.procedure + 1]
    ret$init.method <- .init.method[EMC$init.method + 1]
    ret$substitution.model <-
      as.character(.substitution$model[EMC$substitution.model + 1])
    ret$edist.model <- .edist.model[EMC$edist.model + 1]
    ret$QA$identifier <- .identifier[EMC$identifier + 1]
    ret$code.type <- .code.type[EMC$code.type + 1]
    ret$em.method <- .em.method[EMC$em.method + 1]
    ret$boundary.method <- .boundary.method[EMC$boundary.method + 1]
  }

  if(ret$QA$identifier == "EE"){
    ret$QA$pi <- matrix(ret$QA$pi[1,], nrow = 1)
    ret$QA$kappa <- ret$QA$kappa[1]
    ret$QA$Tt <- ret$QA$Tt[1]
  } else if(ret$QA$identifier == "EV"){
    ret$QA$pi <- matrix(ret$QA$pi[1,], nrow = 1)
    ret$QA$kappa <- ret$QA$kappa[1]
  } else if(ret$QA$identifier == "VE"){
    ret$QA$Tt <- ret$QA$Tt[1]
  }

  if(ret$code.type == "NUCLEOTIDE"){
    colnames(ret$QA$pi) <- as.character(.nucleotide$code[1:ncol(ret$QA$pi)])
  } else if(ret$code.type == "SNP"){
    colnames(ret$QA$pi) <- as.character(.snp$code[1:ncol(ret$QA$pi)])
  }
  rownames(ret$QA$pi) <- paste("k=", 1:nrow(ret$QA$pi), sep = "")

  if(ret$substitution.model %in% c("JC69", "K80", "SNP_JC69")){
    ret$QA$pi <- NULL
  }
  if(ret$substitution.model %in%
     c("JC69", "F81", "SNP_JC69", "SNP_F81", "E_F81", "E_SNP_F81")){
    ret$QA$kappa <- NULL
  }

  ret
} # End of translate.ret().

translate.EMC <- function(EMC){
  EMC$exhaust.iter <- as.integer(EMC$exhaust.iter)
  EMC$fixed.iter <- as.integer(EMC$fixed.iter)
  EMC$short.iter <- as.integer(EMC$short.iter)
  EMC$EM.iter <- as.integer(EMC$EM.iter)
  EMC$short.eps <- as.double(EMC$short.eps)
  EMC$EM.eps <- as.double(EMC$EM.eps)

  EMC$cm.reltol <- as.double(EMC$cm.reltol)
  EMC$cm.maxit <- as.integer(EMC$cm.maxit)

  EMC$nm.abstol.Mu.given.QA <- as.double(EMC$nm.abstol.Mu.given.QA)
  EMC$nm.abstol.QA.given.Mu <- as.double(EMC$nm.abstol.QA.given.Mu)
  EMC$nm.reltol.Mu.given.QA <- as.double(EMC$nm.reltol.Mu.given.QA)
  EMC$nm.reltol.QA.given.Mu <- as.double(EMC$nm.reltol.QA.given.Mu)
  EMC$nm.maxit.Mu.given.QA <- as.integer(EMC$nm.maxit.Mu.given.QA)
  EMC$nm.maxit.QA.given.Mu <- as.integer(EMC$nm.maxit.QA.given.Mu)
  EMC$est.non.seg.site <- as.integer(EMC$est.non.seg.site)

  EMC$min.n.class <- as.integer(EMC$min.n.class)

  if(EMC$init.procedure[1] %in% .init.procedure){
    EMC$init.procedure <- as.integer(which(EMC$init.procedure ==
                                           .init.procedure) - 1)
  } else{
    stop("The initial procedure is not found.")
  }

  if(EMC$init.method[1] %in% .init.method){
    EMC$init.method <- as.integer(which(EMC$init.method == .init.method) - 1)
  } else{
    stop("The initial method is not found.")
  }

  if(EMC$substitution.model[1] %in% as.character(.substitution$model)){
    EMC$substitution.model <-
      as.integer(which(EMC$substitution.model ==
                       as.character(.substitution$model)) - 1)
  } else{
    stop("The substitution model is not found.")
  }

  if(EMC$edist.model[1] %in% .edist.model){
    EMC$edist.model <- as.integer(which(EMC$edist.model == .edist.model) - 1)
  } else{
    stop("The edist model is not found.")
  }

  if(EMC$identifier[1] %in% .identifier){
    EMC$identifier <- as.integer(which(EMC$identifier == .identifier) - 1)
  } else{
    stop("The identifier is not found.")
  }

  if(EMC$code.type[1] %in% .code.type){
    EMC$code.type <- as.integer(which(EMC$code.type == .code.type) - 1)
  } else{
    stop("The code type is not found.")
  }

  if(EMC$em.method[1] %in% .em.method){
    EMC$em.method <- as.integer(which(EMC$em.method == .em.method) - 1)
  } else{
    stop("The em method is not found.")
  }

  if(EMC$boundary.method[1] %in% .boundary.method){
    EMC$boundary.method <- as.integer(which(EMC$boundary.method ==
                                            .boundary.method) - 1)
  } else{
    stop("The boundary method is not found.")
  }

  EMC$max.init.iter <- as.integer(EMC$max.init.iter)

  EMC
} # End of translate.EMC().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_phyclust_Pt.r"
### This file contains functions for computing P(t) = e^{Qt}.
# NUCLEOTIDE: Q <- list(pi = c(0.25, 0.25, 0.25, 0.25), kappa = 0.5)
# SNP: Q <- list(pi = c(0.5, 0.5), kappa = 0.5)

phyclust.Pt <- function(Q, Tt, substitution.model = .substitution$model[1],
    log = FALSE){
  sub.model <- which(substitution.model[1] == as.character(.substitution$model))
  code.type <- which(as.character(.substitution$code.type[sub.model]) ==
                                  .code.type)
  ret <- .Call("R_phyclust_logPt",
               as.double(Q$pi),
               as.double(Q$kappa),
               as.double(Tt),
               as.integer(code.type - 1),
               as.integer(sub.model - 1),
               PACKAGE = "phyclust")
  if(!log){
    ret <- exp(ret)
  }
  if(as.character(.substitution$code.type[sub.model]) == "NUCLEOTIDE"){
    ret <- matrix(ret, nrow = 4, ncol = 4, byrow = TRUE)
  } else{
    ret <- matrix(ret, nrow = 2, ncol = 2, byrow = TRUE)
  }
  attr(ret, "class") <- "Pt"
  attr(ret, "code.type") <- as.character(.substitution$code.type[sub.model])
  attr(ret, "log") <- log
  ret
}

### Print Pt matrix.
print.Pt <- function(x, ...){
  Pt <- x
  if(attr(Pt, "code.type") == "NUCLEOTIDE"){
    colnames(Pt) <- .nucleotide$code[1:4]
    rownames(Pt) <- .nucleotide$code[1:4]
  } else{
    colnames(Pt) <- .snp$code[1:2]
    rownames(Pt) <- .snp$code[1:2]
  }
  if(attr(Pt, "log") == TRUE){
    cat("logPt\n")
  } else{
    cat("Pt\n")
  }
  attr(Pt, "class") <- NULL
  attr(Pt, "code.type") <- NULL
  attr(Pt, "log") <- NULL
  my.print(Pt)
}

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_phyclust_edist.r"
### This file contains functions for compute evolution distances.

phyclust.edist <- function(X, edist.model = .edist.model[1]){
  if(is.vector(X)) stop("X should be a matrix and nrow > 1.")

  if(edist.model[1] %in% .edist.model){
    edist.model <- which(edist.model[1] == .edist.model) - 1
  } else{
    stop("The distance model is not found.")
  }
  N <- nrow(X)
  L <- ncol(X)

  d <- .Call("R_phyclust_edist",
             as.integer(edist.model),
             as.integer(N),
             as.integer(L),
             as.integer(t(X)),
             PACKAGE = "phyclust")
  class(d) <- "dist"
  attr(d, "Size") <- N
  attr(d, "Diag") <- FALSE
  attr(d, "Upper") <- FALSE
  attr(d, "method") <- .edist.model[edist.model[1] + 1]
  d
} # End of edist().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_phyclust_em_step.r"
### This file contains functions to call phyclust_em_step, phyclust_e_step, and
### phyclust_m_step in C.

### EM step.
phyclust.em.step <- function(X, ret.phyclust = NULL,
    K = NULL, Eta = NULL, Mu = NULL, pi = NULL, kappa = NULL, Tt = NULL,
    substitution.model = NULL, identifier = NULL, code.type = NULL){
  if(is.null(ret.phyclust)){
    if(is.null(K) || is.null(Eta) || is.null(Mu) || is.null(Tt) ||
       is.null(substitution.model) ||
       is.null(identifier) || is.null(code.type)){
      stop("Parameters are not specified correctly.")
    } else{
      ret.phyclust <- list(K = K, Eta = Eta, Mu = Mu,
                           QA = list(pi = pi, kappa = kappa, Tt = Tt),
                           substitution.model = substitution.model,
                           identifier = identifier,
                           code.type = code.type)
    }
  } else{
    if(class(ret.phyclust) != "phyclust"){
      stop("ret.phyclust should be in a phyclust class.")
    }
  }

  vect <- convert.QA.to.vect(ret.phyclust)

  ret <- .Call("R_phyclust_em_step",
               as.integer(nrow(X)),
               as.integer(ncol(X)),
               as.integer(t(X)),
               as.integer(ret.phyclust$K),
               as.double(ret.phyclust$Eta),
               as.integer(t(ret.phyclust$Mu)),
               as.double(vect),
               as.integer(which(ret.phyclust$substitution.model ==
                                as.character(.substitution$model)) - 1),
               as.integer(which(ret.phyclust$QA$identifier == .identifier) - 1),
               as.integer(which(ret.phyclust$code.type == .code.type) - 1),
               PACKAGE = "phyclust")

  ret$Z.normalized <- ret$bic <- ret$aic <- ret$icl <-
    ret$class.id <- ret$n.class <- NULL
  ret$substitution.model <- ret.phyclust$substitution.model
  ret$QA$identifier <- ret.phyclust$QA$identifier
  ret$code.type <- ret.phyclust$code.type

  ret <- translate.ret(ret)
  class(ret) <- "phyclust"
  ret
} # End of phyclust.em.step().


### E-step: return a matrix Z.normalized with dim=NxK.
### Z.state = 1 return Z.normalized
###           0 return logPt
phyclust.e.step <- function(X, ret.phyclust = NULL,
    K = NULL, Eta = NULL, Mu = NULL, pi = NULL, kappa = NULL, Tt = NULL,
    substitution.model = NULL, identifier = NULL, code.type = NULL,
    Z.state = TRUE){
  if(is.null(ret.phyclust)){
    if(is.null(K) || is.null(Eta) || is.null(Mu) || is.null(Tt) ||
       is.null(substitution.model) ||
       is.null(identifier) || is.null(code.type)){
      stop("Parameters are not specified correctly.")
    } else{
      ret.phyclust <- list(K = K, Eta = Eta, Mu = Mu,
                           QA = list(pi = pi, kappa = kappa, Tt = Tt,
                                     identifier = identifier),
                           substitution.model = substitution.model,
                           code.type = code.type)
    }
  } else{
    if(class(ret.phyclust) != "phyclust"){
      stop("ret.phyclust should be in a phyclust class.")
    }
  }

  vect <- convert.QA.to.vect(ret.phyclust)

  ret <- .Call("R_phyclust_e_step",
               as.integer(nrow(X)),
               as.integer(ncol(X)),
               as.integer(t(X)),
               as.integer(ret.phyclust$K),
               as.double(ret.phyclust$Eta),
               as.integer(t(ret.phyclust$Mu)),
               as.double(vect),
               as.integer(which(ret.phyclust$substitution.model ==
                                as.character(.substitution$model)) - 1),
               as.integer(which(ret.phyclust$QA$identifier == .identifier) - 1),
               as.integer(which(ret.phyclust$code.type == .code.type) - 1),
               as.integer(Z.state),
               PACKAGE = "phyclust")
  ret <- matrix(ret, nrow = nrow(X), byrow = TRUE)

  ret
} # End of phyclust.e.step().


### M-step: return a object with phyclust class.
phyclust.m.step <- function(X, ret.phyclust = NULL,
    K = NULL, Z.normalized = NULL,
    substitution.model = NULL, identifier = NULL, code.type = NULL){
  if(is.null(ret.phyclust)){
    if(is.null(K) || is.null(Z.normalized) ||
       is.null(substitution.model) ||
       is.null(identifier) || is.null(code.type)){
      stop("Parameters are not specified correctly.")
    } else{
      ret.phyclust <- list(K = K,
                           QA = list(identifier = identifier),
                           Z.normalized = Z.normalized,
                           substitution.model = substitution.model,
                           code.type = code.type)
    }
  } else{
    if(class(ret.phyclust) != "phyclust"){
      stop("ret.phyclust should be in a phyclust class.")
    }
  }

  ret <- .Call("R_phyclust_m_step",
               as.integer(nrow(X)),
               as.integer(ncol(X)),
               as.integer(t(X)),
               as.integer(ret.phyclust$K),
               as.double(t(ret.phyclust$Z.normalized)),
               as.integer(which(ret.phyclust$substitution.model ==
                                as.character(.substitution$model)) - 1),
               as.integer(which(ret.phyclust$QA$identifier == .identifier) - 1),
               as.integer(which(ret.phyclust$code.type == .code.type) - 1),
               PACKAGE = "phyclust")

  ret$Z.normalized <- ret$bic <- ret$aic <- ret$icl <-
    ret$class.id <- ret$n.class <- NULL
  ret$substitution.model <- ret.phyclust$substitution.model
  ret$QA$identifier <- ret.phyclust$QA$identifier
  ret$code.type <- ret.phyclust$code.type

  ret <- translate.ret(ret)
  class(ret) <- "phyclust"
  ret
} # End of phyclust.m.step().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_phyclust_logL.r"
# This file contains functions to call phyclust_logL in C.

phyclust.logL <- function(X, ret.phyclust = NULL,
    K = NULL, Eta = NULL, Mu = NULL, pi = NULL, kappa = NULL, Tt = NULL,
    substitution.model = NULL, identifier = NULL, code.type = NULL){
  if(is.null(ret.phyclust)){
    if(is.null(K) || is.null(Eta) || is.null(Mu) || is.null(Tt) ||
       is.null(substitution.model) ||
       is.null(identifier) || is.null(code.type)){
      stop("Parameters are not specified correctly.")
    } else{
      ret.phyclust <- list(K = K, Eta = Eta, Mu = Mu,
                           QA = list(pi = pi, kappa = kappa, Tt = Tt),
                           substitution.model = substitution.model,
                           identifier = identifier,
                           code.type = code.type)
    }
  } else{
    if(class(ret.phyclust) != "phyclust"){
      stop("ret.phyclust should be in a phyclust class.")
    }
  }

  if(ret.phyclust$substitution.model == "E_F81"){
    ret.phyclust$substitution.model <- "F81"
  } else if(ret.phyclust$substitution.model == "E_HKY85"){
    ret.phyclust$substitution.model <- "HKY85"
  } else if(ret.phyclust$substitution.model == "E_SNP_F81"){
    ret.phyclust$substitution.model <- "SNP_F81"
  }

  vect <- convert.QA.to.vect(ret.phyclust)

  ret <- .Call("R_phyclust_logL",
               as.integer(nrow(X)),
               as.integer(ncol(X)),
               as.integer(t(X)),
               as.integer(ret.phyclust$K),
               as.double(ret.phyclust$Eta),
               as.integer(t(ret.phyclust$Mu)),
               as.double(vect),
               as.integer(which(ret.phyclust$substitution.model ==
                                as.character(.substitution$model)) - 1),
               as.integer(which(ret.phyclust$QA$identifier == .identifier) - 1),
               as.integer(which(ret.phyclust$code.type == .code.type) - 1),
               PACKAGE = "phyclust")

  ret
} # End of phyclust.logL().


### For internal used.
# Storage of vect for EE:
# 	0 to (n_param - 1):				for k = 0, 1, 2,..., K-1.
# Storage of vect for EV:
# 	0 to (n_param - 2):				for k = 0, 1, 2,..., K-1.
#	(n_param - 1) to (n_param - 1 + K):		for Tt[0] to Tt[K-1].
# Storage of vect for VE:
# 	0 to (n_param - 2):				for k = 0,
#	(n_param - 1) to (2*(n_param - 1) - 1):		for k = 1,
#	(2*(n_param - 1)) to (3*(n_param - 1) - 1):	for k = 2,
#	... until k = K-1.
#	((K-1)*(n_param - 1))				for Tt.
# Storage of vect for VV:
#	0 to (n_param - 1):				for k = 0,
#	(n_param) to (2*n_param - 1):			for k = 1,
#	(2*n_param) to (3*n_param - 1):			for k = 2,
# 	... until k = K-1.
convert.QA.to.vect <- function(ret.phyclust){
  if(ret.phyclust$code.type == "NUCLEOTIDE"){
    ncode <- 3
  } else if(ret.phyclust$code.type == "SNP"){
    ncode <- 1
  }
  model.with.pi <- c("F81", "HKY85", "SNP_F81")
  model.with.kappa <- c("K80", "F81", "HKY85", "SNP_F81")

  K <- ret.phyclust$K
  vect <- NULL

  if(ret.phyclust$QA$identifier == "EE"){
    if(ret.phyclust$substitution.model %in% model.with.pi){
      vect <- c(vect, ret.phyclust$QA$pi[1:ncode])
    }
    if(ret.phyclust$substitution.model %in% model.with.kappa){
      vect <- c(vect, ret.phyclust$QA$kappa)
    }
    vect <- c(vect, ret.phyclust$QA$Tt)
  } else if(ret.phyclust$QA$identifier == "EV"){
    if(ret.phyclust$substitution.model %in% model.with.pi){
      vect <- c(vect, ret.phyclust$QA$pi[1:ncode])
    }
    if(ret.phyclust$substitution.model %in% model.with.kappa){
      vect <- c(vect, ret.phyclust$QA$kappa)
    }
    vect <- c(vect, ret.phyclust$QA$Tt)
  } else if(ret.phyclust$QA$identifier == "VE"){
    for(i in 1:K){
      if(ret.phyclust$substitution.model %in% model.with.pi){
        vect <- c(vect, ret.phyclust$QA$pi[i, 1:ncode])
      }
      if(ret.phyclust$substitution.model %in% model.with.kappa){
        vect <- c(vect, ret.phyclust$QA$kappa[i])
      }
    }
    vect <- c(vect, ret.phyclust$QA$Tt)
  } else{
    for(i in 1:K){
      if(ret.phyclust$substitution.model %in% model.with.pi){
        vect <- c(vect, ret.phyclust$QA$pi[i, 1:ncode])
      }
      if(ret.phyclust$substitution.model %in% model.with.kappa){
        vect <- c(vect, ret.phyclust$QA$kappa[i])
      }
      vect <- c(vect, ret.phyclust$QA$Tt[i])
    }
  }

  vect
} # End of convert.QA.to.vect().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_plotdots.r"
### This file contains the function for the dots plot.
### Plot the difference for all sequences based on the first sequence.

plotdots <- function(X, X.class = NULL,
    diff.only = TRUE, fill = FALSE, label = TRUE, xlim = NULL,
    main = "Dots Plot", xlab = "Sites", ylab = "Sequences", ...){
  if(! is.null(X.class)){
    id.X <- order(X.class)
    X <- X[id.X,]
    K <- max(X.class)
    X.nc.cum <- rep(0, K)
    for(k in 1:K) X.nc.cum[k] <- sum(X.class %in% 1:k)
    header <- which(id.X == 1)
  } else{
    header <- 1
  }
  if(diff.only){
    tl.diff <- apply(X, 2, function(x) length(unique(x)))
    X <- X[, tl.diff > 1]
  }

  N <- nrow(X)
  L <- ncol(X)

  if(is.null(xlim)) xlim <- c(1, L + 1)
  if(label){
    ylim <- c(N + 3, 1)
  } else{
    ylim <- c(N + 1, 1)
  }

  plot(NULL, NULL, type = "n", xlim = xlim, ylim = ylim,
       main = main, xlab = xlab, ylab = ylab)

  my.col <- c("green3", "blue2", "#CC00CC", "red2", "white")
  flag.diff <- rowSums(t(X[-header,]) != X[header,]) > 0

  for(i in 1:N){
    for(j in 1:L){
      if(diff.only){
         if(! flag.diff[j]) next
      }
      x.left <- j
      y.top <- i
      if(X[i, j] != X[header, j] || i == header){
        rect(x.left, y.top + 1, x.left + 1, y.top,
             col = my.col[X[i, j] + 1], border = NA)
      } else{
        if(fill & flag.diff[j]){
          rect(x.left, y.top + 1, x.left + 1, y.top,
               col = my.col[X[i, j] + 1], border = NA)
        }
      }
    }
  }

  if(label){
    i <- N + 1
    for(j in 1:L){
      x.left <- j
      y.top <- i
      if(flag.diff[j]){
        rect(x.left, y.top + 2, x.left + 1, y.top, col = "orange", border = NA)
      }
    }
    abline(h = y.top, lty = 3, lwd = 0.5)
  }

  if(! is.null(X.class)) abline(h = X.nc.cum + 1, lty = 3, lwd = 0.5)
} # End of plotdots().
#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_plotfd.r"
### This file contains the function for fitness distirbution plot. 
### Plot the posterior distribution by each mixture component with
### sorted cluster.

plotfd <- function(X, pcobj, merge = TRUE,
    main = NULL, xlab = NULL, ylab = NULL, ...){
  N <- pcobj$N
  K <- pcobj$K

  post.class <- sort(pcobj$class.id)
  post <- pcobj$Z.normalized[order(pcobj$class.id),]
  for(k in 1:K){
    id <- which(post.class == k)
    id.new <- order(post[id, k], decreasing = TRUE)
    post[id,] <- post[id[id.new],]
  }

  xlim <- c(1, N)
  ylim <- c(0, 1)

  if(is.null(main)) main <- "Distribution Plot"
  if(is.null(xlab)) xlab <- "Sequence (sorted)"
  if(is.null(ylab)) ylab <- "Posterior"

  color <- .Color[1:K %% length(.Color) + 1]

  x <- 1:N
  if(merge){
    if(K == 1){
      plot(NULL, NULL, xlim = xlim, ylim = ylim,
           main = main, xlab = xlab, ylab = "K = 1")
      y <- post[, 1]
      lines(list(x = x, y = y), col = color[1])
    } else{
      plot(NULL, NULL, xlim = xlim, ylim = ylim,
           main = main, xlab = xlab, ylab = ylab)
      for(k in 1:K){
        lines(list(x = x, y = post[, k]), col = color[k], lwd = 3, lty = k)
      }
    }
  } else{
    v <- cumsum(pcobj$n.class)[1:(K-1)] + 0.5
    if(K == 1){
      plot(NULL, NULL, xlim = xlim, ylim = ylim,
           main = main, xlab = xlab, ylab = "K = 1")
      y <- post[, 1]
      lines(list(x = x, y = y), col = color[1])
      abline(v = v, lty = 3)
    } else{
      par(mfrow = c(K, 1))
      plot(NULL, NULL, xlim = xlim, ylim = ylim,
           main = main, xlab = "", ylab = "K = 1", axes = FALSE)
      lines(list(x = x, y = post[, 1]), col = color[1], lwd = 2)
      axis(2)
      abline(v = v, lty = 3)
      for(k in 2:(K - 1)){
        plot(NULL, NULL, xlim = xlim, ylim = ylim,
            xlab = "", ylab = paste("K =", k), axes = FALSE)
        lines(list(x = x, y = post[, k]), col = color[k], lwd = 2)
        axis(2)
        abline(v = v, lty = 3)
      }
      plot(NULL, NULL, xlim = xlim, ylim = ylim,
           main = "", xlab = xlab, ylab = paste("K =", K), axes = FALSE)
      lines(list(x = x, y = post[, K]), col = color[K], lwd = 2)
      axis(1)
      axis(2)
      abline(v = v, lty = 3)
    }
  }
} # End of plotfd().
#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_plothist.r"
### This file contains the functions for histograms.

plothist <- function(X, X.class = NULL, Mu = NULL, fill.color = .Color,
    draw.all = TRUE, main = "Mutation counts",
    xlab = "Difference", ylab = "Counts", append = FALSE){
  X.tmp <- t(X)

  if(is.null(Mu)){
    if(is.null(X.class)){
      common <- X.tmp[, 1]
    } else{
      common <- X.tmp[, X.class == 1][, 1]
    }
  } else{
    common <- Mu[1,]
  }

  if(is.null(X.class)){
    color <- .Color[1]

    X.table <- table(colSums(X.tmp != common))
    X.ht <- cbind(as.numeric(names(X.table)), as.numeric(X.table))
    plothist.my(X.ht, main = main, xlab = xlab, ylab = ylab, fill.color = color)
  } else{
    K <- length(unique(X.class))
    color <- .Color[1:K %% length(.Color) + 1]

    X.ht <- NULL
    xlim <- NULL
    ylim <- NULL
    for(k in 1:K){
      X.k.tmp <- X.tmp[, X.class == k]
      X.table <- table(colSums(X.k.tmp != common))
      X.ht[[k]] <- cbind(as.numeric(names(X.table)), as.numeric(X.table))
      xlim <- c(xlim, X.ht[[k]][, 1])
      ylim <- c(ylim, X.ht[[k]][, 2])
    }
    xlim <- range(xlim) + c(-0.6, 0.6)
    ylim <- c(0, max(ylim) + 0.1)

    if(draw.all){
      if(! append){
        par(mfrow = c(K + 1, 1), mar = c(2, 4, 0.5, 2))
      }
      X.table <- table(colSums(X.tmp != common))
      X.ht.all <- cbind(as.numeric(names(X.table)), as.numeric(X.table))
      ylim <- c(0, max(X.ht.all[, 2]) + 0.1)
      plothist.my(X.ht.all, xlim = xlim, ylim = ylim,
                 main = "", xlab = "", ylab = "")
    } else{
      if(! append){
        par(mfrow = c(K, 1), mar = c(2, 4, 0.5, 2))
      }
    }

    for(k in 1:K){
      plothist.my(X.ht[[k]], xlim = xlim, ylim = ylim, fill.color = color[k],
                 main = "", xlab = "", ylab = "")
    }
  }
} # End of plothist().


plothist.my <- function(X.ht, xlim = NULL, ylim = NULL, fill.color = NULL,
  main = NULL, xlab = NULL, ylab = NULL){
  if(is.null(xlim)) xlim <- range(X.ht[, 1]) + c(-0.6, 0.6)
  if(is.null(ylim)) ylim <- c(0, max(X.ht[, 2]) + 0.1)

  plot(NULL, NULL, type = "n", xlim = xlim, ylim = ylim, axes = FALSE,
       main = main, ylab = ylab, xlab = xlab)

  for(i in 1:nrow(X.ht)){
    x.left <- X.ht[i, 1] - 0.5
    x.right <- X.ht[i, 1] + 0.5
    y.top <- X.ht[i, 2]
    y.bottom <- 0
    rect(x.left, y.bottom, x.right, y.top, col = fill.color, border = 1)
  }

  axis(1)
  axis(2)
} # End of plothist.my().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_plotnj.r"
### This file contains functions for plot the unrooted tree.

plotnj <- function(unrooted.tree, X.class = NULL, type = "u", main = NULL,
    show.tip.label = FALSE, show.node.label = FALSE,
    edge.width = 1, edge.width.class = edge.width, ...){
  if(class(unrooted.tree) != "phylo"){
    stop("This is not a phylo tree.")
  }
  if(any(unrooted.tree$edge < 0)){
    stop("This is not a unrooted tree.")
  }
  if(is.null(X.class)){
    plot(unrooted.tree, type = type, show.tip.label = show.tip.label,
         main = main, ...)
  } else{
    unrooted.tree.edge.color <- rep("#000000", length(unrooted.tree$edge.length))
    unrooted.tree.edge.width <- rep(edge.width, length(unrooted.tree$edge.length))
    for(k in sort(unique(X.class))){
      unrooted.tree.edge.color[unrooted.tree$edge[, 2] %in% which(X.class == k)] <-
        .Color[(k - 1) %% length(.Color) + 1]
      unrooted.tree.edge.width[unrooted.tree$edge[, 2] %in% which(X.class == k)] <-
        edge.width.class
    }
    plot(unrooted.tree, type = type,
         edge.color = unrooted.tree.edge.color,
         edge.width = unrooted.tree.edge.width,
         show.tip.label = show.tip.label, main = main, ...)
  }
  if(show.node.label) nodelabels(frame = "n")
} # End of plotnj().
#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_show_option.r"
### This file contains function to show all possible options in phyclust.

.show.option <- function(){
  my.cat("Boundary method: ", paste(.boundary.method, collapse = ", "), "\n")
  my.cat("Code type: ", paste(.code.type, collapse = ", "), "\n")
  my.cat("Edist model: ", paste(.edist.model, collapse = ", "), "\n")
  my.cat("EM method: ", paste(.em.method, collapse = ", "), "\n")
  my.cat("Identifier: ", paste(.identifier, collapse = ", "), "\n")
  my.cat("Init method: ", paste(.init.method, collapse = ", "), "\n")
  my.cat("Init procedure: ", paste(.init.procedure, collapse = ", "), "\n")
  my.cat("Standard code: \n")
  my.print(as.matrix(.nucleotide))
  my.print(as.matrix(.snp))
  my.cat("Substitution model: \n")
  my.print(as.matrix(.substitution))
  invisible()
} # End of .show.option()
#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/f_summary.r"
### This file contains functions for print and summary.

### Print phyclust
print.phyclust <- function(x, digits = max(4, getOption("digits") - 3), ...){
  pcobj <- x

  options(digits = digits)
  init <- NULL

  if(!is.null(pcobj$conv)){
    my.cat("Phyclust Results:\n",
           "code type: ", pcobj$code.type,
           ", em method: ", pcobj$em.method,
           ", boundary method: ", pcobj$boundary.metho, ".\n",
           "init procedure: ", pcobj$init.procedure,
           ", method: ", pcobj$init.method, ".\n",
           "model substitution: ", pcobj$substitution.model,
           ", distance: ", pcobj$edist.model, ".\n",
           "iter: ", pcobj$conv$iter,
           " ", pcobj$conv$inner.iter,
           " ", pcobj$conv$cm.iter,
           ", convergence: ", pcobj$conv$flag,
           ", check.param: ", pcobj$conv$check.param, ".\n",
           "eps: ", pcobj$conv$eps,
           ", error: ", pcobj$conv$error, ".\n")
  }
      
  my.cat("N.X.org: ", pcobj$N.X.org,
         ", N.X.unique: ", pcobj$N.X.unique,
         ", L: ", pcobj$L,
         ", K: ", pcobj$K,
         ", p: ", pcobj$p,
         ", N.seg.site: ", pcobj$N.seg.site, ".\n",
         "logL: ", pcobj$logL)
  if(!is.null(pcobj$bic)) my.cat(", bic: ", pcobj$bic)
  if(!is.null(pcobj$aic)) my.cat(", aic: ", pcobj$aic)
  if(!is.null(pcobj$icl)) my.cat(", icl: ", pcobj$icl)
  my.cat("\n")
  my.cat("identifier: ", pcobj$QA$identifier, "\n")
  cat("  Eta:", pcobj$Eta, "\n")
  if(!is.null(pcobj$Q$pi)){
    my.cat("  pi:\n")
    my.print(pcobj$Q$pi)
  }
  if(!is.null(pcobj$Q$kappa)) cat("  kappa:", pcobj$Q$kappa, "\n")
  cat("  Tt:", pcobj$Q$Tt, "\n")
  if(!is.null(pcobj$n.class)) cat("  n.class:", pcobj$n.class, "\n")
} # End of print.phyclust().

summary.phyclust <- function(object, ...){
  pcobj <- object

  print.phyclust(pcobj)
  cat("Mu:\n")
  for(i in 1:pcobj$K){
    cat("   ", as.character(.nucleotide$code[pcobj$Mu[i,] + 1]), "\n")
  }
  cat("Class id:", pcobj$class.id, "\n")
} # End of summary.phyclust().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/fcn_RRand.r"
### This file contains a wrap to call C function in "src/R_RRand.c".
### Written: Wei-Chen Chen on 2008/10/27.


# Call:
#   SEXP R_RRand(SEXP N, SEXP TRUK, SEXP PREDK, SEXP trcl, SEXP prcl)
# Input:
#   N: SEXP[1], number of observations.
#   TRUK: SEXP[1], number of true clusters.
#   PREDK: SEXP[1], number of predicted clusters.
#   trcl: SEXP[N], true cluster ids.
#   prcl: SEXP[N], predicted cluster ids.
# Output in C:
#   ret: a list contains
#      Rand: SEXP[1], Rand index.
#      adjRand: SEXP[1], adjust Rand index.
#      Eindex: SEXP[1], Eindex.
RRand <- function(trcl, prcl, lab = NULL){
  if(! is.null(lab)){
    trcl <- trcl[lab == 0]
    prcl <- prcl[lab == 0]
  }

  N <- length(trcl)

  if(length(trcl) != N || length(prcl) != N){
    stop("The lengths of trcl and prcl do not agree!")
  } 

  tmp.TRUK <- unique(trcl)
  tmp.PREDK <- unique(prcl)
  TRUK <- max(tmp.TRUK)
  PREDK <- max(tmp.PREDK)

  if(min(tmp.TRUK) < 1 || min(tmp.PREDK) < 1){
    stop("The minimum ID is smaller than 1!")
  }
  if(min(tmp.TRUK) != 1 || min(tmp.PREDK) != 1){
    warnings("The minimum ID is not 1!")
  }

  ret <- .Call("R_RRand",
               as.integer(N),
               as.integer(TRUK),
               as.integer(PREDK),
               as.integer(trcl - 1),
               as.integer(prcl - 1),
               PACKAGE = "phyclust")

  class(ret) <- "RRand"
  ret
} # End of RRand().

print.RRand <- function(x, digits = max(4, getOption("digits") - 3), ...){
  RRand <- x
  my.print(unlist(RRand), digits = digits)
} # End of print.RRand().
#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/p_ms.r"
### This file contains functions for "ms".

ms <- function(nsam = NULL, nreps = 1, opts = NULL){
  temp.file <- tempfile("ms.")

  if(! is.null(opts) && ! is.null(nsam)){
    if(nsam >= 2){
      nsam <- as.character(nsam)
      nreps <- as.character(nreps)
      argv <- c("ms", nsam, nreps, unlist(strsplit(opts, " ")))

      .Call("R_ms_main", argv, temp.file, PACKAGE = "phyclust")
      ret <- scan(file = temp.file, what = "character", sep = "\n", quiet = TRUE)
      unlink(temp.file)

      class(ret) <- "ms"
      return(ret)
    }
  }

  argv <- c("ms", "-h")
  try(.Call("R_ms_main", argv, temp.file, PACKAGE = "phyclust"), silent = TRUE)
  unlink(temp.file)
  invisible()
} # End of ms().

print.ms <- function(x, ...){
  ms <- x
  cat(ms, sep = "\n")
} # End of print.ms().
#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/p_seqgen.r"
### This file contains functions for "seq-gen".

seqgen <- function(opts = NULL, rooted.tree = NULL, newick.tree = NULL,
    input = NULL){
  argv <- "seq-gen"
  temp.file.seqgen <- tempfile("seqgen.")
  temp.file.ms <- tempfile("ms.")

  if(! is.null(opts)){
    if((! is.null(rooted.tree)) && (class(rooted.tree) == "phylo")){
      newick.tree <- write.tree(rooted.tree, digits = 12)
    }
    if(! is.null(newick.tree)){
      write(newick.tree, file = temp.file.ms, sep = "")
    } else if(! is.null(input)){
      write(input, file = temp.file.ms, sep = "\n")
    } else{
      stop("A newick or rooted/phylo tree is required.")
    }

    argv <- c(argv, unlist(strsplit(opts, " ")), temp.file.ms)
    .Call("R_seq_gen_main", argv, temp.file.seqgen, PACKAGE = "phyclust")

    ret <- scan(file = temp.file.seqgen,
                what = "character", sep = "\n", quiet = TRUE)
    class(ret) <- "seqgen"

    unlink(temp.file.ms)
    unlink(temp.file.seqgen)
    return(ret)
  }

  argv <- c(argv, "-h")
  try(.Call("R_seq_gen_main", argv, temp.file.ms, PACKAGE = "phyclust"),
      silent = TRUE)
  unlink(temp.file.ms)
  unlink(temp.file.seqgen)
  invisible()
} # End of seqgen().

print.seqgen <- function(x, ...){
  seqgen <- x
  cat(seqgen, sep = "\n")
} # Enf of print.seqgen().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/p_tool_seq.r"
### Generate sequences by the HKY85 model.
gen.seq.HKY <- function(rooted.tree, pi, kappa, L, rate.scale = 1,
    anc.seq = NULL){
  if(!is.rooted(rooted.tree)){
    stop("A rooted tree is required.")
  }
  if(length(pi) != 4 || sum(pi) != 1 || any(pi <= 0) || any(pi >= 1)){
    stop("pi is not correct.")
  }
  if(length(kappa) != 1 || kappa <= 0){
    stop("kappa is not correct.")
  }
  if(rate.scale <= 0){
    stop("rate.scale is not correct.")
  }
  if(!is.null(anc.seq)){
    if(any(!anc.seq %in% .nucleotide$nid[1:4]) || length(anc.seq) != L){
      stop("anc.seq is not correct.")
    }
  }

  ttips <- rooted.tree$Nnode + 1

  ### seqgen uses the order of A, C, G, T for pi.
  ### seqgen uses ts/tv = (kappa*(freqA*freqG + freqC*freqT))/(freqR*freqY)
  ###              freqR = freqA + freqG and freqY = freqC + freqT
  ts.tv <- kappa * (pi[1] * pi[2] + pi[3] * pi[4]) /
                   ((pi[1] + pi[2]) * (pi[3] + pi[4]))

  if(is.null(anc.seq)){
    opts <- paste("-mHKY",
                  " -t", ts.tv,
                  " -f", paste(pi[c(1, 3, 2, 4)], collapse = ","),
                  " -l", L,
                  " -s", rate.scale,
                  " -u", ttips + 1,
                  " -q",
                  sep = "")
    ret <- seqgen(opts = opts, rooted.tree = rooted.tree)
  } else{
    L <- length(anc.seq)
    mu <- paste(nid2code(anc.seq, lower.case = FALSE), collapse = "")
    seqname <- paste("Ancestor  ", collapse = "")
    input <- c(paste(" 1", length(anc.seq), sep = " "),
               paste(seqname, mu, sep = ""),
               1,
               write.tree(rooted.tree, digits = 12))
    opts <- paste("-mHKY",
                  " -t", ts.tv,
                  " -f", paste(pi[c(1, 3, 2, 4)], collapse = ","),
                  " -l", L,
                  " -s", rate.scale,
                  " -u", ttips + 1,
                  " -k1",
                  " -q",
                  sep = "")
    ret <- seqgen(opts, input = input)
  }

  ret
} # End of gen.seq.HKY().


### Generate sequences by the SNP model.
gen.seq.SNP <- function(rooted.tree, pi, L, rate.scale = 1,
    anc.seq = NULL){
  if(!is.rooted(rooted.tree)){
    stop("A rooted tree is required")
  }
  if(length(pi) != 2 || sum(pi) != 1 || any(pi <= 0) || any(pi >= 1)){
    stop("pi is not correct.")
  }
  if(rate.scale <= 0){
    stop("rate.scale is not correct.")
  }
  if(!is.null(anc.seq)){
    if(any(!anc.seq %in% .snp$sid[1:2]) || length(anc.seq) != L){
      stop("anc.seq is not correct.")
    }
  }

  ttips <- rooted.tree$Nnode + 1

  ### Input pi is for "1" and "2".
  ### transfer to "1" = A + G, "2" = C + T
  pi <- c(pi[1], pi[1], pi[2], pi[2]) / 2

  ### seqgen uses the order of A, C, G, T for pi.
  if(is.null(anc.seq)){
    opts <- paste("-mHKY",
                  " -f", paste(pi[c(1, 3, 2, 4)], collapse = ","),
                  " -l", L,
                  " -s", rate.scale,
                  " -u", ttips + 1,
                  " -q",
                  sep = "")
    ret <- seqgen(opts = opts, rooted.tree = rooted.tree)
  } else{
    L <- length(anc.seq)

    ### Replace half of "1" as "A" and the other half as "G"
    ### Replace half of "2" as "C" and the other half as "T"
    ### To avoid bad state of ancestor states in seqgen.
    mu <- paste(snp2code(sid2snp(anc.seq)), collapse = "")
    seqname <- paste("Ancestor  ", collapse = "")
    input <- c(paste(" 1", length(anc.seq), sep = " "),
               paste(seqname, mu, sep = ""),
               1,
               write.tree(rooted.tree, digits = 12))
    opts <- paste("-mHKY",
                  " -f", paste(pi[c(1, 3, 2, 4)], collapse = ","),
                  " -l", L,
                  " -s", rate.scale,
                  " -u", ttips + 1,
                  " -k1",
                  " -q",
                  sep = "")
    ret <- seqgen(opts, input = input)
  }

  ### Replace "A" and "G" by "1", "C" and "T" by "2".
  nseq <- length(ret) - 1
  org <- do.call("rbind", lapply(ret[2:length(ret)], unstrsplit, ""))
  org <- as.matrix(org, nrow = nseq)
  org[, 10 + (1:L)] <- code2snp(org[, 10 + (1:L)])
  org <- apply(org, 1, paste, sep = "", collapse = "")
  ret <- c(ret[1], org)
  attr(ret, "class") <- "seqgen"

  ret
} # End of gen.seq.SNP().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/p_tool_star_tree.r"
### Compute the height of the rooted tree.
get.rooted.tree.height <- function(rooted.tree, tol = .Machine$double.eps^0.5){
  if(!is.rooted(rooted.tree)){
    stop("A rooted tree is required")
  }
  if(!is.ultrametric(rooted.tree, tol = tol)){
    stop("A rooted tree is not ultrametric")
  }

  edge.from <- rooted.tree$edge[, 1]
  edge.to <- rooted.tree$edge[, 2]

  edge.id <- which(edge.from == (length(rooted.tree$tip.label) + 1))
  n.from <- edge.to[edge.id[1]]
  ret <- rooted.tree$edge.length[edge.id[1]]
  repeat{
    edge.id <- which(edge.from == n.from)
    if(length(edge.id) == 2){
      n.from <- edge.to[edge.id[1]]
      ret <- ret + rooted.tree$edge.length[edge.id[1]]
    } else{
      break
    }
  }
  ret
} # End of get.rooted.tree.height().


### Rescale the rooted tree branche lengths by a constant.
rescale.rooted.tree <- function(rooted.tree, scale.height = 1){
  if(!is.rooted(rooted.tree)){
    stop("A rooted tree is required")
  }
  if(scale.height <= 0){
    stop("scale.height is not correct.")
  }

  rooted.tree$edge.length <- rooted.tree$edge.length * scale.height
  rooted.tree
} # End of rescale.tree().


### Reshap a rooted tree to a star shape.
as.star.tree <- function(rooted.tree, keep.bifurcation = TRUE){
  if(!is.rooted(rooted.tree)){
    stop("A rooted tree is required")
  }

  n.tip <- length(rooted.tree$tip.label)
  tree.height <- get.rooted.tree.height(rooted.tree)
  if(keep.bifurcation){
    rooted.tree$edge.length[rooted.tree$edge[, 2] > n.tip] <- 0
    rooted.tree$edge.length[rooted.tree$edge[, 2] <= n.tip] <- tree.height
  } else{
    rooted.tree$edge <- rooted.tree$edge[rooted.tree$edge[, 2] <= n.tip,]
    rooted.tree$edge[, 1] <- max(rooted.tree$edge[, 2]) + 1
    rooted.tree$Nnode <- 1
    rooted.tree$edge.length <- rep(tree.height, n.tip)
  }
  rooted.tree
} # End of as.star.tree().


### Generate a star tree.
gen.star.tree <- function(N, total.height = 1){
  if(N <= 1){
    stop("N > 1.")
  }

  ms.star <- ms(N, 1, opts = paste("-T", sep = " "))
  tree.star <- read.tree(text = ms.star[3])
  tree.star$tip.label <- paste("s", 1:N, sep = "")
  tree.star <- as.star.tree(tree.star)
  th <- get.rooted.tree.height(tree.star)
  ret <- rescale.rooted.tree(tree.star, total.height / th)

  ret
} # End of gen.star.tree().

### Generate a tree with unit heights and K clusters.
gen.unit.K <- function(K, N.K, rate.anc = 10, rate.dec = 10){
  if(K <= 1){
    stop("K > 1.")
  }
  if(any(N.K <= 1)){
    stop("All N.K > 1.")
  }
  if(length(N.K) != K){
    stop("The length of N.K is not equal to K.")
  }
  if(rate.anc <= 0 || rate.dec <= 0){
    stop("rate.anc or rate.dec <= 0.")
  }

  ms.anc <- ms(K, 1, opts = paste("-T -G", rate.anc, sep = " "))
  tree.anc <- read.tree(text = ms.anc[3])
  tree.anc$tip.label <- paste("a", 1:K, sep = "")

  tree.dec <- NULL
  tree.org <- tree.anc
  tree.star <- tree.anc
  for(k in 1:K){
    ms.dec <- ms(N.K[k], 1, opts = paste("-T -G", rate.dec, sep = " "))
    tree.dec[[k]] <- read.tree(text = ms.dec[3])
    tree.dec[[k]]$tip.label <- paste("d", k, ".", 1:N.K[k], sep = "")
    tree.org <- bind.tree(tree.org, tree.dec[[k]],
                          where = which(tree.org$tip.label ==
                                        paste("a", k, sep = "")))
    tree.star <- bind.tree(tree.star, as.star.tree(tree.dec[[k]]),
                           where = which(tree.star$tip.label ==
                                         paste("a", k, sep = "")))
  }

  tree.anc.height <- get.rooted.tree.height(tree.anc)
  tree.dec.height <- unlist(lapply(tree.dec, get.rooted.tree.height))
  tree.dec.mean.height <- mean(tree.dec.height)
  total.height <- tree.anc.height + tree.dec.mean.height
  max.height <- tree.anc.height + max(tree.dec.height)

  tree.equal <- rescale.rooted.tree(tree.anc, 1 / total.height)
  tmp.scale <- tree.dec.mean.height / total.height
  for(k in 1:K){
    tree.dec.equal <- rescale.rooted.tree(tree.dec[[k]],
                                          tmp.scale / tree.dec.height[k])
    tree.equal <- bind.tree(tree.equal, tree.dec.equal,
                            where = which(tree.equal$tip.label ==
                                          paste("a", k, sep = "")))
  }

  tree.max <- rescale.rooted.tree(tree.org, 1 / max.height)
  tree.star <- rescale.rooted.tree(tree.star, 1 / max.height)

  ret <- list(K = K, N.K = N.K, rate.anc = rate.anc, rate.dec = rate.dec,
              height.anc = tree.anc.height,
              height.dec = tree.dec.height,
              anc = tree.anc, dec = tree.dec,
              org = tree.org,
              equal = tree.equal,
              max = tree.max,
              star = tree.star)
  ret
} # End of gen.unit.K().


### Generate star ancestor tree binded with equal star descendant trees.
gen.equal.star.anc.dec <- function(K, N.K, rate.f = 0.5){
  if(length(N.K) != K){
    stop("The length of N.K is not equal to K.")
  }
  if(rate.f <= 0 || rate.f >= 1){
    stop("rate.f <= 0 || >= 1.")
  }
  if(K <= 0 || length(N.K) != K || any(N.K <= 0)){
    stop("K or N.K is not consistant.")
  }

  ms.anc <- ms(K, 1, opts = paste("-T -G", 1, sep = " "))
  tree.anc <- read.tree(text = ms.anc[3])
  tree.anc$tip.label <- paste("a", 1:K, sep = "")

  tree.dec <- NULL
  for(k in 1:K){
    ms.dec <- ms(N.K[k], 1, opts = paste("-T -G", 1, sep = " "))
    tree.dec[[k]] <- read.tree(text = ms.dec[3])
    tree.dec[[k]]$tip.label <- paste("d", k, ".", 1:N.K[k], sep = "")
  }

  tmp.tree <- as.star.tree(tree.anc)
  tmp.height <- get.rooted.tree.height(tmp.tree)
  tree.equal.star <- rescale.rooted.tree(tmp.tree, rate.f / tmp.height)

  for(k in 1:K){
    tmp.tree <- as.star.tree(tree.dec[[k]])
    tmp.height <- get.rooted.tree.height(tmp.tree)
    tmp.tree <- rescale.rooted.tree(tmp.tree, (1 - rate.f) / tmp.height)
    tree.equal.star <- bind.tree(tree.equal.star, tmp.tree,
                                 where = which(tree.equal.star$tip.label ==
                                               paste("a", k, sep = "")))
  }

  ret <- list(K = K, N.K = N.K, rate.f = rate.f,
              anc = tree.anc, dec = tree.dec,
              equal.star = tree.equal.star)
  ret
} # End of gen.equal.anc.dec.star().

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/tzeng-code.getPIstar.getBigMatB.r"


#getcut.fun<-function(pp.org,nn=2*ngeno, plot=0){
#      pp     <-rev(sort(pp.org))
#      ct     <-round(pp*nn)
#      dimen  <-log((1:length(pp)), base=2 )/ct
#      info   <-cumsum(pp*log(1/pp))
#      netinfo<-info-dimen
#      cutpos <-netinfo==max(netinfo)
#      if(plot==1){
#        plot( netinfo, pp,type="b",cex=0.5)
#        barplot(pp.org)
#        abline(h=pp[cutpos])
#      }
#      return( (1:length(pp))[cutpos])
#    }


    final.BigMatB.fun<-function(mmle.cscn, hapreserv, loci, digit){
      ini.cscn     <-get.ini.possHap.hapReserv.fun(mmle.cscn, loci, digit, hapreserv)
      mut.mat.cscn <-ini.cscn$prepsi
      subPI.cscn   <-get.subPI.fun(mmle.cscn,  mut.mat=mut.mat.cscn)
      ## get BB matrix
      jnk.cscn     <-get.preBB.skip.fun(subPI.cscn, loci, digit, mmle.cscn)
      preBB.cscn   <-jnk.cscn$preBB.lst
      indexBB.cscn <-jnk.cscn$indexBB
      medBB.cscn   <-get.medBB.skip.fun(subPI.lst=subPI.cscn, preBB.lst=preBB.cscn, indexBB.lst=indexBB.cscn,loci, digit, mmle.cscn)
      fullBB.cscn  <-get.fullBB.fun(medBB.lst=medBB.cscn,  subPI.lst=subPI.cscn)
      BB.cscn      <-get.BB.fun(fullBB.cscn)
      BigMatB.cscn <-getBigMatB.fun(BB.cscn, subPI.cscn)
      return(list(fullBB.cscn=fullBB.cscn,
                  ##BB.cscn=BB.cscn,
                  BigMatB.cscn =BigMatB.cscn ,
                  subPI.cscn =  subPI.cscn
                  ))
  }



##--------------------------------------------------------------------------------------
## modified in "get.fullBB.fun":    if(aa+bb==2 && is.null(medBB.lst[[rr]][[mm]])==0)
## 
## on 03/10/05: modified in
##   get.preBB.skip.fun: use mmle.cscn to construct matB
##   get.medBB.skip.fun: use mmle.cscn to construct matB
##   getPIstar.fun:      remove the "if" statement
## on 03/11/05:  add a new function:
##   getBigMatB.fun: BB.lst=BB.cs, subPI.lst=subPI.cs
##
## on 04/06/05: add dimensionality check in function
##   getBigMatB.fun, by using chkdim.Multiply.fun, and chkdim.Add.fun
##   modified part is annotated by ******
##
##--------------------------------------------------------------------------------------


##----
get.step.mut.fun<-function(loci, digit,  Ncol, Nrow,NAMEcol, NAMErow){
  ## Ncol  & NAMEcol = col variable
  ## Nrow  & NAMErow = row variable
  UU<-matrix(0, ncol=loci, nrow=Ncol)
  for(i in 1:loci){    UU[,i]<-as.numeric(substring(NAMEcol,   (i-1)*digit+1,i*digit))}
  VV<-matrix(0, ncol=loci, nrow=Nrow)
  for(i in 1:loci){  VV[,i]<-as.numeric(substring(NAMErow,(i-1)*digit+1,i*digit))}
  step.mut<-matrix(0, nrow=Nrow, ncol=Ncol)
  for(v in 1:Nrow){
    tmp<-rep(0,Ncol); for(i in 1:loci){tmp<-tmp+abs(VV[v,i]-UU[,i])}
    step.mut[v,]<-tmp}
  row.names(step.mut)<-NAMErow
  dimnames(step.mut)[[2]]<-NAMEcol
  return(step.mut)}



get.ini.possHap.hapReserv.fun<-function(em, loci, digit, prehapReserv){
  possHap  <-names(em)
  RR       <-length(possHap)
  tmp      <-prehapReserv[  match(possHap, prehapReserv)]
  hapReserv<-tmp[is.na(tmp)==0]
  Rstar     <-length(hapReserv)
  ## cat("RR=", RR, " Rstar=", Rstar,"\n")
  prepsi    <-get.step.mut.fun(loci, digit,  RR, Rstar, possHap, hapReserv )
  return(list(possHap=possHap,
              hapReserv=hapReserv,
              RR=RR,
              Rstar=Rstar,
              prepsi=prepsi))
}

#################################################
## get matrix BB
#################################################
get.subPI.fun<-function(em, mut.mat){
  subPI.lst<-list()
  allhap <-names(em)
  if(nrow(mut.mat)==ncol(mut.mat)){    subPI.lst[[1]]<-em  }else{
    mmin   <-apply(mut.mat, 2, min)
    tab    <-table(mmin)
    subPI.lst<-list()
    subPI.lst[[1]]<-tmp<-em[mmin==0]
    for(rr in 2:(1+max(mmin))){
      del  <-match(names(tmp), names(em))
      subPI.lst[[rr]]<-tmp1<- (em[-del])[apply(matrix(mut.mat[,-del], ncol=(length(em)-length(del)))==(rr-1), 2, sum)>0]
      tmp<-c(tmp1, tmp)
    }
  }
  return(subPI.lst)
}
get.preBB.skip.fun<-function(subPI.lst=subPI.cs, loci, digit, mmle.cscn){
  preBB.lst<-list(); indexBB.lst<-list()
  len      <-length(subPI.lst)
  preBB.lst[[1]]  <-diag(rep(1, length(subPI.lst[[1]])));  dimnames(preBB.lst[[1]])<-list(names(subPI.lst[[1]]),names(subPI.lst[[1]]))
  indexBB.lst[[1]]<-0
  if(len>1)
    {
      sum.subPI.lst<-sapply(subPI.lst, sum)  
      for(rr in 2:len){
        if(sum.subPI.lst[rr]==0){## no "rr"-step descendants
          preBB.lst[[rr]]<-NA
        } else
        if(sum.subPI.lst[rr]> 0)
          {
            index.all<-0; ss<-0
            while(index.all==0)
              {## in order to find the right "ss"-step ancestors to obtain mut.step
                ss<-ss+1
                if(sum.subPI.lst[rr-ss]>0)
                  { 
                    mut.step<-as.data.frame(get.step.mut.fun(loci, digit, length(subPI.lst[[rr-ss]]), length(subPI.lst[[rr]]),
                                                             names(subPI.lst[[rr-ss]]), names(subPI.lst[[rr]])))
                    mmin      <-apply(mut.step, 1, min)
                    index.all<-max(mmin==ss)
                  } else {index.all<-0}
              }
            names(mmin)<-names(subPI.lst[[rr]])
            indexBB.lst[[rr]]<-mmin
            theta<-0
            ## changed on 03/10/2005, should use mmle.cscn instead of pcs or pcn to construct matB
            poolPI<-mmle.cscn[match(names(subPI.lst[[rr-ss]]), names(mmle.cscn))]
            if(length(poolPI)>1){
                tmp <-(theta^(mut.step-ss))%*% diag(poolPI) }else{
                tmp <-(theta^(mut.step-ss)) *       poolPI  }
            preBB.lst[[rr]]<- tmp/    apply(tmp, 1,sum)
            dimnames(preBB.lst[[rr]])<-dimnames(mut.step)
          }
      }## end of rr-loop
    }## end of if(len>1)
  return(list(preBB.lst=preBB.lst,                indexBB.lst=indexBB.lst))
}


##  subPI.lst<-subPI.cs[[j]][[ii]];preBB.lst<-preBB.cs[[j]][[ii]];indexBB.lst<-indexBB.cs[[j]][[ii]]
##  rm(medBB.lst, subPI.lst, preBB.lst, indexBB.lst, len, rr, ss,tt,tmp, tmp1, tmp0, rowsum.tmp,pos,sstop, chk, prechk)

get.medBB.skip.fun<-function(subPI.lst=subPI.cn, preBB.lst=preBB.cn, indexBB.lst=indexBB.cn,loci, digit, mmle.cscn){
  ## if some hap couldn't find a direct 1-step ancestor, then need to use medBB to construct matB
  medBB.lst<-list()
  medBB.lst[[1]]<-list()
  len      <-length(subPI.lst)
  medBB.lst[[1]][[1]]<-preBB.lst[[1]]
  if(len>1)
    {
      for(rr in 2:len)
        {
          medBB.lst[[rr]]<-list()
          if(prod(indexBB.lst[[rr]]==1))## i.e., each hap has its 1-step mutant ancestor
            {
              medBB.lst[[rr]][[1]]<-preBB.lst[[rr]]
            } else
          if(sum(preBB.lst[[rr]]=='NaN')==0)## i.e., from preBB we've found its closet mutant ancestor
            {
              medBB.lst[[rr]][[1]]<-preBB.lst[[rr]]
            }else{## i.e., some hap dont have 1-step ancestors, need to prepare "BB-1step", "BB-2step", ...etc
              ##--- BB-1step ---
              tmp0<-preBB.lst[[rr]];tmp0[tmp0=='NaN']<-0
              medBB.lst[[rr]][[1]]<-tmp0
              chk<-1; ss<-(min(indexBB.lst[[rr]])+1)
              ##for(ss in (min(indexBB.lst[[rr]])+1):(rr-1)){ ##--- BB-"ss"step ---
              while(chk>0 & ss<rr)
                {
                  pos<-(indexBB.lst[[rr]]>=ss )
                  if(sum(pos)==0) ## i.e., all hap have no "ss" step ancetors
                    {
                      medBB.lst[[rr]][[ss]]<-NULL
                    }else{## i.e., some hap have "ss" step ancestors
                      ## need to get the BBmatrix for rr to rr-ss
                      ## note rr=# of mutations needed for current group of hap to become the reserved hap
                      sstop<-0; tt<-ss
                      while(sstop==0)
                        {## in order to find the right "tt"-step ancestors to obtain mut.step
                          if(sum(subPI.lst[[rr-tt]])==0)
                            {
                              tt<-tt+1
                            }else{
                              mut.step<-as.data.frame(get.step.mut.fun(loci, digit, Nrow=length(pos),NAMErow=names(subPI.lst[[rr]]),
                                                                       Ncol=length(subPI.lst[[rr-tt]]),NAMEcol = names(subPI.lst[[rr-tt]])))
                              mmin     <-apply(mut.step, 1, min)
                              sstop   <-prod(mmin[pos]==tt); if(sstop==0){tt<-tt+1}
                              ## here I assume all hap will find their ancestor at the same tt-step, this need to be modified later
                            }
                        }## end of while
                      theta<-0
                      ## changed on 03/10/2005, should use mmle.cscn instead of pcs or pcn to construct matB
                      poolPI<-mmle.cscn[match(names(subPI.lst[[rr-tt]]), names(mmle.cscn))]
                      if(length(poolPI)>1)
                        {
                          tmp<-(theta^(mut.step-tt))%*% diag(poolPI)
                        }else{
                          tmp<-(theta^(mut.step-tt))  *      poolPI
                        }
                      rowsum.tmp<-apply(tmp, 1,sum); rowsum.tmp[rowsum.tmp==0]<-1
                      tmp1<- tmp/rowsum.tmp
                      medBB.lst[[rr]][[ss]]<-tmp1*0
                      medBB.lst[[rr]][[ss]][pos,]<-tmp1[pos,]
                      dimnames(medBB.lst[[rr]][[ss]])<-dimnames(mut.step)
                    }
                  ## check if each hap finds its ancestor
                  if(length(medBB.lst[[rr]])>1){
                    tmp<-lapply(medBB.lst[[rr]],function(mat){if(1-is.null(mat)){apply(mat, 1, sum)}else{NA}})
                    prechk<-sapply(tmp, sum)
                  }else{  prechk<-apply(medBB.lst[[rr]][[1]],1, sum)}
                  chk<-sum(prechk==0, na.rm=T)   ## chk=0 means all category can find its closest ancestor; chk>0 means some havent found yet.
                  if(chk>0){ss<-ss+1}
                }## end of while
            }
        }## end of rr-loop
    }
  return(medBB.lst)
}




get.fullBB.fun<-function(medBB.lst=medBB.cs,subPI.lst=subPI.cs)
{
  name.lst <-lapply(subPI.lst, names)
  dim.vec  <-sapply(subPI.lst, length)
  len      <-length(dim.vec)

  fullBB.lst<-list()
  fullBB.lst[[1]]<-medBB.lst[[1]]
  if(len>1)
    {
      fullBB.lst[[2]]<-medBB.lst[[2]]
      if(len>2)
        {
          for(rr in 3:length(medBB.lst))
            {
              fullBB.lst[[rr]]<-list()
              for(kk in (1:(rr-1)))
                {
                  mat<-matrix(0, nrow=dim.vec[rr], ncol=dim.vec[rr-kk])
                  row.names(mat)    <-name.lst[[rr]]
                  dimnames(mat)[[2]]<-name.lst[[rr-kk]]
                  fullBB.lst[[rr]][[kk]]<-mat
                }
              
              for(kk in (1:(rr-1)))
                {
                  for(mm in 1:length(medBB.lst[[rr]]))
                    {
                      aa<-sum(is.na(match(dimnames(fullBB.lst[[rr]][[kk]])[[2]], dimnames( medBB.lst[[rr]][[mm]])[[2]] )))==0
                      bb<-sum(is.na(match(dimnames( medBB.lst[[rr]][[mm]])[[2]], dimnames(fullBB.lst[[rr]][[kk]])[[2]] )))==0
                      if(aa+bb==2 && is.null(medBB.lst[[rr]][[mm]])==0)
	              ##   if(aa+bb==2)
                        {
                          fullBB.lst[[rr]][[kk]]<-medBB.lst[[rr]][[mm]]
                        }
                    }## end of mm-loop
                }## end of kk-loop

              
            }## end of rr-loop
        }## end of if(len>2)
    }## end of if (len>1)
  return(fullBB.lst)

}


get.BB.fun<-function(fullBB.lst){
  BB.lst  <-list()
  len     <-c(1,(1:100))
  cumlen  <-cumsum(len)
  BB.lst[[1]]<-fullBB.lst[[1]][[1]]
  if(length(fullBB.lst)>1)
    {
      BB.lst[[ 2 ]]<-fullBB.lst[[2]][[1]]
      if(length(fullBB.lst)>2)
        {
          for(rr in 3:length(fullBB.lst))
            {
              for(kk in 1:length(fullBB.lst[[rr]]))
                {
                  BB.lst[[ cumlen[rr-1]+kk ]]<-fullBB.lst[[rr]][[kk]]
                }
            }
        }
    }
  return(BB.lst)
}


  
  

  
get.allpossROUTE.fun<-function()
{
  ### when this file is modified, function "getPIstar.fun" need to be changed as well!!!!!!
  m1 <-list(1)
  m2 <-  lapply(m1, 2, FUN = c)
  m3 <-c(lapply(m2, 3, FUN = c),
         lapply(m1, 4, FUN = c))
  m4 <-c(lapply(m3, 5, FUN = c),
         lapply(m2, 6, FUN = c),
         lapply(m1, 7, FUN = c))
  m5 <-c(lapply(m4, 8, FUN = c),
         lapply(m3, 9, FUN = c),
         lapply(m2, 10, FUN = c),
         lapply(m1, 11, FUN = c))
  m6 <-c(lapply(m5, 12, FUN = c),
         lapply(m4, 13, FUN = c),
         lapply(m3, 14, FUN = c),
         lapply(m2, 15, FUN = c),
         lapply(m1, 16, FUN = c))
  m7 <-c(lapply(m6, 17, FUN = c),
         lapply(m5, 18, FUN = c),
         lapply(m4, 19, FUN = c),
         lapply(m3, 20, FUN = c),
         lapply(m2, 21, FUN = c),
         lapply(m1, 22, FUN = c))
  m8 <-c(lapply(m7, 23, FUN = c),
         lapply(m6, 24, FUN = c),
         lapply(m5, 25, FUN = c),
         lapply(m4, 26, FUN = c),
         lapply(m3, 27, FUN = c),
         lapply(m2, 28, FUN = c),
         lapply(m1, 29, FUN = c))
  m9 <-c(lapply(m8, 30, FUN = c),
         lapply(m7, 31, FUN = c),
         lapply(m6, 32, FUN = c),
         lapply(m5, 33, FUN = c),
         lapply(m4, 34, FUN = c),
         lapply(m3, 35, FUN = c),
         lapply(m2, 36, FUN = c),
         lapply(m1, 37, FUN = c))
  m10<-c(lapply(m9, 38, FUN = c),
         lapply(m8, 39, FUN = c),
         lapply(m7, 30, FUN = c),
         lapply(m6, 41, FUN = c),
         lapply(m5, 42, FUN = c),
         lapply(m4, 43, FUN = c),
         lapply(m3, 44, FUN = c),
         lapply(m2, 45, FUN = c),
         lapply(m1, 46, FUN = c))
  
  m = c(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10)
  d = 10
  r = c(1, rep(2, d-1))^(c(0,0:(d-2)))
  di = d * (d-1) / 2 + 1
  
  l = rep(1:d, r)
  mat = matrix(0, sum(r), di)
  for(i in 1:sum(r)) {
    mat[i,m[[i]]] = m[[i]]
  }
  row.names(mat)<-paste("Step ",
                        c(1,2,rep(3,2), rep(4,4), rep(5,8), rep(6,16), rep(7,32), rep(8,64), rep(9,128), rep(10, 256)),
                        ".",
                        c(1,1,   (1:2),    (1:4),    (1:8),    (1:16),    (1:32),    (1:64),    (1:128), (1:256)),   sep="")
  dimnames(mat)[[2]]<-paste("mat",c(1,2,rep(3,2), rep(4,3), rep(5,4), rep(6,5),rep(7,6)),".(",(1:di),")",sep="")
  return(mat)
}


    

      chkdim.Multiply.fun<-function(mat){
        chk1<-length(mat)>0
        chk2<-all(dim(mat)>0)
        chk3<-sum(is.na(mat))==0
        return(list("chk1"=chk1,"chk2"=chk2, "chk3"=chk3))
      }


      chkdim.Add.fun<-function(mat, ncolShouldBe=RRstar){
        if(length(mat)==0)
          {
            chk1=chk2=0
          }else{
            chk1<-all(dim(mat)>0)
            chk2<-(ncol(mat) == ncolShouldBe)
          }
        return(chk1*chk2)
      }

getBigMatB.fun<-function(BB.lst=BB.cs, subPI.lst=subPI.cs){
   ## modified on Apr 6, 2005; now check for the compatibility of dimensionality before adding or mulitplying	
   ## see ****** part. functions used:  chkdim.Multiply.fun, and  chkdim.Add.fun
    mmat         <-get.allpossROUTE.fun()
    r.ele        <-c( 1,1,2,4,8,16,32,64,128,256)
    cumsum.r     <-cumsum(r.ele)
    lenPI        <-length(subPI.lst)  # 4
    lenBB        <-length(BB.lst)     # 7=1+1+2+3
    lenROUTE     <-cumsum.r[lenPI]    # 8=1+1+2+4
    RRstar       <-length(subPI.lst[[1]]) ##*********
    
    ## changed on 03/11/05 (remove "if" condition)
    mat.route<-as.matrix(mmat[(1:lenROUTE),(1:lenBB)])
    pi.ele   <-rep((1:10), r.ele)
    j.index  <-pi.ele[1:cumsum.r[lenPI]] ## indexing a mat B is for the level of H^{(j)}

    preBigBB.lst<-list()
    for(kk in 1:lenROUTE)## it's fine to have kk start from "1"; when lenROUTE=1, BB=bb[[1]]=Identity Matrix
      {
        bb<-BB.lst[mat.route[kk,]]
        BB<-bb[[length(bb)]]
        if(length(bb)>1)
          {
            ##***** check the compability of dimenationality  before multiplying
            chk<-sapply(bb, chkdim.Multiply.fun)
            if(all(chk==T))
              {
                for(jj in max(1,(length(bb)-1)):1){BB<-BB%*%bb[[jj]]}### range of jj is modified on 03/12/05 as well
              }
            else
              {
                BB<-matrix(0, nrow=length(subPI.lst[[ j.index[[kk]]  ]]), ncol=RRstar)
                colnames(BB) = names(subPI.lst[[1]])
                rownames(BB) = names(subPI.lst[[ j.index[[kk]]  ]])
              }
          }
        preBigBB.lst[[kk]]<-BB
      }
    
    indx<-split((1:length(preBigBB.lst)), j.index)
    BigMatB<-NULL
    for(jj in 1:length(indx))
      {
        tmpB<-preBigBB.lst[ indx[[jj]] ]
        if(length(tmpB)==1)
          {
            BigMatB<-rbind(BigMatB, tmpB[[1]])
          }else{
            ##***** check the compability of dimenationality  before adding
           key<-sapply(tmpB,chkdim.Add.fun, ncolShouldBe=RRstar); pos<-(1:length(tmpB))[key==1]
            if(all(pos==0))
              {
                sumB = NULL
              } else {
                Add.tmpB<-tmpB[pos]
                sumB    <-Add.tmpB[[1]]
                for(ll in 2:length(Add.tmpB)){sumB<-sumB+Add.tmpB[[ll]]}
              }
            BigMatB<-rbind(BigMatB, sumB)
          }
      }
    return(BigMatB)
  }


##   getBigMatB.fun<-function(BB.lst=BB.cs, subPI.lst=subPI.cs)
##     {
##       mmat         <-get.allpossROUTE.fun()
##       r.ele        <-c( 1,1,2,4,8,16,32,64,128,256)
##       cumsum.r     <-cumsum(r.ele)
##       lenPI        <-length(subPI.lst)  # 4
##       lenBB        <-length(BB.lst)     # 7=1+1+2+3
##       lenROUTE     <-cumsum.r[lenPI]    # 8=1+1+2+4
##   
##       ## changed on 03/11/05 (remove "if" condition)
##       mat.route<-as.matrix(mmat[(1:lenROUTE),(1:lenBB)])
##       pi.ele   <-rep((1:10), r.ele)
##       j.index  <-pi.ele[1:cumsum.r[lenPI]] ## indexing a mat B is for the level of H^{(j)}
##   
##       preBigBB.lst<-list()
##       for(kk in 1:lenROUTE)
##         {
##           bb<-BB.lst[mat.route[kk,]]
##   	## in "getPIstar.fun" there is such if, but it seems to be useless so i removed
##   	## on 03/10/05
##           ## if(sum(is.na(bb))>0){
##           ##   ##          prePIstar.lst[[kk]]<-rep(0, length(subPI.lst[[1]]))
##           ##   cat("dd=",dd); print(" is.na(bb)>0")
##           ## }else{ below...}
##           BB<-bb[[length(bb)]]
##           for(jj in max(1,(length(bb)-1)):1){BB<-BB%*%bb[[jj]]}
##           preBigBB.lst[[kk]]<-BB
##         }
##       indx<-split((1:length(preBigBB.lst)), j.index)
##       BigMatB<-NULL
##       for(jj in 1:length(indx))
##         {
##           tmpB<-preBigBB.lst[ indx[[jj]] ]
##           if(length(tmpB)==1)
##             {
##               BigMatB<-rbind(BigMatB, tmpB[[1]])
##             }else{
##               sumB<-0
##               for(ll in 1:length(tmpB)){sumB<-sumB+tmpB[[ll]]}
##               BigMatB<-rbind(BigMatB, sumB)
##             }
##         }
##       return(BigMatB)
##     }
##   
##   
##   



##rm(BB.lst, subPI.lst, prePIstar.lst, mmat, r.ele, cumsum.r, lenPI, lenBB, lenROUTE, mat.route, pi.ele, PI,bb, BB, kk, jj,PIstar)
getPIstar.fun<-function(BB.lst=BB.cs, subPI.lst=subPI.cs)
  {
    mmat         <-get.allpossROUTE.fun()
    r.ele        <-c( 1,1,2,4,8,16,32,64,128,256)
    cumsum.r     <-cumsum(r.ele)
    lenPI        <-length(subPI.lst)  # 4
    lenBB        <-length(BB.lst)     # 7=1+1+2+3
    lenROUTE     <-cumsum.r[lenPI] # 8=1+1+2+4

    ## changed on 03/11/05 (remove "if" condition)
    mat.route<-as.matrix(mmat[(1:lenROUTE),(1:lenBB)])
    pi.ele   <-rep((1:10), r.ele)
    PI       <-subPI.lst[pi.ele[1:cumsum.r[lenPI]]]

    prePIstar.lst<-list()
    for(kk in 1:lenROUTE)
      {
        bb<-BB.lst[mat.route[kk,]]
        if(sum(is.na(bb))>0){
          prePIstar.lst[[kk]]<-rep(0, length(subPI.lst[[1]]))
        }else{
          BB<-bb[[length(bb)]]
          for(jj in max(1,(length(bb)-1)):1){BB<-BB%*%bb[[jj]]}
          prePIstar.lst[[kk]]<-PI[[kk]] %*% BB
        }
      }

    PIstar<-colSums(matrix(unlist(prePIstar.lst), nrow=length(prePIstar.lst), byrow=T))
    names(PIstar)<-names(subPI.lst[[1]])
    return(PIstar)
  }

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/tzeng-functions.r"
##-------------------------------------------------------------
## modified getI13.fun on Apr3, 2005; use it to replace I13.fun; in the new version, it can be
##                       used for both phased or unphased data
##-------------------------------------------------------------
## modified get.deriBigMatB.fun on Apr8, 2005; add dimensionality check
##   using chkdim.Multiply.fun
##         chkdim.Add.fun
##-------------------------------------------------------------


get.I.fun = function(y, Xt.mat, bigB.dPI, subPI.lst) {
  ## this function calculates I13 %*% I33i %*% t(I13), with input being:
  ## y:          data, length [nsim]
  ## Xt.mat:     data; dim [nsim, RR]
  ## bigB.dPI:   derivative info, calculated using construct.bigB.dPI,
  ##               and see generalform.tex for explaination;
  ## subPI.lst:  list of PIs with frequencies and names listed
  I13       = getI13.fun(y, Xt.mat, bigB.dPI)
  I13       = I13[-1,]  ## the first row is taken off, the base-line row. 
  I33i      = I33Inv.fun(subPI.lst, Xt.mat)
  I13t      = t(I13)
  Imat      = I13 %*% I33i %*% I13t
  ans       = list(I13 = I13, I33i = I33i, Imat = Imat)
  return(ans)
}


I33Inv.fun = function(subPI.lst, Xt.mat)
{
  ## this calculates I33-inverse matrix, and the names follow the same
  ## order as subPI.lst, and not Xt.mat. So Xt.mat's names have to agree
  ## with subPI.lst;
  
  ## subPI.lst: list of PIs with frequencies and names listed
  ## Xt.mat:    design matrix, of [nsim, RR]

  a = unlist(subPI.lst)
  a = a^2
  
  a.name = names(a)
  x.name = colnames(Xt.mat)
  b      = match(a.name, x.name)
  RR     = length(a.name)
  if(sum(b == (1:RR)) != RR)
    print("colnames(Xt.mat) and names(unlist(subPI.lst))) do not match - I33Inv.fun");
  Xt.mat = Xt.mat[,b]

  X = colSums(Xt.mat)
  ans = diag(a / (2 * X))
  rownames(ans) = a.name
  colnames(ans) = a.name

  return(ans)
}

       
getI13.fun = function(t1=y, Xt.mat, bigB.dPI) {
  ## This function returns the I_{13} matrix
  ## Modified Apr3,2005; change y-ybar to t1, where t1= (y - mean(y))/a for phased data
  ##                                                t1= rep((y-mu)/a, nreps) * post for unphased data
  ## IMPORTANT!!!!!!!
  ## colnames(Xt.mat) should agree with rownames(bigB.dPI[[1]]), and
  ##  Xt.mat is sorted according to the rownames of bigB.dPI, so the
  ##  return matrix, I13's dimension names agree with bigB.dPI[[i]], and
  ##  not with colnames(Xt.mat)
  
  ## y is of length(n);
  ## Xt.mat is of dimension [n, RR]
  ## bigB.dPI is a list of length [RR], with each element being a
  ##   matrix of dimension [RR, RRstar]
  
  n      = length(t1)
  RR     = nrow(bigB.dPI[[1]])     ## original PI dimension;
  RRstar = ncol(bigB.dPI[[1]])     ## reduced PI dimension;
  I13    = matrix(0, RRstar, RR)   ## dimension of I13;
  o.name = rownames(bigB.dPI[[1]]) ## original PI names
  r.name = colnames(bigB.dPI[[1]]) ## reduced PI names
  
  ## check if o.name and Xt.mat names agree, and this sorts Xt.mat regardlessly
  x.name = colnames(Xt.mat)
  b      = match(o.name, x.name)
  if(sum(b == (1:RR)) != RR){    print("colnames(Xt.mat) and rownames(bigB.dPI[[1]]) do not match - getI13.fun")}
  Xt.mat = Xt.mat[,b]
  yX    = colSums( t1 * Xt.mat)  ## this should be the same as using ymat
  for(p in 1:RRstar)
    {
      for(k in 1:RR)
        {
          ## bigB.dPI[[k]][h,p] ==>
          ## k: pi that bigB.dPI is taking derivative w.r.t. to;
          ## h: pi that is in the original PI configuration;
          ## p: pi that is in the reduced PI configuration;
          aa = yX %*% bigB.dPI[[k]][,p];
          I13[p, k] = aa
        }
    }
  colnames(I13) = o.name;
  rownames(I13) = r.name;
  return(I13)
}


 ## I13.fun = function(y, Xt.mat, bigB.dPI) {
 ##   ## This function returns the I_{13} matrix
 ## 
 ##   ## IMPORTANT!!!!!!!
 ##   ## colnames(Xt.mat) should agree with rownames(bigB.dPI[[1]]), and
 ##   ##  Xt.mat is sorted according to the rownames of bigB.dPI, so the
 ##   ##  return matrix, I13's dimension names agree with bigB.dPI[[i]], and
 ##   ##  not with colnames(Xt.mat)
 ## 
 ##   ## y is of length(n);
 ##   ## Xt.mat is of dimension [n, RR]
 ##   ## bigB.dPI is a list of length [RR], with each element being a
 ##   ##   matrix of dimension [RR, RRstar]
 ##   
 ##   n      = length(y)
 ##   RR     = nrow(bigB.dPI[[1]])     ## original PI dimension;
 ##   RRstar = ncol(bigB.dPI[[1]])     ## reduced PI dimension;
 ##   I13    = matrix(0, RRstar, RR)   ## dimension of I13;
 ##   ybar   = mean(y)
 ##   o.name = rownames(bigB.dPI[[1]]) ## original PI names
 ##   r.name = colnames(bigB.dPI[[1]]) ## reduced PI names
 ## 
 ##   ## check if o.name and Xt.mat names agree, and this sorts Xt.mat
 ##   ## regardlessly
 ##   x.name = colnames(Xt.mat)
 ##   b      = match(o.name, x.name)
 ##   if(sum(b == (1:RR)) != RR)
 ##     print("colnames(Xt.mat) and rownames(bigB.dPI[[1]]) do not match - I13.fun");
 ##   Xt.mat = Xt.mat[,b]
 ##   
 ##   ## ymat = matrix( rep( y - ybar, RR), ncol = RR)
 ##   ## yX   = colSums(ymat * Xt.mat)
 ##   yX    = colSums( (y - ybar) * Xt.mat)  ## this should be the same as using ymat
 ## 
 ##   for(p in 1:RRstar)
 ##     {
 ##       for(k in 1:RR)
 ##         {
 ##           ## bigB.dPI[[k]][h,p] ==>
 ##           ## k: pi that bigB.dPI is taking derivative w.r.t. to;
 ##           ## h: pi that is in the original PI configuration;
 ##           ## p: pi that is in the reduced PI configuration;
 ##           a = yX %*% bigB.dPI[[k]][,p];
 ##           I13[p, k] = a
 ##         }
 ##     }
 ##   colnames(I13) = o.name;
 ##   rownames(I13) = r.name;
 ##   return(I13)
 ## }



zero.mat.fun = function(mat) {
  ## making the input matrix zero in every cell;
  nr = nrow(mat);
  nc = ncol(mat);
  rn = rownames(mat);
  cn = colnames(mat);
  mat = matrix(0, nr, nc);
  rownames(mat) = rn;
  colnames(mat) = cn;
  return(mat)
}

index.k.fun<-function(len.subPI = 6) {
  ## this function creats a bit strange index that keeps track of
  ## B matrices in the BB-list, with each B indicates how it transforms
  ## layers of haplotypes to their appropriate ancestors;

  ##             __[[1]]  
  ##            |   
  ##  pi(0)     1   __[[2]]
  ##  pi(1)     2  |   __[[3]]
  ##  pi(2)     4  3  |   __[[4]]
  ##  pi(3)     7  6  5  |   __[[5]] 
  ##  pi(4)    11 10  9  8  |   __[[6]]
  ##  pi(5)    16 15 14 13 12  |   __[[7]]
  ##  pi(6)    22 21 20 19 18 17  |   __[[8]]
  ##  pi(7)    29 28 27 26 25 24 23  |
  ##  pi(8)    37 36 35 34 33 32 31 30

  ## the numbers indicate which member of the BB-list belongs to which layer of pi;
  ## [[num]] is the return list in which the vertical numbers are sorted;
  ## NOTE: because pi(num) starts at 0, so len.subPI = 6 ==> pi(5)!!!!
  ## NOTE: the last [[num]] is always set to be NA!!!
  
  sseq <- list()
  KK <- len.subPI - 1
  if(KK == 0)
    {
      sseq[[1]] <- 1
    } else {
      sseq[[1]] <- 1 + cumsum(1:KK)
      if(KK > 1)
        {
          for(k in 2:KK) {
            sseq[[k]] <- sseq[[k-1]][-1]-1
          }
        }
      sseq[[1]] = c(1, sseq[[1]])
      sseq[[(KK+1)]] = NA
    }
  return(sseq)
}


##--------
      chkdim.Multiply.fun<-function(mat){
        chk1<-length(mat)>0
        chk2<-all(dim(mat)>0)
        chk3<-sum(is.na(mat))==0
        return(list("chk1"=chk1,"chk2"=chk2, "chk3"=chk3))
      }

      chkdim.Add.fun<-function(mat, ncolShouldBe=RRstar){
        if(length(mat)==0)
          {
            chk1=chk2=0
          }else{
            chk1<-all(dim(mat)>0)
            chk2<-(ncol(mat) == ncolShouldBe)
          }
        return(chk1*chk2)
      }


get.deriBigMatB.fun<-function(BB.deri, subPI.lst, subPi.num=k)
  { 
    ## modified on Apr 8 2005, add the compatability check of dimensionalith when getting
    ## preBigBB.deri and BigMatB, see "*******".
    ##
    ## this calculate one bigB.dPI (w.r.t. to one pi). For reference,
    ## see generalform.tex;

    ## NOTE: this function is not easy to use, because it does not say,
    ##  specifically which pi that the function is calculating w.r.t.;
    ##  Instead, one has to know how to construct the correct BB.deri,
    ##  which should have the appropriate cells filled out (see
    ##  construct.bigB.dPI) by derivative information. 

    ## BB.deri: It's a BB-list, but for this function, it should have
    ##         appropriate cells filled out with derivative calculated
    ##         information.
    ## subPI.lst:  list of PIs with frequencies and names listed
    ## subPi.num:  the number of "sub-group of PIs" that the derivative
    ##              pi belongs. The number should not exceed length(subPI.lst); 
    
    mmat         <-get.allpossROUTE.fun()
    r.ele        <-c( 1,1,2,4,8,16,32,64,128,256)
    cumsum.r     <-cumsum(r.ele)
    lenPI        <-length(subPI.lst)  ## 4
    lenBB        <-length(BB.deri)     ## 7=1+1+2+3
    lenROUTE     <-cumsum.r[lenPI]    ## 8=1+1+2+4
    RRstar       <-length(subPI.lst[[1]])
    ### changed on 03/11/05 (remove "if" condition)
    mat.route<-as.matrix(mmat[(1:lenROUTE),(1:lenBB)]);
    
    mat.route[,1]<-0 ### this adjustment (made on 03/12/05 is for calculating
                     ### deri.B only, because BB.deri[[1]] is alwyas a 0 matrix.
                     ### Hence has to be exclude from the product calculation,
                     ### o.w. making the final product=0
    pi.ele  <-rep((1:length(r.ele)), r.ele)
    j.index <-pi.ele[1:cumsum.r[lenPI]] ###  indexing a mat B is for the level of H^{(j)}

    preBigBB.deri<-list()
    preBigBB.deri[[1]]<-BB.deri[[1]]
    if(lenROUTE>1)
      {
        for(kk in 2:lenROUTE)
          {
            bb<-BB.deri[mat.route[kk,]]
            BB<-bb[[length(bb)]]
            if(length(bb)>1)
              {
                ##***** check the compability of dimenationality  before multiplying
                chk<-sapply(bb, chkdim.Multiply.fun)
                if(all(chk==T))
                  {
                    for(jj in max(1,(length(bb)-1)):1){BB<-BB%*%bb[[jj]]}### range of jj is modified on 03/12/05 as well
                  }
                else
                  {
                    BB<-matrix(0, nrow=length(subPI.lst[[ j.index[[kk]]  ]]), ncol=RRstar)
                    colnames(BB) = names(subPI.lst[[1]])
                    rownames(BB) = names(subPI.lst[[ j.index[[kk]]  ]])
                  }
            } 
            preBigBB.deri[[kk]]<-BB
          }
      }

    ## If subPi.num is not in the first (reduced group),
    ##  then some of the preBigBB.deri matrices have to be zeroed.     
    index = index.k.fun(lenPI)
    if(subPi.num > 1)
      {
        for(kk in 1:(subPi.num-1))
          {
            index.k = index[[kk]]
            for(jj in index.k)
              {
                preBigBB.deri[[jj]] = zero.mat.fun(preBigBB.deri[[jj]])
              }
          }
      }
    
    indx<-split((1:length(preBigBB.deri)), j.index)
    BigMatB<-NULL
    for(jj in 1:length(indx))
      {
        tmpB<-preBigBB.deri[ indx[[jj]] ]
        if(length(tmpB)==1)
          {
            BigMatB<-rbind(BigMatB, tmpB[[1]])
          }else {
            ##***** check the compability of dimenationality  before adding
            key<-sapply(tmpB,chkdim.Add.fun, ncolShouldBe=RRstar); pos<-(1:length(tmpB))[key==1]
            if(all(pos==0))
              {
                sumB = NULL
              } else {
                Add.tmpB<-tmpB[pos]
                sumB    <-Add.tmpB[[1]]
                for(ll in 2:length(Add.tmpB)){sumB<-sumB+Add.tmpB[[ll]]}
              }
            BigMatB<-rbind(BigMatB, sumB)
          }
      }
    return(BigMatB)
  }
            
## old version 
## get.deriBigMatB.fun<-function(BB.lst, subPI.lst, subPi.num=k)
##   { 
##     ## this calculate one bigB.dPI (w.r.t. to one pi). For reference,
##     ## see generalform.tex;
## 
##     ## NOTE: this function is not easy to use, because it does not say,
##     ##  specifically which pi that the function is calculating w.r.t.;
##     ##  Instead, one has to know how to construct the correct BB.lst,
##     ##  which should have the appropriate cells filled out (see
##     ##  construct.bigB.dPI) by derivative information. 
## 
##     ## BB.lst: It's a BB-list, but for this function, it should have
##     ##         appropriate cells filled out with derivative calculated
##     ##         information.
##     ## subPI.lst:  list of PIs with frequencies and names listed
##     ## subPi.num:  the number of "sub-group of PIs" that the derivative
##     ##              pi belongs. The number should not exceed length(subPI.lst); 
##     
##     mmat         <-get.allpossROUTE.fun()
##     r.ele        <-c( 1,1,2,4,8,16,32,64,128,256)
##     cumsum.r     <-cumsum(r.ele)
##     lenPI        <-length(subPI.lst)  ## 4
##     lenBB        <-length(BB.lst)     ## 7=1+1+2+3
##     lenROUTE     <-cumsum.r[lenPI]    ## 8=1+1+2+4
## 
##     ### changed on 03/11/05 (remove "if" condition)
##     mat.route<-as.matrix(mmat[(1:lenROUTE),(1:lenBB)]);
## 
##     
##     mat.route[,1]<-0 ### this adjustment (made on 03/12/05 is for calculating
##                      ### deri.B only, because BB.lst[[1]] is alwyas a 0 matrix.
##                      ### Hence has to be exclude from the product calculation,
##                      ### o.w. making the final product=0
##     pi.ele<-rep((1:10), r.ele)
##     j.index <-pi.ele[1:cumsum.r[lenPI]] ###  indexing a mat B is for the level of H^{(j)}
## 
##     preBigBB.lst<-list()
##     preBigBB.lst[[1]]<-BB.lst[[1]]
##     if(lenROUTE>1)
##       {
##         for(kk in 2:lenROUTE)
##           {
##             bb<-BB.lst[mat.route[kk,]]
##             BB<-bb[[length(bb)]]
##             if(length(bb)>1){ for(jj in max(1,(length(bb)-1)):1){BB<-BB%*%bb[[jj]]} } ### range of jj is modified on 03/12/05 as well
##             preBigBB.lst[[kk]]<-BB
##           }
##       }
## 
##     ## If subPi.num is not in the first (reduced group),
##     ##  then some of the preBigBB.lst matrices have to be zeroed.     
##     index = index.k.fun(lenPI)
##     if(subPi.num > 1)
##       {
##        for(kk in 1:(subPi.num-1))
##          {
##            index.k = index[[kk]]
##            for(jj in index.k)
##              {
##                preBigBB.lst[[jj]] = zero.mat.fun(preBigBB.lst[[jj]])
##              }
##          }
##       }
## 
## 
##     indx<-split((1:length(preBigBB.lst)), j.index)
##     BigMatB<-NULL
##     for(jj in 1:length(indx))
##       {
##         tmpB<-preBigBB.lst[ indx[[jj]] ]
##         if(length(tmpB)==1)
##           {
##             BigMatB<-rbind(BigMatB, tmpB[[1]])
##           }else{
##             sumB<-0
##             for(ll in 1:length(tmpB)){sumB<-sumB+tmpB[[ll]]}
##             BigMatB<-rbind(BigMatB, sumB)
##           }
##       }
##     return(BigMatB)
##   }


  construct.dB.list.fun = function(CC.mat, subPI.lst) {
    ## modified on Apr 10, 2005, see **** parts to take care of "NA" and "numeric(0)" in CC.mat
    ##
    ## this constructs the dB.list, see reference in generalform.tex for
    ## every CC.mat, calculate, w.r.t. every PI in that CC.mat,
    ## partial(B)/partial(pi). See pt1 -- pt6 calculations in
    ## generalform.tex. The answer it returns is a list of RR (length of
    ## subPI.lst).
    
    ## Note: Because the input CC.mat is of different dimension, the
    ##        output is a list of RR matrices, each with the same dim
    ##        as the input CC.mat matrix.
    ## Note: The derivative is taken with respect to ALL pi's, and if
    ##        any pi is not part of the CC.mat, then the answers are
    ##        automatically zero.
    ## Note: If the CC.mat is all zero, the outcome is automatically
    ##        all ZERO. 
    
    ## CC.mat is ONE CC matrix;
    ## subPI.lst: list of PIs with frequencies and names listed

    ##********* begin modificiation on 04.10.2005 ****************
    if(sum(is.na(CC.mat)) > 0  | length(CC.mat)==0 )
      {
        ## when CC.mat is NA or numeric(0), make dB.list=NULL
          full.PI  = unlist(lapply(subPI.lst, names))
          l        = length(full.PI)
          dB.list  = list()
          for(i in 1:l)
            {
              dB.list[[i]] = numeric(0)
              attr(dB.list[[i]], "name") = full.PI[i]
            }
          ##********* begin modificiation on 04.10.2005 ****************
      } else if(all(CC.mat == 0))
        {
          full.PI  = unlist(lapply(subPI.lst, names))
          l        = length(full.PI)
          dB.list  = list()
          for(i in 1:l)
            {
              dB.list[[i]] = CC.mat * 0
              attr(dB.list[[i]], "name") = full.PI[i]
            }
        } else {
          target.PI = colnames(CC.mat)  #pi^C in the doc
          full.PI = unlist(lapply(subPI.lst, names))
          
          for(i in 1:length(subPI.lst))
            {
              if(!is.na(sum(match(target.PI, names(subPI.lst[[i]])))))
                {
                  k = i;
                  break;
                }
            }
          subPI = subPI.lst[[k]]         ##freq of target.PI
          
          ##pt1 = pt2 = pt3 = pt4 = pt5 = pt6 = CC.mat;
          
          nr  = nrow(CC.mat)
          nc  = ncol(CC.mat)
          
          pt5 = matrix(rep(subPI, nr), nrow = nr, byrow=T)
          
          deno = rowSums(CC.mat * pt5)
          a    = 1/deno;       a[deno==0] = 0
          pt1 = matrix(rep(a, nc), ncol = nc)
          
          pt3 = pt1^2
          pt4 = CC.mat
          
          l = length(full.PI)
          dB.list = list()
          
          ## dB.list is a list of length l, with each element being a
          ## derivative taken with respect to the name attributes. Therefore,
          ## if the name is not an element of target.PI, that dB.list is left
          ## to be zero. 
          
          for(i in 1:l)
            {
              dB.list[[i]] = matrix(0, nr, nc)
              dB.list[[i]] = CC.mat * 0
              attr(dB.list[[i]], "name") = full.PI[i]
            }
          
          no.zero.list = match(full.PI, target.PI)
          
          for(j in 1:l)
            {
              if(!is.na(no.zero.list[j]))
                {
                  i = no.zero.list[j]
                  pt2 = matrix(0, nr, nc)
                  pt2[,i] = CC.mat[,i]
                  
                  pt6 = matrix(rep(pt2[,i], nc), ncol = nc)
                  dB.list[[j]] = pt1 * pt2 - pt3 * pt4 * pt5 * pt6
                  attr(dB.list[[j]], "name") = full.PI[j]
                }
            }
        }
    
    return(dB.list)
  }


construct.bigB.dPI = function(dB.all, subPI.lst, BB.lst) 
{
  ## this function constructs the (bigB.dPI) expression in
  ## generalform.tex. In reality, it calculates derivative of B
  ## w.r.t. to ALL PI's. So for RR-pi, the function returns a list of
  ## length(RR), each is a matrix of the same dimension as BB.

  ## Note: this function calls get.deriBigMatB.fun(), the usage of
  ##        which is not very intuitive. Basically it refills a structure
  ##        that is the same as BB.lst, and use a similar way to construct
  ##        the derivative matrices for every pi. This function can be
  ##        better written later. 
  
  #dB.all: is returned as functions of CC, such as 
  ##    dB.all = list()
  ##    for(i in 1:length(CC)) {
  ##      dB.all[[i]] = construct.dB.list.fun(CC[[i]], subPI.lst)
  ##    }

  #subPI.lst: list of PIs with frequencies and names listed, the
  ##            names have to agree with CC!!!!! 
  #BB.lst: BB.list 

  full.PI = unlist(lapply(subPI.lst, names))
  reduced.PI = names(subPI.lst[[1]])
  l = length(full.PI)
  m = length(reduced.PI)
  bigB.dPI = list()
  for(i in 1:l) {
    bigB.dPI[[i]] = matrix(0, l, m)
    rownames(bigB.dPI[[i]]) = full.PI
    colnames(bigB.dPI[[i]]) = reduced.PI
    attr(bigB.dPI[[i]], "name") = full.PI[i]
  }

  len.subPI = length(subPI.lst)
  index = index.k.fun(len.subPI)  

  for(i in 1:l)
    {
      ##let's do full.PI[i] one at a time;
      current.PI = full.PI[i];
      ##find out which subPI that current.PI belongs to;
      for(j in 1:len.subPI)
        {
          if(!is.na(match(current.PI, names(subPI.lst[[j]]))))
            {
              k = j;
              break;
            }
        }
      ## current.PI belongs to subPI.lst[[k]];
      
      ## make a copy of BB, and change components according to build
      ##   bigB.dPI[[i]]
      BB.deri = BB.lst
      
      ## index associated with k: 
      index.k = index[[k]]
      if(!is.na(sum(index.k)))
        {
          for (jj in index.k)
            {
              BB.deri[[jj]] = dB.all[[jj]][[i]]
            }
          ## Since BB.deri has the same structure as the BB.lst, we are
          ##  using a modified old function to find the derivative matrix. 
          bigB.dPI[[i]] <-get.deriBigMatB.fun(BB.deri, subPI.lst, k)
          attr(bigB.dPI[[i]], "name") = full.PI[i]
        }
    }
  return(bigB.dPI)
}




get.allpossROUTE.fun<-function()
{
  ## when this file is modified, function "getPIstar.fun" need to be changed as well!!!!!!
  m1 <-list(1)
  m2 <-  lapply(m1, 2, FUN = c)
  m3 <-c(lapply(m2, 3, FUN = c),
         lapply(m1, 4, FUN = c))
  m4 <-c(lapply(m3, 5, FUN = c),
         lapply(m2, 6, FUN = c),
         lapply(m1, 7, FUN = c))
  m5 <-c(lapply(m4, 8, FUN = c),
         lapply(m3, 9, FUN = c),
         lapply(m2, 10, FUN = c),
         lapply(m1, 11, FUN = c))
  m6 <-c(lapply(m5, 12, FUN = c),
         lapply(m4, 13, FUN = c),
         lapply(m3, 14, FUN = c),
         lapply(m2, 15, FUN = c),
         lapply(m1, 16, FUN = c))
  m7 <-c(lapply(m6, 17, FUN = c),
         lapply(m5, 18, FUN = c),
         lapply(m4, 19, FUN = c),
         lapply(m3, 20, FUN = c),
         lapply(m2, 21, FUN = c),
         lapply(m1, 22, FUN = c))
  m8 <-c(lapply(m7, 23, FUN = c),
         lapply(m6, 24, FUN = c),
         lapply(m5, 25, FUN = c),
         lapply(m4, 26, FUN = c),
         lapply(m3, 27, FUN = c),
         lapply(m2, 28, FUN = c),
         lapply(m1, 29, FUN = c))
  m9 <-c(lapply(m8, 30, FUN = c),
         lapply(m7, 31, FUN = c),
         lapply(m6, 32, FUN = c),
         lapply(m5, 33, FUN = c),
         lapply(m4, 34, FUN = c),
         lapply(m3, 35, FUN = c),
         lapply(m2, 36, FUN = c),
         lapply(m1, 37, FUN = c))
  m10<-c(lapply(m9, 38, FUN = c),
         lapply(m8, 39, FUN = c),
         lapply(m7, 30, FUN = c),
         lapply(m6, 41, FUN = c),
         lapply(m5, 42, FUN = c),
         lapply(m4, 43, FUN = c),
         lapply(m3, 44, FUN = c),
         lapply(m2, 45, FUN = c),
         lapply(m1, 46, FUN = c))
  
  m = c(m1, m2, m3, m4, m5, m6, m7, m8, m9, m10)
  d = 10
  r = c(1, rep(2, d-1))^(c(0,0:(d-2)))
  di = d * (d-1) / 2 + 1
  
  l = rep(1:d, r)
  mat = matrix(0, sum(r), di)
  for(i in 1:sum(r)) {
    mat[i,m[[i]]] = m[[i]]
  }
  row.names(mat)<-paste("Step ",
                        c(1,2,rep(3,2), rep(4,4), rep(5,8), rep(6,16), rep(7,32), rep(8,64), rep(9,128), rep(10, 256)),
                        ".",
                        c(1,1,   (1:2),    (1:4),    (1:8),    (1:16),    (1:32),    (1:64),    (1:128), (1:256)),   sep="")
  dimnames(mat)[[2]]<-paste("mat",c(1,2,rep(3,2), rep(4,3), rep(5,4), rep(6,5),rep(7,6)),".(",(1:di),")",sep="")
  return(mat)
}


    
#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/tzeng-main.functions.forRDscoreTest.public.r"
## modified on Sep 22, 2006. changes in "getcut.fun"
## add "base=2 for calculating "info"

############################################
## FUNCTIONS needed for Hap-RD
############################################

getcut.fun<-function(pp.org,nn=2*nhap, plot=0){
      pp     <-rev(sort(pp.org))
      ct     <-round(pp*nn)
      dimen  <-log((1:length(pp)), base=2 )/ct
      info   <-cumsum(pp*log(1/pp,base=2))
      netinfo<-info-dimen
      cutpos <-netinfo==max(netinfo)
      if(plot==1){
        plot( netinfo, pp,type="b",cex=0.5)
        barplot(pp.org)
        abline(h=pp[cutpos])
      }
      return( (1:length(pp))[cutpos])
    }

# source("code.getPIstar.getBigMatB")	### Remarked by Wei-Chen Chen 2009-09-16

    final.BigMatB.matC.fun<-function(mmle.cscn, hapreserv, loci, digit){
      ## HERE mmle.cs, mmle.cn MUST NOT CONTAIN "0" CATEGORIES!!!
      ## hapreserv <-sort(unique(c(names(mmle.cs[mmle.cs>=cut.cs]),names(mmle.cn[mmle.cn>=cut.cn]))))
      ini.cscn     <-get.ini.possHap.hapReserv.fun(mmle.cscn, loci, digit, hapreserv)
      mut.mat.cscn <-ini.cscn$prepsi
      subPI.cscn   <-get.subPI.fun(mmle.cscn,  mut.mat=mut.mat.cscn)
      ## get BB matrix
      jnk.cscn     <-get.preBB.skip.fun(subPI.cscn, loci, digit, mmle.cscn)
      preBB.cscn   <-jnk.cscn$preBB.lst
      indexBB.cscn <-jnk.cscn$indexBB
      medBB.cscn   <-get.medBB.skip.fun(subPI.lst=subPI.cscn, preBB.lst=preBB.cscn, indexBB.lst=indexBB.cscn,loci, digit, mmle.cscn)
      fullBB.cscn  <-get.fullBB.fun(medBB.lst=medBB.cscn,  subPI.lst=subPI.cscn)
      BB.cscn      <-get.BB.fun(fullBB.cscn)
      CC.cscn      <-lapply(BB.cscn, function(mat){1*(mat>0)})
      BigMatB.cscn <-getBigMatB.fun(BB.cscn, subPI.cscn)
      return(list(## fullBB.cscn=fullBB.cscn,
                  BB.cscn=BB.cscn,
                  CC.cscn =CC.cscn,
                  BigMatB.cscn =BigMatB.cscn ,
                  subPI.cscn =  subPI.cscn
                  ))
  }


# source("functions.r")	### Remarked by Wei-Chen Chen 2009-09-16
## call  all fun related to calcuating deri.BigB and I13.I33i.I13t
## description is in /home/jytzeng/Research/Hap-RD/Score/MatB.cluster/Fang/generalform.ps



get.matD.matEI.fun<-function(yy=y, mmu=mu, aa=a,
                            XRD.post=xRD.post, XFD.post=x.post, XRD=xRD, XFD = x,
                            post.prob=post, V=v, Nreps=nreps,
                            CC.lst=CC, subPI.lst=subPI, BB.lst=BB,
                            pFD=p.cscn.FD, xx.adj=x.adj
                            ){
  ## this use a hybird I matrix: that is, to use obs I_G for all nonzefo element of E(I_G). So Iaz.2=Iaz.3=0,
  ## and Izz.12= Izz.13=Izz.23 =0 
  ## pFD must match with nrow(matB)
  res    <- ((yy - mmu)/aa)^2
  uRD.mtx<-  (yy - mmu)/aa * XRD.post
  RR     <- ncol(XFD)
  ##-----
  ## Daa
  ##-----
  Daa    <- t(uRD.mtx) %*% uRD.mtx  ## sum_i res * E(Xi*|G) * E(Xi*'|G)
  ## .... double check
  ## tmp.fun<-function(i, mat1=uRD.mtx, mat2=uRD.mtx){mat1[i,]%*%t(mat2[i,])};
  ## .....Daa should = matrix(apply(sapply(1:nhap, tmp.fun, mat1=uRD.mtx, mat2=uRD.mtx), 1,sum), 6,6)
  ##-----
  ## Daz
  ##-----
  Daz.1  <-t(res * XRD.post) %*% xx.adj  ## sum_i res * E(Xi*|G) * t(xx.adj)   ### modified on May 18 for allowing xx.adj
  Daz    <-Daz.1
  ##-----
  ## Dzz 
  ##-----
  Dzz.11 <-t(res * xx.adj) %*% xx.adj  ## sum_i res * xx.adj * t(xx.adj)   ### modified on May 18, 2005 for allowing xx.adj
  Dzz    <-Dzz.11
  ##-----
  ## Iaz
  ##-----
  Iaz.1   <- t(XRD.post) %*% (xx.adj * V) ## it is equal to apply(v * XRD.post,2, sum) ## modified on May 18, 2005 for xx.adj
  Iaz     <- Iaz.1
  ## Izz
  Izz.11 <- t(xx.adj * V) %*% xx.adj    ### modified on May 18 for allowing xx.adj
  Izz    <- Izz.11
  return(list("Daa"=Daa, "Daz"=Daz, "Dzz"=Dzz, "Iaz"=Iaz, "Izz"=Izz))
}





### the function below is for data analysis (allowing for missing data). Will calculate eme and determine cluster base within the fn

haplo.score.RD.unphased.fun<-function(y=Y.vec, geno=Geno, trait.type = "binomial",
                                      miss.val=NA, locus.label=NA, offset = NA, x.adj= NA, skip.haplo = 1e-7){
      ### !!!!! warning: this program can do binary and continuous Y but does not allow ordinal Y yet...
      trait.int <- charmatch(trait.type, c("gaussian", "binomial",  "poisson", "ordinal"))
      if (is.na(trait.int))        stop("Invalid trait type")
      if (trait.int == 0)          stop("Ambiguous trait type")
      if (length(y) != nrow(geno)) stop("Dims of y and geno are not compatible")
      n.loci <- ncol(geno)/2
      if (n.loci != (floor(ncol(geno)/2)))  stop("Odd number of cols of geno")
      adjusted <- TRUE
      if (all(is.na(x.adj)))       adjusted <- FALSE
      if (adjusted)
        {
          x.adj <- as.matrix(x.adj)
          if (nrow(x.adj) != length(y)) stop("Dims of y and x.adj are not compatible")
        }
      miss <- is.na(y)
      if (adjusted)       miss <- miss | apply(is.na(x.adj), 1, any)
      if (trait.int == 3)
        {
          if (all(is.na(offset))) stop("Missing offset")
          miss <- miss | is.na(offset)
        }
      y    <- as.numeric(y[!miss])
      geno <- geno[!miss, ]
      if (adjusted)       x.adj <- x.adj[!miss, , drop = FALSE]
      if (trait.int == 3) offset <- offset[!miss]
      ##-----------------------------------------------------
      ## get EME of hap freq and determine clustering base
      ##-----------------------------------------------------
      haplo <- haplo.em(geno,  locus.label=locus.label, miss.val=miss.val, control = haplo.em.control())
      H.names<-apply(haplo$haplotyp, 1, paste,  collapse="")
      Pcscn  <-haplo$hap.prob; names(Pcscn)<-H.names
      pcscn  <-rev(sort(Pcscn))
#      pos    <-getcut.fun(pcscn,nhap, plot=1)
      pos    <-getcut.fun(pcscn,nhap, plot=0)	### Modified by Wei-Chen Chen 2009-09-16
      reserv <-names(pcscn)[1:pos]
      ## modified on May 18 2005; unmark to allow for missing data

      if (!haplo$converge) stop("EM for haplo failed to converge")
      rows.rem <- haplo$rows.rem
      if (length(rows.rem) > 0)
        {
          keep <- !apply(outer((1:length(y)), rows.rem, "=="), 1, any)
          y <- y[keep]
          if (adjusted)          x.adj  <- x.adj[keep, , drop = FALSE]
          if (trait.int == 3)    offset <- offset[keep]
        }
      if (trait.int == 2)
        {
          if (!all(y == 1 | y == 0))     stop("Invalid y values")
          if (all(y == 1) | all(y == 0)) stop("No variation in y values")
        }

       ##   if (trait.int == 4) {
       ##     y     <- factor(y)
       ##     y.lev <- levels(y)
       ##     y     <- as.numeric(y)
       ##     if (max(y) < 3)             stop("Less than 3 levels for y values")
       ##   }

      ##---- prepare to get X, X*,  E(X*|G), E(X|G) ----
      n.subj     <- length(y)
      hap1       <- haplo$hap1code
      hap2       <- haplo$hap2code
      indx       <- haplo$indx.subj
      post       <- haplo$post
      nreps      <- as.vector(haplo$nreps)
      uhap       <- sort(unique(c(hap1, hap2)))
      which.haplo<- haplo$hap.prob >= skip.haplo
      p.cscn     <- Pcscn[which.haplo]
      h.names    <- names(p.cscn)
      uhap       <- uhap[which.haplo]
      ##---- obtain bigMatB for Reducing Dim ----
      BigMatB    <- final.BigMatB.matC.fun(p.cscn, reserv, n.loci, digit) ##"BB.cscn""CC.cscn""BigMatB.cscn""subPI.cscn"  
      matB       <- BigMatB$BigMatB.cscn
      CC         <- BigMatB$CC.cscn
      BB         <- BigMatB$BB.cscn
      subPI      <- BigMatB$subPI.cscn
      p.cscn.FD  <- p.cscn[match(rownames(matB), h.names)]
      p.cscn.RD  <- p.cscn.FD %*% matB
      RR         <- nrow(matB)
      RRstar     <- ncol(matB)
      ##---- get E(X*|G) =  xRD.post ----
      x          <- 0.5* (outer(hap1, uhap, "==") + outer(hap2, uhap, "==")); colnames(x)<-h.names
      x          <- x[,match(rownames(matB), h.names)]
      ## the "0.5" is because in Schaid's coding the rowsum of X vec is 2 while mine is 1
      ## xRD        <- x %*% matB[,-1]  ## modified on May 17, 2005
      xRD        <- x %*% matB
      n.xRD      <- ncol(xRD)
      xRD.post   <- matrix(rep(NA, n.subj * n.xRD), ncol = n.xRD);for(j in 1:n.xRD){ xRD.post[,j]<-tapply(xRD[,j]*post,indx,sum)}
      ##---- get E(X|G) =  x.post ----
      n.x        <- ncol(x)
      x.post     <- matrix(rep(NA, n.subj * n.x), ncol = n.x);for(j in 1:n.x){ x.post[,j]<-tapply(x[,j]*post,indx,sum)}; 
      ##-----------------------------------
      ## Note E(X*|G) should = E(X|G) %*% matB[,-1], i.e., x.post %*% matB[,-1] - xRD.post should be a 0 mat
      ##-----------------------------------
      ##     get score and Var(score)     
      ##-----------------------------------
    if (trait.int == 4) {
      stop("sorry, the program does not allow ordinal Y yet...")
       ## if (adjusted) {
       ##     library("Design")
       ##     library("Hmisc")
       ##     reg.out <- lrm(y ~ x.adj)
       ##     K <- max(y)
       ##     n.xadj <- ncol(x.adj)
       ##     alpha <- reg.out$coef[1:(K - 1)]
       ##     beta <- reg.out$coeff[K:(K - 1 + n.xadj)]
       ##     tmp <- haplo.score.podds(y, alpha, beta, x.adj, nreps, 
       ##         x.post, post, x)
       ## }
       ## if (!adjusted) {
       ##     tbl <- table(y)
       ##     s <- 1 - (cumsum(tbl) - tbl)/n.subj
       ##     alpha <- -log((1 - s[-1])/s[-1])
       ##     tmp <- haplo.score.podds(y, alpha, beta = NA, x.adj = NA, 
       ##         nreps, x.post, post, x)
       ## }
       ## u.score <- tmp$u.score
       ## v.score <- tmp$v.score
    } else  if (trait.int <= 3) {
      if (!adjusted)
        {
          mu    <- switch(trait.int, mean(y), mean(y), sum(y)/sum(offset))
          a     <- switch(trait.int, var(y), 1, 1)
          x.adj <- matrix(rep(1, n.subj), ncol = 1)
        }
      if (adjusted)
        {
          reg.out <- glm(y ~ x.adj, family = trait.type)
          x.adj   <- cbind(rep(1, n.subj), x.adj)
          mu      <- reg.out$fitted.values
          a       <- switch(trait.int, sum(reg.out$residuals^2)/reg.out$df.residual, 1, 1)
        }
      v           <- switch(trait.int, 1/a, mu * (1 - mu), mu)   ## schaid's "v"=b"/a instead of b" * a
      ##----- get score -------------       
      u.mtx     <- (y - mu) * xRD.post/a
      u.score   <- apply(u.mtx, 2, sum)  ##<-- by sum_i (yi-mu)/a * E(Xi*|G)
      ## --------------------------
      ## the result here should match with using  sum_i (yi-mu)/a * E(Xi|G) %*% matB[,-1]
      ## --------------------------
      ##  uFD.mtx  <-(y-mu)/a * x.post; uFD.score<-apply(uFD.mtx, 2, sum ;   uRD.score<- uFD.score %*% matB[,-1] 
      ## --------------------------
      ##---- get variance ----
      ## Info.mat<-get.matD.matI.fun(yy=y, mmu=mu, aa=a, XRD.post=xRD.post, XFD.post=x.post, XRD=xRD, XFD = x,
      ##                             post.prob=post, V=v, Nreps=nreps, CC.lst=CC, subPI.lst=subPI, BB.lst=BB, pFD=p.cscn.FD, xx.adj=x.adj)
      ## modified on May 19 2005. Use hybrid matEI instead of matI to allow continuous and binary Y:
      Info.mat<-get.matD.matEI.fun(yy=y, mmu=mu, aa=a, XRD.post=xRD.post, XFD.post=x.post, XRD=xRD, XFD = x,
                                  post.prob=post, V=v, Nreps=nreps, CC.lst=CC, subPI.lst=subPI, BB.lst=BB, pFD=p.cscn.FD, xx.adj=x.adj)
      ##  "Daa" "Daz" "Dzz" "Iaz" "Izz"
      if(length(Info.mat$Izz)==1){ I22i = 1/Info.mat$Izz }else{      I22i   <-Ginv(Info.mat$Izz)$Ginv }
      v.score<-(Info.mat$Daa - Info.mat$Iaz %*% I22i %*% t(Info.mat$Daz) - Info.mat$Daz %*% I22i %*% t(Info.mat$Iaz)
                + Info.mat$Iaz %*% I22i %*% Info.mat$Dzz %*% I22i %*% t(Info.mat$Iaz))
    }
      tmp   <- Ginv(v.score)
      df    <- tmp$rank
      g.inv <- tmp$Ginv
      score.global <- u.score %*% g.inv %*% u.score
      score.haplo  <- u.score/sqrt(diag(v.score))
      score.max    <- max(score.haplo^2, na.rm = TRUE)
      
      score.global.p <- 1 - pchisq(score.global, df)
      score.haplo.p  <- 1 - pchisq(score.haplo^2, 1)
      if (all(is.na(locus.label))) {   locus.label <- paste("loc-", 1:n.loci, sep = "")      }
      obj <- (list(score.global = score.global,
                   df = df,
                   score.global.p = score.global.p, 
                   score.haplo = score.haplo, 
                   score.haplo.p = score.haplo.p,
                   ## haplotype = colnames(matB)[-1],  ## modified on May 17, 2005
                   haplotype = colnames(matB),
                   matB = matB,
                   hap.prob.RD = p.cscn.RD,
                   hap.prob.FD = p.cscn.FD,
                   locus.label = locus.label
                   ))
    }




#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/tzeng-modified.r"
### This file contains codes modified from
### "tzeng-main.function.forRDscoreTest.public.r".
### Reference: J.-Y. Tzeng (2005) Genetic Epidemiology.
###            J.-Y. Tzeng, C.-H. Wang, J.-T. K, and C.-H. K. Hsiao (2006) Am. J. Hum. Genet.
### Modified: Wei-Chen Chen 2009-09-17

### X: array[N, L], N haplotypes, L loci.
### ploidy: 1 or 2 => index of X = rep(1:(N/ploidy), each = ploidy)

### This function returns a fake object "haplo" to replace the output of haplo.em().
haplo.phase.known <- function(X, ploidy = 2){
  X.names <- apply(X, 1, paste, collapse="")
  unique.X <- unique(X)
  unique.X.names <- apply(unique.X, 1, paste, collapse="")
  X.summary <- as.data.frame(table(X.names))
  X.summary$prob <- X.summary$Freq / sum(X.summary$Freq)

  haplo <- NULL
  haplo$haplotype <- unique.X
  haplo$hap.prob <- X.summary$prob[match(unique.X.names, as.character(X.summary$X.names))]
  haplo$post <- rep(1, nrow(X))
  if(ploidy == 1){
    haplo$hap1code <- match(X.names, unique.X.names)
    haplo$hap2code <- haplo$hap1code
    haplo$indx.subj <- 1:nrow(X)
  } else if(ploidy == 2){
    tmp <- matrix(match(X.names, unique.X.names), nrow = 2)
    haplo$hap1code <- tmp[1,]
    haplo$hap2code <- tmp[2,]
    haplo$indx.subj <- rep(1:(nrow(X)/2), each = 2)
  } else{
    stop("ploidy = 1 or 2")
  }
  haplo
} # End of haplo.phase.know().



### This function returns posterior probabilities based on the fake "haplo" object.
haplo.post.prob <- function(X, ploidy = 2, skip.haplo = 1e-7, K = NULL){
### Add pos for arbitrarily choosing K. Use default getbut.fun() if K = NULL.
### Modified: Wei-Chen Chen 2009-12-05.
  haplo <- haplo.phase.known(X, ploidy)	# The fake "haplo" object.
  nhap <- nrow(X)
  n.loci <- ncol(X)
  n.subj <- nrow(X)
  digit <- 1

      H.names<-apply(haplo$haplotype, 1, paste,  collapse="")
      Pcscn  <-haplo$hap.prob; names(Pcscn)<-H.names
      pcscn  <-rev(sort(Pcscn))

      if(is.null(K)){
        pos    <-getcut.fun(pcscn,nhap, plot=0)
      } else{
        tl.K <- length(pcscn)
        if(K < 1 || K > (tl.K - 1)){
          stop(paste("K should be in [", 1, ", ", tl.K - 1, "]", sep = ""))
        }
        pos <- K
      }
      reserv <-names(pcscn)[1:pos]

      ##---- prepare to get X, X*,  E(X*|G), E(X|G) ----
#      n.subj     <- length(y)
      hap1       <- haplo$hap1code
      hap2       <- haplo$hap2code
      indx       <- haplo$indx.subj
      post       <- haplo$post
      uhap       <- sort(unique(c(hap1, hap2)))
      which.haplo<- haplo$hap.prob >= skip.haplo
      p.cscn     <- Pcscn[which.haplo]
      h.names    <- names(p.cscn)
      uhap       <- uhap[which.haplo]
      ##---- obtain bigMatB for Reducing Dim ----
      BigMatB    <- final.BigMatB.matC.fun(p.cscn, reserv, n.loci, digit) ##"BB.cscn""CC.cscn""BigMatB.cscn""subPI.cscn"  
      matB       <- BigMatB$BigMatB.cscn
#      CC         <- BigMatB$CC.cscn
#      BB         <- BigMatB$BB.cscn
#      subPI      <- BigMatB$subPI.cscn
#      p.cscn.FD  <- p.cscn[match(rownames(matB), h.names)]
#      p.cscn.RD  <- p.cscn.FD %*% matB
#      RR         <- nrow(matB)
#      RRstar     <- ncol(matB)
      ##---- get E(X*|G) =  xRD.post ----
      x          <- 0.5* (outer(hap1, uhap, "==") + outer(hap2, uhap, "==")); colnames(x)<-h.names
      x          <- x[,match(rownames(matB), h.names)]
      ## the "0.5" is because in Schaid's coding the rowsum of X vec is 2 while mine is 1
      ## xRD        <- x %*% matB[,-1]  ## modified on May 17, 2005
      xRD        <- x %*% matB
      n.xRD      <- ncol(xRD)
      xRD.post   <- matrix(rep(NA, n.subj * n.xRD), ncol = n.xRD);for(j in 1:n.xRD){ xRD.post[,j]<-tapply(xRD[,j]*post,indx,sum)}
      ##---- get E(X|G) =  x.post ----
      n.x        <- ncol(x)
      x.post     <- matrix(rep(NA, n.subj * n.x), ncol = n.x);for(j in 1:n.x){ x.post[,j]<-tapply(x[,j]*post,indx,sum)}; 

  list(haplo = haplo,
       FD.id = match(rownames(matB), H.names),
       RD.id = match(colnames(matB), H.names),
       FD.post = x.post,
       RD.post = xRD.post,
       g.truncate = ncol(matB))
} # End of haplo.post.prob()

#line 1 "/private/tmp/Rtmpvlsw7d/R.INSTALL10d63af1/phyclust/R/zzz.r"
.First.lib <- function(lib, pkg)
{
  library(ape)
  library.dynam("phyclust", pkg, lib)
} # End of .First.lib()

.Last.lib <- function(libpath)
{
#  ret <- .dynLibs()
#  for(dynlib in ret){
#    if(dynlib[1] == "phyclust"){
#      phyclust.libpath <- gsub("/libs/phyclust\\..*", "", dynlib[2])
#      library.dynam.unload("phyclust", phyclust.libpath)
#    }
#  }

  library.dynam.unload("phyclust", libpath)
} # End of .Last.lib()
