def analyze(self, ploidy=2, K=23, covmax=1000000):
    '\n        Analyze Kmer spectrum, calculations derived from\n        allpathslg/src/kmers/KmerSpectra.cc\n        '
    from math import sqrt
    data = self.data
    kf_ceil = max((K for (K, c) in data))
    if (kf_ceil > covmax):
        exceeds = sum((1 for (K, c) in data if (K > covmax)))
        logging.debug('A total of {0} distinct K-mers appear > {1} times. Ignored ...'.format(exceeds, covmax))
        kf_ceil = covmax
    nkf = (kf_ceil + 1)
    a = ([0] * nkf)
    for (kf, c) in data:
        if (kf > kf_ceil):
            continue
        a[kf] = c
    ndk = a
    nk = [(k * c) for (k, c) in enumerate(a)]
    cndk = ([0] * nkf)
    cnk = ([0] * nkf)
    for kf in range(1, nkf):
        cndk[kf] = (cndk[(kf - 1)] + (0.5 * (ndk[(kf - 1)] + ndk[kf])))
        cnk[kf] = (cnk[(kf - 1)] + (0.5 * (nk[(kf - 1)] + nk[kf])))
    _kf_min1 = 10
    while (((_kf_min1 - 1) >= 2) and (nk[(_kf_min1 - 1)] < nk[_kf_min1])):
        _kf_min1 -= 1
    while ((_kf_min1 <= kf_ceil) and (nk[(_kf_min1 + 1)] < nk[_kf_min1])):
        _kf_min1 += 1
    _kf_max2 = _kf_min1
    for kf in range((_kf_min1 + 1), int((0.8 * kf_ceil))):
        if (nk[kf] > nk[_kf_max2]):
            _kf_max2 = kf
    if (ploidy == 2):
        ndk_half = ndk[(_kf_max2 / 2)]
        ndk_double = ndk[(_kf_max2 * 2)]
        if (ndk_double > ndk_half):
            _kf_max2 *= 2
    _kf_max1 = (_kf_max2 / 2)
    _kf_min2 = ((_kf_max1 * ((2 * ndk[_kf_max1]) + ndk[_kf_max2])) / (ndk[_kf_max1] + ndk[_kf_max2]))
    for kf in range((_kf_min1 + 1), _kf_max1):
        if (nk[kf] < nk[_kf_min1]):
            _kf_min1 = kf
    _kf_min3 = ((_kf_max2 * 3) / 2)
    print('kfs:', _kf_min1, _kf_max1, _kf_min2, _kf_max2, _kf_min3, file=sys.stderr)
    self.min1 = _kf_min1
    self.max1 = _kf_max1
    self.min2 = _kf_min2
    self.max2 = _kf_max2
    self.min3 = _kf_min3
    _kf_hi = ((_kf_max2 * sqrt(((4 * ndk[(2 * _kf_max2)]) * _kf_max2))) if ((2 * _kf_max2) < len(ndk)) else (_kf_max2 * sqrt(((4 * ndk[(len(ndk) - 1)]) * _kf_max2))))
    _kf_hi = int(_kf_hi)
    if (_kf_hi > kf_ceil):
        _kf_hi = kf_ceil
    _nk_total = cnk[(len(cnk) - 1)]
    _nk_bad_low_kf = cnk[_kf_min1]
    _nk_good_uniq = (cnk[_kf_min3] - cnk[_kf_min2])
    _nk_bad_high_kf = (_nk_total - cnk[_kf_hi])
    _ndk_good_snp = (cndk[_kf_min2] - cndk[_kf_min1])
    _ndk_good_uniq = (cndk[_kf_min3] - cndk[_kf_min2])
    _kf_ave_uniq = ((_nk_good_uniq * 1.0) / _ndk_good_uniq)
    _genome_size = (((_nk_total - _nk_bad_low_kf) - _nk_bad_high_kf) / _kf_ave_uniq)
    _genome_size_unique = (_ndk_good_uniq + (_ndk_good_snp / 2))
    _genome_size_repetitive = (_genome_size - _genome_size_unique)
    _coverage = ((_nk_total / _genome_size) if _genome_size else 0)
    if (ploidy == 2):
        _d_SNP = ((1.0 / (1.0 - ((1.0 - ((0.5 * _ndk_good_snp) / _genome_size)) ** (1.0 / K)))) if (_ndk_good_snp > 0) else 1000000)
    G = int(_genome_size)
    G1 = int(_genome_size_unique)
    GR = int(_genome_size_repetitive)
    coverage = int(_coverage)
    m = 'Kmer (K={0}) Spectrum Analysis\n'.format(K)
    m += 'Genome size estimate = {0}\n'.format(thousands(G))
    m += 'Genome size estimate CN = 1 = {0} ({1})\n'.format(thousands(G1), percentage(G1, G))
    m += 'Genome size estimate CN > 1 = {0} ({1})\n'.format(thousands(GR), percentage(GR, G))
    m += 'Coverage estimate: {0} x\n'.format(coverage)
    self.repetitive = 'Repeats: {0} percent'.format(((GR * 100) / G))
    if (ploidy == 2):
        d_SNP = int(_d_SNP)
        self.snprate = 'SNP rate ~= 1/{0}'.format(d_SNP)
    else:
        self.snprate = 'SNP rate not computed (Ploidy = {0})'.format(ploidy)
    m += (self.snprate + '\n')
    self.genomesize = int(round(((self.totalKmers * 1.0) / self.max2)))
    print(m, file=sys.stderr)
