{
"cells": [
{
"cell_type": "markdown",
"id": "f0435d52-4276-4026-bfeb-8a741b01fc31",
"metadata": {},
"source": [
"# Buena Muerte Social Club"
]
},
{
"cell_type": "markdown",
"id": "b0f77714-e9ad-4a32-8382-9527016422b1",
"metadata": {},
"source": [
":::note\n",
"This is the companion notebook to the article *Buena Muerte Social Club*, submitted to [FUN 2026](https://fun2026.limos.fr/). It is NOT intended to be self-contained. Please refer to the paper for context and explanations.\n",
":::"
]
},
{
"cell_type": "markdown",
"id": "922af91d-3fd2-4696-aa90-e49f3c1f0d64",
"metadata": {},
"source": [
"## Preliminaries"
]
},
{
"cell_type": "markdown",
"id": "99b65834-2028-469a-8681-6db920df9d70",
"metadata": {},
"source": [
"### Basic Blood Types"
]
},
{
"cell_type": "markdown",
"id": "42dec34d-9cf2-4aa1-a814-ccddeaa2c464",
"metadata": {},
"source": [
"First, we define the different blood types. In addition to the $2^3$ main types depending on the presence of $A$, $B$, and $Rhesus$ markers, we add a void type that will be use later to model non-ghoul humans.\n",
"\n",
"An important property is *can_feed*, that tells if a vampire of given type can feed on a ghoul of other type."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "3ccbc307-90c6-4a63-87d3-1984672fa51d",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:29.186602Z",
"iopub.status.busy": "2026-01-12T15:27:29.185106Z",
"iopub.status.idle": "2026-01-12T15:27:29.316289Z",
"shell.execute_reply": "2026-01-12T15:27:29.316289Z",
"shell.execute_reply.started": "2026-01-12T15:27:29.186602Z"
}
},
"outputs": [],
"source": [
"import numpy as np \n",
"\n",
"class BloodType:\n",
" A = [\"\", \"A\"]\n",
" B = [\"\", \"B\"]\n",
" R = [\"-\", \"+\"]\n",
"\n",
" def __init__(self, i, j, k):\n",
" self.markers = np.array([i, j, k])\n",
" self.x = False\n",
"\n",
" def __repr__(self):\n",
" if self.x:\n",
" return 'X'\n",
" i, j, k = self.markers\n",
" result = f\"{self.A[i]}{self.B[j]}{self.R[k]}\"\n",
" if len(result) == 1:\n",
" return f\"O{result}\"\n",
" else:\n",
" return result\n",
"\n",
" @property\n",
" def sanitized_name(self):\n",
" res = self.__repr__()\n",
" return res.replace('-','m').replace('+','p')\n",
"\n",
" def can_feed_on(self, other):\n",
" return np.all(other.markers <= self.markers)\n",
"\n",
" @classmethod\n",
" def void(cls):\n",
" res = cls(1, 1, 1)\n",
" res.x = True\n",
" return res"
]
},
{
"cell_type": "markdown",
"id": "f7ff914a-2295-4e1f-894a-0b7747d8f4d0",
"metadata": {},
"source": [
"Let us compute all possible groups. We'll store them in a list, a class, and a dict for easy access."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "f6a15fa9-0e85-4410-a9da-2070519d2cfe",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:29.317305Z",
"iopub.status.busy": "2026-01-12T15:27:29.317305Z",
"iopub.status.idle": "2026-01-12T15:27:29.328238Z",
"shell.execute_reply": "2026-01-12T15:27:29.328238Z",
"shell.execute_reply.started": "2026-01-12T15:27:29.317305Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"[O-, O+, B-, B+, A-, A+, AB-, AB+]"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"class Groups:\n",
" pass\n",
"\n",
"groups = [BloodType(i, j, k) for i in range(2) for j in range(2) for k in range(2)]\n",
"g = Groups\n",
"for t in groups:\n",
" setattr(g, t.sanitized_name, t)\n",
"groups_ = {str(g): g for g in groups}\n",
"\n",
"groups"
]
},
{
"cell_type": "markdown",
"id": "8c3d6104-7410-4463-9bff-466f7cf4d8bd",
"metadata": {},
"source": [
"### Vampire/Ghoul Types"
]
},
{
"cell_type": "markdown",
"id": "96ef629a-334e-4e29-b28c-f649ece2c40e",
"metadata": {},
"source": [
"We now define Vampire/Ghoul types, which are defined by two blood types. The following properties are defined :\n",
"- *autarkic*: the vampire can feed on their ghoul.\n",
"- *compatible_with*: each vampire of two pairs can feed on the other one ghoul.\n",
"- *swinger*: a non plug-and-drink pair compatible with another non plug-and-drink pair.\n",
"- *helpless*: a non plug-and-drink pair that is only compatible with plug-and-drink-pairs.\n",
"- *three_compatible*: each vampire of three pairs can feed on a distinct available ghoul, in a cyclic fashion."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "245429b6-8b16-41d8-9495-3275928ce083",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:29.329584Z",
"iopub.status.busy": "2026-01-12T15:27:29.329584Z",
"iopub.status.idle": "2026-01-12T15:27:29.337543Z",
"shell.execute_reply": "2026-01-12T15:27:29.336537Z",
"shell.execute_reply.started": "2026-01-12T15:27:29.329584Z"
}
},
"outputs": [],
"source": [
"class VGType:\n",
" def __init__(self, v, g):\n",
" self.v = v\n",
" self.g = g\n",
" self.compatible_types = []\n",
"\n",
" def __repr__(self):\n",
" return f\"{self.v}/{self.g}\"\n",
"\n",
" @property\n",
" def category(self):\n",
" if self.auto_compatible:\n",
" return \"autarkic\"\n",
" if all(c.auto_compatible for c in self.compatible_types):\n",
" return \"helpless\"\n",
" return \"swinger\"\n",
" \n",
" @property\n",
" def auto_compatible(self):\n",
" return self.v.can_feed_on(self.g)\n",
"\n",
" def compatible_with(self, other):\n",
" return self.v.can_feed_on(other.g) and other.v.can_feed_on(self.g)\n",
"\n",
" def three_cycle(self, other1, other2):\n",
" return self.v.can_feed_on(other1.g) and other1.v.can_feed_on(other2.g) and other2.v.can_feed_on(self.g)\n",
" \n",
" def three_compatible_with(self, other1, other2):\n",
" return self.three_cycle(other1, other2) or self.three_cycle(other2, other1)"
]
},
{
"cell_type": "markdown",
"id": "2d9aad78-0beb-492d-ab80-622214fd2aa8",
"metadata": {},
"source": [
"We now prepare all possible VG pairs and observe their repartition."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "2643225d-c14d-49da-b1c6-8e201312f1b9",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:29.338543Z",
"iopub.status.busy": "2026-01-12T15:27:29.338543Z",
"iopub.status.idle": "2026-01-12T15:27:29.369275Z",
"shell.execute_reply": "2026-01-12T15:27:29.369275Z",
"shell.execute_reply.started": "2026-01-12T15:27:29.338543Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Counter({'autarkic': 27, 'helpless': 19, 'swinger': 18})"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from collections import Counter\n",
"vgs = [VGType(v, g) for v in groups for g in groups]\n",
"for v1 in vgs:\n",
" for v2 in vgs:\n",
" if v1.compatible_with(v2):\n",
" v1.compatible_types.append(v2)\n",
"vgs_ = {str(vg): i for i, vg in enumerate(vgs)}\n",
"Counter(vg.category for vg in vgs)"
]
},
{
"cell_type": "markdown",
"id": "3651463f-d2d3-4aa1-8183-2cd4becbcce9",
"metadata": {},
"source": [
"### Human types"
]
},
{
"cell_type": "markdown",
"id": "2c4c2aec-db78-43fb-a94c-b2de6bfeadfd",
"metadata": {},
"source": [
"We call *human* a potential ghoul not affiliated yet to any vampire. Humans can be modeled as a VG pair where the vampire is a placeholder, denoted $X$, compatible with everyone, e.g. *X/O+*"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "079e1f8e-76c9-4ad3-b8d9-c2b9ce0bc6fe",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:29.370281Z",
"iopub.status.busy": "2026-01-12T15:27:29.370281Z",
"iopub.status.idle": "2026-01-12T15:27:29.375041Z",
"shell.execute_reply": "2026-01-12T15:27:29.375041Z",
"shell.execute_reply.started": "2026-01-12T15:27:29.370281Z"
}
},
"outputs": [],
"source": [
"class Human(VGType):\n",
" category = \"human\"\n",
" def __init__(self, g):\n",
" self.v = BloodType.void() # alway compatible\n",
" self.g = g"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "f7cbb332-3ff7-4156-8fcf-f41877e400b1",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:29.378417Z",
"iopub.status.busy": "2026-01-12T15:27:29.377417Z",
"iopub.status.idle": "2026-01-12T15:27:29.383434Z",
"shell.execute_reply": "2026-01-12T15:27:29.382430Z",
"shell.execute_reply.started": "2026-01-12T15:27:29.378417Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"[X/O-, X/O+, X/B-, X/B+, X/A-, X/A+, X/AB-, X/AB+]"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"humans = [Human(g) for g in groups]\n",
"humans"
]
},
{
"cell_type": "markdown",
"id": "eead2b8c-db17-40a7-aef9-3942ddc9696c",
"metadata": {},
"source": [
"### Blood type distribution"
]
},
{
"cell_type": "markdown",
"id": "0173162c-7e8f-44f3-926e-b07b359af486",
"metadata": {},
"source": [
"For numeric computation, we need statistics on the blood type probabilities. We leverage Wikipedia to access recent data."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "aae163b5-45ee-4cb0-a05d-b686dccee2f3",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:29.384446Z",
"iopub.status.busy": "2026-01-12T15:27:29.384446Z",
"iopub.status.idle": "2026-01-12T15:27:30.631895Z",
"shell.execute_reply": "2026-01-12T15:27:30.631895Z",
"shell.execute_reply.started": "2026-01-12T15:27:29.384446Z"
}
},
"outputs": [],
"source": [
"import requests\n",
"import re\n",
"from bs4 import BeautifulSoup as bs\n",
"from fake_useragent import UserAgent\n",
"ua = UserAgent()\n",
"soup = bs(requests.get(\"https://en.wikipedia.org/wiki/Blood_type_distribution_by_country\", headers={'User-Agent': ua.random}).content)\n",
"\n",
"table = soup('table')[1]('tr')\n",
"names = [s.text.strip().replace(\"−\", \"-\") for s in table[0]('th')[2:]]\n",
"blood_type_by_country = dict()\n",
"for line in table[1:]:\n",
" country = line.th.text.strip().split('[')[0]\n",
" stats = {n: float(s.text.strip()[:-1]) if s.text[0].isdigit() else 0.0 for n, s in zip(names, line('td')[1:])}\n",
" total = sum(stats.values())\n",
" blood_type_by_country[country] = {k: v/total for k, v in stats.items()} \n",
" blood_type_by_country[country]['X'] = 1.0"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "54ea290e-2722-427f-867b-7f7ff511b361",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:30.632900Z",
"iopub.status.busy": "2026-01-12T15:27:30.632900Z",
"iopub.status.idle": "2026-01-12T15:27:30.638315Z",
"shell.execute_reply": "2026-01-12T15:27:30.638315Z",
"shell.execute_reply.started": "2026-01-12T15:27:30.632900Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"{'O+': 0.374,\n",
" 'A+': 0.35700000000000004,\n",
" 'B+': 0.085,\n",
" 'AB+': 0.034,\n",
" 'O-': 0.066,\n",
" 'A-': 0.063,\n",
" 'B-': 0.015,\n",
" 'AB-': 0.006,\n",
" 'X': 1.0}"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"blood_type_by_country['United States']"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "56001340-90a6-4b9e-816d-64a25d322e96",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:30.639320Z",
"iopub.status.busy": "2026-01-12T15:27:30.639320Z",
"iopub.status.idle": "2026-01-12T15:27:30.644580Z",
"shell.execute_reply": "2026-01-12T15:27:30.644580Z",
"shell.execute_reply.started": "2026-01-12T15:27:30.639320Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"{'O+': 0.365,\n",
" 'A+': 0.382,\n",
" 'B+': 0.077,\n",
" 'AB+': 0.025,\n",
" 'O-': 0.065,\n",
" 'A-': 0.068,\n",
" 'B-': 0.013999999999999999,\n",
" 'AB-': 0.004,\n",
" 'X': 1.0}"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"blood_type_by_country['France']"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "65d30348-180a-4508-823f-f30970696912",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:30.645585Z",
"iopub.status.busy": "2026-01-12T15:27:30.645585Z",
"iopub.status.idle": "2026-01-12T15:27:30.649754Z",
"shell.execute_reply": "2026-01-12T15:27:30.649754Z",
"shell.execute_reply.started": "2026-01-12T15:27:30.645585Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"{'O+': 0.3878396121603878,\n",
" 'A+': 0.2757297242702757,\n",
" 'B+': 0.0818099181900818,\n",
" 'AB+': 0.0201999798000202,\n",
" 'O-': 0.1323098676901323,\n",
" 'A-': 0.0818099181900818,\n",
" 'B-': 0.0201999798000202,\n",
" 'AB-': 0.00010099989900010099,\n",
" 'X': 1.0}"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"blood_type_by_country['World']"
]
},
{
"cell_type": "markdown",
"id": "1670d114-d53f-416f-b8e3-aefec6fb4a68",
"metadata": {},
"source": [
"Given a list of VG types and a country, issue normalized arrival rates."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "b33926a1-cd56-4b27-b759-7c891c1f57d6",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:30.650760Z",
"iopub.status.busy": "2026-01-12T15:27:30.650760Z",
"iopub.status.idle": "2026-01-12T15:27:30.656265Z",
"shell.execute_reply": "2026-01-12T15:27:30.655759Z",
"shell.execute_reply.started": "2026-01-12T15:27:30.650760Z"
}
},
"outputs": [],
"source": [
"def arrival_rates(vgs, country=\"United States\"):\n",
" probs = blood_type_by_country[country]\n",
" res = [probs[str(vg.v)]*probs[str(vg.g)] for vg in vgs ]\n",
" total = sum(res)\n",
" return [r/total for r in res]"
]
},
{
"cell_type": "markdown",
"id": "069f6607-4dd2-482b-9f0e-170d9c5c48c9",
"metadata": {},
"source": [
"For the record, one can compute distribution per VG category:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "91c7059b-f913-4595-8218-5576831828fb",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:30.657268Z",
"iopub.status.busy": "2026-01-12T15:27:30.657268Z",
"iopub.status.idle": "2026-01-12T15:27:30.663331Z",
"shell.execute_reply": "2026-01-12T15:27:30.663331Z",
"shell.execute_reply.started": "2026-01-12T15:27:30.657268Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"{'autarkic': 56.6078, 'helpless': 28.1786, 'swinger': 15.2136}"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"from collections import defaultdict\n",
"shares = defaultdict(float)\n",
"rates = arrival_rates(vgs)\n",
"for i, vg in enumerate(vgs):\n",
" shares[vg.category] += float(100*rates[i])\n",
"dict(shares)"
]
},
{
"cell_type": "markdown",
"id": "0c81aba7-689a-4894-88a7-313f5184d5cf",
"metadata": {},
"source": [
"## Instability study"
]
},
{
"cell_type": "markdown",
"id": "773494dc-126b-44b2-8187-00d0ab512e5a",
"metadata": {},
"source": [
"We know that the base model is not stabilisable. Yet, we still can investigate through simulations the practical performance of robust matching algorithm in that settings.\n",
"\n",
"Towards that goal, we write a function that converts a list of VG types into a matching problem."
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "1d0549fc-6fb6-4717-a697-4de6adfd52aa",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:30.665601Z",
"iopub.status.busy": "2026-01-12T15:27:30.664590Z",
"iopub.status.idle": "2026-01-12T15:27:32.332518Z",
"shell.execute_reply": "2026-01-12T15:27:32.332011Z",
"shell.execute_reply.started": "2026-01-12T15:27:30.665601Z"
}
},
"outputs": [],
"source": [
"import stochastic_matching as sm\n",
"\n",
"def build_model(vgs, country=\"United States\", self=True, relaxed=False, threesomes=False):\n",
" edges = []\n",
" n = len(vgs)\n",
" for i, vg in enumerate(vgs):\n",
" if relaxed or (self and vg.auto_compatible):\n",
" edges.append([i])\n",
" for j in range(i+1, n):\n",
" if vg.compatible_with(vgs[j]):\n",
" edges.append([i, j])\n",
" if not threesomes:\n",
" continue\n",
" for k in range(j+1, n):\n",
" if vg.three_compatible_with(vgs[j], vgs[k]):\n",
" if any(str(p.v)==str(p.g) for p in [vg, vgs[j], vgs[k]]):\n",
" continue\n",
" edges.append([i, j, k])\n",
" incidence = np.zeros((n, len(edges)), dtype=int)\n",
" for i, e in enumerate(edges):\n",
" incidence[e, i] = 1\n",
" model = sm.Model(incidence=incidence, names=[str(vg) for vg in vgs], \n",
" rates=arrival_rates(vgs, country=country))\n",
" model.cats = [vg.category for vg in vgs]\n",
" model.edges = edges\n",
" try:\n",
" model.adjacency = sm.model.incidence_to_adjacency(model.incidence)\n",
" except ValueError:\n",
" pass\n",
" return model"
]
},
{
"cell_type": "markdown",
"id": "34e0619c-661f-4f5b-83dd-77cfe1a5bf34",
"metadata": {},
"source": [
"### Swingers-only"
]
},
{
"cell_type": "markdown",
"id": "9c58b2df-3662-4e0d-91aa-832a914cb920",
"metadata": {},
"source": [
"We first focus on the inter-dependent types, because it is smaller and makes sense from a selfish perspective."
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "ec596efc-7003-4171-93a9-dc195556e9a6",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:32.333523Z",
"iopub.status.busy": "2026-01-12T15:27:32.333523Z",
"iopub.status.idle": "2026-01-12T15:27:32.345541Z",
"shell.execute_reply": "2026-01-12T15:27:32.345541Z",
"shell.execute_reply.started": "2026-01-12T15:27:32.333523Z"
}
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"swingers = [vg for vg in vgs if vg.category == 'swinger']\n",
"model = build_model(swingers)\n",
"model.show_graph()"
]
},
{
"cell_type": "markdown",
"id": "55c992ec-c1b6-4926-b63c-f2cc57e4afb9",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-09T16:09:23.999121Z",
"iopub.status.busy": "2026-01-09T16:09:23.999121Z",
"iopub.status.idle": "2026-01-09T16:09:24.003741Z",
"shell.execute_reply": "2026-01-09T16:09:24.003741Z",
"shell.execute_reply.started": "2026-01-09T16:09:23.999121Z"
}
},
"source": [
"Nice structure."
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "a495b1d8-d045-41e4-8185-c60d89494edd",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:32.347008Z",
"iopub.status.busy": "2026-01-12T15:27:32.347008Z",
"iopub.status.idle": "2026-01-12T15:27:32.353226Z",
"shell.execute_reply": "2026-01-12T15:27:32.353226Z",
"shell.execute_reply.started": "2026-01-12T15:27:32.347008Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Kernels of a graph with 18 nodes and 24 edges.\n",
"Node dimension is 0.\n",
"Edge dimension is 6\n",
"Type: Surjective-only\n",
"Node kernel:\n",
"[]\n",
"Edge kernel:\n",
"[[ 0.17028994 -0.17028994 -0.07412584 0.07412584 0. -0.17028994\n",
" -0.12829445 0.12829445 0. 0.07412584 0.12829445 0.\n",
" -0.32576712 0.19747267 -0.44229877 0.36817293 0.19747267 -0.19747267\n",
" 0.36817293 -0.36817293 0.12771746 0.04257249 0.04257249 -0.04257249]\n",
" [ 0.04439801 -0.04439801 -0.26529933 0.26529933 0. -0.04439801\n",
" -0.33257364 0.33257364 0. 0.26529933 0.33257364 0.\n",
" -0.38079546 0.04822182 0.0484483 -0.31374763 0.04822182 -0.04822182\n",
" -0.31374763 0.31374763 0.03329851 0.0110995 0.0110995 -0.0110995 ]\n",
" [ 0.35072029 -0.35072029 0.12045588 -0.12045588 0. -0.35072029\n",
" -0.113703 0.113703 0. -0.12045588 0.113703 0.\n",
" 0.12788313 -0.24158613 0.12232013 -0.00186425 -0.24158613 0.24158613\n",
" -0.00186425 0.00186425 0.51304022 -0.16231993 -0.16231993 0.16231993]\n",
" [ 0.02262008 -0.02262008 -0.20995298 0.20995298 0. -0.02262008\n",
" -0.08260737 0.08260737 0. 0.20995298 0.08260737 0.\n",
" 0.28638173 -0.3689891 -0.29147499 0.08152201 -0.3689891 0.3689891\n",
" 0.08152201 -0.08152201 -0.23303494 0.25565502 0.25565502 -0.25565502]\n",
" [ 0.33486514 -0.33486514 -0.03178975 0.03178975 0. -0.33486514\n",
" 0.22875339 -0.22875339 0. 0.03178975 -0.22875339 0.\n",
" 0.05848433 0.17026906 0.11609151 -0.14788126 0.17026906 -0.17026906\n",
" -0.14788126 0.14788126 0.00114885 0.33371628 0.33371628 -0.33371628]\n",
" [ 0.00676492 -0.00676492 -0.36219861 0.36219861 0. -0.00676492\n",
" 0.25984901 -0.25984901 0. 0.36219861 -0.25984901 0.\n",
" 0.21698293 0.04286608 -0.29770361 -0.064495 0.04286608 -0.04286608\n",
" -0.064495 0.064495 0.25507369 -0.24830877 -0.24830877 0.24830877]]"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model.kernel"
]
},
{
"cell_type": "markdown",
"id": "9e3b6e88-f0a7-40ed-81da-7d385bb1bff2",
"metadata": {},
"source": [
"As predicted, it is not stabilizable."
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "a1677bbb-b2a8-441a-b664-d620b7da1101",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:32.354564Z",
"iopub.status.busy": "2026-01-12T15:27:32.354564Z",
"iopub.status.idle": "2026-01-12T15:27:32.363837Z",
"shell.execute_reply": "2026-01-12T15:27:32.363837Z",
"shell.execute_reply.started": "2026-01-12T15:27:32.354564Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"False"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model.stabilizable"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "1563ffb7-1b71-421b-a7ab-ef8f51f8877d",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:32.364844Z",
"iopub.status.busy": "2026-01-12T15:27:32.364844Z",
"iopub.status.idle": "2026-01-12T15:27:32.368364Z",
"shell.execute_reply": "2026-01-12T15:27:32.368364Z",
"shell.execute_reply.started": "2026-01-12T15:27:32.364844Z"
}
},
"outputs": [],
"source": [
"def grouped_delay(simu):\n",
" from collections import defaultdict\n",
" import numpy as np\n",
" wait = defaultdict(float)\n",
" for i, cat in enumerate(simu.model.cats):\n",
" wait[cat] += float(simu.avg_queues[i])\n",
" return dict(wait)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "b90e8149-e15a-4d4a-96ab-4a5891b9caa1",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:32.369378Z",
"iopub.status.busy": "2026-01-12T15:27:32.369378Z",
"iopub.status.idle": "2026-01-12T15:27:32.374886Z",
"shell.execute_reply": "2026-01-12T15:27:32.374380Z",
"shell.execute_reply.started": "2026-01-12T15:27:32.369378Z"
}
},
"outputs": [],
"source": [
"base_params = {\"simulator\": \"longest\", \"n_steps\": 10**9, \"max_queue\": 100000, \"seed\": 42}"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "7741fa26-425b-408d-9c09-52755a155b7d",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:27:32.375888Z",
"iopub.status.busy": "2026-01-12T15:27:32.374886Z",
"iopub.status.idle": "2026-01-12T15:30:05.070814Z",
"shell.execute_reply": "2026-01-12T15:30:05.070814Z",
"shell.execute_reply.started": "2026-01-12T15:27:32.375888Z"
}
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "84792cdc57734e6d91fb91dcdec1cbd6",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/1 [00:00\n",
"\n",
" Reload\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"helpfuls = [vg for vg in vgs if vg.category != 'helpless']\n",
"model = build_model(helpfuls, self=False)\n",
"\n",
"colors = {\"autarkic\": \"#00FF00\", \"swinger\": \"#FF8888\", \"helpless\": \"#FF0000\"}\n",
"physics = {\n",
" \"solver\": \"barnesHut\",\n",
" \"barnesHut\": {\n",
" \"gravitationalConstant\": -15000, \n",
" \"centralGravity\": 0.25, \n",
" \"springLength\": 200, \n",
" \"springConstant\": 0.0015, \n",
" \"damping\": 0.095, \n",
" \"avoidOverlap\": 0.12 \n",
" },\n",
" \"maxVelocity\": 15,\n",
" \"stabilization\": {\n",
" \"iterations\": 800,\n",
" \"updateInterval\": 25\n",
" }\n",
" }\n",
"nodes_info = [{\"color\": {\"background\": colors[vg.category]}} for vg in helpfuls]\n",
"options = {\"width\": \"1000px\", \"nodes\": {\"font\": {\"size\": 60}}, \"physics\": physics}\n",
"model.show_graph(vis_options=options, nodes_info=nodes_info)"
]
},
{
"cell_type": "markdown",
"id": "64f2baaa-acd0-4f03-b5e7-90d6f9dba8f4",
"metadata": {},
"source": [
"Back to the real thing."
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "48f362c4-28d0-44ff-9256-e3d05b9daf40",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:30:05.158443Z",
"iopub.status.busy": "2026-01-12T15:30:05.158443Z",
"iopub.status.idle": "2026-01-12T15:30:05.255549Z",
"shell.execute_reply": "2026-01-12T15:30:05.255549Z",
"shell.execute_reply.started": "2026-01-12T15:30:05.158443Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Kernels of a graph with 45 nodes and 317 edges.\n",
"Node dimension is 0.\n",
"Edge dimension is 272\n",
"Type: Surjective-only\n",
"Node kernel:\n",
"[]\n",
"Edge kernel:\n",
"[[ 1.85577352e-01 4.60120460e-02 7.29821545e-02 ... 4.78644736e-04\n",
" -1.19916425e-03 -1.67780899e-03]\n",
" [ 1.70012321e-01 2.91401361e-02 6.69110172e-02 ... 2.85079976e-03\n",
" 9.21133624e-04 -1.92966614e-03]\n",
" [ 1.51967077e-01 2.44349889e-02 7.03749366e-02 ... -1.77408140e-03\n",
" -2.70700952e-04 1.50338044e-03]\n",
" ...\n",
" [-2.89699459e-02 -1.53198499e-02 -1.07370205e-02 ... 9.29353881e-01\n",
" -6.42094593e-02 6.43665994e-03]\n",
" [ 8.74057828e-03 8.55991905e-04 3.83245243e-04 ... -6.30747725e-02\n",
" 8.15041946e-01 -1.21883281e-01]\n",
" [ 3.77105242e-02 1.61758418e-02 1.11202657e-02 ... 7.57134673e-03\n",
" -1.20748594e-01 8.71680059e-01]]"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model = build_model(helpfuls)\n",
"model.kernel"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "0a96a988-7f28-4f3e-9819-4ca755dbe321",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:30:05.256774Z",
"iopub.status.busy": "2026-01-12T15:30:05.256774Z",
"iopub.status.idle": "2026-01-12T15:30:05.425368Z",
"shell.execute_reply": "2026-01-12T15:30:05.425368Z",
"shell.execute_reply.started": "2026-01-12T15:30:05.256774Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"True"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model.stabilizable"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "7e59c499-9cae-4f88-979d-5a4def7bb4ff",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:30:05.426730Z",
"iopub.status.busy": "2026-01-12T15:30:05.426730Z",
"iopub.status.idle": "2026-01-12T15:43:47.726564Z",
"shell.execute_reply": "2026-01-12T15:43:47.726051Z",
"shell.execute_reply.started": "2026-01-12T15:30:05.426730Z"
}
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "469d3f258d0d4987bbde88f30c803201",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/1 [00:00\n",
"\n",
" Reload\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"model = build_model(vgs, self=False)\n",
"nodes_info = [{\"color\": {\"background\": colors[vg.category]}} for vg in vgs]\n",
"options = {\"width\": \"1000px\", \"nodes\": {\"font\": {\"size\": 60}}, \"physics\": physics}\n",
"model.show_graph(vis_options=options, nodes_info=nodes_info)"
]
},
{
"cell_type": "markdown",
"id": "9a8ed4dc-d399-49d0-88ec-031c169e8765",
"metadata": {},
"source": [
"Back to the real thing."
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "a675fc45-e3fd-4745-ab29-3ca112f46c68",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:43:47.804906Z",
"iopub.status.busy": "2026-01-12T15:43:47.803906Z",
"iopub.status.idle": "2026-01-12T15:43:47.914094Z",
"shell.execute_reply": "2026-01-12T15:43:47.914094Z",
"shell.execute_reply.started": "2026-01-12T15:43:47.804906Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Kernels of a graph with 64 nodes and 378 edges.\n",
"Node dimension is 0.\n",
"Edge dimension is 314\n",
"Type: Surjective-only\n",
"Node kernel:\n",
"[]\n",
"Edge kernel:\n",
"[[-8.07660786e-02 -5.14815963e-03 1.55185186e-01 ... 1.23681204e-03\n",
" 1.07651019e-03 -1.60301853e-04]\n",
" [-5.27865346e-02 1.35874032e-02 1.66497191e-01 ... 4.31476721e-03\n",
" 3.36397446e-03 -9.50792748e-04]\n",
" [-7.31481008e-02 -1.48678109e-02 1.68269311e-01 ... -8.36046328e-04\n",
" 1.31250595e-03 2.14855228e-03]\n",
" ...\n",
" [ 1.05249784e-02 1.28211539e-02 -4.16977264e-02 ... 9.33214398e-01\n",
" -5.95732615e-02 7.21234047e-03]\n",
" [ 3.51189245e-02 1.37009886e-02 -1.00579058e-01 ... -5.94140329e-02\n",
" 8.18611345e-01 -1.21974622e-01]\n",
" [ 2.45939461e-02 8.79834759e-04 -5.88813319e-02 ... 7.37156906e-03\n",
" -1.21815393e-01 8.70813038e-01]]"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model = build_model(vgs)\n",
"model.kernel"
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "18138ff8-6511-4eba-9f9b-51bcfad3f201",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:43:47.915316Z",
"iopub.status.busy": "2026-01-12T15:43:47.915316Z",
"iopub.status.idle": "2026-01-12T15:43:48.247537Z",
"shell.execute_reply": "2026-01-12T15:43:48.247537Z",
"shell.execute_reply.started": "2026-01-12T15:43:47.915316Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"False"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model.stabilizable"
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "fb457646-8e41-4295-9b48-dd0e037548ce",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T15:43:48.249548Z",
"iopub.status.busy": "2026-01-12T15:43:48.248547Z",
"iopub.status.idle": "2026-01-12T15:56:06.505513Z",
"shell.execute_reply": "2026-01-12T15:56:06.505513Z",
"shell.execute_reply.started": "2026-01-12T15:43:48.249548Z"
}
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "50b1aae3860d43809c3eb134763222e4",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/1 [00:00"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from matplotlib import pyplot as plt\n",
"\n",
"r = res['discarding']\n",
"x = r['regret']\n",
"ys = r['grouped_delay']\n",
"for l in ys[0].keys():\n",
" y = [yy[l] for yy in ys]\n",
" plt.loglog(x, y, label=l)\n",
"plt.xlabel('Regret (wasted ghouls indicator)')\n",
"plt.ylabel('Delay (average queue size)')\n",
"plt.xlim([.00001, 1])\n",
"plt.legend()\n",
"import tikzplotlib\n",
"tikzplotlib.save(\"waste.tex\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "2ddd5185-3d86-4141-9912-bfc96fe53ab1",
"metadata": {},
"source": [
"### Humans"
]
},
{
"cell_type": "markdown",
"id": "58ff9776-3bdf-4788-a505-ace63e4646ed",
"metadata": {},
"source": [
"We now incoporate some amount of *humans* (vampire-free pre-ghouls) into the mix."
]
},
{
"cell_type": "code",
"execution_count": 36,
"id": "dceb9a99-dd0c-49f8-b4a1-5966c7e80171",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T19:10:51.447415Z",
"iopub.status.busy": "2026-01-12T19:10:51.447415Z",
"iopub.status.idle": "2026-01-12T19:10:51.452254Z",
"shell.execute_reply": "2026-01-12T19:10:51.452254Z",
"shell.execute_reply.started": "2026-01-12T19:10:51.447415Z"
}
},
"outputs": [],
"source": [
"def mix_model(α = .001):\n",
" μc = arrival_rates(vgs)\n",
" μd = arrival_rates(humans)\n",
" μ = [(1-α)*p/np.sum(μc) for p in μc]+[α*p/np.sum(μd) for p in μd]\n",
" model = build_model(vgs+humans)\n",
" model.rates = μ\n",
" return model"
]
},
{
"cell_type": "markdown",
"id": "a5f5d14a-7425-4b5e-aeb2-20d1e8bdd47b",
"metadata": {},
"source": [
":::note\n",
"If we try to numerically compute the minimum proportion of humans that allows the system to be stable, one gets a positive value. This is due to the numerical margin of error of the best positive position. For highly skewed problem, it can be difficult to distinguish a zeroed edge from a barely positive one.\n",
"\n",
"Yet, the value found has its one interest: it gives us a bound below which stability may be difficult to observe.\n",
":::"
]
},
{
"cell_type": "code",
"execution_count": 37,
"id": "88973aca-7541-407f-9cff-ea36409932c3",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T19:10:51.453260Z",
"iopub.status.busy": "2026-01-12T19:10:51.453260Z",
"iopub.status.idle": "2026-01-12T19:11:15.050401Z",
"shell.execute_reply": "2026-01-12T19:11:15.050401Z",
"shell.execute_reply.started": "2026-01-12T19:10:51.453260Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"0.0004916191101074219"
]
},
"execution_count": 37,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"def find_threshold(steps=20):\n",
" u = 1.0\n",
" l = 0.0\n",
" for _ in range(steps):\n",
" c = (u+l)/2\n",
" m = mix_model(c)\n",
" if m.stabilizable:\n",
" u = c\n",
" else:\n",
" l = c\n",
" return (u+l)/2\n",
"find_threshold()"
]
},
{
"cell_type": "code",
"execution_count": 38,
"id": "4aea5564-3c62-43c9-a7f6-e616e2d3afa3",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T19:11:15.051415Z",
"iopub.status.busy": "2026-01-12T19:11:15.051415Z",
"iopub.status.idle": "2026-01-12T19:11:15.055633Z",
"shell.execute_reply": "2026-01-12T19:11:15.055633Z",
"shell.execute_reply.started": "2026-01-12T19:11:15.051415Z"
}
},
"outputs": [],
"source": [
"base_params['simulator'] = 'virtual_queue'\n",
"xp = sm.XP('humanities', iterator=sm.Iterator('model', [2**i*.0005 for i in range(11)], 'alpha', lambda a: mix_model(a)), **base_params)"
]
},
{
"cell_type": "code",
"execution_count": 39,
"id": "a1453535-43c5-41c4-8986-ff84d8a16105",
"metadata": {
"execution": {
"iopub.execute_input": "2026-01-12T19:11:15.056662Z",
"iopub.status.busy": "2026-01-12T19:11:15.056662Z",
"iopub.status.idle": "2026-01-12T19:37:16.301614Z",
"shell.execute_reply": "2026-01-12T19:37:16.301614Z",
"shell.execute_reply.started": "2026-01-12T19:11:15.056662Z"
}
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "9d086249e4c24db59bcc26f55b6f09ef",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/11 [00:00"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"r = res['humanities']\n",
"x = r['alpha']\n",
"ys = r['grouped_delay']\n",
"for l in ys[0].keys():\n",
" y = [e[l] for e in ys]\n",
" plt.loglog(x, y, label=l)\n",
"plt.xlabel('Proportion $\\\\alpha$ of participating human volunteers')\n",
"plt.ylabel('Delay (average queue size)')\n",
"plt.xlim([.0005, .5])\n",
"plt.legend()\n",
"import tikzplotlib\n",
"tikzplotlib.save(\"humans.tex\")\n",
"plt.show()"
]
}
],
"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.12.8"
}
},
"nbformat": 4,
"nbformat_minor": 5
}