# Path: covest/constants.py
# MAX_EXP = 200
# GRID_DEPTH = 3
# INF = float('inf')
# VERBOSE = True
# PLOT_LOG_SCALE = True
# USE_BIGFLOAT = False
# STEP = 1.1
# OPTIMIZATION_METHOD = 'L-BFGS-B'
# INITIAL_GRID_COUNT = 20
# INITIAL_GRID_STEP = 3
# DEFAULT_ERR_SCALE = 1
# DEFAULT_K = 21
# DEFAULT_READ_LENGTH = 100
# DEFAULT_REPEAT_MODEL = 0
# DEFAULT_MIN_SINGLECOPY_RATIO = 0.3
# AUTO_SAMPLE_TARGET_COVERAGE = 12
# AUTO_TRIM_PRECISION = 6
# NOISE_THRESHOLD = 10**-6
# MAX_ERRORS = 8
# MAX_NOTRIM = 25
#     DEFAULT_THREAD_COUNT = cpu_count()
#     DEFAULT_THREAD_COUNT = 2
# 
# Path: covest/utils.py
# def estimate_p(cc, alpha):
#     return (cc * (alpha - 1)) / (alpha * cc - alpha - cc)
# 
# def kmer_to_read_coverage(coverage, k, r):
#     return coverage * r / (r - k + 1)
# 
# def fix_coverage(coverage):
#     return inverse(lambda c: (c - c * exp(-c)) / (1 - exp(-c) - c * exp(-c)))(coverage)
# 
# def verbose_print(message):
#     if not constants.VERBOSE:
#         return
#     sys.stderr.write(message + "\n")


import random
from collections import defaultdict
from covest_poisson import poisson_dist
from math import exp, floor, ceil
from scipy.stats import binom
from covest import constants
from .utils import estimate_p, kmer_to_read_coverage, fix_coverage, verbose_print




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

# Path: covest/constants.py
# MAX_EXP = 200
# GRID_DEPTH = 3
# INF = float('inf')
# VERBOSE = True
# PLOT_LOG_SCALE = True
# USE_BIGFLOAT = False
# STEP = 1.1
# OPTIMIZATION_METHOD = 'L-BFGS-B'
# INITIAL_GRID_COUNT = 20
# INITIAL_GRID_STEP = 3
# DEFAULT_ERR_SCALE = 1
# DEFAULT_K = 21
# DEFAULT_READ_LENGTH = 100
# DEFAULT_REPEAT_MODEL = 0
# DEFAULT_MIN_SINGLECOPY_RATIO = 0.3
# AUTO_SAMPLE_TARGET_COVERAGE = 12
# AUTO_TRIM_PRECISION = 6
# NOISE_THRESHOLD = 10**-6
# MAX_ERRORS = 8
# MAX_NOTRIM = 25
#     DEFAULT_THREAD_COUNT = cpu_count()
#     DEFAULT_THREAD_COUNT = 2
# 
# Path: covest/utils.py
# def estimate_p(cc, alpha):
#     return (cc * (alpha - 1)) / (alpha * cc - alpha - cc)
# 
# def kmer_to_read_coverage(coverage, k, r):
#     return coverage * r / (r - k + 1)
# 
# def fix_coverage(coverage):
#     return inverse(lambda c: (c - c * exp(-c)) / (1 - exp(-c) - c * exp(-c)))(coverage)
# 
# def verbose_print(message):
#     if not constants.VERBOSE:
#         return
#     sys.stderr.write(message + "\n")


# Path: covest/histogram.py
import random
from collections import defaultdict
from covest_poisson import poisson_dist
from math import exp, floor, ceil
from scipy.stats import binom
from covest import constants
from .utils import estimate_p, kmer_to_read_coverage, fix_coverage, verbose_print





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))