# Path: covest/histogram.py
# def compute_coverage_apx(hist, k, r):
#     observed_ones = hist.get(1, 0)
#     all_kmers = sum(i * h for i, h in hist.items())
#     total_unique_kmers = sum(h for h in hist.values())
# 
#     if total_unique_kmers == 0:
#         return 0.0, 1.0
# 
#     # discard first column
#     all_kmers -= observed_ones
#     unique_kmers = total_unique_kmers - observed_ones
#     # compute coverage from hist >=2
#     try:
#         cov = all_kmers / unique_kmers
#         cov = fix_coverage(cov)
#         # fix unique kmers
#         unique_kmers /= (1.0 - exp(-cov) - cov * exp(-cov))
#         # compute alpha (error read ratio)
#         estimated_ones = unique_kmers * cov * exp(-cov)
#         estimated_zeros = unique_kmers * exp(-cov)
#         error_ones = max(0.0, observed_ones - estimated_ones)
#         alpha = error_ones / (total_unique_kmers + estimated_zeros)
#         # estimate probability of correct kmer and error rate
#         estimated_p = max(0.0, estimate_p(cov, alpha))
#         e = 1 - estimated_p ** (1.0 / k)
#         # return corrected coverage and error estimate
#         if estimated_p > 0:
#             # function for conversion between kmer and base coverage
#             return float(kmer_to_read_coverage(cov / estimated_p, k, r)), float(e)
#         else:
#             return 0.0, float(e)
#     except ZeroDivisionError:
#         return 0.0, 1.0
# 
# def sample_histogram(hist, factor=2, trim=None):
#     if trim is None:
#         if len(hist) > 300:
#             trim = get_trim(hist)
#         else:
#             trim = max(hist)
#     else:
#         trim = min(max(hist), trim * factor)
# 
#     hist = {k: v for k, v in hist.items() if k < trim}
#     h = defaultdict(int)
#     prob = 1.0 / factor
# 
#     for i, v in hist.items():
#         if i < 100:
#             b = binom(i, prob)
#             probs = [b.pmf(j) for j in range(1, i + 1)]
#         else:
#             probs = poisson_dist(i * prob, i)
# 
#         for j, p in enumerate(probs):
#             h[j + 1] += v * p
# 
#     h = dict(h)
#     for i, v in h.items():
#         d = v - round(v)
#         h[i] = ceil(v) if random.random() < d else floor(v)
#     return {k: v for k, v in h.items() if v > 0}   # remove 0 elements
# 
# def trim_hist(hist, threshold):
#     if threshold >= max(hist):
#         return hist, 0
#     h = {k: v for k, v in hist.items() if k < threshold}
#     tail = sum(v for k, v in hist.items() if k >= threshold)
#     # remove 0 elements
#     return {k: v for k, v in h.items() if v > 0}, tail


import unittest
    import sys
from covest.histogram import compute_coverage_apx, sample_histogram, trim_hist
        histogram = {
            1: 5311035,
            2: 310620,
            3: 223796,
            4: 150045,
            5: 81833,
            6: 37443,
            7: 14628,
            8: 5075,
            9: 1505,
            10: 386,
            11: 133,
            12: 22,
            13: 5,
            14: 9,
            15: 1,
        }
        true_hist = {
            1: 2949548,
            2: 255275,
            3: 108173,
            4: 36802,
            5: 10079,
            6: 2303,
            7: 452,
            8: 79,
            9: 12,
            10: 3,
            11: 1,
        }

# Path: covest/histogram.py
# def compute_coverage_apx(hist, k, r):
#     observed_ones = hist.get(1, 0)
#     all_kmers = sum(i * h for i, h in hist.items())
#     total_unique_kmers = sum(h for h in hist.values())
# 
#     if total_unique_kmers == 0:
#         return 0.0, 1.0
# 
#     # discard first column
#     all_kmers -= observed_ones
#     unique_kmers = total_unique_kmers - observed_ones
#     # compute coverage from hist >=2
#     try:
#         cov = all_kmers / unique_kmers
#         cov = fix_coverage(cov)
#         # fix unique kmers
#         unique_kmers /= (1.0 - exp(-cov) - cov * exp(-cov))
#         # compute alpha (error read ratio)
#         estimated_ones = unique_kmers * cov * exp(-cov)
#         estimated_zeros = unique_kmers * exp(-cov)
#         error_ones = max(0.0, observed_ones - estimated_ones)
#         alpha = error_ones / (total_unique_kmers + estimated_zeros)
#         # estimate probability of correct kmer and error rate
#         estimated_p = max(0.0, estimate_p(cov, alpha))
#         e = 1 - estimated_p ** (1.0 / k)
#         # return corrected coverage and error estimate
#         if estimated_p > 0:
#             # function for conversion between kmer and base coverage
#             return float(kmer_to_read_coverage(cov / estimated_p, k, r)), float(e)
#         else:
#             return 0.0, float(e)
#     except ZeroDivisionError:
#         return 0.0, 1.0
# 
# def sample_histogram(hist, factor=2, trim=None):
#     if trim is None:
#         if len(hist) > 300:
#             trim = get_trim(hist)
#         else:
#             trim = max(hist)
#     else:
#         trim = min(max(hist), trim * factor)
# 
#     hist = {k: v for k, v in hist.items() if k < trim}
#     h = defaultdict(int)
#     prob = 1.0 / factor
# 
#     for i, v in hist.items():
#         if i < 100:
#             b = binom(i, prob)
#             probs = [b.pmf(j) for j in range(1, i + 1)]
#         else:
#             probs = poisson_dist(i * prob, i)
# 
#         for j, p in enumerate(probs):
#             h[j + 1] += v * p
# 
#     h = dict(h)
#     for i, v in h.items():
#         d = v - round(v)
#         h[i] = ceil(v) if random.random() < d else floor(v)
#     return {k: v for k, v in h.items() if v > 0}   # remove 0 elements
# 
# def trim_hist(hist, threshold):
#     if threshold >= max(hist):
#         return hist, 0
#     h = {k: v for k, v in hist.items() if k < threshold}
#     tail = sum(v for k, v in hist.items() if k >= threshold)
#     # remove 0 elements
#     return {k: v for k, v in h.items() if v > 0}, tail


# Path: tests/test_histograms.py
import unittest
    import sys
from covest.histogram import compute_coverage_apx, sample_histogram, trim_hist

        histogram = {
            1: 5311035,
            2: 310620,
            3: 223796,
            4: 150045,
            5: 81833,
            6: 37443,
            7: 14628,
            8: 5075,
            9: 1505,
            10: 386,
            11: 133,
            12: 22,
            13: 5,
            14: 9,
            15: 1,
        }
        true_hist = {
            1: 2949548,
            2: 255275,
            3: 108173,
            4: 36802,
            5: 10079,
            6: 2303,
            7: 452,
            8: 79,
            9: 12,
            10: 3,
            11: 1,
        }

        sampled_hist = sample_histogram(histogram, 2)