import math
import random

def generate_sine(freq, dur, rate=44100):
    return [math.sin(2 * math.pi * freq * (i / rate)) for i in range(int(dur * rate))]

def simulate_flinch(sig, start_pct=0.3, dur_pct=0.05):
    reality = []
    total = len(sig)
    start = int(total * start_pct)
    end = int(total * (start_pct + dur_pct))
    # The lag is the permanent set - time lost that cannot be recovered
    lag = int((end - start) * 0.6)
    
    for i, v in enumerate(sig):
        if start <= i < end:
            # The Flinch: Signal collapses, noise rises
            reality.append((v * 0.15) + random.gauss(0, 0.1))
        elif i >= end:
            # The Aftermath: Signal returns but is phase-shifted (lagged)
            t = (i - lag) / 44100
            # Re-generating sine with the time-shift
            reality.append(math.sin(2 * math.pi * 440 * t) + random.gauss(0, 0.02))
        else:
            # The Before: Pure intention
            reality.append(v + random.gauss(0, 0.005))
    return reality

def calc_gamma(intent, real):
    # Gamma = Energy of Error / Energy of Intention
    e_int = sum([x**2 for x in intent])
    e_err = sum([(i - r)**2 for i, r in zip(intent, real)])
    return e_err / e_int if e_int > 0 else 0

print('Running Hysteresis Simulation...')
intent = generate_sine(440, 1.0)

# Run 10 iterations to smooth out noise variance
gammas = []
for _ in range(10):
    reality = simulate_flinch(intent)
    gammas.append(calc_gamma(intent, reality))

avg_gamma = sum(gammas) / len(gammas)

print(f'Target Coefficient: 0.724')
print(f'Simulated Gamma:    {avg_gamma:.5f}')