The probabilistic method

"""Random colourings and Ramsey numbers: what the first moment gives.

Colour each edge of K_n red or blue by a fair coin. The expected number of
monochromatic copies of K_k is E(n, k) = C(n, k) 2^(1 - C(k, 2)); when it is
below 1, some colouring has none, so R(k, k) > n (Erdos, 1947). The deletion
method improves this: removing one vertex from each monochromatic K_k leaves
a good colouring on at least n - E(n, k) vertices. We tabulate both bounds,
simulate random colourings of K_6 to see that the probability of having no
monochromatic K_4 really is at least 1 - E(6, 4), and check by exhaustive
search that the Paley graph on 17 vertices (edges between quadratic residue
differences mod 17) has no K_4 and no independent 4-set, an explicit
colouring proving R(4, 4) > 17 - far better than the random bound for such
small k. Exact arithmetic for the bounds, exhaustive search for Paley, and a
seeded simulation (statistical evidence) for the random colourings.
"""

import itertools
import math
import random
from fractions import Fraction

def expected_monochromatic(n, k):
    """E(n, k) = C(n, k) 2^(1 - C(k, 2)) as an exact fraction."""
    return Fraction(math.comb(n, k) * 2, 2 ** math.comb(k, 2))

def first_moment_bound(k):
    """Largest n with E(n, k) < 1, so that R(k, k) > n."""
    n = k
    while expected_monochromatic(n + 1, k) < 1:
        n += 1
    return n

def deletion_bound(k):
    """max over n of n - floor(E(n, k)): R(k, k) exceeds this number.

    Some colouring has at most floor(E) monochromatic K_k; deleting one
    vertex from each leaves at least n - floor(E) vertices and none.
    """
    best, n = 0, k
    while True:
        value = n - math.floor(expected_monochromatic(n, k))
        if value < best - 5:
            return best
        best = max(best, value)
        n += 1

def random_colouring(n, rng):
    """A dict {(i, j): bool} for i < j, each edge red with probability 1/2."""
    return {e: rng.random() < 0.5 for e in itertools.combinations(range(n), 2)}

def count_monochromatic(colouring, n, k):
    """Number of k-subsets all of whose edges have the same colour."""
    count = 0
    for subset in itertools.combinations(range(n), k):
        colours = {colouring[e] for e in itertools.combinations(subset, 2)}
        count += len(colours) == 1
    return count

def paley_colouring(q):
    """Paley colouring of K_q, q prime = 1 mod 4: red iff i - j is a square."""
    squares = {(x * x) % q for x in range(1, q)}
    return {(i, j): (j - i) % q in squares
            for i, j in itertools.combinations(range(q), 2)}

def main():
    print(" k  2^(k/2)  first moment: R(k,k) >  deletion: R(k,k) >")
    for k in range(3, 11):
        print(f"{k:>2}  {2 ** (k / 2):>7.1f}  {first_moment_bound(k):>22}"
              f"  {deletion_bound(k):>18}")
    rng = random.Random(1947)
    n, k, trials = 6, 4, 4000
    counts = [count_monochromatic(random_colouring(n, rng), n, k)
              for _ in range(trials)]
    e = expected_monochromatic(n, k)
    good = counts.count(0) / trials
    margin = 1.96 * math.sqrt(good * (1 - good) / trials)
    print(f"Random colourings of K_{n}, {trials} samples (seed 1947):")
    print(f"  mean number of monochromatic K_{k}: {sum(counts) / trials:.3f}"
          f"  (exact E = {float(e):.3f})")
    print(f"  fraction with none: {good:.3f} +- {margin:.3f}"
          f"  (guaranteed >= 1 - E = {1 - float(e):.3f})")
    paley = paley_colouring(17)
    mono = count_monochromatic(paley, 17, 4)
    print(f"Paley colouring of K_17: {mono} monochromatic K_4 among "
          f"{math.comb(17, 4)} sets,")
    print("  so R(4,4) > 17 (in fact R(4,4) = 18); the random bound gives > 6.")
    print()
    print("Epistemic status: the bounds are exact arithmetic and the Paley")
    print("check is exhaustive, so both are proofs of what they state; the")
    print("simulation is statistical evidence (95% interval shown) about one")
    print("seeded run.")

if __name__ == "__main__":
    main()

