# How Epidemics Spread 7, Exercise 2
# Exercise 2: where the random and hub networks start to take off, compared with the threshold T = 1/excess degree.
import random

def random_network(n, mean_degree):
    """Random network: pick random pairs until there are n * mean_degree / 2 links."""
    nbrs = [set() for _ in range(n)]
    links = 0
    while links < n * mean_degree // 2:
        a, b = random.randrange(n), random.randrange(n)
        if a != b and b not in nbrs[a]:
            nbrs[a].add(b)
            nbrs[b].add(a)
            links += 1
    return nbrs

def clustered_network(n, mean_degree, rewire):
    """Ring of neighbours (everyone knows the people nearby); each link is rewired at random with chance `rewire`."""
    nbrs = [set() for _ in range(n)]
    for i in range(n):
        for step in range(1, mean_degree // 2 + 1):
            j = (i + step) % n
            if random.random() < rewire:
                j = random.randrange(n)
                while j == i or j in nbrs[i]:
                    j = random.randrange(n)
            nbrs[i].add(j)
            nbrs[j].add(i)
    return nbrs

def hub_network(n, m):
    """Growth with preferential attachment: each newcomer links to m people, favouring the already well connected."""
    nbrs = [set() for _ in range(n)]
    pool = []                      # a person appears here once per link they have
    targets = list(range(m))
    for new in range(m, n):
        for t in set(targets):
            nbrs[new].add(t)
            nbrs[t].add(new)
        pool.extend(targets)
        pool.extend([new] * m)
        targets = [random.choice(pool) for _ in range(m)]
    return nbrs

def describe(nbrs):
    n = len(nbrs)
    degrees = [len(x) for x in nbrs]
    mean = sum(degrees) / n
    excess = (sum(d * d for d in degrees) / n - mean) / mean   # average number of *other* contacts a contact has
    tri, trip = 0, 0
    for i in range(n):
        d = degrees[i]
        if d >= 2:
            trip += d * (d - 1) / 2
            tri += sum(len(nbrs[i] & nbrs[j]) for j in nbrs[i]) / 2
    return mean, max(degrees), excess, tri / trip


def outbreak(nbrs, T, blocked=frozenset()):
    """One outbreak from a random unblocked person; returns how many people were ever infected."""
    start = random.randrange(len(nbrs))
    while start in blocked:
        start = random.randrange(len(nbrs))
    reached, frontier = {start}, [start]
    while frontier:
        following = []
        for u in frontier:
            for v in nbrs[u]:
                if v not in reached and v not in blocked and random.random() < T:
                    reached.add(v)
                    following.append(v)
        frontier = following
    return len(reached)

def summarise(nbrs, T, runs, blocked=frozenset()):
    sizes = [outbreak(nbrs, T, blocked) for _ in range(runs)]
    big = [s for s in sizes if s > 20]
    fizzled = 1 - len(big) / runs
    share = sum(big) / len(big) / len(nbrs) if big else 0.0   # 0 if none took off
    return fizzled, share, sum(sizes) / runs / len(nbrs)


random.seed(22)
n, runs = 2000, 300
random_g, hub_g = random_network(n, 6), hub_network(n, 3)
for name, g in [("random", random_g), ("hubs", hub_g)]:
    mean, biggest, excess, clustering = describe(g)
    print(f"{name}: mean degree {mean:.2f}, excess degree {excess:.2f}, predicted threshold T = 1/excess = {1 / excess:.3f}")
print("   T   random: average size   hubs: average size")
for T in [0.03, 0.05, 0.08, 0.12, 0.16, 0.20, 0.25]:
    a = summarise(random_g, T, runs)[2]
    b = summarise(hub_g, T, runs)[2]
    print(f"{T:5.2f} {a:20.1%} {b:20.1%}")

# Output when run:
# random: mean degree 6.00, excess degree 5.99, predicted threshold T = 1/excess = 0.167
# hubs: mean degree 5.96, excess degree 15.02, predicted threshold T = 1/excess = 0.067
#    T   random: average size   hubs: average size
#  0.03                 0.1%                 0.1%
#  0.05                 0.1%                 0.1%
#  0.08                 0.1%                 0.3%
#  0.12                 0.2%                 1.9%
#  0.16                 0.7%                 9.7%
#  0.20                 8.4%                20.3%
#  0.25                33.7%                32.3%

# Reading the output: the random network has a predicted threshold of about 0.167 and the hub network about 0.067 (my reasoning: T = 1 / excess degree). In the runs the hub network starts to grow earlier (1.9% at T = 0.12, 9.7% at T = 0.16) while the random network is still near zero (0.2% and 0.7%) and then catches up (8.4% against 20.3% at T = 0.20; 33.7% against 32.3% at T = 0.25). The takeoff is gradual rather than sharp, because the networks have only 2,000 people, so the thresholds are blurred and only roughly matched. That the hub network still shows a gradual takeoff is not necessarily in conflict with reports of no threshold for idealised infinite scale-free networks: this network is small and has a minimum degree of 3 (an inference).