### -----------------------------------------------------------------
### get the subset of axt based on chr, start, end and from which sequence
### Exported!!

setMethod("subAxt", signature(x="Axt", chr="character",
                              start="numeric", end="numeric"),
          function(x, chr, start, end, select=c("target", "query"), qSize=NULL){
            .subAxtMultiple(x, chr=chr, start=start, end=end,
                            select=select, qSize=qSize)
          }
          )

setMethod("subAxt", signature(x="Axt", chr="character",
                              start="missing", end="missing"),
          function(x, chr, start, end, select=c("target", "query"), qSize=NULL){
            select <- match.arg(select)
            if(select == "target"){
              ans <- x[seqnames(targetRanges(x)) == chr]
            }else{
              ans <- x[seqnames(queryRanges(x)) == chr]
            }
            return(ans)
          }
          )

## This is to fetch the Axts within the specific chrs, starts, ends
## based on target sequences.
.subAxtFull <- function(x, chr, start, end,
                        select=c("target", "query"),
                        #type=c("any", "within"),
                        qSize=NULL){
  #type <- match.arg(type)
  select <- match.arg(select)
  if(select == "query"){
    if(is.null(qSize))
      stop("When selecting on query alignments,
           qSize must be provided.")
           if(!is(qSize, "integer"))
             stop("qSize must be an integer object.")
  }
  if(length(chr) != 1L || length(start) != 1L || length(end) != 1L){
    stop("Currently we only support chr, start, end in length of 1.")
  }
  start = as.integer(start)
  end = as.integer(end)

  if(select == "target"){
    searchGRanges <- GRanges(seqnames=chr,
                             ranges=IRanges(start=start, end=end),
                             strand="+")
    ## First we find the axts totally within the coordinates
    indexWithin <- queryHits(findOverlaps(targetRanges(x),
                                          searchGRanges, type="within",
                                          select="all"))
    indexAny <- queryHits(findOverlaps(targetRanges(x),
                                       searchGRanges, type="any",
                                       select="all"))
    ## Specially take care of the axt just partially overlapped with the range
    indexPartial <- setdiff(indexAny, indexWithin)
    newAxts = list()
    if(length(indexPartial) > 0L){
      newStarts <- pmax(start(targetRanges(x)[indexPartial]), start)
      newEnds <- pmin(end(targetRanges(x)[indexPartial]), end)
      for(i in 1:length(indexPartial)){
        newAxts[[i]] <- subAln(x[indexPartial[i]], newStarts[i], newEnds[i],
                               select="target")
      }
    }
    ans <- c(x[indexWithin], do.call(c, newAxts))
    return(ans)
  }else if(select == "query"){
    ## first search Axts on positive strand
    searchGRanges <- GRanges(seqnames=chr,
                             ranges=IRanges(start=start, end=end),
                             strand="+")
    indexPositiveWithin <- queryHits(findOverlaps(queryRanges(x),
                                                  searchGRanges, type="within",
                                                  select="all"))
    indexPositiveAny <- queryHits(findOverlaps(queryRanges(x),
                                               searchGRanges, type="any",
                                               select="all"))
    indexPartial <- setdiff(indexPositiveAny, indexPositiveWithin)
    newAxts = list()
    if(length(indexPartial) > 0L){
      newStarts <- pmax(start(queryRanges(x))[indexPartial], start)
      newEnds <- pmin(end(queryRanges(x))[indexPartial], end)
      for(i in 1:length(indexPartial)){
        newAxts[[i]] <- subAln(x[indexPartial[i]], newStarts[i], newEnds[i],
                               select="query")
      }
    }
    # then search Axts on negative strand.
    # we need to prepare the searchGRanges on negative strand.
    searchGRanges <- GRanges(seqnames=chr,
                             ranges=IRanges(start=qSize-end+1,
                                            end=qSize-start+1),
                             strand="-")
    indexNegativeWithin <- queryHits(findOverlaps(queryRanges(x),
                                                  searchGRanges, type="within",
                                                  select="all"))
    indexNegativeAny <- queryHits(findOverlaps(queryRanges(x),
                                               searchGRanges, type="any",
                                               select="all"))
    indexPartial <- setdiff(indexNegativeAny, indexNegativeWithin)
    newAxts2 = list()
    if(length(indexPartial) > 0L){
      newStarts <- pmax(start(queryRanges(x))[indexPartial], qSize-end+1)
      newEnds <- pmin(end(queryRanges(x))[indexPartial], qSize-start+1)
      for(i in 1:length(indexPartial)){
        newAxts2[[i]] <- subAln(x[indexPartial[i]], newStarts[i], newEnds[i],
                                select="query")
      }
    }
    index <- sort(c(indexPositiveWithin, indexNegativeWithin))
    ans <- c(x[index], do.call(c, newAxts), do.call(c, newAxts2))
    return(ans)
  }
}

.subAxtMultiple <- function(x, chr, start, end,
                            select=c("target", "query"),
                            qSize=NULL){
  ## We will allow multiple chr, start, end
  select <- match.arg(select)
  if(select == "query"){
    if(is.null(qSize)){
      stop("When selecting on query alignments,
           qSize must be provided.")
    }
    if(!is(qSize, "integer")){
      stop("qSize must be an integer object.")
    }
  }
  searchGRanges <- GRanges(seqnames=chr,
                           ranges=IRanges(start=start, end=end),
                           strand="+")
  searchGRanges <- reduce(searchGRanges)
  start <- start(searchGRanges)
  end <- end(searchGRanges)
  if(select == "target"){
    ## First we find the axts totally within the coordinates
    hitsWithin <- findOverlaps(targetRanges(x),
                               searchGRanges, type="within",
                               select="all")
    indexWithin <- queryHits(hitsWithin)
    hitsAny <- findOverlaps(targetRanges(x),
                            searchGRanges, type="any",
                            select="all")
    indexAny <- queryHits(hitsAny)
    ## Specially take care of the axt just partially overlapped with the range
    hitsPartial <- setdiff(hitsAny, hitsWithin)
    newAxts = list()
    if(length(hitsPartial) > 0L){
      newStarts <- pmax(start(targetRanges(x)[queryHits(hitsPartial)]), 
                        start(searchGRanges)[subjectHits(hitsPartial)])
      newEnds <- pmin(end(targetRanges(x)[queryHits(hitsPartial)]), 
                      end(searchGRanges)[subjectHits(hitsPartial)])
      for(i in 1:length(hitsPartial)){
        newAxts[[i]] <- subAln(x[queryHits(hitsPartial[i])], 
                               newStarts[i], newEnds[i],
                               select="target")
      }
    }
    ans <- c(x[indexWithin], do.call(c, newAxts))
    return(ans)
  }else{
    ## first search Axts on positive strand
    hitsPositiveWithin <- findOverlaps(queryRanges(x),
                                       searchGRanges, type="within",
                                       select="all")
    indexPositiveWithin <- queryHits(hitsPositiveWithin)
    hitsPositiveAny <- findOverlaps(queryRanges(x),
                                    searchGRanges, type="any",
                                    select="all")
    indexPositiveAny <- queryHits(hitsPositiveAny)
    hitsPartial <- setdiff(hitsPositiveAny, hitsPositiveWithin)
    newAxts = list()
    if(length(hitsPartial) > 0L){
      newStarts <- pmax(start(queryRanges(x))[queryHits(hitsPartial)], 
                        start(searchGRanges)[subjectHits(hitsPartial)])
      newEnds <- pmin(end(queryRanges(x))[queryHits(hitsPartial)], 
                      end(searchGRanges)[subjectHits(hitsPartial)])
      for(i in 1:length(hitsPartial)){
        newAxts[[i]] <- subAln(x[queryHits(hitsPartial[i])], 
                               newStarts[i], newEnds[i],
                               select="query")
      }
    }
    # then search Axts on negative strand.
    # we need to prepare the searchGRanges on negative strand.
    searchGRanges <- GRanges(seqnames=chr,
                             ranges=IRanges(start=qSize-end+1,
                                            end=qSize-start+1),
                             strand="-")
    searchGRanges <- reduce(searchGRanges)
    hitsNegativeWithin <- findOverlaps(queryRanges(x),
                                       searchGRanges, type="within",
                                       select="all")
    indexNegativeWithin <- queryHits(hitsNegativeWithin)
    hitsNegativeAny <- findOverlaps(queryRanges(x),
                                    searchGRanges, type="any",
                                    select="all")
    indexNegativeAny <- queryHits(hitsNegativeAny)
    hitsPartial <- setdiff(hitsNegativeAny, hitsNegativeWithin)
    newAxts2 = list()
    if(length(hitsPartial) > 0L){
      newStarts <- pmax(start(queryRanges(x))[queryHits(hitsPartial)], 
                        (qSize-end(searchGRanges)+1)[subjectHits(hitsPartial)])
      newEnds <- pmin(end(queryRanges(x))[queryHits(hitsPartial)], 
                      (qSize-start(searchGRanges)+1)[subjectHits(hitsPartial)])
      for(i in 1:length(hitsPartial)){
        newAxts2[[i]] <- subAln(x[queryHits(hitsPartial[i])], 
                                newStarts[i], newEnds[i],
                                select="query")
      }
    }
    index <- sort(c(indexPositiveWithin, indexNegativeWithin))
    ans <- c(x[index], do.call(c, newAxts), do.call(c, newAxts2))
    return(ans)
  }
}



subAln <- function(x, start, end, select=c("target", "query")){
  if(length(x) != 1L){
    stop("subAln only operates on Axt of length 1")
  }
  select <- match.arg(select)
  alignedSeq1 <- strsplit(as.character(targetSeqs(x)[[1]]), "")[[1]]
  alignedSeq2 <- strsplit(as.character(querySeqs(x)[[1]]), "")[[1]]
  indexGap <- alignedSeq1 == "-" | alignedSeq1 == "." | alignedSeq1 == "_"
  seq12aln <- seq_len(length(alignedSeq1))[!indexGap]
  indexGap <- alignedSeq2 == "-" | alignedSeq2 == "." | alignedSeq1 == "_"
  seq22aln = seq_len(length(alignedSeq2))[!indexGap]
  if(select == "target"){
    alnStart <- seq12aln[start - start(targetRanges(x)) + 1L]
    alnEnd <- seq12aln[end - start(targetRanges(x)) + 1L]
  }else{
    alnStart <- seq22aln[start - start(queryRanges(x)) + 1L]
    alnEnd <- seq22aln[end - start(queryRanges(x)) + 1L]
  }
  while(alignedSeq1[alnStart] %in% c("-", ".", "_") ||
        alignedSeq2[alnStart] %in% c("-", ".", "_")){
    alnStart = alnStart + 1L
  }
  while(alignedSeq1[alnEnd] %in% c("-", ".", "_") ||
        alignedSeq2[alnEnd] %in% c("-", ".", "_")){
    alnEnd = alnEnd - 1L
  }
  ## TO DO: check newStart < newEnd
  newTargetStart <- which(alnStart == seq12aln) + start(targetRanges(x)) - 1L
  newTargetEnd <- which(alnEnd == seq12aln) + start(targetRanges(x)) - 1L
  newQueryStart <- which(alnStart == seq22aln) + start(queryRanges(x)) - 1L
  newQueryEnd <- which(alnEnd == seq22aln) + start(queryRanges(x)) - 1L
  ans <- Axt(targetRanges=GRanges(seqnames=seqnames(targetRanges(x)),
                                  ranges=IRanges(start=newTargetStart,
                                                 end=newTargetEnd),
                                  strand="+"),
             targetSeqs=DNAStringSet(paste0(alignedSeq1[alnStart:alnEnd],
                                            collapse="")),
             queryRanges=GRanges(seqnames=seqnames(queryRanges(x)),
                                 ranges=IRanges(start=newQueryStart,
                                                end=newQueryEnd),
                                 strand=strand(queryRanges(x))),
             querySeqs=DNAStringSet(paste0(alignedSeq2[alnStart:alnEnd],
                                           collapse="")),
             score=score(x), symCount=symCount(x))
  return(ans)
}