"""Concentration: Chebyshev versus Chernoff, and balancing sets.

For a sum S of n independent fair signs +-1, Chebyshev's inequality gives
P(S >= a) <= n / a^2, while the Chernoff bound gives P(S >= a) <=
exp(-a^2 / (2n)) - exponentially small. We compare both with the exact
binomial tail and with a seeded simulation for n = 100. Then we use the
Chernoff bound as the lesson does: colour n points red and blue at random;
for each of n given sets, the imbalance |chi(A)| exceeds sqrt(2 n ln 2n) with
probability below 1/n, so with positive probability no set does. For 64
random sets on 64 points we sample random colourings and record how often all
sets are within the bound and the best maximum imbalance found. Exact tails,
seeded simulation (statistical evidence).
"""

import math
import random

def exact_tail(n, a):
    """P(S >= a) for S a sum of n fair +-1 signs (S = 2H - n, H binomial)."""
    need = math.ceil((n + a) / 2)
    return sum(math.comb(n, h) for h in range(max(need, 0), n + 1)) / 2 ** n

def chebyshev_bound(n, a):
    """Chebyshev: P(S >= a) <= P(|S| >= a) <= n / a^2 (a > 0)."""
    return min(1.0, n / a ** 2)

def chernoff_bound(n, a):
    """Chernoff: P(S >= a) <= exp(-a^2 / (2n))."""
    return math.exp(-a * a / (2 * n))

def simulate_tail(n, a, trials, rng):
    """Fraction of seeded trials with S >= a, using random bit strings."""
    hits = 0
    for _ in range(trials):
        heads = bin(rng.getrandbits(n)).count("1")
        hits += 2 * heads - n >= a
    return hits / trials

def random_sets(n_points, n_sets, rng):
    """Each set contains each point independently with probability 1/2."""
    return [[i for i in range(n_points) if rng.random() < 0.5]
            for _ in range(n_sets)]

def max_imbalance(colouring, sets):
    """max over sets A of |sum of colours (+-1) of the points in A|."""
    return max(abs(sum(colouring[i] for i in a)) for a in sets)

def main():
    n = 100
    rng = random.Random(2024)
    print(f"Sum of n = {n} random signs: P(S >= a)")
    print(f"{'a':>4} {'exact':>10} {'simulated':>10} {'Chebyshev':>10} "
          f"{'Chernoff':>10}")
    for a in (10, 20, 30, 40):
        sim = simulate_tail(n, a, 20000, rng)
        print(f"{a:>4} {exact_tail(n, a):>10.2e} {sim:>10.2e} "
              f"{chebyshev_bound(n, a):>10.2e} {chernoff_bound(n, a):>10.2e}")
    points = 64
    sets = random_sets(points, points, rng)
    bound = math.sqrt(2 * points * math.log(2 * points))
    trials = 500
    maxima = sorted(max_imbalance([rng.choice((-1, 1)) for _ in range(points)],
                                  sets) for _ in range(trials))
    within = sum(m <= bound for m in maxima)
    print(f"Balancing {points} random sets on {points} points:")
    print(f"  bound sqrt(2n ln 2n) = {bound:.1f}; random colourings within "
          f"it: {within}/{trials}")
    print(f"  maximum imbalance: median {maxima[trials // 2]}, "
          f"best {maxima[0]}, worst {maxima[-1]}")
    print()
    print("Epistemic status: the tails are exact; the simulated columns are")
    print("statistical evidence from seeded runs. The bounds themselves are")
    print("theorems (Chebyshev; Chernoff via cosh x <= exp(x^2/2)).")

if __name__ == "__main__":
    main()

"""The algorithmic local lemma: Moser and Tardos's resampling algorithm.

A k-CNF formula has clauses of k literals; under a uniformly random
assignment each clause fails with probability p = 2^-k. If every clause
shares variables with at most d others and e p (d + 1) <= 1, the Lovasz local
lemma says the formula is satisfiable - but its proof gives no algorithm,
and the probability that a random assignment works can be exponentially
small. Moser and Tardos (2010) showed that the obvious algorithm is fast:
start at random, and while some clause is violated, resample its variables.
The expected number of resamplings is at most m / d for m clauses. We build a
random 5-CNF instance with dependency degree at most 10 (so e p (d+1) < 1),
run the algorithm from 30 seeds, check the answers, and compare the number of
resamplings with the bound m / d and with naive rejection sampling. Seeded
simulation: statistical evidence about one instance, with every output
assignment verified exactly.
"""

import math
import random

def neighbours(clauses, index, occurrences):
    """Indices of the other clauses sharing a variable with clause index."""
    out = set()
    for var, _ in clauses[index]:
        out |= occurrences[var]
    out.discard(index)
    return out

def build_instance(n_vars, n_clauses, k, max_degree, rng, max_tries=20000):
    """Random k-CNF; a clause is kept only if all dependency degrees stay <= max_degree.

    A clause is a tuple of (variable, wanted_value) pairs; it is satisfied
    when at least one variable has its wanted value.
    """
    clauses, occurrences = [], [set() for _ in range(n_vars)]
    for _ in range(max_tries):
        if len(clauses) == n_clauses:
            break
        vars_ = rng.sample(range(n_vars), k)
        clause = tuple((v, rng.random() < 0.5) for v in vars_)
        clauses.append(clause)
        new = len(clauses) - 1
        for v in vars_:
            occurrences[v].add(new)
        touched = neighbours(clauses, new, occurrences) | {new}
        if any(len(neighbours(clauses, c, occurrences)) > max_degree
               for c in touched):
            for v in vars_:
                occurrences[v].discard(new)
            clauses.pop()
    return clauses

def max_dependency_degree(clauses, n_vars):
    """The largest number of other clauses any clause shares a variable with."""
    occurrences = [set() for _ in range(n_vars)]
    for i, clause in enumerate(clauses):
        for v, _ in clause:
            occurrences[v].add(i)
    return max(len(neighbours(clauses, i, occurrences))
               for i in range(len(clauses)))

def satisfied(clause, assignment):
    """True if some literal of the clause holds."""
    return any(assignment[v] == want for v, want in clause)

def moser_tardos(clauses, n_vars, rng):
    """Resample violated clauses until none is left; return (assignment, steps)."""
    assignment = [rng.random() < 0.5 for _ in range(n_vars)]
    steps = 0
    while True:
        bad = [c for c in clauses if not satisfied(c, assignment)]
        if not bad:
            return assignment, steps
        for v, _ in rng.choice(bad):
            assignment[v] = rng.random() < 0.5
        steps += 1

def naive_successes(clauses, n_vars, tries, rng):
    """How many of `tries` uniformly random assignments satisfy every clause."""
    hits = 0
    for _ in range(tries):
        bits = rng.getrandbits(n_vars)
        hits += all(any(((bits >> v) & 1) == want for v, want in clause)
                    for clause in clauses)
    return hits

def main():
    k, d, n_vars = 5, 10, 400
    rng = random.Random(2010)
    clauses = build_instance(n_vars, 150, k, d, rng)
    m = len(clauses)
    deg = max_dependency_degree(clauses, n_vars)
    p = 2.0 ** -k
    print(f"Random {k}-CNF: {n_vars} variables, {m} clauses, "
          f"max dependency degree {deg}")
    print(f"  LLL condition e p (d+1) = {math.e * p * (deg + 1):.3f} <= 1")
    runs = [moser_tardos(clauses, n_vars, random.Random(seed))
            for seed in range(30)]
    ok = all(all(satisfied(c, a) for c in clauses) for a, _ in runs)
    steps = sorted(s for _, s in runs)
    print(f"Moser-Tardos from 30 seeds: all outputs satisfy every clause: {ok}")
    print(f"  resamplings: mean {sum(steps) / len(steps):.1f}, "
          f"max {steps[-1]}; bound on the mean m/d = {m / deg:.1f}")
    tries = 2000
    hits = naive_successes(clauses, n_vars, tries, rng)
    print(f"Naive: {hits} of {tries} uniformly random assignments satisfy "
          "the formula;")
    print(f"  the LLL only guarantees probability >= (1 - 1/(d+1))^m = "
          f"{(1 - 1 / (deg + 1)) ** m:.1e}")
    print()
    print("Epistemic status: every output assignment is verified exactly;")
    print("the step counts and the naive success rate are statistical")
    print("evidence about one seeded instance. The bound m/d on the expected")
    print("number of resamplings is Moser and Tardos's theorem.")

if __name__ == "__main__":
    main()