.packageName <- "ecolitk"
.buildBNUM2GENBANK <- function() {
  return(.buildGENEPRODFUN(what="genbank"))
}

.buildBNUM2GENEPRODUCT <- function() {
  return(.buildGENEPRODFUN(what="gene product"))
}

.buildGENEPRODFUN <- function(my.file=NULL, my.url="http://genprotec.mbl.edu/files/geneproductfunctions.txt", what="") {
  choices <- c("bnum", "genbank", "gene", "gene type", "gene product")
  what <- match.arg(what, choices)
  what.i <- match(what, choices)
  env.x2y <- new.env(hash=TRUE)
  ##env.genbank <- new.env(hash=TRUE)
  if (is.null(my.file))
    con <- url(my.url, open="r")
  else
    con <- file(my.file, open="r")

  ## skip comments
  mycomments <- ""
  line <- readLines(con, n=1)
  while (length(grep("^bnum", line)) == 0) {
    line <- readLines(con, n=1)
    mycomments <- paste(mycomments, line, sep="\n")
  }
  mycomments <- paste(mycomments, "\n", sep="")
  line <- readLines(con, n=1)

  ##
  while (length(line) != 0) {
    if (line != "") {
      m <- strsplit(line, "\t", extended=TRUE)[[1]]
      bnum <- m[1]
      genbank <- m[what.i]
      if (genbank == "0")
        genbank <- NA
      if (exists(bnum, envir=env.x2y))
        assign(bnum, unique(c(genbank, get(bnum, envir=env.x2y))), envir=env.x2y)
      else
        assign(bnum, genbank, envir=env.x2y)
    }
    line <- readLines(con, n=1)
  }

  attr(env.x2y, "comments") <- mycomments
  return(env.x2y)

}


.buildMultiFun <- function(my.file=NULL, my.url="http://genprotec.mbl.edu/files/MultiFun.txt") {
  env.multiFun <- new.env(hash=TRUE)

  if (is.null(my.file))
    con <- url(my.url, open="r")
  else
    con <- file(my.file, open="r")

  ## skip comments
  mycomments <- ""
  line <- readLines(con, n=1)
  while (length(grep("^[0-9]", line)) == 0) {
    line <- readLines(con, n=1)
    mycomments <- paste(mycomments, line, sep="\n")
  }
  mycomments <- paste(mycomments, "\n", sep="")

  ##
  while (length(line) != 0) {
    if (line != "") {
      ## trim leading white spaces
      m <- regexpr("(\\w\\.?)+", line, perl=TRUE)
      multiFun <- substr(line, m, m+attr(m, "match.length")-1)
      multiFun <- sub("\\.$", "", multiFun)
      multiFunAnnotation <- substr(line, m + attr(m, "match.length") + 1, nchar(line))
      assign(multiFun, multiFunAnnotation, envir=env.multiFun)
    }
    line <- readLines(con, n=1)
  }
  attr(env.multiFun, "comments") <- mycomments
  return(env.multiFun)
}

.buildBNUM2MULTIFUN <- function(my.file=NULL, my.url="http://genprotec.mbl.edu/files/multifunassignments.txt") {
  env.bnum2multifun <- new.env(hash=TRUE)

  if (is.null(my.file))
    con <- url(my.url, open="r")
  else
    con <- file(my.file, open="r")

  ## skip comments
  mycomments <- ""
  line <- readLines(con, n=1)
  while (length(grep("^Bnum", line)) == 0) {
    line <- readLines(con, n=1)
    mycomments <- paste(mycomments, line, sep="\n")
  }
  mycomments <- paste(mycomments, "\n", sep="")
  line <- readLines(con, n=1)

  ##
  while (length(line) != 0) {
    if (line != "") {
      ## trim leading white spaces
      m <- regexpr("b\\d+", line, perl=TRUE)
      bnum <- substr(line, m, m+attr(m, "match.length")-1)
      m <- regexpr("\\s\\S+", line, perl=TRUE)
      multiFun <- substr(line, m+1, m+attr(m, "match.length")-1)
      if (exists(bnum, envir=env.bnum2multifun))
        assign(bnum, unique(c(multiFun, get(bnum, envir=env.bnum2multifun))), envir=env.bnum2multifun)
      else
        assign(bnum, multiFun, envir=env.bnum2multifun)
    }
    line <- readLines(con, n=1)
  }

  return(env.bnum2multifun)
}

