SilverSight/python/hopf_helical.py

#!/usr/bin/env python3
"""
HopfDNA Helical Manifold — 28 exotic classes as winding numbers

Instead of clustering S⁴ points (fighting discretization), model the
braid manifold as a HELICAL manifold where:

  - Each DNA sequence maps to an angle θ on the helix
  - The pitch is the golden corkscrew angle ψ = 2π/φ²
  - The winding number n = θ / (2π) mod 1
  - The 28 exotic classes emerge as 28 distinct winding number slots

The helix never exactly repeats (ψ/2π = 1/φ² is irrational), but its
continued-fraction convergents give rational approximations. The 28
emerges from the convergent structure of the golden ratio.

Key insight: DNA IS a helix. The double helix has 10.5 bases per turn.
The exotic sphere count 28 may relate to how many full helical turns
fit in the 8-base window before the golden-angle irrationality forces
a new symmetry class.
"""

import math
import itertools
from collections import Counter

PHI = (1 + math.sqrt(5)) / 2
PSI = 2 * math.pi / (PHI ** 2)  # corkscrew angle ≈ 2.417 rad

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_norm_sq(q): return q[0]**2+q[1]**2+q[2]**2+q[3]**2

def quat_angle(q):
    """The rotation angle of a unit quaternion q = (cos θ/2, sin θ/2 * axis)."""
    # θ = 2 * arccos(real part)
    real = max(-1.0, min(1.0, q[0]))
    return 2 * math.acos(real)

def dna_to_helix_angle(seq):
    """
    Map 8-base DNA sequence to a helical angle.
    
    Each base contributes a rotation in the S³ fiber:
      A → 0 (identity, no rotation)
      C → π/2 (90° x-rotation)  
      G → π   (180° x-rotation)
      T → 3π/2 (270° x-rotation)
    
    The helical angle is the total accumulated rotation,
    scaled by the golden corkscrew pitch ψ = 2π/φ².
    """
    base_angles = {'A': 0.0, 'C': math.pi/2, 'G': math.pi, 'T': 3*math.pi/2}
    
    # Accumulate rotation through all 8 bases
    total_angle = sum(base_angles[b] for b in seq)
    
    # Scale by golden corkscrew pitch
    helical_angle = total_angle * PSI / (2 * math.pi)
    
    return helical_angle

def winding_number(helical_angle, period=2*math.pi):
    """The winding number n = floor(angle / period)."""
    return int(math.floor(helical_angle / period))

def winding_slot(helical_angle, n_slots=28):
    """
    Map helical angle to a slot in [0, n_slots).
    The slot is the winding number modulo n_slots.
    """
    normalized = helical_angle / (2 * math.pi)
    return int(normalized * n_slots) % n_slots

def continued_fraction_convergents(x, n_terms=20):
    """Compute continued fraction convergents of x."""
    convergents = []
    a_prev, a_curr = 0, 1
    b_prev, b_curr = 1, 0
    
    val = x
    for _ in range(n_terms):
        if val == 0: break
        int_part = int(math.floor(val))
        frac_part = val - int_part
        a_prev, a_curr = a_curr, int_part * a_curr + a_prev
        b_prev, b_curr = b_curr, int_part * b_curr + b_prev
        convergents.append((a_curr, b_curr))
        if frac_part < 1e-12: break
        val = 1.0 / frac_part
    return convergents

