# Listener synchrony: common drive or mutual coupling? Kuramoto phase oscillators (Kuramoto 1984)
# with a common external drive. Each of 30 listeners has a phase theta_i (think of a slow rhythm
# such as heart rate or an attentional cycle) with its own natural frequency:
#   d theta_i/dt = w_i + F sin(phi_i(t) - theta_i) + (K/N) sum_j sin(theta_j - theta_i)
# phi_i(t) is the phase of the stimulus as listener i receives it (the story's slow envelope).
# F: how strongly the stimulus entrains a listener; K: how strongly listeners entrain one another.
# Synchrony = order parameter R = |mean_j exp(i theta_j)| (0 = none, 1 = perfect), averaged over the
# last half of the run. Two controls separate the mechanisms:
#   'tested alone'      - each listener hears the same stimulus in a separate session (K has no effect)
#   'random-delay story' - listeners sit together but each hears the story with its own random delay
import math, random

def run(F, K, alone=False, shuffled=False, N=30, T=120.0, dt=0.01, seed=5):
    rng = random.Random(seed)
    w = [2 * math.pi * rng.gauss(1.0, 0.08) for _ in range(N)]
    th = [rng.uniform(0, 2 * math.pi) for _ in range(N)]
    lag = [rng.uniform(0, 2 * math.pi) if shuffled else 0.0 for _ in range(N)]
    Keff = 0.0 if alone else K
    acc, cnt, phi = 0.0, 0, 0.0
    for step in range(int(T / dt)):
        t = step * dt
        phi = 2 * math.pi * 1.0 * t + 0.6 * math.sin(0.3 * t)          # stimulus phase with slow wander
        cx = sum(math.cos(x) for x in th) / N; sx = sum(math.sin(x) for x in th) / N
        new = []
        for i in range(N):
            d = w[i] + F * math.sin(phi + lag[i] - th[i]) + Keff * (sx * math.cos(th[i]) - cx * math.sin(th[i]))
            new.append(th[i] + d * dt)
        th = new
        if t > T / 2:
            acc += math.hypot(sum(math.cos(x) for x in th) / N, sum(math.sin(x) for x in th) / N); cnt += 1
    return acc / cnt

print("mechanism                       | together, same story | tested alone | together, random delays")
for label, F, K in (("none", 0.0, 0.0), ("common drive only (F = 1.5)", 1.5, 0.0),
                    ("mutual coupling only (K = 2)", 0.0, 2.0), ("both", 1.5, 2.0)):
    print(f"{label:31s} | {run(F, K):20.2f} | {run(F, K, alone=True):12.2f} | {run(F, K, shuffled=True):23.2f}")
# mechanism                       | together, same story | tested alone | together, random delays
# none                            |                 0.15 |         0.15 |                    0.15
# common drive only (F = 1.5)     |                 0.93 |         0.93 |                    0.26
# mutual coupling only (K = 2)    |                 0.96 |         0.15 |                    0.96
# both                            |                 0.99 |         0.93 |                    0.75

# Inter-subject correlation (ISC) of heart rate, with a circular-shift surrogate test, on SYNTHETIC
# data. 20 listeners hear the same 5-minute story; heart rate is sampled once per second. Each
# listener's heart rate = a shared, story-driven slow component (strength a) + their own slow
# fluctuations (AR(1) noise). ISC = mean pairwise Pearson correlation. The null distribution comes
# from circularly shifting each listener's series by a random offset, which keeps each series'
# own structure but breaks its alignment with the story.
import math, random

def series(a, n_sub=20, T=300, seed=11):
    rng = random.Random(seed)
    story, s = [], 0.0
    for t in range(T):                                   # slow stimulus-driven signal
        s = 0.95 * s + rng.gauss(0, 1); story.append(s)
    subs = []
    for _ in range(n_sub):
        e, x = 0.0, []
        for t in range(T):
            e = 0.9 * e + rng.gauss(0, 1)
            x.append(70 + a * story[t] + 1.5 * e)
        subs.append(x)
    return subs

def pearson(x, y):
    n = len(x); mx = sum(x) / n; my = sum(y) / n
    sxy = sum((a - mx) * (b - my) for a, b in zip(x, y))
    return sxy / math.sqrt(sum((a - mx) ** 2 for a in x) * sum((b - my) ** 2 for b in y))

def isc(subs):
    n = len(subs)
    return sum(pearson(subs[i], subs[j]) for i in range(n) for j in range(i + 1, n)) / (n * (n - 1) / 2)

def surrogate_p(subs, n_surr=200, seed=3):
    rng = random.Random(seed); obs = isc(subs); T = len(subs[0]); hits = 0
    for _ in range(n_surr):
        shifted = []
        for x in subs:
            k = rng.randrange(30, T - 30); shifted.append(x[k:] + x[:k])
        hits += isc(shifted) >= obs
    return obs, (hits + 1) / (n_surr + 1)

for a, label in ((1.0, "attentive listening (shared component strong)"), (0.3, "distracted listening (shared component weak)"),
                 (0.0, "no shared component")):
    obs, p = surrogate_p(series(a))
    print(f"{label:46s}: ISC = {obs:+.3f}, surrogate p = {p:.3f}")
# attentive listening (shared component strong) : ISC = +0.589, surrogate p = 0.005
# distracted listening (shared component weak)  : ISC = +0.117, surrogate p = 0.005
# no shared component                           : ISC = -0.001, surrogate p = 0.443