PS2b.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ac1c0408-d967-4d03-bb15-cc77e52e6fa0",
   "metadata": {},
   "source": [
    "## Part B\n",
    "### 0: \n",
    "John Evans, ...\n",
    "### 1:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "09d4b9a5-18dd-477b-85a0-5aaa98cb18be",
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import permutations\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib_inline.backend_inline\n",
    "import numpy as np\n",
    "import networkx as nx\n",
    "\n",
    "# uncomment this before exporting to pdf\n",
    "#matplotlib_inline.backend_inline.set_matplotlib_formats('svg', 'pdf')\n",
    "\n",
    "rng = np.random.default_rng()\n",
    "\n",
    "n_a = 200\n",
    "k_a = 2\n",
    "n_b = 50\n",
    "k_b = 8\n",
    "realizations = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "619b1ccf-55d0-4b6c-ad6f-cb143800216c",
   "metadata": {},
   "outputs": [],
   "source": [
    "stubs = np.repeat(np.arange(n_b), k_b)\n",
    "\n",
    "# this allows for double links. unsure if that's okay.\n",
    "a_links = rng.permuted(np.tile(stubs, (realizations, 1)), axis=1)\n",
    "a_links = np.reshape(a_links, (realizations, n_a, k_a))\n",
    "\n",
    "# likely a way to do this without loops\n",
    "a_proj = np.empty((realizations, n_a, n_a), dtype=int)\n",
    "for i in range(realizations): \n",
    "    for j, k in permutations(range(n_a), 2):\n",
    "        a_proj[i, j, k] = np.any(np.isin(a_links[i, j], a_links[i, k]))\n",
    "\n",
    "a_proj_k = np.tensordot(a_proj, np.ones(n_a), axes=([2], [0]))\n",
    "a_proj_c = np.array([\n",
    "    list(nx.clustering(nx.from_numpy_array(proj)).values())\n",
    "    for proj in a_proj\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0064667e-762a-4ec9-bb52-3fe1794cd6df",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAANtUlEQVR4nO3cf6zdd13H8eeLNgM3cQN6h9J23EYK0ixjGzcTJSLJpukGaVXUrGEJxGX9a4qyaIozi46YMGc0/lF/TMGRBTbGFKyu0OGcYoxbesd+uLYU6hhrC2xlzhFdpKu+/eN8p4fbe3vPbc/t994Pz0dy0/P9nk/O952tffbb7znfk6pCkrT8vaTvASRJ42HQJakRBl2SGmHQJakRBl2SGrGyrwOvWrWqJicn+zq8JC1LDz744DeramK253oL+uTkJNPT030dXpKWpSRfnes5L7lIUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiN6u1NUkvo0ue3u3o79xIfesSiv6xm6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSIwy6JDXCoEtSI0YKepKNSfYnOZBk2yzPn5fkviQPJXk0yRXjH1WSdCLzBj3JCmA7cDmwAdiSZMOMZb8B3FlVFwFXAn847kElSSc2yhn6JcCBqnq8qo4CdwCbZ6wp4Pu6x2cDXxvfiJKkUYwS9NXAwaHtQ92+Yb8JXJXkELAT+MXZXijJ1iTTSaaPHDlyEuNKkuYyrjdFtwC3VtUa4ArgtiTHvXZV3VJVU1U1NTExMaZDS5JgtKAfBtYOba/p9g27GrgToKr+GXgZsGocA0qSRjNK0HcD65OsS3IGgzc9d8xY8yRwKUCSNzIIutdUJOk0mjfoVXUMuBbYBexj8GmWPUluTLKpW3YdcE2SR4DbgfdWVS3W0JKk460cZVFV7WTwZufwvhuGHu8F3jre0SRJC+GdopLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUiJGCnmRjkv1JDiTZNsean0+yN8meJB8f75iSpPmsnG9BkhXAduAngEPA7iQ7qmrv0Jr1wAeAt1bVs0nOXayBJUmzG+UM/RLgQFU9XlVHgTuAzTPWXANsr6pnAarq6fGOKUmazyhBXw0cHNo+1O0b9nrg9Un+Kcn9STbO9kJJtiaZTjJ95MiRk5tYkjSrcb0puhJYD7wd2AL8aZJzZi6qqluqaqqqpiYmJsZ0aEkSjBb0w8Daoe013b5hh4AdVfVCVX0F+BKDwEuSTpNRgr4bWJ9kXZIzgCuBHTPWfJrB2TlJVjG4BPP4+MaUJM1n3qBX1THgWmAXsA+4s6r2JLkxyaZu2S7gmSR7gfuAX62qZxZraEnS8eb92CJAVe0Eds7Yd8PQ4wLe3/1IknrgnaKS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1AiDLkmNMOiS1IiRgp5kY5L9SQ4k2XaCde9KUkmmxjeiJGkU8wY9yQpgO3A5sAHYkmTDLOteDrwPeGDcQ0qS5jfKGfolwIGqeryqjgJ3AJtnWfdB4Cbgv8Y4nyRpRKMEfTVwcGj7ULfv/yS5GFhbVXef6IWSbE0ynWT6yJEjCx5WkjS3U35TNMlLgN8DrptvbVXdUlVTVTU1MTFxqoeWJA0ZJeiHgbVD22u6fS96OXA+8PdJngDeAuzwjVFJOr1GCfpuYH2SdUnOAK4Edrz4ZFU9V1WrqmqyqiaB+4FNVTW9KBNLkmY1b9Cr6hhwLbAL2AfcWVV7ktyYZNNiDyhJGs3KURZV1U5g54x9N8yx9u2nPpYkaaG8U1SSGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGmHQJakRBl2SGrGy7wEkfXeb3HZ33yM0wzN0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWqEQZekRhh0SWrESEFPsjHJ/iQHkmyb5fn3J9mb5NEk9yZ57fhHlSSdyLxBT7IC2A5cDmwAtiTZMGPZQ8BUVV0A3AX8zrgHlSSd2Chn6JcAB6rq8ao6CtwBbB5eUFX3VdXz3eb9wJrxjilJms8oQV8NHBzaPtTtm8vVwGdOZShJ0sKN9dsWk1wFTAE/PsfzW4GtAOedd944Dy1J3/VGOUM/DKwd2l7T7fsOSS4Drgc2VdW3Z3uhqrqlqqaqampiYuJk5pUkzWGUoO8G1idZl+QM4Epgx/CCJBcBf8Ig5k+Pf0xJ0nzmDXpVHQOuBXYB+4A7q2pPkhuTbOqW3Qx8L/DJJA8n2THHy0mSFslI19Craiewc8a+G4YeXzbmuSRJC+SdopLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUCIMuSY0w6JLUiJV9D3AyJrfd3duxn/jQO3o7trRY+vwzpfHxDF2SGmHQJakRIwU9ycYk+5McSLJtludfmuQT3fMPJJkc+6SSpBOaN+hJVgDbgcuBDcCWJBtmLLsaeLaqXgf8PnDTuAeVJJ3YKGfolwAHqurxqjoK3AFsnrFmM/DR7vFdwKVJMr4xJUnzGeVTLquBg0Pbh4AfnmtNVR1L8hzwKuCbw4uSbAW2dpv/kWT/yQwNrJr52qdLTvxvj97mmodzLcxSnQuW7mzOtQC56ZTmeu1cT5zWjy1W1S3ALaf6Okmmq2pqDCONlXMtjHMt3FKdzbkWZrHmGuWSy2Fg7dD2mm7frGuSrATOBp4Zx4CSpNGMEvTdwPok65KcAVwJ7JixZgfwnu7xzwJ/V1U1vjElSfOZ95JLd038WmAXsAL4SFXtSXIjMF1VO4APA7clOQD8G4PoL6ZTvmyzSJxrYZxr4ZbqbM61MIsyVzyRlqQ2eKeoJDXCoEtSI5Zd0JOck+SuJF9Msi/JjyyBmd6Q5OGhn28l+eW+5wJI8itJ9iR5LMntSV7W90wASd7XzbSnz/9WST6S5Okkjw3te2WSzyX5cvfrK5bIXD/X/ff6nyS9fBRvjrlu7v48PprkU0nOWSJzfbCb6eEk9yR5zVKYa+i565JUklXjOt6yCzrwB8Bnq+qHgDcB+3qeh6raX1UXVtWFwJuB54FP9TsVJFkN/BIwVVXnM3hTe7HfsJ5XkvOBaxjchfwm4J1JXtfTOLcCG2fs2wbcW1XrgXu77dPtVo6f6zHgZ4DPn/Zp/t+tHD/X54Dzq+oC4EvAB073UMw+181VdUH35/JvgBtO91DMPhdJ1gI/CTw5zoMtq6AnORt4G4NP1VBVR6vq33sd6niXAv9aVV/te5DOSuB7uvsDzgS+1vM8AG8EHqiq56vqGPAPDEJ12lXV5xl8MmvY8FdZfBT4qdM5E8w+V1Xtq6qTvbt6LOaY657u/yPA/QzuVVkKc31raPMs4LR/AmSO318w+M6rX2PMMy2roAPrgCPAnyd5KMmfJTmr76FmuBK4ve8hAKrqMPC7DM4Cvg48V1X39DsVMDjT/LEkr0pyJnAF33nzWt9eXVVf7x5/A3h1n8MsM78AfKbvIV6U5LeTHATeTT9n6MdJshk4XFWPjPu1l1vQVwIXA39UVRcB/0k//xyeVXfj1Sbgk33PAtBd+93M4C/C1wBnJbmq36kGZ5oMvpHzHuCzwMPAf/c501y6G+T8bO8IklwPHAM+1vcsL6qq66tqLYOZru17nu4E5tdZpL9cllvQDwGHquqBbvsuBoFfKi4HvlBVT/U9SOcy4CtVdaSqXgD+EvjRnmcCoKo+XFVvrqq3Ac8yuPa6VDyV5AcAul+f7nmeJS/Je4F3Au9eoneJfwx4V99DAD/I4ATrkSRPMLg89YUk3z+OF19WQa+qbwAHk7yh23UpsLfHkWbawhK53NJ5EnhLkjO7rzO+lCXwJjJAknO7X89jcP384/1O9B2Gv8riPcBf9TjLkpdkI4PrwZuq6vm+53lRkvVDm5uBL/Y1y4uq6l+q6tyqmqyqSQYnqRd3bRvLAZbVD3AhMA08CnwaeEXfM3VzncXgC8nO7nuWGXP9FoPfyI8BtwEv7Xumbq5/ZPCX8SPApT3OcTuD9xde6P5wXc3gq5/vBb4M/C3wyiUy1093j78NPAXsWiJzHWDw9dkPdz9/vETm+ovu9/2jwF8Dq5fCXDOefwJYNa7jeeu/JDViWV1ykSTNzaBLUiMMuiQ1wqBLUiMMuiQ1wqBLUiMMuiQ14n8B/o767Y0oo84AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAN60lEQVR4nO3df6xk9V3G8fcjW1KpWKh7WyuLXTSAYlMDXhVtLBTaZgMNq9UYiChb0U0apbXBErCJGP+RamPVtLFZ6QpWXFIR62qtLaEg0QD28rPLj1KkSBeoOxRbfzQpxX78407j9vbunXPnnDt3v7vvV3KzM+ecO9/nm9n75MyZmXNSVUiS2vMt6x1AkjQdC1ySGmWBS1KjLHBJapQFLkmN2jDLwTZu3FibN2+e5ZCS1Ly77rrrmaqaW7p8pgW+efNmFhYWZjmkJDUvyb8tt9xDKJLUKAtckho1scCT7EyyL8meJcsvSfJwkgeS/O7aRZQkLafLHvg1wJb9FyR5LbAV+MGq+gHg3cNHkyStZGKBV9VtwLNLFr8FuKqqvjLeZt8aZJMkrWDaY+AnAT+R5M4k/5jkhw+0YZLtSRaSLIxGoymHkyQtNW2BbwBeApwOvAP4UJIst2FV7aiq+aqan5v7po8xSpKmNG2B7wVurEX/AnwN2DhcLEnSJNMW+IeB1wIkOQk4EnhmoEySpA4mfhMzyS7gTGBjkr3AlcBOYOf4o4XPARfVGl8ZYvPlH1nLh1/R41edu25jS9KBTCzwqrrgAKsuHDiLJGkV/CamJDXKApekRlngktQoC1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ1ygKXpEZZ4JLUKAtckhplgUtSoyxwSWqUBS5JjbLAJalRFrgkNWpigSfZmWTf+PJpS9ddmqSSeEFjSZqxLnvg1wBbli5McjzwBuCJgTNJkjqYWOBVdRvw7DKr3gNcBqzpxYwlScub6hh4kq3Ak1V1X4dttydZSLIwGo2mGU6StIxVF3iSo4DfAH6zy/ZVtaOq5qtqfm5ubrXDSZIOYJo98O8FTgDuS/I4sAm4O8l3DhlMkrSyDav9har6FPDSr98fl/h8VT0zYC5J0gRdPka4C7gdODnJ3iQXr30sSdIkE/fAq+qCCes3D5ZGktSZ38SUpEZZ4JLUKAtckhplgUtSoyxwSWqUBS5JjbLAJalRFrgkNcoCl6RGWeCS1CgLXJIaZYFLUqMscElqlAUuSY2ywCWpURa4JDXKApekRnW5pNrOJPuS7Nlv2e8leTjJ/Un+Oskxa5pSkvRNuuyBXwNsWbLsJuCVVfUq4BHgioFzSZImmFjgVXUb8OySZR+vqufHd+8ANq1BNknSCoY4Bv6LwEcHeBxJ0ir0KvAk7wSeB65bYZvtSRaSLIxGoz7DSZL2M3WBJ9kGvBH4uaqqA21XVTuqar6q5ufm5qYdTpK0xIZpfinJFuAy4Iyq+vKwkSRJXXT5GOEu4Hbg5CR7k1wMvBc4Grgpyb1J3r/GOSVJS0zcA6+qC5ZZ/IE1yCJJWgW/iSlJjbLAJalRFrgkNcoCl6RGWeCS1CgLXJIaZYFLUqMscElqlAUuSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1KjLHBJapQFLkmNssAlqVFdrom5M8m+JHv2W/aSJDcl+cz432PXNqYkaakue+DXAFuWLLscuLmqTgRuHt+XJM3QxAKvqtuAZ5cs3gpcO759LfCTw8aSJE0y7THwl1XV0+PbnwdedqANk2xPspBkYTQaTTmcJGmp3m9iVlUBtcL6HVU1X1Xzc3NzfYeTJI1NW+D/nuTlAON/9w0XSZLUxbQFvhu4aHz7IuBvhokjSeqqy8cIdwG3Aycn2ZvkYuAq4PVJPgO8bnxfkjRDGyZtUFUXHGDV2QNnkSStgt/ElKRGWeCS1CgLXJIaZYFLUqMscElqlAUuSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ1ygKXpEb1KvAkb0/yQJI9SXYleeFQwSRJK5u6wJMcB7wVmK+qVwJHAOcPFUyStLK+h1A2AN+aZANwFPBU/0iSpC6mLvCqehJ4N/AE8DTwpar6+NLtkmxPspBkYTQaTZ9UkvQN+hxCORbYCpwAfBfwoiQXLt2uqnZU1XxVzc/NzU2fVJL0DfocQnkd8NmqGlXVV4EbgR8fJpYkaZI+Bf4EcHqSo5IEOBt4aJhYkqRJ+hwDvxO4Abgb+NT4sXYMlEuSNMGGPr9cVVcCVw6URZK0Cn4TU5IaZYFLUqMscElqlAUuSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ1ygKXpEZZ4JLUKAtckhplgUtSo3oVeJJjktyQ5OEkDyX5saGCSZJW1uuamMAfAv9QVT+T5EjgqAEySZI6mLrAk7wYeA2wDaCqngOeGyaWJGmSPodQTgBGwJ8muSfJ1UletHSjJNuTLCRZGI1GPYaTJO2vT4FvAE4D/riqTgX+B7h86UZVtaOq5qtqfm5ursdwkqT99SnwvcDeqrpzfP8GFgtdkjQDUxd4VX0e+FySk8eLzgYeHCSVJGmivp9CuQS4bvwJlMeAN/ePJEnqoleBV9W9wPwwUSRJq+E3MSWpURa4JDXKApekRlngktQoC1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ1ygKXpEZZ4JLUKAtckhplgUtSoyxwSWqUBS5JjbLAJalRvQs8yRFJ7knyd0MEkiR1M8Qe+NuAhwZ4HEnSKvQq8CSbgHOBq4eJI0nqqu8e+B8AlwFfO9AGSbYnWUiyMBqNeg4nSfq6qQs8yRuBfVV110rbVdWOqpqvqvm5ublph5MkLdFnD/zVwHlJHgeuB85K8ueDpJIkTTR1gVfVFVW1qao2A+cDn6iqCwdLJklakZ8Dl6RGbRjiQarqVuDWIR5LktSNe+CS1CgLXJIaZYFLUqMscElqlAUuSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ1ygKXpEZZ4JLUqD5XpT8+yS1JHkzyQJK3DRlMkrSyPpdUex64tKruTnI0cFeSm6rqwYGySZJW0Oeq9E9X1d3j2/8FPAQcN1QwSdLKBjkGnmQzcCpw5zLrtidZSLIwGo2GGE6SxAAFnuTbgL8Cfq2q/nPp+qraUVXzVTU/NzfXdzhJ0livAk/yAhbL+7qqunGYSJKkLqZ+EzNJgA8AD1XV7w8X6eCz+fKPrMu4j1917rqMK6kNffbAXw38PHBWknvHP+cMlEuSNMHUe+BV9U9ABswiSVoFv4kpSY2ywCWpURa4JDXKApekRlngktQoC1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ1ygKXpEZZ4JLUKAtckhplgUtSoyxwSWrU1Bd0kKTWrNflEWFtLpFogR/EvBanpJX0vSr9liSfTvJoksuHCiVJmqzPVemPAN4HvB7YC3wyye6qenCocFofh9rLTOlQ1WcP/EeAR6vqsap6Drge2DpMLEnSJH2OgR8HfG6/+3uBH126UZLtwPbx3f9O8ukeY87CRuCZ9Q6xRg76ueVdvX79oJ9fT86vYXlXr/m9YrmFa/4mZlXtAHas9ThDSbJQVfPrnWMtHMpzA+fXOue3en0OoTwJHL/f/U3jZZKkGehT4J8ETkxyQpIjgfOB3cPEkiRNMvUhlKp6PsmvAh8DjgB2VtUDgyVbP80c7pnCoTw3cH6tc36rlKoa+jElSTPguVAkqVEWuCQ16rAs8EmnAEiyLckoyb3jn19aj5zT6nKKgyQ/m+TBJA8k+YtZZ+yjw/P3nv2eu0eSfHEdYk6tw/y+O8ktSe5Jcn+Sc9Yj5zQ6zO0VSW4ez+vWJJvWI+e0kuxMsi/JngOsT5I/Gs///iSn9Rqwqg6rHxbfcP1X4HuAI4H7gFOWbLMNeO96Z13D+Z0I3AMcO77/0vXOPeT8lmx/CYtvsK979gGfvx3AW8a3TwEeX+/cA87tL4GLxrfPAj643rlXOcfXAKcBew6w/hzgo0CA04E7+4x3OO6BH+qnAOgyv18G3ldV/wFQVftmnLGP1T5/FwC7ZpJsGF3mV8C3j2+/GHhqhvn66DK3U4BPjG/fssz6g1pV3QY8u8ImW4E/q0V3AMckefm04x2OBb7cKQCOW2a7nx6/xLkhyfHLrD9YdZnfScBJSf45yR1JtswsXX9dnz+SvAI4gf8vhBZ0md9vARcm2Qv8PYuvMlrQZW73AW8a3/4p4Ogk3zGDbLPS+f9vF4djgXfxt8DmqnoVcBNw7TrnGdoGFg+jnMniHuqfJDlmPQOtkfOBG6rqf9c7yMAuAK6pqk0sviT/YJJD5W/514EzktwDnMHit7sPtedvMIfKk74aE08BUFVfqKqvjO9eDfzQjLINocspDvYCu6vqq1X1WeARFgu9Bas5hcP5tHX4BLrN72LgQwBVdTvwQhZPBHWw6/K391RVvamqTgXeOV72xZklXHuDnoLkcCzwiacAWHJM6jzgoRnm66vLKQ4+zOLeN0k2snhI5bEZZuyj0ykcknwfcCxw+4zz9dVlfk8AZwMk+X4WC3w005TT6fK3t3G/VxNXADtnnHGt7QZ+YfxplNOBL1XV09M+2GF3SbU6wCkAkvw2sFBVu4G3JjkPeJ7FNyS2rVvgVeo4v48Bb0jyIIsvT99RVV9Yv9TddZwfLJbD9TV+678VHed3KYuHvd7O4hua21qYZ8e5nQn8TpICbgN+Zd0CTyHJLhbnsHH8HsWVwAsAqur9LL5ncQ7wKPBl4M29xmvgeZckLeNwPIQiSYcEC1ySGmWBS1KjLHBJapQFLkmNssAlqVEWuCQ16v8ACq9ULIbt+ssAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(a_proj_k.flatten(), density=True)\n",
    "plt.show()\n",
    "plt.hist(a_proj_c.flatten(), density=True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9599b8db-abfe-4859-9a90-80fa1edeaa60",
   "metadata": {},
   "source": [
    "### 2:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "937067f8-bbae-4f5b-8f3d-91be48700baa",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_ws = 500\n",
    "ws_nn = 4\n",
    "probs = [0, 0.01, 0.05, 0.1, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3eae029a-f641-4e64-9896-14760eada847",
   "metadata": {},
   "outputs": [],
   "source": [
    "# leave the main diagonal = 1 for now to make rewiring easier later\n",
    "ws_base = np.identity(n_ws)\n",
    "\n",
    "for i in range(ws_nn//2):\n",
    "    ws_base += np.roll(np.identity(n_ws), i + 1, axis=1)\n",
    "    ws_base += np.roll(np.identity(n_ws), -i - 1, axis=1)\n",
    "\n",
    "ws = np.tile(ws_base, (len(probs), 1, 1))\n",
    "\n",
    "ws_rewires = []\n",
    "for i, p in enumerate(probs):\n",
    "    rewired_edges = []\n",
    "    \n",
    "    for j in range(n_ws):\n",
    "        rewires = rng.binomial(1, p, size=ws_nn + 1)\n",
    "        possible_nodes = np.nonzero(ws[i, j, :] == 0)[0]\n",
    "        rewire_nodes = rng.choice(possible_nodes, ws_nn + 1, replace=False)\n",
    "        ws[i, j, j] = 0\n",
    "        \n",
    "        for k, (rewire, node) in enumerate(zip(rewires, rewire_nodes)):\n",
    "            neighbor = (j - k + ws_nn//2 + 1) % n_ws\n",
    "            if rewire and ws[i, j, neighbor] != 0:\n",
    "                rewired_edges.append((j, node))\n",
    "                ws[i, j, neighbor] = 0\n",
    "                ws[i, neighbor, j] = 0\n",
    "                ws[i, j, node] = 1\n",
    "                ws[i, node, j] = 1\n",
    "\n",
    "    ws_rewires.append(rewired_edges)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a21b646d-d083-49d1-884f-5d0d203f5488",
   "metadata": {},
   "outputs": [],
   "source": [
    "edge_bets = [nx.edge_betweenness_centrality(nx.from_numpy_array(net)) for net in ws]\n",
    "\n",
    "neighbor_avg_edge_bets = []\n",
    "shortcut_avg_edge_bets = []\n",
    "total_avg_edge_bets = []\n",
    "\n",
    "for edge_bet, rewires in zip(edge_bets, ws_rewires):\n",
    "    total = list(edge_bet.values())\n",
    "    shortcut = [bet for edge, bet in edge_bet.items() if edge in rewires or edge[::-1] in rewires]\n",
    "    neighbor = [bet for edge, bet in edge_bet.items() if edge not in rewires and edge[::-1] not in rewires]\n",
    "    total_avg_edge_bets.append((np.average(total), np.std(total)))\n",
    "\n",
    "    if len(shortcut):\n",
    "        shortcut_avg_edge_bets.append((np.average(shortcut), np.std(shortcut)))\n",
    "    else:\n",
    "        shortcut_avg_edge_bets.append((0, 0))\n",
    "\n",
    "    if len(neighbor):\n",
    "        neighbor_avg_edge_bets.append((np.average(neighbor), np.std(neighbor)))\n",
    "    else:\n",
    "        neighbor_avg_edge_bets.append((0, 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5731b266-9260-4550-b012-77e1f426254c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "edge betweenness\n",
      "probability  total\t      neighbors\t       shortcuts\n",
      "p = 0.00     0.0629+/-0.0619  0.0629+/-0.0619  0.0000+/-0.0000\n",
      "p = 0.01     0.0166+/-0.0235  0.0152+/-0.0204  0.1039+/-0.0310\n",
      "p = 0.05     0.0081+/-0.0098  0.0065+/-0.0071  0.0279+/-0.0141\n",
      "p = 0.10     0.0068+/-0.0057  0.0054+/-0.0043  0.0158+/-0.0057\n",
      "p = 1.00     0.0048+/-0.0014  0.0000+/-0.0000  0.0048+/-0.0014\n"
     ]
    }
   ],
   "source": [
    "print('edge betweenness')\n",
    "print('probability  total\\t      neighbors\\t       shortcuts')\n",
    "for i in range(len(probs)):\n",
    "    print(''.join([\n",
    "        f'p = {probs[i]:.2f}',\n",
    "        f'     {total_avg_edge_bets[i][0]:.4f}+/-{total_avg_edge_bets[i][1]:.4f}',\n",
    "        f'  {neighbor_avg_edge_bets[i][0]:.4f}+/-{neighbor_avg_edge_bets[i][1]:.4f}',\n",
    "        f'  {shortcut_avg_edge_bets[i][0]:.4f}+/-{shortcut_avg_edge_bets[i][1]:.4f}',\n",
    "    ]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "985ecf41-8d00-407c-8d70-b9c273431cfa",
   "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.10.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}