peeler_for_minimum_feedback.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Background\n",
    "## Using a \"peeling\" method to get layers in a directed, weighted graph. \n",
    "\n",
    "Here, I implement a model based on a probabilistic graph traversal model given in ([Schlegel et. al, eLife, 2021](https://doi.org/10.7554/eLife.66018)), which \"peels\" the graph to assign a numerical \"layer\" from a set of start nodes. \n",
    "\n",
    "The basic idea is that we have layers of nodes that we have visited at the same step, and the next layer is determined probabilistically based on the following arguments.\n",
    "\n",
    "Suppose we have a weighted, directed adjacency matrix entry $A_{ij}$, giving the weight of the directed edge that links node $i$ to node $j$. \n",
    "\n",
    "Given a node $i$ which we've visited, the probability of traversing to $j$ in the next layer $p_{ij}$, should be proportional to the importance of node $i$ relative to all of node $j$'s in-partners:\n",
    "$$\n",
    "\n",
    "p_{ij} \\propto \\frac{A_{ij}}{\\sum_k A_{kj}} = \\frac{A_{ij}}{s_{j}^{in}}\n",
    "$$\n",
    "where $s_{j}^{in} = \\sum_k A_{kj}$ is the in-strength. We can also define:  $s_{j}^{out} = \\sum_k A_{jk}$ as the out-strength. \n",
    "\n",
    "Additionally, we want to bias the probability to select for nodes that are more outward than inward. This is because a given node is more likely to be higher up in the topological sort if the outgoing edge strength is greater than the incoming edge strength. So, this follows:\n",
    "\n",
    "$$\n",
    "p_{ij} \\propto \\frac{s_{j}^{out}}{s_{j}^{out} + s_{j}^{in}}.\n",
    "$$\n",
    "\n",
    "Altogether, we implement the model with the probability of getting picked in the next layer:\n",
    "$$\n",
    "p_{ij} = \\frac{A_{ij}s_{j}^{out}}{s_{j}^{in}(s_{j}^{out} + s_{j}^{in})}\n",
    "$$\n",
    "\n",
    "\n",
    "## The data:\n",
    "\n",
    "I'm using the [FlyWire minimum feedback challenge data set](https://codex.flywire.ai/app/mfas_challenge) as the test case. A topological sorting of a connectome can be useful in studying information flow and distribution. The challenge was to sort the neurons in the brain such that the sum of the edge weights directed backwards is minimised. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "edge_list = pd.read_csv('./data/connectome_graph.csv')\n",
    "edge_list.rename(columns={'Source Node  ID':'pre', 'Target Node ID':'post', 'Edge Weight':'weight'}, inplace=True)\n",
    "edge_list.pre = edge_list.pre.astype(str)\n",
    "edge_list.post = edge_list.post.astype(str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_nodes = sorted(list(set(edge_list.pre.values) | set(edge_list.post.values)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Nodes that need sorting. Nodes that don't have any pre or post synapses do not need sorting.\n",
    "changeable_nodes = sorted(list(set(edge_list.pre.values) & set(edge_list.post.values))) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "node_name_dict = dict(zip(all_nodes, range(len(all_nodes))))\n",
    "\n",
    "edge_list.pre = edge_list.pre.apply(lambda x: node_name_dict[x])\n",
    "edge_list.post = edge_list.post.apply(lambda x: node_name_dict[x])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### benchmark:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "benchmark solution: 29023882\n"
     ]
    }
   ],
   "source": [
    "benchmark = pd.read_csv('./data/benchmark.csv')\n",
    "benchmark.rename(columns={'Node ID':'n_id', 'Order':'order'}, inplace=True)\n",
    "benchmark.n_id = benchmark.n_id.astype('str')\n",
    "benchmark.n_id = benchmark.n_id.apply(lambda x: node_name_dict[x])\n",
    "\n",
    "benchmark_dict = benchmark.set_index('n_id')['order'].to_dict()\n",
    "benchmark_id_order = benchmark.n_id.tolist()\n",
    "\n",
    "# how to calculate feedforwardness:\n",
    "edge_list['pre_order'] = edge_list['pre'].apply(lambda x: benchmark_dict[x])\n",
    "edge_list['post_order'] = edge_list['post'].apply(lambda x: benchmark_dict[x])\n",
    "\n",
    "edge_list['is_forward'] = edge_list['pre_order'] - edge_list['post_order'] < 0\n",
    "\n",
    "print('benchmark solution:', edge_list[edge_list['is_forward']].weight.sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the model:\n",
    "from src.graphpeeler import layer_realisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare probability of selection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "bodyId_posts = edge_list['post'].values\n",
    "bodyId_pres = edge_list['pre'].values\n",
    "out_only = list(set(bodyId_pres) - set(bodyId_posts)) # nodes that only have outedges.\n",
    "\n",
    "in_strengths = edge_list.groupby('post').weight.sum().to_dict()\n",
    "out_strengths = edge_list.groupby('pre').weight.sum().to_dict()\n",
    "\n",
    "weights = edge_list['weight'].values\n",
    "edge_df = edge_list.copy(True)\n",
    "\n",
    "# weight of edges as a proportion of their strength:\n",
    "edge_df['in_prop_weight'] = [weights[i]/in_strengths.get(bodyId_posts[i]) for i in range(len(edge_df))] \n",
    "edge_df['out_prop_weight'] = [weights[i]/out_strengths.get(bodyId_pres[i]) for i in range(len(edge_df))]\n",
    "\n",
    "# out-strength as a proportion of the total in and out strength:\n",
    "out_strength_prop = {i : out_strengths.get(i, 0)/(in_strengths.get(i, 0) + out_strengths.get(i, 0)) for i in list(node_name_dict.values())}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# combine to give: p_{ij}\n",
    "edge_df['in_div_out_prop_weight'] = [weights[i]*out_strength_prop[bodyId_posts[i]]/in_strengths.get(bodyId_posts[i]) for i in range(len(edge_df))]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# run the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 50 # number of realisations\n",
    "prepost = 'pre' # run the model forwards by taking the direction of traversal pre --> post. \n",
    "start_ids = out_only # start the 0th layer with nodes that don't have any preceding edges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_layers = []\n",
    "for i in range(N): \n",
    "    nth_r = layer_realisation(df=edge_df, start_ids=start_ids, prepost=prepost, seed=i, prop_weight='in_div_out_prop_weight', correction=1)\n",
    "    id_to_layer = {}\n",
    "    for key, val in nth_r.items():\n",
    "        for neuron in val:\n",
    "            id_to_layer[neuron] = key\n",
    "    d_layers.append(pd.Series(id_to_layer, name=f'r_{i}'))\n",
    "# 13 mins for N=50"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_layer_df = pd.concat(d_layers, axis=1)\n",
    "d_colname = f'm_ds_layer'\n",
    "d_layer_df[d_colname] = d_layer_df.mean(axis=1).sort_values()#.round(0).astype(int) # take the mean layer across all the runs\n",
    "d_layer_df.index.names=['bodyId']\n",
    "d_layer_df = d_layer_df[[d_colname]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEPCAYAAABV6CMBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAAsTAAALEwEAmpwYAAAynElEQVR4nO3deZxU9Znv8c/TGyiytA0CgjQ0Iio4UUFtlxHXTDQmGjUxyk00xoiZ3Iy5mZvEJHeM2SbJTGaSzMTR4BJNRtTEPY4ajYpLtFEaNSJIxJbGVpYWikUQ6O567h/nNNTWXXW6q6qX+r5fr35Vn985derXTVFP/7bnZ+6OiIhIrsr6ugIiIjKwKHCIiEgkChwiIhKJAoeIiESiwCEiIpFU9HUFCm306NE+efLkvq6GiMiA0tjY+J67j8l0btAHjsmTJ7N48eK+roaIyIBiZs1dnVNXlYiIRKLAISIikShwiIhIJAocIiISiQKHiIhEosAhIiKRKHAMAI3NMa59ciWNzbG+roqIyOBfxzHQNTbHmHtjA7va41RVlHHbZfXMqq3u62qJSAlTi6Ofa2jawK72OHGHtvY4DU0b+rpKIlLiFDj6ufq6Gqoqyig3qKwoo76upq+rJCIlTl1V/dys2mpuu6yehqYN1NfVqJtKRPqcWhwiIhKJWhz9XOrg+NVnzSC2fZdaHyLSZxQ4+rnEwfEdbXG+fe+rmKEZViLSZ9RV1c/V19VgCccOxB12aYaViPQRBY5+blZtNWOGD0krLzPTDCsR6RMKHAPAOYdPSDo2g++dPVPdVCLSJzTGMQBcdeYhAPxu8duUlxlHTKpm+rjhfVwrESlVanEMAI3NMbbsbGfLznZa39/Fo8vWceENDcpdJSJ9os9bHGa2CtgKdADt7j7bzPYF7gQmA6uAT7l7zMwM+AVwJrAduMTdl/RFvYulczrujrZ4Unln+hF1V4lIsfWXFsfJ7n64u88Oj68CHnf3acDj4THAGcC08Oty4Lqi17TIGpo2sDMlaIDSj4hI3+kvgSPV2cCt4fe3AucklP/GAw3AKDMb3wf1K5r6uhrKyyyp7EMTR3L7F7SGQ0T6Rn8IHA48amaNZnZ5WDbW3deE368FxobfTwDeTnhuS1iWxMwuN7PFZra4tbW1UPUumpMP3o8yAwOqyo2rPzZDQUNE+kyvxzjM7GCCLqTtwB3uvjniLU5w93fMbD/gMTN7PfGku7uZeZQbuvt8YD7A7NmzIz23P2lsjnHhDUG6kU5x4J4lLQAKHiLSJ3JucZjZ1Wa2Jhy47iw7DXgJ+CnwX8ASM4vU8e7u74SP64F7gaOBdZ1dUOHj+vDyd4ADEp4+MSwblO5Z0pIUNADaO5wFi1Yz90bNqhKRvhGlq+oM4HV335hQ9iOCrqbvEAxUTwGuzPWGZjbMzIZ3fg98GFgKPABcHF52MXB/+P0DwGctUA9sTujSGnS6aio52tRJRPpOlK6qyQQtAgDMbAIwC/h3d/9BWHYwwUD21TnecyxwbzDLlgpggbs/YmYvAr8zs88DzcCnwusfIpiKu5Kga+xzEeo/4Mzcf2SX5+JA9d5VxauMiEgoSuCoBhJbG8cT/PH7YEJZIzAv1xu6exPwoQzlG4BTM5Q78KVc7z/Qxbbv6vKcO1zzh9eYPm64xjpEpKiidFW1kjyD6WSgDViUUFYV8Z7SjWwtCnVXiUhfiNLieBn4uJnNBHYAFwDPuvsHCddMBgbtmEOxddfigCDZoRYBikixRWkd/AswEngFWBF+/2+dJ82snKD7anE+K1jK6utquv0H8gE70VhEBrKcA4e7PwOcBdxHMEh+vrs/nHDJcQRTY+9Nf7b0xKzaai4/sa7L8w7qqhKRoou0ANDdHwEe6eLcM8AR+aiUBBYsWs38Z5q6PF9Zrs2cRKT4+jw7rmTW2Bzj2/e+2uVaDgDr5pyISKFEmgFlZmVm9mUzazCzzWbWnnDuCDP7LzM7KP/VLD33LGnpNmgAtMddXVUiUnRRUo5UAY8BPwemEuyhkfhH71vApcDcPNavZOUy7m3ad1xE+kCUFsfXCNZufJdgxfeNiSfdfRPwNPB3+apcKetu1Xing/bbR4v/RKToogSOucCf3f177h4n8x/FbwGT8lKzEpdtDQdAVYXWWopI8UX55JkCNGS5ZiOwb5ZrJAfBBk7dX3NsXQ3XPrlSWXJFpKiizKraAYzKcs0kYFNPKyPJPH3H2CQ3PPsW7k5VRRm3XaYdAUWkOKK0OF4GPhwOkqcxs5EE4xsv5KFeJa+haQNZ4gYdcSfuylklIsUVJXDMJ9hE6TYzG5F4wsxGAbcQZNC9Pl+VK2X1dTXkOoRRXl6m2VUiUjQ5d1W5++1mdjpwCfBxIAZgZouBGcAQ4Fp3f6gA9SxJZWVlEN/T7jCgotxo60iel3D+rInqphKRook0LcfdLyVYq7EMGEPwWXYkwcZKn3f3L+e9hiWqoWlD2raxDmlBwwzOO3JiEWsmIqUucsoRd78FuMXM9iLomtrs7tvyXbFSV19XQ5lBPMtKwLM/tL9aGyJSVD1eCODuH7j7uwoahTGrtpqpo4dlve6dTR9kvUZEJJ+0gqyfWrBoNW+0JsfkcSOGpF23eFVM6zhEpKi67Koys67zeXfP3X1qD58roZufTf/1n3P4BOY/05TUfdW5J4e6q0SkWLob4ygjPa1IFTA+/L4DeA8YDZSHZWuA7LkyJDtLTpo+engVV515CJNqhvGte19NOvfGuq3FrJmIlLguu6rcfbK7T+n8Aj5EsMNfA0Gyw6HuPh4YCpwCLAJagL8pfLUHv9MO3i/p+PwjgplTFx0zKa3L6uW3NxWrWiIikcY4fkiQcuQkd3/K3TsA3L3D3RcSBJN9w+ukFxqbY9zy/Krdxwbc8vyq3WMZqYv9Dj9gVPEqJyIlL0rg+ARwv7tn7Ipy9x3A/cC5+ahYKUtdw+EkpxXZvqsj6fq9h2gjRxEpniiBowaozHJNZXid9EL13lVp6zc604o0Nsf40/J1Sefe27qziLUTkVIXJXC8CZwfJjNMY2bVwPlAT2djSSjTXhxzDhrDrNpq7lnSkhZU/rR8nabkikjRRAkc1wP7Ay+Y2WfNbLKZ7RU+XkwwOD4OuLYQFS0l1XunJyDunGOVaSF53FF2XBEpmihJDn9pZtOALwO/znCJAf/p7v+Vr8qVqtfe3ZxW1hkwzjtyIgsWrU47nynYiIgUQtQkh1cCxwM3Ay8RdEu9BNwEnBCej8zMys3sJTN7MDyeYmaLzGylmd3ZuQeImQ0Jj1eG5yf35PX6u9YMYxb7DQ+m4M6qrebDh45NO5/LVrMiIvkQOeWIuz/v7l9w99nuPi18vNzdn+tFPa4Elicc/wT4mbsfSJC+/fNh+eeBWFj+s/C6QWf08OR1GmUG5yZkwJ03ZyrlyesD1eIQkaLp81xVZjYR+ChwY3hsBAsK7wovuRU4J/z+7PCY8Pyp4fWDynlHTqQqjAzlBj8457CklCKzaqv59NGTdh+XmVocIlI8kQOHmdWb2Y1m1mhmb5rZEjO7wcyO62Edfg58HXbvlFoDbHL39vC4BZgQfj8BeBsgPL+ZDNN/zexyM1tsZotbW1t7WK2+M6u2mms+PpO/nTaa759zGBcdMyntmnOPnMjQyjLKDaoqtAOgiBRPpJVjZvYD4JvsmeTT6XDgUjP7ibt/K8L9zgLWu3ujmZ0UpS7dcff5BFvdMnv27Cw7WvQ/jc0xrvnDa7S1x1n01kamjxue1OJobI5x95IWTpw2hjHDh3DukdoBUESKJ+cWh5l9EvgWsBq4DKgD9gofLwvLv2Fmn4rw+scDHzezVcAdBF1UvwBGmVlnUJtIkCOL8PGAsD4VwEhg0M1DvWdJC7va4ziwqz3Or556c/e5xuYYF85/ngWLVvPosnX8vrGl7yoqIiUpSlfVl4F1wFHufrO7r3L3neHjzcBRQCvwpVxv6O7fdPeJ7j4Z+DTwhLvPBZ4kWEwIcDFBKhOAB8JjwvNPuPuAa1Fkk/oDPf76+t0L/BqaNiRtH5uYikREpBiiBI4PAXe5+3uZToblvyfotuqtbwBfNbOVBGMYN4XlNwE1YflXgavy8Fr9zsz9kxfnx+O+OzjU19VQmTClqrzcNL4hIkUVJXBUANuzXLOdHuxjDuDuC939rPD7Jnc/2t0PdPdPuvvOsHxHeHxgeH5QpjdJnSFVXrYnOHQOnJeXBcHDHVas1X4cIlI8UXNVnWVmGZ8Tlp8ZXie9UL131e51GuUG3zt7ZtLgd2z7LuJhwqqOuPNP9y9VrioRKZoogWMBcAhwf5h6ZDczm0qwruLQ8DrpocbmGNc8sJTOYYx4hmvq62ooS5jX1hF37lmiQXIRKY4ogePfgacJFustN7PVYdqPZmAFwSK9P4fXSQ81NG1gV8LgtztpLYpZtdVp02/XK7W6iBRJzoEj3MDpdODbwFsE02SPIpge+1ZYfmpXGz1JburratLSiSQOjnelZWO24ScRkfyImuSwzd1/5O7TgBEEQWNEmLPqR+7eVpBalpBZtdV8/5zDkv5hKjOsDN/ZntyJ9fq6rRrnEJGi6PGeo+7+PvB+Husioc4UI3e+uJqxI4Yyb87UtK6pC46axCstr+4pCPfk0ApyESm0Pk9yKOkam2Ncff+rvNKymSdeX5fxmunjhjMuIYvukErlqxKR4ogUOMxsjpk9aGbrzazNzDoyfLVnv5N05/qn3qSzJ6o9HhwnamyO8anrn2NtwoD4R2aMU2tDRIoi564qM/socB9QTpCXagWgIFEAb723Lel4/ZYdSccNTRvoSMlLorQjIlIsUcY4rgHagI+6+6OFqY40Nsd4673koaMLjkpOq56pS2pIRXlB6yUi0ilKV9VM4E4FjcK6Z0kLHQkTpg4ZNzxtP45M28e2bPpAs6pEpCiiBI73gY2FqogEUjPj/nX9+xkDwrw5U3fnqwJwz77WQ0QkH6IEjseBYwtVEQmcd+TEpHQiXQWEWbXVfP/smVSUGWXaBVBEiihK4PgGMNXM/t9g3Oe7P0n87VaUdx0Qpo8bzskH70fd6GGcOG1MkWonIqUuyuD4d4DXgO8SbBP7MrApw3Xu7p/vfdVKU+oYx0kHjck4zbZzJ8DOvFYrW7ex8K+t3P6Fek3LFZGCihI4Lkn4fnL4lYkDChw91JpjssLUZIiwZzdABQ4RKaQogWNKwWohu41OWA2e6bhTfV0NZSSnXTfLPFVXRCSfcg4c7t5cyIpIIHXb2NTjTrNqq5k+bjjLE3b/mz52uFobIlJwylXVz8S276JzbNxI30Y2UVVFWbfHIiKFoE+afqZ676rdazk8PO7KsSndUqnHIiKFoMDRzyxcsb7b40Rbd7Z3eywiUggKHP1MU+v73R4nSt0u9q/rtnZxpYhI/ihw9DN1Y/bp9jhR6irMxuaY8lWJSMEpcPQz8+ZMpTLcdLyy3Jg3Z2qX16bmtYq70quLSOF1OR3XzDYCP3b3fwmPrwYWuvvTxapcKZpVW80dlx9LQ9MG6utqup1euzllxlWZ1nGISBF01+IYBQxNOL4GOKmAdZHQirVbaWjawIq13Y9Z7GyPJx3XjR6mdRwiUnDdLQBcB0wsVkUksGDRar5176sAPPPGewBp+3F0uuCoSbzS8uru45kTMi8WFBHJp+4CRwPwGTPrANaEZSflkBjX3f37+ahcKXp46Zq0464Cx0XHTOKFtzZw38vvAnDfy+9y9JSaLq8XEcmH7gLH14CDgHkJZSeRvbvKgZwCh5kNBZ4GhoR1ucvdv2NmU4A7gBqgEfiMu+8ysyHAb4BZwAbgAndflctrDRQ1w5IX/M0YP6Lb61P3J7/zxdUKHCJSUF0GDndfaWaHESQ3nAAsBG4Bbs3j6+8ETnH3982sEnjWzB4Gvgr8zN3vMLPrCbLtXhc+xtz9QDP7NPAT4II81qdPNTbHePAvyS2O4XtVdvucsSOGAptTjkVECqfbJIfuHgfeBN4Mu6hWuftT+Xpxd3eCLWkBKsMvB04BLgrLbyUYmL8OODv8HuAu4JdmZuF9BryGpg10xPf8KBVllnWW1Lw5U3lyxXraOjzr9F0RkXyIkh23IGs+zKycoDvqQOBagkC1yd0782e0ELR4CB/fDuvTbmabCbqz3ku55+XA5QCTJg2cbpv6uhrKymx38LjshClZZ0lFmb4rIpIPUfbj2M3MJgJHEEzZ3QwscfeWntzL3TuAw81sFHAvcHBP7pNyz/nAfIDZs2cPmNbIY6+tTWpx3PDsW5w+Y1xOwUMBQ0SKJVIrwsxqzewRoBm4j2DM416g2cweMbPJPa2Iu28CngSOBUaZWWdQmwi8E37/DnBAWJcKYCTBIPmg8Mhra5OOO+Ke00rwxuYY1z65UulGRKQocg4cZjYOeBb4MEHg+C3wL+HjqrD82fC6XO85JmxpYGZ7AacDywkCyPnhZRcD94ffPxAeE55/YrCMbwB8ZEbyr648hzGOxuYYc29s4N8eXcHcGxsUPESk4KK0OP6JYIzhG8A0d7/E3b/p7pcQTNv9OrA/8P8i3HM88KSZ/QV4EXjM3R8MX+OrZraSYAzjpvD6m4CasPyrwFURXqvfO33GOMrDfxEz+P7ZM7N2QTU0bWBnW5y4w862OHcv6VGPoYhIzqKMcXwUeNTd/zX1RDhO8VMzOw04C/jfudzQ3f9CMFaSWt4EHJ2hfAfwyQh1HlAamjbQ2X4qo/vd/zqlbvz0+8Vvc96REzXmISIFE6XFMY5g9lN3GsPrpAfq62qoKC/DgPLyspwSFqYGl/aO3MZFRER6Kkrg2AzUZrlmEomr0SSyuDsePuYidWvZsjJlyBWRwooSOJ4Fzjez4zKdNLNjCLqRns1HxUrRPUtaaO8IAkZ7h3NPDuMVr72bGqez5hITEemVKIHjh+HjU2b2WzO71MzOMLPPmdmtwDPh+X/ObxVLR2obI1ubo7E5xu8Xv538HFdXlYgUVpSV40vM7HyCFCBz2ZMSBII/czcCl7p7tnEQ6cLM/Ud2e5yqoWkD7fHk8FJVkdvYiIhIT0VaOe7uD5rZJIKcUUcSLMDbDLwE3Ofu27p7vnRv4Yr1acfdZbqtr6uhoszYFXZvlZcZV581QzOqRKSgIqccCYPDgvBL8qgpJUV66nGqWbXVfHL2Ady2aDUA8bizNG3MQ0QkvwqSuFB6pi1lK9jKsuwD3TMSurM613Fo9biIFJICRz/R2Bxj9cbtSWVt8exTclNnVWkdh4gUmgJHP9HQtCFtFlXd6GFZn/fGuq1Jx2ZaxyEihaXA0U+kLeQzctqUaWdK91bd6GEaHBeRglLg6Cdi23clLd379NGTcgoAFxyVPOvqtEPG5rlmIiLJFDj6icRkhZB9DUdXbvrzWxocF5GCirIfxxNm9v1CVqaUpQ5yp6cSyezhpWuSjts0OC4iBRalxVEPlBeqIqUudWB8/dadOT3vjJnjk45z2fxJRKQ3ogSONwi3bZX8W/ZOcgvjT8vX5dTldNExk7jixLqEkkGzIaKI9FNRAseNwEfDlCOSZ8vWbEk6difnLqfEFeYdcfjVU2/mtW4iIomipBz5A8Ge4H82s58QbPW6lgx/4rr76vxUr3TsO6yKtVuSu6dy7XJat2VHt8ciIvkUJXA0EQQJA37RzXUe8b4CnHLIWBYs2hNvTz90bM7rMY6tq+GVls1JxyIihRLlA/43qAO9YFKn3548fb+cnzt8r8puj0VE8inKfhyXFLAeJa9zAWBnky51L/HupK46Tz0WEcknLQDsJxIXADrRPvwTV51HDToiIlH1aCzCzA4GDgH2cfff5rdKpSl1wV/UFkdPg46ISFSRWhxmdriZLQZeA+4Cbkk4N8fMtpvZx/JbxcGvsTnGHS8kT0Tb+kFbzs/v6apzEZGeiJJy5CBgITCdYFbVwymXPE2w7/j5+apcqWho2kBHyrSD11LWdXSnp6vORUR6IkqL4ztAFXCMu3+VYB3Hbu7uwPPAUfmrXmmor6uhPGWzv9RUIt1JnZH1xOu5rToXEemJKIHjVOAed1/WzTVvA/v3rkqlZ8XarbtbDQZccWIdFx2T+wL91H3GO+Jw95KW/FVQRCRBlMBRDWT7NDKCVonkqLE5xj/d9yqdu8Q6sHVne6R7ZNqZPPtu5SIiPRMlcKwDDsxyzQyCVkdOzOwAM3vSzJaZ2WtmdmVYvq+ZPWZmb4SP1WG5mdl/mNlKM/uLmR0Zof79UqbxjairLGekdFWVlxnnHjmxdxUTEelClMDxBPAxM5ue6aSZHUXQnfXHCPdsB/7R3Q8lSNv+JTM7FLgKeNzdpwGPh8cAZwDTwq/LgesivFa/VF9Xk/aPMGJItFnSqbsHHjlplLaPFZGCiRI4fkTwQf+0mX2RcCzDzGaEx38AtgI/zfWG7r7G3ZeE328FlgMTgLOBW8PLbgXOCb8/G/iNBxqAUWaW+yhyP5Xawng+4kZM9XU1lCf8S764KpaU90pEJJ9yDhzuvgI4j2AM45fAZQRd6X8Brg3Lz+1pZlwzmwwcASwCxrp759Z2a4HOjbQnkNwV1hKWpd7rcjNbbGaLW1tbe1Kdomlo2pAWOIZURFvQP6u2msmj90kqu/NFBQ4RKYxIfSLu/oiZTQEuJuhaqgE2Aw3Ar919Y08qYWb7AHcDX3H3LWZ7Ol7c3c0sUre/u88H5gPMnj27XydmzLTKe+uOaIPjAG3t8aTjTdtzX0AoIhJF5JQj7r6JYAFgd6nVc2ZmlQRB4zZ3vycsXmdm4919TdgVtT4sf4fkXQgnhmUDVqbUIs0bt0e+z872jm6PRUTypU+THFrQtLgJWO7u/55w6gGCVg3h4/0J5Z8NZ1fVA5sTurQGpEybNZX34F/lnMMndHssIpIvkT+izGyumT1uZhvNrD18fNzM5vbg9Y8HPgOcYmYvh19nAj8GTjezN4DTwmOAhwg2lFoJ3AD8fQ9es98pS1l0Mffo2l7fc1LNsF7fQ0Qkk5y7qsIupbuAswgGxTuAVmA0cDJwkpl9Cjjf3XPqYHf3Z+l6rdqpGa534Eu51nkgSN1X/JBxw7nqzEMi3aOxOcavnm5KKrvzxdWRVp+LiOQqSovjm8DHCGY9nQwMdffxwFDgFOAFgqDyjXxXcjCrr6shYS4Ab763LXKeqUwzs8aOGNr7yomIZBAlcHyWoIvoJHd/yt07ANy9w90XAicRdCNdkuc6Dmor1m6lI2FCVFt7PK0Vkk19XQ1VCVkSywzmzZmaryqKiCSJEjgmAve7e8Ydhtx9J8EgtkZlI3h4afLYfk82YppVW801H5+5O8OuEQQkEZFCiBI43gUqs1xTGV4nOZoxfkRaWU+2fo1t37U7UWKHw9X3L1VqdREpiCiBYwFwvpmlf9IBZjaKYBOn2/JQr5KxZHXyh3uZZZ6im02QdmRPd1XcPXKXl4hILqIEju8Bi4EXzOwiM5toZpXh41yC1eMvAN8vREUHo8bmGC+uSg4c40ft1aMEhbNqq7nshCmUWZjbvqKsRwFIRCSbLqfjmlmczBm+DfhtF+XTgA+6u6/scc+SlrRfcKauq1w0Nse4+c9vEffgH+KSYycrQ66IFER3H/BPE31rCMlRY3OM2zNksD1p+n49ut/dS1rYFW7s4cD1TzcxqWaY1nKISN51GTjc/aQi1qPkXP/Um8QzlPdkYBwyr6L89r2vMn3ccLU8RCSv+jRXVSl7671taWXlZdbjcYlMO/456SvTRUR6S4Gjj9SNTs8lddj+I3rcOphVW520CLCTBshFJN8iD2Kb2ceAwwkWBGZa1+Hu/vle1mvQmzdnKo8uW5dUtqoH6dQTTdp3b1a27mnJTBg1VN1UIpJ3UZIc1hJsDzuDrhMTQtBDosCRxazaaqaNGcYbCR/0+1SV9+qepx0ylpWte5Idfuxv9u/V/UREMonSVfUfwEzg1wR5qaYBUzJ81eW3ioNTY3OMrTuTd/pr2bSjV3uFD98ruQG4ZWf0nQRFRLKJ0lV1CvBHd7+sUJUpFY3NMebe2MCOtvR5Vb1Jh771g+Rs9ne8sJrzjpyo7ioRyasoLY424NVCVaSUNDRtYGeGoAEwpKLn8xVeW7Ml6TjuwSJDEZF8ivIp9WeCrirppfq6mrRd/zodOHZ4j+97xszxaWVawSki+RYlcFwNnGhmny5UZUrFrNpqTj1kbFp5eRmcl2E9Rq4uOmYSV5y4Z4ipopf3ExHJJOcxDnd/ycxOBf7HzOYBS4DNmS91JTrMItM6jguOmtTr8YjTZ4zj5udW0dYep6xMy3REJP+iTMcdCfwzsC8wJ/zKxFGG3G41Nse48dm3ksoqyiwvrYOGpg20tcdxoD3cTVCD4yKST1FmVf2MYK/xPxFkx30X0HzPHmho2kBHPHn04Xtnz8zLB3z13lW7xzXiRN9NUEQkmyiB4yzgOXf/cKEqUyo6B8c7EmLH6g3puat64rfPr0o6vu+lFmXIFZG8itIJvhfwXKEqUkpm1VYzMmWx3n0vv9Pr+zY2x1iestf46l6mMRERSRUlcLyEVoXnRWNzLGmbV4B9h/W+SylTJlwlORSRfIsSOL4PfMzMTihUZUpBY3OMC+Y/T+v7yftuHD6p9+Mb9XU1aUnEVqS0QEREeivKGMd44EHgCTNbADSSeTou7v6bPNRtUPrVU2/S3pG+LC8fM6pm1VZTPaySjdv2pB5ZvnYrjc0xzawSkbyJEjhuIZhqa8Bnw6/UT0ALyxQ4urBuy46M5SvWbs3Lh/uJ08Zw38vvJpX95OHl/O6K43p9bxERiBY4PpfvFzezmwlma61395lh2b7AncBkYBXwKXePmZkBvwDOBLYDl7j7knzXqdAuOGoSr7Skp/zqTXLDRNMypCx5sTmmVoeI5E2UleO3FuD1bwF+SXIL5SrgcXf/sZldFR5/AziDIJX7NOAY4LrwcUDpatrtfiOG5uX+GddtOFoIKCJ506c5Kdz9aWBjSvHZQGeQuhU4J6H8Nx5oAEaZWXpWv37ukdfWZiy/Ys7UvNz/tXfTh53MNLtKRPKnPyYzGuvua8Lv1wKd2QAnAG8nXNcSlqUxs8vNbLGZLW5tbS1cTXvg8ANGpZVNqN4rb62BTNlw60YPU2tDRPImSq6qpuxXAUGSw7z8+ezubmaRM4O7+3xgPsDs2bP7VWbxTGMQI4ZE3vq9S+cdOZE7Fq0mcbePmRNG5u3+IiJRWhxlBLOmUr+qCQayJwNVEe+ZybrOLqjwcX1Y/g5wQMJ1E8OyASXTGERVLzZvSjWrtprTDk1O2f7AK+/S2BzL22uISGnL+RPL3Se7+5QMX9XAQcAjwJvAIb2s0wPAxeH3FwP3J5R/1gL1wOaELq0BY+GK9Wllx+Z5/GHM8CFJx9oJUETyKS9/6rr7SuBcgjGH7+T6PDO7HXgemG5mLWb2eeDHwOlm9gZwWngM8BDQBKwEbgD+Ph91L7YX3kpPCzI8JW9Vb83YP71r6q/rtIJcRPIjb53r7r7DzB4DLgS+meNzLuzi1KkZrnXgSz2vYd9bsGg1mz5IzkRv5H/G030vpbcuXlyltRwikh/5nlXVDozL8z0HjYeXpvesHTW5Ou8f5ivXv5+x/G51V4lIHuQtcJjZaOATJE+ZlQRrN32QVvaNM3o7JJRu3MjMiwlTEyCKiPRElOm4V3dzjwMIFuiNJMduqlKzYNFq3mjNz2ZN2Uys3ptla9LHNDKNfYiIRBVljOOaLOe3AD9w93/peXUGrztfXJ2xvBCpQEanzKrqlGlVuYhIVFECx8ldlMeBGPC6u2sP8i7sao9nLC/EnuCZFgECrN+6M++vJSKlJ0qSw6cKWZHB7t3N6eMbALHtuzKW98as2mp+/8Xj+NJ/N7I2IVg88fp6zawSkV7rj7mqBp3G5hjv78jcGCtU8sFZtdVMG5ec3qQj7hm3lxURiaLbFoeZ9SiwuHvmfpkS1NgcY+6NDWTY9I8Jo4YW9K//mgz7mCtLroj0VrbA0NaDr/z3vQxgDU0b2NGWOY5+6eRpBX3tpe9uSSt7rIu07iIiuco2xvE2mTN1Z7IPoD9nU3Q3+D19XHqm3HxpbI7RlGEh4K+ebuL0GeM0ziEiPdZti6ObxIa7vwgSHP5nwr1WFbjOA0Zjc4xrHlia8ZxBQccbGpo2pM2qguCvAK0gF5He6NXguJl9ElgO/CvBZ+HX6X123EGjoWkDbZkGN4DKirKCjjfU19V0uVL8PU3LFZFe6FHgMLPjzOw54A6CVeP/AUx195+6u8Y4QvV1NVSWZ/74vv0L9QXtLppVW83+1XtlPPf46+u0P4eI9FikwGFmU83sLuAZoB64GzjU3f+Pu+uTKMWs2moOHT8irbzMKMoYw8cOy7wle0ccfvXUmwV/fREZnHJaAGhm+xLsszGPYJe/54F/dPeGAtZtwPvxQ8t5uSU9zccJB44uyut3t89HU2vmDLoiItlkW8dRBXwFuAoYRbDD31XufnfBazbANTbHuP7pzNu0X3naQUWpQ3czuirLtfZTRHom26fHCuBHQAdBADlYQSM33c1cmntjQ1HGGLpLZ7J87VYWLMqceFFEpDvZAkdt+GjA/wXeMrPVWb6aC1vlgWFlN1u1trXHi5L6o76uhopu/oV/8D/LNEguIpHl0l9hwL7AxBy/DihITQeQxuYYf8kwttGp0FNxO82qrebOeccxfGjmHsntuzq48IbitH5EZPDItgCwrCdfxap8f9TYHOOCXz3Hji7SqB89uZrbLivsVNxEs2qrOfXg/bo8X6zWj4gMHiX9IV8IDU0b6CJmMG74EH53xXFFT/exYVvXYx0ObP2gjWufXKmWh4jkJMpGTpJFY3OMB195t8vzpx46toi12eOMmeN55o33ujx//dNNlBlUVZQVtTUkIgOTWhxdWLBoNZ+5aVHOM48am2N88rrnWL6260Hxc4+cmK/qRXLRMZM4enL3wSDu6rYSkdyoxZHBgkWr+da9rwLs/kv9omMmdfucr975csakgp3OOXz/Pv1L/htnHMJ51z3X7TVxgm4rEZHuqMWRwZ0vJrcybn4280K+xuYY37r3Vep/+CeaN27v9p7DhvRtjJ5VW82Hs3SVuQfdVpf/ZrHGO0SkS2pxZLDfiKHAnum0K1u3pe3V3dgc48L5z7Ori+y3qXLd1KSQ5s2ZymPL1mWty6PL1vHosnUcuN8+XHr8FC46ZhKNzTEamjZQX1ejMRCREmfu/eEjrXBmz57tixcvjvScxK6qTrX77s2njjqA6r2rWPruZh546R3e39WR8z3v/mLxZ1Nlkulny6ay3Hanh68qN26//Nh+8bOISOGYWaO7z850Ti2ODDKl6mjeuJ1//eOKHt3vihPr+s0H7UXHTGL1hm1d5tHKJHFPkV0dznnXPceIoRUcfsAoNn/QxtgRQ5k3Zyor1m7lzhdXs6s9TlVFGRccNYnp44Zz95IWVq7bysZtu6gbsw/z5kwF2N2CSfy+v/yeRKRrAy5wmNlHgF8A5cCN7v7jfL9GvlZ1Txw1lL8/eVrWgfViu+rMQ1i7ZQf3vdz11OFstuxo5+ndU3w38+iydWnXvNKS3rJZ2bot6VojuRuvc/+SMjNGDK1gyuhhjNq7ijHDh3DukRP57fOrWPjXVv5mwkiOqatJCjad3WlbP2jj+aYNuwNapmCUqestW3dc1O66BYtW8/DSNZwxc3y/ew+I9MaA6qoys3Lgr8DpQAvwInChuy/r6jk96aoCmPbth7rcvS8X//yJw/rlh0Xih99vn1/Vq+AhpcOA0ftUsXHbLhL/W5QbjBk+hA/aOti2s4OqcmNnR5yOePI1e1WWs7MjTtydEUMrMIzNH7ThHszm27uyjBOmjeGk6fux9N3NvNQcoyW2naGV5YzeZwhVFWVMGT2Mpe9sBjMuPX4K08cNj9RqbWyOJbV+K8vLaOuIs++wqt1/nMzYfySx7bvS/qC4e0kLRjClfsXarVn/IGhsjnH9U2+yfssOLjhqUo8/CxL/v+byM3b13J605AdTV9XRwEp3bwIwszuAs4EuA0dPVZWX0daR+xhGov4cNObe2MCu9jgVZQZmaX/xi2TiQOv76V24HQ5rt+zZirg9nv5u6nCSxgNj29vTrtneFt89KSPR1p0du1/3lYT8b9+691Uqy42OuO9+L7d3xLtcxNrtZJbWbUmHiYthgaTn3fHi6t1Bsaup+p1phzozSHS2vKN+JmT6/9rdz9jVcwuxsHegTcedALydcNwSliUxs8vNbLGZLW5tbe3RC40flXnb1e6MGFrB3V88rl8GDQj+WtnVHg8W+3U4be1xnIH3JhCB4D2c+F7ubhFrQ9OGnHsQEu+T+ryOlMVaDy9dk/G1UtMOZboum0z/X3NdqJv03AIs7B2UnxnuPt/dZ7v77DFjxvToHpcePyXyc64645B+PbhbX1dDVUUZ5RaMJVSG31dVlvHPnziMQ8YN7+sqiuSsstzS3stdZZ6ur6vZPX6WTVnCfVKfl7r/2Rkz07dnzrSdQabrsunq/2su2bWTnluAbNwDbYzjWOAad/+78PibAO7+o66e09MxDggGN3/4P8vYFjazDRhSUcbYEUOpLDfqxuxD3ehhvLZmy4AZAM3WZ5o6oPvjh5bzu8Vvs6s9zrZdHerWKlEa49AYR9K5ARY4KggGx08F3iEYHL/I3V/r6jm9CRwiIqVq0AyOu3u7mf1v4I8E03Fv7i5oiIhI/g2owAHg7g8BD/V1PUREStWgHBwXEZHCUeAQEZFIFDhERCQSBQ4REYlkQE3H7QkzawWaIzxlNND1Bt2lRb+LZPp97KHfxR6D9XdR6+4ZV1AP+sARlZkt7mrucqnR7yKZfh976HexRyn+LtRVJSIikShwiIhIJAoc6eb3dQX6Ef0ukun3sYd+F3uU3O9CYxwiIhKJWhwiIhKJAoeIiESiwBEys4+Y2QozW2lmV/V1fYrNzA4wsyfNbJmZvWZmV4bl+5rZY2b2RvjYf3eqyjMzKzezl8zswfB4ipktCt8jd5pZVV/XsVjMbJSZ3WVmr5vZcjM7tlTfG2b2f8L/I0vN7HYzG1pq7w0FDoIPCOBa4AzgUOBCMzu0b2tVdO3AP7r7oUA98KXwd3AV8Li7TwMeD49LxZXA8oTjnwA/c/cDgRjw+T6pVd/4BfCIux8MfIjg91Jy7w0zmwD8AzDb3WcSbO/waUrsvaHAETgaWOnuTe6+C7gDOLuP61RU7r7G3ZeE328l+GCYQPB7uDW87FbgnD6pYJGZ2UTgo8CN4bEBpwB3hZeU0u9iJHAicBOAu+9y902U6HuDYDuKvcKN5fYG1lBi7w0FjsAE4O2E45awrCSZ2WTgCGARMNbd14Sn1gJj+6peRfZz4OsEO5sC1ACb3L09PC6l98gUoBX4ddh1d6OZDaME3xvu/g7wU2A1QcDYDDRSYu8NBQ5JYmb7AHcDX3H3LYnnPJi7Pejnb5vZWcB6d2/s67r0ExXAkcB17n4EsI2UbqkSem9UE7S0pgD7A8OAj/RppfqAAkfgHeCAhOOJYVlJMbNKgqBxm7vfExavM7Px4fnxwPq+ql8RHQ983MxWEXRbnkLQxz8q7J6A0nqPtAAt7r4oPL6LIJCU4nvjNOAtd2919zbgHoL3S0m9NxQ4Ai8C08KZEVUEg10P9HGdiirsw78JWO7u/55w6gHg4vD7i4H7i123YnP3b7r7RHefTPBeeMLd5wJPAueHl5XE7wLA3dcCb5vZ9LDoVGAZJfjeIOiiqjezvcP/M52/i5J6b2jleMjMziTo1y4Hbnb3H/ZtjYrLzE4AngFeZU+//rcIxjl+B0wiSE//KXff2CeV7ANmdhLwf939LDOrI2iB7Au8BPwvd9/Zh9UrGjM7nGCiQBXQBHyO4A/PkntvmNl3gQsIZiK+BFxGMKZRMu8NBQ4REYlEXVUiIhKJAoeIiESiwCEiIpEocIiISCQKHCIiEokCh4iIRKLAISXBzNzMNPecYB2CmXV0pkE3s/3D38/X+7puMjAocIiUnlOBl9w9lnAM8EQf1UcGGAUOkRISZrU9muQgcRqwCVjSF3WSgUeBQySFmZ1jZv9tZn81s23hV6OZ/YOZlaVce3vYzTOni3udF57/ZUr5vmb2o3A3vQ/MbLOZPW5mH85wj0vCe1wS7lS5MLw+p643M6s1swPN7EDgPKASeCOh7FTgZaAuLBvUKcGl95RyREpC54esu1sO175OkK+rkSDL6UiCDLkHAf/t7p9JuPZE4ClgQZgIMfVejwKnA4e5+9KwrBZYCEwmyA/WSJCe+yxgHDDP3W9IuMclwK+B/yFI4f0wQWK9Wnf/dA4/zyqgNtt1CZ5y95MiXC8lRoFDSkLEwDHV3d9MKSsj+PD+LFCfkGIcM1sKHAhMcPcNCeV1wErgeXc/PqF8IcGOehe5+x0J5aMIAsp0YLK7rwvLLwlf24Ez3f2RCD86ZnYGQWCCIJFnDPhueHwWQTbXLwLvhWWt7v5UlNeQ0qKuKpEUqUEjLIsT7MkB8Hcpp68DhgCXpJR/ATDgV50FZvYhYA5wd2LQCF9jE/AdYChBl1Kq+6MGjfC+D7v7XcCfgPHAH9z9rrCsCljn7td3liloSDYV2S8RKS1mVgN8DTgTqGPPX+udUscAfgP8GLgc+LfwHpUEgSRGkHq807Hh40gzuybDy48JHw/JcO6FnH6Ars0h+GPxyZSyp3t5XykxChwiCcLuohcJtgZ9gSAobCTYe2EUcCVB62I3d99qZv8NXGFmJ7v7k8DHCcYrfu7uOxIurwkfTw+/urJPhrK1Pfh5rkk4PKnztc3seGBvgu1PaxKuW+juC6O+jpQWBQ6RZJcRBI3vuvs1iSfM7FiCwJHJdcAVwDyCv+jnheXzU67bHD5e6e7/EbFuPRmQ/E6Gsq+lHJ8SfnVa2IPXkRKiMQ6RZAeGj3dnOJdxyi2Au/8F+DPwCTM7hmBtxNPuvjzl0obw8W97W9FcuLuFEwJGAR0EAbGz7E5gbedx+HVNMeolA5sCh0iyVeHjSYmFZnYE8M0sz72OYLD5boJB8etTL3D3xQRTcM81s0sz3cTMDjOz/SLVOrsTCbZFXphQNodgKrFIJOqqkpJiZrd0c/rvCcY0vgb83MxOBt4AphFMW72HYK/prvwe+BnB4Pl74fWZXESwcvsmM/sHgn3dNwETgb8BZhIMoq/P5WfK0cnATsIWj5kdTDAGszCPryElQoFDSs3F3Zz7iru/a2Z/SzBL6gSCqbevEwSVP9FN4HD3XWZ2G/AV4BZ339nFdS1mNgv4MsG027kErYG1BAv7/hN4NeLPlc3JQEPCQP1J4ePCPL+OlAAtABTJo4TFfdPd/Y0+ro5IQWiMQyRPzOxognGDPypoyGCmriqRXjKzLxKMa3yOIMdVpimwIoOGuqpEeilMIjgRaAKucfcFfVsjkcJS4BARkUg0xiEiIpEocIiISCQKHCIiEokCh4iIRKLAISIikfx/zp6ZVJBQKWoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(d_layer_df[d_layer_df[d_colname]>0].groupby(d_colname).apply(len), '.')\n",
    "plt.xlabel('Layer #', size=20)\n",
    "plt.ylabel('Number of nodes', size=20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## evaluate the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_layer_df['benchmark'] = [benchmark_dict[x] for x in  d_layer_df.index.values]\n",
    "sorted_d_layer_df = d_layer_df.sort_values([d_colname, 'benchmark'])\n",
    "proposed_id_order = sorted_d_layer_df.index.values\n",
    "proposed_id_order_dict = dict(zip(proposed_id_order, range(len(proposed_id_order))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "disconnected_neurons = set(benchmark_dict.keys()) - set(proposed_id_order_dict.keys())\n",
    "from_benchmark_disconnected = {x:benchmark_dict[x]+len(all_nodes) for x in disconnected_neurons}\n",
    "total_proposed_dict = proposed_id_order_dict | from_benchmark_disconnected\n",
    "\n",
    "\n",
    "eval_df = edge_list.copy(True)\n",
    "eval_df['pre_order'] = eval_df['pre'].apply(lambda x: total_proposed_dict[x])\n",
    "eval_df['post_order'] = eval_df['post'].apply(lambda x: total_proposed_dict[x])\n",
    "eval_df['is_forward'] = eval_df['pre_order'] - eval_df['post_order'] < 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "proposed solution: 29193811\n",
      "benchmark solution: 29023882\n",
      "Score: 0.6965\n",
      "Benchmark comparison: 169929\n",
      "\n",
      "\n",
      "Does better than benchmark :)\n"
     ]
    }
   ],
   "source": [
    "print('proposed solution:', eval_df[eval_df['is_forward']].weight.sum())\n",
    "print('benchmark solution:', edge_list[edge_list['is_forward']].weight.sum())\n",
    "score = eval_df[eval_df['is_forward']].weight.sum()/edge_list.weight.sum()\n",
    "benchmark_comparison = eval_df[eval_df['is_forward']].weight.sum() -edge_list[edge_list['is_forward']].weight.sum()\n",
    "\n",
    "print(f'Score: {score:.4f}')\n",
    "print('Benchmark comparison:', benchmark_comparison)\n",
    "print('\\n')\n",
    "if benchmark_comparison > 0:\n",
    "    print('Does better than benchmark :)')\n",
    "else:\n",
    "    print('Does worse than benchmark :(')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### save as file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "numbername_node_dict = dict(zip(range(len(all_nodes)), all_nodes))\n",
    "proposed_order = {numbername_node_dict[key]:val for key, val in total_proposed_dict.items()}\n",
    "\n",
    "pd.DataFrame.from_dict(proposed_order, orient='index', columns=['Order']).reset_index().sort_values('Order').reset_index(drop=True).reset_index().iloc[:,:2].rename(columns={'level_0':'Order', 'index':'Node ID'}).set_index('Node ID').to_csv('./data/output.csv')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}