

import pysam
from numpy import *
from multiprocessing import Pool
import random as rnd
from .kernel import *
import sys
from . import readscounter


class KeyboardInterruptError(Exception):
    
    pass


def nhnoncontrol(uniqueratio, threshold, kernellength, nthreads=4):
    
    bgscore = sim_replicate_nthreads(run_times=200, uniqueratio=uniqueratio,
                                     nthreads=nthreads, kernellength=kernellength, threshold=threshold)
    
    cutoff = bgscore['mean'] + bgscore['std'] * threshold
    
    return cutoff


def dhnoncontrol(uniqueratio, threshold, kernellength, nthreads=4):
    
    bgscore = sim_replicate_nthreads(run_times=200, uniqueratio=uniqueratio,
                                     nthreads=nthreads, kernellength=kernellength, threshold=threshold)
    
    cutoff = bgscore['mean'] + bgscore['std'] * threshold

    return cutoff


def nhcontrol(bamfile, chromosome, paired, chrlength, ultratio, filted_region,maxinsert, kernellength = 600, threshold = 4):
    
    """
        region: chr:start-end
        ultraio = chrlength * uniqratio / chr_total_reads
        filter region
        
    """
    
    region = chromosome + ':' + str(1) + '-' + str(chrlength)
    
    readscount = readscounter.nhreadscounter(bamfile, region, paired,  maxinsert=maxinsert)
    
    kernel = smooth_kernel(kernellength)
    
    kernel_score = list()

    for i in sorted(kernel):
        
        kernel_score.append(kernel[i])

    threshold = filted_region.threshold/100

    for site in readscount:
        
        if readscount[site] > threshold:
            
            readscount[site] = threshold

    smoothed_result = correlate(array(readscount), kernel_score)

    ultratiolist = list()

    ultratiolist.append(ultratio)

    smoothed_result = correlate(smoothed_result, ultratiolist)

    #scores = list()
    
    bg_mean = smoothed_result.mean()
    
    bg_std = smoothed_result.std()
    
    bg_threshold = bg_mean + threshold * bg_std
    
    #bgscore['rand_mean'] = bg_mean
    
    #bgscore['rand_std'] = bg_std
    
    cutoff = bg_threshold

    return cutoff


def nhuniquerate(bamfile, chromosome, paired, fregion, regionstart=1, regionend = -1, maxinsert = 100000):
    
    samfile = pysam.Samfile(bamfile)

    ref_lengths = samfile.lengths

    sam_ref = samfile.references

    refere_ncenumber = samfile.nreferences

    if regionend == -1:
        
        for i in range(refere_ncenumber):
            
            if sam_ref[i] == chromosome:
                
                regionend = ref_lengths[i]

    region = chromosome + ':' + str(regionstart) + '-' + str(regionend)
    
    region_length = regionend - regionstart
        
    nhreadscount = readscounter.nhreadscounter(bamfile, region, paired, maxinsert=maxinsert)

    totaluniq = len(nhreadscount)  + 0.0
    
    uniquerate = totaluniq/region_length
    
    return uniquerate


def dhuniquerate(bamfile, chromosome, regionstart=1, regionend=-1):
    
    """
        Count unique Rate in a region
    
    """
    
    samfile = pysam.Samfile(bamfile)

    ref_lengths = samfile.lengths

    sam_ref = samfile.references

    refere_ncenumber = samfile.nreferences

    if regionend == -1:
        
        for i in range(refere_ncenumber):
            
            if sam_ref[i] == chromosome:
                
                regionend = ref_lengths[i]

    region = chromosome + ':' + str(regionstart) + '-' + str(regionend)
    
    region_length = regionend - regionstart
    
    dhreadscount = readscounter.dhreadscounter(bamfile, region)
    
    totaluniq = len(dhreadscount) + 0.0
    
    uniquerate = totaluniq/region_length
    
    return uniquerate


def ultratio(chrlength, uniqueratio, chrtotalreads, frcount):
    """
    ultraio = chrlength * uniqueratio / chr_total_reads
    """
    ultratio = chrlength * uniqueratio / (chrtotalreads - frcount)

    return ultratio


def sim_replicate_nthreads(run_times=1000, uniqueratio=1, kernellength = 600, threshold = 4, nthreads = 2):
    # randomthresh = list()
    
    pars = list()
    
    for i in range(0,run_times):
        
        par=dict()
        
        par['uniqueratio'] = uniqueratio
        
        par['kernellength'] = kernellength
        
        par['threshold'] = threshold
        
        pars.append(par)
    
    pool=Pool(nthreads)
    
    outscore = dict()
    
    try:
        randomthresh = pool.map(sim_bg_thread_worker, pars)

        summean = 0.0
    
        sumstd = 0.0

        for randscore in randomthresh:

            randmean = randscore['rand_mean']
            
            randstd = randscore['rand_std']
            # print (randmean, randstd)
            summean = summean + randmean
            
            sumstd = sumstd + randstd

        mean_of_mean = summean/run_times
        
        mean_of_std = sumstd/run_times
        # print ('mean_of_mean',mean_of_mean, 'mean_of_std',mean_of_std)

        outscore['mean'] = mean_of_mean
        
        outscore['std'] = mean_of_std
        #return (mean_of_mean, mean_of_std)

        pool.close()
        
        return outscore

    except KeyboardInterrupt:
        
        pool.terminate()
        
        print ("You cancelled the program!")
        
        sys.exit(1)
        
    except Exception as e:
        
        print ('got exception: %r, terminating the pool' % (e,))
        
        pool.terminate()
        
        print ('pool is terminated')
        
    finally:
        # print ('joining pool processes')
        pool.join()
    # print ('join complete')
    # pool.join()
    # pool.close()


