1import json, math
  2from pathlib import Path
  3import numpy as np
  4from scipy.special import lambertw
  5
  6# Resource-driven collective attention toy model.
  7# Agents never observe targets; a captured target leaves only a transient motion/action cue.
  8
  9def lambda_star(ell, tau, r, B=0.64):
 10    K = math.pi/(4*r) * (1.0 - B*r/tau)
 11    z = ell*math.exp(K)/tau
 12    return ell/float(lambertw(z).real)
 13
 14def formula_check():
 15    # Directly verify logarithmic tau dependence of the claimed collective MFPT.
 16    b, D, ell, lam = 10.0, 1.0, 1.0, 1.0
 17    taus = np.array([4., 8., 16., 32., 64., 128., 256.])
 18    Tc = b*b/(2*D)*np.log(taus/(lam*np.exp(-ell/lam)))
 19    # finite differences dTc/d(log tau) should be constant b^2/(2D)
 20    slopes = np.diff(Tc)/np.diff(np.log(taus))
 21    return {"taus": taus.tolist(), "Tc": Tc.tolist(),
 22            "log_tau_slopes": slopes.tolist(),
 23            "expected_slope": b*b/(2*D),
 24            "max_relative_slope_error": float(np.max(np.abs(slopes-b*b/(2*D))/(b*b/(2*D))))}
 25
 26def torus_delta(a, b, L):
 27    d = a-b
 28    return (d + L/2) % L - L/2
 29
 30def simulate(attention_radius, seed, N=8, L=20., steps=300, n_targets=5,
 31             tau=35, target_r=.42, speed=.22, cue_steps=10, social=True):
 32    rng = np.random.default_rng(seed)
 33    x = rng.uniform(0, L, (N, 2))
 34    ang = rng.uniform(-math.pi, math.pi, N)
 35    targets = rng.uniform(0, L, (n_targets, 2))
 36    cooldown = np.zeros(n_targets, dtype=int)
 37    # cue vector and expiry are internal action history, not target information
 38    cues = np.zeros((N, 2)); expiry = np.zeros(N, dtype=int)
 39    captures = 0; duplicate = 0; agg_sum = 0.; agg_count = 0
 40    for t in range(steps):
 41        # local conspecific attention: copy the most recent successful motion cue
 42        new_ang = ang.copy()
 43        if social and attention_radius > 0:
 44            for i in range(N):
 45                d = torus_delta(x, x[i], L)
 46                dist = np.sqrt((d*d).sum(1)); neigh = np.where((dist <= attention_radius) & (dist > 0) & (expiry > t))[0]
 47                if len(neigh):
 48                    # nearest active cue; no proximity or group reward is used
 49                    j = neigh[np.argmin(dist[neigh])]
 50                    v = cues[j]
 51                    if np.linalg.norm(v) > 1e-8:
 52                        desired = math.atan2(v[1], v[0])
 53                        # social influence is modest, preserving exploration
 54                        diff = (desired-new_ang[i]+math.pi)%(2*math.pi)-math.pi
 55                        new_ang[i] += .65*diff
 56        # persistent random exploration
 57        new_ang += rng.normal(0, .20, N)
 58        ang = (new_ang + math.pi)%(2*math.pi)-math.pi
 59        x = (x + speed*np.c_[np.cos(ang), np.sin(ang)]) % L
 60        # target replenishment
 61        cooldown = np.maximum(0, cooldown-1)
 62        for k in range(n_targets):
 63            if cooldown[k] == 0 and np.linalg.norm(torus_delta(x[0], targets[k], L)) < -1: # no-op, keeps target hidden
 64                pass
 65        # each target can be captured once, then replenishes after tau
 66        for i in range(N):
 67            if not np.any(cooldown == 0): break
 68            ds = torus_delta(targets, x[i], L)
 69            hit = np.where((cooldown == 0) & ((ds*ds).sum(1) <= target_r**2))[0]
 70            if len(hit):
 71                k = hit[0]; captures += 1
 72                # capture cue is the agent's recent action, without exposing target location/reward
 73                cues[i] = np.array([math.cos(ang[i]), math.sin(ang[i])])
 74                expiry[i] = t + cue_steps
 75                cooldown[k] = tau
 76                # target is hidden while depleted, then respawns at an unobserved location
 77                targets[k] = rng.uniform(0, L, 2)
 78        # aggregation order parameter from the stated neighbor statistic
 79        for i in range(N):
 80            d = torus_delta(x, x[i], L); dist=np.sqrt((d*d).sum(1))
 81            neigh=np.where((dist>0)&(dist<1.5))[0]
 82            agg_sum += (np.sum(dist[neigh]<1.5)/len(neigh)) if len(neigh) else 0.
 83            agg_count += 1
 84        # duplicate-search proxy: close pair during depleted-target intervals
 85        if np.sum(cooldown > 0) > 0:
 86            for i in range(N):
 87                d=torus_delta(x, x[i], L); duplicate += int(np.sum(((d*d).sum(1)>0)&((d*d).sum(1)<target_r**2))>0)
 88    return {"capture_rate": captures/(steps/1000), "captures": captures,
 89            "aggregation_A": agg_sum/agg_count, "duplicate_fraction": duplicate/(steps*N)}
 90
 91def run():
 92    ell, tau, r = 1.0, 35.0, .42
 93    ls = lambda_star(ell, tau, r)
 94    multipliers = [.25,.5,.75,1.,1.5,2.,4.]
 95    seeds = [3, 11]
 96    rows=[]
 97    for m in multipliers:
 98        lam=m*ls
 99        vals=[simulate(lam,s) for s in seeds]
100        rows.append({"multiplier":m,"lambda":lam,
101          "capture_rate":float(np.mean([v['capture_rate'] for v in vals])),
102          "aggregation_A":float(np.mean([v['aggregation_A'] for v in vals])),
103          "duplicate_fraction":float(np.mean([v['duplicate_fraction'] for v in vals]))})
104    base=[simulate(0,s,social=False) for s in seeds]
105    out={"formula_check":formula_check(),"parameters":{"ell":ell,"tau":tau,"r":r,"lambda_star":ls},
106         "sweep":rows,"independent_baseline":{
107          "capture_rate":float(np.mean([v['capture_rate'] for v in base])),
108          "aggregation_A":float(np.mean([v['aggregation_A'] for v in base])),
109          "duplicate_fraction":float(np.mean([v['duplicate_fraction'] for v in base]))}}
110    Path('results.json').write_text(json.dumps(out,indent=2))
111    print(json.dumps(out,indent=2))
112
113if __name__ == '__main__': run()