# Path: seqio.py
# def fastaIterator(fh):
# 
#     l = fh.readline()
#     if(not l or not l.startswith(">")):
#         raise Exception("No \">\" at start of Fasta File")
#     name = l.strip()[1:]
#     seq = ""
#     while True:
#         l = fh.readline()
#         if not l or l.startswith(">"):
#             yield FastaRecord(name,seq)
#             if not l:
#                 break
#             name = l.strip()[1:]
#             seq = ""
#         else:
#             seq += l.strip()
# 
# Path: cov.py
# 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)
# 
# Path: misc.py
# def create_enum(**enums):
#     return type('Enum', (), enums)


import sys
import copy
from nucio import *
from seqio import fastaIterator
from operator import itemgetter
from itertools import groupby, repeat, izip_longest, imap, count, chain
from collections import namedtuple
from cov import getMarkedRanges
from misc import create_enum
snp_it = lineRecordIterator(sfh, NucSNPRecord, NucSNPRecordTypes)


reads = dict(map(lambda r : (str(r.name), str(r.seq)), fastaIterator(rfh)))
alignments = dict(map(lambda (n,a): (n,list(a)), 
                      groupby(alignment_it, lambda x: x.sname)))

for pbname, snp_entries in groupby(snp_it, lambda x: x.sname):
    warnings = []
    pblen = len(reads[pbname])

    ##no alignments for this pb read
    if pbname not in alignments:
        continue
    ##create ranges of accepted alignments
    accept_alignment_ranges = [None] * pblen
    #alignments[pbname].sort(key=lambda a: (a.send-a.sstart) * pow(a.pctid/100.0,2))
    alignments[pbname].sort(key=lambda a: (a.send-a.sstart))
    for alignment in alignments[pbname]:
        for p in range(alignment.sstart-1,alignment.send):
            accept_alignment_ranges[p] = alignment.qname

    ##
    ##find clr ranges
    ##

    #find ranges
    covered_ranges = map(lambda (s,e): CoverageRange(s,e,1.0,CovStat["COVERED"]),
                         getMarkedRanges(map(lambda c: 1 if not c == None else 0 , accept_alignment_ranges)))
    uncovered_ranges = map(lambda (s,e): CoverageRange(s,e,0.7,CovStat["UNCOVERED"]),

# Path: seqio.py
# def fastaIterator(fh):
# 
#     l = fh.readline()
#     if(not l or not l.startswith(">")):
#         raise Exception("No \">\" at start of Fasta File")
#     name = l.strip()[1:]
#     seq = ""
#     while True:
#         l = fh.readline()
#         if not l or l.startswith(">"):
#             yield FastaRecord(name,seq)
#             if not l:
#                 break
#             name = l.strip()[1:]
#             seq = ""
#         else:
#             seq += l.strip()
# 
# Path: cov.py
# 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)
# 
# Path: misc.py
# def create_enum(**enums):
#     return type('Enum', (), enums)


# Path: pb_correct.py
import sys
import copy
from nucio import *
from seqio import fastaIterator
from operator import itemgetter
from itertools import groupby, repeat, izip_longest, imap, count, chain
from collections import namedtuple
from cov import getMarkedRanges
from misc import create_enum

snp_it = lineRecordIterator(sfh, NucSNPRecord, NucSNPRecordTypes)


reads = dict(map(lambda r : (str(r.name), str(r.seq)), fastaIterator(rfh)))
alignments = dict(map(lambda (n,a): (n,list(a)), 
                      groupby(alignment_it, lambda x: x.sname)))

for pbname, snp_entries in groupby(snp_it, lambda x: x.sname):
    warnings = []
    pblen = len(reads[pbname])

    ##no alignments for this pb read
    if pbname not in alignments:
        continue
    ##create ranges of accepted alignments
    accept_alignment_ranges = [None] * pblen
    #alignments[pbname].sort(key=lambda a: (a.send-a.sstart) * pow(a.pctid/100.0,2))
    alignments[pbname].sort(key=lambda a: (a.send-a.sstart))
    for alignment in alignments[pbname]:
        for p in range(alignment.sstart-1,alignment.send):
            accept_alignment_ranges[p] = alignment.qname

    ##
    ##find clr ranges
    ##

    #find ranges
    covered_ranges = map(lambda (s,e): CoverageRange(s,e,1.0,CovStat["COVERED"]),
                         getMarkedRanges(map(lambda c: 1 if not c == None else 0 , accept_alignment_ranges)))
    uncovered_ranges = map(lambda (s,e): CoverageRange(s,e,0.7,CovStat["UNCOVERED"]),

                           getMarkedRanges(map(lambda c: 1 if c == None else 0 , accept_alignment_ranges)))