def sim_bg_thread_worker(par):

    try:
        
        uniqueratio=par['uniqueratio']
        
        kernellength = par['kernellength']
        
        threshold = par['threshold']

        kernel = smooth_kernel(length=kernellength)

        sim_genome_size = int(1e5)

        total_reads = int(sim_genome_size * uniqueratio)

        region_site = list(range(0,sim_genome_size))

        sim_uniqsite = rnd.sample(region_site, total_reads)

        rand_reads_count = list()

        for i in range(0,sim_genome_size):
            
            rand_reads_count.append(0)

        kernel_score = list()

        for i in sorted(kernel):
            kernel_score.append(kernel[i])



        kdesmooth_result = dict()

        for i in range(0,total_reads):
            
            rand_number = int(rnd.uniform(0,total_reads))
            
            rand_reads = sim_uniqsite[rand_number]

            rand_reads_count[rand_reads] = rand_reads_count[rand_reads] + 1.0

        smoothed_result = correlate(array(rand_reads_count), kernel_score)

        scores = list()

        rand_mean = smoothed_result.mean()
        
        rand_std = smoothed_result.std()
        
        total_sum = smoothed_result.sum()
        
        rand_threshhold = rand_mean + threshold * rand_std

        higher_count = 0

        for now_site in kdesmooth_result:
            
            if kdesmooth_result[now_site] > rand_threshhold:
                
                higher_count = higher_count + 1

        # print (total_sum, rand_mean, rand_std, rand_threshhold, higher_count, total_reads)

        randscore = dict()
        
        randscore['rand_mean'] = rand_mean
        
        randscore['rand_std'] = rand_std

        return randscore

    except KeyboardInterrupt:
        
        raise KeyboardInterruptError()


def get_bpc(bamfile,  hotspots, jobtype, filted_region, nthreads, maxinsert = 100000):

    #bpc average readscount per basepare

    # total_reads = 0
    #
    total_length = 0
    #
    # samfile = pysam.Samfile(bamfile)

    pars = list()

    for hotspot_now in hotspots:

        par = dict()

        par['bamfile'] = bamfile

        par['hotspot'] = hotspot_now

        par['jobtype'] = jobtype

        par['filted_region'] = filted_region

        par['maxinsert'] = maxinsert

        pars.append(par)

        total_length = hotspot_now.end - hotspot_now.start + 1 + total_length


    # print ("total length %s" % total_length)

    pool = Pool(nthreads)

    try:

        reads_count = pool.map(bpc_runner, pars)

        total_reads = 0.0

        for count_now in reads_count:

            total_reads = total_reads + count_now

        bpc = (total_reads+0.0)/total_length

        pool.close()

        return bpc

    except KeyboardInterrupt:

        pool.terminate()

        print ("You cancelled the program!")

        sys.exit(1)

    except Exception as e:

        print ('got exception: %r, terminating the pool' % (e,))

        pool.terminate()

        print ('pool is terminated')

    finally:
        # print ('joining pool processes')
        pool.join()
        # print ('join complete')


def bpc_runner(par):

    try:

        bamfile = par['bamfile']

        hotspot = par['hotspot']

        jobtype = par['jobtype']

        filted_region = par['filted_region']

        maxinsert = par['maxinsert']

        start_site = hotspot.start

        end_site = hotspot.end

        whether_in_fr = 0

        chromosome = hotspot.chromosome

        hotspotregio = chromosome + ':' + str(start_site) + '-' + str(end_site)

        hotspotreads = 0

        for i in range(start_site, end_site + 1):

            parentscare = int(i/100)

            if chromosome in filted_region:

                if parentscare in filted_region[chromosome]:

                    whether_in_fr = 1

        if whether_in_fr == 0:

            readscount = dict()

            if jobtype == 'nhsingle':

                readscount = readscounter.nhreadscounter(bamfile = bamfile, region=hotspotregio, paired=False, maxinsert=maxinsert)

            elif jobtype == 'nhpaired':

                readscount = readscounter.nhreadscounter(bamfile = bamfile, region=hotspotregio, paired=True, maxinsert=maxinsert)

            elif jobtype == 'dh':

                readscount = readscounter.dhreadscounter(bamfile = bamfile, region = hotspotregio)

            else:

                print ("%s count type error!!!!" % jobtype)

                sys.exit(1)

            for i in readscount:

                hotspotreads = hotspotreads+readscount[i]

        # print (hotspotreads)
        return hotspotreads

    except KeyboardInterrupt:

        raise KeyboardInterruptError()
b016830f1077a5202b7b6d035736707de4888da3
<mask token>


class KeyboardInterruptError(Exception):
    pass


def nhnoncontrol(uniqueratio, threshold, kernellength, nthreads=4):
    bgscore = sim_replicate_nthreads(run_times=200, uniqueratio=uniqueratio,
        nthreads=nthreads, kernellength=kernellength, threshold=threshold)
    cutoff = bgscore['mean'] + bgscore['std'] * threshold
    return cutoff


def dhnoncontrol(uniqueratio, threshold, kernellength, nthreads=4):
    bgscore = sim_replicate_nthreads(run_times=200, uniqueratio=uniqueratio,
        nthreads=nthreads, kernellength=kernellength, threshold=threshold)
    cutoff = bgscore['mean'] + bgscore['std'] * threshold
    return cutoff


<mask token>


def ultratio(chrlength, uniqueratio, chrtotalreads, frcount):
    """
    ultraio = chrlength * uniqueratio / chr_total_reads
    """
    ultratio = chrlength * uniqueratio / (chrtotalreads - frcount)
    return ultratio