.buildBNUM2SYMBOL <- function(my.url="http://genprotec.mbl.edu/files/MultiFun.txt") {

}

.buildMULTIFUN2GO <- function(filename) {
  env.multiFun2GO <- new.env(hash=TRUE)
  con <- file(filename, open="r")

  ## skip comments
  mycomments <- ""
  line <- readLines(con, n=1)
  while (length(grep("^!", line)) > 0) {
      line <- readLines(con, n=1)
      mycomments <- paste(mycomments, line, sep="\n")
    }
  mycomments <- paste(mycomments, "\n", sep="")
  ##

  while (length(line) != 0) {
    ## get MultiFun number
    m <- regexpr("^MultiFun:(\\.+?) ", line, perl=TRUE)
    multiFun <- substr(line, 10, attr(m, "match.length") - 1)
    ## get the GOs
    m <- regexpr("> \\.+$", line, perl=TRUE)
    GOs <- substr(line, m+2, m+attr(m, "match.length"))
    if (GOs == "GO:.") {
      ## no go ;)
      GOs <- NA
    } else {
      GOs <- strsplit(GOs, " > ")[[1]]
      GOs <- strsplit(GOs, " ; ")
      GOs <- unlist(lapply(GOs, function(x) x[2]))
    }
    assign(multiFun, GOs, envir=env.multiFun2GO)

    line <- readLines(con, n=1)
  }
  attr(env.multiFun2GO, "comments") <- mycomments
  return(env.multiFun2GO)
}


.buildMultiFunGraphNEL <- function(filename) {
  require("graph") || stop("The graph packages is needed for this operation")

  nid.i <- 1
  nname.i <- 2

  linesInFile <- readLines(filename)

  tmp <- strsplit(linesInFile, ";")

  r <- lapply(tmp, function(x, y) strsplit(x[nid.i], y), "\\.")
  r.names <- unlist(lapply(tmp, function(x, y) x[nname.i]))

  mFunNodenames <- unique(unlist(lapply(tmp, function(x) x[nname.i])))
  nodename2i <- new.env(hash=TRUE)
  nodename2nodeid <- new.env(hash=TRUE)
  nodeid2nodename <- new.env(hash=TRUE)

  multiassign(mFunNodenames, seq(along=mFunNodenames), nodename2i)
  mFunEdges <- vector("list", length=length(mFunNodenames))
  names(mFunEdges) <- mFunNodenames
  multiassign(unlist(lapply(tmp, function(x) x[nid.i])), unlist(lapply(tmp, function(x) x[nname.i])),
              envir=nodeid2nodename)

  for (i in seq(along=r)) {
    n <- length(r[[i]][[1]])
    if (n == 1)
      next
    parent <- paste(r[[i]][[1]][seq(1, n-1, length=n-1)], collapse=".")
    if (! exists(parent, nodeid2nodename))
      next
    parent.i <- get(get(parent, nodeid2nodename), nodename2i)
    child.i <- get(tmp[[i]][nname.i], nodename2i)
    ##parent.i <- get(get(parent, nodeid2nodename), nodename2i)
    ##child.i <- i
    if (is.null(mFunEdges[[parent.i]])) {
      mFunEdges[[parent.i]]$edges = child.i
      mFunEdges[[parent.i]]$weights = 1
    } else {
      mFunEdges[[parent.i]]$edges = c(mFunEdges[[parent.i]]$edges, child.i)
      mFunEdges[[parent.i]]$weights = c(mFunEdges[[parent.i]]$weights, 1)
    }

    if (is.null(mFunEdges[[child.i]])) {
      ##meshedges[[child.i]]$edges = parent.i
###meshedges[[parent.i]]$edges = paste(r[[i]][[1]], collapse=".")
                                        #meshedges[[child.i]]$weights = 1
    } else {
###meshedges[[parent.i]]$edges = c(meshedges[[parent.i]]$edges, paste(r[[i]][[1]], collapse="."))
      ##meshedges[[child.i]]$edges = c(meshedges[[child.i]]$edges, parent.i)
      ##meshedges[[child.i]]$weights = c(meshedges[[child.i]]$weights, 1)
    }
  }

  ##meshnodenames <- seq(along=meshnodenames)
  ##names(meshedges) <- meshnodenames

  gmesh <- new("graphNEL", nodes=mFunNodenames, edgeL=mFunEdges, edgemode="directed")

  return(gmesh)
}


