def reader(working_dirs, references, use_reference, annotations={}):
    for name, sequence in references:
        features = annotations.get(sequence, [])

        if use_reference:
            readers = [reference_reader(sequence)]
        else:
            readers = []

        readers.extend(evidence_reader(working_dir, name) for working_dir in working_dirs)

        active_features = []
        feature_pos = 0

        for i in xrange(len(sequence)):
            if i % 10000 == 0:
                grace.status("%s %s" % (name, grace.pretty_number(i)))

            active_features = [item for item in active_features if item.location.nofuzzy_end > i]
            while feature_pos < len(features) and features[feature_pos].location.nofuzzy_start <= i:
                active_features.append(features[feature_pos])
                feature_pos += 1

            for is_insertion in (True, False):
                yield Calls(name, i, is_insertion, [item.next() for item in readers], active_features)

        for reader in readers:
            for item in reader:
                raise grace.Error("Unexpected extra data in evidence file")

    grace.status("")

def unmill(rast, mill_points):
    padx = max(item[0] for item in mill_points)
    pady = max(item[1] for item in mill_points)
    sy,sx = rast.shape
    padded = numpy.zeros((sy+pady*2,sx+padx*2), rast.dtype)
    padded[pady:pady+sy,padx:padx+sx] = rast
    
    result = rast.copy()
    
    mill_points = sorted(mill_points, key=lambda item:item[1])
    
    for y in xrange(sy):
        row = result[y]
        grace.status('%d' % (sy-y))
        for ox,oy,oheight in mill_points:
            numpy.maximum(row,padded[pady+y-oy,padx-ox:padx+sx-ox]-oheight,row)
    
    grace.status('')
    
    #old_height = 0
    #for x,y,height in sorted(mill_points, key=lambda item:item[2]):
    #    if height != old_height:
    #        numpy.subtract(padded, height-old_height, padded)
    #        old_height = height
    #    print x,y,height
    #    numpy.maximum(
    #        result,
    #        padded[pady-y:pady+sy-y,padx-x:padx+sx-x],
    #        result
    #    )
    return result

    def make_plot(self, plots_name, plot_names, iterator, maximum, color='0,0,0', scale_type='log', windowing='maximum'):
        grace.status('Write '+plots_name)
        filename = self.prefix + plots_name + '.igv'
        f = open(filename, 'wb')
        
        height = max(10,int(100.0/math.sqrt(len(plot_names)))) #...
        
        print >> f, '#track viewLimits=0:%(maximum)f autoScale=off scaleType=%(scale_type)s windowingFunction=%(windowing)s maxHeightPixels=200:%(height)d:1 color=%(color)s' % locals()
        print >> f, '\t'.join(
            [ 'Chromosome', 'Start', 'End', 'Feature'] + [ self.label_prefix + item for item in plot_names ]
        )
        for name, pos, depths in iterator:
            print >> f, '\t'.join(
                [ name, str(pos), str(pos+1), 'F' ] +
                [ str(item) for item in depths ]
            )
        
        f.close()
        grace.status('')
        
        if self.genome:
            #One igvtools process at a time
            self.wait_for_igv()

            p = io.run([
                'igvtools', 'toTDF',
                filename, 
                self.prefix + plots_name + '.tdf',
                self.genome,
                '-f', 'max,mean'
            ], stdin=None, stdout=None)

            self.processes.append((p, filename))

def main():
    with Timer("100 status updates"):
        for i in xrange(100):
            grace.status(str(i))
        grace.status("")

    with Timer("100 parallel processes"):

        @parallel_for(xrange(100))
        def loop(i):
            pass

    with Timer("100 parallel processes updating statuses"):

        @parallel_for(xrange(100))
        def loop(i):
            grace.status(str(i))
            grace.status("")

    with Timer("Nested processes"):

        @parallel_for(xrange(3))
        def _(item):
            print item

            @parallel_for(xrange(3))
            def _(item2):
                print item, item2

def guess_quality_offset(*filenames):
    grace.status('Guessing quality offset')
    try:
        min_value = chr(255)
        #max_value = chr(0)
        any_reads = False
        for filename in filenames:
            for i, item in enumerate(read_sequences(filename, qualities=True)):
                if i > 100000: break
                if len(item) == 2: continue
                
                any_reads = True
                min_value = min(min_value, min(item[2]))
                #max_value = max(max_value, max(item[2]))

        if not any_reads: 
            return 33
        
        low = ord(min_value)
        #high = ord(max_value)
        #print 'Quality chars in range %d-%d in %s' % (low,high,filename)
        if low < 59: return 33 #Sanger and Illumina 1.8+
        return 64 #Illumina pre 1.8
    finally:
        grace.status('')