def sim_replicate_nthreads(run_times=1000, uniqueratio=1, kernellength=600,
    threshold=4, nthreads=2):
    pars = list()
    for i in range(0, run_times):
        par = dict()
        par['uniqueratio'] = uniqueratio
        par['kernellength'] = kernellength
        par['threshold'] = threshold
        pars.append(par)
    pool = Pool(nthreads)
    outscore = dict()
    try:
        randomthresh = pool.map(sim_bg_thread_worker, pars)
        summean = 0.0
        sumstd = 0.0
        for randscore in randomthresh:
            randmean = randscore['rand_mean']
            randstd = randscore['rand_std']
            summean = summean + randmean
            sumstd = sumstd + randstd
        mean_of_mean = summean / run_times
        mean_of_std = sumstd / run_times
        outscore['mean'] = mean_of_mean
        outscore['std'] = mean_of_std
        pool.close()
        return outscore
    except KeyboardInterrupt:
        pool.terminate()
        print('You cancelled the program!')
        sys.exit(1)
    except Exception as e:
        print('got exception: %r, terminating the pool' % (e,))
        pool.terminate()
        print('pool is terminated')
    finally:
        pool.join()


def sim_bg_thread_worker(par):
    try:
        uniqueratio = par['uniqueratio']
        kernellength = par['kernellength']
        threshold = par['threshold']
        kernel = smooth_kernel(length=kernellength)
        sim_genome_size = int(100000.0)
        total_reads = int(sim_genome_size * uniqueratio)
        region_site = list(range(0, sim_genome_size))
        sim_uniqsite = rnd.sample(region_site, total_reads)
        rand_reads_count = list()
        for i in range(0, sim_genome_size):
            rand_reads_count.append(0)
        kernel_score = list()
        for i in sorted(kernel):
            kernel_score.append(kernel[i])
        kdesmooth_result = dict()
        for i in range(0, total_reads):
            rand_number = int(rnd.uniform(0, total_reads))
            rand_reads = sim_uniqsite[rand_number]
            rand_reads_count[rand_reads] = rand_reads_count[rand_reads] + 1.0
        smoothed_result = correlate(array(rand_reads_count), kernel_score)
        scores = list()
        rand_mean = smoothed_result.mean()
        rand_std = smoothed_result.std()
        total_sum = smoothed_result.sum()
        rand_threshhold = rand_mean + threshold * rand_std
        higher_count = 0
        for now_site in kdesmooth_result:
            if kdesmooth_result[now_site] > rand_threshhold:
                higher_count = higher_count + 1
        randscore = dict()
        randscore['rand_mean'] = rand_mean
        randscore['rand_std'] = rand_std
        return randscore
    except KeyboardInterrupt:
        raise KeyboardInterruptError()


<mask token>

<mask token>


class KeyboardInterruptError(Exception):
    pass


def nhnoncontrol(uniqueratio, threshold, kernellength, nthreads=4):
    bgscore = sim_replicate_nthreads(run_times=200, uniqueratio=uniqueratio,
        nthreads=nthreads, kernellength=kernellength, threshold=threshold)
    cutoff = bgscore['mean'] + bgscore['std'] * threshold
    return cutoff


def dhnoncontrol(uniqueratio, threshold, kernellength, nthreads=4):
    bgscore = sim_replicate_nthreads(run_times=200, uniqueratio=uniqueratio,
        nthreads=nthreads, kernellength=kernellength, threshold=threshold)
    cutoff = bgscore['mean'] + bgscore['std'] * threshold
    return cutoff


def nhcontrol(bamfile, chromosome, paired, chrlength, ultratio,
    filted_region, maxinsert, kernellength=600, threshold=4):
    """
        region: chr:start-end
        ultraio = chrlength * uniqratio / chr_total_reads
        filter region
        
    """
    region = chromosome + ':' + str(1) + '-' + str(chrlength)
    readscount = readscounter.nhreadscounter(bamfile, region, paired,
        maxinsert=maxinsert)
    kernel = smooth_kernel(kernellength)
    kernel_score = list()
    for i in sorted(kernel):
        kernel_score.append(kernel[i])
    threshold = filted_region.threshold / 100
    for site in readscount:
        if readscount[site] > threshold:
            readscount[site] = threshold
    smoothed_result = correlate(array(readscount), kernel_score)
    ultratiolist = list()
    ultratiolist.append(ultratio)
    smoothed_result = correlate(smoothed_result, ultratiolist)
    bg_mean = smoothed_result.mean()
    bg_std = smoothed_result.std()
    bg_threshold = bg_mean + threshold * bg_std
    cutoff = bg_threshold
    return cutoff


<mask token>


def ultratio(chrlength, uniqueratio, chrtotalreads, frcount):
    """
    ultraio = chrlength * uniqueratio / chr_total_reads
    """
    ultratio = chrlength * uniqueratio / (chrtotalreads - frcount)
    return ultratio


def sim_replicate_nthreads(run_times=1000, uniqueratio=1, kernellength=600,
    threshold=4, nthreads=2):
    pars = list()
    for i in range(0, run_times):
        par = dict()
        par['uniqueratio'] = uniqueratio
        par['kernellength'] = kernellength
        par['threshold'] = threshold
        pars.append(par)
    pool = Pool(nthreads)
    outscore = dict()
    try:
        randomthresh = pool.map(sim_bg_thread_worker, pars)
        summean = 0.0
        sumstd = 0.0
        for randscore in randomthresh:
            randmean = randscore['rand_mean']
            randstd = randscore['rand_std']
            summean = summean + randmean
            sumstd = sumstd + randstd
        mean_of_mean = summean / run_times
        mean_of_std = sumstd / run_times
        outscore['mean'] = mean_of_mean
        outscore['std'] = mean_of_std
        pool.close()
        return outscore
    except KeyboardInterrupt:
        pool.terminate()
        print('You cancelled the program!')
        sys.exit(1)
    except Exception as e:
        print('got exception: %r, terminating the pool' % (e,))
        pool.terminate()
        print('pool is terminated')
    finally:
        pool.join()


