submodular_function_greedy_maximization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Submodular Functions and Greedy Maximization\n",
    "\n",
    "This Jupyter Notebook is meant to complete [this](https://adrienvdb.com/2020/11/09/submodular-functions-and-greedy-maximization/) blog article. \n",
    "\n",
    "In this notebook, we will achieve the following:\n",
    "\n",
    "- Create some simple submodular functions\n",
    "- Implement both the *naive* and the *greedy* algorithms to maximize submodular functions under cardinality constraint\n",
    "- Re-create the examples mentioned in the blog post and run the algorithms on them\n",
    "- Observe that greedy does not always gives the best solution\n",
    "\n",
    "Let's start, first we import some useful Python packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import networkx as nx\n",
    "import random\n",
    "import time\n",
    "\n",
    "from collections import defaultdict\n",
    "from itertools import combinations, product\n",
    "from matplotlib import pyplot as plt\n",
    "from tqdm import trange\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discovering Submodular Functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Hiking Example\n",
    "\n",
    "In that example, we had the ground set being $N = \\{bottle, compass, sandwich, sunglasses, jacket, shoes\\}$. We could define $f$ to assign some predefined value for every subset of items as seen in the blog post. \n",
    "\n",
    "Here we will, give a value to each item, and we will construct $f$ as being the sum of the values of the items selected. That is, if for each $i \\in N$, we assign it value $v_i$, we have that\n",
    "\n",
    "$$f(S) = \\sum_{i \\in S} v_i.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The ground set together with a value for each item\n",
    "hiking_values = {'bottle': 4, 'compass': 1, 'sandwich': 4, 'sunglasses': 2, 'jacket': 3, 'shoes': 5}\n",
    "\n",
    "# The submodular function computing the sum of value of the items selected\n",
    "def f_hiking(items):\n",
    "    items_values = [hiking_values[item] for item in items]\n",
    "    return sum(items_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see the output for a few possible subsets of items."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Picking the sandwich and the bottle has value 8.\n",
      "Picking the sandwich and the compass has value 5.\n",
      "Picking the sandwich, the sunglasses, the bottle and the jacket has value 13.\n"
     ]
    }
   ],
   "source": [
    "arg1 = {'bottle', 'sandwich'}\n",
    "print('Picking the {} and the {} has value {}.'.format(*arg1, f_hiking(arg1)))\n",
    "\n",
    "arg2 = {'sandwich', 'compass'}\n",
    "print('Picking the {} and the {} has value {}.'.format(*arg2, f_hiking(arg2)))\n",
    "\n",
    "arg3 = {'bottle', 'sandwich', 'jacket', 'sunglasses'}\n",
    "print('Picking the {}, the {}, the {} and the {} has value {}.'.format(*arg3, f_hiking(arg3)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the items have different values, different combinations will obviously give different values. Note also that since all the values are positive here, the more items we select, the higher the value. So selecting all items in this case would give the maximum value $f$ can take. That is, $\\max_{S \\subseteq V}f(S)$ is achieved for $S=\\{bottle, compass, sandwich, sunglasses, jacket, shoes\\}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now as a little exercise, try to prove that $f$ submodular. Recall that a function $f$ is submodular if $f(A) + f(B) \\geq f(A \\cap B) + f(A \\cup B)$ is true for all $A, B \\subseteq N$. You can also verify it using Python if you prefer.\n",
    "\n",
    "\n",
    "<details>\n",
    "<summary><a class=\"btnfire small stroke\">>&nbsp;&nbsp;Show proof</a></summary>\n",
    "    \n",
    "Let $A, B \\subseteq N$. We will show that $f(A) + f(B) \\geq f(A \\cap B) + f(A \\cup B)$. We have\n",
    "\n",
    "$$f(A) + f(B) = \\sum_{i \\in A} v_i + \\sum_{i \\in B} v_i = \\sum_{i \\in A \\setminus B} v_i + \\sum_{i \\in B \\setminus A} v_i + 2\\sum_{i \\in A \\cap B} v_i= \\sum_{i \\in A \\cap B} v_i + \\sum_{i \\in A \\cup B} v_i = f(A \\cap B) + f(A \\cup B).$$\n",
    "\n",
    "In fact, we see here that there is equality, the function is called modular in that case.\n",
    "</details>  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Maximization of Submodular Functions under Cardinality Constraint\n",
    "\n",
    "We now turn to the maximization of submodular function, and more specifically the case where we have a cardinality constraint. That is, given a submodular function $f$ on ground set $N$ and a cardinality constraint $k$, we aim to find\n",
    "\n",
    "$$\\max_{S \\subseteq N: |S| \\leq k} f(S)$$\n",
    "\n",
    "and also (one of) the subset $S$ achieving this maximum value.\n",
    "\n",
    "### Naive Maximization\n",
    "\n",
    "Let's first implement the naive maximization. For that, we basically need to try all possible subsets $S$ of the ground set $N$, with here the restriction that the size of $S$ must be smaller than or equal to $k$. In Python, there is the very useful `combinations` function from the `itertools` package in order for us to get all possible combinations. Given a sequence `seq` and a number `size`, it returns all possible combinations of size `size` of the elements of `seq.`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [],
   "source": [
    "def naive_maximization_under_cc(submodular_function, ground_set, k=5):\n",
    "    \"\"\"Performs the naive maximization under cardinality constraint k given a\n",
    "    set function and its ground set\"\"\"\n",
    "    \n",
    "    # Variable that will hold the maximum value of f\n",
    "    max_value = -float('inf')\n",
    "    \n",
    "    # Variable that will contain (one of) the subset S achieving the maximum\n",
    "    max_subset = None\n",
    "    \n",
    "    # Try every possible size S can take, from 1 to k\n",
    "    for size in trange(1, k + 1, desc='Processing subsets of size'):\n",
    "        \n",
    "        # For each possible combinations of that size\n",
    "        for S in combinations(ground_set, size):\n",
    "            \n",
    "            # Compute the value of the set function for this subset\n",
    "            f_value = submodular_function(S)\n",
    "    \n",
    "            # If the value is bigger than the current maximum one, update the max\n",
    "            if f_value > max_value:\n",
    "                max_value = f_value\n",
    "                max_subset = S\n",
    "    \n",
    "    # Return the max value and the subset achieving it\n",
    "    return max_value, set(max_subset)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a mini-example, let's run the algorithm on our hiking example, with $k=6$ (so basically no cardinality constraint since $|N|=6$ in this case). We already know that the maximum value is achieved when we pick all the items, so let's verify it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Processing subsets of size: 100%|██████████| 6/6 [00:00<00:00, 5025.12it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best subset to pick is composed of {'bottle', 'compass', 'sunglasses', 'shoes', 'sandwich', 'jacket'} and has value 19.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "max_value_hiking, max_set_hiking = naive_maximization_under_cc(f_hiking, hiking_values.keys(), k=len(hiking_values))\n",
    "\n",
    "print('The best subset to pick is composed of {} and has value {}.'.format(set(max_set_hiking), max_value_hiking))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All good! Note how the naive algorithm is still fast in this case as we've seen in the blog post. Indeed, the size of the ground set is small here, so trying every possible combination is still not too costly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Greedy Maximization\n",
    "\n",
    "Let's now implement the greedy algorithm for maximizing a submodular function under cardinality constraint. We recall that the pseudocode is as follows:\n",
    "\n",
    "![title](https://adrienvdb.com/wp-content/uploads/2020/11/Screenshot-from-2020-11-05-11-11-41.png)\n",
    "\n",
    "So the idea is to start with an empty subset $S$, and iteratively add elements from the ground set which give highest marginal contribution at each step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 246,
   "metadata": {},
   "outputs": [],
   "source": [
    "def greedy_maximization_under_cc(submodular_function, ground_set, k=5):\n",
    "    \"\"\"Performs the greedy maximization under cardinality constraint k given a\n",
    "    set function and its ground set\"\"\"\n",
    "    \n",
    "    # Variable that will hold the maximum value greedy can find\n",
    "    greedy_value = 0\n",
    "    \n",
    "    # Variable that will contain the subset S found by greedy\n",
    "    greedy_subset = set()\n",
    "    \n",
    "    # k times, we add a new element to our current set S\n",
    "    for i in range(k):\n",
    "        \n",
    "        # Init some variables\n",
    "        best_contribution = -float('inf')\n",
    "        best_element = None\n",
    "        \n",
    "        # For every element not yet in our subset S\n",
    "        for element in set(ground_set) - greedy_subset:\n",
    "            \n",
    "            # Compute the marginal contribution of the element to S\n",
    "            marginal_contribution = submodular_function(greedy_subset | {element}) - greedy_value\n",
    "            \n",
    "            # If the marginal contribution is better then the current best, update\n",
    "            if marginal_contribution > best_contribution:\n",
    "                best_contribution = marginal_contribution\n",
    "                best_element = element\n",
    "        \n",
    "        # Add the new element to the set S\n",
    "        greedy_subset = greedy_subset | {best_element}\n",
    "        \n",
    "        # Update the current function value achieved by S\n",
    "        greedy_value += best_contribution\n",
    "    \n",
    "    # Return what greedy found\n",
    "    return greedy_value, greedy_subset"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we compute the result for our hiking example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The best subset to pick is composed of {'sunglasses', 'bottle', 'compass', 'shoes', 'sandwich', 'jacket'} and has value 19.\n"
     ]
    }
   ],
   "source": [
    "max_value_hiking, max_set_hiking = greedy_maximization_under_cc(f_hiking, hiking_values.keys(), k=len(hiking_values))\n",
    "\n",
    "print('The best subset to pick is composed of {} and has value {}.'.format(set(max_set_hiking), max_value_hiking))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in the case of *modular* functions (which is the case in the hiking example), the greedy algorithm actually returns the optimal value, as we can see here. This can be proved mathematically."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other Examples \n",
    "\n",
    "We revisit now some of the examples presented in the blog post, and see how the two algorithms perform on them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Graph cuts\n",
    "\n",
    "In order to work with graphs, we use the NetworkX package. Given a graph $G(V, E)$, the submodular function we will consider here is the one on ground set $N=V$, and which compute the cut size of the cut given by the selected subset of nodes of $G$. I.e., we have $f(S) = |\\{(u, v) \\in E: u \\in S, v \\in V \\setminus S\\}|$ for all $S \\subseteq V$.\n",
    "\n",
    "For the purpose of this example, we will work with a random graph on $|V|=20$ vertices, generated using the Barabasi-Albert procedure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/adirlou/.local/lib/python3.6/site-packages/networkx/drawing/nx_pylab.py:611: MatplotlibDeprecationWarning: isinstance(..., numbers.Number)\n",
      "  if cb.is_numlike(alpha):\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXdYVEf3x7/XCliRYqMpsQB2EMQGiWLBCmokGiyJ8Y3RN/r6Myr2GBtqEmMsMRqjiT2KRLFiiS1KsSuKYsGCFQGVvrvf3x8XEGGBXdhCuZ/nmYe9986dOfeye2bmzJkzAklISEhISJQtyulbAAkJCQkJ3SMpfwkJCYkyiKT8JSQkJMogkvKXkJCQKINIyl9CQkKiDCIpfwkJCYkyiKT8JSQkJMogkvKXkJCQKINIyl9CQkKiDFJB3wLkhampKW1sbPQthoSEhESJ4vz58y9JmhWUr9gqfxsbG4SHh+tbDAkJCYkShSAI0arkk8w+EhISEmUQSflLSEhIlEEk5S8hISFRBpGUv4SEhEQZRFL+EhISEmWQYuvtIyEhoQGePwc2bACuXAESEoAaNYAWLYCRIwGzAr0BJUoxkvKXkCiNhIUBCxcCBw6Ixykp764FBACzZwM9ewJ+fkDbtvqRUUKvSGYfCYnSxurVgLs7EBgoKv3sih8AkpPFc4GBYr7Vq/UhpYSekZS/hERpYvVqrBg/Hk5JSahMYkS2S5sBVM2WjEgISUk4P3Gi1ACUQSTlLyFRWggLAyZNQr30dMwA8FmOy0MBvM2WVgFoCKBNSgowaRIgragvU0jKX0KitLBwIZCcDG8A/QGYFJB9I4BhAARANAUtXKhlASWKE5Lyl5AoDTx/Lk7ukipljwZwEqLyByDet38/8OKFlgSUKG5Iyl9CojSwYYNa2f8A0AlAg+wnBUHtciRKLhpR/oIg9BAEIVIQhChBEKbmkedjQRAiBEG4LgjCFk3UKyEhkcGVK7m9evLhDwDDc55MTgauXtWkVBLFmCL7+QuCUB7ASgAeAB4BCBMEYQ/JiGx5GgHwA9CBZJwgCOZFrVdCQiIbCQkqZz0DIAbAQGUX4+I0JJBEcUcTPX9nAFEk75JMA7ANQL8ceb4AsJJkHACQfK6BeiUkJDKpUSProwxACgB5RkrJOJfJRgADAFRTVo6xsdZElCheaEL51wfwMNvxo4xz2WkMoLEgCGcEQTgnCEIPZQUJgjBaEIRwQRDCX0gTTxISKiN3cICsYkUAwDwAhgAWAdiU8XleRr4UADugxOQDAIaGQPPmWpdVonigqwnfCgAaAXAH8AmAtYIg1MyZieSvJJ1IOplJcUckJAokKSkJP//8M9quXAm5XA4AmAOAOdKcjPwGAOIBdFFWGAmMGKFdgSWKDZpQ/o8BWGY7tsg4l51HAPaQTCd5D8AtiI2BhIREIYiPj8f8+fPRsGFDHD9+HL8EBKByv36ix05hEATA01MK9laG0ITyDwPQSBCEBoIgVALgA2BPjjyBEHv9EATBFKIZ6K4G6paQKFM8efIEkydPhq2tLW7fvo3jx48jICAAzs7OYpA2Q8PCFWxoKN4vUWYosvInKQMwDsAhADcA7CB5XRCEuYIg9M3IdghArCAIEQCOA/iGZGxR65aQKCvcuXMHX375JRwcHJCSkoILFy5gw4YNsLOze5epbVtg6VLAyEi9wo2MxPucnDQrtESxRiMhnUnuB7A/x7lZ2T4TwMSMJKFPpPjuJYorV65g0aJFCA4OxpdffonIyEjkOx82Zoz4d9Ik0W8/vxW/giD2+JcufXefRNmBZLFMjo6OlNAgoaGklxdpYCAmUS2IydBQPOflJeaT0DunTp1ir169WLduXS5evJgJCQnqFRAWRnp7i/9XQ8P3/t9pFSuK5729xXwSpQoA4VRBxwpUMRaIrnFycmK4FGVQM6xeLfUESwAkceDAASxcuBAxMTGYPHkyhg8fDgMDg8IX+uKFONK7ehWIi0Pk8+e4pFBg8P790kivlCIIwnmSBdrwpJ28SjuZij8pqeC8pJhv0iTxWGoAdIJMJsPOnTuxaNEikMTUqVMxaNAgVKiggZ+nmRnwzTdZh2/CwzH/s88wWFL8ZR4psFtpJizs3cYewHsbe0QAcAJgnJG6ZpwD8K4BkEZeWiUlJQVr1qxB06ZNsXLlSixYsACXLl3CJ598ohnFr4TmzZsjKioKSap0BiRKNZLyL80sXJjnxh71AOwE8ArASwB9IfroZiHFd9cab968wZIlS9CwYUPs2bMHGzZswKlTp+Dp6QmhsH76KlK5cmXY2dnh8uXLWq1HovgjKf/SSkZ897w29qgJwAbiRh4EUB5AVPYMUnx3jfPixQvMnDkTDRs2xIULF3DgwAHs27cPHTt21KkcTk5OkObTJCTlX1pRMS57TYhL/v8LYFrOi1J8d43w4MEDjB8/Hk2aNMHz589x9uxZbN26FS1bttSLPI6Ojjh//rxe6pYoPkjKv7SiYnz3eAAJAFYAaJ3zohTfvUjcuHEDI0aMQOvWrVG5cmVcv34da9aswQcffKBXuaSevwQgefuUXtSI714FwJcAzCAu0X5vswUpvrvahIWFYeHChThz5gz++9//IioqCsbFKFRys2bNcPfuXSQmJqJKlSr6FkdCT0g9/9JKtvjuqqAAkITcEfmk+O6qQRJHjx5F165dMWDAALi7u+PevXuYMWNGsVL8AFCpUiU4ODjg0qVL+hZFQo9Iyr+00qIFYGCQ58YewQAuZpx7DTHuhjEAu+xlSPHdC0ShUGD37t1wcXHBuHHj4Ovrizt37uDrr7+GkboxdnSIZPqRkMw+pRQOHw7FjBmYB+DbbOc3AZgNwAHiJO8jiJt9OAM4CHHy910hpTC+u4ZiG6Wnp2Pz5s3w9/dHtWrV4Ofnh379+qFcuZLRn3J0dMTJkyf1LYaEPlElBoQ+khTbp/CEh4ezU6dOPFK9OhWC8H4cH1WTIIixX0oLGoptlJiYyJ9++omWlpbs0qULjxw5QoVCoaOH0BwXL16knZ2dvsVQj2fPSH9/cuhQsndv8a+/P/n8ub4lK1ZAxdg+elfyeSVJ+atPTEwMR44cyTp16nDt2rWUnT1LGhkVTvkbGZWeoF+rVonPU1BDKAhivlWrchXx6tUrfvfddzQ3N6e3tzdDS3gAvLS0NBoZGfHNmzf6FqVgpKCEaiEp/zJEUlIS582bRxMTE06ZMuX9CJCZik9dxa9EAZZIivj8MTEx/Oabb1irVi2OGDGCERERen4gzeHs7MyTJ0/qW4z80UDDXdaQlH8ZQKFQcNu2bbS2tuaAAQN4584d5RlV/AEpBIGJAGNmzdLtg2iL0NA8Ff9WgE0BGgFsCPBkjutyAwPO9/KisbExv/76a0ZHR+v7aTTOV199xR9++EHfYuRNWe+4FBJJ+ZdywsLC2KFDB7Zu3Zr//POPKjfkGd89a+js7c3AGTNoZ2dXMswBBeHlpbTBOwzQCuBZgHKAjzJS9jwygNeaNOHzUmxPXr9+PYcOHapvMZST0XD/DNARYCWAw/NQ+N9C3Kc+uDSaLAuBqspf8vYpYcTExGDatGk4fPgw5s2bh+HDh6N8+fIF3+jkBOzalSu+O4yNRXfOESMAMzP0A7D3yRN8/vnn2LZtm9YDjWmNjNhGYO79C2YDmAWgXcZxfSW3lwfgEB2tPfmKAY6OjvD399e3GMpZuBBITkY9ADMg7gObrCTbHQB/Aaib/WRmUMJdu7QvZ0lGlRaioASgB4BIiLHBpuaTbwDEVtqpoDKlnv/7JCUl8bvvvqOJiQn9/Pz4+vVrrdWVnJzMNm3a8Mcff9RaHVrH3z/35GBGj74iwIUAbQHWBzgWYJKyXqWhIbl4sb6fRGukp6ezSpUq6u8Spm2ePcv1v5ueR8+/O8B9AK2z9/wB8f5SPGrLD6jY8y+yU7IgCOUBrATQE4A9gE8EQbBXkq8agPEAQopaZ1mCJLZt24amTZviypUrCAsLw4IFC1CtWjWt1WlgYICdO3di4cKFOH36tNbq0Sp5xDZ6BiAdYjjrUwAuQVzsNk9ZGaU8tlGFChXQvHlzXLx4Ud+ivI+KwQT/AlAZgKeyi1JQwgLRhNnHGUAUybsAIAjCNgD9kG1vkAy+A+AP4BuUFYq4oCg0NBQTJkxAamoqNm3ahE6dOmlf5gwaNGiADRs2wMfHB+Hh4ahTp47O6tYIecQ2Msz4+1+8MxVMhKj85yu7oZTHNspc6evm5qZvUd6hQlDCNxCj0AbnlaGUN9yaQBPLEesDeJjt+BFymFEFQWgDwJLkvvwKEgRhtCAI4YIghL8oyXHkw8IAb2/A2hqYPRvYvBkIChL/zpkDWFmJ18PClN7+6NEjDBs2DF5eXhg9ejTCwsJ0qvgz6dmzJ0aNGoXBgwcjPT1d5/UXiTxiGxkDsIC4j0Em+c5qFLO4PJrGycmp+IV3ViEo4RwAvhD3pMiTUt5wFxWtr0UXBKEcgB8A/F9BeUn+StKJpJNZSd1jdPVqwN0dCAwUey85ezDJyeK5wEAx3+rVWZeSkpIwd+5ctGzZElZWVoiMjMSIESP0GjJg1qxZMDIygp+fn95kKBQZsY2UMRLAzwCeA4gD8COA3soyloHYRo6OjsUvxo8KQQmPAlgOoE5GegjgY4imhSxKecNdVDShVR4DsMx2bIH3g0NWA9AMwD+CINyH6GSxRxCEAneXL3Fk3yxdiZfJe5BZe+Vy9Wps2bIFTZs2RUREBC5cuIB58+ahatWqupE7H8qVK4dNmzZh165d2Llzp77FUZ18YhLNBNAWQGOIgexaA5iuLCNZ+mIb5aBp06aIiYlBghohwLVOtoY7r8CERwFcgzhncwnitqRrAIzNLKMMNNxFRpVZ4fwSxHmDuwAaAKgE4DIAh3zy/4PS6O2Tx4KiWID9MxYTWQHcrMRjIalcOfra2fHUqVP6foo8CQ8Pp6mpKW/cuKFvUVRG1rcv5eosEMqeSltso3xo3749jx07pm8xslA8fUpZxYokwNkZPvzZ02wl/y/J2+cd0JW3D0kZgHEQXXFvANhB8rogCHMFQehb1PJLDBl+yTkZC7FFfAZgM4AxAK7nyFOZxMamTXW+l6s6ODo6YuHChfD29sbbt2/1LU6BPH78GF/cvYs0VdZAKEFeqRJQ0kxdhaS4hHcmib1798LJ0xP/GBpCIQiYg9zaf46Se+8D6Jp5IAiAp6daUVrLIhoxJpPcT7IxSVuS8zPOzSK5R0led5L6/6ZpkjwWFCUC2AXRzakqgI4A+gL4M8ft5UgIBw4U+83SR40aBVdXV4waNSpzFFcsOXXqFNq2bQtbHx9UWr4cUDOuvrxyZUyrVAk/nTlTrJ9TU+h70pckgoKC0LZtW8yYMQPTp0+H28GDkBWy4YahYZlpuIuEKsMDfaQSZfbJY0HRBYCGOc4tAdi7BC8oSkpKYuvWrbls2TJ9i5ILhULB5cuX09zcnAcOHHh3oRDBwe7du8fmzZvz888/Z0pKiv4eSgdcv36dtra2Oq9XoVAwKCiITk5ObN68OXft2kW5XM7Xr1/T09OTPzZpQkXOUCQFJSm2jxTbR6cMHar0i3gSYO0c534F6JbXF9fXV99PohJ3796lubk5T58+rW9RskhMTKSvry9btmypPMCdirGNsseEefPmDfv3788OHTrw6dOnOnwa3SKTyVi1alW+evVKJ/VlV/rNmjXjzp07KZfLSZLR0dFs3rw5R48ezbS0NCmqZyGQlL8u6d1b6RdSWc9/aV49f0Asp4Swb98+1q9fn0+ePNG3KLx79y5btWrFIUOGMDExMf/Mz5+LIyxfX/F9+/qKx3lMDsrlcs6cOZNWVla8ePGiFqQvHnTs2JFHjhzRah0KhYL79u1j27Zt2axZM/71119ZSp8kQ0JCWK9ePX7//ffvb5BTiIa7LCMpf12SR8//LcQ4MreynfMFOKWE9/wzmTVrFt3c3Jienq43GQ4dOkRzc3MuW7ZMqztqbd++naampvzrr7+0Voc+mTBhAhctWqSVshUKBffv309nZ2c6ODhwx44d7yl9kvzrr79oamrKv//+O++C1Gy4yyqS8tcledj8CXAwQJ+MhuA0wOoAr5Vgm392ZDIZu3fvzkmTJum8boVCwQULFrBu3bqqhbTWAOfPn6elpSVnz56dS3mVdDZt2sRBgwZptEyFQsEDBw7QxcWF9vb23L59e673lvl/tLCw4IULFzRaf1lFUv66REkUwswUC7AfRD9/Syj38y/JfskvX76kjY0Nd+7cqbM6ExIS6OXlRRcXFz58+FBn9ZLkkydP6OrqygEDBvDt27c6rVub3Lhxgw0aNNBIWQqFggcPHmS7du3yVPokmZqayhEjRrBNmzZ89OiRRuqWkJS/7slj4xCVUglfUBQWFkYzMzPevHlT63XduHGDTZs25ejRo/XmhZOSksIRI0awZcuWvH//vl5k0DRyuZzVqlVjbGxsoctQKBQ8dOgQXV1daWdnx23btlEmkynN+/LlS7q5ubF///6lqhEtDkjKX9fks2VgQSmtUqUSP1n166+/0t7eXqs7gO3evZtmZmZcu3at1upQFYVCwR9++IF16tQp1iuz1aFz5848fPiw2vflVPpbt27NU+mTZGRkJBs1asRvvvmm1JnPigOS8tcHhdhzVG5gwDl16nDMmDGia1sJRaFQcOTIkfTx8dH4xKtMJuO0adNoaWnJkJAQjZZdVA4cOFBsGqSiMnHiRC5YsEDl/AqFgocPH2b79u3ZtGlTbtmyJV+lT5LHjx+nubl5qXhfxRVJ+esLFf2SZYC4gGXVKsbHx7NXr150d3fnixcv9P0EhSZzAdhPP/2ksTJjY2PZvXt3urm58dmzZxorV5PcvHmTjRs35tdff61Xz6eismXLFg4YMKDAfAqFgsHBwezQoYPKSp8U9ww2NzfXuktpWUdS/vpEBb/kUEtL/jx8eNYtMpmMU6ZMYYMGDXjlyhX9yV5E7ty5o7EFYJcuXWLDhg05ceLEYq9UX716xW7durFr165Fspvrk8jISFpbW+d5XaFQ8MiRI+zYsSObNGnCzZs3q6T05XI5p06dSltb2xIVGLCkIin/4kA+fskxMTE0MTHhrVu33rtl8+bNNDU1ZUBAgJ6ELjpBQUGsX79+kVbFbtq0iaamptyyZYsGJdMu6enpnDhxIj/44ANGREToWxy1kcvlrF69eq7Rp0Kh4NGjR9mpUyc2btyYmzZtUknpk+LK6wEDBrBjx44lelRbkpCUfwlg8eLF9PT0zHU+LCyMFhYWnDt3rlYXLmmTmTNn0t3dXe0ee1paGsePH09bW1tevnxZS9Jpl/Xr19PMzIz79u3Ttyhq4+7uzoMHD5IUlf6xY8fYuXNnNmrUiH/++ada/8+YmBg6OTnR19e31MdHKk5Iyr8EkJqayiZNmnDv3r25rsXExLBdu3YcOHBgiXSFk8lk7NatG7/55huV73n69Ck7d+7Mnj176izOjLY4c+YM69atyyVLlpSoBnzSpEmcN28ejx8/nqX0//jjD7Ub8UuXLtHKyorfffddiXr+0oCk/EsIhw4dYsOGDZmcnJzrWnJycon2J3/58iWtra25a9euAvOePXuWFhYWnDlzZqlx/4uOjmbr1q3p6+ur9P9bHJk1axZNTU35wQcfFErpk6LZz8zMjNu3b9eChBIFISn/EkT//v05b948pdey+5OfPHlSx5IVndDQ0AIXgK1Zs4ZmZmb5x3Upobx9+5aDBg2ii4sLY2Ji9C1Onvzzzz90d3enlZUVa9WqVSilr1AouGzZMtatW5dnz57VgpQSqiAp/xLE3bt3aWJiwgcPHuSZ5+DBgzQ3N+evv/6qQ8k0w5o1a+jg4JDLfJWcnMxRo0bR3t5eJ6uD9YVCoeDcuXNpYWHBsGK2mO/EiRP88MMPaWtryw0bNjAtLY01a9ZU2602PT2dX331FR0cHHjv3j3tCCuhEpLyL2HMnj2bH3/8cb55IiMj2aRJE44dO7ZELQhTKBQcMWIEP/nkkyz774MHD9i2bVsOHDiQr1+/1rOEuiEgIICmpqbcunWrvkXhyZMn+dFHH7Fhw4b8/fff3+vpf/TRR9y/f7/KZcXHx7N79+7s0aMHExIStCGuhBroVPkD6AEgEkAUgKlKrk8EEAHgCoCjAKwLKrOsKf+kpCRaW1vz6NGj+eaLj4+np6cnP/zwwxLlOpeUlMRWrVpx+fLlPH78OOvUqUN/f/8yNxl4+fJl2tjYcNq0aXqZ2zh16hS7dOnChg0bcv369Uo7Ed+OHctgDw8xVHnv3uJff3+lgQfv3btHe3t7jh07ttivxSgr6Ez5AygP4A6AhhD3Kr8MwD5Hng8BGGV8HgNge0HlljXlT5K7du2ig4NDgb16mUzGyZMns0GDBrx69aqOpCs6UVFRrFq1KmvVqsXg4GB9i6M3nj9/zk6dOrFv3746G/WcPn2aXbp0YYMGDfJU+gwNJb28KKtYkSnlylHppileXmI+kv/++y/r1q3L5cuX6+QZJFRDl8rfFcChbMd+APzyyd8awJmCyi2Lyl+hULBr164q74/7559/0tTUlIGBgVqWrOi8ffuWPj4+tLW1ZZ06dUr1toiqkJqayi+++IIODg7Kt50kxVDh/v4q9cDz4syZM+zatSsbNGjA3377Le+OhZrbJYZ+9lmJXctQ2tGl8h8IYF22Y18AK/LJvwLAjDyujQYQDiDcyspKm++n2BIREUFTU1OVlWNoaCgtLCyKtT91VFQUmzdvzuHDhzMpKYkzZswo1AKw0oZCoeDPP//M2rVr8/jx4+8uZPTAaWCQe58IJT1wZZw5c4YeHh60sbHhunXr8h9NFiIgYaIg8OH06Zp7GRIao1gqfwCfAjgHoHJB5ZbFnn8m//d//8eRI0eqnD8mJoYuLi4cNGhQsVsQtm/fPpqbm3PlypVZjVPmArDJkyfrWbriwZEjR2hubs5Vq1YVecPyf//9l926daONjQ3Xrl1bsGNAaCh/rliRjgArARyerY6zALsCNAZoCnAgwJjsMhgZlfhQ5KWRYmf2AdAVwA0A5qqUW5aVf0JCAuvWrctz586pfE9ycjKHDRvGVq1aMTo6WovSqYZcLufcuXNZr149pUHeXrx4QSsrqxIdw0iT3L59m3Pq1GFKhQrvKfmfAaWK+T0FvGoVz549y+7du9Pa2pq//vorU1NTVavYy4u7AO4G+GWOOvYD3AEwAWAiwJEAu5eiTYhKK7pU/hUA3AXQINuEr0OOPK0zJoUbqVpuWVb+JLlx40a2bdtWLY8QhULB77//nnXr1tXrBiPx8fHs27cv27dvn+/CpswFYJGRkTqUrpgSGiqG+M6h3PNSzNlTcrly7F2njnpKn8y1/ej0fOogwPMAq+Y8X0K3Hy0WaGBORxm6dvX0BHArQ8FPzzg3F0DfjM9HADwDcCkj7SmoTI0ofy29XF0gl8vp6urKdevWqX1v5gYj+lgQdu3aNTZq1Ihjx45VSRH98ssvbNasWbEzV+mcArYBzU8xywWBsv791a/T318t5f8jQJec5w0NxUi1EqqjgTmd/Cjbi7y0/HJ1xfnz51m7du1CBTm7efMmmzRpwnHjxulsQdiOHTtoamrKDRs2qHyPQqHg8OHDOWTIkGI7Ya11cvTA1VX+he6BDx2qch2XM2z/J5Vd9/XVznspjRRxTkcVVFX+FVDaWL0amDQJSE4WX2NOkpPFv4GBwKFDwNKlwJgxupVRRdq0aYP+/ftj9uzZWL58uVr3NmnSBOfOncMnn3yCHj16YMeOHTAxMdGKnDKZDNOnT8f27dtx6NAhtGnTRuV7BUHAqlWr4OrqipUrV2LcuHFakbFYs2FD0csQBLGcb75R/Z6EBJWyRQHoCeAnAJ2UXD+zbx+W9O8PY2Nj1KpVC8bGxnl+rlmzJipUKH1qRyUydVNSUsF5STHfpEnisRZ0VOn6LxSzl6sJ5s+fDzs7O3zxxRdo3ry5WvfWrFkTQUFB8PPzg7OzM/7++280a9ZMo/K9fPkSn3zyCUgiPDwcpqamapdhZGSEgIAAuLq6wtHREa6urhqVsdhz5QqQklK0MpKTgatXC8yWmJiI8+fPIyQkBM6XL8OtgPzRED01ZkJ041NGQ0dHDBs2DHFxcXj16hXi4uLw+PHjrM+Z6dWrV0hISECVKlXybCTya0Bq1KiBcuXKqfdeigthYVgxfjw2pKfjKoBPAGzIdnkdgEUAngLoCGA9gHrAOx3Vti3g5KRRkUqP8g8LAyZNwoqkJGwAlL7gJACTAOwAkA6gJYCTWny5msDExARz5szBf//7Xxw/fhyCIKh1f/ny5bF48WI0b94cH374IdatW4d+/fppRLbz589jwIAB8PHxwbx584rUo7O1tcVvv/2Gjz/+GOfPn4e5ublGZCyuvHnzBjdv3kRERATanjsHe00UGhf33qFcLkdERARCQkIQGhqKkJAQREVFoXnz5nBxcUHXDh2g2L0bitRUyADIM1IKRMXwDMBHAMYB+DKvOg0NUdfDA97e3iqJqFAo8Pr1a6UNQ+bnu3fvKr329u1b1KhRQ+1Go1atWqhataravx2NsnAh6qWnYwaAQwCSs136B8A0AMcBNAIwHqLuOpGZITkZWLgQ2LVLoyIJVGYaKQY4OTkxPDxc9Ru8vYHAQASQKId3L3hDtiyfApAB+BlALYgzz46AOGT28tL4y9UUcrkcjo6O8PPzw+DBgwtdTmhoKLy9vTFmzBhMmzatSD+GjRs3YtKkSVi9ejUGDhxY6HJyMn36dJw9exaHDx8uFeaB2NhY3LhxAxEREe/9jY2NRZMmTWBjY4Ovw8Lg/uiR0vtlGelbAI8ArIWomJW9mbdeXjj86adZyv78+fOoW7cuXFxc4OzsDBcXF7Rs2RKVKlUSb3j+HLC2xpyUFHybo6zZAAQAcwBUyVlP9gMDA+DBA8DMTJ3XUihkMhni4+PzbDTya1BSU1NRs2bNfBuJvK4ZGhoWreHIeM+Zo7sZEP+XGzIuT4Koq1ZmHMcAqA/R3GabWYYa71kQhPMkC+zJlg7ln+PlArlf8E0AzhnnqisrQ4df4sJw+vRpfPLJJ7hx4waqVq0IxmTWAAAgAElEQVRa6HJiYmLg5eWFBg0aYP369TAyMlLr/rS0NEycOBGHDx/G7t274eDgUGhZlCGXy9GjRw84Ojpi0aJFGi1bW5DEkydPshR7diWfkpICe3t72NnZwd7eHk2bNgVJXLhwAUFBQYiKisLPlpbwuXkTFdLTc5U9B1CqmOfkOJeUcf4va2t8+OGHGDRoENq1a4datWrlK3vchx+i+j//oHxhHryYd5qyk5aWplIjoewayUI1GsbGxqhcuTKweDEwe3a+yj8JwKqM48cALAAEAsgaoxsaAt9+q9KcjqrKv+R3rQCVJsxCAVhD/IH8CaAuxB/QgMwMhZkw0yEdO3aEm5sbFixYgAULFhS6nHr16uHEiRMYPXo0OnbsiMDAQFhZWal0b0xMDAYNGgRTU1OEhYWhRo0ahZYjL8qXL48tW7bAyckJ7dq1Q//+/TVeR2FRKBSIjo7O1Yu/ceMGKleunKXg7ezsMGDAANjZ2aFu3bpITk7G0aNHERQUhB9//BFGRkbo06cP/P390alTJ1SMixM7L3ko/zkqyGZQuTK+OHoUtU6eREBAAPbu3Ys+ffrA29sbHh4eMDAweC8/SSxfvhz7Ll3C/kqVgLQ09V+IoSHg56f+fXqgUqVKqF27NmrXrq32vcnJyfk2GlFRUXleq1ixIv4EMCCfOZ0eAHwgmtYaQfSRFyA2CNmEUGlORy1UcQnSR1LL1TOHy5oyt7X5AAFwNsBUgP8ArAIwogS5rD1+/JgmJia8fft2kctSKBRcsmSJygvCTp06xXr16nHu3Lk6CUV87tw5mpmZ8datW1qvKydpaWm8ceMGd+3axXnz5nHIkCFs3bo1jYyMaGlpyW7dunHChAlcs2YNT506xZcvX+Yq4+HDh1y9ejV79erFatWq0d3dnUuXLs17QVsBfv4FugXmWGkbHR3Nn376iW5ubqxRowY//vhjbt26lQkJCYyNjWXfvn3p5OQkBpUrRGyfwrohliUUCgXfvHnDxC5dCnSpXQHwA4DmABcArK7MrbZ3b5XqRZny8+/du0Dl/wPAigDTs53rDXBZtuPHjo48d+4cHz9+XGz3kfX392evXr00Vt7+/ftpZmaW52IyhULBFStW6CWC46pVq9i8eXOtLQBLSkrixYsXuWXLFs6YMYMDBgygnZ0dDQwM+MEHH7BPnz6cMmUKN2zYwNDQ0HzDL8vlcoaEhHDmzJls1aoVa9WqxaFDh3Lr1q35rtNISEjgkSNHuH7MGCaXL1845V9AjJ1nz55x7dq19PT0pJGREQ0NDenh4cFHjx69y6QD//MyixrrKQgwEqARwFeFXE+hqvIvHWYfFcwPLZScyzmFc+v5c0z+73/x4MEDxMXFoX79+rCysoKlpSWsrKyyUuZx9epKZw+0yoQJE/Dbb78hKCgIvXv3LnJ5PXv2xKlTp9C3b19cvnwZP/zwQ9ZEa3JyMsaMGYMLFy4gZO9eNDhxAtiyRfQPr1EDaNECGDlSa/MkX375Jc6ePYsvv/wSf/zxR6En3V6/fp3LFh8REYGYmBjY2tpmmWoGDRoEOzs7NG7cOJeZRBlv375FcHAwgoKCsG/fPtSqVQt9+vTB8uXL4erqmmvCOj09HVevXs3yvAkJCcGDBw/QunVruLi44NqIEWizZQvKJSfnUaMSjIzEtSr5eKqZm5vjs88+w6tXrxAWFgZfX188ePAA9vb2aNOmDby9vdG/f39Ytm0repXs3y+aQbPLYWgoqiBPT9HUUww944otLVoAu3ZBlpKi1KtKBnFy1wHAQ4ihjccDMM5ehqEhoKard0GUjgnfbBMqeXlHEIAdgOEQI8+FQLS1hQFoCiBZEOBvZITb/frBw8MDnTp1Akk8fPgQDx48yErZjytUqJCrQch+XL9+fVSsWFHTrwaHDh3C2LFjce3aNZWUlCrEx8fDx8cHMpkM27dvx5s3b+Dt7Y2epqb41sAAFYKDxYzZbZeZCqFnT1EhtG2rEVmyk5SUBFdXV/znP//BV199lW/eFy9e5LLFR0REID4+Hk2aNHlv4tXOzg62trZqexTdv38fQUFBCAoKwpkzZ9CuXTv06dMHvXr1gq1tlm8GSCI6OjpLyYeGhuLSpUuwsbHJ8rxxcXGBg4PD+9+RghYpZpYvCEgmIfP3R/XJkwt8L8OHD0d8fDy2bt0Ka2trAGLjHhwcnDVH8MEHH8Db2xsD3dxge+qUaGOOiwOMjUXFM2JEsXWIKM48v3YNxq1aYb5crnTyfgKAzhBj41QDMBLAPOD9SXjJ2ycPsnn7zEHe3hHXAYyCuJekNYD5ALwyMxkYIPrUKRy6cAHBwcE4duwYateuja5du8LDwwPu7u6oVq1aVpkkERcXl2/j8PTpU5ibmyttHDLPmZiYFKpH279/fzg7O2PatGlq35sXcrkcU6ZMwZYtW5CWloZt7u7ocuAAhAIUEQRBbAi0tFo6KioK7du3x549e+Di4oLHjx8j6t9/UX7TJlS6cQOyV68Qk5iIq4KA8GbNUK9ly/eUvJWVVaEXB8nlcpw7dy5L4T979gyenp7o3bs3unXrljX6i4+PR2hoaFavPjQ0FOXLl89S8s7OznByclJttBgerlIPfHmVKjiRmIidO3fm+R06efIkhgwZgqFDh2LevHl5dkbS09Nx4sQJ7N69G7t370atWrXg7e0NLy8vtGrVSr8+8iWU169f4/vvv8eKFStw3NgYze/ehVAYfaumV1XZUv5Alp9/vkoqL5S8XLlcjosXLyI4OBjBwcEIDQ1Fq1at4OHhAQ8PDzg7OxfYa5TJZIiJicm3gUhJScl39GBpaQlDQ8NcZd+9exfOzs64ePEiLC0t1X9mJZDEkiVLsGDBAoxMScESABVSU1UvINMEoYEGQC6X4/79+1m9+IMHD+LMmTNwKVcOU0l0SU+HIAioJJe/k9/QUPxxFXEkkpCQgMOHD2Pv3r04cOAA6tWrh969e6N3795wdnaGXC7HlStX3uvVP378GG3atMlS9i4uLqhfv37RlOaLF6IHWh498JSUFDg5OcHPzw9Dhw7N9f4WLlyIFStW4Pfff0fPnj1VrlahUCAkJAQBAQEICAiAQqGAt7c3vL294erqWnJX2eqI1NRUrFmzBgsWLEC3bt0wd+5c2Lx4Abi7qxZ9ICdGRsCJEyqb2sqe8g8L0+rLTUpKwqlTpxAcHIwjR47g/v37cHNzy2oMGjduXKgf+tu3b3M1DtkbiEePHqF69epKG4eDBw/ixYsXCAwMLPIP8s2bNxg5ciQePHiAvbNmwXTgQJTPpvhTAXwFMTzrK4iLTxZCjPnyHmp+UdPS0hAVFZXLffLWrVswMzN7z1RTfu1afBIejspk/j2oQoxEoqKiEBQUhL179yI0NBSdOnVC79694enpCblc/p6d/sqVK7C1tX1v8ZS9vb1eFqWdP38ePXv2xMWLF1G/fn0AwNOnT/Hpp58iPT0dW7ZsyTpfGEjiypUr2L17NwICAvDixQv0798fXl5e+PDDD7Vi1iypKBQKbNmyBTNnzoS9vT0WLlyIFi2yzTaqE34mk0J0qMqe8gd09nIB4NmzZzh69GjWyKBcuXLw8PBA165d0bVrV5hpyDaqUCjw/PlzpaOH+/fv4+LFiyhfvjwsLCwKPTl969YteHl5wdXVFStWrIDBkCG5RlGJAJYAGAHACsB+iEvQrwKwyV5YHkPUpKQk3Lx5M9fE671792BlZZXLHt+0adP3F7OtXg1OmgQh2//2UwBHM2SrA2AyRLNeFvn8b2UyGc6cOZOl8BMSEtC7d2+4ubmhSpUquHLlSpYZx8DA4D07vaOjY5EW2mmab7/9FufOncP+/ftx9OhRDBs2DKNGjcKsWbM03iDdvn07qyG4desWevfuDW9vb3Tr1k3tBYOlBZI4ePAg/Pz8YGBgAH9/f7i55RE1ScU5naKYUsum8gd08nJzQhKRkZFZo4ITJ06gQYMGWfMFnTp1Umq60QS7du3CrFmzxF7Z9esw2L4dlW7eBOPj8TI9HRdlMqxITER8xYpKzUsPHjzA999/j++++w5jxoxRulo6L1pAnE8ZkOO8olIlbFu8GBcePsxS9E+fPkWjRo3eU/D29vZo1KiRuAoyP/IY1V0H8AGAyhBXcLsD2IeMkB2ZZBuJxMXF4cCBAwgKCsLBgwdhbW0NJycn1KxZEzExMQgNDcWzZ8/g5OT0nrKvV69ege9Cn6Snp6Ndu3aoW7cuLly4gD/++ANdu3bVer2PHz9GYGAgAgICEB4ejq5du8Lb2xu9evVCzZo1tV5/cSAkJARTp07F06dPsWDBAvTv379gC4CKczqF9aoqu8of0PrLLYj09HSEhoZmjQouX76Mdu3aZTUGrVu31pjdlCTGubjgfykp+OD2bfFkDo8ckkjr0gXRn3yCyOrV8eDBA9y/fx8HDhzA7du3UbNmTcTFxcHc3Bx+5ctj1OPHqJzNlq6MZxAnzS9B9JbKToogYHfr1ojOcJ20t7dHgwYNCt8LVWE+JxKi8v8JwMfZzlMQcNvBAaOMjXHhwgU0btwYxsbGePXqFW7duoVGjRq9NylrZ2eH8uULFexAbzx+/Bj9+/fH5cuXcfr0aTg7O+tchtjYWOzduxcBAQH4559/0KFDB3h5eaFfv36FWlVb3ImMjMS0adMQEhKCOXPmYMSIEep/vwuY0yksqip/vS/myitpZCev58/FXYZ8fcWFYL6+4rGOd/JKSEjg33//zXHjxrFp06Y0MTHhxx9/zF9//ZX37t0rWuGrVlFuYECZKqtAMxbnvHr1ij179mTnzp359OlTkmR6ejofPHjAZ926FbioKA1gF4Cj88nz6KOPeOrUKV6+fJn37t1jbGxs4TaVKWCjkzEADUVPXrYG+EZJnmSANlWq0NLSkoMGDeKSJUt48uTJUrF72L59+1i7dm3OmzeP/v7+dHNz0/sCxdevX3P79u0cPHgwa9SowU6dOvHHH3/k/fv39SqXJnj06BG/+OILmpqa0t/fn0lJSfoWKRdQcZGXRnr+giD0gNjpKg9gHclFOa5XBvAHxBF5LIDBJO/nV2aRev7FnIcPH+LIkSNZZqIaNWpkzRd89NFHqg+ZCzHHoTAwwOyqVfH200+xePHi3BN2ffoAQUF53w9gCIDXAP4GkNd03zkzM/xfo0Z4/fp1Vnrz5g0qVqyI6tWrZ6Vq1arle9zm6FHYb9+O8vnEnpEDOAsxNO4UJTLJKlVC0pQpqD53br7vpiSRnp6O6dOnY+vWrdi8eTM6d+4MuVwONzc3DBw4EBMmTNC3iACAlJQUHD16FAEBAdizZw+sra3h5eUFb29v2NnZ6Vs8lYmLi4O/vz/Wrl2LUaNGYerUqTA2Ni74Rj2gM7OPIAjlIe7f6wFxXVUYgE9IRmTL8xWAFiS/FATBB4AXyXxjE5dm5Z8dhUKBq1evZjUEZ86cgYODQ1Zj4Orq+i4Eb3by8W5yB3AO76L21YdoFslEVqkSKpw5o9zk9emnwObNSmUlgM8A3Ic44ZvvLIavL/DHH+/fTyI5OTlXg6DsODY2Fi9fvsR/Tp9Gt+fP86spiy8B2AP4WkV5SirR0dHw8fGBsbExNm7c+J5zQVRUFNq1a4fTp0+jadOcBjn9IpPJcOrUKQQEBGD37t2oVq1algtpmzZtNL+W4Plz0axy5UqhV6UnJydjxYoVWLJkCfr27Ys5c+bAwsJCs3JqGJ2ZfQC4AjiU7dgPgF+OPIcAuGZ8rgDgJTIanrySRsw+JZDk5GQePXqUU6dOpaOjI6tVq0ZPT0/++OOPvHr16rt9bvMJBOYGcG1BJqAcgcCyyLGpd/b0H4gbeCszrbyXcmzqLZPJ+Pz5c16/fp0nTpzgrl27+Msvv3DevHkcP348hw4dyu7du7NNmza0srKikZERK1euzPr16/NY1aoFmqEy0+cAv87ruopBsYo7gYGBNDc35+LFi/M076xcuZLOzs5MT0/XsXSqkxkLacqUKWzUqBGtrKw4fvx4njhxgjKZrGiFa2AP7/T0dP7222+0tLSkl5cXIyIiiiaTDoGuzD6CIAwE0IPkqIxjXwAuJMdly3MtI8+jjOM7GXle5lVuWen5F0RsbCyOHTuWNXmcmpoKrw4dsPzvv1FeSQhgQOz5f4ocbo85yWu5eB7ePtEQXTor4/044GsAvL+8CEgrVw7eTk6ISkjAy5cvER8fj5o1a8LMzAympqYwNTXN+pz9HEk8evQIt2/fxuXLlxEeHo51KSnwVhLr5jmAYwB6QxyBHAHgDWArgL5KHvdyixZ4tWwZ2rVrpzXPK22SmpqKKVOmIDAwEFu3bs13q0uFQoHu3bvjww8/1OgKcG1BEhEREVmLymJiYtCvXz94eXnho48+KtgbLDtF9PYjiT179sDPzw9mZmZYtGhRidtWVJc9/4EQ7fyZx74AVuTIcw2ARbbjOwBMlZQ1GkA4gHArKyuttIolnaioKP7r5cWUcuXy7AG7ATQFaAKwPcDj+fTO09LS+OTJE169epXHjh3jjh07eKdlS8oLGV5YLgh85OLCkydPMiIigi9evFDak3v58iUPHjzI7777jn379mXdunVpYmLCHj16cObMmdyzZw+fPHmS50jkOcDOAGsArAawGcBf85ApvVIlBrm7s127dqxSpQo7derEGTNmMDg4mImJiXr4L6pHVFQUHR0d2a9fP8bGxqp0z4MHD2hmZsZLly5pWTrNc+fOHS5dupQdOnRgzZo1OWTIEO7cubPgCfoihqY+efIk27dvz+bNm3Pfvn3vRtklDOgqpDMks4/uUbJ/QfZ0DuBrgCkANwCsCjBKSb6tFSuyQoUKNDc3p729PTt37swBAwZwgZcXUytUKJTyVxZe+M2bNzxx4gSXLl3KwYMHs2HDhlkx7r/55hvu2LGD9+7dU/5jK8DbR6VkYJDl4fXmzRseOnSI06ZNY4cOHVilShW2b9+efn5+PHToEN+8eaOL/6DK7Nixg2ZmZvzpp5/UVkYbNmxgixYtmJqaqiXptE9MTAxXr15NDw8PVqtWjf369ePGjRtzh8kODSWNjPgzQEeAlfB+2OR7ED3CqmRLczM7LAYGnNCxI62trfnHH38U3eykZ3Sp/CsAuAugAYBKAC4DcMiRZyyAXzI++wDYUVC5kvLPByX7F+SXugNcruR8avfuebsFFrIXlbZ8OUNDQ7ly5UqOGDGCDg4ONDIyoouLC8eNG8eNGzcyIiJCPXdEDW90kp3ExEQeOXKEM2fOZKdOnVilShW6uLhwypQp3L9/PxMSEtT852iG5ORkjhkzhg0bNmRYPrH680OhULBv376cPn26hqXTD69eveIff/zB/v37s1q1avTw8ODq1asZExOT9R3ZBXA3wC/zUP7pSr4jMoC3W7RgSkqKvh9RI+hM+Yt1wROix88dANMzzs0F0DfjswGAvyCGrQ4F0LCgMiXlnw8F9Pxzph4Af1J2raDNIVTc4EMuCEypUIELrKxoaGjIFi1a8PPPP+cvv/zC8+fPF73nmdGr09RIJD+SkpJ4/Phxzp49m+7u7qxatSrbtm3LSZMmce/evYyPjy/as6jAzZs32bJlSw4aNKjI9T158oS1a9fmuXPnNCRd8eDt27fcuXMnhwwZwg+qV2dKju9ozg1T8lP+OUeHJR2dKn9tJEn550M+HjlxAA9CXNiUDnATxF2BInPmzeGRkydhYaS3NxUGBpTnqDMRYLIgMMTCgpv/9z+ePn1aazb0lGXLmKRu718DO04lJyfzxIkTnDt3Lrt06cKqVauyTZs2nDhxIv/+++98d+kqDH/++SdNTU35yy+/aMzmvGPHDjZp0qRYLkjSBOkLFlBWqZJKyr8ewPoARwB8UZjfQwlAUv6llLS0NO5esyZXTyczPQfoBNHOXwOia+ZhZb31ypXz7ek8evSIu3fv5rRp0+jh4cEPatTgd9Wr87ilJc+YmPCkjQ0Tv/1WZ72lUaNG8TcnJyr0vNVgamoqT58+zfnz52fZoVu1asXx48czICBA6X6+qpCYmMjPPvuMjRs31sokrY+PDydMmKDxcosFKuzh/QZgWEaH6CnAAQC7qTsSLiFIyr+UERsby0WLFtHCwoJubm587OJCRWE9cgDuqViRY8eOZUxMDGNjY1XzvMlg/fr19NXhD2Xz5s1s3LixuIduxkiEBgZiby3bc6VXrCie9/ZWy9RTFNLS0nj27FkuXLiQPXr0YPXq1dm8eXOOGzeOO3fu5HMVGsdr167R3t6en376qdYmnF++fMl69erx+PHjWilfr6iwh3fO9CRjJPA6+/lSshZEUv764Nkz0SQzdKj4RRo6VDwuQu84IiKCX375JWvWrMlhw4bx/Pnz4oUi2MFllStzzRdfsHHjxixXrhwrVqzIDh06cNKkSdy+fTvv3r2br8nh6NGj7NSpU6GfSR0iIyNpamqauzecI27TtTZtGNC+vd7ttunp6QwJCeHixYvp6enJGjVq0N7enl999RW3b9+eFUuJFCdkf/vtN5qamnL9+vVady0MCgqijY1NvhvRlyQUCgWjoqJ4t0MHtZX/0wzlHy/1/ItfKlHKXwMrCrOjUCh48OBB9ujRg+bm5pw1a9Z7Pe8sCuGR8xbgIhubLM+b48eP8z//+Q9NTEw4e/Zslbxb7ty5Q12sw0hKSmLLli25evXqAvOePn2axfE7I5PJGB4ezqVLl7JPnz6sWbMmmzZtypEjR7JDhw5s0qQJr127pjN5Pv/8c44ePVpn9WmSuLg4Hjp0iN9++y09PT1pampKCwsL/tGsGdMyXJPTIc53TQX4Kd7NfZ0DeDNj1PsS4McA3XP+TiWbf/FIxfGHrBQVPWJUsUMnJibyl19+oZ2dHVu0aMH169czOTlZpfoLMgHJBYEyAwOm//yz0mLu3LnDYcOG0czMjP7+/vlO3KamprJixYpaDx/wn//8hx9//LFKPeKkpCQaGRkV+0VbMpmM27Zto5mZGa2trWlsbMxGjRpx1KhR3LRpEx8+fKjV+hMSEmhtbc0DBw68O6mFEWtRSUtL4/nz57lq1SoOGzaMTZo0YdWqVenm5sbJkydz165dfPTo0Tv5MzpdszN69NnTbIBbANpAdH6oA9A3w/ST9RuRvH2KTyoRyr+IKwozefjwIadOnUpTU1P27duXx44dy1fhZQ53t27dyokTJ3JUq1YMLF+eSUDulb+ZIw8V7eAREREcOHAg69aty59//jlP32cLC4uih6POh23bttHW1lYtP/u2bdvy5MmTWpOpqCgUCq5evZqmpqbctGkTSTHGzeXLl7l8+XJ6e3vT1NSUtra2/Oyzz7hx40athEE+evQoLSwsmHDkiEZHrIVFoVAwOjqaO3bs4MSJE7MW3zk4OPCzzz7jmjVrePny5fw7G1pcC1LSkJS/tilgRWEqRI8C64zex/GcDUBYGM+dO0cfHx8aGxvz66+/5u3bt5VW9fjxYwYGBnL69Ons1q0ba9WqRQsLC3p5eXHBggUMDg7mvn376GZvr7H9Cy5cuEBPT09aWVlx3bp1uX54HTp04D///FOoV1cQt2/fpqmp6bv5DRX5+uuvubiYDt3j4+M5aNAgtmzZkpGRkXnmk8vlvHbtGlesWMFBgwbR3NycNjY2HD58OH///fcC52NUZZu7O1PKl9eL59Tr16957NgxLly4kP369WOdOnVYu3Zt9u3bl/Pnz+eRI0fUX1ynw7UgxR1J+WubAlYUpgL8EeCpjGFmduUvFwQeMzamjY0Nf/jhh/cW8sTGxvLQoUOcN28e+/Xrx3r16uXreZPJV199xQULFmj8Mc+cOUN3d3c2atSIW7ZsyVqZO3ToUG7YsEHj9SUnJ7N169ZcsWKF2vdu2bKFXl5eGpepqISFhbFhw4b86quvCjbj5UChUDAiIoKrV6+mj48P69SpQ0tLS/r6+nLdunW8ffu2+o3BqlVU5PCUKsyIVRVkMhmvXLnCtWvX8vPPP2ezZs2yQmr873//47Zt2/IO7aEuGhqJl3Qk5a9NlMSbyc+7oH7Onj9AWcWKfH3nDk+ePMnvv/+ePj4+tLW1ZbVq1ejm5qay5w0pepiYm5vzzp07WnlchULBI0eO0MXFhc2bN2dgYCD9/Pz47bffaryusWPHcuDAgYVSBvfu3WOdOnWKTUAuhULBZcuW0czMjDt27NBYmZGRkfz11185ZMgQ1q9fn/Xr1+eQIUO4Zs0aRkZG5v/8+fSQbwGsDHBoEXrIMTEx3L17N6dOnUp3d3dWq1aNjRs3pq+vL1esWMGwsDDtxhrS4BxcSUVS/tpEyQpbdZV/kiDQr0IFOjs7c+zYsdywYQOvX79eqKBShw8fZtu2bbXwoO+jUCi4Z88etmjRgjY2NuzWrZtGFe1ff/3Fhg0bFjqkgUKhYO3atYvFdoGxsbHs168fnZyctNYok+/mf9atW0dfX19aWlqyTp06HDx4MFevXs2IiIj3/0f52MY9AHbMT/nnsI0nJiby1KlTXLp0KQcOHEhLS0vWqlWLPXv25Jw5c3jgwAGVo5BqlHzWgqg7B1YSkZS/NlFhRWFByp8AZUOGaESczz77jN9//71GylIFuVxOPz8/GhkZ0c3NjadPny5ymXfu3KGZmVmhg5hl0r9/f27btq3I8hSFf//9l9bW1pwwYYLOg4UpFArevXuXv//+O4cPH04bGxuam5tz0KBB/G3hQspzhEHITFsBDoLoGZOn8oe4kO6bESPYunVrGhkZsW3bthw3bhz//PNP3rp1q9iMukgWmz28dY2k/LWJmisK81L+mlhRmJKSwlq1amndRTAnN2/eZMOGDbl+/XpaW1uzZ8+eak/QZpKSkkInJycuW7asyHItWrSI48ePL3I5hUEul9Pf35/m5ub8+++/9SKDMu7fv8+NGzdyu5OT0vhICQAbAXyogvJPKV+eJ/r04b///qv2/IWEbtIRTK4AACAASURBVFBV+ZcrcLcXidzUqKGZcjSwAfThw4fh4OCg831Frays8OjRIwwfPhyRkZHo1asX+vTpg4EDByIiIqLgArIxZcoUWFhY4Ouvle6+qxaurq44d+5ckctRlxcvXqB3794IDAxEaGgo+vZVtp+YfrC2tsawYcPwcZMmMCRzXZ8J4HMAqnyDKsvl6FyzJlxdXWFgYKBpUSV0iKT8C0OLFuI2iABkAFIAyDNSSsY5AEjNOAaAtIzPWT89Q0OgefMii7Jt2zb4+PgUuRx1MTQ0hLGxMZ48eYLKlStj7NixuH37NlxcXODu7o5hw4bhzp07BZYTGBiIv//+G+vXr9fIBt6Ojo64evUqUnJsQ6lNTp48iTZt2qBFixY4ceIErK2tdVa3WiQk5Dp1CeIWmP9Tp5y4OA0JJKFXVBke6CMVa7OPCisKiXc+/tnTvcwhtAZWFCYmJrJGjRp89uyZhh5MPVxcXJTa+xMSEjhnzhyamJhw9OjReZqk7t27R3Nzc43Hmm/Tpg3PnDmj0TKVIZPJ+N1337FOnTrvr5gtriiZq/oR4qrX2hmpCkADgK3zMf2Ulhg4pRVIZh8tYm4O9OwJCALmILeGn5OR7b6SazaAuHm0p2fuzdPVZN++fXB2doa5uXmRyiks1tbWiI6OznW+evXqmD17NiIjI1GzZk20bNkS//vf//D8+fOsPGlpafDx8cHkyZPh4uKiUbnatWunddPP06dP0b17dwQHByM8PBw9evTQan0aIduINZPREHdgupSRvgTQC+K+q0rR0IhVQv9Iyr+w+PmJP4TCYGgo3l9E9GXyySQv5Z+JiYkJ/P39cf36dcjlctjZ2WH69OmIi4vDtGnTYGZmhokTJ2pcLldXV5w9e1bj5WZy5MgRtGnTBu3bt8fRo0dRv359rdWlUUaMyHXKCECdbKkqxG338uyWkErLkSiBqDI80Ecq1mafTPS4ojAhIYHVq1fX+E5S6rBy5Uq1IkTev3+fn3/+OatVq8YaNWpoLTbQ7du3aWFhofFy09PTOWPGDNarV49HjhzRePk6QYqBU+qBLsw+giDUEgQhWBCE2xl/c7mvCILQShCEs4IgXBcE4YogCIOLUmexYswYYOlSwMhINOXkhyCI+ZYuFe8rInv27EHnzp1hrAGPocJSUM9fWf5Zs2ahYsWKcHZ2Rrt27fDDDz8gOTlZo3LZ2toiJSUFjx490liZjx8/RpcuXXDu3DlcuHABXbp00VjZOqUYjFgligdFNftMBXCUZCMARzOOc5IEYBhJBwA9ACwTBKFmEestPowZA5w4AXh5ifbUnD8sQ0PxvJeXmE8Dih/Qv8kHAGxsbHD//n2V86enp8PHxwdTpkzB4cOHERwcjJMnT6JRo0b45ZdfkJaWphG5BEHQqN3/wIEDcHR0RLdu3XDw4EHUrl1bI+XqhbZt33VY1CGz4+LkpB25JHSPKsODvBKASAB1Mz7XBRCpwj2XATQqKF+JMPvkREcrCmNjY1m9enW978j05s0bGhoaqryqc/LkyezZs2dWcLhMQkJC6OHhwQYNGnDjxo2FCnGRk/nz53PixIlFKiMtLY2TJ0+mhYUFT5w4UWSZihVSDJxSC3SxwhdAfLbPQvbjPPI7A7gBoFxBZZdI5a8j1q5dy4EDB+pbDJKkiYnJe1sT5sW+fftoYWHBFy9e5Jnnn3/+YYcOHWhnZ8e//vorVyOhDseOHWP79u0Lff/9+/fZrl079uzZM1+ZSzRlPAZOaUVjyh/iGpBrSlK/nMoeQFw+5dTNGCm0yyfPaADhAMJ1sU1gSaVLly7cuXOnvsUgKfrUh4SE5Jvn4cOHrF27tkobrSgUCh44cICOjo5s3bo19+3bV6h4MW/evKGRkVGhIkgGBgbS3NycixcvLlIDVGIoozFwSiu66vmrZPYBUB3ABQADVS1b6vkr58mTJ6xRowaTkpL0LQpJ0svLK99wxenp6ezYsSPnz5+vVrkKhYK7du2ivb0927dvz+PHj6stm7u9PaPHjVN5e8LU1FSOHz+e1tbW/Pfff9WuT0KiOKCq8i/qhO8eAMMzPg8H8HfODIIgVAKwG8AfJHcWsb4yz86dO9GnTx8YFtZjQ8MUNOk7e/ZsGBkZYepUZb4AeSMIAry9vXHlyhWMGTMGo0aNgoeHB0JCQgq+OSwM8PbG4chI1F2zBti8GQgKEv/OmQNYWQHe3mK+DO7cuYMOHTogOjoaFy9ehKurq1rySkiUOFRpIfJKAEwgevnchmgeqpVx3gnAuozPnwJIx7tFhJcAtCqobKnnr5wOHTowKChI32JksWzZMo4dO1bptUOHDrF+/foaCT+RlpbGNWvW0MLCgn369OHly5eVZyzEROaOHTtoZmbGn376qXiFJJaQKASQQjqXPqKjo1mrVi3t7oSkJoGBgezVq1eu848fP2adOnUKZa7Jj+TkZP7444+sXbs2Bw8ezJs3b767WIhFdykVKnCGyf+3d+fhUVTp4se/L4SQsIgQUEEEFUQWWYSgXlxAwRtFnRgU1BFHFPHeKAiOwAjooAIjCjqKIIqAUe8IKAKCv8QAIQPquIRFWYWgCKJIcCAIk2ZJcn5/nAo2SXfS0GvS7+d5+kl31amqt4rmrepTp85J8HscAaUiha/JX7t3qETee+89+vTpQ2xsbLhDOcHTg15FRUXcfffdpKam0qNHj4BuLy4ujmHDhrF9+3Y6duzIVVddxf3338+exYth+HAoKDhRtk6pV3VgSKn11Sws5JmCArT1uoo2mvwrkXnz5nHHHZH1gHRJnb+94LCeeeYZqlWrxpgxY4K23Tp16jBq1Chyc3M599xzWXP77RS7JX6Aw26vX4B4oK+HdcmRI/Dss0GLValIpMm/kti+fTu7du0K+JW0v84880yqVavGAaeP96ysLN544w3+8Y9/UL169ZBsf9yQIdxUrVq5X+YPgLOAqz3NNAbS02HfvqDEqFQk0uRfScybN4++ffsSExMT7lDKOP/889m5cye//PIL99xzD++88w7nnHNO6AJIS6twIJi3gD9hn0T0SATS0gIbl1IRTJN/JREJffl407x5c7777jvuvvtuBg0aFPpOz9avh3JG7toJrOT3NskeuVywYUOAA1MqcmnyrwQ2btxIfn4+3bp1C3coHjVv3pxZs2ZRXFzMX//619AH4GF4QnfvAFcBF1S0Hh2eUEWRyKtDUGWU3OitVi0yz9WFhYWsWrWK3NzckNTzl/afGjWoXc78t/Hc3WwZYeweW6lQ0+Qf4YwxzJ07lzlz5oQ7FI/y8vJ477336NSpE02aNAnJNl0uF5988gmZmZlkZmZy2/ffM7p6dWoWFZUp+y/gJzy38jmJDk+ooowm/wi3du3akofewh1KGcXFxdxzzz0kJyezbt26oG3HGMOWLVtOJPvPPvuMDh06kJSUxKxZs0hs1ozqF14IHpL/W0AfoG7FG9HhCVVU0eQf4Upu9FbUmiUcJk6ciMvlYsKECbRr1y6g6z5w4ABZWVknEr6IkJSUxKBBg5g7dy5nnllqPKAbb4RFi2wSd/O6LxsTgd69oZHXkWuVqnI0+Uew4uJi5s2bR3p6eug3npdnmz6uX29vqNarBx06wH33QaNGfPLJJ0yZMoXVq1dzzjnncPToUX777TfOOOOM09pcUVEROTk5J5L9hg0buOqqq0hKSuLPf/4zrVu3Lv8EOGoUZGae9ISvz3R4QhWNfOkDIhwv7dvHmE8//dS0bds2tBv96is7yHdcnH15GODjSO/e5qazzjLp6eknFmvTpo1Zv379KW1q9+7dZtasWaZfv36mQYMG5pJLLjGPPfaYWbp0qXG5XKce+6uvmuJT7NtHR6lSVQ0+9u2jV/4RLORt+6dPt/3juFxlqk8AOx2okZ7OgpgYYt26ci7p46d9OTdNjxw5wqpVq05c3e/Zs4devXqRlJTECy+8QNOmTf2LPzWVnJwc2qelEQeIp30oIWKv+CdPDti4ykpVJpr8I1RRURHz589n5cqVodlgSeJ3qzb5AXgI+ByoCdwOvIT90sQWFtryAKmpHvv1N8bw7bffnkj2n376Ke3bt+eGG26wN2oTEwPaNPTQoUOkZGaS/sYbdExPt102iJw4aQE24Rtj6/hHjdIByVXU0uQfoVauXEmTJk1o1apV8DeWk1Mm8YNN/GcBe4B84HrgVeCRkgIFBXa5rl1PXPnn5+eTlZXFxx9/TGZmJgBJSUkMHDiQd999l/pBbEv/zDPPcP3119Nx4EAYOND21ZOWZp/cPXDAtuNv39626tGbuyrKafKPUCGt8nn22ZOvjh07gMFAHHAOcAOwqVQZ43JxYMQI1jRsSHZ2Nq+99tqp3agNkM2bN5OWlsbGjRt/n9ioEYwYEfRtK1UZafKPQMeOHWPBggWsWbMm+BvLy4OMDI91/MOAuUAP4ACQAYwrVUaMofbKldRPTiYhIYGvv/465ENMGmMYMmQITz75JGeffXZIt61UZRWZ/QVEueXLl3PxxRfTvHnz4G+snJ4sr8Fe6Z8BNMWOzXmrh3I14+KYfMkl5Ofnh2Vs4fnz57Nv3z4eeuihkG9bqcrKr+QvIg1EZJmI5Dp/vVboisgZIrJbRKb6s81oENIqHy89YhZjq3n6AP8BfsVe/f/F0zpcLur+8AMHDx6k4HTa2fvh8OHDPPbYY0ydOjUiu7tWKlL5e+X/OJBljLkIO5B7ef1njQNW+bm9Ks/lcrFkyRL69q2wN5rA8NIj5n5gF7bOvyaQANwHeHvcTPLzOe+889i1a1cwovRqwoQJXHPNNVxzzTUh3a5SlZ2/yT8Z230Kzl9PtQKISBfgbGCpn9ur8jIyMujcuXPoBkOpV8/j5IbYLpCnA4XY1j5vAR28rad+fY/NPYNp69atvPHGG0yaNClk21SqqvA3+Z9tjNnjvP8Fm+BPIiLVgBeA4X5uKyrMnTs3tOP0duiAiYvzOGsB8DHQCGgJ1AD+7qmg0yOmp8Hcg8UYwyOPPMLo0aNp3LhxSLapVFVSYSWpiCzHtvQr7aTRuY0xRkQ8PVL5EJBujNldUZM/EXkQeBCgWbNmFYVW5Rw6dIjMzExeffXVkGyvuLiYhXXrctPRo3hK/52Af/qyIqdHzPNffz1kyX/RokXs3r2bIUOGhGR7SlU1FSZ/Y0wvb/NEZK+INDbG7BGRxkCeh2L/BVwtIg8BdYBYETlsjClzf8AYMwOYAZCYmFjOs/lV05IlS7jyyitp2LBhULdjjGHx4sU88cQT1KlTh6u7dSPuX//y3KVDRdx6xGzevDkZGRmBD7iUgoICHn30Ud58801q1KgR9O0pVRX52zxiMXZo1InO3w9LFzDG3F3yXkQGAImeEr+yI3YFu5XPihUrGD16NAUFBfztb3/j5ptvRlavhh49/O4Rs2Qg92B79tlnueKKK7j22muDvi2lqixfen/z9sI2AskCcoHlQANneiIw00P5AcBUX9Ydbb167t+/35xxxhkmPz8/KOv/4osvTM+ePU3Lli3Nu+++a4qKik4u8OqrtodLP3rE3Llzp2nSpElQ4i+Rm5trEhISzI8//hjU7ShVWeFjr55h77rZ2yvakv/s2bNNSkpKwNe7fv16k5ycbJo2bWpmzJhhjh075r1wyQlApPykL+KxK+Tjx4+b2NhYc+TIkYDvR4mbbrrJPPfcc0Fbv1KVna/JX5/wjRCBfrBr+/bt9O/fn169etG9e3dyc3MZNGhQ+XXkqamwciWkpEBcnK3ScRcfb6enpNhypbpCjomJoUmTJvz4448B2w93S5YsYfv27QwbNiwo61cqmugjkREgLy+PL7/8kgULFvi9rp9++olx48Yxf/58hg4dyvTp06lbt8IRbH+XmAgffHDaPWKWNPds2bKl3/vizuVyMXToUF5//XViY2MDum6lopEm/wjwwQcf0Lt3b2rXrn3a6/j111+ZOHEis2fPZtCgQWzdupWEhITTD+o0e8QM1k3fSZMm0blzZ66//vqAr1upaKTJPwLMnTuXxx577LSW/e2333jxxRd55ZVXuPPOO9m4cSNNmjQJcIS+a968ecCf8t2xYwdTpkxh7dq1AV2vUtFM6/zDbPfu3WzYsIGkpKRTWs7lcjF58mRatmzJjh07yMnJYdq0aWFN/BCcK/9HH32URx99NCof/FMqWPTKP5Ty8mw9+vr1tkO1evX4/uBB+iclUbNmTZ9WcezYMWbPns24ceO44ooryM7Opl27dsGN+xQE+so/IyODTZs2MW/evICtUymlyT80cnLsaFklT7+6daF8mQhXxsRAnz72YamuXT2uoqioiDlz5jB27FhatmzJokWL6OqlbNjk5dFx6VKGrVkDt9xiO43r0AHuu++0hk08evQojzzyCFOmTPH55KiU8pEv7UHD8aoy7fz9bDtfXFxsFi5caNq1a2e6detmsrOzw7Mf5fnqK2NSUoyJizPFcXEn71d8vDFxcXb+V1+d0monTJhgkpOTgxS0UlUT+pBXBPDjqdni4mKzbNky07VrV9OxY0fz0UcfmeLi4nDvUVl+nty82blzp0lISDDff/99kHdAqapFk3+4ffWVx8S/Gcy1YM4A0wLMAg+JsjAuzvxPly7moosuMnPnzi3bFUOkCECXEN7cdttt5qmnngrBTihVtfia/LW1T7A8+yy4XCdNKsSOfnMzdqSsGUB/YFvpZY8c4S/GsHnzZu644w6qVYvAf6acHBg+nKkFBSRiR/saUKrIe0AboC7QFlgEtvO44cNh9Wqvq162bBnr1q1j5MiRQQldKaVNPYMjL8/e3DUnd5H8LfAz8ChQHbgOuBJ4p9Ti1YELNm8m5sCBEAR7mpyTWxPgCeD+UrN/wp7YXgR+AyYBf8Tp89vlsst7cOzYMYYMGcJLL70UlsHglYoWmvyDIS3N56IG2OhphsgprSek3E5ufbBjd5Z+lng3cCZwIyDATUBt4DuwJ8X0dNuFRCkvvfQSLVq04Oabbw7mHigV9TT5B8P69Sc15yxxMXAW9ir4OHZA45WAx170XS7br04k8uGklIit8lkMFGGrfGriNgawh5Pb7t27ef7553n55ZepaNQ3pZR/tJ1/MBw86HFyDWwSHAI8h02Q/bBJ0aNIrfbxcnJzVx34E7aq5wgQC7yPvfoHPJ7cRowYQWpqasA7hVNKlaXJPxjq1fM6qwP2ar9EN+wQaB7Vrx+4mALJy8nN3XJgJHYM4M7AGuAPQAZ2bGDgpJNbdnY2n3/+ObNmzQpsrEopj7TaJxg6dLD93nuwHnslXABMBvZQtpUMYPvOb98+SAH6qZyTW4mvgWuwv26qAV2By7EnhROck9vx48cZPHgwL774IrVq1Qp0tEopDzT5B8OAAV5nvQM0xtb9ZwHL8FLtY0y56wkrt5NbIfZkVuS8jjjTugKfYE8CAOuczyfq/N1OblOnTqVp06akpKSEaAeUUn49iAU0wOavXOdvfS/lmmHvb24BNgPnV7TuSv+QV0pKxU+9enkVgtnWvn35Qy6G0969tssGMGNtg6WTXmOd/XjFeZCtDpgLwEx238+4OGPy8szPP/9sGjZsaLZu3RruvVKqSiBED3k9DmQZYy7CXsg+7qXc28AkY0wb4DKc5t5V2qhRZYdB9JHExTGldm06duzI0qVLAxxYAJx1Ftx4I4jwFGWz/1NOscHAduAQ8D1wYsQCEejdGxo1YuTIkTzwwAO0atUqlHuglPLlDOHtBWwFGjvvGwNbPZRpC3x6quuu9Ff+xvjdt8/ixYtNixYtTHJysvnuu+/CvTcn89J9hc/7mJNjVq1aZZo2bWoOHToU7r1RqsogRFf+Zxtj9jjvfwHO9lCmFZAvIgtEZJ2ITBKR6n5ut3JITYXJk6FWLXu1Wx4RW27yZEhNRUS45ZZb2LhxI5dffjmXXXYZY8aM4fDhw6GJvSJdu/6+b6fC2cfCTp0YPHgwL7zwAnXq1AlOjEoprypM/iKyXEQ2englu5dzzjjGwypigKuB4dj7gBfipYGLiDwoIqtFZPU+D09/VkqpqbByJaSk2JukpauC4uPt9JQUWy419aTZcXFxjBo1im+++YadO3fSpk0b5syZU/KrKrz8OLlNnz6dhg0b0rdv39DEqpQ6ifiTRERkK9DDGLNHRBoD/zTGXFyqzBXAc8aY7s7ne4ArjDEPl7fuxMREs7qczr8qpX377FOtGzbYNu7169sWLwMG+DzYyWeffcaQIUOoXbs2U6ZM4dJLLw1qyD5Zvdr21ZOeTmFxMTHHjv0+Lz7eVvb07m3vgyQmsnfvXi655BJWrlxJ27Ztwxe3UlWQiKwxxiRWWM7P5D8J+LcxZqKIPA40MMaMLFWmOrAW6GWM2Scib2LrpKaVt+4qmfwDpKioiNmzZ/Pkk09y6623Mn78eBo2bBjusGDfPpbefTeN9uzh0vPP93pyu++++0hISGDy5MlhC1WpqsrX5O9vnf9E4HoRyQV6OZ8RkUQRmQlgjCnCVvlkicgGbD9fb/i53ahWvXp1Bg0axJYtW6hZsyZt27bllVdeobCwMLyBNWpEWsOGbBw5EpYsgbffhhEjTkr8n3/+OUuXLmXs2LFhDFQp5VfyN8b82xjT0xhzkTGmlzFmvzN9tTHmAbdyy4wxHYwx7Y0xA4wxx7yvVfmqfv36vPzyy2RnZ/Phhx9y6aWXsmLFirDGtHXrVq/NNouKihg8eDCTJk2ibt26IY5MKeVOn/CtAtq1a8eyZct4+umnGThwILfffjs//PBDyOMwxrBt2zavyX/GjBnUqVOHu+66K8SRKaVK0+RfRYgIffr0YfPmzXTs2JEuXbrw1FNPUVDgscPooNizZw/x8fHU99Ah3a+//srYsWOZOnWqdtesVATQ5F/FxMfH8+STT7Ju3Tq2bNlCmzZteP/990PSNHTbtm1cfPHFHueNHj2aP/7xj7SP1M7qlIoymvyrqGbNmjFv3jzefvttxo8fz3XXXceGIA8O462+Pycnh48++oinn346qNtXSvlOk38V1717d9asWUO/fv3o2bMngwcPZv/+/UHZlqcr/+LiYh5++GEmTpxIPR+6glZKhYYm/ygQExNDamoqW7ZswRhDmzZtmD59OkVFRQHdjqcr/1mzZhEbG8s999wT0G0ppfyjyT+KJCQkMG3aNJYtW8a8efPo0qULq1atCtj6S1/579+/nyeeeEJv8ioVgTT5R6EOHTqQnZ3N6NGj6d+/P3fddRc//vijX+s8duwYu3btokWLFiemPfHEE/Tt25dOnTqVs6RSKhw0+UcpEaFfv358++23tGrVik6dOjF+/HiOVDAwuzc7duygadOmxMbGArB27VoWLlzIuHHjAhm2UipANPlHuVq1avH000+zevVq1q1bR9u2bVm4cKFvTUPz8uD556F/f+r1788Mlwuef57ivXt5+OGHmTBhgsc2/0qp8POrY7dg0o7dwiMrK4uhQ4fSuHFjXn75Zc+9bubk2F48MzLsZ/dfC/HxFBYW8kmdOnTPyKDa5ZeHJnClFBC6jt1UFdOzZ0/WrVvHLbfcQvfu3Rk2bBj5+fm/F5g+HXr0gEWLbNIvXU3kchFz/Dg98vOpdt11trxSKuJo8ldl1KhRg0ceeYTNmzfjcrlo3bo1M2fOpHjaNKYOHUpiQQE1jSkzIk8W0BqoBVxnDDsLCmD4cD0BKBWBtNpHVWjt2rVMGzCAqZs2kVFcTDUgE3ABaU6ZX4EWwEzgFuBJ4BPgC7AjeK1cCYkV/hJVSvlJq31UwHTu3JmZLVpQ0xj6ALcCCaXKLADaAX2BOOAp4BvgWwCXy94jUEpFDE3+qmJ5ecjHH1OtnF+Jm4CObp9rY38JbAI7jGN6uh3GUikVETT5q4qlpVVY5DBQuueeesChkg8iPq1HKRUamvxVxdavL9uqp5Q6wG+lpv0GnBivy+WyA9crpSKCX8lfRBqIyDIRyXX+enyiR0SeF5FNIrJFRKaIdvRSuRw8WGGRdtg6/hL/Ab5zpp9w4EBAw1JKnT5/r/wfB7KMMRdhW/o9XrqAiHQDrgQ6AJcAXYHufm5XhZJbV8yFwBGgyHkdcaalABuBD5xpz2D/wVu7r0ef9lUqYvib/JOBt5z3b2EbgpRmsA1AYoGaQA1gr5/bVaHUoQPExQEwHogHJgL/57wfDzTCJv4xQH3gS2Cu+zri40FH8VIqYvjVzl9E8o0xZzrvBThQ8rlUucnAA4AAU40xY7ys70HgQYBmzZp12blz52nHpgIoLw+aN6+w3r9ccXGwaxc0ahS4uJRSZQSsnb+ILBeRjR5eye7ljD2LlDmTiEhLoA3QFDgXuE5Erva0LWPMDGNMojEmsZEmichx1llw4422xc7pEIHevTXxKxVBYioqYIzp5W2eiOwVkcbGmD0i0hjI81AsBfjCGHPYWSYD+C/sA6Cqshg1CjIzoaDg1JeNj7fLK6Uihr91/ouBe5339wIfeiizC+guIjEiUgN7s3eLn9tVoda1K0yebLtqOBW1atnltGsHpSKKv8l/InC9iOQCvZzPiEiiiMx0yszHtvrbgG0N+I0xZomf21XhkJr6+wmgoiogkd8Tf2pqaOJTSvlMO3ZTp271attXT3q6TfIu1+/z4uNtdw69e9uqHr3iVyqkfL3hW2Gdv1JlJCbCBx/YvnrS0uyTuwcO2Hb87dvDgAF6c1epCKfJX52+Ro1gxIhwR6GUOg3at49SSkUhTf5KKRWFNPkrpVQU0uSvlFJRKGKbeorIPiBQnfs0xA4zG+kqS5xQeWLVOAOvssRaWeKEwMba3BhTYXO7iE3+gSQiq31p9xpulSVOqDyxapyBV1lirSxxQnhi1WofpZSKQpr8lVIqCkVL8p8R7gB8VFnihMoTq8YZeJUl1soSJ4Qh1qio81dKKXWyaLnyV0op5abKJH8RaSAiy0Qk1/lbZrRwEblWRL52ex0RukgRaAAABTBJREFUkVudeWkissNtXqdwxemUK3KLZbHb9AtE5EsR2S4i80QkNlxxikgnEflcRDaJyHoRucNtXtCPp4jcICJbnWPxuIf5NZ1jtN05Zue7zRvlTN8qIkmBju0U4/yziGx2jmGWiDR3m+fxexCmOAeIyD63eB5wm3ev813JFZF7Sy8bhlj/7hbnNhHJd5sXymM6W0TyRGSjl/kiIlOc/VgvIp3d5gX3mBpjqsQLeB543Hn/OPBcBeUbAPuBWs7nNOD2SIkTOOxl+nvAnc7714DUcMUJtAIuct43AfYAZ4bieALVseNEXAjEYseKaFuqzEPAa877O4F5zvu2TvmawAXOeqqHMc5r3b6HqSVxlvc9CFOcA7BjcJdetgHwvfO3vvO+fjhjLVV+CDA71MfU2dY1QGdgo5f5vYEM7PjmVwBfhuqYVpkrfyAZeMt5/xZwawXlbwcyjDGnMS6hX041zhNERIDrsAPknPLyp6jCOI0x24wxuc77n7HDeIaqL+fLgO3GmO+NMceAudiY3bnvw3ygp3MMk4G5xpijxpgdwHZnfWGJ0xiT7fY9/AI73nWo+XI8vUkClhlj9htjDgDLgBuCFCeceqx3AXOCGI9XxphV2ItMb5KBt431BXCm2CFxg35Mq1LyP9sYs8d5/wtwdgXl76TsF2KC89Pr7yJSM+ARWr7GGSciq0Xki5KqKSAByDfGFDqfdwPnhjlOAETkMuxV2Hduk4N5PM8FfnT77OlYnCjjHLOD2GPoy7KhjNPdQOyVYAlP34Ng8DXO25x/0/kict4pLhsoPm/PqUK7AFjhNjlUx9QX3vYl6Me0UvXnLyLLgXM8zBrj/sEYY0TEazMm58zaHsh0mzwKm+Risc2u/gI8E8Y4mxtjfhKRC4EVIrIBm7wCJsDH8x3gXmNMsTM5YMczWohIfyARO851iTLfA2PMd57XEHRLgDnGmKMi8j/YX1XXhSkWX90JzDfGFLlNi6RjGjaVKvkbY3p5mycie0WksTFmj5OM8spZVT9goTHmuNu6S65yj4rIm8DwcMZpjPnJ+fu9iPwTuBT4APuzMMa5km0K/BTOOEXkDOD/AWOcn60l6w7Y8fTiJ+A8t8+ejkVJmd0iEgPUA/7t47KhjBMR6YU96XY3xhwtme7lexCMRFVhnMaYf7t9nIm9L1SybI9Sy/4z4BH+7lT+/e4EHnafEMJj6gtv+xL0Y1qVqn0WAyV3xO8FPiynbJk6QCfBldSr3wp4vDsfABXGKSL1S6pJRKQhcCWw2dg7QdnY+xVelw9hnLHAQmyd5fxS84J9PHOAi8S2forF/icv3XLDfR9uB1Y4x3AxcKfY1kAXABcBXwU4Pp/jFJFLgdeBPxhj8tyme/wehDHOxm4f/wBscd5nAv/txFsf+G9O/lUd8lideFtjb5Z+7jYtlMfUF4uBPzmtfq4ADjoXTsE/pqG66x3sF7YuNwvIBZYDDZzpicBMt3LnY8+q1UotvwLYgE1S/wfUCVecQDcnlm+cvwPdlr8Qm6i2A+8DNcMYZ3/gOPC126tTqI4ntqXENuxV2xhn2jPYJAoQ5xyj7c4xu9Bt2THOcluBG4P83awozuXAXrdjuLii70GY4nwW2OTEkw20dlv2fuc4bwfuC2acvsTqfH4KmFhquVAf0znYVnDHsfX2A4H/Bf7XmS/ANGc/NgCJoTqm+oSvUkpFoapU7aOUUspHmvyVUioKafJXSqkopMlfKaWikCZ/pZSKQpr8lVIqCmnyV0qpKKTJXymlotD/B1GiOalnKtqYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Number of vertices in the graph\n",
    "n_vertices = 20\n",
    "\n",
    "# Generate the graph using the barabasi-albert procedure\n",
    "G = nx.random_graphs.barabasi_albert_graph(n_vertices, 2, seed=3)\n",
    "\n",
    "# Draw the graph\n",
    "nx.draw_networkx(G, pos=nx.spring_layout(G, seed = 12))\n",
    "\n",
    "# Define the function to compute the cut size\n",
    "f_cut = lambda S: nx.cuts.cut_size(G, S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now looking at the graph, if we were to pick the cut defined by $S=\\{6, 14, 17\\}$, we see from the drawing of the graph that the cut value in this case would be 5 (upper right hand side of the figure you will find the nodes 6, 17 and 14). Let's just check this with our function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The cut size for cut {6, 14, 17} is 5\n"
     ]
    }
   ],
   "source": [
    "S = {6, 14, 17}\n",
    "print('The cut size for cut {6, 14, 17} is', f_cut(S))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function is particular because it is submodular. Recall that submodularity of $f$ implies that for all $A, B \\subseteq N$, we should have that $f(A)+f(B) \\geq f(A \\cap B) + f(A \\cup B)$. Let us check this with a simple example, and then with various random subsets $A$ and $B$ of size 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The submodular property is respected since 22 >= 18 as expected.\n",
      "\n",
      "It is True that the submodular property is respected for all the random subsets.\n"
     ]
    }
   ],
   "source": [
    "# Check the submodular property for a simple example\n",
    "A = {0, 2, 13}\n",
    "B = {4, 5, 7, 8, 18}\n",
    "\n",
    "print('The submodular property is respected since', f_cut(A) + f_cut(B), '>=', f_cut(A & B) + f_cut(A | B), 'as expected.')\n",
    "print()\n",
    "\n",
    "# Check submodular property for randomly 50 generated subsets A and B\n",
    "submodular_property_check = []\n",
    "\n",
    "for _ in range(50):\n",
    "    \n",
    "    # Generate random subsets A and B\n",
    "    A = set(random.sample(set(range(20)), k=5))\n",
    "    B = set(random.sample(set(range(20)), k=5))\n",
    "    \n",
    "    # Check if the submodular propery is verified and append result to list\n",
    "    submodular_property_check.append(f_cut(A) + f_cut(B) >= f_cut(A & B) + f_cut(A | B))\n",
    "    \n",
    "print('It is', all(submodular_property_check), 'that the submodular property is respected for all the random subsets.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's use our algorithm and find some maximums! We first start by running the naive maximization with cardinality constraint $k=8$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Processing subsets of size: 100%|██████████| 8/8 [00:07<00:00,  1.03it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The maximum value is 30 and it took 7.791031837463379 seconds to compute it\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# Record starting time (to compute the running time)\n",
    "start_time = time.time()\n",
    "\n",
    "# Run the greedy algorithm\n",
    "\n",
    "max_value_cut, max_set_cut = naive_maximization_under_cc(f_cut, range(len(G)), k=8)\n",
    "\n",
    "# Get the running time of the aglorithm\n",
    "duration = time.time() - start_time\n",
    "\n",
    "print('The maximum value is {} and it took {} seconds to compute it'.format(max_value_cut, duration))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And as comparison, we run the greedy algorithm below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The maximum value is 27 and it took 0.0364840030670166 seconds to compute it\n"
     ]
    }
   ],
   "source": [
    "# Record starting time (to compute the running time)\n",
    "start_time = time.time()\n",
    "\n",
    "# Run the greedy algorithm\n",
    "\n",
    "max_value_cut, max_set_cut = greedy_maximization_under_cc(f_cut, range(len(G)), k=8)\n",
    "\n",
    "# Get the running time of the aglorithm\n",
    "duration = time.time() - start_time\n",
    "\n",
    "print('The maximum value is {} and it took {} seconds to compute it'.format(max_value_cut, duration))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe two things here:\n",
    "   1. The value found by greedy is not the maximal one, so greedy is not always the best.\n",
    "   2. Greedy is much faster than the naive method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coverage\n",
    "\n",
    "In the coverage problem, we are given a universe of elements $U$ and subsets $T_1, T_2, \\dots, T_n \\subseteq U$ of this universe. The submodular function we consider is defined on ground set $N=\\{1, 2, \\dots, n\\}$ (the indices of the subsets $T_i$'s) as follows: $f(S) = |\\cup_{i \\in S} T_i|$ for $S \\subseteq N$, that is, it counts the number of elements of $U$ we cover by taking the union of some specific subset of $T_i$'s.\n",
    "\n",
    "As a first example, we consider $U=\\{1, 2, \\dots, 10000\\}$. We will generate the $T_i$'s as being random subsets of $U$, with their size picked at random too. We first create the Python functions to generate all of this and to create the submodular function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [],
   "source": [
    "U = set(range(10000))\n",
    "\n",
    "def generate_random_subset(U):\n",
    "    \"\"\"Function to generate a random subset of universe U\"\"\"\n",
    "    \n",
    "    # Pick a random size (smaller than the size of U divided\n",
    "    # by 5 so that the T_i's are not too big)\n",
    "    size = random.randrange(1, math.ceil((len(U) + 1) / 5))\n",
    "    \n",
    "    # Get a random subset of U of the specified size\n",
    "    subset = random.sample(range(len(U)), size)\n",
    "    \n",
    "    # Return the random set\n",
    "    return set(subset)\n",
    "\n",
    "def generate_subsets(U, n):\n",
    "    \"\"\"Generate a sequence of random subsets of U\"\"\"\n",
    "    subsets = []\n",
    "    \n",
    "    # Generate as many random subsets as necessary\n",
    "    for _ in range(n):\n",
    "        subsets.append(generate_random_subset(U))\n",
    "    \n",
    "    # Return the list of generated subsets\n",
    "    return subsets\n",
    "\n",
    "def f_coverage(S, all_subsets):\n",
    "    \"\"\"Function that computes the coverage given indices of subsets.\n",
    "    The subsets are given in the argument all_subsets\"\"\"\n",
    "    \n",
    "    # Retrieve only the selected subsets\n",
    "    selected_subsets = [all_subsets[i] for i in S]\n",
    "    \n",
    "    # Compute the number of elements in the union\n",
    "    return len(set.union(*selected_subsets))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we generate 40 random subset of $U$. And we use the function to see how much of the universe all these elements actually cover.\n",
    "\n",
    "We also see how much a few of these cover $U$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Coverage of all subsets: 9805 / 10000\n",
      "Coverage of subsets with indices {0, 4, 8}: 2495 / 10000\n"
     ]
    }
   ],
   "source": [
    "# Generate subsets\n",
    "all_subsets = generate_subsets(U, 40)\n",
    "\n",
    "# Take all the indices (1 to 40 basically) and computes how much all our subsets cover the universe\n",
    "S = set(range(40))\n",
    "print('Coverage of all subsets: {} / {}'.format(f_coverage(S, all_subsets), len(U)))\n",
    "\n",
    "# Take indices {0, 4, 8} and computes how much the corresponding subsets cover the universe\n",
    "S = {0, 4, 8}\n",
    "print('Coverage of subsets with indices {{0, 4, 8}}: {} / {}'.format(f_coverage(S, all_subsets), len(U)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I invite you to try and pick some subsets $A$ and $B$ of indices and then check that the submodularity property holds in this case too.\n",
    "\n",
    "Now, let's use our algorithm and find some maximums! We first start by running the naive maximization with cardinality constraint $k=5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Processing subsets of size: 100%|██████████| 5/5 [03:09<00:00, 37.91s/it]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The maximum value is 6200 and it took 189.56952023506165 seconds to compute it\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# Record starting time (to compute the running time)\n",
    "start_time = time.time()\n",
    "\n",
    "# Run the greedy algorithm\n",
    "max_value_coverage, max_set_coverage = naive_maximization_under_cc(lambda S: f_coverage(S, all_subsets), range(len(all_subsets)), k=5)\n",
    "\n",
    "# Get the running time of the aglorithm\n",
    "duration = time.time() - start_time\n",
    "\n",
    "print('The maximum value is {} and it took {} seconds to compute it'.format(max_value_coverage, duration))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Great, we now know the maximum value when $k=5$, but observe it took basically 3 minutes to compute it!\n",
    "\n",
    "Greedy should definitely do better here, but will not necessarily find the best value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The maximum value is 6200 and it took 0.08582663536071777 seconds to compute it\n"
     ]
    }
   ],
   "source": [
    "# Record starting time (to compute the running time)\n",
    "start_time = time.time()\n",
    "\n",
    "# Run the greedy algorithm\n",
    "greedy_maximization(lambda S: f_coverage(S, all_subsets), range(len(all_subsets)), k=5)\n",
    "\n",
    "# Get the running time of the aglorithm\n",
    "duration = time.time() - start_time\n",
    "\n",
    "print('The maximum value is {} and it took {} seconds to compute it'.format(max_value_coverage, duration))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well, it was able to find the best value in that case, but note the running time, less than 0.1 seconds. It's a huge improvement compared to the naive method. That is basically the reason the field of algorithms exist. One tries to find a method to solve a problem as fast and as reliably as possible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Going further - When greedy does not perform well\n",
    "\n",
    "We have seen from the previous examples how the two algorithms perform. We see that there is the need to develop better algorithms! Indeed, the naive method finds the maximum value, but is extremely slow. On the other hand, greedy is extremely fast, but can sometimes give poor solutions.\n",
    "\n",
    "If you remember, I stated an interesting theorem at the end of the blog post, which says that when maximizing a monotone submodular function with cardinality constraint $k$, we have the following. The set $S$ returned by greedy is such that $f(S) \\geq (1-(1-1/k)^k)f(O) \\geq (1-1/e)f(O)$, where $O$ is the optimal value we wish to find.\n",
    "\n",
    "Well we will see below what we call a \"tight\" example, that is an example where the subset found by greedy is such that $f(S) = (1-(1-1/k)^k)f(O)$, so it is as worst as it can get."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A tight example\n",
    "\n",
    "This example is about coverage functions. Here, the universe of elements is $U = \\{1, 2, \\dots, k\\}^n$, where $k$ is the cardinality constraint, and $n$ is some integer. So $U$ is the universe containing vectors of dimension $n$ whose coordinates can take values in $\\{1, 2, \\dots, k\\}$.\n",
    "\n",
    "In order to get the special case, we \"craft\" some specific subset from which we will need to choose later on. First, we define the sets $S_1, S_2, \\dots, S_n$ as follows:\n",
    "\n",
    "$$S_i = \\{x \\in U: x_i = 1\\}$$ for $i \\in \\{1, 2, \\dots, n\\}$. We also define subsets $T_1, T_2, \\dots, T_k$ as\n",
    "\n",
    "$$T_j = \\{x \\in U: x_1 = j\\}$$ for $j \\in \\{1, 2, \\dots, k\\}$.\n",
    "\n",
    "Note that under cardinality constraint $k$, the optimal choice is to pick all the sets $T_i$'s since their union cover all the universe. However, it can be the case that greedy chooses the $k$ sets only from $S_1, S_2, \\dots, S_n$, and the way they were constructed makes it that they cover only a fraction $1-(1-1/k)^k$ of the universe.\n",
    "\n",
    "Let us implement this for the case $n=7$ and $k=6$. To generate all elements in the universe, the Python function `product` from the package `itertools` is really useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 267,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set the values for n and k\n",
    "n = 7\n",
    "k = 6\n",
    "\n",
    "# Generate all elements of the universe\n",
    "universe_tight = set(product(range(k), repeat=n))\n",
    "\n",
    "# Generate the sets S_1, S_2, ... S_n\n",
    "Ss = [{elem for elem in universe_tight if elem[i] == 1} for i in reversed(range(n))]\n",
    "\n",
    "# Generate the sets T_1, T_2, ..., T_k\n",
    "Ts = [{elem for elem in universe_tight if elem[0] == j} for j in range(k)]\n",
    "\n",
    "# Collect all subset in one big sequence\n",
    "all_subsets = Ss + Ts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned before, the optimal value $f(O)$ of $f$ is attained when we pick subsets $T_1, \\dots, T_k$ (so $O=\\{1, 2, \\dots, k\\}$) since it covers the universe. Obviously in this case, $f(O)$ is equal the size of the universe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 268,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The greedy algorithm gives us f(S) = 186186\n",
      "Also, we have that (1-(1-1/k)**k) * f(O) = 186186.0\n"
     ]
    }
   ],
   "source": [
    "max_value_coverage_tight, _ = greedy_maximization_under_cc(lambda x: f_coverage(x, all_subsets), range(len(all_subsets)), k=k)\n",
    "\n",
    "print('The greedy algorithm gives us f(S) =', max_value_coverage_tight)\n",
    "print('Also, we have that (1-(1-1/k)**k) * f(O) =', (1-(1-1/k)**k) * len(universe_tight))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So in the specific case of monotone submodular functions and maximization under cardinality constraint, this is as bad as greedy can get."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion \n",
    "\n",
    "In this Jupyter Notebook, we have seen basic set functions, and in particular submodular functions. We saw how using the naive and greedy algorithms, we were able to sometimes maximize these functions while having a not too long running time.\n",
    "\n",
    "I invite you to go further and try to define your own set function. How long does the naive algorithm take to solve your maximization problem? What about greedy? You can check also online for more advanced algorithms.\n",
    "\n",
    "For example, when a submodular function is not monotone, one can use an algorithm called the *double greedy*, which gives not too bad results and has the advantage of being very simple to implement. You can have some explanation [here](http://people.csail.mit.edu/stefje/fall15/notes_lecture17.pdf)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "The notation and examples shown are greatly inspired from the lecture notes of the Advanced Algorithms class taught at EPFL in 2019-2020."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:anaconda3]",
   "language": "python",
   "name": "conda-env-anaconda3-py"
  },
  "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}