def main():
    #print dir(hello)
    import sys
    print sys.modules['__main__'].__file__
    
    a = nesoni.future(func,'a',[])
    b = nesoni.future(func,'b',[])
    c = nesoni.future(func,'c',[a,b])
    d = nesoni.future(func,'d',[a,b])
    
    c()
    d()

    #print legion.coordinator().get_cores()
    #legion.coordinator().job(hello)

    with Timer('100 status updates'):
        for i in xrange(100):
            grace.status(str(i))
        grace.status('')

    with Timer('100 parallel threads'):
        @nesoni.thread_for(xrange(100))
        def loop(i):
            pass

    with Timer('100 parallel processes'):
        nesoni.parallel_for(xrange(100))(do_nothing)

    with Timer('100 parallel processes updating statuses'):
        nesoni.parallel_for(xrange(100))(do_status)

 def wait_for_igv(self):
     while self.processes:
         p, filename = self.processes.pop()
         grace.status('igvtools processing '+filename)
         assert p.wait() == 0, 'igvtools tile failed'
         grace.status('')
         if self.delete_igv:
             os.unlink(filename)

 def cut_raster(self, raster):
     if self.bit_ball:
         bit = ball_mill(self.res_bit_radius)
     else:
         bit = end_mill(self.res_bit_radius)
     print 'Unmilling'
     inraster = unmill(raster, bit)
     print 'Again'
     inagain = unmill(inraster, ball_mill(self.res_bit_radius))
     print 'Done'
     
     def point_score(x,y,z):
         return 1.0 / ((inagain[y,x]-inraster[y,x])+self.res_bit_radius)
     
     spin = 0.0
     
     cut_z = 0
     min_z = numpy.minimum.reduce(raster.flatten())         
     while True:
         cut_z -= self.res_cutting_depth
         grace.status('%d %d  %f' % (cut_z, min_z, spin))
         inmask = inraster <= cut_z
         self.cut_inmask(cut_z, inmask, 
                         #in_first =self.res_bit_radius/3.0+self.res_finishing_clearance,
                         #out_first=self.res_bit_radius/3.0, 
                         #in_step  =self.res_bit_radius,
                         #out_step =self.res_bit_radius/3.0,
                         in_first  = self.res_finishing_clearance + self.res_cutting_step*spin,
                         out_first = 0.0,
                         in_step   = self.res_cutting_step,
                         out_step  = 0.0,
                         point_score=point_score)
         
         if self.finish:
             infinish_mask = inraster <= cut_z                 
             infinish_mask_lower = inraster <= (cut_z-self.res_cutting_depth)
             self.cut_inmask(cut_z, 
                             infinish_mask & ~infinish_mask_lower, #erode(infinish_mask_lower,self.res_horizontal_step * 1.5), 
                             #in_first =self.res_horizontal_step*2,
                             #out_first=self.res_horizontal_step,
                             #in_step  =self.res_horizontal_step*2,
                             #out_step =self.res_horizontal_step,
                             in_first =self.res_horizontal_step,
                             out_first=0.0,
                             in_step  =self.res_horizontal_step,
                             out_step =0.0,
                             point_score=point_score,
                             down_cut=False
                             )
             self.cut_inmask(cut_z, infinish_mask, down_cut=False)
         
         if cut_z <= min_z: 
             break
         
         spin = (spin-GOLDEN)%1.0
     grace.status('')

 def eat(f):
     for line in f:
         if line.startswith('@'):
             if sam_header_sent[0]: continue
         else:
             n_seen[0] += 1
             if n_seen[0] % 100000 == 0:
                 grace.status('%s alignments produced' % grace.pretty_number(n_seen[0]))
         sam_eater.write_raw(line)
     sam_header_sent[0] = True

 def build_shrimp_mmap(self, cs=False):
     suffix = '-cs' if cs else '-ls'
     
     grace.status('Building SHRiMP mmap')
     io.execute([
         'gmapper' + suffix,
         '--save', self.object_filename('reference' + suffix),
         self.reference_fasta_filename(),
     ])
     grace.status('')

 def eat(process):
     for line in process.stdout:
         if line.startswith('@'):
             if sam_header_sent[0]: continue
         else:
             n_seen[0] += 1
             if n_seen[0] % 100000 == 0:
                 grace.status('%s alignments produced' % grace.pretty_number(n_seen[0]))
         sam_eater.write_raw(line)
         
     assert process.wait() == 0, 'shrimp failed'
     sam_header_sent[0] = True

 def setup(self):
     grace.status('Load depths')
     self.sample_names = [ os.path.split(dirname)[1] for dirname in self.working_dirs ]
     self.workspaces = [ working_directory.Working(dirname, must_exist=True) for dirname in self.working_dirs ]
     
     self.depths = [ item.get_depths() for item in self.workspaces ]
     #self.depths = list(legion.imap(lambda item: item.get_object('depths.pickle.gz'), self.workspaces, local=True))
     
     self.any_pairs = any(item.param['any_pairs'] for item in self.workspaces)
     grace.status('')
     
     lengths = self.workspaces[0].get_reference().get_lengths()
     self.chromosome_names = [ name for name, length in lengths ]
     self.lengths = dict(lengths)
     
     self.processes = [ ]