def sim_bg_thread_worker(par):
    try:
        uniqueratio = par['uniqueratio']
        kernellength = par['kernellength']
        threshold = par['threshold']
        kernel = smooth_kernel(length=kernellength)
        sim_genome_size = int(100000.0)
        total_reads = int(sim_genome_size * uniqueratio)
        region_site = list(range(0, sim_genome_size))
        sim_uniqsite = rnd.sample(region_site, total_reads)
        rand_reads_count = list()
        for i in range(0, sim_genome_size):
            rand_reads_count.append(0)
        kernel_score = list()
        for i in sorted(kernel):
            kernel_score.append(kernel[i])
        kdesmooth_result = dict()
        for i in range(0, total_reads):
            rand_number = int(rnd.uniform(0, total_reads))
            rand_reads = sim_uniqsite[rand_number]
            rand_reads_count[rand_reads] = rand_reads_count[rand_reads] + 1.0
        smoothed_result = correlate(array(rand_reads_count), kernel_score)
        scores = list()
        rand_mean = smoothed_result.mean()
        rand_std = smoothed_result.std()
        total_sum = smoothed_result.sum()
        rand_threshhold = rand_mean + threshold * rand_std
        higher_count = 0
        for now_site in kdesmooth_result:
            if kdesmooth_result[now_site] > rand_threshhold:
                higher_count = higher_count + 1
        randscore = dict()
        randscore['rand_mean'] = rand_mean
        randscore['rand_std'] = rand_std
        return randscore
    except KeyboardInterrupt:
        raise KeyboardInterruptError()


<mask token>

<mask token>


class KeyboardInterruptError(Exception):
    pass


def nhnoncontrol(uniqueratio, threshold, kernellength, nthreads=4):
    bgscore = sim_replicate_nthreads(run_times=200, uniqueratio=uniqueratio,
        nthreads=nthreads, kernellength=kernellength, threshold=threshold)
    cutoff = bgscore['mean'] + bgscore['std'] * threshold
    return cutoff


def dhnoncontrol(uniqueratio, threshold, kernellength, nthreads=4):
    bgscore = sim_replicate_nthreads(run_times=200, uniqueratio=uniqueratio,
        nthreads=nthreads, kernellength=kernellength, threshold=threshold)
    cutoff = bgscore['mean'] + bgscore['std'] * threshold
    return cutoff


def nhcontrol(bamfile, chromosome, paired, chrlength, ultratio,
    filted_region, maxinsert, kernellength=600, threshold=4):
    """
        region: chr:start-end
        ultraio = chrlength * uniqratio / chr_total_reads
        filter region
        
    """
    region = chromosome + ':' + str(1) + '-' + str(chrlength)
    readscount = readscounter.nhreadscounter(bamfile, region, paired,
        maxinsert=maxinsert)
    kernel = smooth_kernel(kernellength)
    kernel_score = list()
    for i in sorted(kernel):
        kernel_score.append(kernel[i])
    threshold = filted_region.threshold / 100
    for site in readscount:
        if readscount[site] > threshold:
            readscount[site] = threshold
    smoothed_result = correlate(array(readscount), kernel_score)
    ultratiolist = list()
    ultratiolist.append(ultratio)
    smoothed_result = correlate(smoothed_result, ultratiolist)
    bg_mean = smoothed_result.mean()
    bg_std = smoothed_result.std()
    bg_threshold = bg_mean + threshold * bg_std
    cutoff = bg_threshold
    return cutoff


def nhuniquerate(bamfile, chromosome, paired, fregion, regionstart=1,
    regionend=-1, maxinsert=100000):
    samfile = pysam.Samfile(bamfile)
    ref_lengths = samfile.lengths
    sam_ref = samfile.references
    refere_ncenumber = samfile.nreferences
    if regionend == -1:
        for i in range(refere_ncenumber):
            if sam_ref[i] == chromosome:
                regionend = ref_lengths[i]
    region = chromosome + ':' + str(regionstart) + '-' + str(regionend)
    region_length = regionend - regionstart
    nhreadscount = readscounter.nhreadscounter(bamfile, region, paired,
        maxinsert=maxinsert)
    totaluniq = len(nhreadscount) + 0.0
    uniquerate = totaluniq / region_length
    return uniquerate


def dhuniquerate(bamfile, chromosome, regionstart=1, regionend=-1):
    """
        Count unique Rate in a region
    
    """
    samfile = pysam.Samfile(bamfile)
    ref_lengths = samfile.lengths
    sam_ref = samfile.references
    refere_ncenumber = samfile.nreferences
    if regionend == -1:
        for i in range(refere_ncenumber):
            if sam_ref[i] == chromosome:
                regionend = ref_lengths[i]
    region = chromosome + ':' + str(regionstart) + '-' + str(regionend)
    region_length = regionend - regionstart
    dhreadscount = readscounter.dhreadscounter(bamfile, region)
    totaluniq = len(dhreadscount) + 0.0
    uniquerate = totaluniq / region_length
    return uniquerate


def ultratio(chrlength, uniqueratio, chrtotalreads, frcount):
    """
    ultraio = chrlength * uniqueratio / chr_total_reads
    """
    ultratio = chrlength * uniqueratio / (chrtotalreads - frcount)
    return ultratio


def sim_replicate_nthreads(run_times=1000, uniqueratio=1, kernellength=600,
    threshold=4, nthreads=2):
    pars = list()
    for i in range(0, run_times):
        par = dict()
        par['uniqueratio'] = uniqueratio
        par['kernellength'] = kernellength
        par['threshold'] = threshold
        pars.append(par)
    pool = Pool(nthreads)
    outscore = dict()
    try:
        randomthresh = pool.map(sim_bg_thread_worker, pars)
        summean = 0.0
        sumstd = 0.0
        for randscore in randomthresh:
            randmean = randscore['rand_mean']
            randstd = randscore['rand_std']
            summean = summean + randmean
            sumstd = sumstd + randstd
        mean_of_mean = summean / run_times
        mean_of_std = sumstd / run_times
        outscore['mean'] = mean_of_mean
        outscore['std'] = mean_of_std
        pool.close()
        return outscore
    except KeyboardInterrupt:
        pool.terminate()
        print('You cancelled the program!')
        sys.exit(1)
    except Exception as e:
        print('got exception: %r, terminating the pool' % (e,))
        pool.terminate()
        print('pool is terminated')
    finally:
        pool.join()


