Step_size_inversion.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Step-size-inversion.ipynb",
      "provenance": [],
      "collapsed_sections": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "S5iIy9s1_0ja"
      },
      "source": [
        "This notebook demonstrates the step-size inversion phenomenon discussed in the paper \"Asymptotic Network Independence and Step-Size for a Distributed Subgradient Method.\" The code provided here will generated Figure 2 of that paper. \n",
        "\n",
        "First, we set the various parameters. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "T2p56WUQ_09a"
      },
      "source": [
        "import numpy as np \n",
        "N = 9; #number of agents\n",
        "m = 1;  #number of data points per agent \n",
        "d = 2;  #dimension "
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "NbrSqQnhAQbF"
      },
      "source": [
        "times = 1000;                                                              #number of times we will run the optimization method to determine each data point\n",
        "stepsizes = [0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85]; #we will use step-sizes that go as 1/k^a, where a is an element of this array\n",
        "L = len(stepsizes)\n",
        "results = np.zeros(L);                                                 #since there are 11 step-sizes, the final result will be an array with 11 entries \n",
        "lambda1 = 1\n",
        "lambda2 = 1/20\n",
        "ones_vector = np.ones(d)\n"
      ],
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-TkdMpvyqvu5"
      },
      "source": [
        "Since the optimization problem is done over the [0,1]^d, it will be useful to have a couple of functions projecting points onto the unit cube. The two functions below do this. The first one projects a numpy array of dimension dim,\n",
        "and the second one projects a scalar."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8pkmYiJGG9yZ"
      },
      "source": [
        "def projcubenp(x,dim):\n",
        "  for i in range(dim):\n",
        "    if x[i] > 1:\n",
        "      x[i]=1\n",
        "\n",
        "    if x[i] < -1:\n",
        "      x[i]=-1\n",
        "  \n",
        "  return x"
      ],
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Sc3FpoPtpeGZ"
      },
      "source": [
        "def projcubescalar(x):\n",
        "  if x > 1:\n",
        "    x=1\n",
        "  if x < -1:\n",
        "    x=-1\n",
        "  return x"
      ],
      "execution_count": 4,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3fDdUyvIrAjc"
      },
      "source": [
        "This function computes the subgradient of the function \n",
        "\n",
        "f(point) = (1/40)*(point^T x - y)^4 + lambda_1||point||_2^2 + lambda_2 ||point||_1 \n",
        "\n",
        "that subgradient is clearly equal to \n",
        "\n",
        "(1/4) (point^T x - y)^3 point + 2 lambda_1 point + lambda_2 sign(point) "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dLqG3ujvRMYi"
      },
      "source": [
        "def compute_grad_elastic(point, x, y, lambda1, lambda2):\n",
        "  # point is theta\n",
        "  #x is a_i\n",
        "  #y is b_i\n",
        "  #lambda1, lambda2 are themselves\n",
        "  sign  = np.sign(point)\n",
        "  diff = np.matmul(point,x) - y\n",
        "  y = (1/10)*(diff ** 3)*x + 2*lambda1*point + lambda2*sign\n",
        "\n",
        "  return y"
      ],
      "execution_count": 5,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6AxmdNbst1UX"
      },
      "source": [
        "Finally, it will be useful to have the following \"double sided relu\" function which is zero when the absolute value of x is below N, and outside that range equals |x|-N"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GhO2SLPbrNJN"
      },
      "source": [
        "def reluN(x,U): \n",
        "  \n",
        "  if np.absolute(x) < U:\n",
        "    return 0\n",
        "  else:\n",
        "    return np.absolute(x) - U"
      ],
      "execution_count": 6,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eGiYr9wHB_lS"
      },
      "source": [
        "Finally, the following package will be used to track progress as we execute the method."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "-v2fkMdqAgZM",
        "outputId": "1ac178bd-68d2-4ed0-8afd-e6784c4e9fec"
      },
      "source": [
        "!pip install enlighten"
      ],
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "Requirement already satisfied: enlighten in /usr/local/lib/python3.7/dist-packages (1.10.1)\n",
            "Requirement already satisfied: blessed>=1.17.7 in /usr/local/lib/python3.7/dist-packages (from enlighten) (1.18.1)\n",
            "Requirement already satisfied: prefixed>=0.3.2 in /usr/local/lib/python3.7/dist-packages (from enlighten) (0.3.2)\n",
            "Requirement already satisfied: wcwidth>=0.1.4 in /usr/local/lib/python3.7/dist-packages (from blessed>=1.17.7->enlighten) (0.2.5)\n",
            "Requirement already satisfied: six>=1.9.0 in /usr/local/lib/python3.7/dist-packages (from blessed>=1.17.7->enlighten) (1.15.0)\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bzAiJKiPuBUS"
      },
      "source": [
        "The following is the main loop of the \"centralized\" method. That is, each iteration of the method will just loop over the data of every agent. This only makes sense under the assumption that this is possible, i.e., a centralized agent can simply access the data of each agent. This method serves as a benchmark to what the distributed method is later compared. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 170
        },
        "id": "4KUWU48_A4mz",
        "outputId": "090ba67d-7d6e-4621-8972-25edd72666e3"
      },
      "source": [
        "import enlighten\n",
        "\n",
        "manager = enlighten.get_manager()\n",
        "ticks = manager.counter(total=L, desc=\"Outer Loop\", unit=\"Iterations\", color=\"red\")\n",
        "\n",
        "\n",
        "for a in range(L):  #there are 11 step-sizes and we loop over them. \n",
        "    ticks.update()\n",
        "    record = np.zeros(times) #this array will contain an entry for every time we run the method \n",
        "    tocks = manager.counter(total=times, desc=\"Inner Loop\", unit=\"Iterations\", color=\"blue\")\n",
        "    for b in range(times): #this loop will fill the \"record\" array\n",
        "        \n",
        "        tocks.update()\n",
        "\n",
        "        X = np.random.rand(m*d,N) #generate random data \n",
        "                                  #this is the a_i in the paper\n",
        "                                  #each column of this array corresponds to information held by an agent\n",
        "                                  #thus there are N columns\n",
        "                                  #each column holds m data points of dimension d\n",
        "\n",
        "        \n",
        "        #the next loop generates the quantities b_i from the paper, which we stack into the vector y\n",
        "        y = np.zeros((m,N)) #there are N agents, each has m data points, so y is m x N\n",
        "        for i in range(N):\n",
        "          for j in range(m):\n",
        "            Xextracted = X[j*d:(j+1)*d,i] #this is x_i\n",
        "            y[j,i] = np.matmul(ones_vector, Xextracted) + 2*np.random.standard_normal() #y_i = 1^T x_i + w_i \n",
        "                                                                                      #this is the all-ones vector is going to be very close to the optimal solution\n",
        "\n",
        "\n",
        "\n",
        "        #now we setup the main loop\n",
        "\n",
        "        seqerror = 1 #this quantity will measure the norm of the subgradient\n",
        "        k=1          #variable counting how many iterations we have done\n",
        "\n",
        "        xprev = np.zeros(d) #initial point\n",
        "        runavg = xprev      #running average of iterates\n",
        "\n",
        "        while seqerror > 3e-2: #as long as the error is large, we keep iterating \n",
        "\n",
        "              #first step is to compute the gradient at the current point, which is xprev\n",
        "              grad = np.zeros(d) #initialize to zero\n",
        "\n",
        "              #sum over all the agents and all the data held by the agents\n",
        "              for i in range(N):\n",
        "                for j in range(m):\n",
        "                  Xextracted = X[j*d:(j+1)*d,i]\n",
        "                  temp_grad = compute_grad_elastic(xprev, Xextracted, y[j,i], lambda1, lambda2)\n",
        "                  grad += temp_grad \n",
        "                \n",
        "              #step-size at this iteration\n",
        "\n",
        "              alpha = 1/( k ** stepsizes[a])\n",
        "\n",
        "              #the main update: projected gradient descent\n",
        "              x = projcubenp(xprev - alpha*(1/(N*m))*grad, d) #note the slightly lopsided way we computed the gradient:\n",
        "                                                              #we sum up over K points, then divide by K\n",
        "                                                              #thus the ||theta||_2^2 and ||theta||_1 terms get added K times\n",
        "                                                              #and then get divided by K\n",
        "\n",
        "              runavg = (1-alpha)*runavg + alpha*x            #update the running average \n",
        "\n",
        "\n",
        "              #next, we will compute the gradient at the running average \n",
        "              #we will stop when this is close to zero \n",
        "              #the rule we use will be explained in the text box after the code \n",
        "              #for now note that this code below computes \n",
        "              #gradrun0: subgradient at the running average with lambda_2 = 0 \n",
        "              #gradrun1: subgradient at the running average \n",
        "              gradrun0 = np.zeros(d)\n",
        "              gradrun1 = np.zeros(d)\n",
        "\n",
        "              for i in range(N):\n",
        "                for j in range(m):\n",
        "                  Xextracted = X[j*d:(j+1)*d,i]\n",
        "                  temp_grad0 = compute_grad_elastic(runavg, Xextracted, y[j,i], lambda1, 0)\n",
        "                  temp_grad1 = compute_grad_elastic(runavg, Xextracted, y[j,i], lambda1, lambda2)\n",
        "                  gradrun0 += temp_grad0\n",
        "                  gradrun1 += temp_grad1\n",
        "\n",
        "              #because the above loop summed over N*m data points, we need to now divide by N*m\n",
        "              gradrun0 = gradrun0/(N*m)\n",
        "              gradrun1 = gradrun1/(N*m)\n",
        "              #compute the norm of the subgradient at the running average \n",
        "              error_cum = 0 \n",
        "\n",
        "              for i in range(d):\n",
        "\n",
        "                if np.absolute(runavg[i]) > 1e-2 and np.absolute(runavg[i]) < 1-1e-2:\n",
        "                  error_cum += np.absolute(gradrun1[i])\n",
        "\n",
        "                if np.absolute(runavg[i]) < 1e-2: \n",
        "                  error_cum += reluN(gradrun0[i],lambda2)\n",
        "\n",
        "                if np.absolute(runavg[i]) > 1-1e-2:\n",
        "                  error_cum += (runavg[i] - projcubescalar(runavg[i]-alpha*gradrun1[i]))/alpha\n",
        "\n",
        "\n",
        "              seqerror = error_cum #set the cumulative error which will determine whether the while loop takes another iteration \n",
        "              xprev = x #update the iterate \n",
        "              k=k+1     #update the number of iterations\n",
        "        record[b]=k     #once the while loop has completed, count how many iterations it took\n",
        "    tocks.close()       #close the inner loop bar once the for loop over b is complete (it will get started again the next time this loop is passed through)\n",
        "\n",
        "    results[a] = np.average(record)\n",
        "\n",
        "\n",
        "manager.stop()\n"
      ],
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "<style>\n",
              ".enlighten-fg-red {\n",
              "  color: #cd0000;\n",
              "}\n",
              ".enlighten-fg-blue {\n",
              "  color: #0000ee;\n",
              "}\n",
              "</style>\n",
              "<div class=\"enlighten\">\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Outer Loop 100%|<span class=\"enlighten-fg-red\">██████████████████████████████████████████████</span>| 8/8 [00:24&lt;00:00, 0.33 Iterations/s]</pre>\n",
              "  </div>\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Inner Loop 100%|<span class=\"enlighten-fg-blue\">███████████████████████████████████████</span>| 1000/1000 [00:24&lt;00:00, 41.81 Iterations/s]</pre>\n",
              "  </div>\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Inner Loop 100%|<span class=\"enlighten-fg-blue\">███████████████████████████████████████</span>| 1000/1000 [00:21&lt;00:00, 47.03 Iterations/s]</pre>\n",
              "  </div>\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Inner Loop 100%|<span class=\"enlighten-fg-blue\">███████████████████████████████████████</span>| 1000/1000 [00:18&lt;00:00, 54.23 Iterations/s]</pre>\n",
              "  </div>\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Inner Loop 100%|<span class=\"enlighten-fg-blue\">███████████████████████████████████████</span>| 1000/1000 [00:16&lt;00:00, 64.29 Iterations/s]</pre>\n",
              "  </div>\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Inner Loop 100%|<span class=\"enlighten-fg-blue\">███████████████████████████████████████</span>| 1000/1000 [00:13&lt;00:00, 79.84 Iterations/s]</pre>\n",
              "  </div>\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Inner Loop 100%|<span class=\"enlighten-fg-blue\">██████████████████████████████████████</span>| 1000/1000 [00:09&lt;00:00, 110.25 Iterations/s]</pre>\n",
              "  </div>\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Inner Loop 100%|<span class=\"enlighten-fg-blue\">██████████████████████████████████████</span>| 1000/1000 [00:05&lt;00:00, 206.55 Iterations/s]</pre>\n",
              "  </div>\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Inner Loop 100%|<span class=\"enlighten-fg-blue\">██████████████████████████████████████</span>| 1000/1000 [00:06&lt;00:00, 171.50 Iterations/s]</pre>\n",
              "  </div>\n",
              "</div>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5FnTpyxRaPuk"
      },
      "source": [
        "We now explain the way we computed the subgradinet. First, the motivation: as explained in the paper, we want to choose a subgradient of minimum norm. Otherwise, we might be at the optimal point and move away from it with a nonzero subgradient. This is particularly important since we use the 1-norm of the subgradient as a stopping criterion. \n",
        "\n",
        "Clearly, the subgradient of our objective decomposes component by component. We thus consider three cases:\n",
        "\n",
        "(1) theta_i is not 0 or 1. \n",
        "\n",
        "In this case, the situation is simple. There is only one choice for subgradient. This is exactly what is computed by gradrun1. \n",
        "\n",
        "\n",
        "(2) theta_i is zero. In this case, there is considerable freedom to choose the subgradient of the last term: it can be chosen to be anything between -lambda2 and lambda2 \n",
        "\n",
        "So we choose it to make that entry of the subgradient the smallest. What this means is that:\n",
        "\n",
        "  (a) we compute the subgradient with lambda_2=0, for which there is only one choice; this is gradrun0\n",
        "\n",
        "  (b) the final answre is zero with |gradrun0[i]|<lambda2. Otherwise, it is |gradrun0[i]|-lambda_2. This is computed by the reluN function. \n",
        "\n",
        "\n",
        "(3) theta_i=1. In this case, you compute the norm of the \"gradient mapping\" "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 312
        },
        "id": "5emFTGavypos",
        "outputId": "846c9a5b-c56f-4cc5-cc8f-3c7b44ceb7e9"
      },
      "source": [
        "import matplotlib.pyplot as plt\n",
        "plt.scatter(stepsizes,results)\n",
        "plt.xlabel('beta')\n",
        "plt.title('Iterations vs beta')"
      ],
      "execution_count": 15,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Text(0.5, 1.0, 'Iterations vs beta')"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 15
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEWCAYAAAB/tMx4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAVNUlEQVR4nO3dfZBdd33f8fcHSeA1GNauBFhrZPFkGQxp1CjlYUoCA0GGaWsVaMGEFgeIW6ahnSYRoNAG2g41RUzJA8m4pjgmCZh6qKKSpqnsQoDC2CQyMpEdIzDmwV5hLGMvhWRrFPHtH/esuVr2ee/u3p/0fs3cuef+ztP3/rTz0bm/c+65qSokSe152FoXIElaGgNckhplgEtSowxwSWqUAS5JjTLAJalRBrialeR7SZ601nUsRpKtSSrJ+rWuRe0zwLUkSb6W5EXd9GVJPrPC+/tkkjf0t1XVo6rqzpXc7zAx/DWdAa41ZyBJS2OAa1mSPA24EnhON6Qx0bU/Isl7knwjybeSXJlkpJv3/CR3J3lLknuA30lydpL/keRYkge66fO65d8JPA94X7eP93XtleQp3fRjkvxut/7Xk/zrJA/r5l2W5DNdPQ8k+WqSl/S9h8uS3Jnku928n53hfW5OMpnknL627UnuS7IhyVOSfCrJd7q2/zpP170uydEk30zyy33bfFiStyb5SpJvJ7mub5+f7p4nun54TpInJ/lEt+x9ST6UZHTh/4JqmQGuZamq24F/BtzYDWlMhce7gAuAHweeAowBv9q36uOBc4Dzgcvp/S3+Tvd6CzAJvK/bx9uA/wP8QrePX5ihlN8EHgM8Cfhp4J8AP9c3/1nAEWAj8G7gA+l5JPAbwEuq6izgucAtM7zPo8CNwMv7ml8NfLSqjgP/HrgeOBs4r6tnLi8Angq8GHjL1HAU8CZgV/ceNgMPAL/Vzfup7nm064cbgQBXdMs+DXgC8I559q1TRVX58LHoB/A14EXd9GXAZ/rmBfhL4Ml9bc8BvtpNPx/4PnDGHNv/ceCBvtefBN4wbZmi95/Dum57T++b90+BT/bVd0ffvDO7dR8PPBKYoBfMI/O85zcAn+h7j3cBP9W9/l3gKuC8ebaxtdv3hX1t7wY+0E3fDrywb965wHFgfd+66+fY/i7g0Fr/ffhYnYdH4FoJm+iF5M1JJrphlf/VtU85VlX/b+pFkjOT/Odu+OP/0hsuGE2ybgH72whsAL7e1/Z1ekf9U+6Zmqiqv+omH1VVfwm8kt6niG8m+aMkF86yn/9Gb6joXHpHwz+g98kA4M30Qv1Pk9yW5HXz1HzXtFo3d9PnA3/Q12+3AyeAx820kSSPS/KRJONdv/0+vf7QacAA1yBMv6XlffSGQC6qqtHu8ZiqetQc6/wSsA14VlU9mh8OF2SW5afv7zi98JuyBRhfUPFVB6rqZ+gd7X4ReP8syz1Ab5jklfSGTz5S1Tvsrap7qurnq2ozvaP/354an5/FE6bVerSbvovecM5o3+OMqhpn5j74D137M7t+ew0/7DOd4gxwDcK3gPOSPBygqn5ALwTfm+SxAEnGkuycYxtn0Qv9ie6k3dtn2MeM13xX1QngOuCdSc5Kcj7wi/SORufUHcFe0o2FPwh8j96R9Ww+TG98/RXd9NR2/uHUSVd649Y1z3b+Tfep4yJ6Y/VTJz2v7N7H+d12NyW5pJt3rNtmfz+c1dX8nSRjwO753rNOHQa4BuETwG3APUnu69reAtwB3NR9tP/f9I6wZ/NrwAi9o+mb6A259Pt14BXdVSS/McP6b6I37n4n8Bl64Xr1Amp/GL2wPwrcT+/k4RvnWP5j9E4+3lNVX+hr/0ngc0m+1y3zL2vua9Q/Ra9/Pg68p6qu79p/vVv/+iTfpdcXz4KHhn7eCXy2G2J5NvBvgb8FfAf4I2DfAt6zThHpPgFKkhrjEbgkNcoAl6RGGeCS1CgDXJIatao3Edq4cWNt3bp1NXcpSc27+eab76uqTdPbVzXAt27dysGDB1dzl5LUvCRfn6ndIRRJapQBLkmNMsAlqVEGuCQ1ygCXpEb5W4SStIL2Hxpn74EjHJ2YZPPoCLt3bmPX9rH5V1wAA1ySVsj+Q+Ps2XeYyeMnABifmGTPvsMAAwlxh1AkaYXsPXDkofCeMnn8BHsPHBnI9g1wSVohRycmF9W+WAa4JK2QzaMji2pfLANcklbI7p3bGNlw8u9yj2xYx+6dc/041cJ5ElOSVsjUiUqvQpGkBu3aPjawwJ7OIRRJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSo+YN8CRXJ7k3ya3T2t+U5ItJbkvy7pUrUZI0k4UcgV8DXNzfkOQFwCXA36yqi4D3DL40SdJc5g3wqvo0cP+05jcC76qqB7tl7l2B2iRJc1jqGPgFwPOSfC7Jp5L85GwLJrk8ycEkB48dO7bE3UmSpltqgK8HzgGeDewGrkuSmRasqquqakdV7di0adMSdydJmm6pAX43sK96/hT4AbBxcGVJkuaz1ADfD7wAIMkFwMOB+wZVlCRpfvPeDzzJtcDzgY1J7gbeDlwNXN1dWvh94LVVVStZqCTpZPMGeFVdOsus1wy4FknSIvhNTElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVHzBniSq5Pcm+TWGeb9UpJKsnFlypMkzWYhR+DXABdPb0zyBODFwDcGXJMkaQHmDfCq+jRw/wyz3gu8GahBFyVJmt+SxsCTXAKMV9UXBlyPJGmB1i92hSRnAr9Cb/hkIctfDlwOsGXLlsXuTpI0i6UcgT8ZeCLwhSRfA84DPp/k8TMtXFVXVdWOqtqxadOmpVcqSTrJoo/Aq+ow8Nip112I76iq+wZYlyRpHvMGeJJrgecDG5PcDby9qj6w0oVJ0kz2Hxpn74EjHJ2YZPPoCLt3bmPX9rG1LmtNzBvgVXXpPPO3DqwaSZrD/kPj7Nl3mMnjJwAYn5hkz77DAKdliPtNTEnN2HvgyEPhPWXy+An2HjiyRhWtLQNcUjOOTkwuqv1UZ4BLasbm0ZFFtZ/qDHBJzdi9cxsjG9ad1DayYR27d25bo4rW1qIvI5SktTJ1otKrUHoMcElN2bV97LQN7OkcQpGkRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNmjfAk1yd5N4kt/a17U3yxSR/nuQPkoyubJmSpOkWcgR+DXDxtLYbgGdU1Y8BXwL2DLguSdI85g3wqvo0cP+0tuur6q+7lzcB561AbZKkOQxiDPx1wB8PYDuSpEVYVoAneRvw18CH5ljm8iQHkxw8duzYcnYnSeqz5ABPchnwd4GfraqabbmquqqqdlTVjk2bNi11d5KkadYvZaUkFwNvBn66qv5qsCVJkhZiIZcRXgvcCGxLcneS1wPvA84CbkhyS5IrV7hOSdI08x6BV9WlMzR/YAVqkSQtgt/ElKRGLWkMXNKpZf+hcfYeOMLRiUk2j46we+c2dm0fW+uyNA8DXDrN7T80zp59h5k8fgKA8YlJ9uw7DGCIDzmHUKTT3N4DRx4K7ymTx0+w98CRNapIC2WAS6e5oxOTi2rX8DDApdPc5tGRRbVreBjg0mlu985tjGxYd1LbyIZ17N65bY0q0kJ5ElM6zU2dqPQqlPYY4JLYtX3MwG6QQyiS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIatX6tC5BORfsPjbP3wBGOTkyyeXSE3Tu3sWv72FqXpVPMvEfgSa5Ocm+SW/vazklyQ5Ivd89nr2yZUjv2Hxpnz77DjE9MUsD4xCR79h1m/6HxtS5Np5iFDKFcA1w8re2twMer6qnAx7vXkoC9B44wefzESW2Tx0+w98CRNapIp6p5A7yqPg3cP635EuCD3fQHgV0Drktq1tGJyUW1S0u11JOYj6uqb3bT9wCPG1A9UvM2j44sql1aqmVfhVJVBdRs85NcnuRgkoPHjh1b7u6kobd75zZGNqw7qW1kwzp279y2RhXpVLXUAP9WknMBuud7Z1uwqq6qqh1VtWPTpk1L3J3Ujl3bx7jiZc9kbHSEAGOjI1zxsmd6FYoGbqmXEX4MeC3wru75vw+sIukUsGv7mIGtFbeQywivBW4EtiW5O8nr6QX3zyT5MvCi7rUkaRXNewReVZfOMuuFA65FkrQIfpVekhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVHz/iq9NCz2Hxpn74EjHJ2YZPPoCLt3bmPX9rG1LktaMwa4mrD/0Dh79h1m8vgJAMYnJtmz7zCAIa7TlkMoasLeA0ceCu8pk8dPsPfAkTWqSFp7BriacHRiclHt0unAAFcTNo+OLKpdOh0Y4GrC7p3bGNmw7qS2kQ3r2L1z2xpVJK09T2Kexlq6qmOqrlbqlVaDAX6aavGqjl3bx4a2NmktOIRymvKqDql9Bvhpyqs6pPY5hDJgrYwrbx4dYXyGsPaqDqkdHoEP0NS48vjEJMUPx5X3Hxpf69J+hFd1SO0zwAeopXHlXdvHuOJlz2RsdIQAY6MjXPGyZw7lpwVJM3MIZYBaG1f2qg6pbR6BD5DfFpS0mpYV4En+VZLbktya5NokZwyqsBY5rixpNS05wJOMAf8C2FFVzwDWAa8aVGEtclxZ0mpa7hj4emAkyXHgTODo8ktqm+PKklbLkgO8qsaTvAf4BjAJXF9V109fLsnlwOUAW7ZsWfR+WrmuWpJW23KGUM4GLgGeCGwGHpnkNdOXq6qrqmpHVe3YtGnTovbR0nXVkrTalnMS80XAV6vqWFUdB/YBzx1MWT0tXVctSattOQH+DeDZSc5MEuCFwO2DKaunteuqJWk1LTnAq+pzwEeBzwOHu21dNaC6AK+rlqS5LOs68Kp6e1VdWFXPqKp/XFUPDqow8LpqSZrLUH+V3l9hkaTZDXWAg9dVS9JsvBeKJDXKAJekRhngktQoA1ySGmWAS1KjUlWrt7PkGPD1Ja6+EbhvgOWstJbqbalWaKvelmqFtuptqVZYXr3nV9WP3ExqVQN8OZIcrKoda13HQrVUb0u1Qlv1tlQrtFVvS7XCytTrEIokNcoAl6RGtRTgA71R1ipoqd6WaoW26m2pVmir3pZqhRWot5kxcEnSyVo6Apck9THAJalRQxHgSS5OciTJHUneOsP8y5IcS3JL93hD37zXJvly93jtkNd6oq/9Yytd60Lq7Zb5R0n+IsltST7c1z5UfTtPrUPXt0ne21fTl5JM9M0bqr6dp9Zh7NstSf4kyaEkf57kpX3z9nTrHUmyc1hrTbI1yWRf31656J1X1Zo+gHXAV4AnAQ8HvgA8fdoylwHvm2Hdc4A7u+ezu+mzh7HWbt73hrBvnwocmuo34LFD3Lcz1jqsfTtt+TcBVw9r385W67D2Lb0Tgm/spp8OfK1v+gvAI+j94PpXgHVDWutW4Nbl7H8YjsD/NnBHVd1ZVd8HPkLv1+4XYidwQ1XdX1UPADcAF69QnbC8WtfCQur9eeC3uv6jqu7t2oexb2erdS0s9m/hUuDabnoY+3a2WtfCQuot4NHd9GOAo930JcBHqurBqvoqcEe3vWGsddmGIcDHgLv6Xt/dtU338u7jx0eTPGGR6w7KcmoFOCPJwSQ3Jdm1gnVOWUi9FwAXJPlsV9fFi1h3kJZTKwxn3wKQ5Hx6R4OfWOy6A7KcWmE4+/YdwGuS3A38T3qfGha67iAtp1aAJ3ZDK59K8rzF7nwYAnwh/hDYWlU/Ru9o5YNrXM9c5qr1/Op9lfbVwK8lefJaFDjNenpDE8+nd+T1/iSja1rR7OaqdRj7dsqrgI9W1Ym1LmQBZqp1GPv2UuCaqjoPeCnwe0mGNc9mq/WbwJaq2g78IvDhJI+eYzs/Yhje8DjQf5R6Xtf2kKr6dv3wB5P/C/ATC113wJZTK1U13j3fCXwS2L6CtcLC+udu4GNVdbz7yPkleiE5dH07R63D2rdTXsXJQxLD2LdTptc6rH37euC6rq4bgTPo3SxqGPt2xlq7YZ5vd+030xtLv2BRe1+pwf1FnARYT+8kzhP54UmAi6Ytc27f9D8AbuqmzwG+Su9E0Nnd9DlDWuvZwCO66Y3Al5njRNIq1nsx8MG+uu4C/saQ9u1stQ5l33bLXQh8je5Lc8P6dztHrUPZt8AfA5d100+jN64c4CJOPol5Jyt7EnM5tW6aqo3eSdDxxf4drNg/wiI74aX0jqa+Aryta/t3wN/vpq8Abus650+AC/vWfR29ExV3AD83rLUCzwUOd+2HgdcPSd8G+E/AX3R1vWqI+3bGWoe1b7vX7wDeNcO6Q9W3s9U6rH1L72qOz3Z13QK8uG/dt3XrHQFeMqy1Ai/vsuIW4PPA31vsvv0qvSQ1ahjGwCVJS2CAS1KjDHBJapQBLkmNMsAlqVEGuE5p3R3fbl3E8pcl2bySNUmDYoBLJ7sMMMDVBANcp4P1ST6U5PbuBmNnJvmJ7gZCNyc5kOTcJK8AdgAf6u7PPJLkV5P8WZJbk1yVJGv9ZqQpfpFHp7QkW+l9Vf3vVNVnk1wN3E7vNgeXVNWxJK8EdlbV65J8EvjlqjrYrX9OVd3fTf8ecF1V/eEavBXpR6xf6wKkVXBXVX22m/594FeAZwA3dAfU6+jdGW4mL0jyZuBMevcwuY3eHSelNWeA63Qw/WPmd4Hbquo5c62U5Azgt4EdVXVXknfQu5OcNBQcA9fpYEuSqbB+NXATsGmqLcmGJBd1878LnNVNT4X1fUkeBbxitQqWFsIA1+ngCPDPk9xO7/aov0kvjP9jkqk7xD23W/Ya4MoktwAPAu8HbgUOAH+2ynVLc/IkpiQ1yiNwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIa9f8BNSCBMYP6ijgAAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8BCFSyBIg5sq"
      },
      "source": [
        "We now code the distributed method. \n",
        "\n",
        "\n",
        "The first step is given in the following chunk of code, which returns P to be the Metropolis matrix on the path. \n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "m3-WaNBXoVFE"
      },
      "source": [
        "import networkx as nx\n",
        "\n",
        "def adj_to_met(U, dimen):\n",
        "  #converts adjacency matrix to metropolis matrix\n",
        "\n",
        "  degs = U.sum(axis=0) #vector of degrees\n",
        "  \n",
        "\n",
        "  X = np.zeros((dimen, dimen)) #initialize a zeros dimen x dimen matrix\n",
        "  \n",
        "  #loop below sets X[i,j] to be U[i,j] divided by the maximum degrees of i and j, plus one \n",
        "  #thus zero entries in U translate to zero entries in X\n",
        "  #if necessary, can try to rewrite all this in the future so that it does not take quadratic time\n",
        "  for i in range(dimen):\n",
        "    for j in range(dimen):\n",
        "      X[i,j] = float(U[i,j])/float((max(degs[i]+1, degs[j]+1)))\n",
        "  \n",
        "  \n",
        "  #finally, we need to take care of the diagonal entries\n",
        "\n",
        "  s = X.sum(axis=0) #sum of the rows of X\n",
        "  d = np.ones(dimen) - s #what remains to make things add up to one \n",
        "\n",
        "  X += np.diag(d) #put that remainder on the diagonal and add it to X\n",
        "\n",
        "  return X\n",
        "\n",
        "G = nx.grid_graph(dim=(3,3)) #create an nxn path\n",
        "Q = nx.adjacency_matrix(G)   #adjacency matrix of this path\n",
        "U = Q.A                      #convert to a numpy array\n",
        "P = adj_to_met(U, N) \n"
      ],
      "execution_count": 11,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eGR_9sNC2srw"
      },
      "source": [
        "Now we do the main iteration loop. This is similar to the centralized case above, except that now the \"iterate\" is a matrix (one column for each agent)\n",
        "and is updated with a consensus term. \n",
        "\n",
        "We stop when the gradient at the running average across network + time is sufficiently close. "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gJxUt5a7g5DW",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 170
        },
        "outputId": "bd6ed8d6-ae54-466c-ac1e-16b10fcaea7d"
      },
      "source": [
        "import enlighten\n",
        "\n",
        "manager = enlighten.get_manager()\n",
        "ticks = manager.counter(total=L, desc=\"Outer Loop\", unit=\"Iterations\", color=\"red\")\n",
        "\n",
        "results_dist = np.zeros(L)\n",
        "\n",
        "\n",
        "for a in range(L):  #there are 11 step-sizes and we loop over them. \n",
        "    ticks.update()\n",
        "    record = np.zeros(times) #this array will contain an entry for every time we run the method \n",
        "    tocks = manager.counter(total=times, desc=\"Inner Loop\", unit=\"Iterations\", color=\"blue\")\n",
        "    for b in range(times): #this loop will fill the \"record\" array\n",
        "        \n",
        "        tocks.update()\n",
        "\n",
        "        X = np.random.rand(m*d,N) #generate random data \n",
        "                                  #this is the a_i in the paper\n",
        "                                  #each column of this array corresponds to information held by an agent\n",
        "                                  #thus there are N columns\n",
        "                                  #each column holds m data points of dimension d\n",
        "\n",
        "        \n",
        "        #the next loop generates the quantities b_i from the paper, which we stack into the vector y\n",
        "        y = np.zeros((m,N)) #there are N agents, each has m data points, so y is m x N\n",
        "        for i in range(N):\n",
        "          for j in range(m):\n",
        "            Xextracted = X[j*d:(j+1)*d,i] #this is x_i\n",
        "            y[j,i] = np.matmul(ones_vector, Xextracted) + 2*np.random.standard_normal() #y_i = 1^T x_i + w_i \n",
        "                                                                                      #this is the all-ones vector is going to be very close to the optimal solution\n",
        "\n",
        "\n",
        "\n",
        "        #now we setup the main loop\n",
        "\n",
        "        seqerror = 1 #this quantity will measure the norm of the subgradient\n",
        "        k=1          #variable counting how many iterations we have done\n",
        "\n",
        "        xprev = np.zeros((d, N)) #initial point. \n",
        "                                 #note the difference from the centralized case\n",
        "                                 #now the initial point is a dxN matrix\n",
        "                                 #the i'th column of this matrix is the initial point at node i\n",
        "\n",
        "        runavg = np.zeros(d)      #running average of iterates\n",
        "                                  #this is a vector in R^d\n",
        "\n",
        "        while seqerror > 3e-2: #as long as the error is large, we keep iterating \n",
        "\n",
        "              #first step is to compute the gradient at the current point, which is xprev\n",
        "              grad = np.zeros((d,N)) #initialize to zero\n",
        "                                     #we will have different gradients at different nodes\n",
        "                                     #the gradient at node i will be the ith column of this\n",
        "\n",
        "              #compute the gradients\n",
        "              for i in range(N):\n",
        "                for j in range(m):\n",
        "                  Xextracted = X[j*d:(j+1)*d,i]\n",
        "                  temp_grad = compute_grad_elastic(xprev[:,i], Xextracted, y[j,i], lambda1, lambda2)\n",
        "                  grad[:,i] += temp_grad \n",
        "                \n",
        "              #step-size at this iteration\n",
        "\n",
        "              alpha = 1/( k ** stepsizes[a])\n",
        "              x = np.zeros((d,N))\n",
        "              #the main update: first, projected gradient descent\n",
        "              for i in range(N):\n",
        "                x[:,i] = projcubenp(xprev[:,i] - alpha*(1/(N*m))*grad[:,i], d) #note the slightly lopsided way we computed the gradient:\n",
        "                                                              #we sum up over K points, then divide by K\n",
        "                                                              #thus the ||theta||_2^2 and ||theta||_1 terms get added K times\n",
        "                                                              #and then get divided by K\n",
        "              \n",
        "              \n",
        "              #second part of main update: the consensus multiplication\n",
        "              # maybe more helpful to think of this as x <--- (P x^T)^T\n",
        "              #which works out to the same thing since P is symmetric\n",
        "              x = np.matmul(x,P)\n",
        "\n",
        "              xavg = x.sum(axis=1)/N                            #average across the network\n",
        "              runavg = (1-alpha)*runavg + alpha*xavg            #update the running average \n",
        "\n",
        "\n",
        "              #the next part is identical to the centralized case\n",
        "              gradrun0 = np.zeros(d)\n",
        "              gradrun1 = np.zeros(d)\n",
        "\n",
        "              for i in range(N):\n",
        "                for j in range(m):\n",
        "                  Xextracted = X[j*d:(j+1)*d,i]\n",
        "                  temp_grad0 = compute_grad_elastic(runavg, Xextracted, y[j,i], lambda1, 0)\n",
        "                  temp_grad1 = compute_grad_elastic(runavg, Xextracted, y[j,i], lambda1, lambda2)\n",
        "                  gradrun0 += temp_grad0\n",
        "                  gradrun1 += temp_grad1\n",
        "\n",
        "              #because the above loop summed over N*m data points, we need to now divide by N*m\n",
        "              gradrun0 = gradrun0/(N*m)\n",
        "              gradrun1 = gradrun1/(N*m)\n",
        "              #compute the norm of the subgradient at the running average \n",
        "              error_cum = 0 \n",
        "\n",
        "              for i in range(d):\n",
        "\n",
        "                if np.absolute(runavg[i]) > 1e-2 and np.absolute(runavg[i]) < 1-1e-2:\n",
        "                  error_cum += np.absolute(gradrun1[i])\n",
        "\n",
        "                if np.absolute(runavg[i]) < 1e-2: \n",
        "                  error_cum += reluN(gradrun0[i],lambda2)\n",
        "\n",
        "                if np.absolute(runavg[i]) > 1-1e-2:\n",
        "                  error_cum += (runavg[i] - projcubescalar(runavg[i]-alpha*gradrun1[i]))/alpha\n",
        "\n",
        "\n",
        "              seqerror = error_cum #set the cumulative error which will determine whether the while loop takes another iteration \n",
        "              xprev = x #update the iterate \n",
        "              k=k+1     #update the number of iterations\n",
        "        record[b]=k     #once the while loop has completed, count how many iterations it took\n",
        "    tocks.close()       #close the inner loop bar once the for loop over b is complete (it will get started again the next time this loop is passed through)\n",
        "\n",
        "    results_dist[a] = np.average(record)\n",
        "\n",
        "\n",
        "manager.stop()\n"
      ],
      "execution_count": 12,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/html": [
              "<style>\n",
              ".enlighten-fg-red {\n",
              "  color: #cd0000;\n",
              "}\n",
              ".enlighten-fg-blue {\n",
              "  color: #0000ee;\n",
              "}\n",
              "</style>\n",
              "<div class=\"enlighten\">\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Outer Loop 100%|<span class=\"enlighten-fg-red\">██████████████████████████████████████████████</span>| 8/8 [35:44&lt;00:00, 0.00 Iterations/s]</pre>\n",
              "  </div>\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Inner Loop 100%|<span class=\"enlighten-fg-blue\">████████████████████████████████████████</span>| 1000/1000 [35:44&lt;00:00, 0.47 Iterations/s]</pre>\n",
              "  </div>\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Inner Loop 100%|<span class=\"enlighten-fg-blue\">████████████████████████████████████████</span>| 1000/1000 [21:45&lt;00:00, 0.77 Iterations/s]</pre>\n",
              "  </div>\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Inner Loop 100%|<span class=\"enlighten-fg-blue\">████████████████████████████████████████</span>| 1000/1000 [15:41&lt;00:00, 1.06 Iterations/s]</pre>\n",
              "  </div>\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Inner Loop 100%|<span class=\"enlighten-fg-blue\">████████████████████████████████████████</span>| 1000/1000 [12:06&lt;00:00, 1.38 Iterations/s]</pre>\n",
              "  </div>\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Inner Loop 100%|<span class=\"enlighten-fg-blue\">████████████████████████████████████████</span>| 1000/1000 [09:35&lt;00:00, 1.74 Iterations/s]</pre>\n",
              "  </div>\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Inner Loop 100%|<span class=\"enlighten-fg-blue\">████████████████████████████████████████</span>| 1000/1000 [07:18&lt;00:00, 2.28 Iterations/s]</pre>\n",
              "  </div>\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Inner Loop 100%|<span class=\"enlighten-fg-blue\">████████████████████████████████████████</span>| 1000/1000 [04:30&lt;00:00, 3.71 Iterations/s]</pre>\n",
              "  </div>\n",
              "  <div class=\"enlighten-bar\">\n",
              "    <pre>Inner Loop 100%|<span class=\"enlighten-fg-blue\">████████████████████████████████████████</span>| 1000/1000 [08:43&lt;00:00, 1.91 Iterations/s]</pre>\n",
              "  </div>\n",
              "</div>\n"
            ],
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ]
          },
          "metadata": {
            "tags": []
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 312
        },
        "id": "hEPrAEEi2Dsp",
        "outputId": "2e89ce3f-d8cc-4df3-dcb2-7a8a75033ba7"
      },
      "source": [
        "plt.scatter(stepsizes,results_dist)\n",
        "plt.xlabel('beta')\n",
        "plt.title('Iterations vs beta')"
      ],
      "execution_count": 14,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Text(0.5, 1.0, 'Iterations vs beta')"
            ]
          },
          "metadata": {
            "tags": []
          },
          "execution_count": 14
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEWCAYAAACKSkfIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAbwklEQVR4nO3df5RcZZ3n8feHJEIrYMLQk0k6CQEMURhnE6gF1IHFo0MCZ8dEZTG4IwSQ6Ay44+pEie4OrB5W1ugwMihs0Awwww8zGEJmxAkRBBYOUTokEgJGmgCbdEJoCA0ovUyI3/3jPgU3TfXPqu6u6vt5nVOnbz3317ducj5163lu3VJEYGZmxbDfSBdgZmbDx6FvZlYgDn0zswJx6JuZFYhD38ysQBz6ZmYF4tC3QpH0G0lHjHQdAyFpuqSQNHaka7HG59C3YSPpaUkfTtMLJd0/xPu7R9Kn820RcWBEbB3K/dYTv2FYdw59a0gOMbPBcejbsJP0HuAa4H2pu6Uzte8v6VuS/q+kXZKukdSU5p0iabukL0t6Fvh7SRMk/YukDkkvpukpafnLgJOAq9I+rkrtIeldafqdkm5I6z8j6b9J2i/NWyjp/lTPi5KeknRa7jUslLRV0itp3n+u8DonS+qSdEiubbak5yWNk/QuSfdKeim1/bCPQ3eepB2Sdkr6q9w295N0saQnJb0gaUVun/elv53pOLxP0pGS7k7LPi/pRknj+/8vaI3MoW/DLiIeBz4LPJi6W8qBczlwFDALeBfQAvx1btU/AA4BDgMWkf3//fv0fBrQBVyV9vFV4P8AF6V9XFShlL8D3gkcAfwH4Gzg3Nz8E4AtwKHAN4EfKPMO4ErgtIg4CHg/sLHC69wBPAh8PNf8SeDWiNgDfB24E5gATEn19OaDwAzgVODL5a4y4HPA/PQaJgMvAt9N805Of8en4/AgIOAbadn3AFOBS/vYt40WEeGHH8PyAJ4GPpymFwL35+YJ+C1wZK7tfcBTafoU4N+AA3rZ/izgxdzze4BPd1smyN5QxqTtHZ2b9xngnlx9bbl5b0/r/gHwDqCTLMyb+njNnwbuzr3GbcDJ6fkNwDJgSh/bmJ72/e5c2zeBH6Tpx4EP5eZNAvYAY3Prju1l+/OBDSP9/8OP4Xn4TN/qRTNZsK6X1Jm6fP41tZd1RMT/Kz+R9HZJ/zt1zbxM1pUxXtKYfuzvUGAc8Eyu7RmyTxdlz5YnIuLVNHlgRPwW+ATZp5Wdkn4s6d097OdHZN1Yk8jOun9H9gkE4EtkbwS/kLRZ0nl91LytW62T0/RhwG254/Y4sBeYWGkjkiZKukVSezpu/0h2PKwAHPo2Urrf3vV5su6ZYyJifHq8MyIO7GWdLwIzgRMi4mDe7MpQD8t3398essAsmwa096v4iDUR8SdkZ9W/Aq7tYbkXybpwPkHWtXNLRHZ6HRHPRsQFETGZ7FPG98rjDT2Y2q3WHWl6G1lX0/jc44CIaKfyMfifqf296bj9GW8eMxvlHPo2UnYBUyS9DSAifkcWnFdI+n0ASS2S5vSyjYPI3ig608DlJRX2UfGa/IjYC6wALpN0kKTDgC+QnfX2Kp0pz0t9+68BvyE7g+/JTWTjBWek6fJ2/lN54JmsHz762M5/T59ujiEbeygP/F6TXsdhabvNkualeR1pm/njcFCq+SVJLcDivl6zjR4OfRspdwObgWclPZ/avgy0AetSt8NPyc7ke/K3QBPZWfs6su6gvO8AZ6Srb66ssP7nyMYRtgL3kwXy8n7Uvh/ZG8QOYDfZAOqf97L8arIB2Gcj4pe59n8P/FzSb9Iyfxm9f4fgXrLjcxfwrYi4M7V/J61/p6RXyI7FCfBGt9RlwAOp++dE4H8AxwIvAT8GVvbjNdsoofRJ08zMCsBn+mZmBeLQNzMrEIe+mVmBOPTNzAqk7m9adeihh8b06dNHugwzs4axfv365yOiudK8ug/96dOn09raOtJlmJk1DEnP9DTP3TtmZgXi0DczKxCHvplZgTj0zcwKxKFvZlYgdX/1zmCs2tDO0jVb2NHZxeTxTSyeM5P5s1v6XtHMbJQbdaG/akM7S1ZuomvPXgDaO7tYsnITgIPfzApv1HXvLF2z5Y3AL+vas5ela7aMUEVmZvWjz9CXNFXSzyQ9ln7S7S9T+yGS1kp6Iv2dkNol6UpJbZIekXRsblvnpOWfkHTOULygHZ1dA2o3MyuS/pzpvw58MSKOBk4ELpR0NHAxcFdEzCD7UYeL0/Knkf1gxAxgEXA1ZG8SZL9sdAJwPHBJ+Y2iliaPbxpQu5lZkfQZ+hGxMyIeTtOvkP3ocgswD7g+LXY9MD9NzwNuiMw6sh+qngTMAdZGxO70u6Frgbk1fTXA4jkzaRq37+9iN40bw+I5vf0Ak5lZMQxoIFfSdGA28HNgYkTsTLOeBSam6RayH2ou257aemqvtJ9FZJ8SmDZt2kBKfGOw1lfvmJm9Vb9DX9KBwI+Az0fEy5LemBcRIalmv7sYEcuAZQClUmnA250/u8Uhb2ZWQb+u3pE0jizwb4yI8o8o70rdNqS/z6X2dmBqbvUpqa2ndjMzGyb9uXpHwA+AxyPib3KzVgPlK3DOAW7PtZ+druI5EXgpdQOtAU6VNCEN4J6a2szMbJj0p3vnA8CngE2SNqa2rwCXAysknQ88A5yZ5t0BnA60Aa8C5wJExG5JXwceSst9LSJ21+RVmJlZvyiiZl3xQ6JUKoV/RMXMrP8krY+IUqV5o+4buWZm1jOHvplZgTj0zcwKxKFvZlYgDn0zswJx6JuZFYhD38ysQBz6ZmYF4tA3MysQh76ZWYE49M3MCsShb2ZWIA59M7MCceibmRWIQ9/MrEAc+mZmBeLQNzMrEIe+mVmBOPTNzAqkz9CXtFzSc5IezbX9UNLG9Hi6/IPpkqZL6srNuya3znGSNklqk3SlJA3NSzIzs56M7ccy1wFXATeUGyLiE+VpSd8GXsot/2REzKqwnauBC4CfA3cAc4GfDLxkMzMbrD7P9CPiPmB3pXnpbP1M4ObetiFpEnBwRKyLiCB7A5k/8HLNzKwa1fbpnwTsiogncm2HS9og6V5JJ6W2FmB7bpntqc3MzIZRf7p3enMW+57l7wSmRcQLko4DVkk6ZqAblbQIWAQwbdq0Kks0M7OyQZ/pSxoLfAz4YbktIl6LiBfS9HrgSeAooB2Yklt9SmqrKCKWRUQpIkrNzc2DLdHMzLqppnvnw8CvIuKNbhtJzZLGpOkjgBnA1ojYCbws6cQ0DnA2cHsV+zYzs0HozyWbNwMPAjMlbZd0fpq1gLcO4J4MPJIu4bwV+GxElAeB/wL4PtBG9gnAV+6YmQ0zZRfT1K9SqRStra0jXYaZWcOQtD4iSpXm+Ru5ZmYF4tA3MysQh76ZWYE49M3MCsShb2ZWIA59M7MCceibmRWIQ9/MrEAc+mZmBeLQNzMrEIe+mVmBOPTNzArEoW9mViAOfTOzAnHom5kViEPfzKxAHPpmZgXi0DczKxCHvplZgTj0zcwKpM/Ql7Rc0nOSHs21XSqpXdLG9Dg9N2+JpDZJWyTNybXPTW1tki6u/UsxM7O+9OdM/zpgboX2KyJiVnrcASDpaGABcExa53uSxkgaA3wXOA04GjgrLWtmZsNobF8LRMR9kqb3c3vzgFsi4jXgKUltwPFpXltEbAWQdEta9rEBV2xmZoNWTZ/+RZIeSd0/E1JbC7Att8z21NZTe0WSFklqldTa0dFRRYlmZpY32NC/GjgSmAXsBL5ds4qAiFgWEaWIKDU3N9dy02ZmhdZn904lEbGrPC3pWuBf0tN2YGpu0SmpjV7azcxsmAzqTF/SpNzTjwLlK3tWAwsk7S/pcGAG8AvgIWCGpMMlvY1ssHf14Ms2M7PB6PNMX9LNwCnAoZK2A5cAp0iaBQTwNPAZgIjYLGkF2QDt68CFEbE3beciYA0wBlgeEZtr/mrMzKxXioiRrqFXpVIpWltbR7oMM7OGIWl9RJQqzfM3cs3MCsShb2ZWIA59M7MCceibmRWIQ9/MrEAc+mZmBeLQNzMrEIe+mVmBOPTNzArEoW9mViAOfTOzAnHom5kViEPfzKxAHPpmZgXi0DczKxCHvplZgTj0zcwKxKFvZlYgDn0zswJx6JuZFUifoS9puaTnJD2aa1sq6VeSHpF0m6TxqX26pC5JG9Pjmtw6x0naJKlN0pWSNDQvyczMetKfM/3rgLnd2tYCfxgRfwT8GliSm/dkRMxKj8/m2q8GLgBmpEf3bZqZ2RAb29cCEXGfpOnd2u7MPV0HnNHbNiRNAg6OiHXp+Q3AfOAnA6zXzGxUW7WhnaVrtrCjs4vJ45tYPGcm82e31Gz7tejTP499w/twSRsk3SvppNTWAmzPLbM9tVUkaZGkVkmtHR0dNSjRzKz+rdrQzpKVm2jv7CKA9s4ulqzcxKoN7TXbR1WhL+mrwOvAjalpJzAtImYDXwBuknTwQLcbEcsiohQRpebm5mpKNDNrGEvXbKFrz9592rr27GXpmi0120ef3Ts9kbQQ+I/AhyIiACLiNeC1NL1e0pPAUUA7MCW3+pTUZmZmyY7OrgG1D8agzvQlzQW+BHwkIl7NtTdLGpOmjyAbsN0aETuBlyWdmK7aORu4verqzcxGkcnjmwbUPhj9uWTzZuBBYKak7ZLOB64CDgLWdrs082TgEUkbgVuBz0bE7jTvL4DvA23Ak3gQ18xsH4vnzKRp3Jh92prGjWHxnJk124dSz0zdKpVK0draOtJlmJkNi1pcvSNpfUSUKs0bdJ++mZnV3vzZLTW9RLM734bBzKxAHPpmZgXi0DczKxCHvplZgTj0zcwKxKFvZlYgDn0zswJx6JuZFYhD38ysQBz6ZmYF4tA3MysQh76ZWYE49M3MCsShb2ZWIA59M7MCceibmRWIQ9/MrEAc+mZmBeLQNzMrkH6FvqTlkp6T9Giu7RBJayU9kf5OSO2SdKWkNkmPSDo2t845afknJJ1T+5djZma96e+Z/nXA3G5tFwN3RcQM4K70HOA0YEZ6LAKuhuxNArgEOAE4Hrik/EZhZmbDo1+hHxH3Abu7Nc8Drk/T1wPzc+03RGYdMF7SJGAOsDYidkfEi8Ba3vpGYmZmQ6iaPv2JEbEzTT8LTEzTLcC23HLbU1tP7W8haZGkVkmtHR0dVZRoZmZ5NRnIjYgAohbbSttbFhGliCg1NzfXarNmZoVXTejvSt02pL/PpfZ2YGpuuSmprad2MzMbJtWE/mqgfAXOOcDtufaz01U8JwIvpW6gNcCpkiakAdxTU5uZmQ2Tsf1ZSNLNwCnAoZK2k12FczmwQtL5wDPAmWnxO4DTgTbgVeBcgIjYLenrwENpua9FRPfBYTMzG0LKuuPrV6lUitbW1pEuw8ysYUhaHxGlSvP8jVwzswJx6JuZFYhD38ysQBz6ZmYF4tA3MysQh76ZWYH06zp9GzqrNrSzdM0WdnR2MXl8E4vnzGT+7Iq3JDIzq5pDfwSt2tDOkpWb6NqzF4D2zi6WrNwE4OA3syHh7p0RtHTNljcCv6xrz16WrtkyQhWZ2Wjn0B9BOzq7BtRuZlYth/4Imjy+aUDtZmbVcuiPoMVzZtI0bsw+bU3jxrB4zswRqsjMRjsP5I6g8mCtr94xs+Hi0B9h82e3OOTNbNi4e8fMrEAc+mZmBeLQNzMrEIe+mVmBOPTNzArEoW9mViCDDn1JMyVtzD1elvR5SZdKas+1n55bZ4mkNklbJM2pzUswM7P+GvR1+hGxBZgFIGkM0A7cBpwLXBER38ovL+loYAFwDDAZ+KmkoyJi3zuOmZnZkKlV986HgCcj4plelpkH3BIRr0XEU0AbcHyN9m9mZv1Qq9BfANyce36RpEckLZc0IbW1ANtyy2xPbW8haZGkVkmtHR0dNSrRzMyqDn1JbwM+AvxTaroaOJKs62cn8O2BbjMilkVEKSJKzc3N1ZZoZmZJLc70TwMejohdABGxKyL2RsTvgGt5swunHZiaW29KajMzs2FSi9A/i1zXjqRJuXkfBR5N06uBBZL2l3Q4MAP4RQ32b2Zm/VTVXTYlvQP4E+AzueZvSpoFBPB0eV5EbJa0AngMeB240FfumJkNr6pCPyJ+C/xet7ZP9bL8ZcBl1ezTzMwGz9/INTMrEIe+mVmBOPTNzArEoW9mViAOfTOzAnHom5kVSFWXbJqZNYJVG9pZumYLOzq7mDy+icVzZjJ/dsVbf416Dn0zG9VWbWhnycpNdO3Jvgva3tnFkpWbAAoZ/O7eMbNRbemaLW8EflnXnr0sXbNlhCoaWQ59MxvVdnR2Dah9tHPom9moNnl804DaRzuHvpmNaovnzKRp3Jh92prGjWHxnJkjVNHI8kCumY1q5cFaX72Tceib2ag3f3ZLYUO+O3fvmJkViEPfzKxAHPpmZgXi0DczKxCHvplZgfjqHRsQ37jKrLFVfaYv6WlJmyRtlNSa2g6RtFbSE+nvhNQuSVdKapP0iKRjq92/DZ/yjavaO7sI3rxx1aoN7SNdmpn1U626dz4YEbMiopSeXwzcFREzgLvSc4DTgBnpsQi4ukb7t2HgG1eZNb6h6tOfB1yfpq8H5ufab4jMOmC8pElDVIPVmG9cZdb4ahH6Adwpab2kRaltYkTsTNPPAhPTdAuwLbfu9tS2D0mLJLVKau3o6KhBiVYLvnGVWeOrRej/cUQcS9Z1c6Gkk/MzIyLI3hj6LSKWRUQpIkrNzc01KNFqwTeuMmt8VV+9ExHt6e9zkm4Djgd2SZoUETtT981zafF2YGpu9SmpzRqAb1xl1viqCn1J7wD2i4hX0vSpwNeA1cA5wOXp7+1pldXARZJuAU4AXsp1A1kD8I2rzBpbtWf6E4HbJJW3dVNE/Kukh4AVks4HngHOTMvfAZwOtAGvAudWuX8zMxuAqkI/IrYC/65C+wvAhyq0B3BhNfs0M7PB820YzMwKxKFvZlYgDn0zswJx6JuZFYjvsmlmA+a7rTYuh76ZDUj5bqvlm++V77YKOPgbgLt3zGxAfLfVxubQN7MB8d1WG5tD38wGxHdbbWwOfRu1Vm1o5wOX383hF/+YD1x+t3/hq0Z8t9XG5oFcG5U82Dh0fLfVxubQt1Gpt8FGh1P1fLfVxuXuHRuVPNhoVpnP9G1Umjy+ifYKAV/Pg43+wpMNB5/p26jUaION5TGI9s4ugjfHIDz4bLXm0LdRaf7sFr7xsffSMr4JAS3jm/jGx95bt2fO/sKTDRd379io1UiDjR6DsOHiM32zOuAvPNlwceib1YFGG4OwxuXuHbM64C882XAZdOhLmgrcAEwEAlgWEd+RdClwAdCRFv1KRNyR1lkCnA/sBf5LRKyponazUaWRxiCscVVzpv868MWIeFjSQcB6SWvTvCsi4lv5hSUdDSwAjgEmAz+VdFRE7HvJgpmZDZlB9+lHxM6IeDhNvwI8DvR2mjIPuCUiXouIp4A24PjB7t/MzAauJgO5kqYDs4Gfp6aLJD0iabmkCamtBdiWW207PbxJSFokqVVSa0dHR6VFzMxsEKoOfUkHAj8CPh8RLwNXA0cCs4CdwLcHus2IWBYRpYgoNTc3V1uimZklVYW+pHFkgX9jRKwEiIhdEbE3In4HXMubXTjtwNTc6lNSm5mZDRNFxOBWlARcD+yOiM/n2idFxM40/V+BEyJigaRjgJvI3gQmA3cBM/oayJXUATwzqCLhUOD5Qa473BqpVmisehupVmisehupVmisequp9bCIqNhNUs3VOx8APgVskrQxtX0FOEvSLLLLOJ8GPgMQEZslrQAeI7vy58L+XLnTU+H9Iak1IkqDXX84NVKt0Fj1NlKt0Fj1NlKt0Fj1DlWtgw79iLgfUIVZd/SyzmXAZYPdp5mZVce3YTAzK5DRHvrLRrqAAWikWqGx6m2kWqGx6m2kWqGx6h2SWgc9kGtmZo1ntJ/pm5lZjkPfzKxAGjL0Jc2VtEVSm6SLK8xfKKlD0sb0+HRu3jmSnkiPcxqg3r259tUjXWta5kxJj0naLOmmXHvdHds+6q2rYyvpilw9v5bUmZtXd8e2j3rr7dhOk/QzSRvSLWJOz81bktbbImnOUNdaTb2Spkvqyh3bawa884hoqAcwBngSOAJ4G/BL4OhuyywErqqw7iHA1vR3QpqeUK/1pnm/qbNjOwPYUD5uwO/X+bGtWG89Httuy38OWF7Px7aneuvx2JINiv55mj4aeDo3/Utgf+DwtJ0xdVzvdODRavbfiGf6xwNtEbE1Iv4NuIXsDp79MQdYGxG7I+JFYC0wd4jqLKum3uHWn1ovAL6bjh8R8Vxqr9dj21O9w22g/w/OAm5O0/V6bPPy9Q63/tQawMFp+p3AjjQ9Enf/rabeqjVi6Pf3bp0fTx+LblX2gy8DWbeWqqkX4ABldxxdJ2n+kFbav1qPAo6S9ECqae4A1q21auqF+ju2AEg6jOys8+6BrltD1dQL9XdsLwX+TNJ2si+Qfm4A69ZaNfUCHJ66fe6VdNJAd96Iod8f/wxMj4g/Ijsrun6E6+lLb/UeFtlXsT8J/K2kI0eiwJyxZF0mp5Cd3V0rafyIVtS73uqtt2NbtgC4NRrnB4Yq1Vtvx/Ys4LqImAKcDvyDpHrOv57q3QlMi4jZwBeAmyQd3Mt23qKeX3RP+rxbZ0S8EBGvpaffB47r77pDoJp6iYj29HcrcA/Z7xaMWK1kZyWrI2JP+jj8a7JQrctjS8/11uOxLVvAvl0l9Xpsy7rXW4/H9nxgRarpQeAAshua1euxrVhv6oZ6IbWvJxsbOGpAex/KAYshGgQZSzaQdThvDoIc022ZSbnpjwLr0vQhwFNkg2ET0vQhdVzvBGD/NH0o8AS9DKYNU61zgetzNW0Dfq+Oj21P9dbdsU3LvZvsRoXKtdXlse2l3ro7tsBPgIVp+j1kfeQi+/nW/EDuVoZ+ILeaepvL9ZENBLcP9P/CkL2wIT5op5OdsT0JfDW1fQ34SJr+BrA5HcyfAe/OrXse2WBNG3BuPdcLvB/YlNo3AefXQa0C/obsbqmbgAV1fmwr1luPxzY9vxS4vMK6dXdse6q3Ho8t2RUwD6SaNgKn5tb9alpvC3BaPRzbnuoFPp6yYiPwMPCnA923b8NgZlYgjdinb2Zmg+TQNzMrEIe+mVmBOPTNzArEoW9mViAOfbNu0p0MHx3A8gslTR7KmsxqxaFvVr2FgEPfGoJD36yysZJulPR4ugne2yUdl25ytV7SGkmTJJ0BlIAb0/3NmyT9taSHJD0qaZkkjfSLMSvzl7PMupE0nexWB38cEQ9IWg48TnaLjHkR0SHpE8CciDhP0j3AX0VEa1r/kIjYnab/AVgREf88Ai/F7C3GjnQBZnVqW0Q8kKb/EfgK8IfA2nTiPobsjoeVfFDSl4C3k903ZzPZnVTNRpxD36yy7h+BXwE2R8T7eltJ0gHA94BSRGyTdCnZHRLN6oL79M0qmyapHPCfBNYBzeU2SeMkHZPmvwIclKbLAf+8pAOBM4arYLP+cOibVbYFuFDS42S3Cv47sgD/X5LKdz58f1r2OuAaSRuB14BrgUeBNcBDw1y3Wa88kGtmViA+0zczKxCHvplZgTj0zcwKxKFvZlYgDn0zswJx6JuZFYhD38ysQP4/sPPd2+hYAB8AAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "dyOVO30kCArA"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}