def read_userplot(filename):
    grace.status('Load '+filename)
    
    f = open(filename,'rb')
    lines = f.readlines()
    f.close()
    
    n = 0
    while n < len(lines) and lines[n].startswith('#'):
        n += 1
    headers = lines[:n]

    is_multiplot = (n > 0)

    data = [ tuple([ int(item) for item in lines[i].strip().split() ])
             for i in xrange(n,len(lines)) ]

    grace.status('')
    return headers, is_multiplot, data

def guess_quality_offset(filename):
    grace.status('Guessing quality offset')
    try:
        min_value = chr(255)
        #max_value = chr(0)
        for i, item in enumerate(read_sequences(filename, qualities=True)):
            if len(item) == 2: return 33 #Not fastq
            
            min_value = min(min_value, min(item[2]))
            #max_value = max(max_value, max(item[2]))
            
            if i >= 100000: break
        
        low = ord(min_value)
        #high = ord(max_value)
        #print 'Quality chars in range %d-%d in %s' % (low,high,filename)
        if low < 59: return 33 #Sanger and Illumina 1.8+
        return 64 #Illumina pre 1.8
    finally:
        grace.status('')

    def run(self):
        workspace = working_directory.Working(self.output_dir)
        workspace.setup_reference(self.reference)
        workspace.update_param(snp_cost=self.snp_cost)

        # assert os.path.exists(self.reference), 'Reference file does not exist'
        # reference_filename = workspace._object_filename('reference.fa')
        # if os.path.exists(reference_filename):
        #   os.unlink(reference_filename)
        # os.symlink(os.path.relpath(self.reference, self.output_dir), reference_filename)

        bam_filename = io.abspath(self.output_dir, "alignments.bam")
        bam_prefix = io.abspath(self.output_dir, "alignments")

        if sam.is_bam(self.input):
            sort_input_filename = self.input
            temp_filename = None
        else:
            temp_filename = io.abspath(self.output_dir, "temp.bam")
            sort_input_filename = temp_filename
            writer = io.Pipe_writer(temp_filename, ["samtools", "view", "-S", "-b", "-"])
            f = open(self.input, "rb")
            while True:
                data = f.read(1 << 20)
                if not data:
                    break
                writer.write(data)
            writer.close()
            f.close()

        grace.status("Sort")

        # io.execute([
        #    'samtools', 'sort', '-n', sort_input_filename, bam_prefix
        # ])
        sam.sort_bam(sort_input_filename, bam_prefix, by_name=True)

        if temp_filename is not None:
            os.unlink(temp_filename)

        grace.status("")

def _make_inner(action):
    timestamp = coordinator().time()
    assert timestamp > LOCAL.time, 'Time running in reverse.'
    
    cores = action.cores_required()
    if cores > 1:
        coordinator().trade_cores(1, cores)
    try:        
        config.write_colored_text(sys.stderr, '\n'+action.describe()+'\n')
        
        if LOCAL.abort_make and not selection.matches(LOCAL.do_selection, [action.shell_name()]):
            raise grace.Error('%s would be run. Stopping here.' % action.ident())
        
        old_status = grace.status(action.shell_name())
        try:
            _run_and_save_state(action, timestamp)
        finally:
            grace.status(old_status)
    finally:
        if cores > 1:
            coordinator().trade_cores(cores, 1)

    def find_maximum_depth(self):
        grace.status('Finding maximum depth')        
        #if self.strand_specific:
        #    iterator = itertools.chain(
        #        self.iter_over(lambda item: item.ambiguous_depths[0], zeros=False),
        #        self.iter_over(lambda item: item.ambiguous_depths[1], zeros=False),
        #    )
        #else:
        #    iterator = self.iter_over_unstranded(lambda item: item.ambiguous_depths, zeros=False)

        maximum = 1
        norm_maximum = 1.0
        normalize = self.normalizer()
        n = 0
        #for name, pos, depths in iterator:
        #    if n % 1000000 == 0:
        #        grace.status('Finding maximum depth %s %s' % (name, grace.pretty_number(pos)))
        #    n += 1
        #    
        #    maximum = max(maximum,max(depths))
        #    norm_maximum = max(norm_maximum,max(normalize(depths)))

        #futures = [ ]        
        #for name in self.chromosome_names:
        #    for i, depth in enumerate(self.depths):
        #        if self.strand_specific:
        #            #this = futures.append(max(
        #            #    max(depth[name].ambiguous_depths[0]),
        #            #    max(depth[name].ambiguous_depths[1])
        #            #)
        #            
        #            futures.append( (name,i,legion.future(max, depth[name].ambiguous_depths[0])) )
        #            futures.append( (name,i,legion.future(max, depth[name].ambiguous_depths[1])) )
        #        else:
        #            #this = max(iter_add(depth[name].ambiguous_depths[0],depth[name].ambiguous_depths[1]))
        #            futures.append( (name,i,legion.future(lambda item: max(iter_add(item[0],item[1])), depth[name].ambiguous_depths)) )
        #
        #for name, i, future in futures:
        #    grace.status('Finding maximum depth %s %s' % (name, self.sample_names[i]))
        #    this = future()
        #    maximum = max(maximum, this)
        #    norm_maximum = max(norm_maximum, self.norm_mult[i] * this)

        for name in self.chromosome_names:
            for i, depth in enumerate(self.depths):
                grace.status('Finding maximum depth %s %s' % (name, self.sample_names[i]))
                if self.strand_specific:
                    this = max(
                        depth[name].ambiguous_depths[0].maximum(),
                        depth[name].ambiguous_depths[1].maximum()
                        )
                else:
                    this = (depth[name].ambiguous_depths[0] + depth[name].ambiguous_depths[1]).maximum()
                maximum = max(maximum, this)
                norm_maximum = max(norm_maximum, self.norm_mult[i] * this)
        
        self.maximum = maximum
        self.norm_maximum = norm_maximum
        grace.status('')

    def run(self):
        workspace = working_directory.Working(self.output_dir)        
        workspace.setup_reference(self.reference)
        workspace.update_param(snp_cost = self.snp_cost)
        
        #assert os.path.exists(self.reference), 'Reference file does not exist'
        #reference_filename = workspace._object_filename('reference.fa')
        #if os.path.exists(reference_filename):
        #   os.unlink(reference_filename)
        #os.symlink(os.path.relpath(self.reference, self.output_dir), reference_filename)
        
        bam_filename = io.abspath(self.output_dir, 'alignments.bam')
        bam_prefix = io.abspath(self.output_dir, 'alignments')

        if sam.is_bam(self.input):
            sort_input_filename = self.input
            temp_filename = None
        else:
            temp_filename = io.abspath(self.output_dir, 'temp.bam')
            sort_input_filename = temp_filename
            writer = io.Pipe_writer(temp_filename, ['samtools', 'view', '-S', '-b', '-'])
            f = open(self.input, 'rb')
            while True:
                data = f.read(1<<20)
                if not data: break
                writer.write(data)
            writer.close()
            f.close()
        
        grace.status('Sort')
        
        io.execute([
            'samtools', 'sort', '-n', sort_input_filename, bam_prefix
        ])
        
        if temp_filename is not None:
            os.unlink(temp_filename)
        
        grace.status('')

 def calculate_norm_mult(self):
     grace.status('Calculating normalization')
     
     totals = [ 0 ] * len(self.working_dirs)
     
     for i in xrange(len(self.workspaces)):
         for name in self.lengths:
             totals[i] += self.depths[i][name].depths[0].total() + self.depths[i][name].depths[1].total()
     
     #for name, pos, depths in self.iter_over_unstranded(lambda item: item.depths):
     #    if pos % 1000000 == 0:
     #        grace.status('Calculating normalization %s %s' % (name, grace.pretty_number(pos)))
     #    for i, depth in enumerate(depths):
     #        totals[i] += depth
     
     grace.status('')
     
     nonzero = [ item for item in totals if item ]
     geomean = math.exp(sum( math.log(item) for item in nonzero ) / len(nonzero))
     self.norm_mult = [
         1.0 if not item else geomean / item
         for item in totals
     ]

 def run(self):
     working = io.Workspace(self.output_dir, must_exist=False)
 
     for filename in self.files:
         reader = io.Table_reader(filename)
         
         name = os.path.splitext(os.path.split(filename)[1])[0]
         
         rname = None
         files = None
         for record in reader:
             if record['Chromosome'] != rname:
                 if files: 
                     for item in files: 
                         item.close()
                 rname = record['Chromosome']
                 grace.status('Convert '+name+' '+rname)
                 files = [
                     open(working / (
                         name + 
                         '-' + grace.filesystem_friendly_name(rname) + 
                         '-' + grace.filesystem_friendly_name(item) + '.userplot'
                     ), 'wb')
                     for item in reader.headings[4:]
                 ]
                 pos = 0
             assert int(record['Start']) == pos and int(record['End']) == pos + 1
             
             for val, f in zip(record.values()[4:], files):
                 print >> f, val
             
             pos += 1
         
         if files: 
             for item in files: 
                 item.close()
         grace.status('')