"""Common full-read scoring for native and graph-proposed fragment placements.

Candidate retrieval scores are not used as likelihoods. Each mate is realigned
with edlib and scored from its original bases/qualities. Scores are approximate
alignment likelihood ratios (not calibrated genotype confidence). The same
procedure applies to IPD, observed paths and competing HLA loci.
"""
import math
import re

COMPLEMENT = str.maketrans('ACGTN', 'TGCAN')


def reverse_complement(sequence):
    return sequence.translate(COMPLEMENT)[::-1]


def alignment_score(query, qualities, target, cigar):
    """Likelihood ratio to iid uniform bases; affine indel error approximation."""
    if len(query) != len(qualities):
        raise ValueError('Quality length mismatch')
    operations = re.findall(r'(\d+)([=XID])', cigar)
    if not operations or ''.join(n + op for n, op in operations) != cigar:
        raise ValueError('Invalid extended CIGAR')
    i = j = 0
    score = 0.0
    for count, op in operations:
        count = int(count)
        if count < 1:
            raise ValueError('Empty operation')
        if op in '=X':
            if i + count > len(query) or j + count > len(target):
                raise ValueError('CIGAR exceeds sequence')
            for offset in range(count):
                a, b = query[i + offset], target[j + offset]
                if a == 'N' or b == 'N':
                    continue
                if (a == b) != (op == '='):
                    raise ValueError('CIGAR differs from bases')
                q = min(40, max(3, ord(qualities[i + offset]) - 33))
                error = 10 ** (-q / 10)
                score += math.log(4 * ((1 - error) if a == b else error / 3))
            i += count
            j += count
        else:
            score += math.log(.001) + (count - 1) * math.log(.1)
            if op == 'I':
                i += count
            else:
                j += count
    if i != len(query) or j != len(target):
        raise ValueError('Alignment did not consume its complete sequences')
    return score


def gene_penalty(query, qualities, target, cigar, offset, body):
    """Within-gene error penalty relative to a perfect match, without length bonus.

    This separates allele-discriminating gene evidence from a linked flank's
    mapping advantage. A perfectly matching conserved gene segment scores zero
    for every allele; extra flanking matches cannot create gene discrimination.
    """
    i = j = covered = 0
    score = 0.0
    for count, op in re.findall(r'(\d+)([=XID])', cigar):
        count = int(count)
        if op in '=X':
            for n in range(count):
                if body[0] <= offset + j + n < body[1]:
                    a, b = query[i + n], target[j + n]
                    if a == 'N' or b == 'N':
                        continue
                    covered += 1
                    if a != b:
                        q = min(40, max(3, ord(qualities[i + n]) - 33))
                        error = 10 ** (-q / 10)
                        score += math.log((error / 3) / (1 - error))
            i += count
            j += count
        else:
            overlap = (body[0] <= offset + j < body[1]) if op == 'I' else (
                max(body[0], offset + j) < min(body[1], offset + j + count))
            if overlap:
                score += math.log(.001) + (count - 1) * math.log(.1)
            if op == 'I':
                i += count
            else:
                j += count
    return score, covered


def realign(query, qualities, reference, start, end, reverse=False, padding=12, body=None):
    import edlib
    if not query or not 0 <= start < end <= len(reference):
        raise ValueError('Invalid alignment window')
    if reverse:
        query, qualities = reverse_complement(query), qualities[::-1]
    low, high = max(0, start - padding), min(len(reference), end + padding)
    window = reference[low:high]
    result = edlib.align(query, window, mode='HW', task='path',
                         k=max(8, len(query) // 5))
    if result['editDistance'] < 0:
        return None
    # edlib's CIGAR corresponds to its first optimal location. Alternative seed
    # windows are retained by the caller; this is not exhaustive alignment.
    begin, last = result['locations'][0]
    if begin is None or last is None or not result.get('cigar'):
        return None
    result = dict(start=low + begin, end=low + last + 1, reverse=reverse,
                score=alignment_score(query, qualities, window[begin:last + 1], result['cigar']),
                edit_distance=result['editDistance'], cigar=result['cigar'])
    if body is not None:
        penalty, covered = gene_penalty(query, qualities, window[begin:last + 1], result['cigar'], low + begin, body)
        result.update(gene_penalty=penalty, gene_bases=covered)
    return result


def pair_score(first, second, mean, stdev, maximum_insert):
    if mean <= 0 or stdev <= 0 or maximum_insert <= 0:
        raise ValueError('Invalid insert parameters')
    if first is None or second is None or first['reverse'] == second['reverse']:
        return None
    forward, reverse = (second, first) if first['reverse'] else (first, second)
    span = reverse['end'] - forward['start']
    if not (forward['start'] <= reverse['start'] and forward['end'] <= reverse['end']
            and 0 < span <= maximum_insert):
        return None
    return first['score'] + second['score'] - .5 * ((span - mean) / stdev) ** 2


def overlaps_gene(placement, body):
    """At least one sequenced base overlaps the body; insert span alone is insufficient."""
    return any(max(placement[a], body[0]) < min(placement[b], body[1])
               for a, b in [('first_start', 'first_end'), ('second_start', 'second_end')])


def internal_pair(placement, body):
    return all(body[0] <= placement[a] < placement[b] <= body[1]
               for a, b in [('first_start', 'first_end'), ('second_start', 'second_end')])