def sim_bg_thread_worker(par):
    try:
        uniqueratio = par['uniqueratio']
        kernellength = par['kernellength']
        threshold = par['threshold']
        kernel = smooth_kernel(length=kernellength)
        sim_genome_size = int(100000.0)
        total_reads = int(sim_genome_size * uniqueratio)
        region_site = list(range(0, sim_genome_size))
        sim_uniqsite = rnd.sample(region_site, total_reads)
        rand_reads_count = list()
        for i in range(0, sim_genome_size):
            rand_reads_count.append(0)
        kernel_score = list()
        for i in sorted(kernel):
            kernel_score.append(kernel[i])
        kdesmooth_result = dict()
        for i in range(0, total_reads):
            rand_number = int(rnd.uniform(0, total_reads))
            rand_reads = sim_uniqsite[rand_number]
            rand_reads_count[rand_reads] = rand_reads_count[rand_reads] + 1.0
        smoothed_result = correlate(array(rand_reads_count), kernel_score)
        scores = list()
        rand_mean = smoothed_result.mean()
        rand_std = smoothed_result.std()
        total_sum = smoothed_result.sum()
        rand_threshhold = rand_mean + threshold * rand_std
        higher_count = 0
        for now_site in kdesmooth_result:
            if kdesmooth_result[now_site] > rand_threshhold:
                higher_count = higher_count + 1
        randscore = dict()
        randscore['rand_mean'] = rand_mean
        randscore['rand_std'] = rand_std
        return randscore
    except KeyboardInterrupt:
        raise KeyboardInterruptError()


def get_bpc(bamfile, hotspots, jobtype, filted_region, nthreads, maxinsert=
    100000):
    total_length = 0
    pars = list()
    for hotspot_now in hotspots:
        par = dict()
        par['bamfile'] = bamfile
        par['hotspot'] = hotspot_now
        par['jobtype'] = jobtype
        par['filted_region'] = filted_region
        par['maxinsert'] = maxinsert
        pars.append(par)
        total_length = hotspot_now.end - hotspot_now.start + 1 + total_length
    pool = Pool(nthreads)
    try:
        reads_count = pool.map(bpc_runner, pars)
        total_reads = 0.0
        for count_now in reads_count:
            total_reads = total_reads + count_now
        bpc = (total_reads + 0.0) / total_length
        pool.close()
        return bpc
    except KeyboardInterrupt:
        pool.terminate()
        print('You cancelled the program!')
        sys.exit(1)
    except Exception as e:
        print('got exception: %r, terminating the pool' % (e,))
        pool.terminate()
        print('pool is terminated')
    finally:
        pool.join()


def bpc_runner(par):
    try:
        bamfile = par['bamfile']
        hotspot = par['hotspot']
        jobtype = par['jobtype']
        filted_region = par['filted_region']
        maxinsert = par['maxinsert']
        start_site = hotspot.start
        end_site = hotspot.end
        whether_in_fr = 0
        chromosome = hotspot.chromosome
        hotspotregio = chromosome + ':' + str(start_site) + '-' + str(end_site)
        hotspotreads = 0
        for i in range(start_site, end_site + 1):
            parentscare = int(i / 100)
            if chromosome in filted_region:
                if parentscare in filted_region[chromosome]:
                    whether_in_fr = 1
        if whether_in_fr == 0:
            readscount = dict()
            if jobtype == 'nhsingle':
                readscount = readscounter.nhreadscounter(bamfile=bamfile,
                    region=hotspotregio, paired=False, maxinsert=maxinsert)
            elif jobtype == 'nhpaired':
                readscount = readscounter.nhreadscounter(bamfile=bamfile,
                    region=hotspotregio, paired=True, maxinsert=maxinsert)
            elif jobtype == 'dh':
                readscount = readscounter.dhreadscounter(bamfile=bamfile,
                    region=hotspotregio)
            else:
                print('%s count type error!!!!' % jobtype)
                sys.exit(1)
            for i in readscount:
                hotspotreads = hotspotreads + readscount[i]
        return hotspotreads
    except KeyboardInterrupt:
        raise KeyboardInterruptError()

import pysam
from numpy import *
from multiprocessing import Pool
import random as rnd
from .kernel import *
import sys
from . import readscounter


class KeyboardInterruptError(Exception):
    pass


def nhnoncontrol(uniqueratio, threshold, kernellength, nthreads=4):
    bgscore = sim_replicate_nthreads(run_times=200, uniqueratio=uniqueratio,
        nthreads=nthreads, kernellength=kernellength, threshold=threshold)
    cutoff = bgscore['mean'] + bgscore['std'] * threshold
    return cutoff


def dhnoncontrol(uniqueratio, threshold, kernellength, nthreads=4):
    bgscore = sim_replicate_nthreads(run_times=200, uniqueratio=uniqueratio,
        nthreads=nthreads, kernellength=kernellength, threshold=threshold)
    cutoff = bgscore['mean'] + bgscore['std'] * threshold
    return cutoff


def nhcontrol(bamfile, chromosome, paired, chrlength, ultratio,
    filted_region, maxinsert, kernellength=600, threshold=4):
    """
        region: chr:start-end
        ultraio = chrlength * uniqratio / chr_total_reads
        filter region
        
    """
    region = chromosome + ':' + str(1) + '-' + str(chrlength)
    readscount = readscounter.nhreadscounter(bamfile, region, paired,
        maxinsert=maxinsert)
    kernel = smooth_kernel(kernellength)
    kernel_score = list()
    for i in sorted(kernel):
        kernel_score.append(kernel[i])
    threshold = filted_region.threshold / 100
    for site in readscount:
        if readscount[site] > threshold:
            readscount[site] = threshold
    smoothed_result = correlate(array(readscount), kernel_score)
    ultratiolist = list()
    ultratiolist.append(ultratio)
    smoothed_result = correlate(smoothed_result, ultratiolist)
    bg_mean = smoothed_result.mean()
    bg_std = smoothed_result.std()
    bg_threshold = bg_mean + threshold * bg_std
    cutoff = bg_threshold
    return cutoff


def nhuniquerate(bamfile, chromosome, paired, fregion, regionstart=1,
    regionend=-1, maxinsert=100000):
    samfile = pysam.Samfile(bamfile)
    ref_lengths = samfile.lengths
    sam_ref = samfile.references
    refere_ncenumber = samfile.nreferences
    if regionend == -1:
        for i in range(refere_ncenumber):
            if sam_ref[i] == chromosome:
                regionend = ref_lengths[i]
    region = chromosome + ':' + str(regionstart) + '-' + str(regionend)
    region_length = regionend - regionstart
    nhreadscount = readscounter.nhreadscounter(bamfile, region, paired,
        maxinsert=maxinsert)
    totaluniq = len(nhreadscount) + 0.0
    uniquerate = totaluniq / region_length
    return uniquerate


def dhuniquerate(bamfile, chromosome, regionstart=1, regionend=-1):
    """
        Count unique Rate in a region
    
    """
    samfile = pysam.Samfile(bamfile)
    ref_lengths = samfile.lengths
    sam_ref = samfile.references
    refere_ncenumber = samfile.nreferences
    if regionend == -1:
        for i in range(refere_ncenumber):
            if sam_ref[i] == chromosome:
                regionend = ref_lengths[i]
    region = chromosome + ':' + str(regionstart) + '-' + str(regionend)
    region_length = regionend - regionstart
    dhreadscount = readscounter.dhreadscounter(bamfile, region)
    totaluniq = len(dhreadscount) + 0.0
    uniquerate = totaluniq / region_length
    return uniquerate


def ultratio(chrlength, uniqueratio, chrtotalreads, frcount):
    """
    ultraio = chrlength * uniqueratio / chr_total_reads
    """
    ultratio = chrlength * uniqueratio / (chrtotalreads - frcount)
    return ultratio


def sim_replicate_nthreads(run_times=1000, uniqueratio=1, kernellength=600,
    threshold=4, nthreads=2):
    pars = list()
    for i in range(0, run_times):
        par = dict()
        par['uniqueratio'] = uniqueratio
        par['kernellength'] = kernellength
        par['threshold'] = threshold
        pars.append(par)
    pool = Pool(nthreads)
    outscore = dict()
    try:
        randomthresh = pool.map(sim_bg_thread_worker, pars)
        summean = 0.0
        sumstd = 0.0
        for randscore in randomthresh:
            randmean = randscore['rand_mean']
            randstd = randscore['rand_std']
            summean = summean + randmean
            sumstd = sumstd + randstd
        mean_of_mean = summean / run_times
        mean_of_std = sumstd / run_times
        outscore['mean'] = mean_of_mean
        outscore['std'] = mean_of_std
        pool.close()
        return outscore
    except KeyboardInterrupt:
        pool.terminate()
        print('You cancelled the program!')
        sys.exit(1)
    except Exception as e:
        print('got exception: %r, terminating the pool' % (e,))
        pool.terminate()
        print('pool is terminated')
    finally:
        pool.join()


def sim_bg_thread_worker(par):
    try:
        uniqueratio = par['uniqueratio']
        kernellength = par['kernellength']
        threshold = par['threshold']
        kernel = smooth_kernel(length=kernellength)
        sim_genome_size = int(100000.0)
        total_reads = int(sim_genome_size * uniqueratio)
        region_site = list(range(0, sim_genome_size))
        sim_uniqsite = rnd.sample(region_site, total_reads)
        rand_reads_count = list()
        for i in range(0, sim_genome_size):
            rand_reads_count.append(0)
        kernel_score = list()
        for i in sorted(kernel):
            kernel_score.append(kernel[i])
        kdesmooth_result = dict()
        for i in range(0, total_reads):
            rand_number = int(rnd.uniform(0, total_reads))
            rand_reads = sim_uniqsite[rand_number]
            rand_reads_count[rand_reads] = rand_reads_count[rand_reads] + 1.0
        smoothed_result = correlate(array(rand_reads_count), kernel_score)
        scores = list()
        rand_mean = smoothed_result.mean()
        rand_std = smoothed_result.std()
        total_sum = smoothed_result.sum()
        rand_threshhold = rand_mean + threshold * rand_std
        higher_count = 0
        for now_site in kdesmooth_result:
            if kdesmooth_result[now_site] > rand_threshhold:
                higher_count = higher_count + 1
        randscore = dict()
        randscore['rand_mean'] = rand_mean
        randscore['rand_std'] = rand_std
        return randscore
    except KeyboardInterrupt:
        raise KeyboardInterruptError()


def get_bpc(bamfile, hotspots, jobtype, filted_region, nthreads, maxinsert=
    100000):
    total_length = 0
    pars = list()
    for hotspot_now in hotspots:
        par = dict()
        par['bamfile'] = bamfile
        par['hotspot'] = hotspot_now
        par['jobtype'] = jobtype
        par['filted_region'] = filted_region
        par['maxinsert'] = maxinsert
        pars.append(par)
        total_length = hotspot_now.end - hotspot_now.start + 1 + total_length
    pool = Pool(nthreads)
    try:
        reads_count = pool.map(bpc_runner, pars)
        total_reads = 0.0
        for count_now in reads_count:
            total_reads = total_reads + count_now
        bpc = (total_reads + 0.0) / total_length
        pool.close()
        return bpc
    except KeyboardInterrupt:
        pool.terminate()
        print('You cancelled the program!')
        sys.exit(1)
    except Exception as e:
        print('got exception: %r, terminating the pool' % (e,))
        pool.terminate()
        print('pool is terminated')
    finally:
        pool.join()


def bpc_runner(par):
    try:
        bamfile = par['bamfile']
        hotspot = par['hotspot']
        jobtype = par['jobtype']
        filted_region = par['filted_region']
        maxinsert = par['maxinsert']
        start_site = hotspot.start
        end_site = hotspot.end
        whether_in_fr = 0
        chromosome = hotspot.chromosome
        hotspotregio = chromosome + ':' + str(start_site) + '-' + str(end_site)
        hotspotreads = 0
        for i in range(start_site, end_site + 1):
            parentscare = int(i / 100)
            if chromosome in filted_region:
                if parentscare in filted_region[chromosome]:
                    whether_in_fr = 1
        if whether_in_fr == 0:
            readscount = dict()
            if jobtype == 'nhsingle':
                readscount = readscounter.nhreadscounter(bamfile=bamfile,
                    region=hotspotregio, paired=False, maxinsert=maxinsert)
            elif jobtype == 'nhpaired':
                readscount = readscounter.nhreadscounter(bamfile=bamfile,
                    region=hotspotregio, paired=True, maxinsert=maxinsert)
            elif jobtype == 'dh':
                readscount = readscounter.dhreadscounter(bamfile=bamfile,
                    region=hotspotregio)
            else:
                print('%s count type error!!!!' % jobtype)
                sys.exit(1)
            for i in readscount:
                hotspotreads = hotspotreads + readscount[i]
        return hotspotreads
    except KeyboardInterrupt:
        raise KeyboardInterruptError()



import pysam
from numpy import *
from multiprocessing import Pool
import random as rnd
from .kernel import *
import sys
from . import readscounter


class KeyboardInterruptError(Exception):
    
    pass


def nhnoncontrol(uniqueratio, threshold, kernellength, nthreads=4):
    
    bgscore = sim_replicate_nthreads(run_times=200, uniqueratio=uniqueratio,
                                     nthreads=nthreads, kernellength=kernellength, threshold=threshold)
    
    cutoff = bgscore['mean'] + bgscore['std'] * threshold
    
    return cutoff


def dhnoncontrol(uniqueratio, threshold, kernellength, nthreads=4):
    
    bgscore = sim_replicate_nthreads(run_times=200, uniqueratio=uniqueratio,
                                     nthreads=nthreads, kernellength=kernellength, threshold=threshold)
    
    cutoff = bgscore['mean'] + bgscore['std'] * threshold

    return cutoff


def nhcontrol(bamfile, chromosome, paired, chrlength, ultratio, filted_region,maxinsert, kernellength = 600, threshold = 4):
    
    """
        region: chr:start-end
        ultraio = chrlength * uniqratio / chr_total_reads
        filter region
        
    """
    
    region = chromosome + ':' + str(1) + '-' + str(chrlength)
    
    readscount = readscounter.nhreadscounter(bamfile, region, paired,  maxinsert=maxinsert)
    
    kernel = smooth_kernel(kernellength)
    
    kernel_score = list()

    for i in sorted(kernel):
        
        kernel_score.append(kernel[i])

    threshold = filted_region.threshold/100

    for site in readscount:
        
        if readscount[site] > threshold:
            
            readscount[site] = threshold

    smoothed_result = correlate(array(readscount), kernel_score)

    ultratiolist = list()

    ultratiolist.append(ultratio)

    smoothed_result = correlate(smoothed_result, ultratiolist)

    #scores = list()
    
    bg_mean = smoothed_result.mean()
    
    bg_std = smoothed_result.std()
    
    bg_threshold = bg_mean + threshold * bg_std
    
    #bgscore['rand_mean'] = bg_mean
    
    #bgscore['rand_std'] = bg_std
    
    cutoff = bg_threshold

    return cutoff


def nhuniquerate(bamfile, chromosome, paired, fregion, regionstart=1, regionend = -1, maxinsert = 100000):
    
    samfile = pysam.Samfile(bamfile)

    ref_lengths = samfile.lengths

    sam_ref = samfile.references

    refere_ncenumber = samfile.nreferences

    if regionend == -1:
        
        for i in range(refere_ncenumber):
            
            if sam_ref[i] == chromosome:
                
                regionend = ref_lengths[i]

    region = chromosome + ':' + str(regionstart) + '-' + str(regionend)
    
    region_length = regionend - regionstart
        
    nhreadscount = readscounter.nhreadscounter(bamfile, region, paired, maxinsert=maxinsert)

    totaluniq = len(nhreadscount)  + 0.0
    
    uniquerate = totaluniq/region_length
    
    return uniquerate


def dhuniquerate(bamfile, chromosome, regionstart=1, regionend=-1):
    
    """
        Count unique Rate in a region
    
    """
    
    samfile = pysam.Samfile(bamfile)

    ref_lengths = samfile.lengths

    sam_ref = samfile.references

    refere_ncenumber = samfile.nreferences

    if regionend == -1:
        
        for i in range(refere_ncenumber):
            
            if sam_ref[i] == chromosome:
                
                regionend = ref_lengths[i]

    region = chromosome + ':' + str(regionstart) + '-' + str(regionend)
    
    region_length = regionend - regionstart
    
    dhreadscount = readscounter.dhreadscounter(bamfile, region)
    
    totaluniq = len(dhreadscount) + 0.0
    
    uniquerate = totaluniq/region_length
    
    return uniquerate


def ultratio(chrlength, uniqueratio, chrtotalreads, frcount):
    """
    ultraio = chrlength * uniqueratio / chr_total_reads
    """
    ultratio = chrlength * uniqueratio / (chrtotalreads - frcount)

    return ultratio


