# Exact overlap and food-entry bounds: code and saved outputs

Reactions can depend on overlapping pools of possible catalysts. This program calculates how that overlap affects the probability that all reactions receive catalytic support. Finite probability checks use exact integers and fractions. Curves numerically evaluate
the paper's finite formulas; they are not sampled RAF frequencies.

## Catalyst overlap and joint support probability

![Two catalyst-set families and their exact probability difference](overlap_comparison.svg)

Each row requires at least one catalyst from its marked molecules. Both families
have two molecules per row and one shared molecule per pair of rows. When each molecule is independently available with probability p=1/2,
all three rows are satisfied in 9 of 16 assignments for Family I, and 8 of 16
for Family II. This is a joint event. The different event "at least one row is
satisfied" has probabilities 15/16 and 7/8. The small example is an abstract
requirements example from Theorem 6.1, not a complete binary-polymer RAF network.

## Fixed food constrains the growth of the model

![Gateway probability for fixed and linearly growing catalysis](gateway_scaling.svg)

The food set remains A, B, AA, AB, BA, BB. A channel is a reversible reaction;
one catalysis assignment applies to both directions. There are exactly 34
channels with a reaction side made entirely of food once n is at least 4.
Every RAF must include a catalysed such channel. The plots give the probability
that at least one of these channels has any catalyst in the whole molecule set.
This is an upper bound on RAF probability, since that catalyst might not be
reachable and the remaining required reactions might not have catalytic support.

Here f is the mean number of channels catalysed by one molecule. At fixed f, the bound tends to zero. At f=λn, it tends to 1−exp(−34λ), which
is below one for every finite λ. This paper does not determine the RAF limit.
The n axis measures maximum word length, not number of molecule types.

![Map of the gateway upper bound](gateway_bound_map.svg)

This is an upper-bound map, not a phase diagram of actual RAF existence. For
example, a 10% contour says that RAF probability is at most 10% there. A bound
near one does not establish that RAFs are common. Full numerical values are in gateway_map.csv.

## Actual RAF check in a small restricted network

Keep the six food molecules and all 30 ambient words of length at most four,
but allow only A + AA ⇌ AAA and A + AAA ⇌ AAAA. The two nonempty
food-generated supports are the first channel alone and both channels together.
The script compares the closure-based RAF definition, the exhaustive support
catalogue, and iterative pruning for all 65,536 assignments involving reachable
molecules. The 44 remaining catalysis coordinates cannot support this restricted
network; their possible assignments are nevertheless counted in fixed-Q results.

If only AAAA catalyses the first channel, an entry is catalysed but no RAF
exists. If AAA also catalyses the second channel, both channels form a RAF.
Closure is computed without demanding a catalyst at every construction step;
the RAF condition is not a claim of an entirely catalysed startup sequence.

For q=1−p the exact RAF probability of this restricted network is
1−q⁷+q⁷p(1−q⁸). At p=1/2 it is 65279/65536. These are not probabilities
for the full 64-channel network at n=4.

## Saved diagnostic output

```text
Family I: all requirements at p=1/2 = 9/16; at least one = 15/16.
Family II: all requirements at p=1/2 = 1/2; at least one = 7/8.
Exact polynomial difference: q^3 - q^4 = p(1-p)^3.
Its maximum is 27/256 at p=1/4 (differentiate p(1-p)^3).
Dedicated channels: Bernoulli union = 63/64; at fixed Q=2, union = 17/22, product of fixed-Q marginals would give 7273/10648 (not equal).
Both inclusion-exclusion levels agree with literal enumeration; fixed-Q coefficients agree for every Q in each checked catalogue.
Channel counts, n=2..10: 4, 18, 64, 188, 502, 1264, 3048, 7130, 16334
Food-accessible counts, n=2..10: 4, 18, 34, 34, 34, 34, 34, 34, 34
Restricted n=4 system: checked all 65,536 relevant assignments; RAF definition, catalogue activation, and pruning agree.
Its exact RAF probability at p=1/2 is 65279/65536; all 61 fixed-Q coefficients on the 60 ambient coordinates agree.
A catalysed entry without a RAF, and a two-channel RAF, both verified.
Gateway bound at f=1.3:
  n=   6: 0.999985014
  n=  10: 0.996060576
  n=  64: 0.509781087
  n= 128: 0.295870965
  n= 512: 0.083017294
  n=1024: 0.042326652
  n=4096: 0.010738217
Linear-scale example λ=0.03:
  n=   6: 0.784839513
  n=  64: 0.651076726
  n=4096: 0.639584696
  limiting gateway bound: 0.639405060
All exact finite checks passed. General asymptotic claims rely on the paper.
```

