Greedy influence maximization in Python
Greedy seed selection with a (1-1/e) approximation guarantee.
- File
greedy_influence_maximization.py, 138 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/greedy-influence-maximization/. The authors are listed in prov/DDTH3Z.jsonl.
import random
from collections.abc import Callable, Hashable, Iterable
def greedy_influence_maximization(
vertices: Iterable[int], k: int, sigma: Callable[[frozenset[int]], float]
) -> list[int]:
"""Return k seeds chosen greedily for the spread function sigma, in the order picked.
sigma maps a frozenset of vertices to its (expected or estimated) spread.
Each round adds the vertex with the largest marginal gain
sigma(S + v) - sigma(S); ties go to the smallest vertex. If k exceeds the
number of vertices, every vertex is returned. The function makes at most
k * n + 1 calls to sigma, and is a (1 - 1/e)-approximation whenever sigma
is monotone and submodular with sigma(empty set) = 0 and the calls are exact.
"""
pool = sorted(set(vertices))
if k < 0:
raise ValueError("k must be non-negative")
seeds: list[int] = [] # @line L1
chosen: set[int] = set()
base = sigma(frozenset())
for _ in range(min(k, len(pool))): # @line L2
best, best_gain, best_value = -1, 0.0, base
for v in pool: # @line L3
if v in chosen:
continue
value = sigma(frozenset(chosen | {v}))
gain = value - base # @line L4
if best < 0 or gain > best_gain: # @line L5
best, best_gain, best_value = v, gain, value
seeds.append(best) # @line L6
chosen.add(best)
base = best_value
return seeds # @line L7
def _out_lists(n: int, arcs: Iterable[tuple[int, int, float]]) -> list[list[tuple[int, float]]]:
out: list[list[tuple[int, float]]] = [[] for _ in range(n)]
for u, v, p in arcs:
if not (0 <= u < n and 0 <= v < n):
raise ValueError("arc end out of range")
if not 0 <= p <= 1:
raise ValueError("probability or weight outside [0, 1]")
out[u].append((v, p))
return out
def ic_spread(
n: int, arcs: list[tuple[int, int, float]], seeds: Iterable[int], runs: int, rng: random.Random
) -> float:
"""Monte Carlo estimate of the expected spread under the independent cascade model.
arcs holds (u, v, p): when u becomes active it activates v with probability p,
once. Each run is a breadth-first cascade; the result is the mean number of
active vertices over the runs, seeds included. Time O(runs * (n + m)).
"""
out = _out_lists(n, arcs)
start = sorted(set(seeds))
if runs <= 0 or not start:
return 0.0
total = 0
for _ in range(runs):
active = [False] * n
for s in start:
active[s] = True
frontier = list(start)
count = len(start)
while frontier:
nxt = []
for u in frontier:
for v, p in out[u]:
if not active[v] and rng.random() < p:
active[v] = True
nxt.append(v)
count += 1
frontier = nxt
total += count
return total / runs
def lt_spread(
n: int, arcs: list[tuple[int, int, float]], seeds: Iterable[int], runs: int, rng: random.Random
) -> float:
"""Monte Carlo estimate of the expected spread under the linear threshold model.
arcs holds (u, v, b): u pushes v with weight b, and the weights into each
vertex sum to at most 1. Each run draws a uniform threshold in [0, 1] for
every vertex; a vertex turns active once the weight from its active
in-neighbours reaches its threshold. Time O(runs * (n + m)).
"""
out = _out_lists(n, arcs)
into = [0.0] * n
for _, v, b in arcs:
into[v] += b
if any(x > 1 + 1e-9 for x in into):
raise ValueError("the weights into a vertex sum to more than 1")
start = sorted(set(seeds))
if runs <= 0 or not start:
return 0.0
total = 0
for _ in range(runs):
theta = [rng.random() for _ in range(n)]
pressure = [0.0] * n
active = [False] * n
for s in start:
active[s] = True
frontier = list(start)
count = len(start)
while frontier:
nxt = []
for u in frontier:
for v, b in out[u]:
if active[v]:
continue
pressure[v] += b
if pressure[v] >= theta[v]:
active[v] = True
nxt.append(v)
count += 1
frontier = nxt
total += count
return total / runs
def memoized(sigma: Callable[[frozenset[int]], float]) -> Callable[[frozenset[int]], float]:
"""Wrap a spread function so that each seed set is evaluated once."""
cache: dict[Hashable, float] = {}
def f(s: frozenset[int]) -> float:
if s not in cache:
cache[s] = sigma(s)
return cache[s]
return f
test_greedy_influence_maximization.py 183 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/greedy-influence-maximization/. The authors are listed in prov/DDTH3Z.jsonl.
import math
import random
from fractions import Fraction
from itertools import combinations, product
from greedy_influence_maximization import greedy_influence_maximization, ic_spread, lt_spread, memoized
HALF = Fraction(1, 2)
# The example graph of the page: arcs 0->1, 0->2, 1->3, 2->3, 3->4, 5->4, all with 1/2.
EXAMPLE = [(0, 1, HALF), (0, 2, HALF), (1, 3, HALF), (2, 3, HALF), (3, 4, HALF), (5, 4, HALF)]
def reach(n, live, seeds):
out = [[] for _ in range(n)]
for u, v in live:
out[u].append(v)
seen = set(seeds)
stack = list(seeds)
while stack:
u = stack.pop()
for v in out[u]:
if v not in seen:
seen.add(v)
stack.append(v)
return len(seen)
def ic_exact(n, arcs):
"""sigma for the independent cascade model, by summing over all 2^m live-edge graphs."""
worlds = []
for mask in product([False, True], repeat=len(arcs)):
pr = Fraction(1)
live = []
for keep, (u, v, p) in zip(mask, arcs):
pr *= p if keep else 1 - p
if keep:
live.append((u, v))
if pr:
worlds.append((pr, live))
return lambda s: sum((pr * reach(n, live, s) for pr, live in worlds), Fraction(0))
def lt_exact(n, arcs):
"""sigma for the linear threshold model: each vertex keeps at most one in-arc, arc (u, v, b) with probability b."""
options = []
for v in range(n):
ins = [(u, b) for u, w, b in arcs if w == v]
opts = [(b, (u, v)) for u, b in ins if b]
rest = 1 - sum((b for u, b in ins), Fraction(0))
if rest:
opts.append((rest, None))
options.append(opts)
worlds = []
for choice in product(*options):
pr = Fraction(1)
live = []
for b, arc in choice:
pr *= b
if arc is not None:
live.append(arc)
worlds.append((pr, live))
return lambda s: sum((pr * reach(n, live, s) for pr, live in worlds), Fraction(0))
def best_value(n, k, sigma):
return max(sigma(frozenset(c)) for c in combinations(range(n), min(k, n)))
def random_ic(rng, n):
arcs = []
for u in range(n):
for v in range(n):
if u != v and rng.random() < 0.35:
arcs.append((u, v, rng.choice([Fraction(1, 4), HALF, Fraction(3, 4), Fraction(1)])))
return arcs[:9]
def random_lt(rng, n):
arcs = []
for v in range(n):
ins = [u for u in range(n) if u != v and rng.random() < 0.4]
if ins:
share = Fraction(rng.choice([2, 3, 4]), 4 * len(ins))
arcs += [(u, v, share) for u in ins]
return arcs
def test_example_ic():
sigma = ic_exact(6, EXAMPLE)
assert sigma(frozenset({0})) == Fraction(85, 32)
assert [sigma(frozenset({v})) for v in range(1, 6)] == [
Fraction(7, 4), Fraction(7, 4), Fraction(3, 2), Fraction(1), Fraction(3, 2)
]
seeds = greedy_influence_maximization(range(6), 2, sigma)
assert seeds == [0, 5]
assert sigma(frozenset(seeds)) == Fraction(259, 64)
assert sigma(frozenset({0, 5})) - sigma(frozenset({0})) == Fraction(89, 64)
assert sigma(frozenset({0, 3})) - sigma(frozenset({0})) == Fraction(27, 32)
assert sigma(frozenset({0, 1})) - sigma(frozenset({0})) == Fraction(25, 32)
assert best_value(6, 2, sigma) == Fraction(259, 64)
def test_example_lt():
sigma = lt_exact(6, EXAMPLE)
assert sigma(frozenset({0})) == Fraction(11, 4)
seeds = greedy_influence_maximization(range(6), 2, sigma)
assert seeds == [0, 5]
assert sigma(frozenset(seeds)) == Fraction(17, 4)
def test_coverage_gap():
# p = 1 turns the cascade into coverage: vertex 0 is a decoy set {6, 7, 9, 10},
# vertices 1 and 2 the sets {6, 7, 8} and {9, 10, 11}; vertices 3 to 5 are unused.
arcs = [(0, x, 1) for x in (6, 7, 9, 10)] + [(1, x, 1) for x in (6, 7, 8)] + [(2, x, 1) for x in (9, 10, 11)]
sigma = ic_exact(12, arcs)
seeds = greedy_influence_maximization(range(12), 2, sigma)
assert seeds == [0, 1]
assert sigma(frozenset(seeds)) == 7
assert best_value(12, 2, sigma) == 8
def test_edge_cases():
sigma = ic_exact(3, [])
assert greedy_influence_maximization(range(3), 0, sigma) == []
assert greedy_influence_maximization(range(3), 2, sigma) == [0, 1]
assert greedy_influence_maximization(range(3), 9, sigma) == [0, 1, 2]
assert greedy_influence_maximization([], 2, sigma) == []
try:
greedy_influence_maximization(range(3), -1, sigma)
except ValueError:
pass
else:
raise AssertionError("negative k accepted")
def test_approximation_bound_and_submodularity():
rng = random.Random(2003)
bound = 1 - 1 / math.e
for trial in range(60):
n = rng.randint(1, 6)
if trial % 2:
arcs = random_ic(rng, n)
sigma = memoized(ic_exact(n, arcs))
else:
arcs = random_lt(rng, n)
sigma = memoized(lt_exact(n, arcs))
for k in range(1, n + 1):
got = sigma(frozenset(greedy_influence_maximization(range(n), k, sigma)))
assert got >= bound * best_value(n, k, sigma)
for a in range(1 << n):
s = frozenset(i for i in range(n) if a >> i & 1)
for v in range(n):
if v in s:
continue
gain_s = sigma(s | {v}) - sigma(s)
assert gain_s >= 0
for w in range(n):
if w not in s and w != v:
assert sigma(s | {v, w}) - sigma(s | {w}) <= gain_s
def test_monte_carlo_estimates():
floats = [(u, v, float(p)) for u, v, p in EXAMPLE]
est = ic_spread(6, floats, [0], 20000, random.Random(7))
assert abs(est - 85 / 32) < 0.05
est = lt_spread(6, floats, [0], 20000, random.Random(7))
assert abs(est - 11 / 4) < 0.05
assert ic_spread(6, floats, [], 10, random.Random(1)) == 0.0
sigma = memoized(lambda s: ic_spread(6, floats, s, 4000, random.Random(11)))
assert greedy_influence_maximization(range(6), 2, sigma) == [0, 5]
if __name__ == "__main__":
test_example_ic()
test_example_lt()
test_coverage_gap()
test_edge_cases()
test_approximation_bound_and_submodularity()
test_monte_carlo_estimates()
print("ok")
What the next level needs
- the tests did not run in this build