Analyze Kmer spectrum, calculations derived from
allpathslg/src/kmers/KmerSpectra.cc
def analyze(self, ploidy=2, K=23, covmax=1000000):
    '\n        Analyze Kmer spectrum, calculations derived from\n        allpathslg/src/kmers/KmerSpectra.cc\n        '
    from math import sqrt
    data = self.data
    kf_ceil = max((K for (K, c) in data))
    if (kf_ceil > covmax):
        exceeds = sum((1 for (K, c) in data if (K > covmax)))
        logging.debug('A total of {0} distinct K-mers appear > {1} times. Ignored ...'.format(exceeds, covmax))
        kf_ceil = covmax
    nkf = (kf_ceil + 1)
    a = ([0] * nkf)
    for (kf, c) in data:
        if (kf > kf_ceil):
            continue
        a[kf] = c
    ndk = a
    nk = [(k * c) for (k, c) in enumerate(a)]
    cndk = ([0] * nkf)
    cnk = ([0] * nkf)
    for kf in range(1, nkf):
        cndk[kf] = (cndk[(kf - 1)] + (0.5 * (ndk[(kf - 1)] + ndk[kf])))
        cnk[kf] = (cnk[(kf - 1)] + (0.5 * (nk[(kf - 1)] + nk[kf])))
    _kf_min1 = 10
    while (((_kf_min1 - 1) >= 2) and (nk[(_kf_min1 - 1)] < nk[_kf_min1])):
        _kf_min1 -= 1
    while ((_kf_min1 <= kf_ceil) and (nk[(_kf_min1 + 1)] < nk[_kf_min1])):
        _kf_min1 += 1
    _kf_max2 = _kf_min1
    for kf in range((_kf_min1 + 1), int((0.8 * kf_ceil))):
        if (nk[kf] > nk[_kf_max2]):
            _kf_max2 = kf
    if (ploidy == 2):
        ndk_half = ndk[(_kf_max2 / 2)]
        ndk_double = ndk[(_kf_max2 * 2)]
        if (ndk_double > ndk_half):
            _kf_max2 *= 2
    _kf_max1 = (_kf_max2 / 2)
    _kf_min2 = ((_kf_max1 * ((2 * ndk[_kf_max1]) + ndk[_kf_max2])) / (ndk[_kf_max1] + ndk[_kf_max2]))
    for kf in range((_kf_min1 + 1), _kf_max1):
        if (nk[kf] < nk[_kf_min1]):
            _kf_min1 = kf
    _kf_min3 = ((_kf_max2 * 3) / 2)
    print('kfs:', _kf_min1, _kf_max1, _kf_min2, _kf_max2, _kf_min3, file=sys.stderr)
    self.min1 = _kf_min1
    self.max1 = _kf_max1
    self.min2 = _kf_min2
    self.max2 = _kf_max2
    self.min3 = _kf_min3
    _kf_hi = ((_kf_max2 * sqrt(((4 * ndk[(2 * _kf_max2)]) * _kf_max2))) if ((2 * _kf_max2) < len(ndk)) else (_kf_max2 * sqrt(((4 * ndk[(len(ndk) - 1)]) * _kf_max2))))
    _kf_hi = int(_kf_hi)
    if (_kf_hi > kf_ceil):
        _kf_hi = kf_ceil
    _nk_total = cnk[(len(cnk) - 1)]
    _nk_bad_low_kf = cnk[_kf_min1]
    _nk_good_uniq = (cnk[_kf_min3] - cnk[_kf_min2])
    _nk_bad_high_kf = (_nk_total - cnk[_kf_hi])
    _ndk_good_snp = (cndk[_kf_min2] - cndk[_kf_min1])
    _ndk_good_uniq = (cndk[_kf_min3] - cndk[_kf_min2])
    _kf_ave_uniq = ((_nk_good_uniq * 1.0) / _ndk_good_uniq)
    _genome_size = (((_nk_total - _nk_bad_low_kf) - _nk_bad_high_kf) / _kf_ave_uniq)
    _genome_size_unique = (_ndk_good_uniq + (_ndk_good_snp / 2))
    _genome_size_repetitive = (_genome_size - _genome_size_unique)
    _coverage = ((_nk_total / _genome_size) if _genome_size else 0)
    if (ploidy == 2):
        _d_SNP = ((1.0 / (1.0 - ((1.0 - ((0.5 * _ndk_good_snp) / _genome_size)) ** (1.0 / K)))) if (_ndk_good_snp > 0) else 1000000)
    G = int(_genome_size)
    G1 = int(_genome_size_unique)
    GR = int(_genome_size_repetitive)
    coverage = int(_coverage)
    m = 'Kmer (K={0}) Spectrum Analysis\n'.format(K)
    m += 'Genome size estimate = {0}\n'.format(thousands(G))
    m += 'Genome size estimate CN = 1 = {0} ({1})\n'.format(thousands(G1), percentage(G1, G))
    m += 'Genome size estimate CN > 1 = {0} ({1})\n'.format(thousands(GR), percentage(GR, G))
    m += 'Coverage estimate: {0} x\n'.format(coverage)
    self.repetitive = 'Repeats: {0} percent'.format(((GR * 100) / G))
    if (ploidy == 2):
        d_SNP = int(_d_SNP)
        self.snprate = 'SNP rate ~= 1/{0}'.format(d_SNP)
    else:
        self.snprate = 'SNP rate not computed (Ploidy = {0})'.format(ploidy)
    m += (self.snprate + '\n')
    self.genomesize = int(round(((self.totalKmers * 1.0) / self.max2)))
    print(m, file=sys.stderr)