def main():
    print("=" * 60)
    print("HopfDNA Helical Manifold Analysis")
    print("=" * 60)
    
    print(f"\nGolden ratio φ = {PHI:.10f}")
    print(f"Corkscrew angle ψ = 2π/φ² = {PSI:.10f} rad")
    print(f"ψ / (2π) = 1/φ² = {1/PHI**2:.10f}")
    print(f"This is irrational → helix never exactly repeats")
    
    # Continued fraction of 1/φ²
    print(f"\nContinued fraction convergents of 1/φ²:")
    convs = continued_fraction_convergents(1/PHI**2)
    for i, (p, q) in enumerate(convs[:15]):
        ratio = p/q if q else 0
        error = abs(ratio - 1/PHI**2)
        print(f"  n={i}: {p}/{q} = {ratio:.10f}  (error: {error:.2e})")
    
    # Generate all 4⁸ sequences
    bases = ['A', 'C', 'G', 'T']
    all_seqs = [''.join(p) for p in itertools.product(bases, repeat=8)]
    
    # Map each to helical angle and winding slot
    print(f"\n--- 8-base DNA helical analysis ---")
    print(f"Total sequences: {len(all_seqs)}")
    
    angles = []
    slots_28 = []
    winding_numbers = []
    
    for seq in all_seqs:
        theta = dna_to_helix_angle(seq)
        angles.append(theta)
        slots_28.append(winding_slot(theta, 28))
        winding_numbers.append(winding_number(theta))
    
    slot_dist = Counter(slots_28)
    wn_dist = Counter(winding_numbers)
    
    print(f"\nWinding number distribution (mod 2π):")
    print(f"  Distinct winding numbers: {len(wn_dist)}")
    for wn, count in sorted(wn_dist.items()):
        print(f"  n={wn}: {count} sequences")
    
    print(f"\nWinding slots (mod 28):")
    print(f"  Non-empty slots: {len(slot_dist)}")
    print(f"  Empty slots: {28 - len(slot_dist)}")
    
    if len(slot_dist) == 28:
        print(f"\n  ✓ ALL 28 SLOTS POPULATED — exact 28 with no intervention!")
    else:
        print(f"\n  Slots populated: {len(slot_dist)}/28")
        # Try different slot counts
        print(f"\n  Testing slot counts:")
        for n_slots in [4, 7, 8, 14, 16, 20, 24, 28, 32, 35, 40, 56]:
            slots = [winding_slot(dna_to_helix_angle(s), n_slots) for s in all_seqs]
            dist = Counter(slots)
            print(f"    n={n_slots:3d}: {len(dist):3d}/{n_slots} slots filled", end="")
            if len(dist) == n_slots: print(" ✓ FULL", end="")
            if n_slots == 28: print(" ← 28", end="")
            print()
    
    # ── The key: use the golden angle directly ──────────────────────
    print(f"\n--- Golden angle slot analysis ---")
    golden_angle = 2 * math.pi * (1 - 1/PHI)  # ≈ 137.5°
    print(f"Golden angle = 2π(1-1/φ) = {golden_angle:.6f} rad = {math.degrees(golden_angle):.1f}°")
    
    # Map each sequence to a point on the golden-angle circle
    golden_slots = []
    for seq in all_seqs:
        theta = dna_to_helix_angle(seq)
        # The golden angle distributes points uniformly on S¹
        slot = int((theta / golden_angle) % 1.0 * 28)
        golden_slots.append(slot)
    
    g_dist = Counter(golden_slots)
    print(f"Golden-angle slots (mod 28): {len(g_dist)}/28 filled")
    
    # ── The REAL approach: use S³ fiber winding ─────────────────────
    print(f"\n--- S³ fiber winding (the correct approach) ---")
    
    def dna_to_s3_fiber(seq):
        """Map DNA → S³ fiber element (unit quaternion via composition)."""
        q = (1.0, 0.0, 0.0, 0.0)  # identity
        for b in seq:
            # Each base = 90° rotation about its axis
            angle = math.pi / 2
            half = angle / 2
            cos_h = math.cos(half)
            sin_h = math.sin(half)
            base_q = DNA_TO_Q[b]
            rot = (cos_h, sin_h * base_q[1], sin_h * base_q[2], sin_h * base_q[3])
            q = quat_mul(q, rot)
        return q
    
    # Compute S³ fiber for all sequences
    fiber_angles = []
    for seq in all_seqs:
        q = dna_to_s3_fiber(seq)
        # The rotation angle of q = 2*arccos(q[0])
        theta = quat_angle(q)
        fiber_angles.append(theta)
    
    # The winding number in S³ is θ/(2π)
    # Count distinct winding numbers
    s3_windings = [int(round(theta / (2 * math.pi))) for theta in fiber_angles]
    s3_dist = Counter(s3_windings)
    
    print(f"S³ fiber winding numbers (θ/2π):")
    print(f"  Distinct windings: {len(s3_dist)}")
    for wn, count in sorted(s3_dist.items()):
        bar = "█" * min(40, count // 50)
        print(f"  k={wn:3d}: {count:5d} {bar}")
    
    print(f"\n{'═' * 60}")
    if len(s3_dist) == 28:
        print(f"✓ EXACT 28: S³ fiber winding gives {len(s3_dist)} classes")
    else:
        # Use mod 28
        mod28 = [wn % 28 for wn in s3_windings]
        mod28_dist = Counter(mod28)
        print(f"  S³ windings mod 28: {len(mod28_dist)} distinct classes")
        if len(mod28_dist) == 28:
            print(f"  ✓ EXACT 28 via modular winding!")
        else:
            # The CORRECT count comes from the S³ fiber
            # Count up to winding 13 (since σ²⁸ = id, period is 28)
            # But the DNA grid can only produce windings 0..7
            # (8 bases × π/2 each = max 4π = 2 full turns)
            print(f"  Max winding from 8 bases: {max(s3_dist)} (= 4π/2π = 2 turns)")
            print(f"  Need longer sequences for 28 windings")
            print(f"  56-base DNA → 56×π/2 = 28π → 14 turns → still need 2× → 14 windings")
            print(f"  112-base DNA → 28 full S³ turns → 28 exotic classes")
            print(f"\n  OR: use golden-angle pitch (non-π/2) spacing")
    print(f"{'═' * 60}")

if __name__ == "__main__":
    main()