<gh_stars>0
#!/usr/bin/env python3

#
# Copyright 2015 Dovetail Genomics LLC
#
#

import os
import re
import pysam
from bamtags import BamTags
import gc
import random
import bisect
import numpy as np

def parse_pline(l,scaffold_lengths):
    c=l.strip().split()
    m=re.match("(.*)_(\d+)\.([35])",c[2])
    scaffold_lengths[c[1]]=int(c[6])
    return( c[1],m.group(1),int(m.group(2)),m.group(3),int(c[3]),int(c[7]) )

def parse_pline2(l,scaffold_lengths):
    c=l.strip().split()
    base=int(c[3])
    x=int(c[5])
    scaffold_lengths[c[1]]=max(scaffold_lengths.get(c[1],0),x)
    return( c[1],c[2],base,c[4],x,int(c[6]) )

def aln2bins(aln,binsize):        
    return( ( aln.tid,int(aln.pos/binsize) ),( aln.rnext,int(aln.pnext/binsize)))


class GrowingNPArray:
    def __init__(self,dtype):
        self.data = np.zeros((10000,),dtype=dtype)
        self.N=10000
        self.n=0
        self.dtype = dtype

    def append(self,d):
        if self.n==self.N:
            self.N *=2
            dd = np.zeros((self.N,),dtype=self.dtype)
            dd[:self.n] = self.data
            self.data = dd
        self.data[self.n]=d
        self.n+=1

    def finalize(self):
        dd = np.zeros((self.n,),dtype=self.dtype)
        dd[:self.n] = self.data[:self.n]
        return dd


def pairs_overlap(x,y):
    a=min(x[0],x[1])
    b=max(x[0],x[1])
    c=min(y[0],y[1])
    d=max(y[0],y[1])

    if a<=c and c<=b: return True
    if a<=d and d<=b: return True
    if c<=a and a<=d: return True
    if c<=b and b<=d: return True
    return False
                    
class HiriseAssembly:
    """This class is designed to encapsulate a Hirise Assembly."""

    def __init__(self,options={}):
        """Options currently used are:  bams, bam, datamodel, layout"""
        self.bams=[]
        if options.get("bams"): self.bams+=list(options["bams"])
        if options.get("bam") : self.bams +=[options["bam"]]
        self.model_params=False

        self.binsize=10000
        self.bam_objects=False
        self.shotgun_bam_objects=False

        if "bam" in options or "bams" in options:
            self.load_ocontig_info()
        
        if "datamodel" in options and options["datamodel"]: 
            self.set_modelparams(options["datamodel"])

        self.layout_lines=[]
        if options.get("layout"):
            self.load_playout(options["layout"])
        self.masked_regions={}
        self.check_bams()
        self.ocontig_fa=False
        self.ocontig_fasta=False
        if "shotgun" in options: 
            self.shotgun = options["shotgun"]
        else:
            self.shotgun = []
        self.contigx=False

    def get_seq(self,contig):
        if not self.ocontig_fasta:
            raise Exception 
        
        if not os.path.exists(self.ocontig_fasta+".fai"):
            pysam.faidx(self.ocontig_fasta)
        if not self.ocontig_fa:
            self.ocontig_fa = pysam.Fastafile(self.ocontig_fasta)

        m=re.match("(.*)_(\d+)",contig)
        ocontig=m.group(1)
        x=int(m.group(2))
        y=x+self.contig_length(contig)
        leng=self.ocontig_lengths[ocontig]
        #return self.ocontig_lengths[m.group(1)]


        #dummy,scaffold,x,y,leng = l.strip().split()
        #x=int(x)
        #y=int(y)
        #leng=int(leng)
#        y=min(y,leng)
        seq = self.ocontig_fa.fetch(ocontig,x-1,y)
        return(seq)

    def contig_length(self,contig):
        if contig in self.contig_lengths:
            return self.contig_lengths[contig]
        else:
            m=re.match("(.*)_(\d+)",contig)
            return self.ocontig_lengths[m.group(1)]
        return 0

    def contig_end_to_region(self,contig_end_string,w):
        ((this_contig,this_base,this_end),) = re.findall("^(.*)_(\d+)\.([53])$",contig_end_string)
        print(this_contig,this_base,this_end)
        this_base = int(this_base)
#        this_end = int(this_end)
        for scaffold,contig,base,end,x,l in self.layout_lines:
            slen = self.scaffold_lengths[scaffold]
            if this_contig == contig and this_base == base and this_end == end:
                return(scaffold,max(0,x-w),min(slen,x+w),(x<slen-x))
        print("so such end found",contig_end_string)
        raise Exception

    def index_ocontigs(self):
        self.ocontig_bases  = {}
        self.contig_lengths = {}
        self.contig_scaffold= {}
        self.load_ocontig_info()
#        self.scaffold_ocontigs={}
        
        ocb=self.ocontig_bases
        for scaffold,contig,base,end,x,l in self.layout_lines:

#            if not scaffold in self.scaffold_ocontigs:
#                self.scaffold_ocontigs[scaffold]={}
#            self.scaffold_ocontigs[scaffold][self.contig_ocontig[contig]]=1

            self.contig_scaffold[contig,base]=scaffold
            if end in [5,"5"] and not base==1:
                ocb[contig] = ocb.get(contig,[])+[ base ]
            if end in [5,"5"]:  
                self.contig_lengths[ contig+"_"+str(base) ]=l
                self.contig_lengths[ contig,base ]=l
        for b in ocb.keys():
            ocb[b].sort()
            #print("ocb:",b,ocb[b])
        self.total_ocontig_length = sum(self.ocontig_lengths.values())
        self.rc_buffer=[]
        x=0
        for i in self.ocontigs_iter():
            self.rc_buffer.append((x,i))
            x+=self.ocontig_lengths[i]

    def broken_coord(self,ocontig,x,tuples=False):
        if tuples:
            if not ocontig in self.ocontig_bases:
                if (not (ocontig,1) in self.contig_lengths ) or x < self.contig_lengths[(ocontig,1)]: 
#                    print(ocontig,"1")
                    return((ocontig,1), x)
                else:
#                    print(ocontig,"2")
                    return((ocontig, self.contig_lengths[(ocontig,1)]), x)
            for b in self.ocontig_bases[ocontig][::-1]:
                if b<x:
                    new_contig = (ocontig,b)
                    #print("zxzx",ocontig,x,b,new_contig,self.contig_lengths.get(new_contig,False))
                    if (not new_contig in self.contig_lengths) or (x-b)<self.contig_lengths[new_contig]:
#                        print(ocontig,"3")
                        return((ocontig,b), x-b)
                    else:
                        return (False,False)
            new_contig = (ocontig,1)
            if (not new_contig in self.contig_lengths) or x < self.contig_lengths[new_contig]:
#                print(ocontig,"4",new_contig,new_contig in self.contig_lengths)
                return(new_contig, x)
            else:
                return (False,False)
            

        else:
            if not ocontig in self.ocontig_bases:
                if (not ocontig+"_1" in self.contig_lengths ) or x < self.contig_lengths[ocontig+"_1"]: 
                    return(ocontig+"_1", x)
                else:
                    return(ocontig+"_"+str( self.contig_lengths[ocontig+"_1"] ), x)
            for b in self.ocontig_bases[ocontig][::-1]:
                if b<x:
                    new_contig = ocontig+"_"+str(b)
                    #print("zxzx",ocontig,x,b,new_contig,self.contig_lengths.get(new_contig,False))
                    if (not new_contig in self.contig_lengths) or (x-b)<self.contig_lengths[new_contig]:
                        return(ocontig+"_"+str(b), x-b)
                    else:
                        return (False,False)
            new_contig = ocontig+"_1"
            if (not new_contig in self.contig_lengths) or x < self.contig_lengths[new_contig]:
                return(new_contig, x)
            else:
                return (False,False)

    def merge_masked_regions(self,debug=False):
        mr = self.masked_regions
        for ocontig in mr.keys():
            new_regions=[]
            edges = []
            for x,y in mr[ocontig]:
                edges.append((x,1))
                edges.append((y,-1))
            edges.sort(key=lambda x:(x[0],-x[1]))
            rs=0
            state=0
            for x,b in edges:
                rs+=b
                if state==0 and rs>0:
                    start=x
                    state=1
                elif state==1 and rs==0:
                    end=x
                    new_regions.append((start,end))
                    state=0
            if debug and not (len(mr[ocontig])==len(new_regions)):
                print("#mmr",ocontig,len(mr[ocontig]),len(new_regions),mr[ocontig],new_regions)
            mr[ocontig]=new_regions

    def add_mask_regions(self,segments=[],filename=False):
        masked_regions = self.masked_regions
        if filename:
            for line in open(filename):
                c=line.strip().split()
                ocontig,x,y = c[0],int(c[1]),int(c[2])
                if not ocontig in masked_regions: masked_regions[ocontig]=[]
                masked_regions[ocontig].append((x,y))
        for ocontig,x,y in segments:
            if not ocontig in masked_regions: masked_regions[ocontig]=[]
            masked_regions[ocontig].append((x,y))
            
        for oc in masked_regions.keys():
            masked_regions[oc].sort()

    def setup_mapper(self,debug=False):
#        debug=True
        self.contigx={}
        self.contigst={}
        cex={}
        for scaffold,contig,base,end,x,l in self.layout_lines:
            if end in [5,"5"]:
                self.contigx[contig,base] = x
            cex[contig,base,end]=x
        for contig,base in self.contigx.keys():
            if cex[contig,base,"5"] < cex[contig,base,"3"]: 
                self.contigst[contig,base]=1
            else: 
                self.contigst[contig,base]=-1
#        if debug:
#            for contig,base in self.contigx.keys():
#                print("m:",contig,base,self.contigx[contig,base],self.contigst[contig,base])


    def scaffolds_iter(self,debug=False):
        if not self.layout_lines or len(self.layout_lines)==0:
            self.make_trivial_layout_lines(debug=debug)
            self.validate()

        last_scaffold=False
        for scaffold,contig,base,end,x,l in self.layout_lines:
            if (not last_scaffold) or not scaffold==last_scaffold:
                last_scaffold=scaffold
#                print("#",scaffold)
                yield(scaffold)

    def scaffold_coord(self,c,x,debug=False):
        if not self.contigx: self.setup_mapper(debug=debug)
        y=self.contigx[c]+self.contigst[c]*x
        if debug: print("sc:",c,x,y)
        return y

    def chicago_support_scan(self,scaffold,minsupport=0.0,mapq=10,debug=False,minx=0,maxx=1e6,gpf=False,logfile=False):

        import chicago_edge_scores as ces
        ces.set_exp_insert_size_dist_fit_params(self.model_params)
        model=ces.model

        links={}
        self.get_scaffold_links(scaffold,skipI=False,mapq=mapq,links=links,debug=debug,tuples=True)
#        model = ces.
        pairs=[]
        edges=[]
        buff=100
        nb=[]
        joins=[]
        t1=50.0
        t2=minsupport
        slen=self.ocontig_lengths[scaffold]
        for c1,c2 in links.keys():
            if debug: print(c1,c2,links[c1,c2])
            for x,y in links[c1,c2]:
                pairs.append( ( self.scaffold_coord(c1,x,debug=debug),self.scaffold_coord(c2,y,debug=debug)) )
        for a,b in pairs:               
            a,b=min(a,b),max(a,b)
            if debug:
                print("z",a,b)


            if (b-a)>minx and (b-a)<maxx:
                if gpf: gpf.write("{}\t{}\t{}\t{}\n".format(0.5*(a+b),0.5*(b-a),w,z))
                ll=model.lnF(b-a)

                edges.append( tuple([a+buff         , ll,c1,"a"]) )
                edges.append( tuple([b-buff         ,-ll,c2,"b"]) )

                nb.append( tuple([a+buff         , 1]) )
                nb.append( tuple([b-buff         ,-1]) )

        if gpf: gpf.write("\n\n\n")
        edges.sort()
        nb.sort()
        rs=0
        n=0
        tripped=False
        last=0
        state=0
        stretch_start=0
        low_point=0.0
        ji=0
        gap_scores={}
        gap_coverages={}
        max_t1=0
        min_t1=slen
        breakers=[]
        
        for i in range(len(edges)):
            rs+=edges[i][1]
            n+=nb[i][1]
            x=edges[i][0]

            try:
                score=model.cutScore(slen,x,n,rs,rangeCutoff=maxx,minSep=minx)
            except Exception as e:
                print(edges)
                print(i,edges[i],nb[i],rs,n,x,scaffold,slen,len(edges))
                raise e


            if score>t1:
                tripped=True
                last_t1=x
                min_t1 = min(min_t1,x)
                max_t1 = max(max_t1,x)
            if tripped and score<t2 and state==0:
                stretch_start=x
                state=1
                low_point =score
                minx = x
            if state==1 and score>t2:
                breakers.append((scaffold,stretch_start,x,minx,low_point,slen))
#                if logfile: logfile.write("{} {} {} {} {} {}\n".format(scaffold,stretch_start,x,minx,low_point,slen))
                state=0
            if state==1:
                if score<low_point:
                    low_point=score
                    minx=x
            if debug: print("dd:",scaffold,edges[i][0],score,rs,state,x,stretch_start,score,n,edges[i][2])
            if gpf: gpf.write("{}\t{}\t{}\t{}\n".format(edges[i][0],score,rs,n))

            last=edges[i][0]
            
        for scaffold,stretch_start,x,minx,low_point,slen in breakers:
            if stretch_start < last_t1:
#                if logfile: logfile.write("{} {} {} {} {} {} {} {}\n".format(scaffold,stretch_start,x,minx,low_point,slen,min_t1,max_t1))
                if logfile: logfile.write("{} {} {} {} {} {}\n".format(scaffold,stretch_start,x,minx,low_point,slen))

    def make_trivial_layout_lines(self,debug=False):
        self.layout_lines=[]
        #print(self.ocontig_lengths)
        if len(self.ocontig_lengths)==0:
            print("#load_contig_info")
            self.load_ocontig_info()
        i=1
        for ocontig in self.ocontigs_iter():
            scaffold=ocontig
#            if debug: print(scaffold,ocontig,1,"5",0,self.ocontig_lengths[ocontig],sep="\t")
            self.layout_lines.append((scaffold,ocontig,1,"5",0,self.ocontig_lengths[ocontig])) #
            self.layout_lines.append((scaffold,ocontig,1,"3",self.ocontig_lengths[ocontig],self.ocontig_lengths[ocontig])) #
            self.scaffold_lengths[scaffold] = self.ocontig_lengths[ocontig]
            self.scaffold_ocontigs[scaffold] = self.scaffold_ocontigs.get(scaffold,[])+[ocontig]
            i+=1
        self.index_ocontigs()
        self.setup_mapper()

    def get_scaffold_links(self,Tscaffold,skipI=False,mapq=10,links={},debug=False,tuples=False):
        if not self.layout_lines or len(self.layout_lines)==0:
            self.make_trivial_layout_lines(debug=debug)
            self.validate()
        contigs={}
        ocontigs={}
        if not tuples:
            for scaffold,contig,base,end,x,l in self.layout_lines:
                if scaffold==Tscaffold and end in [5,"5"]:
                    contigs[ contig+"_"+str(base) ]=1  
                    ocontigs[contig]=1
        else:
            for scaffold,contig,base,end,x,l in self.layout_lines:
                if scaffold==Tscaffold and end in [5,"5"]:
                    contigs[ contig,base ]=1  
                    ocontigs[contig]=1
                
        if debug: print(ocontigs,contigs)
        self.get_links(list(ocontigs.keys()),skipI,mapq,links,contigs,tuples=tuples)

    def segments_iter(self):
        for l in self.layout_lines:
            scaffold1,contig1,base1,end1,x1,l21 =l
            if end1 in [5 , "5"]:
                #cc=self.contig_length(contig1+"_"+str(base1))
                yield( contig1,base1,base1+l21,self.ocontig_lengths[contig1] )

    def contigs_iter(self):
        for l in self.layout_lines:
            scaffold1,contig1,base1,end1,x1,l21 =l
            if end1 in [5 , "5"]:
                yield( contig1 +"_"+str(base1) )

    def ocontigs_iter(self):
        for ocontig in self.ocontig_lengths.keys():
            yield(ocontig)

    def load_ocontig_info(self):
        import pysam

        self.ocontig_lengths={}
        if not self.bam_objects:
            self.bam_objects={}
            for b in self.bams:
                self.bam_objects[b] =  pysam.Samfile(b,"rb") 
        for bam in self.bam_objects.values():
            h=bam.header
            seqs=h['SQ']
            for s in seqs:
                if s['SN'] in self.ocontig_lengths:
                    if not  self.ocontig_lengths[s['SN']]==s['LN'] :
                        raise Exception 
                else:
                    self.ocontig_lengths[s['SN']]=s['LN']

    def random_window(self,wlen=1000):
        r=int(random.random()*self.total_ocontig_length)
        x=bisect.bisect_right(self.rc_buffer,(r,"a"))
        while r+x > self.total_ocontig_length or (x<len(self.rc_buffer)-1 and r+wlen>self.rc_buffer[x+1][0]):
            r=int(random.random()*self.total_ocontig_length)
            x=bisect.bisect_right(self.rc_buffer,(r,"a"))
#            print( r,self.rc_buffer[x],self.rc_buffer[x+1],r-self.rc_buffer[x][0])
        x-=1
