diff --git a/ipynb/Paint.ipynb b/ipynb/Paint.ipynb new file mode 100644 index 0000000..fd56852 --- /dev/null +++ b/ipynb/Paint.ipynb @@ -0,0 +1,965 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "d3bc52a7-6c7d-49b5-b6d4-b43a74e0ef34", + "metadata": {}, + "source": [ + "
Peter Norvig, April 2024
\n", + "\n", + "# Counting Cluster Sizes\n", + "\n", + "Zach Wissner-Gross's *Fiddler on the Proof!* column [poses a question](https://thefiddler.substack.com/p/can-you-paint-by-number) that I will restate as:\n", + "\n", + "Consider a two-dimensional grid of squares, where each square is colored either red or blue. A **cluster** is a group of contiguous squares of the same color (contiguity must be horizontal or vertical, not diagonal). For example, in the following 10 × 2 grid there are four red clusters and five blue clusters:\n", + "\n", + "![](https://substackcdn.com/image/fetch/w_1456,c_limit,f_webp,q_auto:good,fl_progressive:steep/https%3A%2F%2Fsubstack-post-media.s3.amazonaws.com%2Fpublic%2Fimages%2F079ab505-66c7-427b-ad6a-a2cf6a1794a6_1600x384.png)\n", + "\n", + "Altogether there are 20 squares and 9 clusters, so the average cluster size is about 2.22 squares. Under the assumption that every configuration of red and blue squares is equally likely, Zach poses two questions (and I'll add two more):\n", + "\n", + "- What is the average cluster size for a grid consisting of a single infinitely long row?\n", + "- What is the average cluster size for a grid consisting of two infinitely long rows?\n", + "- What is the average cluster size for a grid of any given size *w* × *h*?\n", + "- What if there are three or more colors rather than just two?\n", + "\n", + "# Code to Make Grids and to Count Clusters\n", + "\n", + "My strategy will be to try consider every possible grid of a given size, compute the average cluster size for each grid, and then average across all of them. I'll define two data types and five functions:\n", + "- `Square`: a pair of `(x, y)` coordinates, e.g. `(2, 1)`.\n", + "- `Grid`: data type for a grid; a dict of `{square: contents}`, where `contents` initially starts as a color (e.g. `'R'` or `'B'` for red or blue).\n", + "- `grids(width, height)`: returns a list of all possible grids of the given size. Each grid has a different combination of colors. With *c* colors, there are *c*(*w* × *h*) possible grids of size *w* × *h*.\n", + "- `one_grid(width, height, colorseq)`: returns one grid of given size, filling in squares with entries from colorseq.\n", + "- `mean_cluster_size(grid)`: returns the average size of the clusters in a grid.\n", + "- `cluster(grid)`: Mutates the grid so that each square's contents is changed from a color (a string) to a cluster number, where 1 is the number of the first cluster, 2 of the next cluster, and so on. Uses a [flood fill](https://en.wikipedia.org/wiki/Flood_fill) algorithm.\n", + "- `neighbors(square)`: the four squares surrounding the given square.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "d3991ae8-4c2c-4e4f-8884-95fd74aec5be", + "metadata": {}, + "outputs": [], + "source": [ + "import itertools \n", + "from statistics import mean, stdev\n", + "from typing import *\n", + "\n", + "Square = Tuple[int, int]\n", + "Grid = Dict[Square, Union[str, int]] # e.g. {(1, 1): 1, (2, 1): 'B'}\n", + "\n", + "def grids(width: int, height: int, colors='RB') -> List[Grid]:\n", + " \"\"\"All possible grids with given width, height, and color set.\"\"\"\n", + " return [one_grid(width, height, colorseq) \n", + " for colorseq in itertools.product(colors, repeat=(width * height))]\n", + "\n", + "def one_grid(width: int, height: int, colorseq: Sequence[str]) -> Grid: \n", + " \"\"\"A grid of given size made from the sequence of colors.\"\"\"\n", + " squares = [(x, y) for y in range(height) for x in range(width)]\n", + " return dict(zip(squares, colorseq))\n", + "\n", + "def mean_cluster_size(grid: Grid) -> float: \n", + " \"\"\"Mean size of clusters in a grid.\"\"\"\n", + " return len(grid) / max(cluster(grid).values())\n", + " \n", + "def cluster(grid: Grid) -> Grid:\n", + " \"\"\"Do a flood fill, replacing one cluster of adjacent characters with an integer,\n", + " then incrementing the integer for the next cluster and continuing.\"\"\"\n", + " cluster_number = 0 \n", + " for square in grid:\n", + " c = grid[square]\n", + " if isinstance(c, str):\n", + " cluster_number += 1\n", + " # Assign cluster number to square and all its neighbors with the same contents\n", + " Q = [square] # queue of squares in cluster `i`\n", + " while Q: \n", + " sq = Q.pop()\n", + " if grid.get(sq) == c: \n", + " grid[sq] = cluster_number\n", + " Q.extend(neighbors(sq))\n", + " return grid\n", + "\n", + "def neighbors(square: Square) -> List[Square]:\n", + " \"\"\"The four neighbors of a square.\"\"\"\n", + " (x, y) = square\n", + " return [(x + 1, y), (x - 1, y), (x, y + 1), (x, y - 1)]" + ] + }, + { + "cell_type": "markdown", + "id": "1bcebc7b-d59e-4bea-b735-c372fd91cfb1", + "metadata": {}, + "source": [ + "# Answering the Questions\n", + "\n", + "Here's one more function to help solve the questions:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "816f854a-7c99-4022-b539-98f87bc20369", + "metadata": {}, + "outputs": [], + "source": [ + "def do(n, h=1, colors='RB') -> dict: \n", + " \"\"\"A dict of {(w, h): mean-cluster-size} for grids with width w from 1 to n, and with height h.\"\"\"\n", + " return {(w, h): mean(map(mean_cluster_size, grids(w, h, colors)))\n", + " for w in range(1, n + 1)}" + ] + }, + { + "cell_type": "markdown", + "id": "0c5a9d5a-58b2-4cf9-b2ff-33078461bb97", + "metadata": {}, + "source": [ + "# One-Row Grids\n", + "\n", + "I can't make an infinitely long row (a grid with *n* squares has 2*n* possible color arrangments, it will be slow to investigate grids with more than about 20 squares (a million arrangements)). However, I can look at successively wider rows and see if the average cluster size seems to be converging:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "ee733740-c0b2-490b-9541-1b07e37a6056", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 1.29 s, sys: 43 ms, total: 1.33 s\n", + "Wall time: 1.34 s\n" + ] + }, + { + "data": { + "text/plain": [ + "{(1, 1): 1.0,\n", + " (2, 1): 1.5,\n", + " (3, 1): 1.75,\n", + " (4, 1): 1.875,\n", + " (5, 1): 1.9375,\n", + " (6, 1): 1.96875,\n", + " (7, 1): 1.984375,\n", + " (8, 1): 1.9921875,\n", + " (9, 1): 1.99609375,\n", + " (10, 1): 1.998046875,\n", + " (11, 1): 1.9990234375,\n", + " (12, 1): 1.99951171875,\n", + " (13, 1): 1.999755859375,\n", + " (14, 1): 1.9998779296875,\n", + " (15, 1): 1.99993896484375}" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "%time do(15)" + ] + }, + { + "cell_type": "markdown", + "id": "13dd529d-519e-487e-9430-e782c8b23e0f", + "metadata": {}, + "source": [ + "For rows that are at least 11 squares wide we're getting very close (within 0.001) to an average cluster size of 2.\n", + "\n", + "Now that I see the result, I can come up with a justification: Every cluster starts with one square. The next square in the row will be the same color half the time. Continuing, we would get three of the same color in a row a quarter of the time, four in a row an eigth of the time, and so on. So we get:\n", + "\n", + "mean cluster size   =   Σi∈{0,1,...∞} (1/2)n   =   2\n", + "\n", + "# Two-Row Grids\n", + "\n", + "Now consider grids that are two rows tall:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "e2f0c3ca-acac-4dd8-9970-c846e3389ad5", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 8.12 s, sys: 269 ms, total: 8.39 s\n", + "Wall time: 8.71 s\n" + ] + }, + { + "data": { + "text/plain": [ + "{(1, 2): 1.5,\n", + " (2, 2): 2.125,\n", + " (3, 2): 2.46875,\n", + " (4, 2): 2.6682291666666664,\n", + " (5, 2): 2.7920386904761907,\n", + " (6, 2): 2.873721168154762,\n", + " (7, 2): 2.9304973563762626,\n", + " (8, 2): 2.971709969311288,\n", + " (9, 2): 3.0027141140414764}" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "%time do(9, 2)" + ] + }, + { + "cell_type": "markdown", + "id": "8e1ac7a1-f947-4b76-8c5d-d1b26974d018", + "metadata": {}, + "source": [ + "It looks like the answer is converging on 3, but it is not as clear as with the first question.\n", + "\n", + "Here's a justification: consider a cluster that extends a certain length all along the top row, with possibly some additional squares below it in the bootom row. You can call this a π-shape, because there is a horizontal bar at the top, and ero or more vertical bars extending down. By the argument above, the average length of the bar is 2. Now for each column in the bar there is a 50% chance that the square below will be the same color and thus join the cluster. So the average cluster size in total is 2 + 2 × 1/2 = 3. By a similar argument, clusters that extend along the bottom row with possibly some squares above also have an average size of 3. But what about an Z-shape or S-shape cluster, which starts along the top (or bottom) row, then switches to the bottom (or top)? The trick is that you can take any Z- or S-shape and turn it into a π-shape by just swapping top and bottom squares when the cluster is on the bottom. This swap doesn't change any probabilities, so a Z-, S-, and π-shaped clusters all have a mean cluster size of 3.\n", + "\n", + "# Arbitrary *w* × *h* Grids\n", + "\n", + "Let's consider grids of size *n* × 3 and *n* × 4:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "b29ddec0-5fd2-4f8e-96af-f4b380c6eb59", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 6.52 s, sys: 170 ms, total: 6.69 s\n", + "Wall time: 6.72 s\n" + ] + }, + { + "data": { + "text/plain": [ + "{(1, 3): 1.75,\n", + " (2, 3): 2.46875,\n", + " (3, 3): 2.8986049107142855,\n", + " (4, 3): 3.1591145833333334,\n", + " (5, 3): 3.326567058836346,\n", + " (6, 3): 3.4407984463961334}" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "%time do(6, 3)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "cb493f9f-734e-4767-b245-24cdf906a06a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 27.9 s, sys: 952 ms, total: 28.9 s\n", + "Wall time: 29.1 s\n" + ] + }, + { + "data": { + "text/plain": [ + "{(1, 4): 1.875,\n", + " (2, 4): 2.6682291666666664,\n", + " (3, 4): 3.1591145833333334,\n", + " (4, 4): 3.461628834083668,\n", + " (5, 4): 3.6608226758600444}" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "%time do(5, 4)" + ] + }, + { + "cell_type": "markdown", + "id": "26567fd7-dbf5-44b1-98a6-de8d4f95a20f", + "metadata": {}, + "source": [ + "# Adding Colors\n", + "\n", + "What if we add a third color? That should make the average cluster size smaller (since there is a 2/3 rather than 1/2 chance that the next square will be a different color and start a new cluster). I'll start with a single row:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "dc4dc1f5-f5f4-4d69-8bdf-c00d21d978cd", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 13.1 s, sys: 510 ms, total: 13.6 s\n", + "Wall time: 13.7 s\n" + ] + }, + { + "data": { + "text/plain": [ + "{(1, 1): 1.0,\n", + " (2, 1): 1.3333333333333333,\n", + " (3, 1): 1.4444444444444444,\n", + " (4, 1): 1.4814814814814814,\n", + " (5, 1): 1.4938271604938271,\n", + " (6, 1): 1.4979423868312758,\n", + " (7, 1): 1.4993141289437586,\n", + " (8, 1): 1.4997713763145861,\n", + " (9, 1): 1.499923792104862,\n", + " (10, 1): 1.4999745973682874,\n", + " (11, 1): 1.4999915324560957,\n", + " (12, 1): 1.4999971774853653}" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "%time do(12, 1, 'RGB')" + ] + }, + { + "cell_type": "markdown", + "id": "7afd16c2-99bd-410b-b5a8-d20862b7b1f5", + "metadata": {}, + "source": [ + "This is straightforward: the mean cluster size converges to 3/2. We can update the equation to deal with a single row with *c* equiprobable colors:\n", + "\n", + "mean cluster size   =   Σi∈{0,1,...∞} (1/*c*)n   =   *c* / (*c* - 1)\n", + "\n", + "Now for a two-row grid with three colors:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "37f4a825-d75c-4275-b2c4-6dd3de080b3f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 10.5 s, sys: 412 ms, total: 11 s\n", + "Wall time: 11.3 s\n" + ] + }, + { + "data": { + "text/plain": [ + "{(1, 2): 1.3333333333333333,\n", + " (2, 2): 1.654320987654321,\n", + " (3, 2): 1.7736625514403292,\n", + " (4, 2): 1.8261240664532845,\n", + " (5, 2): 1.8535078615096905,\n", + " (6, 2): 1.8697928838761029}" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "%time do(6, 2, 'RGB')" + ] + }, + { + "cell_type": "markdown", + "id": "3c18f723-447f-4cd6-ba27-59db37e87659", + "metadata": {}, + "source": [ + "By the same reasoning as we had last time, this should converge to 3/2 + 3/2 × 1/3 = 2. " + ] + }, + { + "cell_type": "markdown", + "id": "6a36b07b-8e5c-4048-8046-7619b34f6b26", + "metadata": {}, + "source": [ + "# Random Sampling for Larger Grids\n", + "\n", + "I can't analyze an infinite two-dimensional grid.\n", + "\n", + "And I can't even enumerate all the, say, 100 × 100 grids, because there are 210,000 of them.\n", + "\n", + "But I can **randomly sample** a bunch of 100 × 100 grids, and that should give me a good estimate of the full distribution of grids.\n", + "\n", + "First, a function to make a random grid of a given size:" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "51cc5dca-f3b4-4a2b-96af-479683216d4b", + "metadata": {}, + "outputs": [], + "source": [ + "import random\n", + "\n", + "def random_grid(width: int, height: int, colors='RB') -> Grid:\n", + " \"\"\"A width × height grid filled with random colors.\"\"\"\n", + " return one_grid(width, height, [random.choice(colors) for _ in range(width * height)])" + ] + }, + { + "cell_type": "markdown", + "id": "50ac4938-fd46-4e27-89d0-82bca289f37d", + "metadata": {}, + "source": [ + "Here's the mean cluster size of one random grid:" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "6abd1ae5-852f-459c-ad6e-1e91d91566b0", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "7.58150113722517" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mean_cluster_size(random_grid(100, 100))" + ] + }, + { + "cell_type": "markdown", + "id": "197a9c3e-58e4-4a4d-a37f-c793c98f36e8", + "metadata": {}, + "source": [ + "But that's just one grid. Better to sample *N* different grids (say, *N* = 200) and report on the mean cluster sizes for them:" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "39fd3321-907d-45eb-b06b-70d590427a09", + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "\n", + "def report(numbers, bins=30) -> dict:\n", + " \"\"\"Plot a histogram and return statistics on these numbers.\"\"\"\n", + " plt.hist(numbers, bins=bins, rwidth=0.8, align='left')\n", + " return dict(mean=mean(numbers), max=max(numbers), min=min(numbers), \n", + " stdev=stdev(numbers), N=len(numbers))" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "083dc3ea-a395-4eae-b751-347e4384a373", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'mean': 7.246114781372947,\n", + " 'max': 7.716049382716049,\n", + " 'min': 6.830601092896175,\n", + " 'stdev': 0.181119344027548,\n", + " 'N': 200}" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAdfUlEQVR4nO3dfZBV9X348c/KygVSHgoJD1vARYtBxQIVkygkwkSwiCRN6kMaBSoxgyMRldaw1PhT0pFFO7U0UnB0UkwHIU6NUg3xgUkCiHkSDDYxGRCEsEGRqTW7gOkq7Pn90XEnKw/u1XO/u3f39Zo5f9xzv2fPhzmuvOfce7kVWZZlAQCQyEltPQAA0LmIDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASKqyrQd4t6ampnjllVeiZ8+eUVFR0dbjAACtkGVZHDhwIKqqquKkk058b6Pdxccrr7wSQ4YMaesxAID3oa6uLgYPHnzCNe0uPnr27BkR/zd8r1692ngaAKA1GhoaYsiQIc1/j59Iu4uPd15q6dWrl/gAgDLTmrdMeMMpAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACCpyrYeAGhb1TVri1q/e/HUEk0CdBbufAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASKro+Ni4cWNMmzYtqqqqoqKiItasWXPUml//+tfxmc98Jnr37h09e/aMT3ziE7Fnz5485gUAylzR8XHo0KEYNWpULF269JjP79y5M8aPHx8jRoyI9evXxwsvvBC33nprdOvW7QMPCwCUv8piD5gyZUpMmTLluM/fcsstcfHFF8ddd93VvO/UU099f9MBAB1Oru/5aGpqirVr18bpp58eF110UfTv3z8+/vGPH/OlmXc0NjZGQ0NDiw0A6LhyjY/9+/fHwYMHY/HixfEXf/EX8fTTT8fnPve5+PznPx8bNmw45jG1tbXRu3fv5m3IkCF5jgQAtDO53/mIiPjsZz8bN910U4wePTpqamrikksuiXvvvfeYxyxYsCDq6+ubt7q6ujxHAgDamaLf83EiH/7wh6OysjLOPPPMFvvPOOOM2LRp0zGPKRQKUSgU8hwDAGjHcr3z0bVr1zj33HNj27ZtLfZv3749TjnllDxPBQCUqaLvfBw8eDB27NjR/HjXrl2xdevW6Nu3bwwdOjRuvvnmuOKKK+JTn/pUTJw4MZ588sl4/PHHY/369XnODQCUqaLjY/PmzTFx4sTmx/PmzYuIiJkzZ8YDDzwQn/vc5+Lee++N2tramDt3bnz0ox+N73znOzF+/Pj8pgYAylbR8TFhwoTIsuyEa2bNmhWzZs1630MBAB2X73YBAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJFV0fGzcuDGmTZsWVVVVUVFREWvWrDnu2tmzZ0dFRUUsWbLkA4wIAHQkRcfHoUOHYtSoUbF06dITrluzZk389Kc/jaqqqvc9HADQ8VQWe8CUKVNiypQpJ1yzd+/e+MpXvhJPPfVUTJ069X0PBwB0PEXHx3tpamqK6dOnx8033xxnnXXWe65vbGyMxsbG5scNDQ15jwQAtCO5x8edd94ZlZWVMXfu3Fatr62tjYULF+Y9BpSV6pq1RR+ze3Hb31Usdu72MDPQ9nL9tMuWLVviX/7lX+KBBx6IioqKVh2zYMGCqK+vb97q6uryHAkAaGdyjY9nnnkm9u/fH0OHDo3KysqorKyM3/zmN/G3f/u3UV1dfcxjCoVC9OrVq8UGAHRcub7sMn369Ljwwgtb7Lvoooti+vTpcfXVV+d5KgCgTBUdHwcPHowdO3Y0P961a1ds3bo1+vbtG0OHDo1+/fq1WH/yySfHwIED46Mf/egHnxYAKHtFx8fmzZtj4sSJzY/nzZsXEREzZ86MBx54ILfBAICOqej4mDBhQmRZ1ur1u3fvLvYUAEAH5rtdAICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJVbb1AAApVdesLWr97sVTSzQJdF7ufAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASKro+Ni4cWNMmzYtqqqqoqKiItasWdP83Ntvvx3z58+Ps88+Oz70oQ9FVVVVzJgxI1555ZU8ZwYAyljR8XHo0KEYNWpULF269Kjn3nzzzXj++efj1ltvjeeffz4eeeSR2L59e3zmM5/JZVgAoPxVFnvAlClTYsqUKcd8rnfv3rFu3boW++6555742Mc+Fnv27ImhQ4e+vykBgA6j6PgoVn19fVRUVESfPn2O+XxjY2M0NjY2P25oaCj1SABAGyppfPzv//5v1NTUxBe/+MXo1avXMdfU1tbGwoULSzkGtFp1zdqi1u9ePLVEk3AirhOUt5J92uXtt9+OL3zhC9HU1BTLli077roFCxZEfX1981ZXV1eqkQCAdqAkdz7efvvtuPzyy2PXrl3xgx/84Lh3PSIiCoVCFAqFUowBALRDucfHO+Hx0ksvxQ9/+MPo169f3qcAAMpY0fFx8ODB2LFjR/PjXbt2xdatW6Nv375RVVUVl156aTz//PPx3e9+N44cORL79u2LiIi+fftG165d85scAChLRcfH5s2bY+LEic2P582bFxERM2fOjNtvvz0ee+yxiIgYPXp0i+N++MMfxoQJE97/pABAh1B0fEyYMCGyLDvu8yd6DgDAd7sAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkio6PjZu3BjTpk2LqqqqqKioiDVr1rR4PsuyuP3226Oqqiq6d+8eEyZMiBdffDGveQGAMld0fBw6dChGjRoVS5cuPebzd911V9x9992xdOnSeO6552LgwIExadKkOHDgwAceFgAof5XFHjBlypSYMmXKMZ/LsiyWLFkSt9xyS3z+85+PiIhvfetbMWDAgFi1alXMnj37g00LAJS9XN/zsWvXrti3b19Mnjy5eV+hUIgLLrggfvSjHx3zmMbGxmhoaGixAQAdV9F3Pk5k3759ERExYMCAFvsHDBgQv/nNb455TG1tbSxcuDDPMaBNVNesLWr97sVTSzQJQPtWkk+7VFRUtHicZdlR+96xYMGCqK+vb97q6upKMRIA0E7keudj4MCBEfF/d0AGDRrUvH///v1H3Q15R6FQiEKhkOcYAEA7luudj2HDhsXAgQNj3bp1zfveeuut2LBhQ5x//vl5ngoAKFNF3/k4ePBg7Nixo/nxrl27YuvWrdG3b98YOnRo3HjjjbFo0aIYPnx4DB8+PBYtWhQ9evSIL37xi7kODgCUp6LjY/PmzTFx4sTmx/PmzYuIiJkzZ8YDDzwQX/3qV+P3v/99XHfddfHGG2/Exz/+8Xj66aejZ8+e+U0NAJStouNjwoQJkWXZcZ+vqKiI22+/PW6//fYPMhcA0EH5bhcAICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFKVbT0AHVd1zdqi1u9ePLVEkwDQnrjzAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEnlHh+HDx+Or33tazFs2LDo3r17nHrqqfH1r389mpqa8j4VAFCGKvP+gXfeeWfce++98a1vfSvOOuus2Lx5c1x99dXRu3fvuOGGG/I+HQBQZnKPjx//+Mfx2c9+NqZOnRoREdXV1bF69erYvHlz3qcCAMpQ7i+7jB8/Pr7//e/H9u3bIyLihRdeiE2bNsXFF1+c96kAgDKU+52P+fPnR319fYwYMSK6dOkSR44ciTvuuCP++q//+pjrGxsbo7GxsflxQ0ND3iMBAO1I7vHx0EMPxcqVK2PVqlVx1llnxdatW+PGG2+MqqqqmDlz5lHra2trY+HChXmPQU6qa9YWtX734qklmqT1ynHmD6oz/pnL0Qe5Tq4xHUnuL7vcfPPNUVNTE1/4whfi7LPPjunTp8dNN90UtbW1x1y/YMGCqK+vb97q6uryHgkAaEdyv/Px5ptvxkkntWyaLl26HPejtoVCIQqFQt5jAADtVO7xMW3atLjjjjti6NChcdZZZ8XPf/7zuPvuu2PWrFl5nwoAKEO5x8c999wTt956a1x33XWxf//+qKqqitmzZ8f/+3//L+9TAQBlKPf46NmzZyxZsiSWLFmS948GADoA3+0CACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEiqsq0HgGOprllb1Prdi6eWaBJKxTWmVPy31f658wEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS4gMASEp8AABJiQ8AIKmSxMfevXvjqquuin79+kWPHj1i9OjRsWXLllKcCgAoM5V5/8A33ngjxo0bFxMnTownnngi+vfvHzt37ow+ffrkfSoAoAzlHh933nlnDBkyJFasWNG8r7q6Ou/TAABlKveXXR577LEYO3ZsXHbZZdG/f/8YM2ZM3H///cdd39jYGA0NDS02AKDjyv3Ox8svvxzLly+PefPmxd///d/Hz372s5g7d24UCoWYMWPGUetra2tj4cKFeY/BH6iuWVvU+t2Lp5ZoEihvfpcgH7nf+Whqaoo///M/j0WLFsWYMWNi9uzZ8eUvfzmWL19+zPULFiyI+vr65q2uri7vkQCAdiT3+Bg0aFCceeaZLfadccYZsWfPnmOuLxQK0atXrxYbANBx5R4f48aNi23btrXYt3379jjllFPyPhUAUIZyj4+bbropfvKTn8SiRYtix44dsWrVqrjvvvtizpw5eZ8KAChDucfHueeeG48++misXr06Ro4cGf/wD/8QS5YsiSuvvDLvUwEAZSj3T7tERFxyySVxySWXlOJHAwBlzne7AABJiQ8AICnxAQAkJT4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSqmzrAQAoreqatUWt3714aokmab1ynJnWc+cDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkip5fNTW1kZFRUXceOONpT4VAFAGShofzz33XNx3333xZ3/2Z6U8DQBQRkoWHwcPHowrr7wy7r///vjjP/7jUp0GACgzJYuPOXPmxNSpU+PCCy884brGxsZoaGhosQEAHVdlKX7ot7/97Xj++efjueeee8+1tbW1sXDhwlKMAUAbqq5ZW9T63YunlmgS2pvc73zU1dXFDTfcECtXroxu3bq95/oFCxZEfX1981ZXV5f3SABAO5L7nY8tW7bE/v3745xzzmned+TIkdi4cWMsXbo0Ghsbo0uXLs3PFQqFKBQKeY8BALRTucfHpz/96fjFL37RYt/VV18dI0aMiPnz57cIDwCg88k9Pnr27BkjR45sse9DH/pQ9OvX76j9AEDn4184BQCSKsmnXd5t/fr1KU4DAJQBdz4AgKTEBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAElVtvUAALRf1TVri1q/e/HUEk1CR+LOBwCQlPgAAJISHwBAUuIDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApHKPj9ra2jj33HOjZ8+e0b9///jLv/zL2LZtW96nAQDKVO7xsWHDhpgzZ0785Cc/iXXr1sXhw4dj8uTJcejQobxPBQCUocq8f+CTTz7Z4vGKFSuif//+sWXLlvjUpz6V9+kAgDKTe3y8W319fURE9O3b95jPNzY2RmNjY/PjhoaGUo8EALShksZHlmUxb968GD9+fIwcOfKYa2pra2PhwoWlHKOF6pq1Ra3fvXhqiSZpvWJnjmgfcwOUm3L9/225/d1W0k+7fOUrX4n/+q//itWrVx93zYIFC6K+vr55q6urK+VIAEAbK9mdj+uvvz4ee+yx2LhxYwwePPi46wqFQhQKhVKNAQC0M7nHR5Zlcf3118ejjz4a69evj2HDhuV9CgCgjOUeH3PmzIlVq1bFf/7nf0bPnj1j3759ERHRu3fv6N69e96nAwDKTO7v+Vi+fHnU19fHhAkTYtCgQc3bQw89lPepAIAyVJKXXQAAjsd3uwAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkqps6wHKSXXN2qLW7148tUSTAED5cucDAEhKfAAASYkPACAp8QEAJCU+AICkxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkipZfCxbtiyGDRsW3bp1i3POOSeeeeaZUp0KACgjJYmPhx56KG688ca45ZZb4uc//3l88pOfjClTpsSePXtKcToAoIyUJD7uvvvu+NKXvhTXXHNNnHHGGbFkyZIYMmRILF++vBSnAwDKSGXeP/Ctt96KLVu2RE1NTYv9kydPjh/96EdHrW9sbIzGxsbmx/X19RER0dDQkPdoERHR1PhmUev/cI4PcuwHUex5333utvozO7Z0x7bluR37/o5ty3M7tnTHvvv4ttJWfz8d62dmWfbei7Oc7d27N4uI7Nlnn22x/4477shOP/30o9bfdtttWUTYbDabzWbrAFtdXd17tkLudz7eUVFR0eJxlmVH7YuIWLBgQcybN6/5cVNTU/zP//xP9OvX75jrj6WhoSGGDBkSdXV10atXrw82OLlybdon16X9cm3aL9fmxLIsiwMHDkRVVdV7rs09Pj784Q9Hly5dYt++fS3279+/PwYMGHDU+kKhEIVCocW+Pn36vK9z9+rVy38Q7ZRr0z65Lu2Xa9N+uTbH17t371aty/0Np127do1zzjkn1q1b12L/unXr4vzzz8/7dABAmSnJyy7z5s2L6dOnx9ixY+O8886L++67L/bs2RPXXnttKU4HAJSRksTHFVdcEa+//np8/etfj1dffTVGjhwZ3/ve9+KUU04pxemiUCjEbbfddtTLN7Q916Z9cl3aL9em/XJt8lORZa35TAwAQD58twsAkJT4AACSEh8AQFLiAwBIqiziY+/evXHVVVdFv379okePHjF69OjYsmXLCY958MEHY9SoUdGjR48YNGhQXH311fH6668nmrhzqK6ujoqKiqO2OXPmHPeYDRs2xDnnnBPdunWLU089Ne69996EE3cexV6bRx55JCZNmhQf+chHolevXnHeeefFU089lXjqju/9/M6849lnn43KysoYPXp06QfthN7PtWlsbIxbbrklTjnllCgUCnHaaafFv/3bvyWcunyV7J9Xz8sbb7wR48aNi4kTJ8YTTzwR/fv3j507d57wX0HdtGlTzJgxI/75n/85pk2bFnv37o1rr702rrnmmnj00UfTDd/BPffcc3HkyJHmx7/85S9j0qRJcdlllx1z/a5du+Liiy+OL3/5y7Fy5cp49tln47rrrouPfOQj8Vd/9Vepxu4Uir02GzdujEmTJsWiRYuiT58+sWLFipg2bVr89Kc/jTFjxqQau8Mr9rq8o76+PmbMmBGf/vSn47XXXiv1mJ3S+7k2l19+ebz22mvxzW9+M/70T/809u/fH4cPH04xbvnL5dvkSmj+/PnZ+PHjizrmH//xH7NTTz21xb5vfOMb2eDBg/McjXe54YYbstNOOy1ramo65vNf/epXsxEjRrTYN3v27OwTn/hEivE6tfe6Nsdy5plnZgsXLizhVLT2ulxxxRXZ1772tey2227LRo0alWa4Tu69rs0TTzyR9e7dO3v99dcTT9YxtPuXXR577LEYO3ZsXHbZZdG/f/8YM2ZM3H///Sc85vzzz4/f/va38b3vfS+yLIvXXnstHn744Zg6dWqiqTuft956K1auXBmzZs067hcC/vjHP47Jkye32HfRRRfF5s2b4+23304xZqfUmmvzbk1NTXHgwIHo27dviafrvFp7XVasWBE7d+6M2267LeF0nVtrrs07fzfddddd8Sd/8idx+umnx9/93d/F73//+8TTlqd2Hx8vv/xyLF++PIYPHx5PPfVUXHvttTF37tz493//9+Mec/7558eDDz4YV1xxRXTt2jUGDhwYffr0iXvuuSfh5J3LmjVr4ne/+138zd/8zXHX7Nu376gvFxwwYEAcPnw4/vu//7vEE3Zerbk27/ZP//RPcejQobj88stLN1gn15rr8tJLL0VNTU08+OCDUVnZ7l8l7zBac21efvnl2LRpU/zyl7+MRx99NJYsWRIPP/xwq96/Q7T/l11OPvnk7Lzzzmux7/rrrz/hrfoXX3wxGzRoUHbXXXdlL7zwQvbkk09mZ599djZr1qxSj9tpTZ48ObvkkktOuGb48OHZokWLWuzbtGlTFhHZq6++WsrxOrXWXJs/tGrVqqxHjx7ZunXrSjgV73VdDh8+nI0dOzZbvnx58z4vu6TRmt+ZSZMmZd26dct+97vfNe/7zne+k1VUVGRvvvlmqUcse+0+PoYOHZp96UtfarFv2bJlWVVV1XGPueqqq7JLL720xb5nnnkmi4jslVdeKcmcndnu3buzk046KVuzZs0J133yk5/M5s6d22LfI488klVWVmZvvfVWKUfstFp7bd7x7W9/O+vevXv23e9+t8STdW6tuS5vvPFGFhFZly5dmreKiormfd///vcTTtx5tPZ3ZsaMGdlpp53WYt+vfvWrLCKy7du3l3LEDqHd38cbN25cbNu2rcW+7du3n/BL6t58882jblF26dIlIiIyX2WTuxUrVkT//v3f8z015513Xjz++OMt9j399NMxduzYOPnkk0s5YqfV2msTEbF69eqYNWtWrF692vujSqw116VXr17xi1/8osW+ZcuWxQ9+8IN4+OGHY9iwYaUes1Nq7e/MuHHj4j/+4z/i4MGD8Ud/9EcR8X9/N5100kkxePDgFKOWt7aun/fys5/9LKusrMzuuOOO7KWXXsoefPDBrEePHtnKlSub19TU1GTTp09vfrxixYqssrIyW7ZsWbZz585s06ZN2dixY7OPfexjbfFH6NCOHDmSDR06NJs/f/5Rz737urz88stZjx49sptuuin71a9+lX3zm9/MTj755Ozhhx9OOXKnUcy1WbVqVVZZWZn967/+a/bqq682b394S5l8FHNd3s3LLqVVzLU5cOBANnjw4OzSSy/NXnzxxWzDhg3Z8OHDs2uuuSblyGWr3cdHlmXZ448/no0cOTIrFArZiBEjsvvuu6/F8zNnzswuuOCCFvu+8Y1vZGeeeWbWvXv3bNCgQdmVV16Z/fa3v004defw1FNPZRGRbdu27ajnjnVd1q9fn40ZMybr2rVrVl1d3eL1bPJVzLW54IILsog4aps5c2a6gTuJYn9n/pD4KK1ir82vf/3r7MILL8y6d++eDR48OJs3b573e7RSRZZ5HQIASKfdf9QWAOhYxAcAkJT4AACSEh8AQFLiAwBISnwAAEmJDwAgKfEBACQlPgCApMQHAJCU+AAAkhIfAEBS/x970SbjkheWjgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "N = 200\n", + "sizes = [mean_cluster_size(random_grid(100, 100)) for _ in range(N)]\n", + "report(sizes)" + ] + }, + { + "cell_type": "markdown", + "id": "9f7b9a04-31ff-44ed-bbf4-19b56c910859", + "metadata": {}, + "source": [ + "A 100 x 100 grid isn't infinite, so let's look at larger grids:" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "9e1f69b9-870f-4439-ad62-92d91684392d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'mean': 7.402852321356604,\n", + " 'max': 7.635044855888529,\n", + " 'min': 7.126313914127917,\n", + " 'stdev': 0.10047800378358163,\n", + " 'N': 200}" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAilElEQVR4nO3de3CU1f3H8c+G6CYiCROEXLgkSLkI2MiAkkC5lQk0XGoLCEgloYiXEamYoTQRKfCrGrSoKYJQWkhEhos2cqlhFKgE5CIjl3hpEUNJSAqJiEoWUDcgz+8Ph61bksDCPsnZ5f2aOTM+5znnyXcPMfnM2SfPOizLsgQAAGCwkIYuAAAA4HIILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA44U2dAH+cuHCBR0/flxNmjSRw+Fo6HIAAMAVsCxLp0+fVlxcnEJCat9HCZrAcvz4cbVu3bqhywAAAFehvLxcrVq1qvV80ASWJk2aSPr+BUdERDRwNQAA4Eq4XC61bt3a83u8NkETWC6+DRQREUFgAQAgwFzudg5uugUAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwXmhDFwDg+pKQWeDT+NK5Q22qBEAgYYcFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxfA4s27dv1/DhwxUXFyeHw6F169Z5nXc4HDW2P/7xj7VeMy8vr8Y53377rc8vCAAABB+fA8vZs2eVmJioBQsW1Hi+oqLCqy1btkwOh0MjR46s87oRERGXzA0LC/O1PAAAEIR8fg5LamqqUlNTaz0fExPjdbx+/XoNGDBAt956a53XdTgcl8wFAACQbL6H5bPPPlNBQYHuv//+y449c+aM4uPj1apVKw0bNkwHDhyoc7zb7ZbL5fJqAAAgONkaWF555RU1adJEI0aMqHNcp06dlJeXpw0bNmjVqlUKCwtT7969VVxcXOuc7OxsRUZGelrr1q39XT4AADCErYFl2bJl+tWvfnXZe1GSkpJ03333KTExUX369NFrr72mDh066KWXXqp1TlZWlqqqqjytvLzc3+UDAABD2PZZQu+++64OHTqkNWvW+Dw3JCREd955Z507LE6nU06n81pKBAAAAcK2HZalS5eqe/fuSkxM9HmuZVkqKipSbGysDZUBAIBA4/MOy5kzZ3T48GHPcUlJiYqKihQVFaU2bdpIklwul15//XU9//zzNV4jLS1NLVu2VHZ2tiRpzpw5SkpKUvv27eVyuTR//nwVFRVp4cKFV/OaAABAkPE5sOzdu1cDBgzwHGdkZEiS0tPTlZeXJ0lavXq1LMvSvffeW+M1ysrKFBLy382dU6dO6cEHH1RlZaUiIyPVrVs3bd++XXfddZev5QEAgCDksCzLaugi/MHlcikyMlJVVVWKiIho6HIA1CIhs8Cn8aVzh9pUCQATXOnvbz5LCAAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHihDV0AANSHhMwCn8aXzh1qUyUArgY7LAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADj+RxYtm/fruHDhysuLk4Oh0Pr1q3zOj9hwgQ5HA6vlpSUdNnr5ufnq3PnznI6nercubPWrl3ra2kAACBI+RxYzp49q8TERC1YsKDWMT/72c9UUVHhaRs3bqzzmrt379aYMWM0fvx4ffDBBxo/frxGjx6tPXv2+FoeAAAIQqG+TkhNTVVqamqdY5xOp2JiYq74mjk5OUpJSVFWVpYkKSsrS9u2bVNOTo5WrVrla4kAACDI2HIPS2FhoVq0aKEOHTrogQce0IkTJ+ocv3v3bg0aNMirb/Dgwdq1a1etc9xut1wul1cDAADByecdlstJTU3VPffco/j4eJWUlGjmzJn66U9/qn379snpdNY4p7KyUtHR0V590dHRqqysrPXrZGdna86cOX6tHQCCSUJmgU/jS+cOtakS4Nr5PbCMGTPG899du3ZVjx49FB8fr4KCAo0YMaLWeQ6Hw+vYsqxL+n4oKytLGRkZnmOXy6XWrVtfQ+UAAMBUfg8s/ys2Nlbx8fEqLi6udUxMTMwluyknTpy4ZNflh5xOZ607NgAAILjY/hyWL774QuXl5YqNja11THJysjZv3uzVt2nTJvXq1cvu8gAAQADweYflzJkzOnz4sOe4pKRERUVFioqKUlRUlGbPnq2RI0cqNjZWpaWleuKJJ3TLLbfol7/8pWdOWlqaWrZsqezsbEnSY489pr59++rZZ5/V3XffrfXr12vLli3asWOHH14iAAAIdD4Hlr1792rAgAGe44v3kaSnp2vRokX66KOPtHz5cp06dUqxsbEaMGCA1qxZoyZNmnjmlJWVKSTkv5s7vXr10urVq/Xkk09q5syZateundasWaOePXtey2sDAABBwufA0r9/f1mWVev5t99++7LXKCwsvKRv1KhRGjVqlK/lAACA6wCfJQQAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeD5/WjMAXG8SMgt8nlM6d6gNlQDXL3ZYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHihDV0AAFyphMwCn8aXzh1qUyUA6hs7LAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjOdzYNm+fbuGDx+uuLg4ORwOrVu3znPu3Llz+t3vfqfbb79djRs3VlxcnNLS0nT8+PE6r5mXlyeHw3FJ+/bbb31+QQAAIPj4HFjOnj2rxMRELViw4JJzX3/9tfbv36+ZM2dq//79euONN/Tpp5/q5z//+WWvGxERoYqKCq8WFhbma3kAACAI+fwcltTUVKWmptZ4LjIyUps3b/bqe+mll3TXXXeprKxMbdq0qfW6DodDMTExvpYDAACuA7bfw1JVVSWHw6GmTZvWOe7MmTOKj49Xq1atNGzYMB04cKDO8W63Wy6Xy6sBAIDgZGtg+fbbb5WZmalx48YpIiKi1nGdOnVSXl6eNmzYoFWrViksLEy9e/dWcXFxrXOys7MVGRnpaa1bt7bjJQAAAAPYFljOnTunsWPH6sKFC3r55ZfrHJuUlKT77rtPiYmJ6tOnj1577TV16NBBL730Uq1zsrKyVFVV5Wnl5eX+fgkAAMAQtnyW0Llz5zR69GiVlJTonXfeqXN3pSYhISG6884769xhcTqdcjqd11oqAAAIAH7fYbkYVoqLi7VlyxY1a9bM52tYlqWioiLFxsb6uzwAABCAfN5hOXPmjA4fPuw5LikpUVFRkaKiohQXF6dRo0Zp//79evPNN/Xdd9+psrJSkhQVFaUbb7xRkpSWlqaWLVsqOztbkjRnzhwlJSWpffv2crlcmj9/voqKirRw4UJ/vEYAABDgfA4se/fu1YABAzzHGRkZkqT09HTNnj1bGzZskCTdcccdXvO2bt2q/v37S5LKysoUEvLfzZ1Tp07pwQcfVGVlpSIjI9WtWzdt375dd911l6/lAQCAIORzYOnfv78sy6r1fF3nLiosLPQ6fvHFF/Xiiy/6WgoAALhO8FlCAADAeAQWAABgPAILAAAwni3PYQEA+EdCZoFP40vnDrWpEqBhscMCAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMbzObBs375dw4cPV1xcnBwOh9atW+d13rIszZ49W3FxcQoPD1f//v31z3/+87LXzc/PV+fOneV0OtW5c2etXbvW19IAAECQ8jmwnD17VomJiVqwYEGN55977jm98MILWrBggd5//33FxMQoJSVFp0+frvWau3fv1pgxYzR+/Hh98MEHGj9+vEaPHq09e/b4Wh4AAAhCob5OSE1NVWpqao3nLMtSTk6OZsyYoREjRkiSXnnlFUVHR2vlypV66KGHapyXk5OjlJQUZWVlSZKysrK0bds25eTkaNWqVb6WCAAAgoxf72EpKSlRZWWlBg0a5OlzOp3q16+fdu3aVeu83bt3e82RpMGDB9c5x+12y+VyeTUAABCcfN5hqUtlZaUkKTo62qs/OjpaR48erXNeTXMuXq8m2dnZmjNnzjVUC1zfEjILfBpfOneoX+YCwNWw5a+EHA6H17FlWZf0XeucrKwsVVVVeVp5efnVFwwAAIzm1x2WmJgYSd/vmMTGxnr6T5w4cckOyv/O+9/dlMvNcTqdcjqd11gxAAAIBH7dYWnbtq1iYmK0efNmT191dbW2bdumXr161TovOTnZa44kbdq0qc45AADg+uHzDsuZM2d0+PBhz3FJSYmKiooUFRWlNm3aaOrUqXrmmWfUvn17tW/fXs8884xuuukmjRs3zjMnLS1NLVu2VHZ2tiTpscceU9++ffXss8/q7rvv1vr167Vlyxbt2LHDDy8RAAAEOp8Dy969ezVgwADPcUZGhiQpPT1deXl5mj59ur755hs98sgj+uqrr9SzZ09t2rRJTZo08cwpKytTSMh/N3d69eql1atX68knn9TMmTPVrl07rVmzRj179ryW1wYAAIKEz4Glf//+siyr1vMOh0OzZ8/W7Nmzax1TWFh4Sd+oUaM0atQoX8sBAADXAT5LCAAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwXmhDFwAEuoTMAp/nlM4detXz/TUXqAvfWzANOywAAMB4BBYAAGA8AgsAADAegQUAABjP74ElISFBDofjkjZ58uQaxxcWFtY4/pNPPvF3aQAAIED5/a+E3n//fX333Xee448//lgpKSm655576px36NAhRUREeI6bN2/u79IAAECA8ntg+d+gMXfuXLVr1079+vWrc16LFi3UtGlTf5cDAACCgK33sFRXV2vFihWaOHGiHA5HnWO7deum2NhYDRw4UFu3brWzLAAAEGBsfXDcunXrdOrUKU2YMKHWMbGxsVqyZIm6d+8ut9utV199VQMHDlRhYaH69u1b6zy32y232+05drlc/iwdAAAYxNbAsnTpUqWmpiouLq7WMR07dlTHjh09x8nJySovL9e8efPqDCzZ2dmaM2eOX+sFAABmsu0toaNHj2rLli2aNGmSz3OTkpJUXFxc55isrCxVVVV5Wnl5+dWWCgAADGfbDktubq5atGihoUN9/3yJAwcOKDY2ts4xTqdTTqfzassDAAABxJbAcuHCBeXm5io9PV2hod5fIisrS8eOHdPy5cslSTk5OUpISFCXLl08N+nm5+crPz/fjtIAAEAAsiWwbNmyRWVlZZo4ceIl5yoqKlRWVuY5rq6u1rRp03Ts2DGFh4erS5cuKigo0JAhQ+woDQAABCBbAsugQYNkWVaN5/Ly8ryOp0+frunTp9tRBgAACBJ8lhAAADAegQUAABiPwAIAAIxn64PjAABSQmaBT+NL5/r+OAiTXG+vF/WDHRYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYL7ShCwD8JSGzwKfxpXOH+mUuAMB+7LAAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGM/vgWX27NlyOBxeLSYmps4527ZtU/fu3RUWFqZbb71Vixcv9ndZAAAggIXacdEuXbpoy5YtnuNGjRrVOrakpERDhgzRAw88oBUrVmjnzp165JFH1Lx5c40cOdKO8gAAQICxJbCEhoZedlflosWLF6tNmzbKycmRJN12223au3ev5s2bR2ABAACSbLqHpbi4WHFxcWrbtq3Gjh2rI0eO1Dp29+7dGjRokFff4MGDtXfvXp07d86O8gAAQIDxe2Dp2bOnli9frrffflt/+ctfVFlZqV69eumLL76ocXxlZaWio6O9+qKjo3X+/HmdPHmy1q/jdrvlcrm8GgAACE5+f0soNTXV89+33367kpOT1a5dO73yyivKyMiocY7D4fA6tiyrxv4fys7O1pw5c/xQMQDAFAmZBT6NL5079Krn/u98mM32P2tu3Lixbr/9dhUXF9d4PiYmRpWVlV59J06cUGhoqJo1a1brdbOyslRVVeVp5eXlfq0bAACYw5abbn/I7Xbr4MGD6tOnT43nk5OT9fe//92rb9OmTerRo4duuOGGWq/rdDrldDr9WisAADCT33dYpk2bpm3btqmkpER79uzRqFGj5HK5lJ6eLun7nZG0tDTP+IcfflhHjx5VRkaGDh48qGXLlmnp0qWaNm2av0sDAAAByu87LP/5z39077336uTJk2revLmSkpL03nvvKT4+XpJUUVGhsrIyz/i2bdtq48aNevzxx7Vw4ULFxcVp/vz5/EkzAADw8HtgWb16dZ3n8/LyLunr16+f9u/f7+9SAABAkOCzhAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgPAILAAAwHoEFAAAYz++f1gwAQCBKyCzwaXzp3KE2VVI/Au31ssMCAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA44U2dAGoW0JmgU/jS+cOtakSAIAd+Dl/ZdhhAQAAxiOwAAAA4xFYAACA8QgsAADAeH4PLNnZ2brzzjvVpEkTtWjRQr/4xS906NChOucUFhbK4XBc0j755BN/lwcAAAKQ3wPLtm3bNHnyZL333nvavHmzzp8/r0GDBuns2bOXnXvo0CFVVFR4Wvv27f1dHgAACEB+/7Pmt956y+s4NzdXLVq00L59+9S3b98657Zo0UJNmzb1d0kAACDA2X4PS1VVlSQpKirqsmO7deum2NhYDRw4UFu3bq1zrNvtlsvl8moAACA42RpYLMtSRkaGfvKTn6hr1661jouNjdWSJUuUn5+vN954Qx07dtTAgQO1ffv2WudkZ2crMjLS01q3bm3HSwAAAAaw9Um3jz76qD788EPt2LGjznEdO3ZUx44dPcfJyckqLy/XvHnzan0bKSsrSxkZGZ5jl8tFaAEAIEjZtsMyZcoUbdiwQVu3blWrVq18np+UlKTi4uJazzudTkVERHg1AAAQnPy+w2JZlqZMmaK1a9eqsLBQbdu2varrHDhwQLGxsX6uDgAABCK/B5bJkydr5cqVWr9+vZo0aaLKykpJUmRkpMLDwyV9/3bOsWPHtHz5cklSTk6OEhIS1KVLF1VXV2vFihXKz89Xfn6+v8sDAAAByO+BZdGiRZKk/v37e/Xn5uZqwoQJkqSKigqVlZV5zlVXV2vatGk6duyYwsPD1aVLFxUUFGjIkCH+Lg8AAAQgW94Supy8vDyv4+nTp2v69On+LgUAAAQJPksIAAAYj8ACAACMR2ABAADGs/XBccEiIbPAp/Glc4faVAkAwEQN9Xvievr9xA4LAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMYjsAAAAOMRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjBfa0AXAPgmZBT6NL507NKDnAgCCFzssAADAeAQWAABgPAILAAAwHoEFAAAYj8ACAACMR2ABAADGI7AAAADjEVgAAIDxCCwAAMB4BBYAAGA8AgsAADAegQUAABiPwAIAAIxHYAEAAMazLbC8/PLLatu2rcLCwtS9e3e9++67dY7ftm2bunfvrrCwMN16661avHixXaUBAIAAY0tgWbNmjaZOnaoZM2bowIED6tOnj1JTU1VWVlbj+JKSEg0ZMkR9+vTRgQMH9MQTT+g3v/mN8vPz7SgPAAAEGFsCywsvvKD7779fkyZN0m233aacnBy1bt1aixYtqnH84sWL1aZNG+Xk5Oi2227TpEmTNHHiRM2bN8+O8gAAQIAJ9fcFq6urtW/fPmVmZnr1Dxo0SLt27apxzu7duzVo0CCvvsGDB2vp0qU6d+6cbrjhhkvmuN1uud1uz3FVVZUkyeVyXetLuMQF99c+jfdnDdfytZlr5tyG/NrMrZ+5Dfm1meubQKzbhLn+dPG6lmXVPdDys2PHjlmSrJ07d3r1P/3001aHDh1qnNO+fXvr6aef9urbuXOnJck6fvx4jXNmzZplSaLRaDQajRYErby8vM584fcdloscDofXsWVZl/RdbnxN/RdlZWUpIyPDc3zhwgV9+eWXatasmWeOy+VS69atVV5eroiIiKt6Hagba2w/1rh+sM72Y43tF4hrbFmWTp8+rbi4uDrH+T2w3HLLLWrUqJEqKyu9+k+cOKHo6Oga58TExNQ4PjQ0VM2aNatxjtPplNPp9Opr2rRpjWMjIiIC5h8uULHG9mON6wfrbD/W2H6BtsaRkZGXHeP3m25vvPFGde/eXZs3b/bq37x5s3r16lXjnOTk5EvGb9q0ST169Kjx/hUAAHB9seWvhDIyMvTXv/5Vy5Yt08GDB/X444+rrKxMDz/8sKTv385JS0vzjH/44Yd19OhRZWRk6ODBg1q2bJmWLl2qadOm2VEeAAAIMLbcwzJmzBh98cUX+r//+z9VVFSoa9eu2rhxo+Lj4yVJFRUVXs9kadu2rTZu3KjHH39cCxcuVFxcnObPn6+RI0deUx1Op1OzZs265K0j+A9rbD/WuH6wzvZjje0XzGvssKzL/R0RAABAw+KzhAAAgPEILAAAwHgEFgAAYDwCCwAAMF7ABpaEhAQ5HI5L2uTJk2scX1FRoXHjxqljx44KCQnR1KlT67fgAOTrGr/xxhtKSUlR8+bNFRERoeTkZL399tv1XHVg8XWNd+zYod69e6tZs2YKDw9Xp06d9OKLL9Zz1YHH13X+oZ07dyo0NFR33HGH/YUGMF/XuLCwsMbxn3zyST1XHjiu5vvY7XZrxowZio+Pl9PpVLt27bRs2bJ6rNp/bHs0v93ef/99fffdd57jjz/+WCkpKbrnnntqHO92u9W8eXPNmDGDH/BXyNc13r59u1JSUvTMM8+oadOmys3N1fDhw7Vnzx5169atvsoOKL6ucePGjfXoo4/qxz/+sRo3bqwdO3booYceUuPGjfXggw/WV9kBx9d1vqiqqkppaWkaOHCgPvvsM7vLDGhXu8aHDh3yeiJr8+bNbasx0F3NGo8ePVqfffaZli5dqh/96Ec6ceKEzp8/Xx/l+l3Q/Fnz1KlT9eabb6q4uLjOzyySpP79++uOO+5QTk5O/RQXJHxZ44u6dOmiMWPG6Pe//73N1QWHq1njESNGqHHjxnr11Vdtri54XOk6jx07Vu3bt1ejRo20bt06FRUV1V+RAe5ya1xYWKgBAwboq6++qvVjVVC3y63xW2+9pbFjx+rIkSOKiopqgAr9K2DfEvqh6upqrVixQhMnTrziH/LwzdWs8YULF3T69Omg+B+lPlzNGh84cEC7du1Sv379bK4ueFzpOufm5urf//63Zs2aVY/VBQdfvpe7deum2NhYDRw4UFu3bq2nCgPflazxhg0b1KNHDz333HNq2bKlOnTooGnTpumbb76p52r9I2DfEvqhdevW6dSpU5owYUJDlxK0rmaNn3/+eZ09e1ajR4+2r7Ag4ssat2rVSp9//rnOnz+v2bNna9KkSfYXGCSuZJ2Li4uVmZmpd999V6GhQfFjsl5dyRrHxsZqyZIl6t69u9xut1599VUNHDhQhYWF6tu3b/0VG6CuZI2PHDmiHTt2KCwsTGvXrtXJkyf1yCOP6MsvvwzM+1isIDBo0CBr2LBhVzy+X79+1mOPPWZfQUHI1zVeuXKlddNNN1mbN2+2sarg4ssaHzlyxPrwww+tJUuWWFFRUdbKlSttri54XG6dz58/b/Xo0cNatGiRp2/WrFlWYmJiPVQXHHz9eXHRsGHDrOHDh9tQUfC5kjVOSUmxwsLCrFOnTnn68vPzLYfDYX399dd2l+h3AR9YSktLrZCQEGvdunVXPIfA4htf13j16tVWeHi49eabb9pcWfC4mu/ji/7whz9YHTp0sKGq4HMl6/zVV19ZkqxGjRp5msPh8PT94x//qMeKA8+1fC8/9dRTVqdOnWyoKrhc6RqnpaVZ7dq18+r717/+ZUmyPv30UztLtEXA73Xm5uaqRYsWGjp0aEOXErR8WeNVq1Zp4sSJWrVqFf8mPriW72PLsuR2u22oKvhcyTpHREToo48+8up7+eWX9c477+hvf/ub2rZta3eZAe1avpcPHDig2NhYG6oKLle6xr1799brr7+uM2fO6Oabb5YkffrppwoJCVGrVq3qo1S/CujAcuHCBeXm5io9Pf2S95mzsrJ07NgxLV++3NN38Q7/M2fO6PPPP1dRUZFuvPFGde7cuT7LDii+rPGqVauUlpamP/3pT0pKSlJlZaUkKTw8XJGRkfVee6DwZY0XLlyoNm3aqFOnTpK+fy7LvHnzNGXKlHqvO9Bc6TqHhISoa9euXudbtGihsLCwS/rhzZfv5ZycHCUkJKhLly6eG0jz8/OVn5/fEKUHDF/WeNy4cfrDH/6gX//615ozZ45Onjyp3/72t5o4caLCw8MbovxrEtCBZcuWLSorK9PEiRMvOVdRUaGysjKvvh8+C2Tfvn1auXKl4uPjVVpaanepAcuXNf7zn/+s8+fPa/LkyV4PMkpPT1deXl59lBuQfFnjCxcuKCsrSyUlJQoNDVW7du00d+5cPfTQQ/VZckDy9ecFfOfLGldXV2vatGk6duyYwsPD1aVLFxUUFGjIkCH1WXLA8WWNb775Zm3evFlTpkxRjx491KxZM40ePVpPPfVUfZbsN0HzHBYAABC8guI5LAAAILgRWAAAgPEILAAAwHgEFgAAYDwCCwAAMB6BBQAAGI/AAgAAjEdgAQAAxiOwAAAA4xFYAACA8QgsAADAeAQWAABgvP8H6XaZoF5z+kYAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "report([mean_cluster_size(random_grid(200, 200)) for _ in range(N)])" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "b580bc66-9bc2-406c-8ce4-806b9c020dc7", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'mean': 7.469780220918766,\n", + " 'max': 7.719358435543357,\n", + " 'min': 7.246960302761897,\n", + " 'stdev': 0.06989052018543407,\n", + " 'N': 200}" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAW/0lEQVR4nO3df6zVdf3A8deBO4+IYAOFe28gMIXZwKxJw7AUUiliboXmry0htWj+2Bgzx40/vKQBc2m4kW5uRmpRrjSyMBMz8AezodNlVop5kZuCpOK9wNddhpzvH847b6DcA+e8LufweGyfzfM5n3M/L987HJ77nMM9hVKpVAoAgCT9+noAAODwIj4AgFTiAwBIJT4AgFTiAwBIJT4AgFTiAwBIJT4AgFQNfT3A/9qzZ0+8/vrrMWjQoCgUCn09DgDQC6VSKbZv3x7Nzc3Rr9/HX9s45OLj9ddfj5EjR/b1GADAAWhvb48RI0Z87DGHXHwMGjQoIt4ffvDgwX08DQDQG52dnTFy5Mjuv8c/ziEXHx+81TJ48GDxAQA1pjcfmfCBUwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFI19PUAQO0aPX9VWcdvXDKjSpMAtcSVDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAgVVnxsXjx4vjc5z4XgwYNimHDhsXXvva1ePHFF3scUyqVorW1NZqbm2PAgAExZcqUeOGFFyo6NABQu8qKj7Vr18ZVV10VTz31VKxevTp2794d06ZNi507d3Yfc9NNN8Utt9wSy5Yti/Xr10djY2Occ845sX379ooPDwDUnoZyDn7ooYd63F6+fHkMGzYsnnnmmTjjjDOiVCrF0qVLY8GCBTFz5syIiLjrrrti+PDhsWLFipgzZ07lJgcAatJBfeajo6MjIiKGDBkSERFtbW2xZcuWmDZtWvcxxWIxzjzzzFi3bt0+f0ZXV1d0dnb22ACA+lXWlY8PK5VKMW/evPjCF74QEyZMiIiILVu2RETE8OHDexw7fPjwePXVV/f5cxYvXhwLFy480DGgLoyev6rsx2xcMqNPzl2p8wKHrwO+8nH11VfH3/72t/jlL3+5132FQqHH7VKptNe+D7S0tERHR0f31t7efqAjAQA14ICufFxzzTXxwAMPxGOPPRYjRozo3t/Y2BgR718BaWpq6t6/devWva6GfKBYLEaxWDyQMQCAGlTWlY9SqRRXX3113H///fHoo4/GmDFjetw/ZsyYaGxsjNWrV3fv27VrV6xduzYmT55cmYkBgJpW1pWPq666KlasWBG/+93vYtCgQd2f8TjmmGNiwIABUSgUYu7cubFo0aIYO3ZsjB07NhYtWhRHHXVUXHLJJVX5HwAAaktZ8XH77bdHRMSUKVN67F++fHnMnj07IiKuu+66ePfdd+PKK6+Mbdu2xaRJk+Lhhx+OQYMGVWRgAKC2lRUfpVJpv8cUCoVobW2N1tbWA50JAKhjvtsFAEglPgCAVOIDAEglPgCAVOIDAEglPgCAVOIDAEglPgCAVAf0xXLAoWX0/FVlHb9xyYwqTQKwf658AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpxAcAkEp8AACpGvp6AKgXo+evKuv4jUtmVGkSgEObKx8AQCrxAQCkEh8AQCrxAQCkEh8AQCrxAQCkEh8AQCrxAQCkEh8AQCrxAQCkEh8AQKqy4+Oxxx6Lc889N5qbm6NQKMTKlSt73D979uwoFAo9ttNOO61S8wIANa7s+Ni5c2eccsopsWzZso885itf+Ups3ry5e3vwwQcPakgAoH6U/a2206dPj+nTp3/sMcViMRobGw94KACgflXlMx9r1qyJYcOGxbhx4+Lb3/52bN269SOP7erqis7Ozh4bAFC/Kh4f06dPj1/84hfx6KOPxs033xzr16+PL33pS9HV1bXP4xcvXhzHHHNM9zZy5MhKjwQAHELKfttlfy688MLu/54wYUJMnDgxRo0aFatWrYqZM2fudXxLS0vMmzev+3ZnZ6cAAYA6VvH4+F9NTU0xatSo2LBhwz7vLxaLUSwWqz0GAHCIqPrv+Xjrrbeivb09mpqaqn0qAKAGlH3lY8eOHfHyyy93325ra4vnnnsuhgwZEkOGDInW1tY477zzoqmpKTZu3Bjf//7349hjj42vf/3rFR0cAKhNZcfH008/HVOnTu2+/cHnNWbNmhW33357PP/883H33XfHO++8E01NTTF16tS49957Y9CgQZWbGgCoWWXHx5QpU6JUKn3k/X/6058OaiAAoL75bhcAIJX4AABSiQ8AIJX4AABSiQ8AIJX4AABSiQ8AIJX4AABSVf2L5QD2ZfT8VWUdv3HJjIo8Fuh7rnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQquz4eOyxx+Lcc8+N5ubmKBQKsXLlyh73l0qlaG1tjebm5hgwYEBMmTIlXnjhhUrNCwDUuLLjY+fOnXHKKafEsmXL9nn/TTfdFLfcckssW7Ys1q9fH42NjXHOOefE9u3bD3pYAKD2NZT7gOnTp8f06dP3eV+pVIqlS5fGggULYubMmRERcdddd8Xw4cNjxYoVMWfOnIObFgCoeRX9zEdbW1ts2bIlpk2b1r2vWCzGmWeeGevWrdvnY7q6uqKzs7PHBgDUr7KvfHycLVu2RETE8OHDe+wfPnx4vPrqq/t8zOLFi2PhwoWVHAPgI42ev6qs4zcumVGlSfIcjv/PHNqq8q9dCoVCj9ulUmmvfR9oaWmJjo6O7q29vb0aIwEAh4iKXvlobGyMiPevgDQ1NXXv37p1615XQz5QLBajWCxWcgwA4BBW0SsfY8aMicbGxli9enX3vl27dsXatWtj8uTJlTwVAFCjyr7ysWPHjnj55Ze7b7e1tcVzzz0XQ4YMieOPPz7mzp0bixYtirFjx8bYsWNj0aJFcdRRR8Ull1xS0cEBgNpUdnw8/fTTMXXq1O7b8+bNi4iIWbNmxc9+9rO47rrr4t13340rr7wytm3bFpMmTYqHH344Bg0aVLmpAYCaVXZ8TJkyJUql0kfeXygUorW1NVpbWw9mLgCgTvluFwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAglfgAAFKJDwAgVUW/1RYqZfT8VWUdv3HJjCpNAkClufIBAKQSHwBAKvEBAKQSHwBAKvEBAKQSHwBAKvEBAKQSHwBAKvEBAKQSHwBAKvEBAKQSHwBAKvEBAKQSHwBAqoa+HgCA/Rs9f1VZx29cMqNKk8DBc+UDAEglPgCAVOIDAEglPgCAVOIDAEglPgCAVOIDAEglPgCAVOIDAEglPgCAVOIDAEglPgCAVOIDAEglPgCAVA19PQBArfC19lAZrnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQqqGvB4BDia9MB6g+Vz4AgFTiAwBIJT4AgFTiAwBIJT4AgFTiAwBIJT4AgFTiAwBIJT4AgFTiAwBIJT4AgFQVj4/W1tYoFAo9tsbGxkqfBgCoUVX5Yrnx48fHI4880n27f//+1TgNAFCDqhIfDQ0NrnYAAPtUlc98bNiwIZqbm2PMmDFx0UUXxSuvvPKRx3Z1dUVnZ2ePDQCoXxW/8jFp0qS4++67Y9y4cfHGG2/EjTfeGJMnT44XXnghhg4dutfxixcvjoULF1Z6DA5jo+evKuv4jUtmVGkS6MlzE95X8Ssf06dPj/POOy9OPvnkOPvss2PVqvf/sN111137PL6lpSU6Ojq6t/b29kqPBAAcQqrymY8PGzhwYJx88smxYcOGfd5fLBajWCxWewwA4BBR9d/z0dXVFf/85z+jqamp2qcCAGpAxePj2muvjbVr10ZbW1v89a9/jfPPPz86Oztj1qxZlT4VAFCDKv62y3/+85+4+OKL480334zjjjsuTjvttHjqqadi1KhRlT4VAFCDKh4fv/rVryr9IwGAOuK7XQCAVOIDAEglPgCAVOIDAEglPgCAVOIDAEglPgCAVOIDAEglPgCAVOIDAEglPgCAVOIDAEglPgCAVOIDAEglPgCAVOIDAEglPgCAVOIDAEglPgCAVOIDAEglPgCAVOIDAEglPgCAVA19PQAAh67R81eVdfzGJTOqNAn1xJUPACCV+AAAUokPACCV+AAAUokPACCV+AAAUokPACCV+AAAUokPACCV+AAAUokPACCV+AAAUokPACCV+AAAUokPACBVQ18PwP6Nnr+q7MdsXDKjT8794fMezGMBqF+ufAAAqcQHAJBKfAAAqcQHAJBKfAAAqcQHAJBKfAAAqcQHAJBKfAAAqcQHAJBKfAAAqcQHAJBKfAAAqcQHAJCqoa8HyOZr3gFyHMzrbV+9Vpd73kqe+2DU2t9trnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQSnwAAKnEBwCQqmrxcdttt8WYMWPiyCOPjFNPPTUef/zxap0KAKghVYmPe++9N+bOnRsLFiyIZ599Nr74xS/G9OnTY9OmTdU4HQBQQ6oSH7fccktcfvnlccUVV8SnPvWpWLp0aYwcOTJuv/32apwOAKghFf9W2127dsUzzzwT8+fP77F/2rRpsW7dur2O7+rqiq6uru7bHR0dERHR2dlZ6dEiImJP1/+VdXy15ihHuTNHVG7ug1kvj63eY/vy3B57YI/ty3Mfbo89GH35enswDoW/2z74maVSaf8HlyrstddeK0VE6cknn+yx/4c//GFp3Lhxex1//fXXlyLCZrPZbDZbHWzt7e37bYWKX/n4QKFQ6HG7VCrttS8ioqWlJebNm9d9e8+ePfH222/H0KFD93k8H6+zszNGjhwZ7e3tMXjw4L4e57Bj/fuW9e9b1r9v9fX6l0ql2L59ezQ3N+/32IrHx7HHHhv9+/ePLVu29Ni/devWGD58+F7HF4vFKBaLPfZ94hOfqPRYh53Bgwf7w9+HrH/fsv59y/r3rb5c/2OOOaZXx1X8A6dHHHFEnHrqqbF69eoe+1evXh2TJ0+u9OkAgBpTlbdd5s2bF9/85jdj4sSJ8fnPfz7uuOOO2LRpU3z3u9+txukAgBpSlfi48MIL46233oof/OAHsXnz5pgwYUI8+OCDMWrUqGqcjg8pFotx/fXX7/VWFjmsf9+y/n3L+vetWlr/QqnUm38TAwBQGb7bBQBIJT4AgFTiAwBIJT4AgFTio4aMHj06CoXCXttVV121z+OfeOKJOP3002Po0KExYMCAOOmkk+LHP/5x8tT1o9z1/7Ann3wyGhoa4jOf+Uz1B61T5a7/mjVr9nn8v/71r+TJ68OBPP+7urpiwYIFMWrUqCgWi3HCCSfET3/608Sp60e56z979ux9Hj9+/Pjkyfetar9encpbv359vPfee923//73v8c555wT3/jGN/Z5/MCBA+Pqq6+OT3/60zFw4MB44oknYs6cOTFw4MD4zne+kzV23Sh3/T/Q0dERl156aZx11lnxxhtvVHvMunWg6//iiy/2+G2Pxx13XNVmrGcHsv4XXHBBvPHGG3HnnXfGiSeeGFu3bo3du3dnjFt3yl3/W2+9NZYsWdJ9e/fu3XHKKafs989LFv/UtobNnTs3/vCHP8SGDRt6/T04M2fOjIEDB8Y999xT5enqX2/X/6KLLoqxY8dG//79Y+XKlfHcc8/lDVnH9rf+a9asialTp8a2bdt8ZUMV7G/9H3roobjooovilVdeiSFDhvTBhPWt3Nf/lStXxsyZM6Otre2Q+J1b3napUbt27Yqf//zncdlll/U6PJ599tlYt25dnHnmmVWerv71dv2XL18e//73v+P6669PnK7+lfP8/+xnPxtNTU1x1llnxV/+8pekCetbb9b/gQceiIkTJ8ZNN90Un/zkJ2PcuHFx7bXXxrvvvps8bf05kNf/O++8M84+++xDIjwivO1Ss1auXBnvvPNOzJ49e7/HjhgxIv773//G7t27o7W1Na644orqD1jnerP+GzZsiPnz58fjjz8eDQ3+qFVSb9a/qakp7rjjjjj11FOjq6sr7rnnnjjrrLNizZo1ccYZZ+QNW4d6s/6vvPJKPPHEE3HkkUfGb3/723jzzTfjyiuvjLffftvnPg5SOa//ERGbN2+OP/7xj7FixYrqDlYGb7vUqC9/+ctxxBFHxO9///v9HtvW1hY7duyIp556KubPnx/Lli2Liy++OGHK+rW/9X/vvffitNNOi8svv7z7O41aW1u97VIh5Tz/P+zcc8+NQqEQDzzwQJUmOzz0Zv2nTZsWjz/+eGzZsqX7m07vv//+OP/882Pnzp0xYMCArHHrTrnP/8WLF8fNN98cr7/+ehxxxBFVnq6XStScjRs3lvr161dauXJl2Y+94YYbSuPGjavCVIeP3qz/tm3bShFR6t+/f/dWKBS69/35z39OnLi+HMzz/8YbbyyddNJJVZjq8NHb9b/00ktLJ5xwQo99//jHP0oRUXrppZeqOWJdK/f5v2fPntKJJ55Ymjt3bpUnK49rwTVo+fLlMWzYsJgxY0bZjy2VStHV1VWFqQ4fvVn/wYMHx/PPP99j32233RaPPvpo/OY3v4kxY8ZUe8y6dTDP/2effTaampqqMNXho7frf/rpp8evf/3r2LFjRxx99NEREfHSSy9Fv379YsSIERmj1qVyn/9r166Nl19+OS6//PIqT1Ye8VFj9uzZE8uXL49Zs2bt9TmClpaWeO211+Luu++OiIif/OQncfzxx8dJJ50UEe//3o8f/ehHcc0116TPXS96u/79+vWLCRMm9Lh/2LBhceSRR+61n94r5/m/dOnSGD16dIwfP777A3r33Xdf3HfffX0xel0oZ/0vueSSuOGGG+Jb3/pWLFy4MN5888343ve+F5dddpm3XA5QOev/gTvvvDMmTZp0yL3uiI8a88gjj8SmTZvisssu2+u+zZs3x6ZNm7pv79mzJ1paWqKtrS0aGhrihBNOiCVLlsScOXMyR64r5aw/lVfO+u/atSuuvfbaeO2112LAgAExfvz4WLVqVXz1q1/NHLmulLP+Rx99dKxevTquueaamDhxYgwdOjQuuOCCuPHGGzNHrivlvv50dHTEfffdF7feemvWiL3mA6cAQCq/5wMASCU+AIBU4gMASCU+AIBU4gMASCU+AIBU4gMASCU+AIBU4gMASCU+AIBU4gMASCU+AIBU/w80z+BJw/qnVgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "report([mean_cluster_size(random_grid(300, 300)) for _ in range(N)])" + ] + }, + { + "cell_type": "markdown", + "id": "6d20f32a-8c7a-4456-aa41-d2989e33fd57", + "metadata": {}, + "source": [ + "I think that what's happening is that the clusters that are near the edge of the grid get arbitrarily cut off, and since the edge is a smaller percentage of the larger grid, the larger grid has a larger mean cluster size, one that is a better representative of what would happen on an infinite grid. But I can't say what the mean converges to." + ] + }, + { + "cell_type": "markdown", + "id": "f3f21c52-8baf-4e71-a33c-e6f4cefd0877", + "metadata": {}, + "source": [ + "# My First Random Cluster\n", + "\n", + "The great thing about an *n* × 1 grid is that the clusters can't interfere with each other. Going left-to-right, once one cluster ends, the next one starts, and there's no looping back. But in an arbitrary 2D array there is looping back. So I want to get a feel for what the \"first\" cluster on an infinite grid looks like: start by placing the first color at the origin square, then recursively place random colors in neighboring squares, extending the cluster, but stop when the cluster is completely surrounded by squares of the second color. " + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "ee2a29b8-7ec6-4558-b86a-5933ddc716f0", + "metadata": {}, + "outputs": [], + "source": [ + "def random_cluster(colors='@.') -> Grid:\n", + " \"\"\"Create just enough of a grid to cover a single random cluster of the first color.\"\"\"\n", + " c = colors[0]\n", + " cluster = {(0, 0): c}\n", + " Q = neighbors((0, 0))\n", + " while Q:\n", + " sq = Q.pop()\n", + " if sq not in cluster:\n", + " cluster[sq] = random.choice(colors)\n", + " if cluster[sq] == c:\n", + " Q.extend(neighbors(sq))\n", + " return cluster" + ] + }, + { + "cell_type": "markdown", + "id": "58e3e7bc-2576-4332-921c-04e02da6490b", + "metadata": {}, + "source": [ + "I'll define `show` to print a picture of the cluster:" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "ee6683dd-7b4c-4e4c-b28e-1e8075cc07aa", + "metadata": {}, + "outputs": [], + "source": [ + "def show(cluster: Grid, color='@') -> int:\n", + " \"\"\"Print a representation of the grid and return the size of the one cluster.\"\"\"\n", + " xs = sorted({x for (x, y) in cluster})\n", + " ys = sorted({y for (x, y) in cluster})\n", + " for y in ys:\n", + " print(*[cluster.get((x, y), ' ') for x in xs])\n", + " return cluster_size(cluster, color)\n", + "\n", + "def cluster_size(cluster, color='@') -> int: return Counter(cluster.values())[color]" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "1827747e-76ac-496e-ad3b-b7cbf8b5201a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " . \n", + ". @ . \n", + ". @ @ . . \n", + " . @ @ @ .\n", + " . @ . . \n", + " . @ @ . \n", + " . @ @ .\n", + " . @ . \n", + " . \n" + ] + }, + { + "data": { + "text/plain": [ + "12" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "show(random_cluster())" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "fd429e22-b8bd-475e-9bb0-cbdac2d57c0e", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " . . . . \n", + " . @ @ . . @ @ . \n", + ". @ @ @ @ @ @ @ @ . \n", + " . @ @ . . @ . @ @ . \n", + " . . . @ . . @ . \n", + " . @ . . @ @ @ . @ . \n", + " . @ . . @ @ . . \n", + " . @ @ @ . @ @ . \n", + " . @ @ @ @ . @ . \n", + " . @ @ @ . @ . \n", + " . @ @ . @ @ . \n", + " . @ @ @ @ . \n", + " . @ . @ . \n", + " . @ @ . . \n", + " . @ @ @ . @ .\n", + " . . @ @ @ .\n", + " . @ . . \n", + " . \n" + ] + }, + { + "data": { + "text/plain": [ + "61" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "show(random_cluster())" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "13ba443e-c5b9-4897-9075-0833f92081ae", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " . \n", + " . @ . . . . \n", + " . @ @ @ @ @ . \n", + " . . @ @ . @ @ @ . \n", + " . . @ @ @ @ . @ . . \n", + " . @ . . @ . @ @ @ . \n", + " . @ @ @ . . . @ @ . \n", + " . @ @ . . @ . . . . . \n", + " . . @ @ . . . . . @ . . @ @ @ @ .\n", + " . . . @ @ @ . @ @ @ @ . . @ . . . . @ . . . \n", + " . @ @ @ @ @ . . . . @ @ @ @ @ . @ @ @ @ @ . \n", + ". @ @ @ @ @ @ @ @ @ @ @ . . . @ . @ . @ @ . \n", + " . @ @ @ . @ @ @ . @ . . @ @ @ @ @ . . \n", + " . . @ . @ . . . . @ @ @ @ @ @ . \n", + " . . . @ @ @ . . . . \n", + " . @ @ . . \n", + " . . \n" + ] + }, + { + "data": { + "text/plain": [ + "100" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "show(random_cluster())" + ] + }, + { + "cell_type": "markdown", + "id": "faf079b0-8534-4d92-aa8e-acd035dc5420", + "metadata": {}, + "source": [ + "We can see that there is a lof of variation in size and shape of the first clusters. We can also see that the shape will constrain other clusters: there are some small clusters of the second color that are trapped inside the cluster of the first color. Below I sample 10,000 first clusters and report on their sizes:" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "ab6ba895-aa51-41af-baf1-f09ee47e973f", + "metadata": {}, + "outputs": [], + "source": [ + "sizes = [cluster_size(random_cluster()) for _ in range(10_000)]" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "df7f6ae0-32aa-4002-8056-9de0dc2da07a", + "metadata": { + "collapsed": true, + "jupyter": { + "outputs_hidden": true + }, + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "{'mean': 59.3772, 'max': 703, 'min': 1, 'stdev': 75.3564687329362, 'N': 10000}" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjEAAAGdCAYAAADjWSL8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAhkElEQVR4nO3df2xV9f3H8deVtheo7RkF2+sdVap2CCswV1wp8ztQoEConTEZbHUdRgSRnw0QfsgfsCW2yDJQ08kQCajguj8UxyZWSsQqgwIWbiyITGOVIr2Uabkt2LVYPt8/Fk68FNEWsP3cPh/JSew579uez0kjz5zeHx5jjBEAAIBlruvoEwAAAGgPIgYAAFiJiAEAAFYiYgAAgJWIGAAAYCUiBgAAWImIAQAAViJiAACAlaI6+gSulfPnz+vEiROKi4uTx+Pp6NMBAADfgTFGDQ0N8vv9uu66y99ridiIOXHihJKTkzv6NAAAQDtUV1erb9++l52J2IiJi4uT9L+LEB8f38FnAwAAvov6+nolJye7/45fTsRGzIU/IcXHxxMxAABY5rs8FYQn9gIAACsRMQAAwEpEDAAAsBIRAwAArETEAAAAKxExAADASkQMAACwEhEDAACsRMQAAAArETEAAMBKRAwAALASEQMAAKxExAAAACsRMQAAwEpRHX0Ctuq3+LU2zX+yYsI1OhMAALom7sQAAAArETEAAMBKRAwAALASEQMAAKxExAAAACsRMQAAwEpEDAAAsBIRAwAArETEAAAAKxExAADASkQMAACwEhEDAACsRMQAAAArETEAAMBKRAwAALASEQMAAKxExAAAACsRMQAAwEpEDAAAsBIRAwAArETEAAAAKxExAADASkQMAACw0hVFTGFhoTwej/Lz8919xhgtX75cfr9fPXr00MiRI3X48OGwxzU1NWn27Nnq06ePYmNjlZOTo+PHj4fN1NXVKS8vT47jyHEc5eXl6fTp01dyugAAIIK0O2L279+vZ599VoMHDw7bv3LlSq1atUpFRUXav3+/fD6fxowZo4aGBncmPz9fW7ZsUXFxsXbt2qUzZ84oOztbLS0t7kxubq4CgYBKSkpUUlKiQCCgvLy89p4uAACIMO2KmDNnzuiBBx7QunXr1KtXL3e/MUZPPvmkli5dqvvvv19paWl6/vnn9eWXX+qll16SJIVCIa1fv15/+tOfNHr0aN1xxx3atGmTKisrtWPHDknSkSNHVFJSoueee06ZmZnKzMzUunXr9M9//lNHjx69CssGAAC2a1fEzJw5UxMmTNDo0aPD9ldVVSkYDCorK8vd5/V6NWLECO3evVuSVFFRoXPnzoXN+P1+paWluTN79uyR4zjKyMhwZ4YNGybHcdyZizU1Nam+vj5sAwAAkSuqrQ8oLi7WgQMHtH///lbHgsGgJCkpKSlsf1JSkj799FN3JiYmJuwOzoWZC48PBoNKTExs9f0TExPdmYsVFhbq97//fVuXAwAALNWmOzHV1dWaO3euNm3apO7du3/jnMfjCfvaGNNq38UunrnU/OW+z5IlSxQKhdyturr6sj8PAADYrU0RU1FRodraWqWnpysqKkpRUVEqKyvT008/raioKPcOzMV3S2pra91jPp9Pzc3Nqquru+zMyZMnW/38U6dOtbrLc4HX61V8fHzYBgAAIlebImbUqFGqrKxUIBBwt6FDh+qBBx5QIBDQLbfcIp/Pp9LSUvcxzc3NKisr0/DhwyVJ6enpio6ODpupqanRoUOH3JnMzEyFQiHt27fPndm7d69CoZA7AwAAurY2PScmLi5OaWlpYftiY2PVu3dvd39+fr4KCgqUmpqq1NRUFRQUqGfPnsrNzZUkOY6jKVOmaP78+erdu7cSEhK0YMECDRo0yH2i8IABAzRu3DhNnTpVa9eulSRNmzZN2dnZ6t+//xUvGgAA2K/NT+z9NgsXLlRjY6NmzJihuro6ZWRkaPv27YqLi3NnVq9eraioKE2cOFGNjY0aNWqUNm7cqG7durkzmzdv1pw5c9xXMeXk5KioqOhqny4AALCUxxhjOvokroX6+no5jqNQKHRNnh/Tb/FrbZr/ZMWEq34OAABEmrb8+81nJwEAACsRMQAAwEpEDAAAsBIRAwAArETEAAAAKxExAADASkQMAACwEhEDAACsRMQAAAArETEAAMBKRAwAALASEQMAAKxExAAAACsRMQAAwEpEDAAAsBIRAwAArETEAAAAKxExAADASkQMAACwEhEDAACsRMQAAAArETEAAMBKRAwAALASEQMAAKxExAAAACsRMQAAwEpEDAAAsBIRAwAArETEAAAAKxExAADASkQMAACwEhEDAACsRMQAAAArETEAAMBKRAwAALASEQMAAKxExAAAACsRMQAAwEpEDAAAsBIRAwAArETEAAAAKxExAADASkQMAACwEhEDAACsRMQAAAArETEAAMBKRAwAALASEQMAAKxExAAAACsRMQAAwEpEDAAAsBIRAwAArETEAAAAKxExAADASkQMAACwEhEDAACsRMQAAAArETEAAMBKRAwAALASEQMAAKxExAAAACsRMQAAwEpEDAAAsBIRAwAArETEAAAAKxExAADASkQMAACwEhEDAACsRMQAAAArETEAAMBKbYqYNWvWaPDgwYqPj1d8fLwyMzP1+uuvu8eNMVq+fLn8fr969OihkSNH6vDhw2Hfo6mpSbNnz1afPn0UGxurnJwcHT9+PGymrq5OeXl5chxHjuMoLy9Pp0+fbv8qAQBAxGlTxPTt21crVqzQu+++q3fffVf33HOPfvnLX7qhsnLlSq1atUpFRUXav3+/fD6fxowZo4aGBvd75Ofna8uWLSouLtauXbt05swZZWdnq6WlxZ3Jzc1VIBBQSUmJSkpKFAgElJeXd5WWDAAAIoHHGGOu5BskJCToj3/8ox566CH5/X7l5+dr0aJFkv531yUpKUlPPPGEHnnkEYVCId1www168cUXNWnSJEnSiRMnlJycrG3btmns2LE6cuSIBg4cqPLycmVkZEiSysvLlZmZqQ8++ED9+/f/TudVX18vx3EUCoUUHx9/JUu8pH6LX2vT/CcrJlz1cwAAINK05d/vdj8npqWlRcXFxTp79qwyMzNVVVWlYDCorKwsd8br9WrEiBHavXu3JKmiokLnzp0Lm/H7/UpLS3Nn9uzZI8dx3ICRpGHDhslxHHfmUpqamlRfXx+2AQCAyNXmiKmsrNT1118vr9er6dOna8uWLRo4cKCCwaAkKSkpKWw+KSnJPRYMBhUTE6NevXpddiYxMbHVz01MTHRnLqWwsNB9Do3jOEpOTm7r0gAAgEXaHDH9+/dXIBBQeXm5Hn30UU2ePFnvv/++e9zj8YTNG2Na7bvYxTOXmv+277NkyRKFQiF3q66u/q5LAgAAFmpzxMTExOi2227T0KFDVVhYqCFDhuipp56Sz+eTpFZ3S2pra927Mz6fT83Nzaqrq7vszMmTJ1v93FOnTrW6y/N1Xq/XfdXUhQ0AAESuK36fGGOMmpqalJKSIp/Pp9LSUvdYc3OzysrKNHz4cElSenq6oqOjw2Zqamp06NAhdyYzM1OhUEj79u1zZ/bu3atQKOTOAAAARLVl+LHHHtP48eOVnJyshoYGFRcX66233lJJSYk8Ho/y8/NVUFCg1NRUpaamqqCgQD179lRubq4kyXEcTZkyRfPnz1fv3r2VkJCgBQsWaNCgQRo9erQkacCAARo3bpymTp2qtWvXSpKmTZum7Ozs7/zKJAAAEPnaFDEnT55UXl6eampq5DiOBg8erJKSEo0ZM0aStHDhQjU2NmrGjBmqq6tTRkaGtm/frri4OPd7rF69WlFRUZo4caIaGxs1atQobdy4Ud26dXNnNm/erDlz5rivYsrJyVFRUdHVWC8AAIgQV/w+MZ0V7xMDAIB9vpf3iQEAAOhIRAwAALASEQMAAKxExAAAACsRMQAAwEpEDAAAsBIRAwAArETEAAAAKxExAADASkQMAACwEhEDAACsRMQAAAArETEAAMBKRAwAALASEQMAAKxExAAAACsRMQAAwEpEDAAAsBIRAwAArETEAAAAKxExAADASkQMAACwEhEDAACsRMQAAAArETEAAMBKRAwAALASEQMAAKxExAAAACsRMQAAwEpEDAAAsBIRAwAArETEAAAAKxExAADASkQMAACwEhEDAACsRMQAAAArETEAAMBKRAwAALASEQMAAKxExAAAACsRMQAAwEpEDAAAsBIRAwAArETEAAAAKxExAADASkQMAACwEhEDAACsRMQAAAArETEAAMBKRAwAALASEQMAAKxExAAAACsRMQAAwEpEDAAAsBIRAwAArETEAAAAKxExAADASkQMAACwEhEDAACsRMQAAAArETEAAMBKRAwAALASEQMAAKxExAAAACsRMQAAwEpEDAAAsBIRAwAArETEAAAAKxExAADASkQMAACwUpsiprCwUHfeeafi4uKUmJio++67T0ePHg2bMcZo+fLl8vv96tGjh0aOHKnDhw+HzTQ1NWn27Nnq06ePYmNjlZOTo+PHj4fN1NXVKS8vT47jyHEc5eXl6fTp0+1bJQAAiDhtipiysjLNnDlT5eXlKi0t1VdffaWsrCydPXvWnVm5cqVWrVqloqIi7d+/Xz6fT2PGjFFDQ4M7k5+fry1btqi4uFi7du3SmTNnlJ2drZaWFncmNzdXgUBAJSUlKikpUSAQUF5e3lVYMgAAiAQeY4xp74NPnTqlxMRElZWV6Re/+IWMMfL7/crPz9eiRYsk/e+uS1JSkp544gk98sgjCoVCuuGGG/Tiiy9q0qRJkqQTJ04oOTlZ27Zt09ixY3XkyBENHDhQ5eXlysjIkCSVl5crMzNTH3zwgfr37/+t51ZfXy/HcRQKhRQfH9/eJX6jfotfa9P8JysmXPVzAAAg0rTl3+8rek5MKBSSJCUkJEiSqqqqFAwGlZWV5c54vV6NGDFCu3fvliRVVFTo3LlzYTN+v19paWnuzJ49e+Q4jhswkjRs2DA5juPOXKypqUn19fVhGwAAiFztjhhjjObNm6e77rpLaWlpkqRgMChJSkpKCptNSkpyjwWDQcXExKhXr16XnUlMTGz1MxMTE92ZixUWFrrPn3EcR8nJye1dGgAAsEC7I2bWrFl677339Ne//rXVMY/HE/a1MabVvotdPHOp+ct9nyVLligUCrlbdXX1d1kGAACwVLsiZvbs2dq6dat27typvn37uvt9Pp8ktbpbUltb696d8fl8am5uVl1d3WVnTp482ernnjp1qtVdngu8Xq/i4+PDNgAAELnaFDHGGM2aNUuvvPKK3nzzTaWkpIQdT0lJkc/nU2lpqbuvublZZWVlGj58uCQpPT1d0dHRYTM1NTU6dOiQO5OZmalQKKR9+/a5M3v37lUoFHJnAABA1xbVluGZM2fqpZde0t///nfFxcW5d1wcx1GPHj3k8XiUn5+vgoICpaamKjU1VQUFBerZs6dyc3Pd2SlTpmj+/Pnq3bu3EhIStGDBAg0aNEijR4+WJA0YMEDjxo3T1KlTtXbtWknStGnTlJ2d/Z1emQQAACJfmyJmzZo1kqSRI0eG7d+wYYMefPBBSdLChQvV2NioGTNmqK6uThkZGdq+fbvi4uLc+dWrVysqKkoTJ05UY2OjRo0apY0bN6pbt27uzObNmzVnzhz3VUw5OTkqKipqzxoBAEAEuqL3ienMeJ8YAADs8729TwwAAEBHIWIAAICViBgAAGAlIgYAAFiJiAEAAFYiYgAAgJWIGAAAYCUiBgAAWImIAQAAViJiAACAlYgYAABgJSIGAABYiYgBAABWImIAAICViBgAAGAlIgYAAFgpqqNPoCvqt/i1Ns1/smLCNToTAADsxZ0YAABgJSIGAABYiYgBAABWImIAAICViBgAAGAlIgYAAFiJiAEAAFYiYgAAgJWIGAAAYCUiBgAAWImIAQAAViJiAACAlYgYAABgJSIGAABYiYgBAABWImIAAICViBgAAGAlIgYAAFiJiAEAAFYiYgAAgJWIGAAAYCUiBgAAWImIAQAAViJiAACAlYgYAABgJSIGAABYiYgBAABWImIAAICViBgAAGAlIgYAAFiJiAEAAFYiYgAAgJWIGAAAYCUiBgAAWImIAQAAViJiAACAlYgYAABgJSIGAABYiYgBAABWImIAAICViBgAAGAlIgYAAFiJiAEAAFYiYgAAgJWiOvoE0Db9Fr/WpvlPVky4RmcCAEDH4k4MAACwEhEDAACsRMQAAAArETEAAMBKRAwAALASEQMAAKxExAAAACsRMQAAwEptjpi3335b9957r/x+vzwej1599dWw48YYLV++XH6/Xz169NDIkSN1+PDhsJmmpibNnj1bffr0UWxsrHJycnT8+PGwmbq6OuXl5clxHDmOo7y8PJ0+fbrNCwQAAJGpzRFz9uxZDRkyREVFRZc8vnLlSq1atUpFRUXav3+/fD6fxowZo4aGBncmPz9fW7ZsUXFxsXbt2qUzZ84oOztbLS0t7kxubq4CgYBKSkpUUlKiQCCgvLy8diwRAABEojZ/7MD48eM1fvz4Sx4zxujJJ5/U0qVLdf/990uSnn/+eSUlJemll17SI488olAopPXr1+vFF1/U6NGjJUmbNm1ScnKyduzYobFjx+rIkSMqKSlReXm5MjIyJEnr1q1TZmamjh49qv79+7d3vQAAIEJc1efEVFVVKRgMKisry93n9Xo1YsQI7d69W5JUUVGhc+fOhc34/X6lpaW5M3v27JHjOG7ASNKwYcPkOI47c7GmpibV19eHbQAAIHJd1YgJBoOSpKSkpLD9SUlJ7rFgMKiYmBj16tXrsjOJiYmtvn9iYqI7c7HCwkL3+TOO4yg5OfmK1wMAADqva/LqJI/HE/a1MabVvotdPHOp+ct9nyVLligUCrlbdXV1O84cAADY4qpGjM/nk6RWd0tqa2vduzM+n0/Nzc2qq6u77MzJkydbff9Tp061ustzgdfrVXx8fNgGAAAi11WNmJSUFPl8PpWWlrr7mpubVVZWpuHDh0uS0tPTFR0dHTZTU1OjQ4cOuTOZmZkKhULat2+fO7N3716FQiF3BgAAdG1tfnXSmTNn9NFHH7lfV1VVKRAIKCEhQTfddJPy8/NVUFCg1NRUpaamqqCgQD179lRubq4kyXEcTZkyRfPnz1fv3r2VkJCgBQsWaNCgQe6rlQYMGKBx48Zp6tSpWrt2rSRp2rRpys7O5pVJAABAUjsi5t1339Xdd9/tfj1v3jxJ0uTJk7Vx40YtXLhQjY2NmjFjhurq6pSRkaHt27crLi7Ofczq1asVFRWliRMnqrGxUaNGjdLGjRvVrVs3d2bz5s2aM2eO+yqmnJycb3xvGgAA0PV4jDGmo0/iWqivr5fjOAqFQtfk+TH9Fr/WpvlPVkzo8McCANDZteXfbz47CQAAWImIAQAAViJiAACAlYgYAABgJSIGAABYiYgBAABWImIAAICViBgAAGAlIgYAAFiJiAEAAFYiYgAAgJWIGAAAYCUiBgAAWImIAQAAViJiAACAlYgYAABgJSIGAABYiYgBAABWiuroE8D3p9/i19o0/8mKCdfoTAAAuHLciQEAAFYiYgAAgJWIGAAAYCUiBgAAWImIAQAAViJiAACAlYgYAABgJSIGAABYiYgBAABWImIAAICViBgAAGAlPjsJ3wmfuwQA6Gy4EwMAAKxExAAAACsRMQAAwEpEDAAAsBIRAwAArETEAAAAKxExAADASkQMAACwEhEDAACsRMQAAAArETEAAMBKRAwAALASHwCJa66tHx4p8QGSAIBvx50YAABgJSIGAABYiYgBAABWImIAAICViBgAAGAlIgYAAFiJl1ij02vrS7R5eTYAdA3ciQEAAFYiYgAAgJWIGAAAYCUiBgAAWImIAQAAViJiAACAlYgYAABgJSIGAABYiTe7Q0TjjfIAIHJxJwYAAFiJiAEAAFbiz0nAN+BPUQDQuXEnBgAAWImIAQAAViJiAACAlXhODHANXMnzaXguDgB8N9yJAQAAViJiAACAlfhzEhBBOurPWG197MWPB4D26PR3Yp555hmlpKSoe/fuSk9P1zvvvNPRpwQAADqBTh0xf/vb35Sfn6+lS5fq4MGD+r//+z+NHz9ex44d6+hTAwAAHaxTR8yqVas0ZcoUPfzwwxowYICefPJJJScna82aNR19agAAoIN12ufENDc3q6KiQosXLw7bn5WVpd27d7eab2pqUlNTk/t1KBSSJNXX11+T8zvf9GWb5r9+Hjy28/5sHvv9PPbix6cte6NNjz30+7E8th2PBWxw4f8NxphvHzad1GeffWYkmX/9619h+x9//HHzox/9qNX8smXLjCQ2NjY2Nja2CNiqq6u/tRU67Z2YCzweT9jXxphW+yRpyZIlmjdvnvv1+fPn9cUXX6h3796XnL/a6uvrlZycrOrqasXHx1/zn9cZdfVr0NXXL3ENJK5BV1+/xDWQruwaGGPU0NAgv9//rbOdNmL69Omjbt26KRgMhu2vra1VUlJSq3mv1yuv1xu27wc/+MG1PMVLio+P77K/tBd09WvQ1dcvcQ0krkFXX7/ENZDafw0cx/lOc532ib0xMTFKT09XaWlp2P7S0lINHz68g84KAAB0Fp32TowkzZs3T3l5eRo6dKgyMzP17LPP6tixY5o+fXpHnxoAAOhgnTpiJk2apM8//1x/+MMfVFNTo7S0NG3btk0333xzR59aK16vV8uWLWv1J62upKtfg66+folrIHENuvr6Ja6B9P1dA48x3+U1TAAAAJ1Lp31ODAAAwOUQMQAAwEpEDAAAsBIRAwAArETEXCXPPPOMUlJS1L17d6Wnp+udd97p6FO6Kt5++23de++98vv98ng8evXVV8OOG2O0fPly+f1+9ejRQyNHjtThw4fDZpqamjR79mz16dNHsbGxysnJ0fHjx7/HVbRfYWGh7rzzTsXFxSkxMVH33Xefjh49GjYT6ddgzZo1Gjx4sPumVZmZmXr99dfd45G+/osVFhbK4/EoPz/f3Rfp12D58uXyeDxhm8/nc49H+vov+Oyzz/Tb3/5WvXv3Vs+ePfWTn/xEFRUV7vFIvw79+vVr9Xvg8Xg0c+ZMSR20/iv7hCMYY0xxcbGJjo4269atM++//76ZO3euiY2NNZ9++mlHn9oV27Ztm1m6dKl5+eWXjSSzZcuWsOMrVqwwcXFx5uWXXzaVlZVm0qRJ5sYbbzT19fXuzPTp080Pf/hDU1paag4cOGDuvvtuM2TIEPPVV199z6tpu7Fjx5oNGzaYQ4cOmUAgYCZMmGBuuukmc+bMGXcm0q/B1q1bzWuvvWaOHj1qjh49ah577DETHR1tDh06ZIyJ/PV/3b59+0y/fv3M4MGDzdy5c939kX4Nli1bZn784x+bmpoad6utrXWPR/r6jTHmiy++MDfffLN58MEHzd69e01VVZXZsWOH+eijj9yZSL8OtbW1Yb8DpaWlRpLZuXOnMaZj1k/EXAU/+9nPzPTp08P23X777Wbx4sUddEbXxsURc/78eePz+cyKFSvcff/973+N4zjmL3/5izHGmNOnT5vo6GhTXFzsznz22WfmuuuuMyUlJd/buV8ttbW1RpIpKyszxnTNa2CMMb169TLPPfdcl1p/Q0ODSU1NNaWlpWbEiBFuxHSFa7Bs2TIzZMiQSx7rCus3xphFixaZu+666xuPd5Xr8HVz5841t956qzl//nyHrZ8/J12h5uZmVVRUKCsrK2x/VlaWdu/e3UFn9f2oqqpSMBgMW7vX69WIESPctVdUVOjcuXNhM36/X2lpaVZen1AoJElKSEiQ1PWuQUtLi4qLi3X27FllZmZ2qfXPnDlTEyZM0OjRo8P2d5Vr8OGHH8rv9yslJUW//vWv9fHHH0vqOuvfunWrhg4dql/96ldKTEzUHXfcoXXr1rnHu8p1uKC5uVmbNm3SQw89JI/H02HrJ2Ku0H/+8x+1tLS0+lDKpKSkVh9eGWkurO9yaw8Gg4qJiVGvXr2+ccYWxhjNmzdPd911l9LS0iR1nWtQWVmp66+/Xl6vV9OnT9eWLVs0cODALrP+4uJiHThwQIWFha2OdYVrkJGRoRdeeEFvvPGG1q1bp2AwqOHDh+vzzz/vEuuXpI8//lhr1qxRamqq3njjDU2fPl1z5szRCy+8IKlr/B583auvvqrTp0/rwQcflNRx6+/UHztgE4/HE/a1MabVvkjVnrXbeH1mzZql9957T7t27Wp1LNKvQf/+/RUIBHT69Gm9/PLLmjx5ssrKytzjkbz+6upqzZ07V9u3b1f37t2/cS6Sr8H48ePd/x40aJAyMzN166236vnnn9ewYcMkRfb6Jen8+fMaOnSoCgoKJEl33HGHDh8+rDVr1uh3v/udOxfp1+GC9evXa/z48fL7/WH7v+/1cyfmCvXp00fdunVrVZG1tbWtijTSXHh1wuXW7vP51NzcrLq6um+cscHs2bO1detW7dy5U3379nX3d5VrEBMTo9tuu01Dhw5VYWGhhgwZoqeeeqpLrL+iokK1tbVKT09XVFSUoqKiVFZWpqefflpRUVHuGiL5GlwsNjZWgwYN0ocfftglfgck6cYbb9TAgQPD9g0YMEDHjh2T1HX+XyBJn376qXbs2KGHH37Y3ddR6ydirlBMTIzS09NVWloatr+0tFTDhw/voLP6fqSkpMjn84Wtvbm5WWVlZe7a09PTFR0dHTZTU1OjQ4cOWXF9jDGaNWuWXnnlFb355ptKSUkJO94VrsGlGGPU1NTUJdY/atQoVVZWKhAIuNvQoUP1wAMPKBAI6JZbbon4a3CxpqYmHTlyRDfeeGOX+B2QpJ///Oet3l7h3//+t/uBxF3lOkjShg0blJiYqAkTJrj7Omz97Xo6MMJceIn1+vXrzfvvv2/y8/NNbGys+eSTTzr61K5YQ0ODOXjwoDl48KCRZFatWmUOHjzovnx8xYoVxnEc88orr5jKykrzm9/85pIvqevbt6/ZsWOHOXDggLnnnnuseUnho48+ahzHMW+99VbYSwu//PJLdybSr8GSJUvM22+/baqqqsx7771nHnvsMXPdddeZ7du3G2Mif/2X8vVXJxkT+ddg/vz55q233jIff/yxKS8vN9nZ2SYuLs79f1ykr9+Y/728Pioqyjz++OPmww8/NJs3bzY9e/Y0mzZtcme6wnVoaWkxN910k1m0aFGrYx2xfiLmKvnzn/9sbr75ZhMTE2N++tOfui/Btd3OnTuNpFbb5MmTjTH/e1nhsmXLjM/nM16v1/ziF78wlZWVYd+jsbHRzJo1yyQkJJgePXqY7Oxsc+zYsQ5YTdtdau2SzIYNG9yZSL8GDz30kPu7fcMNN5hRo0a5AWNM5K//Ui6OmEi/Bhfe7yM6Otr4/X5z//33m8OHD7vHI339F/zjH/8waWlpxuv1mttvv908++yzYce7wnV44403jCRz9OjRVsc6Yv0eY4xp3z0cAACAjsNzYgAAgJWIGAAAYCUiBgAAWImIAQAAViJiAACAlYgYAABgJSIGAABYiYgBAABWImIAAICViBgAAGAlIgYAAFiJiAEAAFb6f0wZHsSVVhjeAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "report(sizes)" + ] + }, + { + "cell_type": "markdown", + "id": "7d2c2f1a-5b51-4d02-a35c-05d6f90e5a09", + "metadata": {}, + "source": [ + "There's quite a spread, with the standard deviation being larger than the mean, and the maximum being about 10 times the mean. It might help to zoom in on the smaller sizes:" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "2f1dacba-d47c-46ab-be44-b5129e317ef4", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'mean': 11.6556, 'max': 32, 'min': 1, 'stdev': 9.159878236846971, 'N': 5000}" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAdkElEQVR4nO3df3TV9X348VckEIGFlB+SS2aK6UbXH6HOBsegtnAKxDqp9nDOsNXt0FPWg0PZMnAWxjap5xgom8gmLT06j1IZxT8mm+doV+LUKGWeIcWJtMf2HFFiJSdnLkuCskTh8/3D49338sNwEbnvJI/HOZ9zms99X/K+n77bPM/7fnJTlmVZFgAACbmg1BMAADiRQAEAkiNQAIDkCBQAIDkCBQBIjkABAJIjUACA5AgUACA55aWewNk4fvx4vP7661FZWRllZWWlng4AcAayLIuenp6oqamJCy54/z2SARkor7/+etTW1pZ6GgDAWWhra4uLL774fccMyECprKyMiHdf4JgxY0o8GwDgTHR3d0dtbW3+5/j7GZCB8t7bOmPGjBEoADDAnMntGW6SBQCSI1AAgOQIFAAgOQIFAEiOQAEAkiNQAIDkCBQAIDkCBQBIjkABAJIjUACA5AgUACA5AgUASI5AAQCSI1AAgOSUl3oCKbpk5aNFjX9l3dUf0kwAYGiygwIAJEegAADJESgAQHIECgCQHIECACRHoAAAyREoAEByBAoAkByBAgAkR6AAAMkRKABAcgQKAJAcgQIAJEegAADJESgAQHIECgCQHIECACRHoAAAyREoAEByBAoAkByBAgAkR6AAAMkRKABAcgQKAJAcgQIAJEegAADJESgAQHIECgCQHIECACRHoAAAyREoAEByBAoAkByBAgAkR6AAAMkRKABAcgQKAJAcgQIAJEegAADJESgAQHIECgCQHIECACRHoAAAyREoAEByBAoAkByBAgAkR6AAAMkRKABAcgQKAJAcgQIAJEegAADJKSpQ3nnnnfjLv/zLqKuri5EjR8bHPvaxuP322+P48eP5MVmWxZo1a6KmpiZGjhwZs2fPjgMHDhT8O729vbFs2bKYMGFCjB49Oq655pp47bXXzs0rAgAGvKIC5Tvf+U58//vfj02bNsXPf/7zWL9+ffzN3/xN3H333fkx69evjw0bNsSmTZtiz549kcvlYt68edHT05Mf09TUFDt27Ijt27fHrl274siRIzF//vw4duzYuXtlAMCAVV7M4H//93+Pa6+9Nq6++uqIiLjkkkvihz/8YTz33HMR8e7uycaNG2P16tWxYMGCiIjYsmVLVFdXx7Zt22LJkiXR1dUV9913Xzz44IMxd+7ciIjYunVr1NbWxuOPPx5XXnnluXx9AMAAVNQOyhVXXBH/9m//Fr/4xS8iIuI///M/Y9euXfF7v/d7ERFx8ODBaG9vj8bGxvxzKioqYtasWbF79+6IiNi7d2+8/fbbBWNqamqivr4+P+ZEvb290d3dXXAAAINXUTso3/rWt6Krqys+8YlPxLBhw+LYsWNxxx13xNe+9rWIiGhvb4+IiOrq6oLnVVdXx6uvvpofM2LEiBg7duxJY957/onWrl0b3/72t4uZKgAwgBW1g/LQQw/F1q1bY9u2bfHTn/40tmzZEn/7t38bW7ZsKRhXVlZW8HWWZSedO9H7jVm1alV0dXXlj7a2tmKmDQAMMEXtoPz5n/95rFy5Mr761a9GRMTUqVPj1VdfjbVr18aiRYsil8tFxLu7JJMmTco/r6OjI7+rksvloq+vLzo7Owt2UTo6OmLmzJmn/L4VFRVRUVFR3CsDAAasonZQ3nrrrbjggsKnDBs2LP9rxnV1dZHL5aKlpSX/eF9fX7S2tubjo6GhIYYPH14w5vDhw/Hiiy+eNlAAgKGlqB2UL3/5y3HHHXfERz/60fj0pz8d+/btiw0bNsQ3vvGNiHj3rZ2mpqZobm6OKVOmxJQpU6K5uTlGjRoV119/fUREVFVVxeLFi2PFihUxfvz4GDduXNxyyy0xderU/G/1AABDW1GBcvfdd8df/dVfxdKlS6OjoyNqampiyZIl8dd//df5MbfeemscPXo0li5dGp2dnTF9+vTYuXNnVFZW5sfcddddUV5eHgsXLoyjR4/GnDlz4oEHHohhw4adu1cGAAxYZVmWZaWeRLG6u7ujqqoqurq6YsyYMef8379k5aNFjX9l3dXnfA4AMNgU8/Pb3+IBAJIjUACA5AgUACA5AgUASI5AAQCSI1AAgOQIFAAgOQIFAEiOQAEAkiNQAIDkCBQAIDkCBQBIjkABAJIjUACA5AgUACA5AgUASI5AAQCSI1AAgOQIFAAgOQIFAEiOQAEAkiNQAIDkCBQAIDkCBQBIjkABAJIjUACA5AgUACA5AgUASI5AAQCSI1AAgOQIFAAgOQIFAEiOQAEAkiNQAIDkCBQAIDkCBQBIjkABAJIjUACA5AgUACA5AgUASI5AAQCSI1AAgOQIFAAgOQIFAEiOQAEAkiNQAIDkCBQAIDkCBQBIjkABAJIjUACA5AgUACA5AgUASI5AAQCSI1AAgOQIFAAgOQIFAEiOQAEAkiNQAIDkCBQAIDkCBQBIjkABAJIjUACA5AgUACA5RQfKr371q/iDP/iDGD9+fIwaNSp++7d/O/bu3Zt/PMuyWLNmTdTU1MTIkSNj9uzZceDAgYJ/o7e3N5YtWxYTJkyI0aNHxzXXXBOvvfbaB381AMCgUFSgdHZ2xuc+97kYPnx4/OhHP4qf/exnceedd8ZHPvKR/Jj169fHhg0bYtOmTbFnz57I5XIxb9686OnpyY9pamqKHTt2xPbt22PXrl1x5MiRmD9/fhw7duycvTAAYOAqy7IsO9PBK1eujJ/85CfxzDPPnPLxLMuipqYmmpqa4lvf+lZEvLtbUl1dHd/5zndiyZIl0dXVFRdddFE8+OCDcd1110VExOuvvx61tbXx2GOPxZVXXtnvPLq7u6Oqqiq6urpizJgxZzr9M3bJykeLGv/KuqvP+RwAYLAp5ud3UTsojzzySEybNi1+//d/PyZOnBiXXXZZ3HvvvfnHDx48GO3t7dHY2Jg/V1FREbNmzYrdu3dHRMTevXvj7bffLhhTU1MT9fX1+TEn6u3tje7u7oIDABi8igqUl19+OTZv3hxTpkyJH//4x3HjjTfGn/zJn8QPfvCDiIhob2+PiIjq6uqC51VXV+cfa29vjxEjRsTYsWNPO+ZEa9eujaqqqvxRW1tbzLQBgAGmqEA5fvx4fPazn43m5ua47LLLYsmSJfHNb34zNm/eXDCurKys4Ossy046d6L3G7Nq1aro6urKH21tbcVMGwAYYIoKlEmTJsWnPvWpgnOf/OQn49ChQxERkcvlIiJO2gnp6OjI76rkcrno6+uLzs7O0445UUVFRYwZM6bgAAAGr6IC5XOf+1y89NJLBed+8YtfxOTJkyMioq6uLnK5XLS0tOQf7+vri9bW1pg5c2ZERDQ0NMTw4cMLxhw+fDhefPHF/BgAYGgrL2bwn/3Zn8XMmTOjubk5Fi5cGP/xH/8R99xzT9xzzz0R8e5bO01NTdHc3BxTpkyJKVOmRHNzc4waNSquv/76iIioqqqKxYsXx4oVK2L8+PExbty4uOWWW2Lq1Kkxd+7cc/8KAYABp6hAufzyy2PHjh2xatWquP3226Ouri42btwYN9xwQ37MrbfeGkePHo2lS5dGZ2dnTJ8+PXbu3BmVlZX5MXfddVeUl5fHwoUL4+jRozFnzpx44IEHYtiwYefulQEAA1ZRn4OSCp+DAgADz4f2OSgAAOeDQAEAkiNQAIDkCBQAIDkCBQBIjkABAJIjUACA5AgUACA5AgUASI5AAQCSI1AAgOQIFAAgOQIFAEiOQAEAkiNQAIDkCBQAIDkCBQBIjkABAJIjUACA5AgUACA5AgUASI5AAQCSI1AAgOQIFAAgOQIFAEiOQAEAkiNQAIDkCBQAIDkCBQBIjkABAJIjUACA5AgUACA5AgUASI5AAQCSI1AAgOQIFAAgOQIFAEiOQAEAkiNQAIDkCBQAIDkCBQBIjkABAJIjUACA5AgUACA5AgUASI5AAQCSI1AAgOQIFAAgOQIFAEiOQAEAkiNQAIDkCBQAIDkCBQBIjkABAJIjUACA5AgUACA5AgUASI5AAQCSI1AAgOQIFAAgOQIFAEiOQAEAkiNQAIDkCBQAIDkfKFDWrl0bZWVl0dTUlD+XZVmsWbMmampqYuTIkTF79uw4cOBAwfN6e3tj2bJlMWHChBg9enRcc8018dprr32QqQAAg8hZB8qePXvinnvuic985jMF59evXx8bNmyITZs2xZ49eyKXy8W8efOip6cnP6apqSl27NgR27dvj127dsWRI0di/vz5cezYsbN/JQDAoHFWgXLkyJG44YYb4t57742xY8fmz2dZFhs3bozVq1fHggULor6+PrZs2RJvvfVWbNu2LSIiurq64r777os777wz5s6dG5dddlls3bo19u/fH48//vi5eVUAwIB2VoFy0003xdVXXx1z584tOH/w4MFob2+PxsbG/LmKioqYNWtW7N69OyIi9u7dG2+//XbBmJqamqivr8+POVFvb290d3cXHADA4FVe7BO2b98eP/3pT2PPnj0nPdbe3h4REdXV1QXnq6ur49VXX82PGTFiRMHOy3tj3nv+idauXRvf/va3i50qADBAFbWD0tbWFn/6p38aW7dujQsvvPC048rKygq+zrLspHMner8xq1atiq6urvzR1tZWzLQBgAGmqEDZu3dvdHR0RENDQ5SXl0d5eXm0trbG3//930d5eXl+5+TEnZCOjo78Y7lcLvr6+qKzs/O0Y05UUVERY8aMKTgAgMGrqECZM2dO7N+/P55//vn8MW3atLjhhhvi+eefj4997GORy+WipaUl/5y+vr5obW2NmTNnRkREQ0NDDB8+vGDM4cOH48UXX8yPAQCGtqLuQamsrIz6+vqCc6NHj47x48fnzzc1NUVzc3NMmTIlpkyZEs3NzTFq1Ki4/vrrIyKiqqoqFi9eHCtWrIjx48fHuHHj4pZbbompU6eedNMtADA0FX2TbH9uvfXWOHr0aCxdujQ6Oztj+vTpsXPnzqisrMyPueuuu6K8vDwWLlwYR48ejTlz5sQDDzwQw4YNO9fTAQAGoLIsy7JST6JY3d3dUVVVFV1dXR/K/SiXrHy0qPGvrLv6nM8BAAabYn5++1s8AEByBAoAkByBAgAkR6AAAMkRKABAcgQKAJAcgQIAJEegAADJESgAQHIECgCQHIECACRHoAAAyREoAEByBAoAkByBAgAkR6AAAMkRKABAcgQKAJAcgQIAJEegAADJESgAQHIECgCQHIECACRHoAAAyREoAEByBAoAkByBAgAkR6AAAMkRKABAcgQKAJAcgQIAJEegAADJESgAQHIECgCQHIECACRHoAAAyREoAEByyks9Af7PJSsfLWr8K+uu/pBmAgClJVDOMZEBAB+ct3gAgOQIFAAgOQIFAEiOQAEAkiNQAIDkCBQAIDkCBQBIjkABAJIjUACA5AgUACA5AgUASI5AAQCSI1AAgOT4a8aDhL+iDMBgYgcFAEiOQAEAkiNQAIDkCBQAIDkCBQBIjkABAJIjUACA5PgcFIr+DJUIn6MCwIfLDgoAkByBAgAkp6hAWbt2bVx++eVRWVkZEydOjK985Svx0ksvFYzJsizWrFkTNTU1MXLkyJg9e3YcOHCgYExvb28sW7YsJkyYEKNHj45rrrkmXnvttQ/+agCAQaGoQGltbY2bbropnn322WhpaYl33nknGhsb480338yPWb9+fWzYsCE2bdoUe/bsiVwuF/PmzYuenp78mKamptixY0ds3749du3aFUeOHIn58+fHsWPHzt0rAwAGrKJukv3Xf/3Xgq/vv//+mDhxYuzduze+8IUvRJZlsXHjxli9enUsWLAgIiK2bNkS1dXVsW3btliyZEl0dXXFfffdFw8++GDMnTs3IiK2bt0atbW18fjjj8eVV155jl4aADBQfaB7ULq6uiIiYty4cRERcfDgwWhvb4/Gxsb8mIqKipg1a1bs3r07IiL27t0bb7/9dsGYmpqaqK+vz485UW9vb3R3dxccAMDgddaBkmVZLF++PK644oqor6+PiIj29vaIiKiuri4YW11dnX+svb09RowYEWPHjj3tmBOtXbs2qqqq8kdtbe3ZThsAGADOOlBuvvnmeOGFF+KHP/zhSY+VlZUVfJ1l2UnnTvR+Y1atWhVdXV35o62t7WynDQAMAGf1QW3Lli2LRx55JJ5++um4+OKL8+dzuVxEvLtLMmnSpPz5jo6O/K5KLpeLvr6+6OzsLNhF6ejoiJkzZ57y+1VUVERFRcXZTJXzoNgPevMhbwD0p6gdlCzL4uabb46HH344nnjiiairqyt4vK6uLnK5XLS0tOTP9fX1RWtraz4+GhoaYvjw4QVjDh8+HC+++OJpAwUAGFqK2kG56aabYtu2bfEv//IvUVlZmb9npKqqKkaOHBllZWXR1NQUzc3NMWXKlJgyZUo0NzfHqFGj4vrrr8+PXbx4caxYsSLGjx8f48aNi1tuuSWmTp2a/60eAGBoKypQNm/eHBERs2fPLjh///33x9e//vWIiLj11lvj6NGjsXTp0ujs7Izp06fHzp07o7KyMj/+rrvuivLy8li4cGEcPXo05syZEw888EAMGzbsg70aBhxvDwFwKkUFSpZl/Y4pKyuLNWvWxJo1a0475sILL4y777477r777mK+PQAwRPhbPABAcgQKAJAcgQIAJOesPgcFhjI39gJ8+OygAADJESgAQHK8xcOA5a0WgMHLDgoAkByBAgAkx1s8DEneHgJImx0UACA5AgUASI5AAQCSI1AAgOQIFAAgOX6LB+iX33qCgW+g/e9YoMAAMdD+zwXggxAoMAQUGzcRAgcoLYEC55FdEIAz4yZZACA5AgUASI5AAQCS4x4U4ENVqvtu3O8DA5sdFAAgOQIFAEiOQAEAkiNQAIDkuEkWSNZAvNF1IM4ZUiRQABIhborjeg1uAgXgBEPxB99QfM2kzT0oAEByBAoAkBxv8QBQMt5a4nQECgAfiMjgwyBQAAYJoXDmXKv0uQcFAEiOHRQAOE+K3bmJGLq7NwIFAIrg7aHzw1s8AEByBAoAkByBAgAkR6AAAMkRKABAcgQKAJAcgQIAJEegAADJ8UFtADBADKUPibODAgAkR6AAAMkRKABAcgQKAJAcgQIAJEegAADJESgAQHIECgCQHIECACRHoAAAyREoAEByBAoAkByBAgAkR6AAAMkRKABAcgQKAJCckgbK9773vairq4sLL7wwGhoa4plnninldACARJQsUB566KFoamqK1atXx759++Lzn/98XHXVVXHo0KFSTQkASETJAmXDhg2xePHi+KM/+qP45Cc/GRs3boza2trYvHlzqaYEACSivBTftK+vL/bu3RsrV64sON/Y2Bi7d+8+aXxvb2/09vbmv+7q6oqIiO7u7g9lfsd73ypq/P8/j6Hw3FJ+b889P88t5ff23LN7bim/t+d+eM8t5ff+MH7GvvdvZlnW/+CsBH71q19lEZH95Cc/KTh/xx13ZB//+MdPGn/bbbdlEeFwOBwOh2MQHG1tbf22Qkl2UN5TVlZW8HWWZSedi4hYtWpVLF++PP/18ePH47//+79j/Pjxpxx/Ot3d3VFbWxttbW0xZsyYs5/4IOc69c816p9r1D/XqH+u0ZkZKNcpy7Lo6emJmpqafseWJFAmTJgQw4YNi/b29oLzHR0dUV1dfdL4ioqKqKioKDj3kY985Ky//5gxY5L+LzAVrlP/XKP+uUb9c4365xqdmYFwnaqqqs5oXElukh0xYkQ0NDRES0tLwfmWlpaYOXNmKaYEACSkZG/xLF++PP7wD/8wpk2bFjNmzIh77rknDh06FDfeeGOppgQAJKJkgXLdddfFG2+8EbfffnscPnw46uvr47HHHovJkyd/aN+zoqIibrvttpPeLqKQ69Q/16h/rlH/XKP+uUZnZjBep7IsO5Pf9QEAOH/8LR4AIDkCBQBIjkABAJIjUACA5AypQPne974XdXV1ceGFF0ZDQ0M888wzpZ5SMtasWRNlZWUFRy6XK/W0Surpp5+OL3/5y1FTUxNlZWXxz//8zwWPZ1kWa9asiZqamhg5cmTMnj07Dhw4UJrJllB/1+nrX//6SWvrd3/3d0sz2RJYu3ZtXH755VFZWRkTJ06Mr3zlK/HSSy8VjLGWzuw6DfW1tHnz5vjMZz6T/zC2GTNmxI9+9KP844NtHQ2ZQHnooYeiqakpVq9eHfv27YvPf/7zcdVVV8WhQ4dKPbVkfPrTn47Dhw/nj/3795d6SiX15ptvxqWXXhqbNm065ePr16+PDRs2xKZNm2LPnj2Ry+Vi3rx50dPTc55nWlr9XaeIiC996UsFa+uxxx47jzMsrdbW1rjpppvi2WefjZaWlnjnnXeisbEx3nzzzfwYa+nMrlPE0F5LF198caxbty6ee+65eO655+KLX/xiXHvttfkIGXTr6IP/6b+B4Xd+53eyG2+8seDcJz7xiWzlypUlmlFabrvttuzSSy8t9TSSFRHZjh078l8fP348y+Vy2bp16/Ln/vd//zerqqrKvv/975dghmk48TplWZYtWrQou/baa0synxR1dHRkEZG1trZmWWYtnc6J1ynLrKVTGTt2bPYP//APg3IdDYkdlL6+vti7d280NjYWnG9sbIzdu3eXaFbp+eUvfxk1NTVRV1cXX/3qV+Pll18u9ZSSdfDgwWhvby9YUxUVFTFr1ixr6hSeeuqpmDhxYnz84x+Pb37zm9HR0VHqKZVMV1dXRESMGzcuIqyl0znxOr3HWnrXsWPHYvv27fHmm2/GjBkzBuU6GhKB8l//9V9x7Nixk/4QYXV19Ul/sHComj59evzgBz+IH//4x3HvvfdGe3t7zJw5M954441STy1J760ba6p/V111VfzjP/5jPPHEE3HnnXfGnj174otf/GL09vaWemrnXZZlsXz58rjiiiuivr4+IqylUznVdYqwliIi9u/fH7/2a78WFRUVceONN8aOHTviU5/61KBcRyX7qPtSKCsrK/g6y7KTzg1VV111Vf4/T506NWbMmBG/8Ru/EVu2bInly5eXcGZps6b6d9111+X/c319fUybNi0mT54cjz76aCxYsKCEMzv/br755njhhRdi165dJz1mLf2f010naynit37rt+L555+P//mf/4l/+qd/ikWLFkVra2v+8cG0jobEDsqECRNi2LBhJ1VkR0fHSbXJu0aPHh1Tp06NX/7yl6WeSpLe+w0na6p4kyZNismTJw+5tbVs2bJ45JFH4sknn4yLL744f95aKnS663QqQ3EtjRgxIn7zN38zpk2bFmvXro1LL700/u7v/m5QrqMhESgjRoyIhoaGaGlpKTjf0tISM2fOLNGs0tbb2xs///nPY9KkSaWeSpLq6uoil8sVrKm+vr5obW21pvrxxhtvRFtb25BZW1mWxc033xwPP/xwPPHEE1FXV1fwuLX0rv6u06kMtbV0KlmWRW9v7+BcRyW7Pfc82759ezZ8+PDsvvvuy372s59lTU1N2ejRo7NXXnml1FNLwooVK7Knnnoqe/nll7Nnn302mz9/flZZWTmkr09PT0+2b9++bN++fVlEZBs2bMj27duXvfrqq1mWZdm6deuyqqqq7OGHH87279+ffe1rX8smTZqUdXd3l3jm59f7Xaeenp5sxYoV2e7du7ODBw9mTz75ZDZjxozs13/914fMdfrjP/7jrKqqKnvqqaeyw4cP54+33norP8Za6v86WUtZtmrVquzpp5/ODh48mL3wwgvZX/zFX2QXXHBBtnPnzizLBt86GjKBkmVZ9t3vfjebPHlyNmLEiOyzn/1swa+vDXXXXXddNmnSpGz48OFZTU1NtmDBguzAgQOlnlZJPfnkk1lEnHQsWrQoy7J3fz30tttuy3K5XFZRUZF94QtfyPbv31/aSZfA+12nt956K2tsbMwuuuiibPjw4dlHP/rRbNGiRdmhQ4dKPe3z5lTXJiKy+++/Pz/GWur/OllLWfaNb3wj/zPsoosuyubMmZOPkywbfOuoLMuy7Pzt1wAA9G9I3IMCAAwsAgUASI5AAQCSI1AAgOQIFAAgOQIFAEiOQAEAkiNQAIDkCBQAIDkCBQBIjkABAJIjUACA5Pw/9c6vqN/KeusAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "report(sorted(sizes)[:5000])" + ] + }, + { + "cell_type": "markdown", + "id": "0758d4b5-76fa-47e6-85c1-478a772e400e", + "metadata": {}, + "source": [ + "# Tests\n", + "\n", + "Some unit tests to give confidence in the code, and show examples of use:" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "da7a19ef-d47f-4c69-8fe5-cfaf79861dab", + "metadata": {}, + "outputs": [], + "source": [ + "assert grids(1, 3) == [\n", + " {(0, 0): 'R', (0, 1): 'R', (0, 2): 'R'},\n", + " {(0, 0): 'R', (0, 1): 'R', (0, 2): 'B'},\n", + " {(0, 0): 'R', (0, 1): 'B', (0, 2): 'R'},\n", + " {(0, 0): 'R', (0, 1): 'B', (0, 2): 'B'},\n", + " {(0, 0): 'B', (0, 1): 'R', (0, 2): 'R'},\n", + " {(0, 0): 'B', (0, 1): 'R', (0, 2): 'B'},\n", + " {(0, 0): 'B', (0, 1): 'B', (0, 2): 'R'},\n", + " {(0, 0): 'B', (0, 1): 'B', (0, 2): 'B'}]\n", + "\n", + "assert grids(1, 2, 'RGB') == [\n", + " {(0, 0): 'R', (0, 1): 'R'},\n", + " {(0, 0): 'R', (0, 1): 'G'},\n", + " {(0, 0): 'R', (0, 1): 'B'},\n", + " {(0, 0): 'G', (0, 1): 'R'},\n", + " {(0, 0): 'G', (0, 1): 'G'},\n", + " {(0, 0): 'G', (0, 1): 'B'},\n", + " {(0, 0): 'B', (0, 1): 'R'},\n", + " {(0, 0): 'B', (0, 1): 'G'},\n", + " {(0, 0): 'B', (0, 1): 'B'}]\n", + "\n", + "grid6x1 = one_grid(6, 1, 'RRBBBR')\n", + "assert grid6x1 == {(0, 0): 'R', (1, 0): 'R', (2, 0): 'B', (3, 0): 'B', (4, 0): 'B', (5, 0): 'R'}\n", + "assert cluster(grid6x1) == {(0, 0): 1, (1, 0): 1, (2, 0): 2, (3, 0): 2, (4, 0): 2, (5, 0): 3}\n", + "assert mean_cluster_size(grid6x1) == 2\n", + "\n", + "grid5x3 = one_grid(5, 3, 'RR:RR'\n", + " '.R:R.'\n", + " '.RRR.')\n", + "assert cluster(grid5x3) == {\n", + " (0, 0): 1, (1, 0): 1, (2, 0): 2, (3, 0): 1, (4, 0): 1,\n", + " (0, 1): 3, (1, 1): 1, (2, 1): 2, (3, 1): 1, (4, 1): 4,\n", + " (0, 2): 3, (1, 2): 1, (2, 2): 1, (3, 2): 1, (4, 2): 4}\n", + "assert mean_cluster_size(grid5x3) == 3.75\n", + "\n", + "grid10x2 = one_grid(10, 2, 'RBRBRRRBBR' # Example from diagram at top of notebook\n", + " 'RRRRBBRBRB') \n", + "assert mean_cluster_size(grid10x2) == 20/9\n", + "\n", + "assert len(random_grid(3, 4)) == 12" + ] + } + ], + "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.8.15" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/ipynb/Stubborn.ipynb b/ipynb/Stubborn.ipynb new file mode 100644 index 0000000..8405d25 --- /dev/null +++ b/ipynb/Stubborn.ipynb @@ -0,0 +1,596 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "fc5197e8-5f15-4809-8356-8eb70471637a", + "metadata": {}, + "source": [ + "
Peter Norvig
Mar 2024
\n", + "\n", + "# Stubborn Number Endings\n", + "\n", + "[Francis Su](https://www.francissu.com/)'s book *Mathematics for Human Flourishing* mentions the fact that numbers that end in \"5\" have a square that also ends in \"5\". \n", + "\n", + "For example, 5² = 25, 15² = 225, and 25² = 625. \n", + "\n", + "This leads to some questions:\n", + "\n", + "- Is there an easy way to calculate the square of a number ending in \"5\"?\n", + "- What should we call this property of \"square has same ending\"?\n", + "- Are there other digits besides 5 that have this property?\n", + "- Can we prove the property, not just show some examples?\n", + "- Are there multi-digit endings that have this property?\n", + "\n", + "\n", + "## Is there an easy way to calculate the square of a number ending in \"5\"?\n", + "\n", + "Let's make a table of {number: square} pairs and try to see a pattern:" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "fd3c4f25-fb63-43e5-925c-e4b7b8f134ec", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{15: 225,\n", + " 25: 625,\n", + " 35: 1225,\n", + " 45: 2025,\n", + " 55: 3025,\n", + " 65: 4225,\n", + " 75: 5625,\n", + " 85: 7225,\n", + " 95: 9025}" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "{i: i ** 2 for i in range(15, 100, 10)}" + ] + }, + { + "cell_type": "markdown", + "id": "44d2cf09-b397-46d7-ab4f-e068e571895d", + "metadata": {}, + "source": [ + "I see a pattern here: for example 95 squared is 9025, which you get by taking the \"9\", multiplying it by one more than itself to get 9⋅10 = 90, then squaring 5 to get 25, then putting the \"90\" next to the \"25\" to get \"9025\".\n", + "\n", + "Does the trick also work for numbers with more digits? Let's investigate:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "7681cb5f-bf5a-4317-b842-e0be21745d90", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{105: 11025,\n", + " 115: 13225,\n", + " 125: 15625,\n", + " 135: 18225,\n", + " 145: 21025,\n", + " 155: 24025,\n", + " 165: 27225,\n", + " 175: 30625,\n", + " 185: 34225,\n", + " 195: 38025,\n", + " 205: 42025,\n", + " 215: 46225,\n", + " 225: 50625,\n", + " 235: 55225,\n", + " 245: 60025}" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "{i: i ** 2 for i in range(105, 250, 10)}" + ] + }, + { + "cell_type": "markdown", + "id": "89b1e783-0546-4305-989c-b3b83781cf92", + "metadata": {}, + "source": [ + "Yes, it looks like the trick still works: 105 squared is 10⋅11 = \"110\", followed by \"25,\" to make \"11025\". And 245 squared is 24⋅25 = \"600\", followed by \"25\", to make \"60025\". \n", + "\n", + "\n", + "## What should we call this property?\n", + "\n", + "Let's define an **ending** as the rightmost-digits (zero or more) of a number in decimal notation. \n", + "\n", + "Then we can say that an ending is **stubborn** if every number with that ending has a square with the same ending.\n", + "\n", + "## Are there other digits besides 5 that are stubborn?\n", + "\n", + "We could work this out in our heads, or with paper and pencil, or we can compute it with an expression:" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "c86f99ed-48fa-4954-9467-b72d67252d73", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'0', '1', '5', '6'}" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "digits = '0123456789'\n", + "\n", + "{e for e in digits if str(int(e) ** 2).endswith(e)}" + ] + }, + { + "cell_type": "markdown", + "id": "f4cb87cf-2fa8-4d2e-bc24-94c378930792", + "metadata": {}, + "source": [ + "These are the four digits whose square ends in the same digit.\n", + "\n", + "## Can we prove stubborness?\n", + "\n", + "We have seen that 0² ends in 0, 1² ends in 1, 5² ends in 5, and 6² ends in 6. And we have checked some numbers with those endings, for example, 245² ends in 5. But can we **prove** that **every** number ending in 0, 1, 5, or 6 has a square that ends in the same digit?\n", + "\n", + "Some notation: I'll use quote marks, as in \"*se*\" to mean the string of staring digits \"*s*\" followed by the string of ending digits \"*e*\".\n", + "\n", + "With a little bit of algebra we can see that if *s* is any string of digits and *e* is a single ending digit, then:\n", + "\n", + "\"*se*\"² = (10⋅*s* + *e*)² = (10⋅*s*)² + 2⋅(10⋅*s* ⋅ *e*) + *e*² = 10 ⋅ (10⋅*s*² + 2⋅*s*⋅*e*) + *e*²\n", + "\n", + "This is ten times some integer, plus *e*², so \"*se*\"² ends in the digit *e* if and only if *e*² ends in *e*, and we know that is true for 0, 1, 5, and 6, and for no other digits.\n", + "\n", + "## Are there multi-digit endings that are stubborn?\n", + "\n", + "The algebraic argument above can be extended to work with an ending string *e* that is *k* digits long:\n", + "\n", + "\"*se*\"² = (10*k*⋅*s* + *e*)² = (10*k*⋅*s*)² + 2⋅(10*k*⋅*s* ⋅ *e*) + *e*² = 10*k* ⋅ (10*k*⋅*s*² + 2⋅*s*⋅*e*) + *e*²\n", + "\n", + "This is 10*k* times some integer, plus *e*², so again \"*se*\"² ends in *e* if and only if *e*² ends in *e*. To put it another way, to test whether *e* is stubborn, all we have to do is square *e* and see if the result ends in *e*. There is one complication: we'd like to say that \"00\" is stubborn, because any number ending in \"00\", when squared, also ends in \"00\", for example 200² = 40000. But 00² is zero, which we write as \"0\", not as \"00\" or \"0000\". To make sure that \"00\" is considered stubborn, I will define the predicate function `stubborn(ending)` to test the square of \"1\" followed by the ending (I could have chosen any other starting digit string and the result would be the same). \n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "1419b122-257d-475b-90c6-e832c0aae8cf", + "metadata": {}, + "outputs": [], + "source": [ + "def stubborn(ending: str) -> bool:\n", + " \"\"\"Does the square of any number with this ending also end with this ending?\"\"\"\n", + " return str(int(\"1\" + ending) ** 2).endswith(ending)" + ] + }, + { + "cell_type": "markdown", + "id": "9ec3ae19-96cb-4bea-ae6f-075b3ade7d3b", + "metadata": {}, + "source": [ + "Now we can find all two-digit stubborn endings as follows:" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "57e64e31-0674-4849-bf72-f6c6b437009a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['00', '01', '25', '76']" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "[s + e for s in digits for e in digits if stubborn(s + e)]" + ] + }, + { + "cell_type": "markdown", + "id": "61a2b8e2-2bdf-458f-adaf-ffde56125862", + "metadata": {}, + "source": [ + "We could easily continue on in this way, enumerating all three-, four- or even six-digit endings, and checking each one to see if it is stubborn. There are only a million six-digit endings. But there are a quadrillion 15-digit endings, so it would take a *long* time to check all of those. And 100-digit endings? Forget about it. Instead we need to rely on a simplification:\n", + "- Any two-digit stubborn ending ('00', '01', '25', '76') has to end in a one-digit-stubborn ending ('0', '1', '5', '6').\n", + "- In general, any *d*-digit stubborn ending has to end in a (*d*-1)-digit stubborn ending.\n", + "- So, to find the stubborn endings of length 100, I don't need to generate and test all 10100 endings, I only need to consider the stubborn endings of length 99 and check each one of them. (If this simplification is not obvious, stop and convince yourself it is true.)\n", + "\n", + "Using this simplification we can efficiently compute all stubborn-endings of a given length *d* as follows (caching greatly improves efficiency):" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "abeba9cb-1357-4e2b-a661-8616d0d3f768", + "metadata": {}, + "outputs": [], + "source": [ + "from functools import lru_cache\n", + "\n", + "@lru_cache(None)\n", + "def stubborn_endings(d: int) -> list:\n", + " \"\"\"A list of all stubborn endings of length `d` digits.\"\"\"\n", + " if d == 0:\n", + " return [''] # The empty ending is the sole stubborn ending of length 0.\n", + " else:\n", + " return [(s + e) for e in stubborn_endings(d - 1) \n", + " for s in digits \n", + " if stubborn(s + e)]" + ] + }, + { + "cell_type": "markdown", + "id": "05c524e7-6fa8-4fdd-875d-ef5cc0453560", + "metadata": {}, + "source": [ + "For example:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "3943be37-ec8a-4ec0-b946-acb2f449dd72", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['000', '001', '625', '376']" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "stubborn_endings(3)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "d1dc31eb-e5dd-405e-b16d-eceef5c9493c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['0000', '0001', '0625', '9376']" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "stubborn_endings(4)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "81730e6b-7218-4fca-9107-2a9539f62ea3", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['00000', '00001', '90625', '09376']" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "stubborn_endings(5)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "288506c8-6dc1-40d5-842a-eb65ed94ae98", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['000000', '000001', '890625', '109376']" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "stubborn_endings(6)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "73491ea5-a4c8-4221-8fc7-d6ce44d17bb4", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000',\n", + " '0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001',\n", + " '3953007319108169802938509890062166509580863811000557423423230896109004106619977392256259918212890625',\n", + " '6046992680891830197061490109937833490419136188999442576576769103890995893380022607743740081787109376']" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "stubborn_endings(100)" + ] + }, + { + "cell_type": "markdown", + "id": "37a45341-8dac-4720-be7d-388fa47eec1a", + "metadata": {}, + "source": [ + "## More questions!\n", + "\n", + "This leads to a few new questions.\n", + "\n", + "## Can each stubborn ending be extended exactly one way?\n", + "\n", + "We know that each stubborn ending of length *d* has to build on the endings of length *d - 1*, but is it always the case that for any length *d* there will be exactly four endings? Might there be a case where two digits work with one of the endings, or no digits?\n", + "\n", + "I can show that there are always 4 endings up to length *d* = 2000, but I leave it to you to describe a proof that this will always be true." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "654d53af-e58f-4ed2-874b-c4573f3a9a7a", + "metadata": {}, + "outputs": [], + "source": [ + "for d in range(1, 2000):\n", + " assert len(stubborn_endings(d)) == 4 " + ] + }, + { + "cell_type": "markdown", + "id": "537556a3-d570-4a50-93bf-9632b40f9b47", + "metadata": {}, + "source": [ + "One cool implication: if it is true that a stubborn ending of any length can always be extended, that means there is an infinitely long integer whose square ends with itself!\n", + "\n", + "## What digits are used to extend each ending?\n", + "\n", + "There doesn't seem to be a pattern, and all digits seemingly get used roughly evenly. \n", + "\n", + "I don't have a theory of which digit comes next, but I can count how many times each digit appears in the 2000-digit endings that end in \"5\" and \"6\", and see that each of the ten digits appears about 200 times:" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "4dbd7ace-8b41-4ea7-b71f-8f7e318243e1", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Counter({'0': 205,\n", + " '3': 173,\n", + " '2': 208,\n", + " '6': 197,\n", + " '9': 198,\n", + " '5': 197,\n", + " '4': 206,\n", + " '8': 214,\n", + " '7': 205,\n", + " '1': 197})" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from collections import Counter\n", + "\n", + "zero, one, five, six = stubborn_endings(2000)\n", + "\n", + "Counter(five)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "96bdab37-c5a3-4597-b4ca-fb560e4c3163", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Counter({'9': 205,\n", + " '6': 174,\n", + " '7': 208,\n", + " '3': 197,\n", + " '0': 198,\n", + " '4': 196,\n", + " '5': 206,\n", + " '1': 214,\n", + " '2': 205,\n", + " '8': 197})" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Counter(six)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "808e37de-9f5f-4fbf-b55e-7ba59c89e6c6", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 58, + "id": "881e3458-1e14-453e-866c-99e86396daaa", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{(1, 1): 1.0,\n", + " (2, 1): 1.5,\n", + " (3, 1): 1.75,\n", + " (4, 1): 1.875,\n", + " (5, 1): 1.9375,\n", + " (6, 1): 1.96875,\n", + " (7, 1): 1.984375,\n", + " (8, 1): 1.9921875,\n", + " (9, 1): 1.99609375,\n", + " (10, 1): 1.998046875,\n", + " (11, 1): 1.9990234375,\n", + " (12, 1): 1.99951171875,\n", + " (13, 1): 1.999755859375,\n", + " (14, 1): 1.9998779296875}" + ] + }, + "execution_count": 58, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "{(w, 1): mean_mean_size(w, 1) for w in range(1, 15)}" + ] + }, + { + "cell_type": "code", + "execution_count": 59, + "id": "df36e42e-5dd3-44f7-a0e8-b99027408a4a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{(1, 2): 1.5,\n", + " (2, 2): 2.125,\n", + " (3, 2): 2.46875,\n", + " (4, 2): 2.6682291666666664,\n", + " (5, 2): 2.7920386904761907,\n", + " (6, 2): 2.873721168154762,\n", + " (7, 2): 2.9304973563762626,\n", + " (8, 2): 2.971709969311288}" + ] + }, + "execution_count": 59, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "{(w, 2): mean_mean_size(w, 2) for w in range(1, 9)}" + ] + }, + { + "cell_type": "code", + "execution_count": 62, + "id": "47938141-8caa-4304-8575-8c6206c9be96", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{(1, 3): 1.75,\n", + " (2, 3): 2.46875,\n", + " (3, 3): 2.8986049107142855,\n", + " (4, 3): 3.1591145833333334,\n", + " (5, 3): 3.326567058836346,\n", + " (6, 3): 3.4407984463961334,\n", + " (7, 3): 3.522721585152764}" + ] + }, + "execution_count": 62, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "{(w, 3): mean_mean_size(w, 3) for w in range(1, 8)}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "85a991a3-65c9-488b-b50f-799e937493f6", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "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.8.15" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}