linkedmultiget <- function(x, envir.list=list(), unique=TRUE) {

  if (! is.character(x))
    stop("x must be a vector of mode 'character'")

  f <- function(x, y) {
    tmp <- mget(x, envir=y, ifnotfound=NA)
    tmp <- unlist(tmp)
    if (all(is.na(tmp)))
      tmp <- as.character(tmp)
    if (! is.character(tmp))
      stop("Values in environments must be of mode 'character'")
    return(tmp)
  }

  ##r <- vector("list", length=length(x))
  r <- as.list(x)
  for (i in seq(along=envir.list)) {
    r <- lapply(r, f, envir.list[[i]])
    if (unique)
      r <- lapply(r, unique)
  }
  names(r) <- x
  return(r)
}
cPlotCircle <- function(radius=1, xlim=c(-2, 2), ylim=xlim, edges=300, main=NULL, main.inside=NULL, ...) {
  plot.new()
  plot.window(xlim, ylim, ...)
  linesCircle(radius, edges=edges, ...)
  title(main)
  text(0, 0, main.inside)
}

linesCircle <- function(radius, center.x = 0, center.y = 0, edges=300, ...) {
  xy <- polar2xy(radius, seq(0, 2*pi, length=edges))
  ##x <- (radius * cos(seq(0, 2*pi, length=edges))) + center.x
  ##y <- (radius * sin(seq(0, 2*pi, length=edges))) + center.y
  lines(xy$x + center.x, xy$y + center.y, ...)
}

chromPos2angle <- function(pos, len.chrom, rot=pi/2, clockwise=TRUE) {
  if (any(abs(pos) > len.chrom, na.rm=TRUE))
    warning(paste(pos, ">", len.chrom, ": abs(pos) > len.chrom !!!"))
  
  theta <- pos * 2 * pi / len.chrom
  
  if (clockwise)
    theta <- - theta

  theta <- theta + rot

  return(theta)
}

polygonDisk <- function(radius, center.x=0, center.y=0, edges=300, ...) {
  ##x <- (radius * cos(seq(0, 2*pi, length=edges))) + center.x
  ##y <- (radius * sin(seq(0, 2*pi, length=edges))) + center.y
  
  xy <- polar2xy(radius, seq(0, 2*pi, length=edges))
  polygon(xy$x + center.x, xy$y + center.y, ...)
  
}

# linesPolar <- function(theta, radius, center.x = 0, center.y = 0, ...) {
#   xy <- polar2xy(radius, theta0, theta1, length=edges))
#   xy$x <- xy$x + center.x
#   xy$y <- xy$y + center.y
#   lines(xy$x, xy$y, ...)
# }

pointsArc <- function(theta0, theta1, radius, center.x = 0, center.y = 0, ...) {
  xy <- polar2xy(radius, seq(theta0, theta1, length=length(radius)))
  ##x <- (radius * cos(seq(0, 2*pi, length=edges))) + center.x
  ##y <- (radius * sin(seq(0, 2*pi, length=edges))) + center.y
  points(xy$x + center.x, xy$y + center.y, ...)
}

linesArc <- function(theta0, theta1, radius, center.x = 0, center.y = 0, ...) {
  xy <- polar2xy(radius, seq(theta0, theta1, length=length(radius)))
  ##x <- (radius * cos(seq(0, 2*pi, length=edges))) + center.x
  ##y <- (radius * sin(seq(0, 2*pi, length=edges))) + center.y
  lines(xy$x + center.x, xy$y + center.y, ...)
}

arrowsArc <- function(theta0, theta1, radius, center.x = 0, center.y = 0, edges = 10,
                      length = 0.25, angle = 30, code = 2, ...) {
  xy <- polar2xy(radius, seq(theta0, theta1, length=edges))
  xy$x <- xy$x + center.x
  xy$y <- xy$y + center.y
  lines(xy$x, xy$y, ...)
  n <- length(xy$x)
  if (code == 2 | code == 3)
    arrows(xy$x[n-1], xy$y[n-1], xy$x[n], xy$y[n], length=length, angle=angle, ...)
  if (code == 1 | code == 3)
    arrows(xy$x[1], xy$y[1], xy$x[2], xy$y[2], length=length, angle=angle, ...)
  
}

polygonArc <- function(theta0, theta1, radius.in, radius.out,
                       center.x = 0, center.y = 0,
                       edges=10,
                       col="black",
                       border = NA,
                       ...) {
  
  
  if (length(edges) == 1)
    edges <- rep(edges, length=length(theta0))

  col <- rep(col, length = length(theta0))

  ok <- ! (is.na(theta0) | is.na(theta1))
  
  for (i in seq(along=theta0[ok])) {
    theta.seq <- seq(theta0[ok][i], theta1[ok][i], length=edges[ok][i])
    x <- c(radius.in * cos(theta.seq), radius.out * cos(rev(theta.seq))) + center.x 
    y <- c(radius.in * sin(theta.seq), radius.out * sin(rev(theta.seq))) + center.y
    polygon(x, y, col=col[ok][i], border=border, ...)
  }
}

polygonChrom <- function(begin, end, len.chrom,
                         radius.in, radius.out,
                         total.edges=300,
                         edges=max(round(abs(end-begin)/len.chrom* total.edges), 2, na.rm=TRUE),
                         rot=pi/2, clockwise=TRUE,
                         ...) {

  theta0 <- chromPos2angle(begin, len.chrom, rot=rot, clockwise=clockwise)
  theta1 <- chromPos2angle(end, len.chrom, rot=rot, clockwise=clockwise)

  if (any(theta0 == theta1, na.rm=TRUE))
    warning(paste("identical angles for: ", which(theta0 == theta1), collapse=TRUE))
  
  polygonArc(theta0, theta1, radius.in, radius.out, edges=edges, ...)
}

linesChrom <- function(begin, end, len.chrom, radius,
                       total.edges=300,
                       edges=max(round(abs(end-begin)/len.chrom* total.edges), 2, na.rm=TRUE),
                       rot=pi/2, clockwise=TRUE,
                       ...) {

  theta0 <- chromPos2angle(begin, len.chrom, rot=rot, clockwise=clockwise)
  theta1 <- chromPos2angle(end, len.chrom, rot=rot, clockwise=clockwise)

  if (any(theta0 == theta1, na.rm=TRUE))
    warning(paste("identical angles for: ", which(theta0 == theta1), collapse=TRUE))

  if (length(edges) == 1)
    edges <- rep(edges, length=length(theta0))

  ok <- ! (is.na(theta0) | is.na(theta1))
  
  for (i in seq(along=theta0[ok])) {
    linesArc(theta0[ok][i], theta1[ok][i], rep(radius, edges[ok][i]), ...)
  }
}

ecoli.len <- 4639221
xy2polar <- function(x, y) {
  if (missing(y)) {
    y <- x$y
    x <- x$x
  }
  rho <- sqrt(x^2 + y^2)
  theta <- atan(y/x)
  return(list(rho=rho, theta=theta))
}

polar2xy <- function(rho, theta) {
  if (missing(theta)) {
    theta <- rho$theta
    rho <- rho$rho
  }
  x <- rho * cos(theta)
  y <- rho * sin(theta)
  return(list(x=x, y=y))
}

rotate <- function(x, y, alpha) {
  pol <- xy2polar(x, y)
  pol$theta <- pol$theta + alpha
  xy <- polar2xy(pol)
  return(xy)
}
                    
wstringapply <- function(x, SIZE, SLIDE, FUN, ...) {
  n <- (nchar(x) - SIZE) %/% SLIDE
  res <- vector("list", length = n)
  for (i in seq(0, n-1, length=n)) {
    ##x.start <- (i-1) * SIZE + 1
    ##x.stop <- min(nchar(x), i * SIZE)
    w.x <- substring(x, i * SLIDE + 1, i * SLIDE + SIZE)
    res[[i+1]] <- FUN(w.x, ...)
    ##res[[i]] <- FUN(w.x)
  }
  return(res)
}

# wapply <- function(x, SIZE, FUN, ...) {
#   n <- length(x)
#   res <- vector("list", length = n-1)
#   for (i in seq(1, n-1, length=n-1)) {
#     w.x <- x[i, i+SIZE]
#     res[[i]] <- FUN(w.x, ...)
#     ##res[[i]] <- FUN(w.x)
#   }
#   return(res)
# }

gccontent <- function(x) {x <- toupper(x); n <- nchar(x); sum(as.integer(strsplit(x, "")[[1]] %in% c("G", "C")))/n}