#        print( r,self.rc_buffer[x],self.rc_buffer[x+1],r-self.rc_buffer[x][0])
        return( self.rc_buffer[x][1], r-self.rc_buffer[x][0], r-self.rc_buffer[x][0]+wlen )

    

    def top_other_contigs(self,reference,mapq=10,n=2,debug=False,bins=False,binx=0):
        nother={}
        if not self.bam_objects:
            self.bam_objects={}
            for b in self.bams:
                self.bam_objects[b] = pysam.Samfile(b,"rb")
        if debug: print(self.bam_objects)

        if bins:
            region="{}:{}-{}".format(reference,binx*self.binsize,(binx+1)*self.binsize)
            for bam in self.bam_objects.values():
                for aln in bam.fetch(region=region,until_eof=True):
                    if aln.is_duplicate :              continue
                    if aln.mapq < mapq               : continue
                    if BamTags.mate_mapq(aln) < mapq : continue
                    my_bin,rbin = aln2bins(aln,self.binsize)
                    if not my_bin == rbin : 
                        nother[rbin] = nother.get(rbin,0)+1
        else:
            for bam in self.bam_objects.values():
                for aln in bam.fetch(reference=reference,until_eof=True):
                    if aln.is_duplicate :              continue
                    if aln.mapq < mapq               : continue
                    if BamTags.mate_mapq(aln) < mapq : continue
                    if not aln.rnext == aln.tid : 
                        nother[aln.rnext] = nother.get(aln.rnext,0)+1
 
        kk=list(nother.keys())
        kk.sort(key=lambda x: nother[x],reverse=True)
        if debug: print(kk[:2])
        return(tuple(kk[:2]))

    def window_stats(self,reference,xa,xb,mapq=10,debug=False,bins=False):
        if not self.bam_objects:
            self.bam_objects=[]
            for b in self.bams:
                self.bam_objects.append( pysam.Samfile(b,"rb") )

        if not self.shotgun_bam_objects:
            self.shotgun_bam_objects=[]
            for b in self.shotgun:
                self.shotgun_bam_objects.append( pysam.Samfile(b,"rb") )
        
        binx = 0
        if bins:
            binx = int(xa/self.binsize)

        top2 = self.top_other_contigs(reference,mapq=mapq,debug=debug,bins=bins,binx=binx)
        region = "{}:{}-{}".format(reference,max(1,xa-1000),xb+1000)
        n_aligned=0
        n_hiq=0
        templates={}
        n_lowq   =0
        occ ={}
        occ2={}
        for bam in self.bam_objects.values():
            if debug: print("#",region,xa,xb,bam)
            for aln in bam.fetch(region=region):
                if aln.pos > xb:                   continue
                if aln.pos < xa:                   continue
                if aln.is_duplicate :              continue
                n_aligned+=1
                if aln.mapq < mapq               : continue
                if BamTags.mate_mapq(aln) < mapq : continue
                templates[aln.query_name]=templates.get(aln.query_name,0)+1
                if bins:
                    my_bin,rbin = aln2bins(aln,self.binsize)
                    if not my_bin==rbin:
                        occ[rbin] = occ.get(rbin,0)+1
                        if not rbin in top2:
                            occ2[rbin] = occ2.get(rbin,0)+1
                            
                else:
                    if (not aln.rnext == aln.tid):
                        occ[aln.rnext] = occ.get(aln.rnext,0)+1
                        if (not aln.rnext in top2) : 
                            occ2[aln.rnext] = occ2.get(aln.rnext,0)+1
                n_hiq+=1
            if debug: print("#",n_hiq,n_aligned)
        n_both   = len([ i for i in templates.keys() if templates[i]>1  ])
        n_other  = len(occ.keys())
        n_other2 = len([ i for i in occ.keys() if occ[i]>1  ])
        n_other3 = len([ i for i in occ.keys() if occ[i]>2  ])

        n_otherB  = len(occ2.keys())
        n_otherB2 = len([ i for i in occ2.keys() if occ[i]>1  ])
        n_otherB3 = len([ i for i in occ2.keys() if occ[i]>2  ])

        n_shotgun=0
        for bam in self.shotgun_bam_objects:
            for aln in bam.fetch(region=region):
                if aln.pos > xb: continue
                if aln.is_duplicate :              continue
                if aln.mapq < mapq               : continue
                n_shotgun+=1
        
        print("\t".join(map(str,[region,n_aligned,n_hiq,n_both,n_other,n_other2,n_other3,n_otherB,n_otherB2,n_otherB3,n_shotgun])))

    def chicago_promiscuity_mask(self,ocontig,mapq=10,w=1000,minlinks=2,maxothers=3,outfile=False,bins=False,debug=False):
        if not self.bam_objects:
            self.bam_objects=[]
            for b in self.bams:
                self.bam_objects.append( pysam.Samfile(b,"rb") )
#        dtype=[('x',int),('t',int),('y',int),('b',int),('r',int)]
        
        dtype=[('x',int),('t',int),('y',int)]
        buffer = GrowingNPArray(dtype=dtype)

        i=0
        tid=0
        for bam in self.bam_objects:                
            if debug: print("#",ocontig,bam)
            tid = bam.gettid(ocontig)
            for aln in bam.fetch(reference=ocontig):
                if aln.rnext<0: continue
                if aln.mapq < mapq               : continue
                if BamTags.mate_mapq(aln) < mapq : continue
                if aln.is_duplicate : continue
#                buffer.append((aln.pos,aln.rnext,aln.pnext,i,aln.tid))
                buffer.append((aln.pos,aln.rnext,aln.pnext))
            i+=1

        buffer = buffer.finalize()
        buffer = np.sort(buffer,order='x',kind="mergesort")

        bad_edges=[]


        def get_contig(x,b):
            return x[1]

        def get_contig_bin(x,b):
            return(x[1]+(int(x[2]/b)/1.0e6))

        if bins:
            get_bin = get_contig_bin
        else:
            get_bin = get_contig

        bc={}
        i=0
        j=0
        N=len(buffer)
        while i < N:
            mybin = get_bin((0,tid,buffer[i][0]),self.binsize)
            while j<N and buffer[j][0]<buffer[i][0]+w: 
                jbin = get_bin(buffer[j],self.binsize)
                bc[ jbin ] = bc.get(jbin,0)+1
                j+=1
                
            n_in_window = j-i
            
            nw=len([ ii for ii in bc.keys() if bc[ii]>=minlinks and not ii==mybin ])  #n_other_windows(buffer[i:j],minlinks,bins,self.binsize)

#            if debug: print(tid,buffer[i][0],n_in_window,nw,mybin,[ (ii,bc[ii]) for ii in bc.keys() if bc[ii]>=minlinks and not ii==mybin ],[ get_bin(x,self.binsize) for x in buffer[i:j] ])
            if debug: print(tid,buffer[i][0],n_in_window,nw,mybin,[ (ii,bc[ii]) for ii in bc.keys() if bc[ii]>=minlinks and not ii==mybin ])


            if nw > maxothers:
                bad_edges.append((buffer[i  ][0], 1))
                bad_edges.append((buffer[j-1][0],-1))

            ibin = get_bin(buffer[i],self.binsize)
            bc[ ibin ] = bc[ibin]-1
            if bc[ibin]==0: bc.pop(ibin)
            i+=1

        bad_edges.sort()
        if debug:
            for be in bad_edges:
                print(e)

        rs=0
        state=0
        for x,dy in bad_edges:
            rs+=dy
            if state==0 and rs>0:
                startx=x
                state=1
            elif state==1 and rs==0:
                outfile.write("{} {} {}\n".format(ocontig,startx,x))
                if debug: print("#{} {} {}\n".format(ocontig,startx,x))
                state=0
        outfile.flush()

        
    def chicago_read_density_mask(self,ocontigs,mapq=10,outfile=False,w=1000,cutoff=100,debug=False,shotgunT=False):

        if not self.bam_objects:
            self.bam_objects=[]
            for b in self.bams:
                self.bam_objects.append( pysam.Samfile(b,"rb") )

        if not self.shotgun_bam_objects:
            self.shotgun_bam_objects=[]
            for b in self.shotgun:
                self.shotgun_bam_objects.append( pysam.Samfile(b,"rb") )

        for ocontig in ocontigs:
            locations=[]
            shotgun_locations=[]
#            gc.disable()
            for bam in self.bam_objects:                
                if debug: print("#",ocontig,bam)
                for aln in bam.fetch(reference=ocontig):
#                    if debug: print("#",aln.pos,aln.is_duplicate,aln.mapq,BamTags.mate_mapq(aln),aln.query_name,len(locations))
                    if aln.rnext<0: continue
                    if aln.mapq < mapq               : continue
                    if BamTags.mate_mapq(aln) < mapq : continue
                    if aln.is_duplicate : continue
                    locations.append(aln.pos)
                if debug: print("#",len(locations))
            if shotgunT:
                for bam in self.shotgun_bam_objects:                
                    if debug: print("#",ocontig,bam)
                    for aln in bam.fetch(reference=ocontig):
                        if aln.mapq < mapq  : continue
                        if aln.is_duplicate : continue
                        shotgun_locations.append(aln.pos)
                    if debug: print("#",len(locations),len(shotgun_locations))

#            gc.enable()
            locations.sort()
            shotgun_locations.sort()

            i=0
            j=0
            bad_edges=[]

            while i<len(locations):
                while j<len(locations) and locations[j]<locations[i]+w: j+=1
                n_in_window = j-i
                if debug:
                    print("chicago_density",ocontig,locations[i],locations[j-1],j-i,locations[j-1]-locations[i],locations[i:j])
                if n_in_window > cutoff:
    #                outfile.write("{} {} {}\n".format(ocontig,locations[i],locations[j-1]))
                    bad_edges.append((locations[i],1)  )
                    bad_edges.append((locations[j-1],-1)  )
                i+=1

            if shotgunT:
                i=0
                j=0
                while i<len(shotgun_locations):
                    while j<len(shotgun_locations) and shotgun_locations[j]<shotgun_locations[i]+w: j+=1
                    n_in_window = j-i
                    if debug:
                        print("shotgun_density",ocontig,shotgun_locations[i],shotgun_locations[j-1],j-i,shotgun_locations[j-1]-shotgun_locations[i])
                    if n_in_window > shotgunT:
        #                outfile.write("{} {} {}\n".format(ocontig,locations[i],locations[j-1]))
                        bad_edges.append((shotgun_locations[i],1)  )
                        bad_edges.append((shotgun_locations[j-1],-1)  )
                    i+=1


            bad_edges.sort()

            rs=0
            state=0
            for x,dy in bad_edges:
                rs+=dy
                if state==0 and rs>0:
                    startx=x
                    state=1
                elif state==1 and rs==0:
                    outfile.write("{} {} {}\n".format(ocontig,startx,x))
                    state=0
        outfile.flush()


    def make_scaffold_mask(self,scaffold): #,mapper,mask):
        s_mask = {}
        segments={}
        slen = self.scaffold_lengths[scaffold]

        for ocontig in self.scaffold_ocontigs[scaffold]:
#            print(ocontig)
            for a,b in self.masked_regions.get(ocontig,[]):
                
                c1,x = self.broken_coord(ocontig,a,tuples=True)
                c2,y = self.broken_coord(ocontig,b,tuples=True)

                s1=self.contig_scaffold.get(c1)
                s2=self.contig_scaffold.get(c2)

                if c1==c2 and s1==s2 and s1==scaffold:

                    xx=self.scaffold_coord(c1,x)
                    yy=self.scaffold_coord(c2,y)
                    #print("#msm:",scaffold,ocontig,a,b,x,y,xx,yy,s1,s2,scaffold,sep="\t")

                    segments[ max(0,min(xx,yy)) , min( slen, max(xx,yy) )   ]=1
                elif (not c1==c2) and s1==s2 and s1==scaffold:
                    
                    xA = x
                    yA = self.contig_lengths[c1]
                    xx=self.scaffold_coord(c1,xA)
                    yy=self.scaffold_coord(c1,yA)
                    segments[ max(0,min(xx,yy)) , min( slen, max(xx,yy) )   ]=1

                    xB = 1
                    yB = y
                    xx=self.scaffold_coord(c2,xB)
                    yy=self.scaffold_coord(c2,yB)
                    segments[ max(0,min(xx,yy)) , min( slen, max(xx,yy) )   ]=1

                elif (not c1==c2) and s1==scaffold and not s2==scaffold:
                    
                    xA = x
                    yA = self.contig_lengths[c1]
                    xx=self.scaffold_coord(c1,xA)
                    yy=self.scaffold_coord(c1,yA)
                    segments[ max(0,min(xx,yy)) , min( slen, max(xx,yy) )   ]=1

                elif (not c1==c2) and s2==scaffold and not s1==scaffold:
                    
                    xB = 1
                    yB = y
                    xx=self.scaffold_coord(c2,xB)
                    yy=self.scaffold_coord(c2,yB)
                    segments[ max(0,min(xx,yy)) , min( slen, max(xx,yy) )   ]=1
 
        segments=list(segments.keys())
        segments.sort()
        #print("#msm:",scaffold,slen,segments,sep="\t")
        return(segments)

    def read_deserts(self,outfile,mapq=10,min_len=1000):
        if not self.bam_objects:
            self.bam_objects={}
            for b in self.bams:                
                self.bam_objects[b]= pysam.Samfile(b,"rb") 
        desert_edges=[]
        for bam in self.bam_objects.values():
            lastx=0
            lastr=-1
            for aln in bam.fetch(until_eof=True):
                if aln.mapq < mapq               : continue
                ref = aln.tid
                if not ref==lastr:
                    if lastr>-1:
                        lastname =bam.getrname(lastr)
                        lastl = self.ocontig_lengths[ lastname ]
                        if lastl - lastx > min_len:
                            desert_edges.append(( lastname,lastx  , 1 ))
                            desert_edges.append(( lastname,lastl  ,-1 ))
                    lastx=0
                deltax= aln.pos - lastx
                if deltax>min_len:
                    desert_edges.append(( bam.getrname(aln.tid),lastx  , 1 ))
                    desert_edges.append(( bam.getrname(aln.tid),aln.pos,-1 ))

                lastx=aln.pos
                lastr=aln.tid

        
        nbams = len( self.bams)
        state=0
        startx=0
        desert_edges.sort()
        y=0
        for contig,x,dy in desert_edges:

            y+=dy

            if state==0 and y == nbams:
                state=1
                startx=x
            elif state==1 and y<nbams:
                endx = x
                state=0
                if endx-startx>min_len:
                    print(contig,startx,endx,"readDesert",file=outfile)

            last_contig,last_x = contig,x


    def chicago_pairs(self,mapq=10,scaffold_contig_counts={},contig_scaffold_counts={},callback=False,bamfile=False,my_scaffolds=False):
        read_length = 100
        if not self.bam_objects:
            self.bam_objects={}
            for b in self.bams:                
                self.bam_objects[b]= pysam.Samfile(b,"rb") 
                


        if bamfile:
            bams_iter = iter( [ self.bam_objects[bamfile] ] )
        else:
            bams_iter = iter( self.bam_objects.values() )

        for bam in bams_iter:
#            print("#",bam)
            last_contig=False
            mri=0
            mr=[]
            last_scaffold=False
            for aln in bam.fetch(until_eof=True):
                if not aln.is_read1: continue
                track=False
                skip=False
                if aln.rnext<0:      continue
                if aln.mapq < mapq               : continue
                if BamTags.mate_mapq(aln) < mapq : continue
                if aln.is_duplicate : continue
                contig  = bam.getrname(aln.tid)
                contig2 = bam.getrname(aln.rnext)
                
                if not last_contig or not last_contig==contig:
                    if contig in self.masked_regions:
                        mr = self.masked_regions[contig]
                        mri=0
                    else:
                        mr=[]
                        mri=0
                c1,x = self.broken_coord(contig ,aln.pos  ,tuples=True)
                scaffold=self.contig_scaffold.get(c1)
                c2,y = self.broken_coord(contig2,aln.pnext,tuples=True)
                nscaffold=self.contig_scaffold.get(c2)

                if (not my_scaffolds==False) and (not ((scaffold in my_scaffolds) and (nscaffold in my_scaffolds))): continue

                if last_contig and not last_contig==contig:
#                    print("#done with contig",last_contig,list(contig_scaffold_counts.get(last_contig,{}).keys()))
                    for sc in contig_scaffold_counts.get(last_contig,{}).keys():  # for each hirise scaffold that this contig contributes a segment to:
                        if sc in scaffold_contig_counts:                   
                            if last_contig in scaffold_contig_counts[sc]:  # if we are still tracking the number of times we need to read this contig from a bam file for this scaffold
#                                print("##",sc,last_contig,scaffold_contig_counts[sc][last_contig])
                                scaffold_contig_counts[sc][last_contig] -= 1 # decrease the count by one
#                                print("##",sc,last_contig,scaffold_contig_counts[sc][last_contig])
                                if scaffold_contig_counts[sc][last_contig] ==0:
                                    del scaffold_contig_counts[sc][last_contig]
                            if len(scaffold_contig_counts[sc])==0:
#                                print("scaffold {} ready now",sc)
                                if callback: 
                                    callback(sc)
                                else:
                                    yield(sc,"done")
                                
                                del scaffold_contig_counts[sc]
#                                raise ScaffoldReady(sc)
#                            print("#s",last_contig,sc,scaffold_contig_counts.get(sc))
#                    print("#c",last_contig,contig_scaffold_counts.get(last_contig))
                last_contig = contig        
                last_scaffold = scaffold

                while mri<len(mr) and aln.pos > mr[mri][1]: mri+=1                        
                if mri<len(mr) and mr[mri][0]<aln.pos and aln.pos+read_length<mr[mri][1] : continue #this read maps to a masked region 

                if contig2 in self.masked_regions:
                    if [ (aa,bb) for aa,bb in self.masked_regions[contig2] if aa<aln.pnext and aln.pnext+read_length<bb ]: continue  #the pair maps to a masked region
                
                if not scaffold: continue
                if not nscaffold: continue

                xx=self.scaffold_coord(c1,x)
                yy=self.scaffold_coord(c2,y)

                #                yield(s1,self.scaffold_coord(c1,x),s2,self.scaffold_coord(c2,y))

                if self.scaffold_lengths[scaffold] < xx or  self.scaffold_lengths[nscaffold] < yy:
                    print("coordinate out of range error",xx,yy,scaffold,nscaffold,self.scaffold_lengths[scaffold], self.scaffold_lengths[nscaffold],contig,contig2,c1,c2,x,y)
                    raise Exception 


                yield( scaffold, nscaffold, xx, yy, c1, c2, aln.query_name )


    def chicago_pairs_for_scaffolds(self,mapq=10,callback=False,bamfile=False,scaffolds=[],contigs=False,minsep=0):
        read_length = 100
        if not self.bam_objects:
            self.bam_objects={}
            for b in self.bams:                
                self.bam_objects[b]= pysam.Samfile(b,"rb") 

        if bamfile:
            bams_iter = iter( [ self.bam_objects[bamfile] ] )
        else:
            bams_iter = iter( self.bam_objects.values() )


        ocontigs_to_hit = {}
        my_ocontigs=[]
        if contigs:
            my_ocontigs=contigs
        else:
            for scaffold in scaffolds:
                for oc in self.scaffold_ocontigs[scaffold]:
                    ocontigs_to_hit[oc]=1

            my_ocontigs = list(ocontigs_to_hit.keys())

        for bam in bams_iter:
#            print("#",bam)
            last_contig=False
            mri=0
            mr=[]
            last_scaffold=False
            for ocontig in my_ocontigs:
                for aln in bam.fetch(reference=ocontig):
                    if not aln.is_read1: continue
                    track=False
                    skip=False
                    if aln.rnext<0:      continue
                    if aln.mapq < mapq               : continue
                    if BamTags.mate_mapq(aln) < mapq : continue
                    if aln.is_duplicate : continue
                    contig  = bam.getrname(aln.tid)
                    contig2 = bam.getrname(aln.rnext)

                    if not last_contig or not last_contig==contig:
                        if contig in self.masked_regions:
                            mr = self.masked_regions[contig]
                            mri=0
                        else:
                            mr=[]
                            mri=0

                    if last_contig and not last_contig==contig:
                        pass
    #                    print("#done with contig",last_contig,list(contig_scaffold_counts.get(last_contig,{}).keys()))

                    c1,x = self.broken_coord(contig ,aln.pos  ,tuples=True)
                    scaffold=self.contig_scaffold.get(c1)
                    last_contig = contig        
                    last_scaffold = scaffold

                    while mri<len(mr) and aln.pos > mr[mri][1]: mri+=1                        
                    if mri<len(mr) and mr[mri][0]<aln.pos and aln.pos+read_length<mr[mri][1] : continue #this read maps to a masked region 

                    if contig2 in self.masked_regions:
                        if [ (aa,bb) for aa,bb in self.masked_regions[contig2] if aa<aln.pnext and aln.pnext+read_length<bb ]: continue  #the pair maps to a masked region

                    c2,y = self.broken_coord(contig2,aln.pnext,tuples=True)
                    nscaffold=self.contig_scaffold.get(c2)
                    if not scaffold: continue
                    if not nscaffold: continue

                    if not scaffold in scaffolds: continue
                    if not nscaffold in scaffolds: continue

                    xx=self.scaffold_coord(c1,x)
                    yy=self.scaffold_coord(c2,y)
                    if scaffold==nscaffold and abs(xx-yy)<minsep: continue

                    #                yield(s1,self.scaffold_coord(c1,x),s2,self.scaffold_coord(c2,y))

                    if self.scaffold_lengths[scaffold] < xx or  self.scaffold_lengths[nscaffold] < yy:
                        print("coordinate out of range error",xx,yy,scaffold,nscaffold,self.scaffold_lengths[scaffold], self.scaffold_lengths[nscaffold],contig,contig2,c1,c2,x,y)
                        raise Exception 

                    yield( scaffold, nscaffold, xx, yy, c1, c2, aln.query_name )

    def get_links(self,ocontigs,skipI=False,mapq=10,links={},contigs=False,raw=False,tuples=False,debug=False):
#        links={}

        if not self.bam_objects:
            self.bam_objects={}
            for b in self.bams:

                self.bam_objects[b]= pysam.Samfile(b,"rb") 

        for bam in self.bam_objects.values():
#            bam = pysam.Samfile(b,"rb") 

            for ocontig in ocontigs:
#                print("#",ocontig)
                mr=[]
                mri=0
                if ocontig in self.masked_regions: mr=self.masked_regions[ocontig]
                for aln in bam.fetch(reference=ocontig):
#                    if debug: print("#",aln.query_name,aln.pos,aln.pnext,aln.mapq,sep="\t")
                    if aln.rnext<0: continue
                    if aln.mapq < mapq               : continue
                    if BamTags.mate_mapq(aln) < mapq : continue
                    if aln.is_duplicate : continue

                    while mri<len(mr) and aln.pos > mr[mri][1]: mri+=1                        
                    if mri<len(mr) and mr[mri][0]<aln.pos and aln.pos<mr[mri][1] : 
                        if debug: print("#m1",ocontig,bam.getrname(aln.rnext),aln.pos,aln.pnext,mr[mri],aln.query_name,sep="\t")
                        continue #this read maps to a masked region 

                    if raw:
                        c1,x1 = ocontig,aln.pos
                    else:
                        c1,x1 = self.broken_coord(ocontig,aln.pos,tuples=tuples)
                        if not c1: continue

                    if contigs and not c1 in contigs: continue
                    ocontig2 = bam.getrname(aln.rnext)
                    
                    
                    if raw:
                        c2,x2 = ocontig2,aln.pnext
                    else:
                        c2,x2 = self.broken_coord(ocontig2,aln.pnext,tuples=tuples)
                        if not c2: continue
                    if contigs and not c2 in contigs: continue

                    if skipI and c1 == c2: continue

                    if ocontig2 in self.masked_regions:
                        overlapping_masked_segs = [ (aa,bb) for aa,bb in self.masked_regions[ocontig2] if aa<aln.pnext and aln.pnext<bb ]
                        if overlapping_masked_segs: 
                            if debug: print("#m2",ocontig,ocontig2, aln.pos,aln.pnext,overlapping_masked_segs ,aln.query_name,sep="\t")
                            continue  #the pair maps to a masked region

                    if not (c1,c2) in links: 
                        links[c1,c2]=[]
                    links[c1,c2].append((x1,x2))
                    if debug: print("#p",c1,c2,aln.query_name,aln.pos,aln.pnext,sep="\t")
                    #                print(self.broken_coord(ocontig,aln.pos),self.broken_coord(bam.getrname(aln.rnext),aln.pnext))
#        for c1,c2 in links:
#            print(c1,c2,links[c1,c2])
        return(links)

    def validate(self):
        contig_mins={}
        contig_maxs={}
        scaffold_mins={}
        scaffold_maxs={}
        contig_lengths={}
        counts={}
        #        scaffold_lengths={]
        lastx=0
        last_scaffold=False
        last_contig=(False,)
        for scaffold,contig,base,end,x,l in self.layout_lines:
            if last_scaffold and scaffold == last_scaffold and not x>lastx: 
                print(last_scaffold,scaffold,x,lastx)
                raise Exception 
            if x<0: raise Exception
            counts[contig,base]=counts.get((contig,base),0)+1
            contig_mins[contig,base]=min( contig_mins.get((contig,base),1e99) , x )
            contig_maxs[contig,base]=max( contig_maxs.get((contig,base),-1) , x )
            if (contig,base) in contig_lengths:
                if not contig_lengths[contig,base]==l: raise Exception 
            else:
                contig_lengths[contig,base]=l
            if (contig,base) == last_contig:
#                print(scaffold,contig,base,end,x,l,last_contig,lastx,x-lastx,"XXXX",sep="\t")
                if not l == x-lastx :
                    print(scaffold,contig,base,end,x,l,last_contig,lastx,sep="\t")
                    raise Exception 
            scaffold_mins[scaffold]=min( scaffold_mins.get(scaffold,1e99) , x )
            scaffold_maxs[scaffold]=max( scaffold_maxs.get(scaffold,-1)   , x )
            last_contig=(contig,base)
            lastx,last_scaffold=x,scaffold
        for contig,base in counts.keys():
            if not counts[contig,base]==2: raise Exception 
            if not contig_maxs[contig,base]-contig_mins[contig,base]==contig_lengths[contig,base]: 
                print("wtf?",contig,base,contig_mins[contig,base],contig_maxs[contig,base],contig_lengths[contig,base],contig_maxs[contig,base]-contig_mins[contig,base])
                raise Exception 

        self.setup_mapper()
        lastk=(1,0)
        lastl=0
        lastx=0
        for k in sorted(self.contigx.keys()):
            if k[0]==lastk[0]:
                #print("XX",k,contig_lengths.get(k),self.contigx[k],lastl,k[1]-lastk[1],lastl==k[1]-lastk[1],sep="\t")
                if not lastl==k[1]-lastk[1]:
                    print(k,contig_lengths.get(k),self.contigx[k],lastl,k[1]-lastk[1],lastl==k[1]-lastk[1],sep="\t")
                    raise Exception 
            lastk=k
            lastl = contig_lengths.get(k)
            lastx=self.contigx[k]
#        self.contigx={}
#        self.contigst={}


        print("#validation success")
#            contig_maxs[contig,base]=max( contig_maxs.get((contig,base),-1) , x )


    def add_breakpoints_ranges(self,breaks,debug=False,scores={},unlink_threshold=100.0, contig_break_threshold=100.0,score_types={},clipped_likelihood_for_gaps_only=True):
        """This function filters breaks expressed as spans according to several criteria before passing selected break points along to add_breakpoints()"""

        if not self.layout_lines or len(self.layout_lines)==0:
            self.make_trivial_layout_lines(debug=debug)
            self.validate()

        break_points = []
        category={}
        br={}
        scaffold_breaks={}
        for s,a,b,c in breaks:
            if not s in scaffold_breaks: scaffold_breaks[s]=[]

        for s,a,b,c in breaks:
            scaffold_breaks[s].append((a,b,c))

        for s in scaffold_breaks.keys():
             scaffold_breaks[s].sort()
        i=0
        offset={}
        while i< len(self.layout_lines)-1:
            #print(self.layout_lines[i])
            scaffold1,contig1,base1,end1,x1,l21 = self.layout_lines[  i  ]
            scaffold2,contig2,base2,end2,x2,l22 = self.layout_lines[ i+1 ]
            i+=1
            if not scaffold1 in scaffold_breaks: continue
            if not scaffold1 == scaffold2 : continue
            for a,b,c in scaffold_breaks[scaffold1]:
                if pairs_overlap((a,b),(x1,x2)) and (not (contig1,base1)==(contig2,base2)):  # you don't have to break a contig to make this break
                    if x1<=c and c<=x2:
                        this_offset_from_min_score = 0
                    else:
                        this_offset_from_min_score = min(abs(x1-c),abs(x2-c))

                    if (not (scaffold1,a,b,c) in category) or \
                       (category[scaffold1,a,b,c]=="break") or \
                       (category[scaffold1,a,b,c]=="gap" and offset[scaffold1,a,b,c]>this_offset_from_min_score):
                        category[scaffold1,a,b,c]="gap"
                        br[scaffold1,a,b,c]=int((x1+x2)/2)
                        offset[scaffold1,a,b,c] = this_offset_from_min_score
                elif x1<a and b<x2:
                    if not (scaffold1,a,b,c) in category:
                        if scores.get((scaffold1,a,b,c),-1e9) < contig_break_threshold:
                            category[scaffold1,a,b,c]="break"
                            br[scaffold1,a,b,c]=c
        categories={}
        for s,a,b,c in br.keys():
            if (category[s,a,b,c]=="break" and scores.get((s,a,b,c),100)< contig_break_threshold) or \
               (category[s,a,b,c]=="gap"   and scores.get((s,a,b,c),100)< unlink_threshold) :

                if clipped_likelihood_for_gaps_only and (score_types.get((s,a,b,c),"fine")=="clipped") and category[s,a,b,c]=="break": 
                    if debug: print("Skip break because we're only using clipped likelihood for un-scaffolding, not contig breaking",s,a,b,c,scores.get((s,a,b,c),-100),sep="\t")
                    if debug: print("X",s,a,b,c,br[s,a,b,c],category[s,a,b,c],score_types.get((s,a,b,c),"fine"),scores.get((s,a,b,c),-100))
                    continue
                
                break_points.append((s,br[s,a,b,c]))
                if True or debug: print("Q",s,a,b,c,br[s,a,b,c],category[s,a,b,c],score_types.get((s,a,b,c),"fine"),scores.get((s,a,b,c),-100))
                categories[s,br[s,a,b,c]]=category[s,a,b,c]
        self.add_breakpoints(break_points,debug=debug)

    def add_breakpoints(self,breaks,debug=False,categories={}):
        """breaks is a list of (scaffold,x) pairs where scaffold is the id of a scaffold in the assembly, 
and x is a basepair offset from the scaffold 5' end of the scaffold where a break should be added.  If x
falls in a gap between contigs, the linkage is broken.  If x falls in a contig, the contig is broken."""
        import networkx as nx

        g = nx.DiGraph()

        scaffold_breaks={}

        for s,x in breaks:
            if not s in scaffold_breaks: scaffold_breaks[s]=[]

        for s,x in breaks:
            scaffold_breaks[s].append(x)

        for s in scaffold_breaks.keys():
             scaffold_breaks[s].sort()

        if debug: print(scaffold_breaks)

        new_layout_lines=[]

        scaffold_broken=False

        break_indices={}

        i=0
        break_counter=0
#        debug=False

        edge_breaks={}
        elength={}
        scaffold={}
        xcoord={}
#
#  Convert the layout lines into a directed graph
#
        contig_length={}
        while i< len(self.layout_lines)-1:
            #print(self.layout_lines[i])
            scaffold1,contig1,base1,end1,x1,l21 = self.layout_lines[  i  ]
            scaffold2,contig2,base2,end2,x2,l22 = self.layout_lines[ i+1 ]
            scaffold[contig1,base1]=scaffold1
            scaffold[contig2,base2]=scaffold2
            i+=1
            if scaffold1==scaffold2: 
                if scaffold1 in scaffold_breaks and (not (contig1,base1)==(contig2,base2)) and len([z for   z in scaffold_breaks[scaffold1] if x1<z and z<x2])>0:  
                    pass # this is an OO join to remove
                else:
                    g.add_edge( (contig1,base1,end1),(contig2,base2,end2) )
                    elength[ (contig1,base1,end1),(contig2,base2,end2) ]=x2-x1
                    if (contig1,base1)==(contig2,base2): contig_length[contig1,base2]=x2-x1

                    if not ( scaffold1==scaffold2 and scaffold1 in scaffold_breaks):
                        continue

                    else:
                        break_positions = [z for   z in scaffold_breaks[scaffold1] if x1<z and z<x2]

                        if len(break_positions)==0: 
                            continue

                        edge_breaks[(contig1,base1,end1),(contig2,base2,end2)]=[]
                        for x in break_positions: # these's a break in here!

                            edge_breaks[(contig1,base1,end1),(contig2,base2,end2)].append( x-x1 )

        for e1,e2 in edge_breaks.keys():
#            print(e1,e2)
            next_nodes= g.neighbors(e2) 
            prev_nodes= g.predecessors(e1) 
#            print(prev_nodes)
            if len(list(next_nodes))>1: 
                raise Exception 
            if len(list(prev_nodes))>1: 
                raise Exception 
            next_node=False
            prev_node=False
            
            if len(list(next_nodes)) == 0:
                continue
            if len(list(previous_nodes)) == 0:
                continue

            if next_nodes:
                next_node = list(next_nodes)[0]
            if prev_nodes:
                prev_node = list(prev_nodes)[0]
            g.remove_edge(e1,e2)

            #print(e1,e2,edge_breaks[e1,e2],scaffold[e1[0],e1[1]],next_node,prev_node)

            if not (e1[0],e1[1])==(e2[0],e2[1]):
                #print("this is just an oo link to udo")
                #this is just an edge to remove
                continue

            ocontig = e1[0]
            
            if e1[2] in ["5",5]: #forward

                
                
                if next_node: g.remove_edge(e2,next_node)
                g.remove_node(e2)

                #                          a        b
                #         prev    e1     x1          x2       e2   next
                #  ---------->     --------------------------->    ----------------->

                #just implement the first and last break.  (this is not necessarily good if the contig is really, really long. XXX TODO )
                edge_breaks[e1,e2].sort()
                x1=edge_breaks[e1,e2][0]
                x2=edge_breaks[e1,e2][-1]

                for i in range(len(edge_breaks[e1,e2])-1):
                    xa = edge_breaks[e1,e2][i  ]
                    xb = edge_breaks[e1,e2][i+1]
                    new_node_a = ( e1[0] , e1[1]+xa, '5' )
                    new_node_b = ( e1[0] , e1[1]+xa, '3' )
                    g.add_edge(new_node_a,new_node_b)
                    elength[ new_node_a,new_node_b ] = xb-xa
                    scaffold[e1[0] , e1[1]+xa] = scaffold[e1[0],e1[1]]

                new_node_1 = (e1[0],e1[1]   ,'3')
                new_node_2 = (e1[0],e1[1]+x2,'5')
                new_node_3 = (e1[0],e1[1]+x2,'3')
                scaffold[e1[0],e1[1]+x2] = scaffold[e1[0],e1[1]]

                if debug: print("b",e1,e2,prev_node,next_node,x1,x2,edge_breaks[e1,e2],new_node_1,new_node_2,   sep="\t")

                g.add_edge(e1,new_node_1)
                g.add_edge(new_node_2,new_node_3)
                if next_node: g.add_edge(new_node_3,next_node)

                LL = elength[e1,e2]
                if next_node: LL2 = elength[e2,next_node]
                elength[ e1,new_node_1 ] = x1
                elength[ new_node_2,new_node_3 ] = LL - x2
                if next_node: elength[ new_node_3,next_node  ] = LL2


                if debug: print("b",e1,e2,prev_node,next_node,x1,x2,edge_breaks[e1,e2],new_node_1,new_node_2, elength.get(( e1,new_node_1 )),elength.get(( new_node_2,new_node_3 )), elength.get((new_node_3,next_node)), sep="\t")
#                print("zzz",contig_length[e1[0],e1[1]],LL,elength[e1,e2],edge_breaks[e1,e2],x1,x2,elength[e1,e2],elength[ e1,new_node_1 ],elength[ new_node_2,new_node_3 ])

            else:                #reverse
                
                if prev_node: g.remove_edge(prev_node,e1)
                g.remove_node(e1)

                #just implement the first and last break.  (this is not necessarily good if the contig is really, really long. XXX TODO )
                edge_breaks[e1,e2].sort()                
                x1=elength[e1,e2]-edge_breaks[e1,e2][-1]
                x2=elength[e1,e2]-edge_breaks[e1,e2][0]


                #                          b        a
                #         prev    e1     x2          x1       e2   next
                #  ---------->    <----------------------------    ----------------->
                #                  A      B           C

                for i in range(len(edge_breaks[e1,e2])-1):
                    xb = elength[e1,e2]-edge_breaks[e1,e2][i  ]
                    xa = elength[e1,e2]-edge_breaks[e1,e2][i+1]
                    new_node_a = ( e1[0] , e2[1]+xa, '5' )
                    new_node_b = ( e1[0] , e2[1]+xa, '3' )
                    g.add_edge(new_node_a,new_node_b)
                    elength[ new_node_a,new_node_b ] = abs(xb-xa)
                    scaffold[e1[0] , e2[1]+xa] = scaffold[e1[0],e2[1]]

                new_node_C = (e1[0],e2[1]   ,'3') #
                new_node_B = (e1[0],e2[1]+x2,'5') #
                new_node_A = (e1[0],e2[1]+x2,'3') #
                scaffold[e1[0],e2[1]+x2] = scaffold[e1[0],e2[1]]

                g.add_edge(new_node_C,e2) #
                g.add_edge(new_node_A,new_node_B) #
                if prev_node: g.add_edge(prev_node,new_node_A) #

                LL =  elength[e1,e2]
                if prev_node: LL2 = elength[prev_node,e1]
                elength[ new_node_C, e2 ] = x1
                elength[ new_node_A,new_node_B ] = LL - x2
                if prev_node: elength[ prev_node, new_node_A  ] = LL2
#                print("zzW",edge_breaks[e1,e2],x1,x2,elength[e1,e2],elength[ new_node_1, e2 ],elength[ new_node_3,new_node_2 ])
            
#        for e1,e2 in elength.keys():
#            if elength[e1,e2]<0:
#                print("wtf?",e1,e2,elength[e1,e2])
#
#  Compute the new layout coordinates and scaffold IDs
#
        scaffold_incr={}
        node_x={}
        contig_length={}
        new_scaffold={}
        for n in sorted(g.nodes()):
            #print(n,scaffold[n[0],n[1]],g.degree(n),g.in_degree(n),g.out_degree(n))
            if g.in_degree(n)==0:
                scaffold_id=scaffold[n[0],n[1]]
                new_scaffold_id = scaffold_id 
                if scaffold_id in scaffold_incr: 
                    new_scaffold_id = scaffold_id + "."+str(scaffold_incr.get(scaffold_id,0))
                new_scaffold[n]=new_scaffold_id
                scaffold_incr[scaffold_id] = scaffold_incr.get(scaffold_id,0)+1
                ll=list(nx.dfs_preorder_nodes(g,n))
                xx=0
                lastl=False
                for l in ll:
                    if lastl:
                        xx+=elength[ lastl,l ]
                        #print("z",l,lastl,elength[lastl,l],node_x[lastl],xx,new_scaffold[n],sep="\t")
                        if (lastl[0],lastl[1])==(l[0],l[1]):
                            contig_length[lastl[0],lastl[1]]=elength[ lastl,l ]
                    #print(new_scaffold_id,l,xx)
                    node_x[l]=xx
                    lastl = l
#
#  Build the new layout_lines list
#
        self.scaffold_lengths={}
        new_layout = []
        for n in g.nodes():
            if g.in_degree(n)==0:

                new_scaffold_id = new_scaffold[n]
                ll=list(nx.dfs_preorder_nodes(g,n))

                for l in ll:
#                    print("w",new_scaffold_id,l[0],l[1],l[2],node_x[l])
#                    print("z",new_scaffold_id,l[0],l[1],l[2],node_x[l],contig_length[l[0],l[1]])
                    new_layout.append( (new_scaffold_id,l[0],l[1],l[2],node_x[l],contig_length[l[0],l[1]] ) )
#                    new_layout.append( (new_scaffold_id,l[0],l[1],l[2],node_x[l],elength[l[0],l[1]] ) )
                    self.scaffold_lengths[new_scaffold_id] = max( node_x[l] , self.scaffold_lengths.get(new_scaffold_id,0))

                    
        
        self.layout_lines=new_layout
#        self.layout_lines.sort( key=lambda x: (x[0],x[4]) )
        self.index_ocontigs()

    def load_assembly(self,infile):
        self.scaffold_lengths={}
        self.bams=[]
        self.layout_lines=[]
        self.model_params={}
        self.ocontig_lengths={}
        self.contig_lengths={}
        self.scaffold_ocontigs={}
        self.bam_objects=False
#        self.model_params={}
        f=open(infile,"r")
        path = ""
        #if infile.startswith("/"):
        path = os.path.dirname(os.path.normpath(infile))
        for l in f:
            #print(l)
            if l[0]=="M":
                c=l.strip().split()
                self.model_params[c[1]] = eval(" ".join(c[2:]))[0]
            if l[0]=="B":
                #self.bams.append(l[2:].strip())
                possible_path = l[2:].strip()
                if possible_path.startswith("/"):
                   self.bams.append(l[2:].strip())
                else:
                    self.bams.append(os.path.join(path,l[2:].strip()))   
            if l[0]=="S":
                #self.shotgun.append(l[2:].strip())
                possible_path = l[2:].strip()
                if possible_path.startswith("/"):
                    self.shotgun.append(l[2:].strip())
                else:
                    self.shotgun.append(os.path.join(path,l[2:].strip()))
            if l[0]=="P": # Layout
                self.layout_lines.append( parse_pline2(l,self.scaffold_lengths) )
                scaffold,contig,base,end,x,lll = self.layout_lines[-1]
                self.contig_lengths[contig+"_"+str(base)]=lll
                self.contig_lengths[contig,base]=lll
                self.scaffold_ocontigs[scaffold] = self.scaffold_ocontigs.get(scaffold,[])+[contig]
            if l[0]=="D": # Masked regions (from D for shotgun depth)
                dummy,ocontig,a,b = l.strip().split()
                if not ocontig in self.masked_regions: self.masked_regions[ocontig]=[]
                self.masked_regions[ocontig].append((int(a),int(b)))
            if l[0]=="L": # Length of input contigs, same as bam header lengths  
                dummy,ocontig,l = l.strip().split()
                #if not ocontig in self.masked_regions: self.masked_regions[ocontig]=[]
                self.ocontig_lengths[ocontig] = int(l)
        f.close()
        self.index_ocontigs()
#        print("setup mapper")
        self.setup_mapper(debug=True)

    def save_assembly(self,outfile):
        f=open(outfile,"w")
        for b in self.shotgun:
            f.write("S {}\n".format(b))
        for b in self.bams:
            f.write("B {}\n".format(b))
        #f.write("M {}\n".format(self.model_params))
        if self.model_params:
            for k in sorted(self.model_params.keys()):
                f.write("M {} {}\n".format(k,(self.model_params[k],)))
        for l in self.layout_lines:
            f.write("P "+" ".join(map(str,l))+"\n")
        for oc in sorted(self.masked_regions.keys()):
            for a,b in self.masked_regions[oc]:
                f.write("D {} {} {}\n".format(oc,a,b))
        for oc in sorted(self.ocontig_lengths.keys()):
            f.write("L {} {}\n".format(oc,self.ocontig_lengths[oc]))

#        self.masked_regions={}
        f.close()

    def scaffold_mapping():
        for i in range(len(self.layout_lines)):
            scaffold1,contig1,base1,end1,x1,l11 = self.layout_lines[i  ]
            if end1 == "3":
                yield("{}_{}".format(contig1,base1),scaffold1)
        
    def as_edges(self):
        for i in range(len(self.layout_lines)-1):
            scaffold1,contig1,base1,end1,x1,l21 = self.layout_lines[i  ]
            scaffold2,contig2,base2,end2,x2,l22 = self.layout_lines[i+1]
            if scaffold1==scaffold2 :
                if contig1==contig2 and base1==base2: 
                    yield(["#edge:","{}_{}.{}".format(contig1,base1,end1) ,"{}_{}.{}".format(contig2,base2,end2) ,{'length':l21,'contig':True}])
                else:
                    yield(["#edge:","{}_{}.{}".format(contig1,base1,end1) ,"{}_{}.{}".format(contig2,base2,end2) ,{'length':x2-x1,'contig':False}])
#                    yield(["#edge:",contig1+"."+end1,contig2+"."+end2,{'length':x2-x1,'contig':False}])

    def contig_coords(self):
        seen = {}
        for l in self.layout_lines:
            scaffold1,contig1,base1,end1,x1,l21 = l
            if not (contig1,base1) in seen:
                yield( {"contig":contig1,"base":base1,"scaffold":scaffold1,"span":(x1,x1+l21),"flipped":end1 in [3,"3"] })

            seen[contig1,base1]=True

    def dump_layout(self,f):
#        f=open(outfile,"w")
        for s in self.scaffold_lengths.keys():
            f.write("{} slen\n".format(self.scaffold_lengths[s]))
        for l in self.layout_lines:
            scaffold1,contig1,base1,end1,x1,l21 = l
            f.write("p: {} {}_{}.{} {} - -1 {} {} {}\n".format(scaffold1,contig1,base1,end1,x1,self.scaffold_lengths[scaffold1],l21,False))
#        f.close()
        #p: 1 Scaffold105707_1.3 3019 - -1 85132546 3048 False

    def load_playout(self,layoutfile):
        #print(layoutfile)
        self.scaffold_lengths={}
#        self.contig_lengths={}
        if not os.path.exists(layoutfile):
            raise Exception 
        self.layout_lines = [ parse_pline(c,self.scaffold_lengths) for c in open(layoutfile,"rt") if c[:2]=="p:" ]
        #for l in self.layout_lines:
            #print(l)
        self.index_ocontigs()

    def check_bams(self):
        for bamfile in self.bams:
            if not os.path.exists(bamfile):
                raise Exception 

    def set_modelparams(self,datamodel):
        if os.path.exists(datamodel):
            self.model_params = eval(open(datamodel).read())
        else:
            self.model_params=datamodel

if __name__=="__main__":

    import sys
    import argparse

    parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument('-d','--debug',default=False,action="store_true",help="Turn on debugging ouput")
    parser.add_argument('-E','--edges',default=False,action="store_true",help="Print the layout as edges.")
#    parser.add_argument('-E','--edges',default=False,action="store_true",help="Print the layout as edges.")
    parser.add_argument('-D','--dump_layout',default=False,action="store_true",help="Print the layout as edges.")
    #parser.add_argument('-p','--progress',default=False,action="store_true",help="Print progress info")
    parser.add_argument('-M','--datamodel',default=False, help="Name of a file containing the model parameter dictionary")
    parser.add_argument('-m','--mask'   ,default=[], action="append",help="Name of a file containing regions of the input assembly to mask out")
    parser.add_argument('-b','--bamfile',default=[], action="append",help="A chicago bam file")
    parser.add_argument('-S','--shotgunbamfile',default=[], action="append",help="A shotgun bam file")
    #parser.add_argument('-p','--progress',default=False,action="store_true",help="Print progress info")
    parser.add_argument('-L','--layout',default=False,help="File containing a scaffolding layout in 'p:' lines.")
    parser.add_argument('-o','--outfile',default=False,help="Filename for serialization of the assembly file.")
    parser.add_argument('-i','--infile',default=False,help="Filename for serialised assembly input file.")

    args = parser.parse_args()
    
#    print(args)
    if args.infile:
        asf = HiriseAssembly()
        asf.load_assembly(args.infile)
#        asf.validate()
    else:
        asf = HiriseAssembly( {"bams":args.bamfile, "datamodel": args.datamodel, "layout": args.layout , "debug":args.debug, "shotgun": args.shotgunbamfile} )
    if args.debug: print("# done with intital load")
    asf.validate()
#    asf.index_ocontigs()

    if args.datamodel:
        asf.set_modelparams(args.datamodel)

    for segments_file in args.mask:
        asf.add_mask_regions(filename=segments_file)
        asf.merge_masked_regions()

    if args.outfile:
        asf.save_assembly( args.outfile )
        
    if args.edges:
        for e in asf.as_edges():
            print("\t".join(map(str,e)))

#    asf.add_breakpoints( [( '1', 1692 ),('1',10000)] )
#    asf.add_breakpoints( [( '1', 1692 ),('1',2000),('1',10000)] )
#    asf.add_breakpoints( [( '1', 1692 ),('1',2000),('1',2100),('1',10000)] )

#    asf.validate()

#    asf.load_ocontig_info()
#    print("#loaded ocontig info")
#    asf.load_ocontig_info()
#    print("#loaded ocontig info")
#    print(asf.get_links(["Scaffold116325"],links={}))
#    print(asf.get_links(["Scaffold116325"],links={}))
#    asf.get_scaffold_links("1000",skipI=True)

#1691
    if args.dump_layout:
        asf.dump_layout(sys.stdout)

#!/usr/bin/env python3

#
# Copyright 2015 Dovetail Genomics LLC
#
#

import os
import re
import pysam
from bamtags import BamTags
import gc
import random
import bisect
import numpy as np

def parse_pline(l,scaffold_lengths):
    c=l.strip().split()
    m=re.match("(.*)_(\d+)\.([35])",c[2])
    scaffold_lengths[c[1]]=int(c[6])
    return( c[1],m.group(1),int(m.group(2)),m.group(3),int(c[3]),int(c[7]) )

def parse_pline2(l,scaffold_lengths):
    c=l.strip().split()
    base=int(c[3])
    x=int(c[5])
    scaffold_lengths[c[1]]=max(scaffold_lengths.get(c[1],0),x)
    return( c[1],c[2],base,c[4],x,int(c[6]) )

def aln2bins(aln,binsize):        
    return( ( aln.tid,int(aln.pos/binsize) ),( aln.rnext,int(aln.pnext/binsize)))


class GrowingNPArray:
    def __init__(self,dtype):
        self.data = np.zeros((10000,),dtype=dtype)
        self.N=10000
        self.n=0
        self.dtype = dtype

    def append(self,d):
        if self.n==self.N:
            self.N *=2
            dd = np.zeros((self.N,),dtype=self.dtype)
            dd[:self.n] = self.data
            self.data = dd
        self.data[self.n]=d
        self.n+=1

    def finalize(self):
        dd = np.zeros((self.n,),dtype=self.dtype)
        dd[:self.n] = self.data[:self.n]
        return dd


def pairs_overlap(x,y):
    a=min(x[0],x[1])
    b=max(x[0],x[1])
    c=min(y[0],y[1])
    d=max(y[0],y[1])

    if a<=c and c<=b: return True
    if a<=d and d<=b: return True
    if c<=a and a<=d: return True
    if c<=b and b<=d: return True
    return False
                    
class HiriseAssembly:
    """This class is designed to encapsulate a Hirise Assembly."""

    def __init__(self,options={}):
        """Options currently used are:  bams, bam, datamodel, layout"""
        self.bams=[]
        if options.get("bams"): self.bams+=list(options["bams"])
        if options.get("bam") : self.bams +=[options["bam"]]
        self.model_params=False

        self.binsize=10000
        self.bam_objects=False
        self.shotgun_bam_objects=False

        if "bam" in options or "bams" in options:
            self.load_ocontig_info()
        
        if "datamodel" in options and options["datamodel"]: 
            self.set_modelparams(options["datamodel"])

        self.layout_lines=[]
        if options.get("layout"):
            self.load_playout(options["layout"])
        self.masked_regions={}
        self.check_bams()
        self.ocontig_fa=False
        self.ocontig_fasta=False
        if "shotgun" in options: 
            self.shotgun = options["shotgun"]
        else:
            self.shotgun = []
        self.contigx=False

    def get_seq(self,contig):
        if not self.ocontig_fasta:
            raise Exception 
        
        if not os.path.exists(self.ocontig_fasta+".fai"):
            pysam.faidx(self.ocontig_fasta)
        if not self.ocontig_fa:
            self.ocontig_fa = pysam.Fastafile(self.ocontig_fasta)

        m=re.match("(.*)_(\d+)",contig)
        ocontig=m.group(1)
        x=int(m.group(2))
        y=x+self.contig_length(contig)
        leng=self.ocontig_lengths[ocontig]
        #return self.ocontig_lengths[m.group(1)]


        #dummy,scaffold,x,y,leng = l.strip().split()
        #x=int(x)
        #y=int(y)
        #leng=int(leng)
#        y=min(y,leng)
        seq = self.ocontig_fa.fetch(ocontig,x-1,y)
        return(seq)

    def contig_length(self,contig):
        if contig in self.contig_lengths:
            return self.contig_lengths[contig]
        else:
            m=re.match("(.*)_(\d+)",contig)
            return self.ocontig_lengths[m.group(1)]
        return 0

    def contig_end_to_region(self,contig_end_string,w):
        ((this_contig,this_base,this_end),) = re.findall("^(.*)_(\d+)\.([53])$",contig_end_string)
        print(this_contig,this_base,this_end)
        this_base = int(this_base)
#        this_end = int(this_end)
        for scaffold,contig,base,end,x,l in self.layout_lines:
            slen = self.scaffold_lengths[scaffold]
            if this_contig == contig and this_base == base and this_end == end:
                return(scaffold,max(0,x-w),min(slen,x+w),(x<slen-x))
        print("so such end found",contig_end_string)
        raise Exception

    def index_ocontigs(self):
        self.ocontig_bases  = {}
        self.contig_lengths = {}
        self.contig_scaffold= {}
        self.load_ocontig_info()
#        self.scaffold_ocontigs={}
        
        ocb=self.ocontig_bases
        for scaffold,contig,base,end,x,l in self.layout_lines:

#            if not scaffold in self.scaffold_ocontigs:
#                self.scaffold_ocontigs[scaffold]={}
#            self.scaffold_ocontigs[scaffold][self.contig_ocontig[contig]]=1

            self.contig_scaffold[contig,base]=scaffold
            if end in [5,"5"] and not base==1:
                ocb[contig] = ocb.get(contig,[])+[ base ]
            if end in [5,"5"]:  
                self.contig_lengths[ contig+"_"+str(base) ]=l
                self.contig_lengths[ contig,base ]=l
        for b in ocb.keys():
            ocb[b].sort()
            #print("ocb:",b,ocb[b])
        self.total_ocontig_length = sum(self.ocontig_lengths.values())
        self.rc_buffer=[]
        x=0
        for i in self.ocontigs_iter():
            self.rc_buffer.append((x,i))
            x+=self.ocontig_lengths[i]

    def broken_coord(self,ocontig,x,tuples=False):
        if tuples:
            if not ocontig in self.ocontig_bases:
                if (not (ocontig,1) in self.contig_lengths ) or x < self.contig_lengths[(ocontig,1)]: 
#                    print(ocontig,"1")
                    return((ocontig,1), x)
                else:
#                    print(ocontig,"2")
                    return((ocontig, self.contig_lengths[(ocontig,1)]), x)
            for b in self.ocontig_bases[ocontig][::-1]:
                if b<x:
                    new_contig = (ocontig,b)
                    #print("zxzx",ocontig,x,b,new_contig,self.contig_lengths.get(new_contig,False))
                    if (not new_contig in self.contig_lengths) or (x-b)<self.contig_lengths[new_contig]:
#                        print(ocontig,"3")
                        return((ocontig,b), x-b)
                    else:
                        return (False,False)
            new_contig = (ocontig,1)
            if (not new_contig in self.contig_lengths) or x < self.contig_lengths[new_contig]:
#                print(ocontig,"4",new_contig,new_contig in self.contig_lengths)
                return(new_contig, x)
            else:
                return (False,False)
            

        else:
            if not ocontig in self.ocontig_bases:
                if (not ocontig+"_1" in self.contig_lengths ) or x < self.contig_lengths[ocontig+"_1"]: 
                    return(ocontig+"_1", x)
                else:
                    return(ocontig+"_"+str( self.contig_lengths[ocontig+"_1"] ), x)
            for b in self.ocontig_bases[ocontig][::-1]:
                if b<x:
                    new_contig = ocontig+"_"+str(b)
                    #print("zxzx",ocontig,x,b,new_contig,self.contig_lengths.get(new_contig,False))
                    if (not new_contig in self.contig_lengths) or (x-b)<self.contig_lengths[new_contig]:
                        return(ocontig+"_"+str(b), x-b)
                    else:
                        return (False,False)
            new_contig = ocontig+"_1"
            if (not new_contig in self.contig_lengths) or x < self.contig_lengths[new_contig]:
                return(new_contig, x)
            else:
                return (False,False)

    def merge_masked_regions(self,debug=False):
        mr = self.masked_regions
        for ocontig in mr.keys():
            new_regions=[]
            edges = []
            for x,y in mr[ocontig]:
                edges.append((x,1))
                edges.append((y,-1))
            edges.sort(key=lambda x:(x[0],-x[1]))
            rs=0
            state=0
            for x,b in edges:
                rs+=b
                if state==0 and rs>0:
                    start=x
                    state=1
                elif state==1 and rs==0:
                    end=x
                    new_regions.append((start,end))
                    state=0
            if debug and not (len(mr[ocontig])==len(new_regions)):
                print("#mmr",ocontig,len(mr[ocontig]),len(new_regions),mr[ocontig],new_regions)
            mr[ocontig]=new_regions

    def add_mask_regions(self,segments=[],filename=False):
        masked_regions = self.masked_regions
        if filename:
            for line in open(filename):
                c=line.strip().split()
                ocontig,x,y = c[0],int(c[1]),int(c[2])
                if not ocontig in masked_regions: masked_regions[ocontig]=[]
                masked_regions[ocontig].append((x,y))
        for ocontig,x,y in segments:
            if not ocontig in masked_regions: masked_regions[ocontig]=[]
            masked_regions[ocontig].append((x,y))
            
        for oc in masked_regions.keys():
            masked_regions[oc].sort()

    def setup_mapper(self,debug=False):
#        debug=True
        self.contigx={}
        self.contigst={}
        cex={}
        for scaffold,contig,base,end,x,l in self.layout_lines:
            if end in [5,"5"]:
                self.contigx[contig,base] = x
            cex[contig,base,end]=x
        for contig,base in self.contigx.keys():
            if cex[contig,base,"5"] < cex[contig,base,"3"]: 
                self.contigst[contig,base]=1
            else: 
                self.contigst[contig,base]=-1
#        if debug:
#            for contig,base in self.contigx.keys():
#                print("m:",contig,base,self.contigx[contig,base],self.contigst[contig,base])


    def scaffolds_iter(self,debug=False):
        if not self.layout_lines or len(self.layout_lines)==0:
            self.make_trivial_layout_lines(debug=debug)
            self.validate()

        last_scaffold=False
        for scaffold,contig,base,end,x,l in self.layout_lines:
            if (not last_scaffold) or not scaffold==last_scaffold:
                last_scaffold=scaffold
#                print("#",scaffold)
                yield(scaffold)

    def scaffold_coord(self,c,x,debug=False):
        if not self.contigx: self.setup_mapper(debug=debug)
        y=self.contigx[c]+self.contigst[c]*x
        if debug: print("sc:",c,x,y)
        return y

    def chicago_support_scan(self,scaffold,minsupport=0.0,mapq=10,debug=False,minx=0,maxx=1e6,gpf=False,logfile=False):

        import chicago_edge_scores as ces
        ces.set_exp_insert_size_dist_fit_params(self.model_params)
        model=ces.model

        links={}
        self.get_scaffold_links(scaffold,skipI=False,mapq=mapq,links=links,debug=debug,tuples=True)
#        model = ces.
        pairs=[]
        edges=[]
        buff=100
        nb=[]
        joins=[]
        t1=50.0
        t2=minsupport
        slen=self.ocontig_lengths[scaffold]
        for c1,c2 in links.keys():
            if debug: print(c1,c2,links[c1,c2])
            for x,y in links[c1,c2]:
                pairs.append( ( self.scaffold_coord(c1,x,debug=debug),self.scaffold_coord(c2,y,debug=debug)) )
        for a,b in pairs:               
            a,b=min(a,b),max(a,b)
            if debug:
                print("z",a,b)


            if (b-a)>minx and (b-a)<maxx:
                if gpf: gpf.write("{}\t{}\t{}\t{}\n".format(0.5*(a+b),0.5*(b-a),w,z))
                ll=model.lnF(b-a)

                edges.append( tuple([a+buff         , ll,c1,"a"]) )
                edges.append( tuple([b-buff         ,-ll,c2,"b"]) )

                nb.append( tuple([a+buff         , 1]) )
                nb.append( tuple([b-buff         ,-1]) )

        if gpf: gpf.write("\n\n\n")
        edges.sort()
        nb.sort()
        rs=0
        n=0
        tripped=False
        last=0
        state=0
        stretch_start=0
        low_point=0.0
        ji=0
        gap_scores={}
        gap_coverages={}
        max_t1=0
        min_t1=slen
        breakers=[]
        
        for i in range(len(edges)):
            rs+=edges[i][1]
            n+=nb[i][1]
            x=edges[i][0]

            try:
                score=model.cutScore(slen,x,n,rs,rangeCutoff=maxx,minSep=minx)
            except Exception as e:
                print(edges)
                print(i,edges[i],nb[i],rs,n,x,scaffold,slen,len(edges))
                raise e


            if score>t1:
                tripped=True
                last_t1=x
                min_t1 = min(min_t1,x)
                max_t1 = max(max_t1,x)
            if tripped and score<t2 and state==0:
                stretch_start=x
                state=1
                low_point =score
                minx = x
            if state==1 and score>t2:
                breakers.append((scaffold,stretch_start,x,minx,low_point,slen))
#                if logfile: logfile.write("{} {} {} {} {} {}\n".format(scaffold,stretch_start,x,minx,low_point,slen))
                state=0
            if state==1:
                if score<low_point:
                    low_point=score
                    minx=x
            if debug: print("dd:",scaffold,edges[i][0],score,rs,state,x,stretch_start,score,n,edges[i][2])
            if gpf: gpf.write("{}\t{}\t{}\t{}\n".format(edges[i][0],score,rs,n))

            last=edges[i][0]
            
        for scaffold,stretch_start,x,minx,low_point,slen in breakers:
            if stretch_start < last_t1:
#                if logfile: logfile.write("{} {} {} {} {} {} {} {}\n".format(scaffold,stretch_start,x,minx,low_point,slen,min_t1,max_t1))
                if logfile: logfile.write("{} {} {} {} {} {}\n".format(scaffold,stretch_start,x,minx,low_point,slen))

    def make_trivial_layout_lines(self,debug=False):
        self.layout_lines=[]
        #print(self.ocontig_lengths)
        if len(self.ocontig_lengths)==0:
            print("#load_contig_info")
            self.load_ocontig_info()
        i=1
        for ocontig in self.ocontigs_iter():
            scaffold=ocontig
#            if debug: print(scaffold,ocontig,1,"5",0,self.ocontig_lengths[ocontig],sep="\t")
            self.layout_lines.append((scaffold,ocontig,1,"5",0,self.ocontig_lengths[ocontig])) #
            self.layout_lines.append((scaffold,ocontig,1,"3",self.ocontig_lengths[ocontig],self.ocontig_lengths[ocontig])) #
            self.scaffold_lengths[scaffold] = self.ocontig_lengths[ocontig]
            self.scaffold_ocontigs[scaffold] = self.scaffold_ocontigs.get(scaffold,[])+[ocontig]
            i+=1
        self.index_ocontigs()
        self.setup_mapper()

    def get_scaffold_links(self,Tscaffold,skipI=False,mapq=10,links={},debug=False,tuples=False):
        if not self.layout_lines or len(self.layout_lines)==0:
            self.make_trivial_layout_lines(debug=debug)
            self.validate()
        contigs={}
        ocontigs={}
        if not tuples:
            for scaffold,contig,base,end,x,l in self.layout_lines:
                if scaffold==Tscaffold and end in [5,"5"]:
                    contigs[ contig+"_"+str(base) ]=1  
                    ocontigs[contig]=1
        else:
            for scaffold,contig,base,end,x,l in self.layout_lines:
                if scaffold==Tscaffold and end in [5,"5"]:
                    contigs[ contig,base ]=1  
                    ocontigs[contig]=1
                
        if debug: print(ocontigs,contigs)
        self.get_links(list(ocontigs.keys()),skipI,mapq,links,contigs,tuples=tuples)

    def segments_iter(self):
        for l in self.layout_lines:
            scaffold1,contig1,base1,end1,x1,l21 =l
            if end1 in [5 , "5"]:
                #cc=self.contig_length(contig1+"_"+str(base1))
                yield( contig1,base1,base1+l21,self.ocontig_lengths[contig1] )

    def contigs_iter(self):
        for l in self.layout_lines:
            scaffold1,contig1,base1,end1,x1,l21 =l
            if end1 in [5 , "5"]:
                yield( contig1 +"_"+str(base1) )

    def ocontigs_iter(self):
        for ocontig in self.ocontig_lengths.keys():
            yield(ocontig)

    def load_ocontig_info(self):
        import pysam

        self.ocontig_lengths={}
        if not self.bam_objects:
            self.bam_objects={}
            for b in self.bams:
                self.bam_objects[b] =  pysam.Samfile(b,"rb") 
        for bam in self.bam_objects.values():
            h=bam.header
            seqs=h['SQ']
            for s in seqs:
                if s['SN'] in self.ocontig_lengths:
                    if not  self.ocontig_lengths[s['SN']]==s['LN'] :
                        raise Exception 
                else:
                    self.ocontig_lengths[s['SN']]=s['LN']

    def random_window(self,wlen=1000):
        r=int(random.random()*self.total_ocontig_length)
        x=bisect.bisect_right(self.rc_buffer,(r,"a"))
        while r+x > self.total_ocontig_length or (x<len(self.rc_buffer)-1 and r+wlen>self.rc_buffer[x+1][0]):
            r=int(random.random()*self.total_ocontig_length)
            x=bisect.bisect_right(self.rc_buffer,(r,"a"))
#            print( r,self.rc_buffer[x],self.rc_buffer[x+1],r-self.rc_buffer[x][0])
        x-=1
#        print( r,self.rc_buffer[x],self.rc_buffer[x+1],r-self.rc_buffer[x][0])
        return( self.rc_buffer[x][1], r-self.rc_buffer[x][0], r-self.rc_buffer[x][0]+wlen )

    

    def top_other_contigs(self,reference,mapq=10,n=2,debug=False,bins=False,binx=0):
        nother={}
        if not self.bam_objects:
            self.bam_objects={}
            for b in self.bams:
                self.bam_objects[b] = pysam.Samfile(b,"rb")
        if debug: print(self.bam_objects)

        if bins:
            region="{}:{}-{}".format(reference,binx*self.binsize,(binx+1)*self.binsize)
            for bam in self.bam_objects.values():
                for aln in bam.fetch(region=region,until_eof=True):
                    if aln.is_duplicate :              continue
                    if aln.mapq < mapq               : continue
                    if BamTags.mate_mapq(aln) < mapq : continue
                    my_bin,rbin = aln2bins(aln,self.binsize)
                    if not my_bin == rbin : 
                        nother[rbin] = nother.get(rbin,0)+1
        else:
            for bam in self.bam_objects.values():
                for aln in bam.fetch(reference=reference,until_eof=True):
                    if aln.is_duplicate :              continue
                    if aln.mapq < mapq               : continue
                    if BamTags.mate_mapq(aln) < mapq : continue
                    if not aln.rnext == aln.tid : 
                        nother[aln.rnext] = nother.get(aln.rnext,0)+1
 
        kk=list(nother.keys())
        kk.sort(key=lambda x: nother[x],reverse=True)
        if debug: print(kk[:2])
        return(tuple(kk[:2]))

    def window_stats(self,reference,xa,xb,mapq=10,debug=False,bins=False):
        if not self.bam_objects:
            self.bam_objects=[]
            for b in self.bams:
                self.bam_objects.append( pysam.Samfile(b,"rb") )

        if not self.shotgun_bam_objects:
            self.shotgun_bam_objects=[]
            for b in self.shotgun:
                self.shotgun_bam_objects.append( pysam.Samfile(b,"rb") )
        
        binx = 0
        if bins:
            binx = int(xa/self.binsize)

        top2 = self.top_other_contigs(reference,mapq=mapq,debug=debug,bins=bins,binx=binx)
        region = "{}:{}-{}".format(reference,max(1,xa-1000),xb+1000)
        n_aligned=0
        n_hiq=0
        templates={}
        n_lowq   =0
        occ ={}
        occ2={}
        for bam in self.bam_objects.values():
            if debug: print("#",region,xa,xb,bam)
            for aln in bam.fetch(region=region):
                if aln.pos > xb:                   continue
                if aln.pos < xa:                   continue
                if aln.is_duplicate :              continue
                n_aligned+=1
                if aln.mapq < mapq               : continue
                if BamTags.mate_mapq(aln) < mapq : continue
                templates[aln.query_name]=templates.get(aln.query_name,0)+1
                if bins:
                    my_bin,rbin = aln2bins(aln,self.binsize)
                    if not my_bin==rbin:
                        occ[rbin] = occ.get(rbin,0)+1
                        if not rbin in top2:
                            occ2[rbin] = occ2.get(rbin,0)+1
                            
                else:
                    if (not aln.rnext == aln.tid):
                        occ[aln.rnext] = occ.get(aln.rnext,0)+1
                        if (not aln.rnext in top2) : 
                            occ2[aln.rnext] = occ2.get(aln.rnext,0)+1
                n_hiq+=1
            if debug: print("#",n_hiq,n_aligned)
        n_both   = len([ i for i in templates.keys() if templates[i]>1  ])
        n_other  = len(occ.keys())
        n_other2 = len([ i for i in occ.keys() if occ[i]>1  ])
        n_other3 = len([ i for i in occ.keys() if occ[i]>2  ])

        n_otherB  = len(occ2.keys())
        n_otherB2 = len([ i for i in occ2.keys() if occ[i]>1  ])
        n_otherB3 = len([ i for i in occ2.keys() if occ[i]>2  ])

        n_shotgun=0
        for bam in self.shotgun_bam_objects:
            for aln in bam.fetch(region=region):
                if aln.pos > xb: continue
                if aln.is_duplicate :              continue
                if aln.mapq < mapq               : continue
                n_shotgun+=1
        
        print("\t".join(map(str,[region,n_aligned,n_hiq,n_both,n_other,n_other2,n_other3,n_otherB,n_otherB2,n_otherB3,n_shotgun])))

    def chicago_promiscuity_mask(self,ocontig,mapq=10,w=1000,minlinks=2,maxothers=3,outfile=False,bins=False,debug=False):
        if not self.bam_objects:
            self.bam_objects=[]
            for b in self.bams:
                self.bam_objects.append( pysam.Samfile(b,"rb") )
#        dtype=[('x',int),('t',int),('y',int),('b',int),('r',int)]
        
        dtype=[('x',int),('t',int),('y',int)]
        buffer = GrowingNPArray(dtype=dtype)

        i=0
        tid=0
        for bam in self.bam_objects:                
            if debug: print("#",ocontig,bam)
            tid = bam.gettid(ocontig)
            for aln in bam.fetch(reference=ocontig):
                if aln.rnext<0: continue
                if aln.mapq < mapq               : continue
                if BamTags.mate_mapq(aln) < mapq : continue
                if aln.is_duplicate : continue