def sim_replicate_nthreads(run_times=1000, uniqueratio=1, kernellength = 600, threshold = 4, nthreads = 2):
    # randomthresh = list()
    
    pars = list()
    
    for i in range(0,run_times):
        
        par=dict()
        
        par['uniqueratio'] = uniqueratio
        
        par['kernellength'] = kernellength
        
        par['threshold'] = threshold
        
        pars.append(par)
    
    pool=Pool(nthreads)
    
    outscore = dict()
    
    try:
        randomthresh = pool.map(sim_bg_thread_worker, pars)

        summean = 0.0
    
        sumstd = 0.0

        for randscore in randomthresh:

            randmean = randscore['rand_mean']
            
            randstd = randscore['rand_std']
            # print (randmean, randstd)
            summean = summean + randmean
            
            sumstd = sumstd + randstd

        mean_of_mean = summean/run_times
        
        mean_of_std = sumstd/run_times
        # print ('mean_of_mean',mean_of_mean, 'mean_of_std',mean_of_std)

        outscore['mean'] = mean_of_mean
        
        outscore['std'] = mean_of_std
        #return (mean_of_mean, mean_of_std)

        pool.close()
        
        return outscore

    except KeyboardInterrupt:
        
        pool.terminate()
        
        print ("You cancelled the program!")
        
        sys.exit(1)
        
    except Exception as e:
        
        print ('got exception: %r, terminating the pool' % (e,))
        
        pool.terminate()
        
        print ('pool is terminated')
        
    finally:
        # print ('joining pool processes')
        pool.join()
    # print ('join complete')
    # pool.join()
    # pool.close()


def sim_bg_thread_worker(par):

    try:
        
        uniqueratio=par['uniqueratio']
        
        kernellength = par['kernellength']
        
        threshold = par['threshold']

        kernel = smooth_kernel(length=kernellength)

        sim_genome_size = int(1e5)

        total_reads = int(sim_genome_size * uniqueratio)

        region_site = list(range(0,sim_genome_size))

        sim_uniqsite = rnd.sample(region_site, total_reads)

        rand_reads_count = list()

        for i in range(0,sim_genome_size):
            
            rand_reads_count.append(0)

        kernel_score = list()

        for i in sorted(kernel):
            kernel_score.append(kernel[i])



        kdesmooth_result = dict()

        for i in range(0,total_reads):
            
            rand_number = int(rnd.uniform(0,total_reads))
            
            rand_reads = sim_uniqsite[rand_number]

            rand_reads_count[rand_reads] = rand_reads_count[rand_reads] + 1.0

        smoothed_result = correlate(array(rand_reads_count), kernel_score)

        scores = list()

        rand_mean = smoothed_result.mean()
        
        rand_std = smoothed_result.std()
        
        total_sum = smoothed_result.sum()
        
        rand_threshhold = rand_mean + threshold * rand_std

        higher_count = 0

        for now_site in kdesmooth_result:
            
            if kdesmooth_result[now_site] > rand_threshhold:
                
                higher_count = higher_count + 1

        # print (total_sum, rand_mean, rand_std, rand_threshhold, higher_count, total_reads)

        randscore = dict()
        
        randscore['rand_mean'] = rand_mean
        
        randscore['rand_std'] = rand_std

        return randscore

    except KeyboardInterrupt:
        
        raise KeyboardInterruptError()


def get_bpc(bamfile,  hotspots, jobtype, filted_region, nthreads, maxinsert = 100000):

    #bpc average readscount per basepare

    # total_reads = 0
    #
    total_length = 0
    #
    # samfile = pysam.Samfile(bamfile)

    pars = list()

    for hotspot_now in hotspots:

        par = dict()

        par['bamfile'] = bamfile

        par['hotspot'] = hotspot_now

        par['jobtype'] = jobtype

        par['filted_region'] = filted_region

        par['maxinsert'] = maxinsert

        pars.append(par)

        total_length = hotspot_now.end - hotspot_now.start + 1 + total_length


    # print ("total length %s" % total_length)

    pool = Pool(nthreads)

    try:

        reads_count = pool.map(bpc_runner, pars)

        total_reads = 0.0

        for count_now in reads_count:

            total_reads = total_reads + count_now

        bpc = (total_reads+0.0)/total_length

        pool.close()

        return bpc

    except KeyboardInterrupt:

        pool.terminate()

        print ("You cancelled the program!")

        sys.exit(1)

    except Exception as e:

        print ('got exception: %r, terminating the pool' % (e,))

        pool.terminate()

        print ('pool is terminated')

    finally:
        # print ('joining pool processes')
        pool.join()
        # print ('join complete')


def bpc_runner(par):

    try:

        bamfile = par['bamfile']

        hotspot = par['hotspot']

        jobtype = par['jobtype']

        filted_region = par['filted_region']

        maxinsert = par['maxinsert']

        start_site = hotspot.start

        end_site = hotspot.end

        whether_in_fr = 0

        chromosome = hotspot.chromosome

        hotspotregio = chromosome + ':' + str(start_site) + '-' + str(end_site)

        hotspotreads = 0

        for i in range(start_site, end_site + 1):

            parentscare = int(i/100)

            if chromosome in filted_region:

                if parentscare in filted_region[chromosome]:

                    whether_in_fr = 1

        if whether_in_fr == 0:

            readscount = dict()

            if jobtype == 'nhsingle':

                readscount = readscounter.nhreadscounter(bamfile = bamfile, region=hotspotregio, paired=False, maxinsert=maxinsert)

            elif jobtype == 'nhpaired':

                readscount = readscounter.nhreadscounter(bamfile = bamfile, region=hotspotregio, paired=True, maxinsert=maxinsert)

            elif jobtype == 'dh':

                readscount = readscounter.dhreadscounter(bamfile = bamfile, region = hotspotregio)

            else:

                print ("%s count type error!!!!" % jobtype)

                sys.exit(1)

            for i in readscount:

                hotspotreads = hotspotreads+readscount[i]

        # print (hotspotreads)
        return hotspotreads

    except KeyboardInterrupt:

        raise KeyboardInterruptError()