# Path: nucio.py
# def getNucmerAlignmentIterator(fh):
#     '''Get nucmer alignments from show-coords output
#     (Deprecated legacy)
#     '''
#     return lineRecordIterator(fh, NucRecord, NucRecordTypes)
# 
# Path: cov.py
# def fillc(l,s,e):
#     '''l is an array
#     s and e are start and end of
#     a region the coverage should be incremented in'''
# 
#     if(len(l) < e):
#         l += list(repeat(0, e-len(l)))
#     for cr in range(s,e):
#         l[cr] += 1
# 
# def getMarkedRanges(v):
#     '''Takes an array v with wanted elements marked as 1
#     all other elements are 0. This function returns
#     the index ranges of these 1's'''
#     
#     #subtract adjacent elements to find -1
#     #flip these to 1's and change everything else to 0
#     #so that we can use it to find the index
#     breaks = map(lambda x: 1 if x == -1 else 0 , pairwise(v + [0],operator.sub))
# 
#     #create a cumsum of the inverse to know how many elements
#     #are a part of this region
#     lengths = accumulate_mod(v)
# 
#     #zip up the indexes with the lengths
#     z = zip(count(),list(lengths))
# 
#     #use the breaks to select only the indexes we want
#     #and use the cumsum to know how many elements came
#     #previously
#     endAndLength = compress(z,breaks)
# 
#     return map(lambda (e,l): (e-l+1,e), endAndLength)
# 
# def getCoverageFromNucAlignments(alignments):
#     '''Gets coverage from nucmer alignments'''
#     it = iter(alignments)
#     n = it.next()
# 
#     def nucAGetter():
#         yield (n.sstart-1, n.send-1)
#         for g in alignments:
#             yield (n.sstart-1, n.send-1)
# 
#     return getCoverage( nucAGetter(), n.slen)


import sys
from itertools import groupby,ifilter,izip_longest,count
from nucio import getNucmerAlignmentIterator
from cov import fillc, getMarkedRanges, getCoverageFromNucAlignments

fh = open(sys.argv[1])

COV_GAP_MIN = int(sys.argv[2])

fhist = open(sys.argv[3] + ".uncov.hist", "w")
freg = open(sys.argv[3] + ".uncov.regions", "w")
ftbases = open(sys.argv[3] + ".uncov.total.bases", "w")

pcov = [] #pct cov

total_bases = 0
total_uncovered_bases = 0

for pbname,alignments in groupby(getNucmerAlignmentIterator(fh), lambda x: x.sname):
    a = list(alignments)
    cov = getCoverageFromNucAlignments(a)

    #mark the 0 coverage regions
    zcov = map(lambda c: 1 if c == 0 else 0, cov)
    
    #ranges with 0 coverage
    zcov_ranges = getMarkedRanges(zcov)
    
    #only look at the gaps larger than the min gap size
    zcov_ranges_filt = filter(lambda (x,y) : y-x > COV_GAP_MIN ,zcov_ranges )

    #write out the regions that pass filter to region file
    freg.write("\t".join([pbname, " ".join(map(lambda t: "%d,%d" % t, zcov_ranges_filt))]) + "\n")


# Path: nucio.py
# def getNucmerAlignmentIterator(fh):
#     '''Get nucmer alignments from show-coords output
#     (Deprecated legacy)
#     '''
#     return lineRecordIterator(fh, NucRecord, NucRecordTypes)
# 
# Path: cov.py
# def fillc(l,s,e):
#     '''l is an array
#     s and e are start and end of
#     a region the coverage should be incremented in'''
# 
#     if(len(l) < e):
#         l += list(repeat(0, e-len(l)))
#     for cr in range(s,e):
#         l[cr] += 1
# 
# def getMarkedRanges(v):
#     '''Takes an array v with wanted elements marked as 1
#     all other elements are 0. This function returns
#     the index ranges of these 1's'''
#     
#     #subtract adjacent elements to find -1
#     #flip these to 1's and change everything else to 0
#     #so that we can use it to find the index
#     breaks = map(lambda x: 1 if x == -1 else 0 , pairwise(v + [0],operator.sub))
# 
#     #create a cumsum of the inverse to know how many elements
#     #are a part of this region
#     lengths = accumulate_mod(v)
# 
#     #zip up the indexes with the lengths
#     z = zip(count(),list(lengths))
# 
#     #use the breaks to select only the indexes we want
#     #and use the cumsum to know how many elements came
#     #previously
#     endAndLength = compress(z,breaks)
# 
#     return map(lambda (e,l): (e-l+1,e), endAndLength)
# 
# def getCoverageFromNucAlignments(alignments):
#     '''Gets coverage from nucmer alignments'''
#     it = iter(alignments)
#     n = it.next()
# 
#     def nucAGetter():
#         yield (n.sstart-1, n.send-1)
#         for g in alignments:
#             yield (n.sstart-1, n.send-1)
# 
#     return getCoverage( nucAGetter(), n.slen)


# Path: pb_uncovered.py
import sys
from itertools import groupby,ifilter,izip_longest,count
from nucio import getNucmerAlignmentIterator
from cov import fillc, getMarkedRanges, getCoverageFromNucAlignments


fh = open(sys.argv[1])

COV_GAP_MIN = int(sys.argv[2])

fhist = open(sys.argv[3] + ".uncov.hist", "w")
freg = open(sys.argv[3] + ".uncov.regions", "w")
ftbases = open(sys.argv[3] + ".uncov.total.bases", "w")

pcov = [] #pct cov

total_bases = 0
total_uncovered_bases = 0

for pbname,alignments in groupby(getNucmerAlignmentIterator(fh), lambda x: x.sname):
    a = list(alignments)
    cov = getCoverageFromNucAlignments(a)

    #mark the 0 coverage regions
    zcov = map(lambda c: 1 if c == 0 else 0, cov)
    
    #ranges with 0 coverage
    zcov_ranges = getMarkedRanges(zcov)
    
    #only look at the gaps larger than the min gap size
    zcov_ranges_filt = filter(lambda (x,y) : y-x > COV_GAP_MIN ,zcov_ranges )

    #write out the regions that pass filter to region file
    freg.write("\t".join([pbname, " ".join(map(lambda t: "%d,%d" % t, zcov_ranges_filt))]) + "\n")


    total_bases += a[0].slen