#                buffer.append((aln.pos,aln.rnext,aln.pnext,i,aln.tid))
                buffer.append((aln.pos,aln.rnext,aln.pnext))
            i+=1

        buffer = buffer.finalize()
        buffer = np.sort(buffer,order='x',kind="mergesort")

        bad_edges=[]


        def get_contig(x,b):
            return x[1]

        def get_contig_bin(x,b):
            return(x[1]+(int(x[2]/b)/1.0e6))

        if bins:
            get_bin = get_contig_bin
        else:
            get_bin = get_contig

        bc={}
        i=0
        j=0
        N=len(buffer)
        while i < N:
            mybin = get_bin((0,tid,buffer[i][0]),self.binsize)
            while j<N and buffer[j][0]<buffer[i][0]+w: 
                jbin = get_bin(buffer[j],self.binsize)
                bc[ jbin ] = bc.get(jbin,0)+1
                j+=1
                
            n_in_window = j-i
            
            nw=len([ ii for ii in bc.keys() if bc[ii]>=minlinks and not ii==mybin ])  #n_other_windows(buffer[i:j],minlinks,bins,self.binsize)

#            if debug: print(tid,buffer[i][0],n_in_window,nw,mybin,[ (ii,bc[ii]) for ii in bc.keys() if bc[ii]>=minlinks and not ii==mybin ],[ get_bin(x,self.binsize) for x in buffer[i:j] ])
            if debug: print(tid,buffer[i][0],n_in_window,nw,mybin,[ (ii,bc[ii]) for ii in bc.keys() if bc[ii]>=minlinks and not ii==mybin ])


            if nw > maxothers:
                bad_edges.append((buffer[i  ][0], 1))
                bad_edges.append((buffer[j-1][0],-1))

            ibin = get_bin(buffer[i],self.binsize)
            bc[ ibin ] = bc[ibin]-1
            if bc[ibin]==0: bc.pop(ibin)
            i+=1

        bad_edges.sort()
        if debug:
            for be in bad_edges:
                print(e)

        rs=0
        state=0
        for x,dy in bad_edges:
            rs+=dy
            if state==0 and rs>0:
                startx=x
                state=1
            elif state==1 and rs==0:
                outfile.write("{} {} {}\n".format(ocontig,startx,x))
                if debug: print("#{} {} {}\n".format(ocontig,startx,x))
                state=0
        outfile.flush()

        
    def chicago_read_density_mask(self,ocontigs,mapq=10,outfile=False,w=1000,cutoff=100,debug=False,shotgunT=False):

        if not self.bam_objects:
            self.bam_objects=[]
            for b in self.bams:
                self.bam_objects.append( pysam.Samfile(b,"rb") )

        if not self.shotgun_bam_objects:
            self.shotgun_bam_objects=[]
            for b in self.shotgun:
                self.shotgun_bam_objects.append( pysam.Samfile(b,"rb") )

        for ocontig in ocontigs:
            locations=[]
            shotgun_locations=[]
#            gc.disable()
            for bam in self.bam_objects:                
                if debug: print("#",ocontig,bam)
                for aln in bam.fetch(reference=ocontig):
#                    if debug: print("#",aln.pos,aln.is_duplicate,aln.mapq,BamTags.mate_mapq(aln),aln.query_name,len(locations))
                    if aln.rnext<0: continue
                    if aln.mapq < mapq               : continue
                    if BamTags.mate_mapq(aln) < mapq : continue
                    if aln.is_duplicate : continue
                    locations.append(aln.pos)
                if debug: print("#",len(locations))
            if shotgunT:
                for bam in self.shotgun_bam_objects:                
                    if debug: print("#",ocontig,bam)
                    for aln in bam.fetch(reference=ocontig):
                        if aln.mapq < mapq  : continue
                        if aln.is_duplicate : continue
                        shotgun_locations.append(aln.pos)
                    if debug: print("#",len(locations),len(shotgun_locations))

#            gc.enable()
            locations.sort()
            shotgun_locations.sort()

            i=0
            j=0
            bad_edges=[]

            while i<len(locations):
                while j<len(locations) and locations[j]<locations[i]+w: j+=1
                n_in_window = j-i
                if debug:
                    print("chicago_density",ocontig,locations[i],locations[j-1],j-i,locations[j-1]-locations[i],locations[i:j])
                if n_in_window > cutoff:
    #                outfile.write("{} {} {}\n".format(ocontig,locations[i],locations[j-1]))
                    bad_edges.append((locations[i],1)  )
                    bad_edges.append((locations[j-1],-1)  )
                i+=1

            if shotgunT:
                i=0
                j=0
                while i<len(shotgun_locations):
                    while j<len(shotgun_locations) and shotgun_locations[j]<shotgun_locations[i]+w: j+=1
                    n_in_window = j-i
                    if debug:
                        print("shotgun_density",ocontig,shotgun_locations[i],shotgun_locations[j-1],j-i,shotgun_locations[j-1]-shotgun_locations[i])
                    if n_in_window > shotgunT:
        #                outfile.write("{} {} {}\n".format(ocontig,locations[i],locations[j-1]))
                        bad_edges.append((shotgun_locations[i],1)  )
                        bad_edges.append((shotgun_locations[j-1],-1)  )
                    i+=1


            bad_edges.sort()

            rs=0
            state=0
            for x,dy in bad_edges:
                rs+=dy
                if state==0 and rs>0:
                    startx=x
                    state=1
                elif state==1 and rs==0:
                    outfile.write("{} {} {}\n".format(ocontig,startx,x))
                    state=0
        outfile.flush()


    def make_scaffold_mask(self,scaffold): #,mapper,mask):
        s_mask = {}
        segments={}
        slen = self.scaffold_lengths[scaffold]

        for ocontig in self.scaffold_ocontigs[scaffold]:
#            print(ocontig)
            for a,b in self.masked_regions.get(ocontig,[]):
                
                c1,x = self.broken_coord(ocontig,a,tuples=True)
                c2,y = self.broken_coord(ocontig,b,tuples=True)

                s1=self.contig_scaffold.get(c1)
                s2=self.contig_scaffold.get(c2)

                if c1==c2 and s1==s2 and s1==scaffold:

                    xx=self.scaffold_coord(c1,x)
                    yy=self.scaffold_coord(c2,y)
                    #print("#msm:",scaffold,ocontig,a,b,x,y,xx,yy,s1,s2,scaffold,sep="\t")

                    segments[ max(0,min(xx,yy)) , min( slen, max(xx,yy) )   ]=1
                elif (not c1==c2) and s1==s2 and s1==scaffold:
                    
                    xA = x
                    yA = self.contig_lengths[c1]
                    xx=self.scaffold_coord(c1,xA)
                    yy=self.scaffold_coord(c1,yA)
                    segments[ max(0,min(xx,yy)) , min( slen, max(xx,yy) )   ]=1

                    xB = 1
                    yB = y
                    xx=self.scaffold_coord(c2,xB)
                    yy=self.scaffold_coord(c2,yB)
                    segments[ max(0,min(xx,yy)) , min( slen, max(xx,yy) )   ]=1

                elif (not c1==c2) and s1==scaffold and not s2==scaffold:
                    
                    xA = x
                    yA = self.contig_lengths[c1]
                    xx=self.scaffold_coord(c1,xA)
                    yy=self.scaffold_coord(c1,yA)
                    segments[ max(0,min(xx,yy)) , min( slen, max(xx,yy) )   ]=1

                elif (not c1==c2) and s2==scaffold and not s1==scaffold:
                    
                    xB = 1
                    yB = y
                    xx=self.scaffold_coord(c2,xB)
                    yy=self.scaffold_coord(c2,yB)
                    segments[ max(0,min(xx,yy)) , min( slen, max(xx,yy) )   ]=1
 
        segments=list(segments.keys())
        segments.sort()
        #print("#msm:",scaffold,slen,segments,sep="\t")
        return(segments)

    def read_deserts(self,outfile,mapq=10,min_len=1000):
        if not self.bam_objects:
            self.bam_objects={}
            for b in self.bams:                
                self.bam_objects[b]= pysam.Samfile(b,"rb") 
        desert_edges=[]
        for bam in self.bam_objects.values():
            lastx=0
            lastr=-1
            for aln in bam.fetch(until_eof=True):
                if aln.mapq < mapq               : continue
                ref = aln.tid
                if not ref==lastr:
                    if lastr>-1:
                        lastname =bam.getrname(lastr)
                        lastl = self.ocontig_lengths[ lastname ]
                        if lastl - lastx > min_len:
                            desert_edges.append(( lastname,lastx  , 1 ))
                            desert_edges.append(( lastname,lastl  ,-1 ))
                    lastx=0
                deltax= aln.pos - lastx
                if deltax>min_len:
                    desert_edges.append(( bam.getrname(aln.tid),lastx  , 1 ))
                    desert_edges.append(( bam.getrname(aln.tid),aln.pos,-1 ))

                lastx=aln.pos
                lastr=aln.tid

        
        nbams = len( self.bams)
        state=0
        startx=0
        desert_edges.sort()
        y=0
        for contig,x,dy in desert_edges:

            y+=dy

            if state==0 and y == nbams:
                state=1
                startx=x
            elif state==1 and y<nbams:
                endx = x
                state=0
                if endx-startx>min_len:
                    print(contig,startx,endx,"readDesert",file=outfile)

            last_contig,last_x = contig,x


    def chicago_pairs(self,mapq=10,scaffold_contig_counts={},contig_scaffold_counts={},callback=False,bamfile=False,my_scaffolds=False):
        read_length = 100
        if not self.bam_objects:
            self.bam_objects={}
            for b in self.bams:                
                self.bam_objects[b]= pysam.Samfile(b,"rb") 
                


        if bamfile:
            bams_iter = iter( [ self.bam_objects[bamfile] ] )
        else:
            bams_iter = iter( self.bam_objects.values() )

        for bam in bams_iter:
#            print("#",bam)
            last_contig=False
            mri=0
            mr=[]
            last_scaffold=False
            for aln in bam.fetch(until_eof=True):
                if not aln.is_read1: continue
                track=False
                skip=False
                if aln.rnext<0:      continue
                if aln.mapq < mapq               : continue
                if BamTags.mate_mapq(aln) < mapq : continue
                if aln.is_duplicate : continue
                contig  = bam.getrname(aln.tid)
                contig2 = bam.getrname(aln.rnext)
                
                if not last_contig or not last_contig==contig:
                    if contig in self.masked_regions:
                        mr = self.masked_regions[contig]
                        mri=0
                    else:
                        mr=[]
                        mri=0
                c1,x = self.broken_coord(contig ,aln.pos  ,tuples=True)
                scaffold=self.contig_scaffold.get(c1)
                c2,y = self.broken_coord(contig2,aln.pnext,tuples=True)
                nscaffold=self.contig_scaffold.get(c2)

                if (not my_scaffolds==False) and (not ((scaffold in my_scaffolds) and (nscaffold in my_scaffolds))): continue

                if last_contig and not last_contig==contig:
#                    print("#done with contig",last_contig,list(contig_scaffold_counts.get(last_contig,{}).keys()))
                    for sc in contig_scaffold_counts.get(last_contig,{}).keys():  # for each hirise scaffold that this contig contributes a segment to:
                        if sc in scaffold_contig_counts:                   
                            if last_contig in scaffold_contig_counts[sc]:  # if we are still tracking the number of times we need to read this contig from a bam file for this scaffold
#                                print("##",sc,last_contig,scaffold_contig_counts[sc][last_contig])
                                scaffold_contig_counts[sc][last_contig] -= 1 # decrease the count by one
#                                print("##",sc,last_contig,scaffold_contig_counts[sc][last_contig])
                                if scaffold_contig_counts[sc][last_contig] ==0:
                                    del scaffold_contig_counts[sc][last_contig]
                            if len(scaffold_contig_counts[sc])==0:
#                                print("scaffold {} ready now",sc)
                                if callback: 
                                    callback(sc)
                                else:
                                    yield(sc,"done")
                                
                                del scaffold_contig_counts[sc]
#                                raise ScaffoldReady(sc)
#                            print("#s",last_contig,sc,scaffold_contig_counts.get(sc))
#                    print("#c",last_contig,contig_scaffold_counts.get(last_contig))
                last_contig = contig        
                last_scaffold = scaffold

                while mri<len(mr) and aln.pos > mr[mri][1]: mri+=1                        
                if mri<len(mr) and mr[mri][0]<aln.pos and aln.pos+read_length<mr[mri][1] : continue #this read maps to a masked region 

                if contig2 in self.masked_regions:
                    if [ (aa,bb) for aa,bb in self.masked_regions[contig2] if aa<aln.pnext and aln.pnext+read_length<bb ]: continue  #the pair maps to a masked region
                
                if not scaffold: continue
                if not nscaffold: continue

                xx=self.scaffold_coord(c1,x)
                yy=self.scaffold_coord(c2,y)

                #                yield(s1,self.scaffold_coord(c1,x),s2,self.scaffold_coord(c2,y))

                if self.scaffold_lengths[scaffold] < xx or  self.scaffold_lengths[nscaffold] < yy:
                    print("coordinate out of range error",xx,yy,scaffold,nscaffold,self.scaffold_lengths[scaffold], self.scaffold_lengths[nscaffold],contig,contig2,c1,c2,x,y)
                    raise Exception 


                yield( scaffold, nscaffold, xx, yy, c1, c2, aln.query_name )


    def chicago_pairs_for_scaffolds(self,mapq=10,callback=False,bamfile=False,scaffolds=[],contigs=False,minsep=0):
        read_length = 100
        if not self.bam_objects:
            self.bam_objects={}
            for b in self.bams:                
                self.bam_objects[b]= pysam.Samfile(b,"rb") 

        if bamfile:
            bams_iter = iter( [ self.bam_objects[bamfile] ] )
        else:
            bams_iter = iter( self.bam_objects.values() )


        ocontigs_to_hit = {}
        my_ocontigs=[]
        if contigs:
            my_ocontigs=contigs
        else:
            for scaffold in scaffolds:
                for oc in self.scaffold_ocontigs[scaffold]:
                    ocontigs_to_hit[oc]=1

            my_ocontigs = list(ocontigs_to_hit.keys())

        for bam in bams_iter:
#            print("#",bam)
            last_contig=False
            mri=0
            mr=[]
            last_scaffold=False
            for ocontig in my_ocontigs:
                for aln in bam.fetch(reference=ocontig):
                    if not aln.is_read1: continue
                    track=False
                    skip=False
                    if aln.rnext<0:      continue
                    if aln.mapq < mapq               : continue
                    if BamTags.mate_mapq(aln) < mapq : continue
                    if aln.is_duplicate : continue
                    contig  = bam.getrname(aln.tid)
                    contig2 = bam.getrname(aln.rnext)

                    if not last_contig or not last_contig==contig:
                        if contig in self.masked_regions:
                            mr = self.masked_regions[contig]
                            mri=0
                        else:
                            mr=[]
                            mri=0

                    if last_contig and not last_contig==contig:
                        pass
    #                    print("#done with contig",last_contig,list(contig_scaffold_counts.get(last_contig,{}).keys()))

                    c1,x = self.broken_coord(contig ,aln.pos  ,tuples=True)
                    scaffold=self.contig_scaffold.get(c1)
                    last_contig = contig        
                    last_scaffold = scaffold

                    while mri<len(mr) and aln.pos > mr[mri][1]: mri+=1                        
                    if mri<len(mr) and mr[mri][0]<aln.pos and aln.pos+read_length<mr[mri][1] : continue #this read maps to a masked region 

                    if contig2 in self.masked_regions:
                        if [ (aa,bb) for aa,bb in self.masked_regions[contig2] if aa<aln.pnext and aln.pnext+read_length<bb ]: continue  #the pair maps to a masked region

                    c2,y = self.broken_coord(contig2,aln.pnext,tuples=True)
                    nscaffold=self.contig_scaffold.get(c2)
                    if not scaffold: continue
                    if not nscaffold: continue

                    if not scaffold in scaffolds: continue
                    if not nscaffold in scaffolds: continue

                    xx=self.scaffold_coord(c1,x)
                    yy=self.scaffold_coord(c2,y)
                    if scaffold==nscaffold and abs(xx-yy)<minsep: continue

                    #                yield(s1,self.scaffold_coord(c1,x),s2,self.scaffold_coord(c2,y))

                    if self.scaffold_lengths[scaffold] < xx or  self.scaffold_lengths[nscaffold] < yy:
                        print("coordinate out of range error",xx,yy,scaffold,nscaffold,self.scaffold_lengths[scaffold], self.scaffold_lengths[nscaffold],contig,contig2,c1,c2,x,y)
                        raise Exception 

                    yield( scaffold, nscaffold, xx, yy, c1, c2, aln.query_name )

    def get_links(self,ocontigs,skipI=False,mapq=10,links={},contigs=False,raw=False,tuples=False,debug=False):
#        links={}

        if not self.bam_objects:
            self.bam_objects={}
            for b in self.bams:

                self.bam_objects[b]= pysam.Samfile(b,"rb") 

        for bam in self.bam_objects.values():
#            bam = pysam.Samfile(b,"rb") 

            for ocontig in ocontigs:
#                print("#",ocontig)
                mr=[]
                mri=0
                if ocontig in self.masked_regions: mr=self.masked_regions[ocontig]
                for aln in bam.fetch(reference=ocontig):
#                    if debug: print("#",aln.query_name,aln.pos,aln.pnext,aln.mapq,sep="\t")
                    if aln.rnext<0: continue
                    if aln.mapq < mapq               : continue
                    if BamTags.mate_mapq(aln) < mapq : continue
                    if aln.is_duplicate : continue

                    while mri<len(mr) and aln.pos > mr[mri][1]: mri+=1                        
                    if mri<len(mr) and mr[mri][0]<aln.pos and aln.pos<mr[mri][1] : 
                        if debug: print("#m1",ocontig,bam.getrname(aln.rnext),aln.pos,aln.pnext,mr[mri],aln.query_name,sep="\t")
                        continue #this read maps to a masked region 

                    if raw:
                        c1,x1 = ocontig,aln.pos
                    else:
                        c1,x1 = self.broken_coord(ocontig,aln.pos,tuples=tuples)
                        if not c1: continue

                    if contigs and not c1 in contigs: continue
                    ocontig2 = bam.getrname(aln.rnext)
                    
                    
                    if raw:
                        c2,x2 = ocontig2,aln.pnext
                    else:
                        c2,x2 = self.broken_coord(ocontig2,aln.pnext,tuples=tuples)
                        if not c2: continue
                    if contigs and not c2 in contigs: continue

                    if skipI and c1 == c2: continue

                    if ocontig2 in self.masked_regions:
                        overlapping_masked_segs = [ (aa,bb) for aa,bb in self.masked_regions[ocontig2] if aa<aln.pnext and aln.pnext<bb ]
                        if overlapping_masked_segs: 
                            if debug: print("#m2",ocontig,ocontig2, aln.pos,aln.pnext,overlapping_masked_segs ,aln.query_name,sep="\t")
                            continue  #the pair maps to a masked region

                    if not (c1,c2) in links: 
                        links[c1,c2]=[]
                    links[c1,c2].append((x1,x2))
                    if debug: print("#p",c1,c2,aln.query_name,aln.pos,aln.pnext,sep="\t")
                    #                print(self.broken_coord(ocontig,aln.pos),self.broken_coord(bam.getrname(aln.rnext),aln.pnext))
#        for c1,c2 in links:
#            print(c1,c2,links[c1,c2])
        return(links)

    def validate(self):
        contig_mins={}
        contig_maxs={}
        scaffold_mins={}
        scaffold_maxs={}
        contig_lengths={}
        counts={}
        #        scaffold_lengths={]
        lastx=0
        last_scaffold=False
        last_contig=(False,)
        for scaffold,contig,base,end,x,l in self.layout_lines:
            if last_scaffold and scaffold == last_scaffold and not x>lastx: 
                print(last_scaffold,scaffold,x,lastx)
                raise Exception 
            if x<0: raise Exception
            counts[contig,base]=counts.get((contig,base),0)+1
            contig_mins[contig,base]=min( contig_mins.get((contig,base),1e99) , x )
            contig_maxs[contig,base]=max( contig_maxs.get((contig,base),-1) , x )
            if (contig,base) in contig_lengths:
                if not contig_lengths[contig,base]==l: raise Exception 
            else:
                contig_lengths[contig,base]=l
            if (contig,base) == last_contig:
#                print(scaffold,contig,base,end,x,l,last_contig,lastx,x-lastx,"XXXX",sep="\t")
                if not l == x-lastx :
                    print(scaffold,contig,base,end,x,l,last_contig,lastx,sep="\t")
                    raise Exception 
            scaffold_mins[scaffold]=min( scaffold_mins.get(scaffold,1e99) , x )
            scaffold_maxs[scaffold]=max( scaffold_maxs.get(scaffold,-1)   , x )
            last_contig=(contig,base)
            lastx,last_scaffold=x,scaffold
        for contig,base in counts.keys():
            if not counts[contig,base]==2: raise Exception 
            if not contig_maxs[contig,base]-contig_mins[contig,base]==contig_lengths[contig,base]: 
                print("wtf?",contig,base,contig_mins[contig,base],contig_maxs[contig,base],contig_lengths[contig,base],contig_maxs[contig,base]-contig_mins[contig,base])
                raise Exception 

        self.setup_mapper()
        lastk=(1,0)
        lastl=0
        lastx=0
        for k in sorted(self.contigx.keys()):
            if k[0]==lastk[0]:
                #print("XX",k,contig_lengths.get(k),self.contigx[k],lastl,k[1]-lastk[1],lastl==k[1]-lastk[1],sep="\t")
                if not lastl==k[1]-lastk[1]:
                    print(k,contig_lengths.get(k),self.contigx[k],lastl,k[1]-lastk[1],lastl==k[1]-lastk[1],sep="\t")
                    raise Exception 
            lastk=k
            lastl = contig_lengths.get(k)
            lastx=self.contigx[k]
#        self.contigx={}
#        self.contigst={}


        print("#validation success")
#            contig_maxs[contig,base]=max( contig_maxs.get((contig,base),-1) , x )


    def add_breakpoints_ranges(self,breaks,debug=False,scores={},unlink_threshold=100.0, contig_break_threshold=100.0,score_types={},clipped_likelihood_for_gaps_only=True):
        """This function filters breaks expressed as spans according to several criteria before passing selected break points along to add_breakpoints()"""

        if not self.layout_lines or len(self.layout_lines)==0:
            self.make_trivial_layout_lines(debug=debug)
            self.validate()

        break_points = []
        category={}
        br={}
        scaffold_breaks={}
        for s,a,b,c in breaks:
            if not s in scaffold_breaks: scaffold_breaks[s]=[]

        for s,a,b,c in breaks:
            scaffold_breaks[s].append((a,b,c))

        for s in scaffold_breaks.keys():
             scaffold_breaks[s].sort()
        i=0
        offset={}
        while i< len(self.layout_lines)-1:
            #print(self.layout_lines[i])
            scaffold1,contig1,base1,end1,x1,l21 = self.layout_lines[  i  ]
            scaffold2,contig2,base2,end2,x2,l22 = self.layout_lines[ i+1 ]
            i+=1
            if not scaffold1 in scaffold_breaks: continue
            if not scaffold1 == scaffold2 : continue
            for a,b,c in scaffold_breaks[scaffold1]:
                if pairs_overlap((a,b),(x1,x2)) and (not (contig1,base1)==(contig2,base2)):  # you don't have to break a contig to make this break
                    if x1<=c and c<=x2:
                        this_offset_from_min_score = 0
                    else:
                        this_offset_from_min_score = min(abs(x1-c),abs(x2-c))

                    if (not (scaffold1,a,b,c) in category) or \
                       (category[scaffold1,a,b,c]=="break") or \
                       (category[scaffold1,a,b,c]=="gap" and offset[scaffold1,a,b,c]>this_offset_from_min_score):
                        category[scaffold1,a,b,c]="gap"
                        br[scaffold1,a,b,c]=int((x1+x2)/2)
                        offset[scaffold1,a,b,c] = this_offset_from_min_score
                elif x1<a and b<x2:
                    if not (scaffold1,a,b,c) in category:
                        if scores.get((scaffold1,a,b,c),-1e9) < contig_break_threshold:
                            category[scaffold1,a,b,c]="break"
                            br[scaffold1,a,b,c]=c
        categories={}
        for s,a,b,c in br.keys():
            if (category[s,a,b,c]=="break" and scores.get((s,a,b,c),100)< contig_break_threshold) or \
               (category[s,a,b,c]=="gap"   and scores.get((s,a,b,c),100)< unlink_threshold) :

                if clipped_likelihood_for_gaps_only and (score_types.get((s,a,b,c),"fine")=="clipped") and category[s,a,b,c]=="break": 
                    if debug: print("Skip break because we're only using clipped likelihood for un-scaffolding, not contig breaking",s,a,b,c,scores.get((s,a,b,c),-100),sep="\t")
                    if debug: print("X",s,a,b,c,br[s,a,b,c],category[s,a,b,c],score_types.get((s,a,b,c),"fine"),scores.get((s,a,b,c),-100))
                    continue
                
                break_points.append((s,br[s,a,b,c]))
                if True or debug: print("Q",s,a,b,c,br[s,a,b,c],category[s,a,b,c],score_types.get((s,a,b,c),"fine"),scores.get((s,a,b,c),-100))
                categories[s,br[s,a,b,c]]=category[s,a,b,c]
        self.add_breakpoints(break_points,debug=debug)

    def add_breakpoints(self,breaks,debug=False,categories={}):
        """breaks is a list of (scaffold,x) pairs where scaffold is the id of a scaffold in the assembly, 
and x is a basepair offset from the scaffold 5' end of the scaffold where a break should be added.  If x
falls in a gap between contigs, the linkage is broken.  If x falls in a contig, the contig is broken."""
        import networkx as nx

        g = nx.DiGraph()

        scaffold_breaks={}

        for s,x in breaks:
            if not s in scaffold_breaks: scaffold_breaks[s]=[]

        for s,x in breaks:
            scaffold_breaks[s].append(x)

        for s in scaffold_breaks.keys():
             scaffold_breaks[s].sort()

        if debug: print(scaffold_breaks)

        new_layout_lines=[]

        scaffold_broken=False

        break_indices={}

        i=0
        break_counter=0
#        debug=False

        edge_breaks={}
        elength={}
        scaffold={}
        xcoord={}
#
#  Convert the layout lines into a directed graph
#
        contig_length={}
        while i< len(self.layout_lines)-1:
            #print(self.layout_lines[i])
            scaffold1,contig1,base1,end1,x1,l21 = self.layout_lines[  i  ]
            scaffold2,contig2,base2,end2,x2,l22 = self.layout_lines[ i+1 ]
            scaffold[contig1,base1]=scaffold1
            scaffold[contig2,base2]=scaffold2
            i+=1
            if scaffold1==scaffold2: 
                if scaffold1 in scaffold_breaks and (not (contig1,base1)==(contig2,base2)) and len([z for   z in scaffold_breaks[scaffold1] if x1<z and z<x2])>0:  
                    pass # this is an OO join to remove
                else:
                    g.add_edge( (contig1,base1,end1),(contig2,base2,end2) )
                    elength[ (contig1,base1,end1),(contig2,base2,end2) ]=x2-x1
                    if (contig1,base1)==(contig2,base2): contig_length[contig1,base2]=x2-x1

                    if not ( scaffold1==scaffold2 and scaffold1 in scaffold_breaks):
                        continue

                    else:
                        break_positions = [z for   z in scaffold_breaks[scaffold1] if x1<z and z<x2]

                        if len(break_positions)==0: 
                            continue

                        edge_breaks[(contig1,base1,end1),(contig2,base2,end2)]=[]
                        for x in break_positions: # these's a break in here!

                            edge_breaks[(contig1,base1,end1),(contig2,base2,end2)].append( x-x1 )

        for e1,e2 in edge_breaks.keys():
#            print(e1,e2)
            next_nodes= g.neighbors(e2) 
            prev_nodes= g.predecessors(e1) 
#            print(prev_nodes)
            if len(list(next_nodes))>1: 
                raise Exception 
            if len(list(prev_nodes))>1: 
                raise Exception 
            next_node=False
            prev_node=False
            
            if len(list(next_nodes)) == 0:
                continue
            if len(list(previous_nodes)) == 0:
                continue

            if next_nodes:
                next_node = list(next_nodes)[0]
            if prev_nodes:
                prev_node = list(prev_nodes)[0]
            g.remove_edge(e1,e2)

            #print(e1,e2,edge_breaks[e1,e2],scaffold[e1[0],e1[1]],next_node,prev_node)

            if not (e1[0],e1[1])==(e2[0],e2[1]):
                #print("this is just an oo link to udo")
                #this is just an edge to remove
                continue

            ocontig = e1[0]
            
            if e1[2] in ["5",5]: #forward

                
                
                if next_node: g.remove_edge(e2,next_node)
                g.remove_node(e2)

                #                          a        b
                #         prev    e1     x1          x2       e2   next
                #  ---------->     --------------------------->    ----------------->

                #just implement the first and last break.  (this is not necessarily good if the contig is really, really long. XXX TODO )
                edge_breaks[e1,e2].sort()
                x1=edge_breaks[e1,e2][0]
                x2=edge_breaks[e1,e2][-1]

                for i in range(len(edge_breaks[e1,e2])-1):
                    xa = edge_breaks[e1,e2][i  ]
                    xb = edge_breaks[e1,e2][i+1]
                    new_node_a = ( e1[0] , e1[1]+xa, '5' )
                    new_node_b = ( e1[0] , e1[1]+xa, '3' )
                    g.add_edge(new_node_a,new_node_b)
                    elength[ new_node_a,new_node_b ] = xb-xa
                    scaffold[e1[0] , e1[1]+xa] = scaffold[e1[0],e1[1]]

                new_node_1 = (e1[0],e1[1]   ,'3')
                new_node_2 = (e1[0],e1[1]+x2,'5')
                new_node_3 = (e1[0],e1[1]+x2,'3')
                scaffold[e1[0],e1[1]+x2] = scaffold[e1[0],e1[1]]

                if debug: print("b",e1,e2,prev_node,next_node,x1,x2,edge_breaks[e1,e2],new_node_1,new_node_2,   sep="\t")

                g.add_edge(e1,new_node_1)
                g.add_edge(new_node_2,new_node_3)
                if next_node: g.add_edge(new_node_3,next_node)

                LL = elength[e1,e2]
                if next_node: LL2 = elength[e2,next_node]
                elength[ e1,new_node_1 ] = x1
                elength[ new_node_2,new_node_3 ] = LL - x2
                if next_node: elength[ new_node_3,next_node  ] = LL2


                if debug: print("b",e1,e2,prev_node,next_node,x1,x2,edge_breaks[e1,e2],new_node_1,new_node_2, elength.get(( e1,new_node_1 )),elength.get(( new_node_2,new_node_3 )), elength.get((new_node_3,next_node)), sep="\t")
#                print("zzz",contig_length[e1[0],e1[1]],LL,elength[e1,e2],edge_breaks[e1,e2],x1,x2,elength[e1,e2],elength[ e1,new_node_1 ],elength[ new_node_2,new_node_3 ])

            else:                #reverse
                
                if prev_node: g.remove_edge(prev_node,e1)
                g.remove_node(e1)

                #just implement the first and last break.  (this is not necessarily good if the contig is really, really long. XXX TODO )
                edge_breaks[e1,e2].sort()                
                x1=elength[e1,e2]-edge_breaks[e1,e2][-1]
                x2=elength[e1,e2]-edge_breaks[e1,e2][0]


                #                          b        a
                #         prev    e1     x2          x1       e2   next
                #  ---------->    <----------------------------    ----------------->
                #                  A      B           C

                for i in range(len(edge_breaks[e1,e2])-1):
                    xb = elength[e1,e2]-edge_breaks[e1,e2][i  ]
                    xa = elength[e1,e2]-edge_breaks[e1,e2][i+1]
                    new_node_a = ( e1[0] , e2[1]+xa, '5' )
                    new_node_b = ( e1[0] , e2[1]+xa, '3' )
                    g.add_edge(new_node_a,new_node_b)
                    elength[ new_node_a,new_node_b ] = abs(xb-xa)
                    scaffold[e1[0] , e2[1]+xa] = scaffold[e1[0],e2[1]]

                new_node_C = (e1[0],e2[1]   ,'3') #
                new_node_B = (e1[0],e2[1]+x2,'5') #
                new_node_A = (e1[0],e2[1]+x2,'3') #
                scaffold[e1[0],e2[1]+x2] = scaffold[e1[0],e2[1]]

                g.add_edge(new_node_C,e2) #
                g.add_edge(new_node_A,new_node_B) #
                if prev_node: g.add_edge(prev_node,new_node_A) #

                LL =  elength[e1,e2]
                if prev_node: LL2 = elength[prev_node,e1]
                elength[ new_node_C, e2 ] = x1
                elength[ new_node_A,new_node_B ] = LL - x2
                if prev_node: elength[ prev_node, new_node_A  ] = LL2
#                print("zzW",edge_breaks[e1,e2],x1,x2,elength[e1,e2],elength[ new_node_1, e2 ],elength[ new_node_3,new_node_2 ])
            
#        for e1,e2 in elength.keys():
#            if elength[e1,e2]<0:
#                print("wtf?",e1,e2,elength[e1,e2])
#
#  Compute the new layout coordinates and scaffold IDs
#
        scaffold_incr={}
        node_x={}
        contig_length={}
        new_scaffold={}
        for n in sorted(g.nodes()):
            #print(n,scaffold[n[0],n[1]],g.degree(n),g.in_degree(n),g.out_degree(n))
            if g.in_degree(n)==0:
                scaffold_id=scaffold[n[0],n[1]]
                new_scaffold_id = scaffold_id 
                if scaffold_id in scaffold_incr: 
                    new_scaffold_id = scaffold_id + "."+str(scaffold_incr.get(scaffold_id,0))
                new_scaffold[n]=new_scaffold_id
                scaffold_incr[scaffold_id] = scaffold_incr.get(scaffold_id,0)+1
                ll=list(nx.dfs_preorder_nodes(g,n))
                xx=0
                lastl=False
                for l in ll:
                    if lastl:
                        xx+=elength[ lastl,l ]
                        #print("z",l,lastl,elength[lastl,l],node_x[lastl],xx,new_scaffold[n],sep="\t")
                        if (lastl[0],lastl[1])==(l[0],l[1]):
                            contig_length[lastl[0],lastl[1]]=elength[ lastl,l ]
                    #print(new_scaffold_id,l,xx)
                    node_x[l]=xx
                    lastl = l
#
#  Build the new layout_lines list
#
        self.scaffold_lengths={}
        new_layout = []
        for n in g.nodes():
            if g.in_degree(n)==0:

                new_scaffold_id = new_scaffold[n]
                ll=list(nx.dfs_preorder_nodes(g,n))

                for l in ll:
#                    print("w",new_scaffold_id,l[0],l[1],l[2],node_x[l])
#                    print("z",new_scaffold_id,l[0],l[1],l[2],node_x[l],contig_length[l[0],l[1]])
                    new_layout.append( (new_scaffold_id,l[0],l[1],l[2],node_x[l],contig_length[l[0],l[1]] ) )
#                    new_layout.append( (new_scaffold_id,l[0],l[1],l[2],node_x[l],elength[l[0],l[1]] ) )
                    self.scaffold_lengths[new_scaffold_id] = max( node_x[l] , self.scaffold_lengths.get(new_scaffold_id,0))

                    
        
        self.layout_lines=new_layout
#        self.layout_lines.sort( key=lambda x: (x[0],x[4]) )
        self.index_ocontigs()

    def load_assembly(self,infile):
        self.scaffold_lengths={}
        self.bams=[]
        self.layout_lines=[]
        self.model_params={}
        self.ocontig_lengths={}
        self.contig_lengths={}
        self.scaffold_ocontigs={}
        self.bam_objects=False
#        self.model_params={}
        f=open(infile,"r")
        path = ""
        #if infile.startswith("/"):
        path = os.path.dirname(os.path.normpath(infile))
        for l in f:
            #print(l)
            if l[0]=="M":
                c=l.strip().split()
                self.model_params[c[1]] = eval(" ".join(c[2:]))[0]
            if l[0]=="B":
                #self.bams.append(l[2:].strip())
                possible_path = l[2:].strip()
                if possible_path.startswith("/"):
                   self.bams.append(l[2:].strip())
                else:
                    self.bams.append(os.path.join(path,l[2:].strip()))   
            if l[0]=="S":
                #self.shotgun.append(l[2:].strip())
                possible_path = l[2:].strip()
                if possible_path.startswith("/"):
                    self.shotgun.append(l[2:].strip())
                else:
                    self.shotgun.append(os.path.join(path,l[2:].strip()))
            if l[0]=="P": # Layout
                self.layout_lines.append( parse_pline2(l,self.scaffold_lengths) )
                scaffold,contig,base,end,x,lll = self.layout_lines[-1]
                self.contig_lengths[contig+"_"+str(base)]=lll
                self.contig_lengths[contig,base]=lll
                self.scaffold_ocontigs[scaffold] = self.scaffold_ocontigs.get(scaffold,[])+[contig]
            if l[0]=="D": # Masked regions (from D for shotgun depth)
                dummy,ocontig,a,b = l.strip().split()
                if not ocontig in self.masked_regions: self.masked_regions[ocontig]=[]
                self.masked_regions[ocontig].append((int(a),int(b)))
            if l[0]=="L": # Length of input contigs, same as bam header lengths  
                dummy,ocontig,l = l.strip().split()
                #if not ocontig in self.masked_regions: self.masked_regions[ocontig]=[]
                self.ocontig_lengths[ocontig] = int(l)
        f.close()
        self.index_ocontigs()
#        print("setup mapper")
        self.setup_mapper(debug=True)

    def save_assembly(self,outfile):
        f=open(outfile,"w")
        for b in self.shotgun:
            f.write("S {}\n".format(b))
        for b in self.bams:
            f.write("B {}\n".format(b))
        #f.write("M {}\n".format(self.model_params))
        if self.model_params:
            for k in sorted(self.model_params.keys()):
                f.write("M {} {}\n".format(k,(self.model_params[k],)))
        for l in self.layout_lines:
            f.write("P "+" ".join(map(str,l))+"\n")
        for oc in sorted(self.masked_regions.keys()):
            for a,b in self.masked_regions[oc]:
                f.write("D {} {} {}\n".format(oc,a,b))
        for oc in sorted(self.ocontig_lengths.keys()):
            f.write("L {} {}\n".format(oc,self.ocontig_lengths[oc]))

#        self.masked_regions={}
        f.close()

    def scaffold_mapping():
        for i in range(len(self.layout_lines)):
            scaffold1,contig1,base1,end1,x1,l11 = self.layout_lines[i  ]
            if end1 == "3":
                yield("{}_{}".format(contig1,base1),scaffold1)
        
    def as_edges(self):
        for i in range(len(self.layout_lines)-1):
            scaffold1,contig1,base1,end1,x1,l21 = self.layout_lines[i  ]
            scaffold2,contig2,base2,end2,x2,l22 = self.layout_lines[i+1]
            if scaffold1==scaffold2 :
                if contig1==contig2 and base1==base2: 
                    yield(["#edge:","{}_{}.{}".format(contig1,base1,end1) ,"{}_{}.{}".format(contig2,base2,end2) ,{'length':l21,'contig':True}])
                else:
                    yield(["#edge:","{}_{}.{}".format(contig1,base1,end1) ,"{}_{}.{}".format(contig2,base2,end2) ,{'length':x2-x1,'contig':False}])
#                    yield(["#edge:",contig1+"."+end1,contig2+"."+end2,{'length':x2-x1,'contig':False}])

    def contig_coords(self):
        seen = {}
        for l in self.layout_lines:
            scaffold1,contig1,base1,end1,x1,l21 = l
            if not (contig1,base1) in seen:
                yield( {"contig":contig1,"base":base1,"scaffold":scaffold1,"span":(x1,x1+l21),"flipped":end1 in [3,"3"] })

            seen[contig1,base1]=True

    def dump_layout(self,f):
#        f=open(outfile,"w")
        for s in self.scaffold_lengths.keys():
            f.write("{} slen\n".format(self.scaffold_lengths[s]))
        for l in self.layout_lines:
            scaffold1,contig1,base1,end1,x1,l21 = l
            f.write("p: {} {}_{}.{} {} - -1 {} {} {}\n".format(scaffold1,contig1,base1,end1,x1,self.scaffold_lengths[scaffold1],l21,False))
#        f.close()
        #p: 1 Scaffold105707_1.3 3019 - -1 85132546 3048 False

    def load_playout(self,layoutfile):
        #print(layoutfile)
        self.scaffold_lengths={}
#        self.contig_lengths={}
        if not os.path.exists(layoutfile):
            raise Exception 
        self.layout_lines = [ parse_pline(c,self.scaffold_lengths) for c in open(layoutfile,"rt") if c[:2]=="p:" ]
        #for l in self.layout_lines:
            #print(l)
        self.index_ocontigs()

    def check_bams(self):
        for bamfile in self.bams:
            if not os.path.exists(bamfile):
                raise Exception 

    def set_modelparams(self,datamodel):
        if os.path.exists(datamodel):
            self.model_params = eval(open(datamodel).read())
        else:
            self.model_params=datamodel

if __name__=="__main__":

    import sys
    import argparse

    parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument('-d','--debug',default=False,action="store_true",help="Turn on debugging ouput")
    parser.add_argument('-E','--edges',default=False,action="store_true",help="Print the layout as edges.")
#    parser.add_argument('-E','--edges',default=False,action="store_true",help="Print the layout as edges.")
    parser.add_argument('-D','--dump_layout',default=False,action="store_true",help="Print the layout as edges.")
    #parser.add_argument('-p','--progress',default=False,action="store_true",help="Print progress info")
    parser.add_argument('-M','--datamodel',default=False, help="Name of a file containing the model parameter dictionary")
    parser.add_argument('-m','--mask'   ,default=[], action="append",help="Name of a file containing regions of the input assembly to mask out")
    parser.add_argument('-b','--bamfile',default=[], action="append",help="A chicago bam file")
    parser.add_argument('-S','--shotgunbamfile',default=[], action="append",help="A shotgun bam file")
    #parser.add_argument('-p','--progress',default=False,action="store_true",help="Print progress info")
    parser.add_argument('-L','--layout',default=False,help="File containing a scaffolding layout in 'p:' lines.")
    parser.add_argument('-o','--outfile',default=False,help="Filename for serialization of the assembly file.")
    parser.add_argument('-i','--infile',default=False,help="Filename for serialised assembly input file.")

    args = parser.parse_args()
    
#    print(args)
    if args.infile:
        asf = HiriseAssembly()
        asf.load_assembly(args.infile)
#        asf.validate()
    else:
        asf = HiriseAssembly( {"bams":args.bamfile, "datamodel": args.datamodel, "layout": args.layout , "debug":args.debug, "shotgun": args.shotgunbamfile} )
    if args.debug: print("# done with intital load")
    asf.validate()
#    asf.index_ocontigs()

    if args.datamodel:
        asf.set_modelparams(args.datamodel)

    for segments_file in args.mask:
        asf.add_mask_regions(filename=segments_file)
        asf.merge_masked_regions()

    if args.outfile:
        asf.save_assembly( args.outfile )
        
    if args.edges:
        for e in asf.as_edges():
            print("\t".join(map(str,e)))

#    asf.add_breakpoints( [( '1', 1692 ),('1',10000)] )
#    asf.add_breakpoints( [( '1', 1692 ),('1',2000),('1',10000)] )
#    asf.add_breakpoints( [( '1', 1692 ),('1',2000),('1',2100),('1',10000)] )

#    asf.validate()

#    asf.load_ocontig_info()
#    print("#loaded ocontig info")
#    asf.load_ocontig_info()
#    print("#loaded ocontig info")
#    print(asf.get_links(["Scaffold116325"],links={}))
#    print(asf.get_links(["Scaffold116325"],links={}))
#    asf.get_scaffold_links("1000",skipI=True)

#1691
    if args.dump_layout:
        asf.dump_layout(sys.stdout)
#!/usr/bin/env python3 # # Copyright 2015 Dovetail Genomics LLC # # This class is designed to encapsulate a Hirise Assembly. Options currently used are:  bams, bam, datamodel, layout #return self.ocontig_lengths[m.group(1)] #dummy,scaffold,x,y,leng = l.strip().split() #x=int(x) #y=int(y) #leng=int(leng) #        y=min(y,leng) #        this_end = int(this_end) #        self.scaffold_ocontigs={} #            if not scaffold in self.scaffold_ocontigs: #                self.scaffold_ocontigs[scaffold]={} #            self.scaffold_ocontigs[scaffold][self.contig_ocontig[contig]]=1 #print("ocb:",b,ocb[b]) #                    print(ocontig,"1") #                    print(ocontig,"2") #print("zxzx",ocontig,x,b,new_contig,self.contig_lengths.get(new_contig,False)) #                        print(ocontig,"3") #                print(ocontig,"4",new_contig,new_contig in self.contig_lengths) #print("zxzx",ocontig,x,b,new_contig,self.contig_lengths.get(new_contig,False)) #        debug=True #        if debug: #            for contig,base in self.contigx.keys(): #                print("m:",contig,base,self.contigx[contig,base],self.contigst[contig,base]) #                print("#",scaffold) #        model = ces. #                if logfile: logfile.write("{} {} {} {} {} {}\n".format(scaffold,stretch_start,x,minx,low_point,slen)) #                if logfile: logfile.write("{} {} {} {} {} {} {} {}\n".format(scaffold,stretch_start,x,minx,low_point,slen,min_t1,max_t1)) #print(self.ocontig_lengths) #            if debug: print(scaffold,ocontig,1,"5",0,self.ocontig_lengths[ocontig],sep="\t") # # #cc=self.contig_length(contig1+"_"+str(base1)) #            print( r,self.rc_buffer[x],self.rc_buffer[x+1],r-self.rc_buffer[x][0]) #        print( r,self.rc_buffer[x],self.rc_buffer[x+1],r-self.rc_buffer[x][0]) #        dtype=[('x',int),('t',int),('y',int),('b',int),('r',int)] #                buffer.append((aln.pos,aln.rnext,aln.pnext,i,aln.tid)) #n_other_windows(buffer[i:j],minlinks,bins,self.binsize) #            if debug: print(tid,buffer[i][0],n_in_window,nw,mybin,[ (ii,bc[ii]) for ii in bc.keys() if bc[ii]>=minlinks and not ii==mybin ],[ get_bin(x,self.binsize) for x in buffer[i:j] ]) #            gc.disable() #                    if debug: print("#",aln.pos,aln.is_duplicate,aln.mapq,BamTags.mate_mapq(aln),aln.query_name,len(locations)) #            gc.enable() #                outfile.write("{} {} {}\n".format(ocontig,locations[i],locations[j-1])) #                outfile.write("{} {} {}\n".format(ocontig,locations[i],locations[j-1])) #,mapper,mask): #            print(ocontig) #print("#msm:",scaffold,ocontig,a,b,x,y,xx,yy,s1,s2,scaffold,sep="\t") #print("#msm:",scaffold,slen,segments,sep="\t") #            print("#",bam) #                    print("#done with contig",last_contig,list(contig_scaffold_counts.get(last_contig,{}).keys())) # for each hirise scaffold that this contig contributes a segment to: # if we are still tracking the number of times we need to read this contig from a bam file for this scaffold #                                print("##",sc,last_contig,scaffold_contig_counts[sc][last_contig]) # decrease the count by one #                                print("##",sc,last_contig,scaffold_contig_counts[sc][last_contig]) #                                print("scaffold {} ready now",sc) #                                raise ScaffoldReady(sc) #                            print("#s",last_contig,sc,scaffold_contig_counts.get(sc)) #                    print("#c",last_contig,contig_scaffold_counts.get(last_contig)) #this read maps to a masked region #the pair maps to a masked region #                yield(s1,self.scaffold_coord(c1,x),s2,self.scaffold_coord(c2,y)) #            print("#",bam) #                    print("#done with contig",last_contig,list(contig_scaffold_counts.get(last_contig,{}).keys())) #this read maps to a masked region #the pair maps to a masked region #                yield(s1,self.scaffold_coord(c1,x),s2,self.scaffold_coord(c2,y)) #        links={} #            bam = pysam.Samfile(b,"rb") #                print("#",ocontig) #                    if debug: print("#",aln.query_name,aln.pos,aln.pnext,aln.mapq,sep="\t") #this read maps to a masked region #the pair maps to a masked region #                print(self.broken_coord(ocontig,aln.pos),self.broken_coord(bam.getrname(aln.rnext),aln.pnext)) #        for c1,c2 in links: #            print(c1,c2,links[c1,c2]) #        scaffold_lengths={] #                print(scaffold,contig,base,end,x,l,last_contig,lastx,x-lastx,"XXXX",sep="\t") #print("XX",k,contig_lengths.get(k),self.contigx[k],lastl,k[1]-lastk[1],lastl==k[1]-lastk[1],sep="\t") #        self.contigx={} #        self.contigst={} #            contig_maxs[contig,base]=max( contig_maxs.get((contig,base),-1) , x ) This function filters breaks expressed as spans according to several criteria before passing selected break points along to add_breakpoints() #print(self.layout_lines[i]) # you don't have to break a contig to make this break breaks is a list of (scaffold,x) pairs where scaffold is the id of a scaffold in the assembly,  and x is a basepair offset from the scaffold 5' end of the scaffold where a break should be added.  If x falls in a gap between contigs, the linkage is broken.  If x falls in a contig, the contig is broken. #        debug=False # #  Convert the layout lines into a directed graph # #print(self.layout_lines[i]) # this is an OO join to remove # these's a break in here! #            print(e1,e2) #            print(prev_nodes) #print(e1,e2,edge_breaks[e1,e2],scaffold[e1[0],e1[1]],next_node,prev_node) #print("this is just an oo link to udo") #this is just an edge to remove #forward #                          a        b #         prev    e1     x1          x2       e2   next #  ---------->     --------------------------->    -----------------> #just implement the first and last break.  (this is not necessarily good if the contig is really, really long. XXX TODO ) #                print("zzz",contig_length[e1[0],e1[1]],LL,elength[e1,e2],edge_breaks[e1,e2],x1,x2,elength[e1,e2],elength[ e1,new_node_1 ],elength[ new_node_2,new_node_3 ]) #reverse #just implement the first and last break.  (this is not necessarily good if the contig is really, really long. XXX TODO ) #                          b        a #         prev    e1     x2          x1       e2   next #  ---------->    <----------------------------    -----------------> #                  A      B           C # # # # # # #                print("zzW",edge_breaks[e1,e2],x1,x2,elength[e1,e2],elength[ new_node_1, e2 ],elength[ new_node_3,new_node_2 ]) #        for e1,e2 in elength.keys(): #            if elength[e1,e2]<0: #                print("wtf?",e1,e2,elength[e1,e2]) # #  Compute the new layout coordinates and scaffold IDs # #print(n,scaffold[n[0],n[1]],g.degree(n),g.in_degree(n),g.out_degree(n)) #print("z",l,lastl,elength[lastl,l],node_x[lastl],xx,new_scaffold[n],sep="\t") #print(new_scaffold_id,l,xx) # #  Build the new layout_lines list # #                    print("w",new_scaffold_id,l[0],l[1],l[2],node_x[l]) #                    print("z",new_scaffold_id,l[0],l[1],l[2],node_x[l],contig_length[l[0],l[1]]) #                    new_layout.append( (new_scaffold_id,l[0],l[1],l[2],node_x[l],elength[l[0],l[1]] ) ) #        self.layout_lines.sort( key=lambda x: (x[0],x[4]) ) #        self.model_params={} #if infile.startswith("/"): #print(l) #self.bams.append(l[2:].strip()) #self.shotgun.append(l[2:].strip()) # Layout # Masked regions (from D for shotgun depth) # Length of input contigs, same as bam header lengths #if not ocontig in self.masked_regions: self.masked_regions[ocontig]=[] #        print("setup mapper") #f.write("M {}\n".format(self.model_params)) #        self.masked_regions={} #                    yield(["#edge:",contig1+"."+end1,contig2+"."+end2,{'length':x2-x1,'contig':False}]) #        f=open(outfile,"w") #        f.close() #p: 1 Scaffold105707_1.3 3019 - -1 85132546 3048 False #print(layoutfile) #        self.contig_lengths={} #for l in self.layout_lines: #print(l) #    parser.add_argument('-E','--edges',default=False,action="store_true",help="Print the layout as edges.") #parser.add_argument('-p','--progress',default=False,action="store_true",help="Print progress info") #parser.add_argument('-p','--progress',default=False,action="store_true",help="Print progress info") #    print(args) #        asf.validate() #    asf.index_ocontigs() #    asf.add_breakpoints( [( '1', 1692 ),('1',10000)] ) #    asf.add_breakpoints( [( '1', 1692 ),('1',2000),('1',10000)] ) #    asf.add_breakpoints( [( '1', 1692 ),('1',2000),('1',2100),('1',10000)] ) #    asf.validate() #    asf.load_ocontig_info() #    print("#loaded ocontig info") #    asf.load_ocontig_info() #    print("#loaded ocontig info") #    print(asf.get_links(["Scaffold116325"],links={})) #    print(asf.get_links(["Scaffold116325"],links={})) #    asf.get_scaffold_links("1000",skipI=True) #1691