The fixed-Q tables count assignments directly; their probabilities are exact
fractions. Fixed-Q sampling does not make even disjoint channel events
independent. The independent-core product is checked under Bernoulli sampling.

The manuscript reports Lean verification except for Proposition 2.1 and
Corollary 9.6. This program does not compile Lean or audit those proof files.
It checks the stated finite examples and evaluates the gateway formulas.

## Run and inspect

Install numpy and matplotlib, then run:

```bash
python -X utf8 overlap_raf.py --output results
```

The program needs no external data, random seed, account, or API key.
CSV files contain all plotted data, the reaction census, the 34 entry channels,
the 16 catalyst subsets for each family, and every fixed-Q result for the
restricted network. SVG copies of the figures support publication at any size.
diagnostics.json records the source hash. Inclusion-exclusion is exponential;
the script deliberately applies it only to small catalogues.


[Download the complete package](overlap_package.zip). The source below includes the revised chart title; its scientific calculations are unchanged. UTF-8 mode is needed on Windows.

The exact channel count gives the sharper coefficient $1+o(1)$ after Proposition 9.1. The bound $272f_n/n$ and limiting constant 34 are unaffected. Theorem 7.2 supplies a sufficient coordinate-disjointness hypothesis for the Bernoulli product.

[Download the manuscript TeX source](manuscript_source.zip).

## Complete Python source

```python
#!/usr/bin/env python3
"""Exact finite examples and gateway bounds for the supplied RAF paper.

Run: python overlap_raf.py --output results
Dependencies: Python >= 3.10, numpy, matplotlib. No downloads or random samples.
Fractions and integer coefficient lists give exact finite checks. Plots use floats.
The finite examples do not replace the paper's proofs for arbitrary networks or n.
"""

import argparse
import csv
import hashlib
import json
import math
from collections import Counter
from fractions import Fraction as F
from functools import lru_cache
from itertools import combinations, product
from pathlib import Path


def subsets(items):
    items = tuple(items)
    for k in range(len(items) + 1):
        yield from combinations(items, k)


def multiply(a, b):
    """Multiply polynomials represented by coefficients in increasing degree."""
    c = [0] * (len(a) + len(b) - 1)
    for i, x in enumerate(a):
        for j, y in enumerate(b):
            c[i + j] += x * y
    return c


def union_sizes(requirements):
    """Signed counts of union sizes: the inner inclusion-exclusion calculation.

    Empty sets mean unused channels, as in the paper; repeated sets count once.
    """
    family = sorted({frozenset(w) for w in requirements if w}, key=lambda w: sorted(w))
    terms = Counter()
    for chosen in subsets(family):
        union = frozenset().union(*chosen)
        terms[len(union)] += (-1) ** len(chosen)
    return dict(sorted(terms.items()))


def local_probability(requirements, p):
    """Probability that a channel's catalyst set meets every requirement."""
    if not 0 <= p <= 1:
        raise ValueError('p must be between 0 and 1')
    return sum(c * (1 - p) ** d for d, c in union_sizes(requirements).items())


def local_counts(requirements, molecule_count):
    """Coefficient Q counts satisfying assignments with exactly Q catalysts."""
    answer = [0] * (molecule_count + 1)
    for size, sign in union_sizes(requirements).items():
        for q in range(molecule_count - size + 1):
            answer[q] += sign * math.comb(molecule_count - size, q)
    return answer


def catalogue_law(cores, molecule_count, channel_count, p):
    """Theorem 5.1: exact probability of AT LEAST ONE activated core.

    A core is {channel_index: frozenset(eligible_molecule_indices)}.
    Missing/empty requirements mean unused channels. All ambient channels,
    including unused ones, contribute to the fixed-Q count polynomial.
    Runtime is exponential in catalogue size; use only small catalogues here.
    """
    for core in cores:
        for j, w in core.items():
            if not 0 <= j < channel_count or any(not 0 <= x < molecule_count for x in w):
                raise ValueError('Requirement outside the declared molecule/channel sets')
    probability = F(0)
    counts = [0] * (molecule_count * channel_count + 1)
    for selected in subsets(range(len(cores))):
        if not selected:
            continue
        sign = (-1) ** (len(selected) + 1)
        joint = F(1)
        joint_counts = [1]
        for j in range(channel_count):
            requirements = [cores[k].get(j, frozenset()) for k in selected]
            joint *= local_probability(requirements, p)
            joint_counts = multiply(joint_counts, local_counts(requirements, molecule_count))
        probability += sign * joint
        counts = [a + sign * b for a, b in zip(counts, joint_counts)]
    return probability, counts


def activated(core, fibres):
    return all(not w or bool(w & fibres[j]) for j, w in core.items())


def enumerate_catalogue(cores, molecule_count, channel_count):
    """Independent, literal enumeration; no inclusion-exclusion is used."""
    total = molecule_count * channel_count
    if total > 20:
        raise ValueError('Literal enumeration is limited to 20 coordinates')
    counts = [0] * (total + 1)
    for mask in range(1 << total):
        fibres = [{x for x in range(molecule_count)
                   if mask & (1 << (j * molecule_count + x))}
                  for j in range(channel_count)]
        if any(activated(core, fibres) for core in cores):
            counts[mask.bit_count()] += 1
    return counts


def probability_from_counts(counts, p):
    total = len(counts) - 1
    return sum(c * p ** q * (1 - p) ** (total - q) for q, c in enumerate(counts))


def words(n):
    return [''.join(w) for length in range(1, n + 1)
            for w in product('AB', repeat=length)]


def repository_channels(n):
    """Section 2.1: keep both orders unless they display the same product."""
    molecules = words(n)
    order = {w: i for i, w in enumerate(molecules)}
    return [(u, v, u + v) for u in molecules for v in molecules
            if len(u) + len(v) <= n
            and (u + v != v + u or order[u] <= order[v])]


@lru_cache(None)
def primitive_counts(maximum):
    """P[d] counts binary words not equal to a shorter word repeated."""
    p = [0] + [1 << d for d in range(1, maximum + 1)]
    for d in range(1, maximum + 1):
        for multiple in range(2 * d, maximum + 1, d):
            p[multiple] -= p[d]
    return p


@lru_cache(None)
def channel_count(n):
    """Proposition 2.1, evaluated using integers; no network is generated."""
    if n < 1:
        raise ValueError('n must be positive')
    primitive = primitive_counts(n // 3)
    correction = sum(primitive[d] * (((n // d - 1) ** 2) // 4)
                     for d in range(1, n // 3 + 1))
    return (n - 2) * (1 << (n + 1)) + 4 - correction


def gateway_probability(n, f):
    """Numerically evaluate 1-(1-f/J[n])**(34*|X[n]|), n >= 4.

    Form K*p as an integer ratio first, so huge n does not underflow p before
    multiplying by K. For tiny p, log(1-p)/(-p)=1+p/2+p^2/3+... . The
    omitted terms below affect the exponent by less than 3e-25 relatively.
    """
    if not isinstance(n, int) or n < 4 or not math.isfinite(f) or f < 0:
        raise ValueError('Require integer n >= 4 and finite f >= 0')
    j = channel_count(n)
    f_exact = F(f)
    if f_exact > j:
        raise ValueError('f cannot exceed the total number of channels')
    if f == 0:
        return 0.0
    if f_exact == j:
        return 1.0
    k = 34 * ((1 << (n + 1)) - 2)
    p = float(f_exact / j)
    mean = float(f_exact * k / j)
    if p == 1.0:  # Rounding can reach 1 before the exact fraction does.
        return 1.0
    correction = (-math.log1p(-p) / p if p >= 1e-8 else 1 + p / 2 + p * p / 3)
    return -math.expm1(-mean * correction)


def closure(food, channels, support):
    """Build molecules using either direction, without requiring catalysts."""
    available = set(food)
    while True:
        enlarged = set(available)
        for j in support:
            u, v, w = channels[j]
            if u in available and v in available:
                enlarged.add(w)
            if w in available:
                enlarged.update((u, v))
        if enlarged == available:
            return available
        available = enlarged


def is_raf(food, channels, support, fibres):
    available = closure(food, channels, support)
    return bool(support) and all(set(channels[j]) <= available
                                and bool(fibres[j] & available) for j in support)


def max_raf(food, channels, fibres):
    """Repeated closure and deletion returns the largest RAF, or the empty set."""
    support = set(range(len(channels)))
    while support:
        available = closure(food, channels, support)
        smaller = {j for j in support if set(channels[j]) <= available
                   and fibres[j] & available}
        if smaller == support:
            return support
        support = smaller
    return set()


def support_catalogue(food, channels, molecules):
    """Definition 3.1: all nonempty food-generated supports, not just minimal ones."""
    indices = {x: i for i, x in enumerate(molecules)}
    cores, supports = [], []
    for support in subsets(range(len(channels))):
        available = closure(food, channels, support)
        if support and all(set(channels[j]) <= available for j in support):
            eligible = frozenset(indices[x] for x in available)
            cores.append({j: eligible for j in support})
            supports.append(support)
    return cores, supports


def write_csv(path, rows):
    with path.open('w', newline='', encoding='utf-8') as f:
        writer = csv.DictWriter(f, fieldnames=list(rows[0]))
        writer.writeheader()
        writer.writerows(rows)


def exact_checks(out):
    report, diagnostics = [], {}
    first = [frozenset(w) for w in ((0, 1), (0, 2), (0, 3))]
    second = [frozenset(w) for w in ((0, 1), (0, 2), (1, 2))]
    assert union_sizes(first) == {0: 1, 2: -3, 3: 3, 4: -1}
    assert union_sizes(second) == {0: 1, 2: -3, 3: 2}
    fixture_rows = []
    for name, family in [('I', first), ('II', second)]:
        assert [len(w) for w in family] == [2, 2, 2]
        assert [len(a & b) for a, b in combinations(family, 2)] == [1, 1, 1]
        local = local_counts(family, 4)
        direct = [0] * 5
        for mask in range(16):
            catalysts = {x for x in range(4) if mask & (1 << x)}
            all_hit = all(catalysts & w for w in family)
            any_hit = any(catalysts & w for w in family)
            direct[mask.bit_count()] += int(all_hit)
            fixture_rows.append(dict(family=name, mask=mask,
                                     catalysts=','.join(map(str, sorted(catalysts))),
                                     all_requirements_met=int(all_hit),
                                     at_least_one_met=int(any_hit)))
        assert local == direct
        cores = [{0: w} for w in family]
        union, counts = catalogue_law(cores, 4, 1, F(1, 2))
        assert counts == enumerate_catalogue(cores, 4, 1)
        assert union == probability_from_counts(counts, F(1, 2))
        diagnostics[name] = dict(hit_coefficients=union_sizes(family),
                                 all_counts_by_Q=local, any_counts_by_Q=counts,
                                 all_at_half=str(local_probability(family, F(1, 2))),
                                 any_at_half=str(union))
        report.append(f'Family {name}: all requirements at p=1/2 = '
                      f'{local_probability(family, F(1, 2))}; at least one = {union}.')
    write_csv(out / 'catalyst_subsets.csv', fixture_rows)
    report.append('Exact polynomial difference: q^3 - q^4 = p(1-p)^3.')
    report.append('Its maximum is 27/256 at p=1/4 (differentiate p(1-p)^3).')

    # Separate channels give independent events under Bernoulli sampling only.
    dedicated = [{j: w} for j, w in enumerate(first)]
    independent, counts = catalogue_law(dedicated, 4, 3, F(1, 2))
    assert independent == 1 - (1 - F(3, 4)) ** 3 == F(63, 64)
    assert counts == enumerate_catalogue(dedicated, 4, 3)
    fixed_two = F(counts[2], math.comb(12, 2))
    single_fixed = 1 - F(math.comb(10, 2), math.comb(12, 2))
    false_product = 1 - (1 - single_fixed) ** 3
    assert fixed_two != false_product
    report.append(f'Dedicated channels: Bernoulli union = {independent}; '
                  f'at fixed Q=2, union = {fixed_two}, product of fixed-Q marginals '
                  f'would give {false_product} (not equal).')

    mixed = [{0: first[0], 1: first[1]}, {0: first[1]},
             {0: first[0], 1: first[2]}]
    mixed_p, mixed_counts = catalogue_law(mixed, 4, 2, F(1, 3))
    assert mixed_counts == enumerate_catalogue(mixed, 4, 2)
    assert mixed_p == probability_from_counts(mixed_counts, F(1, 3))
    # Empty catalogue, unconstrained core, unused channel, and p endpoints.
    for cores in [[], [{}], [{0: frozenset({0})}], [{0: frozenset()}]]:
        for p in [F(0), F(1, 3), F(1)]:
            probability, c = catalogue_law(cores, 2, 2, p)
            assert c == enumerate_catalogue(cores, 2, 2)
            assert probability == probability_from_counts(c, p)
    report.append('Both inclusion-exclusion levels agree with literal enumeration; '
                  'fixed-Q coefficients agree for every Q in each checked catalogue.')

    census = []
    for n in range(2, 11):
        channels = repository_channels(n)
        # Independent displayed-reaction deduplication checks reaction identity.
        displayed = {(tuple(sorted((w[:s], w[s:]))), w)
                     for w in words(n) for s in range(1, len(w))}
        assert len(channels) == len(displayed) == channel_count(n)
        seeds = [(u, v, w) for u, v, w in channels if len(u) <= 2 and len(v) <= 2]
        assert len(seeds) == (4 if n == 2 else 18 if n == 3 else 34)
        split = (n - 2) * (1 << (n + 1)) + 4
        census.append(dict(n=n, molecules=(1 << (n + 1)) - 2,
                           channels=len(channels), directed_records=2 * len(channels),
                           split_position_channels=split, removed_duplicates=split-len(channels),
                           food_accessible_channels=len(seeds)))
    write_csv(out / 'reaction_counts.csv', census)
    seeds = [(u, v, w) for u, v, w in repository_channels(4) if len(u) <= 2 and len(v) <= 2]
    write_csv(out / 'food_channels.csv', [dict(index=i, left_1=u, left_2=v, product=w)
                                         for i, (u, v, w) in enumerate(seeds)])
    report.append('Channel counts, n=2..10: ' + ', '.join(str(r['channels']) for r in census))
    report.append('Food-accessible counts, n=2..10: ' + ', '.join(str(r['food_accessible_channels']) for r in census))

    # Actual reversible RAF semantics on a restricted n=4 repository network.
    # It has all 30 ambient molecules, but only these two reaction channels.
    channels = [('A', 'AA', 'AAA'), ('A', 'AAA', 'AAAA')]
    food = set(words(2))
    relevant = words(2) + ['AAA', 'AAAA']
    ambient = words(4)
    cores, supports = support_catalogue(food, channels, relevant)
    assert supports == [(0,), (0, 1)]
    counts = [0] * 17
    for mask in range(1 << 16):
        fibres = [{relevant[x] for x in range(8) if mask & (1 << (8 * j + x))}
                  for j in range(2)]
        literal = any(is_raf(food, channels, support, fibres)
                      for support in subsets(range(2)))
        indexed = [{i for i, x in enumerate(relevant) if x in f} for f in fibres]
        assert literal == any(activated(core, indexed) for core in cores)
        assert literal == bool(max_raf(food, channels, fibres))
        if literal:
            assert fibres[0]  # The only food-accessible channel must be catalysed.
            counts[mask.bit_count()] += 1
    # The other 22 molecules cannot be made by these channels. Their 44
    # catalysis assignments are unrestricted and still count in fixed-Q sampling.
    full_counts = multiply(counts, [math.comb(44, q) for q in range(45)])
    ambient_cores, _ = support_catalogue(food, channels, ambient)
    probability, exact_full_counts = catalogue_law(ambient_cores, 30, 2, F(1, 2))
    assert full_counts == exact_full_counts
    assert probability == probability_from_counts(counts, F(1, 2)) == F(65279, 65536)
    for p in [F(0), F(1, 10), F(1, 2), F(1)]:
        exact, _ = catalogue_law(ambient_cores, 30, 2, p)
        q = 1 - p
        assert exact == 1 - q**7 + q**7 * p * (1 - q**8)
    # A catalysed entry can fail to support any RAF; mutual support can rescue it.
    assert not max_raf(food, channels, [{'AAAA'}, set()])
    assert max_raf(food, channels, [{'AAAA'}, {'AAA'}]) == {0, 1}
    report.append('Restricted n=4 system: checked all 65,536 relevant assignments; '
                  'RAF definition, catalogue activation, and pruning agree.')
    report.append(f'Its exact RAF probability at p=1/2 is {probability}; '
                  'all 61 fixed-Q coefficients on the 60 ambient coordinates agree.')
    report.append('A catalysed entry without a RAF, and a two-channel RAF, both verified.')
    write_csv(out / 'restricted_network_fixed_Q.csv',
              [dict(Q=q, successful_assignments=c, total_assignments=math.comb(60, q),
                    exact_probability=str(F(c, math.comb(60, q))))
               for q, c in enumerate(full_counts)])
    diagnostics['dedicated'] = dict(bernoulli_half=str(independent),
                                    fixed_Q2=str(fixed_two), invalid_product=str(false_product))
    diagnostics['restricted_network'] = dict(molecules=ambient, channels=channels,
                                             food=sorted(food), raf_at_half=str(probability),
                                             assignments_checked=65536)
    return report, diagnostics


def make_figures(out):
    import numpy as np
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    from matplotlib.colors import LinearSegmentedColormap

    plt.rcParams.update({'font.family': 'DejaVu Sans', 'font.size': 10,
                         'axes.spines.top': False, 'axes.spines.right': False,
                         'axes.labelcolor': '#263641', 'text.color': '#263641',
                         'figure.facecolor': '#fcfbf7', 'axes.facecolor': '#fcfbf7',
                         'savefig.facecolor': '#fcfbf7', 'axes.titleweight': 'bold'})
    teal, coral, blue = '#007f7a', '#c45440', '#4764a4'

    def save(fig, name):
        fig.savefig(out / f'{name}.png', dpi=190, bbox_inches='tight')
        fig.savefig(out / f'{name}.svg', bbox_inches='tight')
        plt.close(fig)

    fig = plt.figure(figsize=(12, 7.7), layout='constrained')
    grid = fig.add_gridspec(2, 2, height_ratios=[1, 1.7])
    families = [((0, 1), (0, 2), (0, 3)), ((0, 1), (0, 2), (1, 2))]
    for index, (family, color) in enumerate(zip(families, [teal, coral])):
        ax = fig.add_subplot(grid[0, index])
        for row, w in enumerate(family):
            for column in range(4):
                ax.scatter(column, row, s=500, color=color if column in w else '#e7e7e0')
                if column in w:
                    ax.text(column, row, '✓', ha='center', va='center', color='white', fontsize=14)
        ax.set(xticks=range(4), xticklabels=['0', '1', '2', '3'], yticks=range(3),
               yticklabels=['Requirement 1', 'Requirement 2', 'Requirement 3'],
               xlim=(-.65, 3.65), ylim=(2.6, -.7), xlabel='Possible catalyst molecule')
        ax.set_title(f'Family {"I" if index == 0 else "II"}: {"9/16" if index == 0 else "8/16"} at p = 1/2', pad=10)
        for spine in ax.spines.values():
            spine.set_visible(False)
        ax.tick_params(length=0)
    p = np.linspace(0, 1, 401)
    q = 1 - p
    first, second = 1 - 3*q*q + 3*q**3 - q**4, 1 - 3*q*q + 2*q**3
    ax = fig.add_subplot(grid[1, 0])
    ax.plot(p, first, color=teal, label='Family I', lw=2.5)
    ax.plot(p, second, color=coral, label='Family II', lw=2.5, ls='--')
    ax.fill_between(p, second, first, color=teal, alpha=.12)
    ax.set(xlabel='Probability p of each catalysis assignment', ylabel='Probability ALL three requirements are met', ylim=(0, 1.02))
    ax.legend(frameon=False)
    ax.grid(alpha=.15)
    ax = fig.add_subplot(grid[1, 1])
    ax.plot(p, p*q**3, color=blue, lw=2.5)
    ax.scatter([.25, .5], [27/256, 1/16], color=blue)
    ax.annotate('Largest gap: 27/256 ≈ 10.55 percentage points', (.25, 27/256),
                xytext=(.03, .127), fontsize=9, arrowprops={'arrowstyle': '-', 'color': blue})
    ax.annotate('At p = 1/2: 1/16 = 6.25 points', (.5, 1/16), xytext=(.42, .084), fontsize=9)
    ax.set(xlabel='Probability p of each catalysis assignment', ylabel='Family I minus Family II', ylim=(0, .145))
    ax.grid(alpha=.15)
    fig.suptitle('Joint catalyst-support probabilities for two overlap patterns', fontsize=18)
    save(fig, 'overlap_comparison')
    write_csv(out / 'overlap_curves.csv', [dict(p=float(x), all_I=float(a), all_II=float(b),
                                               gap=float(a-b), any_I=float(1-(1-x)**4),
                                               any_II=float(1-(1-x)**3),
                                               independent_union=float(1-(1-x)**6))
                                           for x, a, b in zip(p, first, second)])

    ns = sorted(set(np.rint(np.geomspace(4, 4096, 220)).astype(int).tolist()) | {6, 10, 64, 128, 512, 1024, 4096})
    rows = []
    fig, axes = plt.subplots(1, 2, figsize=(12.5, 5.5), layout='constrained')
    for f, color in zip([.3, 1.3, 3.0], [teal, coral, blue]):
        g = [gateway_probability(n, f) for n in ns]
        axes[0].plot(ns, g, color=color, label=f'f = {f:g}', lw=2.3)
        rows += [dict(mode='fixed_f', value=f, n=n, f=f, gateway_probability=v, limit=0.0) for n, v in zip(ns, g)]
    axes[0].set(xscale='log', xlabel='Maximum molecule length n (log scale)', ylabel='Upper bound on RAF probability',
                title='Fixed f: the bound tends to zero', ylim=(0, 1.03))
    for rate, color in zip([.01, .03, .1], [teal, coral, blue]):
        g = [gateway_probability(n, rate*n) for n in ns]
        limit = -math.expm1(-34*rate)
        axes[1].plot(ns, g, color=color, label=f'λ = {rate:g}; limit {limit:.4f}', lw=2.3)
        axes[1].axhline(limit, color=color, linestyle=':', alpha=.75)
        rows += [dict(mode='linear_f', value=rate, n=n, f=rate*n, gateway_probability=v, limit=limit) for n, v in zip(ns, g)]
    axes[1].set(xscale='log', xlabel='Maximum molecule length n (log scale)',
                title='f = λn: the bound approaches 1 − exp(−34λ)', ylim=(0, 1.03))
    for ax in axes:
        ax.grid(alpha=.16)
        ax.legend(frameon=False, fontsize=9)
    fig.suptitle('Only 34 channels can be used directly from food', fontsize=18)
    save(fig, 'gateway_scaling')
    write_csv(out / 'gateway_curves.csv', rows)

    fs = np.geomspace(.01, 4, 150)
    z = np.array([[gateway_probability(n, float(f)) for n in ns] for f in fs])
    fig, ax = plt.subplots(figsize=(11, 6.5), layout='constrained')
    cmap = LinearSegmentedColormap.from_list('gateway', ['#fcf8ed', '#a8d4c5', '#267f85', '#173c5b'])
    mesh = ax.pcolormesh(ns, fs, z, shading='nearest', cmap=cmap, vmin=0, vmax=1, rasterized=True)
    contour = ax.contour(ns, fs, z, levels=[.1, .5, .9], colors=['#283d45'], linewidths=1)
    labels = ax.clabel(contour, fmt={.1: '10% bound', .5: '50% bound', .9: '90% bound'},
                      fontsize=9, inline=False, manual=[(30, .087), (40, .775), (35, 2.23)])
    for label in labels:
        label.set_rotation(0)
        label.set_bbox(dict(facecolor='#fcfbf7', edgecolor='none', alpha=.85, pad=3))
    ax.set(xscale='log', yscale='log', xlabel='Maximum molecule length n (log scale)',
           ylabel='Expected channels catalysed per molecule f (log scale)',
           title='Upper bound map: colour shows a catalysed entry, not a RAF')
    bar = fig.colorbar(mesh, ax=ax, pad=.025)
    bar.set_label('Gateway-open probability Gₙ')
    save(fig, 'gateway_bound_map')
    write_csv(out / 'gateway_map.csv', [dict(n=n, f=float(f), gateway_probability=float(z[i, j]))
                                       for i, f in enumerate(fs) for j, n in enumerate(ns)])
    return rows


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--output', type=Path, default=Path('results'))
    args = parser.parse_args()
    out = args.output.resolve()
    out.mkdir(parents=True, exist_ok=True)
    report, diagnostics = exact_checks(out)
    make_figures(out)
    report.append('Gateway bound at f=1.3:')
    for n in [6, 10, 64, 128, 512, 1024, 4096]:
        report.append(f'  n={n:4d}: {gateway_probability(n, 1.3):.9f}')
    report.append('Linear-scale example λ=0.03:')
    for n in [6, 64, 4096]:
        report.append(f'  n={n:4d}: {gateway_probability(n, .03*n):.9f}')
    report.append(f'  limiting gateway bound: {-math.expm1(-1.02):.9f}')
    report.append('All exact finite checks passed. General asymptotic claims rely on the paper.')
    source = Path(__file__).read_text(encoding='utf-8')
    diagnostics['source_sha256'] = hashlib.sha256(source.encode()).hexdigest()
    diagnostics['method'] = 'Exact integers/Fractions for finite checks; double precision for figures; no randomness.'
    (out / 'diagnostics.json').write_text(json.dumps(diagnostics, indent=2) + '\n')
    (out / 'diagnostics.txt').write_text('\n'.join(report) + '\n')
    print('\n'.join(report))
    text = '# Exact overlap and food-entry bounds: code and saved outputs\n\nReactions can depend on overlapping pools of possible catalysts. This program calculates how that overlap affects the probability that all reactions receive catalytic support. Finite probability checks use exact integers and fractions. Curves numerically evaluate\nthe paper\'s finite formulas; they are not sampled RAF frequencies.\n\n## Catalyst overlap and joint support probability\n\n![Two catalyst-set families and their exact probability difference](overlap_comparison.svg)\n\nEach row requires at least one catalyst from its marked molecules. Both families\nhave two molecules per row and one shared molecule per pair of rows. When each molecule is independently available with probability p=1/2,\nall three rows are satisfied in 9 of 16 assignments for Family I, and 8 of 16\nfor Family II. This is a joint event. The different event "at least one row is\nsatisfied" has probabilities 15/16 and 7/8. The small example is an abstract\nrequirements example from Theorem 6.1, not a complete binary-polymer RAF network.\n\n## Fixed food constrains the growth of the model\n\n![Gateway probability for fixed and linearly growing catalysis](gateway_scaling.svg)\n\nThe food set remains A, B, AA, AB, BA, BB. A channel is a reversible reaction;\none catalysis assignment applies to both directions. There are exactly 34\nchannels with a reaction side made entirely of food once n is at least 4.\nEvery RAF must include a catalysed such channel. The plots give the probability\nthat at least one of these channels has any catalyst in the whole molecule set.\nThis is an upper bound on RAF probability, since that catalyst might not be\nreachable and the remaining required reactions might not have catalytic support.\n\nHere f is the mean number of channels catalysed by one molecule. At fixed f, the bound tends to zero. At f=λn, it tends to 1−exp(−34λ), which\nis below one for every finite λ. This paper does not determine the RAF limit.\nThe n axis measures maximum word length, not number of molecule types.\n\n![Map of the gateway upper bound](gateway_bound_map.svg)\n\nThis is an upper-bound map, not a phase diagram of actual RAF existence. For\nexample, a 10% contour says that RAF probability is at most 10% there. A bound\nnear one does not establish that RAFs are common. Full numerical values are in gateway_map.csv.\n\n## Actual RAF check in a small restricted network\n\nKeep the six food molecules and all 30 ambient words of length at most four,\nbut allow only A + AA ⇌ AAA and A + AAA ⇌ AAAA. The two nonempty\nfood-generated supports are the first channel alone and both channels together.\nThe script compares the closure-based RAF definition, the exhaustive support\ncatalogue, and iterative pruning for all 65,536 assignments involving reachable\nmolecules. The 44 remaining catalysis coordinates cannot support this restricted\nnetwork; their possible assignments are nevertheless counted in fixed-Q results.\n\nIf only AAAA catalyses the first channel, an entry is catalysed but no RAF\nexists. If AAA also catalyses the second channel, both channels form a RAF.\nClosure is computed without demanding a catalyst at every construction step;\nthe RAF condition is not a claim of an entirely catalysed startup sequence.\n\nFor q=1−p the exact RAF probability of this restricted network is\n1−q⁷+q⁷p(1−q⁸). At p=1/2 it is 65279/65536. These are not probabilities\nfor the full 64-channel network at n=4.\n\n## Saved diagnostic output\n\n```text\nFamily I: all requirements at p=1/2 = 9/16; at least one = 15/16.\nFamily II: all requirements at p=1/2 = 1/2; at least one = 7/8.\nExact polynomial difference: q^3 - q^4 = p(1-p)^3.\nIts maximum is 27/256 at p=1/4 (differentiate p(1-p)^3).\nDedicated channels: Bernoulli union = 63/64; at fixed Q=2, union = 17/22, product of fixed-Q marginals would give 7273/10648 (not equal).\nBoth inclusion-exclusion levels agree with literal enumeration; fixed-Q coefficients agree for every Q in each checked catalogue.\nChannel counts, n=2..10: 4, 18, 64, 188, 502, 1264, 3048, 7130, 16334\nFood-accessible counts, n=2..10: 4, 18, 34, 34, 34, 34, 34, 34, 34\nRestricted n=4 system: checked all 65,536 relevant assignments; RAF definition, catalogue activation, and pruning agree.\nIts exact RAF probability at p=1/2 is 65279/65536; all 61 fixed-Q coefficients on the 60 ambient coordinates agree.\nA catalysed entry without a RAF, and a two-channel RAF, both verified.\nGateway bound at f=1.3:\n  n=   6: 0.999985014\n  n=  10: 0.996060576\n  n=  64: 0.509781087\n  n= 128: 0.295870965\n  n= 512: 0.083017294\n  n=1024: 0.042326652\n  n=4096: 0.010738217\nLinear-scale example λ=0.03:\n  n=   6: 0.784839513\n  n=  64: 0.651076726\n  n=4096: 0.639584696\n  limiting gateway bound: 0.639405060\nAll exact finite checks passed. General asymptotic claims rely on the paper.\n```\n\nThe fixed-Q tables count assignments directly; their probabilities are exact\nfractions. Fixed-Q sampling does not make even disjoint channel events\nindependent. The independent-core product is checked under Bernoulli sampling.\n\nThe manuscript reports Lean verification except for Proposition 2.1 and\nCorollary 9.6. This program does not compile Lean or audit those proof files.\nIt checks the stated finite examples and evaluates the gateway formulas.\n\n## Run and inspect\n\nInstall numpy and matplotlib, then run:\n\n```bash\npython -X utf8 overlap_raf.py --output results\n```\n\nThe program needs no external data, random seed, account, or API key.\nCSV files contain all plotted data, the reaction census, the 34 entry channels,\nthe 16 catalyst subsets for each family, and every fixed-Q result for the\nrestricted network. SVG copies of the figures support publication at any size.\ndiagnostics.json records the source hash. Inclusion-exclusion is exponential;\nthe script deliberately applies it only to small catalogues.\n\n\n[Download the complete package](overlap_package.zip). The source below includes the revised chart title; its scientific calculations are unchanged. UTF-8 mode is needed on Windows.\n\nThe exact channel count gives the sharper coefficient $1+o(1)$ after Proposition 9.1. The bound $272f_n/n$ and limiting constant 34 are unaffected. Theorem 7.2 supplies a sufficient coordinate-disjointness hypothesis for the Bernoulli product.\n\n[Download the manuscript TeX source](manuscript_source.zip).\n\n## Complete Python source\n\n```python\n' + source + '\n```\n'
    (out / 'Code_and_Outputs.md').write_text(text, encoding='utf-8')


if __name__ == '__main__':
    main()

```
