{
"cells": [
{
"cell_type": "markdown",
"id": "6c1a6ba1-4ea4-4400-a4f2-5bf083577896",
"metadata": {},
"source": [
"# Online Stochastic Matchings: Stability on Hypergraphs"
]
},
{
"cell_type": "markdown",
"id": "cell-0",
"metadata": {},
"source": [
"This is the companion notebook of the paper *Online Stochastic Matchings: Stability on Hypergraphs*.\n",
"\n",
"It covers the numerical claims of the paper, focusing on the candy hypergraph:\n",
"\n",
"1. **The characterization of stability.** `model.stabilizable` decides any $(G,\\lambda)$; `model.maximin`\n",
" gives the witness flow $\\mu$; `model.incidence` gives the matrix.\n",
"2. **The candy's closed-form stability region**, checked against `model.stabilizable`.\n",
"3. **Regular greedy ML vs. VQML** on the candy $\\alpha$-family (the paper's figure).\n",
"4. **The provable greedy-instability threshold** $\\alpha<2/21$: the domination and\n",
" throughput-conservation ingredients.\n",
"5. **Accessibility and tie-breaking of the virtual chain**: proof-level combinatorics."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "cell-1",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-18T22:02:30.298352Z",
"iopub.status.busy": "2026-07-18T22:02:30.297878Z",
"iopub.status.idle": "2026-07-18T22:02:32.024621Z",
"shell.execute_reply": "2026-07-18T22:02:32.023566Z",
"shell.execute_reply.started": "2026-07-18T22:02:30.298327Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"stochastic_matching: 0.4.0\n"
]
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import multiprocess as mp\n",
"import stochastic_matching as sm\n",
"from stochastic_matching import XP, Iterator, evaluate\n",
"\n",
"print(\"stochastic_matching:\", sm.__version__)\n",
"np.set_printoptions(precision=4, suppress=True)"
]
},
{
"cell_type": "markdown",
"id": "60101877-3f1a-47a1-85ee-34cfb99a39c1",
"metadata": {},
"source": [
"First, recall the candy structure: two triangles bridged by a hyperedge that includes a seventh node:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "ed3104d6-c2e0-4c29-bf58-3c574a38faa6",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-18T22:02:39.840531Z",
"iopub.status.busy": "2026-07-18T22:02:39.840248Z",
"iopub.status.idle": "2026-07-18T22:02:39.848673Z",
"shell.execute_reply": "2026-07-18T22:02:39.847447Z",
"shell.execute_reply.started": "2026-07-18T22:02:39.840512Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"candy incidence A (from model.incidence):\n",
"[[1 1 0 0 0 0 0]\n",
" [1 0 1 0 0 0 0]\n",
" [0 1 1 0 0 0 1]\n",
" [0 0 0 0 0 0 1]\n",
" [0 0 0 1 1 0 1]\n",
" [0 0 0 1 0 1 0]\n",
" [0 0 0 0 1 1 0]] \n",
"\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"candy = sm.HyperPaddle()\n",
"print(\"candy incidence A (from model.incidence):\")\n",
"print(np.asarray(candy.incidence), \"\\n\")\n",
"candy.show_graph()"
]
},
{
"cell_type": "markdown",
"id": "e4dfd7e2-9f50-4368-bde4-6af2960a677a",
"metadata": {},
"source": [
"## Characterization of stability"
]
},
{
"cell_type": "markdown",
"id": "cell-2",
"metadata": {},
"source": [
"The package decides stabilizability directly: `model.stabilizable` is `True` iff the maximin\n",
"conservation flow is strictly positive **and** the incidence matrix is surjective, i.e.\n",
"conditions (iii)/(iv) of the main theorem. When stabilizable, `model.maximin` *is* a positive\n",
"witness $\\mu$ with $A\\mu=\\lambda$; when not, either the flow has a non-positive coordinate\n",
"(the cone is left) or $A$ is rank-deficient.\n",
"\n",
"> Note: The dual *certificate* $y$ of condition (ii) is stated in the paper; the package works directly with the primal witness."
]
},
{
"cell_type": "markdown",
"id": "f3cf04d4-6d54-4125-9fcf-fbaf346708df",
"metadata": {},
"source": [
"We will use the `report` function to observe a few cases."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "a7f829fb-a736-4852-9ea7-9c2561aa3914",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-18T22:02:44.606857Z",
"iopub.status.busy": "2026-07-18T22:02:44.606547Z",
"iopub.status.idle": "2026-07-18T22:02:44.612981Z",
"shell.execute_reply": "2026-07-18T22:02:44.611161Z",
"shell.execute_reply.started": "2026-07-18T22:02:44.606839Z"
}
},
"outputs": [],
"source": [
"def report(name, model):\n",
" A = np.asarray(model.incidence)\n",
" mu = np.asarray(model.maximin)\n",
" rank, n = np.linalg.matrix_rank(A), A.shape[0]\n",
" print(f\"{name}\")\n",
" print(f\" stabilizable = {model.stabilizable} \"\n",
" f\"(rank A = {rank} of {n}; maximin min = {mu.min():.3f})\")\n",
" print(f\" maximin flow mu = {mu.round(4)}\")"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "55893b82-b05d-4988-8aa2-29e8ca9eebe8",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-18T22:02:47.960013Z",
"iopub.status.busy": "2026-07-18T22:02:47.959636Z",
"iopub.status.idle": "2026-07-18T22:02:47.967683Z",
"shell.execute_reply": "2026-07-18T22:02:47.966665Z",
"shell.execute_reply.started": "2026-07-18T22:02:47.959992Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"candy, alpha-family (alpha=0.3): lambda=(1,1,3a,a,3a,1,1)\n",
" stabilizable = True (rank A = 7 of 7; maximin min = 0.300)\n",
" maximin flow mu = [0.7 0.3 0.3 0.3 0.3 0.7 0.3]\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"candy.rates = [1, 1, 0.9, 0.3, 0.9, 1, 1]\n",
"report(\"candy, alpha-family (alpha=0.3): lambda=(1,1,3a,a,3a,1,1)\", candy)\n",
"candy.show_flow()"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "4f02772e-729a-4d76-8d9c-d5ce8cc2c05a",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-18T22:02:54.138175Z",
"iopub.status.busy": "2026-07-18T22:02:54.137779Z",
"iopub.status.idle": "2026-07-18T22:02:54.145303Z",
"shell.execute_reply": "2026-07-18T22:02:54.144215Z",
"shell.execute_reply.started": "2026-07-18T22:02:54.138157Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"candy, center overloaded λ_4=5\n",
" stabilizable = False (rank A = 7 of 7; maximin min = -2.000)\n",
" maximin flow mu = [ 3. -2. -2. -2. -2. 3. 5.]\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"candy.rates = [1, 1, 1, 5, 1, 1, 1]\n",
"report(\"candy, center overloaded λ_4=5\", candy)\n",
"candy.show_flow()"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "e00bcca7-8b8d-4950-bb06-f32725cecba4",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-18T22:03:01.455931Z",
"iopub.status.busy": "2026-07-18T22:03:01.455649Z",
"iopub.status.idle": "2026-07-18T22:03:01.466339Z",
"shell.execute_reply": "2026-07-18T22:03:01.464955Z",
"shell.execute_reply.started": "2026-07-18T22:03:01.455914Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"rem:wall degenerate (stabilizable; every basic solution degenerate)\n",
" stabilizable = True (rank A = 3 of 3; maximin min = 0.500)\n",
" maximin flow mu = [0.5 0.5 0.5 0.5]\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"model = sm.Model(incidence=[[0, 1, 1, 0], [0, 0, 1, 1], [1, 1, 1, 1]], rates=[1, 1, 2]) \n",
"report(\"rem:wall degenerate (stabilizable; every basic solution degenerate)\", model)\n",
"model.show_kernel(disp_flow=True)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "bf89300c-dbfd-4e98-a1f5-24964f68dbb4",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-18T22:03:05.454129Z",
"iopub.status.busy": "2026-07-18T22:03:05.453801Z",
"iopub.status.idle": "2026-07-18T22:03:05.462073Z",
"shell.execute_reply": "2026-07-18T22:03:05.460778Z",
"shell.execute_reply.started": "2026-07-18T22:03:05.454110Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Moyal single 3-hyperedge {1,2,3} (rank 1 < 3)\n",
" stabilizable = False (rank A = 1 of 3; maximin min = 1.000)\n",
" maximin flow mu = [1.]\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"model = sm.Model(incidence=[[1], [1], [1]], rates=[1, 1, 1])\n",
"report(\"Moyal single 3-hyperedge {1,2,3} (rank 1 < 3)\", model)\n",
"model.show_kernel(disp_flow=True)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "cell-3",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-18T22:03:10.486598Z",
"iopub.status.busy": "2026-07-18T22:03:10.486205Z",
"iopub.status.idle": "2026-07-18T22:03:10.494311Z",
"shell.execute_reply": "2026-07-18T22:03:10.493062Z",
"shell.execute_reply.started": "2026-07-18T22:03:10.486579Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Moyal RM21 Prop-4 flavor: {1,2,3}+{3,4}, classes 1,2 degree one\n",
" stabilizable = False (rank A = 2 of 4; maximin min = 1.000)\n",
" maximin flow mu = [1. 1.]\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"model = sm.Model(incidence=[[1, 0], [1, 0], [1, 1], [0, 1]], rates=[1, 1, 2, 1])\n",
"report(\"Moyal RM21 Prop-4 flavor: {1,2,3}+{3,4}, classes 1,2 degree one\", model)\n",
"model.show_kernel(disp_flow=True)"
]
},
{
"cell_type": "markdown",
"id": "c44a014b-b3cf-4071-9cf3-18ba7a38a5b3",
"metadata": {},
"source": [
"## The candy's closed-form stability region"
]
},
{
"cell_type": "markdown",
"id": "cell-4",
"metadata": {},
"source": [
"For the candy $A$ is square and invertible, so stabilizability $\\iff$ the unique solution of\n",
"$A\\mu=\\lambda$ is strictly positive, which solves to\n",
"$$|\\lambda_1-\\lambda_2| < \\lambda_3-\\lambda_4 < \\lambda_1+\\lambda_2 \\quad\\text{and}\\quad\n",
" |\\lambda_7-\\lambda_6| < \\lambda_5-\\lambda_4 < \\lambda_6+\\lambda_7.$$\n",
"\n",
"Recipe for the formula above:\n",
"- The 3-edge must be $\\lambda_4$, leaving $\\bar{\\lambda_3}=\\lambda_3-\\lambda_4$ and $\\bar{\\lambda_5}=\\lambda_5-\\lambda_4$ for nodes 3 and 5.\n",
"- That leaves the triangles $\\lambda_1$, $\\lambda_2$, $\\bar{\\lambda_3}$, and $\\lambda_6$, $\\lambda_7$, $\\bar{\\lambda_5}$, which must uphold the triangular inequality.\n",
"- Triangular inequality for $a$, $b$, $c$: $|b-c| ϵ:\n",
" alphas = [αmin + (αmax - αmin)*i/(span+1) for i in range(1, span+1)]\n",
" xp = XP('ML', simulator='longest', \n",
" iterator=Iterator('model', alphas, name='alpha', process=candy_at), **common)\n",
" res = evaluate(xp, ['steps_done'], pool)\n",
" if res['ML']['steps_done'][0] == steps_done:\n",
" αmax = alphas[0]\n",
" else:\n",
" for i in range(1, span):\n",
" if res['ML']['steps_done'][i] == steps_done:\n",
" αmin, αmax = alphas[i-1], alphas[i]\n",
" break\n",
" else:\n",
" αmin = alphas[-1]\n",
" α0 = (αmin+αmax)/2\n",
"α0"
]
},
{
"cell_type": "markdown",
"id": "ea674508-689a-495b-b681-adacf16db910",
"metadata": {},
"source": [
"We then sweep $\\alpha$."
]
},
{
"cell_type": "markdown",
"id": "cell-6",
"metadata": {},
"source": [
"We plot two metrics side by side: the mean central queue $\\overline{Q_4}$ (the coordinate the\n",
"instability proof bounds, and the paper's figure) and the package's built-in `delay` (mean\n",
"waiting time by Little's law, $\\sum_i\\overline{Q_i}/\\Lambda$). \n",
"\n",
"Both show the same thing on the left: greedy blows up as $\\alpha\\to \\alpha_0$ while VQML goes to 0, so the separation does not hinge on singling out node~4. \n",
"\n",
"At large $\\alpha$, delay and $\\overline{Q_4}$ differ (the instability comes from outer nodes).\n",
"\n",
"Also note that for delay and large $\\alpha$, the curves cross: VQML carries slightly more global delay than\n",
"greedy. It is understandable, since VQML is a maximally *stable* witness through a reservation mechanism, not a delay optimizer."
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "b3de364b-bcb3-49fe-891e-77cd145e95e6",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-18T22:04:09.362459Z",
"iopub.status.busy": "2026-07-18T22:04:09.362224Z",
"iopub.status.idle": "2026-07-18T22:07:03.662474Z",
"shell.execute_reply": "2026-07-18T22:07:03.661032Z",
"shell.execute_reply.started": "2026-07-18T22:04:09.362442Z"
}
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "4374942f37ea4da09a945e601fc376f1",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/154 [00:00"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig, axes = plt.subplots(1, 2, figsize=(9.4, 3.4), sharex=True)\n",
"panels = [('central_queue', r'mean central queue $\\overline{Q_4}$'),\n",
" ('delay', r'delay (mean waiting time)')]\n",
"\n",
"\n",
"for ax, (key, ylabel) in zip(axes, panels):\n",
" # yml, yvq = np.array(ml[key]), np.array(vq[key])\n",
" ax.semilogy(res['ML']['alpha'], res['ML'][key], ms=4, label='ML (greedy)')\n",
" ax.semilogy(res['VQML']['alpha'], res['VQML'][key], ms=5, label='VQML')\n",
" ax.axvline(2 / 21, color='gray', ls=':', lw=1)\n",
" ax.text(2 / 21, 0.03, ' provable 2/21', rotation=90, va='bottom', fontsize=8,\n",
" transform=ax.get_xaxis_transform())\n",
" ax.axvline(α0, color='gray', ls=':', lw=1)\n",
" ax.text(α0, 0.03, '$\\\\alpha_0$', rotation=90, va='bottom', fontsize=8,\n",
" transform=ax.get_xaxis_transform())\n",
" ax.set_xlabel(r'$\\alpha$'); ax.set_ylabel(ylabel)\n",
" ax.set_xlim(0, 1); ax.grid(True, which='both', alpha=0.2)\n",
"axes[0].legend()\n",
"plt.tight_layout(); plt.show()"
]
},
{
"cell_type": "markdown",
"id": "cell-8",
"metadata": {},
"source": [
"## The provable greedy-instability threshold $\\alpha<2/21$\n",
"\n",
"For $\\alpha<2/21\\approx0.095$ **no** greedy policy is stable. The proof combines a\n",
"triangle/hyperedge invariant, an $M/M/1$ domination giving $\\Pr(Q_3>0)\\le 3\\alpha/2$, and the\n",
"throughput identity $r_h=\\lambda_4$. We corroborate the two ingredients where ML is stable\n",
"($\\alpha\\ge0.5$)."
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "cell-9",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-18T22:07:04.310640Z",
"iopub.status.busy": "2026-07-18T22:07:04.310303Z",
"iopub.status.idle": "2026-07-18T22:07:16.834326Z",
"shell.execute_reply": "2026-07-18T22:07:16.833342Z",
"shell.execute_reply.started": "2026-07-18T22:07:04.310613Z"
}
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "0ab65b62efbf482d904cbd49d9598a94",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/4 [00:000)=0.220 <= 3a/2=0.750? True; rate(h)=0.4960 vs lambda_4=0.5000 (conservation)\n",
"alpha=0.60: P(Qbridge>0)=0.345 <= 3a/2=0.900? True; rate(h)=0.5949 vs lambda_4=0.6000 (conservation)\n",
"alpha=0.70: P(Qbridge>0)=0.489 <= 3a/2=1.050? True; rate(h)=0.6947 vs lambda_4=0.7000 (conservation)\n",
"alpha=0.80: P(Qbridge>0)=0.651 <= 3a/2=1.200? True; rate(h)=0.7928 vs lambda_4=0.8000 (conservation)\n",
"\n",
"Provable threshold 2/21 = 0.0952 ; simulated ML threshold ~0.45.\n"
]
}
],
"source": [
"# Self-contained metrics (run in pool workers). Bridge = node 3 (0-based 2),\n",
"# hyperedge {3,4,5} = edge index 6.\n",
"def p_bridge(simu): # P(Q_bridge > 0) by PASTA\n",
" return 1.0 - simu.logs.queue_log[2, 0] / simu.logs.steps_done\n",
"\n",
"def rate_h_time(simu): # hyperedge firings per unit time\n",
" return float(simu.logs.traffic[6] / simu.logs.steps_done * sum(simu.model.rates))\n",
"\n",
"xp4 = XP('greedy', simulator='longest', model=None, n_steps=1_000_000, seed=1, max_queue=4000,\n",
" iterator=Iterator('model', [0.5, 0.6, 0.7, 0.8], name='alpha', process=candy_at))\n",
"with mp.Pool() as pool:\n",
" r4 = evaluate(xp4, [p_bridge, rate_h_time], pool)['greedy']\n",
"\n",
"for a, pb, rh in zip(r4['alpha'], r4['p_bridge'], r4['rate_h_time']):\n",
" print(f\"alpha={a:.2f}: P(Qbridge>0)={pb:.3f} <= 3a/2={1.5*a:.3f}? {pb <= 1.5*a + 1e-3}; \"\n",
" f\"rate(h)={rh:.4f} vs lambda_4={a:.4f} (conservation)\")\n",
"print(\"\\nProvable threshold 2/21 =\", round(2 / 21, 4), \"; simulated ML threshold ~0.45.\")"
]
},
{
"cell_type": "markdown",
"id": "cell-10",
"metadata": {},
"source": [
"## Accessibility and tie-breaking of the virtual chain\n",
"\n",
"Two structural facts about the signed virtual queue $Q$. These are *proof-level* combinatorics\n",
"(a custom adversarial tie-break rule, and a deterministic steering word) that the package's\n",
"`virtual_queue` simulator does not expose as knobs, so we check them directly with numpy.\n",
"\n",
"**(a) The $I_2$ trap (tie-breaking remark).** Take $A=I_2$ (two mono-edges), $\\lambda=(1,1)$,\n",
"and a legal *deterministic* maximizer rule that violates the idle clause (among maximizers,\n",
"prefer $2e_1, e_1, 2e_2, 0, e_2, e_1+e_2$). The set reachable from $0$ then has **9** states,\n",
"with $0$ visited exactly once and a **7**-state absorbing class. So the idle clause is what the\n",
"accessibility argument needs.\n",
"\n",
"**(b) Steering to $0$ (accessibility lemma).** Under the canonical idle clause, the Appendix-A\n",
"steering drives $Q$ to $0$ from every state within $\\lVert q\\rVert_1 + 2 a_{\\max} V(q)$. We run\n",
"it from many mixed-sign start states over a suite of incidence matrices (from the package where\n",
"available)."
]
},
{
"cell_type": "code",
"execution_count": 28,
"id": "cell-11",
"metadata": {
"execution": {
"iopub.execute_input": "2026-07-18T22:07:16.835319Z",
"iopub.status.busy": "2026-07-18T22:07:16.834876Z",
"iopub.status.idle": "2026-07-18T22:07:17.133608Z",
"shell.execute_reply": "2026-07-18T22:07:17.132635Z",
"shell.execute_reply.started": "2026-07-18T22:07:16.835295Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(a) I2 trap: reachable from 0 = 9 (paper 9); 0 re-enterable = False (paper False); absorbing SCC sizes = [7] (paper [7])\n",
"(b) steering: 1000/1000 start states over the suite reach 0 within the bound\n"
]
}
],
"source": [
"import itertools\n",
"\n",
"# (a) adversarial I2 tie-break trap\n",
"A_I2 = np.eye(2, dtype=int)\n",
"feas = [s for s in itertools.product(range(3), repeat=2) if sum(s) <= 2]\n",
"pref = [(2, 0), (1, 0), (0, 2), (0, 0), (0, 1), (1, 1)] # most preferred maximizer first\n",
"\n",
"def i2_rule(q):\n",
" scores = {s: int(np.array(q) @ (A_I2 @ np.array(s))) for s in feas}\n",
" best = max(scores.values())\n",
" cand = [s for s in feas if scores[s] == best]\n",
" return next(s for s in pref if s in cand)\n",
"\n",
"start = (0, 0); seen = {start}; frontier = [start]; edges = {}\n",
"while frontier:\n",
" q = frontier.pop(); s = np.array(i2_rule(q)); succ = set()\n",
" for a in ([1, 0], [0, 1]):\n",
" qp = tuple(np.array(q) + np.array(a) - (A_I2 @ s)); succ.add(qp)\n",
" if qp not in seen: seen.add(qp); frontier.append(qp)\n",
" edges[q] = succ\n",
"\n",
"def freach(src):\n",
" acc = {src}; st = [src]\n",
" while st:\n",
" for v in edges.get(st.pop(), ()):\n",
" if v not in acc: acc.add(v); st.append(v)\n",
" return acc\n",
"\n",
"reach = {u: freach(u) for u in seen}\n",
"sccs = []; done = set()\n",
"for u in seen:\n",
" if u in done: continue\n",
" comp = {v for v in seen if v in reach[u] and u in reach[v]}\n",
" sccs.append(comp); done |= comp\n",
"bottom = sorted(len(c) for c in sccs if all(edges[u] <= c for u in c))\n",
"reenter = start in {v for x in edges[start] for v in reach[x]}\n",
"print(f\"(a) I2 trap: reachable from 0 = {len(seen)} (paper 9); \"\n",
" f\"0 re-enterable = {reenter} (paper False); absorbing SCC sizes = {bottom} (paper [7])\")\n",
"\n",
"# (b) steering to 0 under the canonical idle clause\n",
"def steer_to_zero(A, q0):\n",
" A = np.asarray(A, int); n, m = A.shape\n",
" a_max = max(int(A[:, k].sum()) for k in range(m))\n",
" q = np.array(q0, int)\n",
" bound = int(np.abs(q).sum()) + 2 * a_max * int(np.maximum(q, 0).sum())\n",
" steps = 0\n",
" while np.any(q != 0):\n",
" if steps > bound: return False, steps, bound\n",
" scores = A.T @ q\n",
" if scores.max() <= 0: # idle: drain a negative coord\n",
" q = q + np.eye(n, dtype=int)[int(np.argmin(q))]\n",
" else: # active: s = 2 e_kstar\n",
" r = q - 2 * A[:, int(np.argmax(scores))]\n",
" neg = np.where(r <= -1)[0]\n",
" q = r + np.eye(n, dtype=int)[int(neg[0]) if neg.size else 0]\n",
" steps += 1\n",
" return True, steps, bound\n",
"\n",
"suite = {\n",
" \"I2\": A_I2,\n",
" \"candy\": np.asarray(sm.HyperPaddle(rates=[1, 1, 1, 1, 1, 1, 1]).incidence), # incidence from package\n",
" \"rem:wall\": np.array([[0, 1, 1, 0], [0, 0, 1, 1], [1, 1, 1, 1]]),\n",
" \"Moyal {1,2,3}\": np.array([[1], [1], [1]]),\n",
" \"Moyal {1,2,3}+{3,4}\": np.array([[1, 0], [1, 0], [1, 1], [0, 1]]),\n",
"}\n",
"rng2 = np.random.default_rng(0); total = good = 0\n",
"for name, A in suite.items():\n",
" for _ in range(200):\n",
" reached, steps, bound = steer_to_zero(A, rng2.integers(-4, 6, size=A.shape[0]))\n",
" total += 1; good += int(reached and steps <= bound)\n",
"print(f\"(b) steering: {good}/{total} start states over the suite reach 0 within the bound\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.14.5"
}
},
"nbformat": 4,
"nbformat_minor": 5
}