#!/usr/bin/env python3
"""
HopfDNA Fiber Analysis — Count exotic classes via S³ fiber structure

The 28 exotic spheres are NOT features of the S⁴ base.
They are TWISTS in the S³ FIBER over each base point.

For each S⁴ projection (base point), collect all DNA sequences
that map to it. Their fiber coordinates (the S³ = SU(2) element)
should show exactly 28 distinct twist patterns.

This is the correct topological approach: the 28 lives in the fiber.
"""

import math
import itertools
from collections import defaultdict

DNA_TO_Q = {
    'A': (1, 0, 0, 0),
    'C': (0, 1, 0, 0),
    'G': (0, 0, 1, 0),
    'T': (0, 0, 0, 1),
}

def quat_mul(a, b):
    a1, a2, a3, a4 = a
    b1, b2, b3, b4 = b
    return (
        a1*b1 - a2*b2 - a3*b3 - a4*b4,
        a1*b2 + a2*b1 + a3*b4 - a4*b3,
        a1*b3 - a2*b4 + a3*b1 + a4*b2,
        a1*b4 + a2*b3 - a3*b2 + a4*b1,
    )

def quat_conj(q):
    return (q[0], -q[1], -q[2], -q[3])

def quat_norm_sq(q):
    return q[0]**2 + q[1]**2 + q[2]**2 + q[3]**2

def quat_inv(q):
    ns = quat_norm_sq(q)
    if ns == 0: return None
    c = quat_conj(q)
    return tuple(x/ns for x in c)

def quat_div(a, b):
    bi = quat_inv(b)
    if bi is None: return None
    return quat_mul(a, bi)

def dna_to_s7(seq):
    def bases_to_quat(bases):
        q = [0.0, 0.0, 0.0, 0.0]
        for b in bases:
            v = DNA_TO_Q[b]
            q[0] += v[0]; q[1] += v[1]; q[2] += v[2]; q[3] += v[3]
        return tuple(q)
    q1 = bases_to_quat(seq[:4])
    q2 = bases_to_quat(seq[4:])
    n1 = quat_norm_sq(q1)
    n2 = quat_norm_sq(q2)
    total = n1 + n2
    if total == 0: return ((1,0,0,0), (0,0,0,1))
    scale = math.sqrt(total)
    return (tuple(x/scale for x in q1), tuple(x/scale for x in q2))

def hopf_map(s7_point):
    q1, q2 = s7_point
    return quat_div(q1, q2)

def fiber_coordinate(s7_point):
    """
    The FIBER coordinate of (q1, q2) ∈ S⁷ over base π(q1,q2) = q1/q2.
    
    Two points (q1,q2) and (q1',q2') are in the same fiber iff
    q1/q2 = q1'/q2', i.e., q1*q2⁻¹ = q1'*q2'⁻¹.
    
    The fiber element is λ = q2/|q2| ∈ S³ (the phase of q2).
    This is the SU(2) rotation that distinguishes exotic classes.
    """
    q1, q2 = s7_point
    n2 = math.sqrt(quat_norm_sq(q2))
    if n2 < 1e-15:
        return (1, 0, 0, 0)  # identity at infinity
    return tuple(x/n2 for x in q2)

def fiber_distance(a, b):
    """S³ = SU(2) distance: chordal distance on the 3-sphere."""
    return math.sqrt(sum((a[i] - b[i])**2 for i in range(4)))

def s4_key(p, precision=4):
    """Canonical S⁴ key for grouping points in the same fiber."""
    if p is None:
        return None
    if p[0] < 0:
        return tuple(round(-x, precision) for x in p)
    return tuple(round(x, precision) for x in p)

def main():
    print("=" * 60)
    print("HopfDNA Fiber Analysis — 28 exotic classes via S³ fiber")
    print("=" * 60)
    
    bases = ['A', 'C', 'G', 'T']
    all_seqs = [''.join(p) for p in itertools.product(bases, repeat=8)]
    print(f"\nTotal sequences: {len(all_seqs)}")
    
    # Group sequences by S⁴ base point (fiber)
    base_to_fibers = defaultdict(list)
    for seq in all_seqs:
        s7 = dna_to_s7(seq)
        s4 = hopf_map(s7)
        fiber = fiber_coordinate(s7)
        key = s4_key(s4)
        base_to_fibers[key].append(fiber)
    
    print(f"Distinct S⁴ base points: {len(base_to_fibers)}")
    
    # For each base point, count distinct fiber elements
    total_fiber_elements = set()
    base_fiber_counts = []
    
    for base_key, fibers in base_to_fibers.items():
        # Deduplicate fiber elements within this base point
        unique_fibers = set()
        for f in fibers:
            # Round fiber coordinates
            fr = tuple(round(x, 3) for x in f)
            unique_fibers.add(fr)
        base_fiber_counts.append(len(unique_fibers))
        for uf in unique_fibers:
            total_fiber_elements.add((base_key, uf))
    
    print(f"\nFiber analysis:")
    print(f"  Total (base, fiber) pairs: {len(total_fiber_elements)}")
    print(f"  Avg fibers per base: {sum(base_fiber_counts)/len(base_fiber_counts):.1f}")
    print(f"  Max fibers per base: {max(base_fiber_counts)}")
    print(f"  Min fibers per base: {min(base_fiber_counts)}")
    
    # Cluster ALL fiber elements globally
    # The 28 exotic classes should appear as clusters in the fiber S³
    all_fibers = list(total_fiber_elements)
    fiber_points = [f for (_, f) in all_fibers]
    
    print(f"\n  Total distinct fiber elements: {len(fiber_points)}")
    
    # Cluster fiber elements via single-linkage on S³
    print(f"\n  S³ fiber clustering (single-linkage):")
    for eps in [0.5, 0.8, 1.0, 1.2, 1.3, 1.4, 1.41, 1.42, 1.5, 1.6, 1.8, 2.0]:
        n = len(fiber_points)
        parent = list(range(n))
        def find(x):
            while parent[x] != x:
                parent[x] = parent[parent[x]]
                x = parent[x]
            return x
        def union(x, y):
            parent[find(x)] = find(y)
        for i in range(min(n, 5000)):  # subsample for speed
            for j in range(i+1, min(n, 5000)):
                d = fiber_distance(fiber_points[i], fiber_points[j])
                if d < eps:
                    union(i, j)
        comps = len(set(find(i) for i in range(min(n, 5000))))
        marker = " ◄═══ 28" if comps == 28 else ""
        print(f"    eps={eps:.2f}: {comps} clusters{marker}")
    
    # Also try: count distinct fiber elements per base, 
    # then look at the DISTRIBUTION of these counts
    count_dist = Counter(base_fiber_counts)
    print(f"\n  Distribution of fiber counts per base:")
    for count, freq in sorted(count_dist.items()):
        print(f"    {count} fibers: {freq} base points")
    
    print(f"\n{'═' * 60}")

if __name__ == "__main__":
    main()