ml/ppo.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "ppo.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyOPG23g9ISoTlhGTCWrBVV1",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/letianzj/QuantResearch/blob/master/ml/ppo.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JlNkxcEuDTVn",
        "colab_type": "text"
      },
      "source": [
        "## PPO CartPole\n",
        "\n",
        "Use PPO (Proximal Policy Optimization) to solve CartPole game.\n",
        "\n",
        "[OpenAI Gym CartPole](https://gym.openai.com/envs/CartPole-v0/) has four states, cart position, cart speed, pole angle, and pole speed. The actions are either going left or right. The objective is to keep pole from falling. Every move that doesn't lead to a fall gets reward 1.\n",
        "\n",
        "PPO is off-policy TD method. It re-uses old experiences by importance sampling. Inherit notations from [reinforce notebook](./reinforce.ipynb),\n",
        "\n",
        "$$\n",
        "\\begin{aligned}\n",
        "J(\\theta) &= E_{\\tau \\sim \\pi_{\\theta}(\\tau)}[R(\\tau)] \\\\\n",
        " &= \\sum_{t=1}^T E_{(s_t, a_t) \\sim \\pi_{\\theta}(s_t, a_t)}\\left[ r(s_t, a_t) \\right] \\\\\n",
        " &= \\sum_{t=1}^T E_{(s_t, a_t) \\sim \\pi_{\\theta'}(s_t, a_t)}\\left[ \\frac{\\pi_{\\theta}(s_t, a_t)}{\\pi_{\\theta'}(s_t, a_t)} r(s_t, a_t) \\right] \\\\\n",
        " &\\approx \\sum_{t=1}^T E_{(s_t, a_t) \\sim \\pi_{\\theta'}(s_t, a_t)}\\left[ \\frac{\\pi_{\\theta}(a_t|s_t)}{\\pi_{\\theta'}(a_t|s_t)} r(s_t, a_t) \\right] \\\\\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "where last equation assumes $\\pi_{\\theta}(s_t) \\approx \\pi_{\\theta'}(s_t)$, or state $s_t$ has similar probability of occurrence under two policies.\n",
        "\n",
        "PPO improves TRPO (Trust Region Policy Optimization) by moving KL-Divergence condition from constraints to Lagrangian penalties in the objective function. KL-Divergence is a better measure than simple parameter distance.\n",
        "\n",
        "$$\n",
        "J_{PPO}^{\\theta'}(\\theta)=E_{(s_t, a_t) \\sim \\pi_{\\theta'}(s_t, a_t)}\\left[ \\frac{\\pi_{\\theta}(a_t|s_t)}{\\pi_{\\theta'}(a_t|s_t)} A^{\\theta'}(s_t, a_t) \\right] - \\beta KL(\\theta, \\theta')\n",
        "$$\n",
        "\n",
        "PPO2 replaces KL-Divergence by eplison clips to keep policies close.\n",
        "\n",
        "This notebook follows closely [here](https://github.com/dragen1860/Deep-Learning-with-TensorFlow-book/blob/master/ch14-%E5%BC%BA%E5%8C%96%E5%AD%A6%E4%B9%A0/ppo_tf_cartpole.py). Just added some comments for my own understanding.\n",
        "\n",
        "\n",
        "[Reference]\n",
        "* Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.\n",
        "* [RL by David Silver](https://www.davidsilver.uk/teaching/)\n",
        "* [OpenAI Spinning Up](https://spinningup.openai.com/en/latest/algorithms/ppo.html)\n",
        "* [OpenAI Baseline](https://github.com/openai/baselines/tree/master/baselines/ppo2)\n",
        "* [Deep Learning with TensorFlow](https://github.com/dragen1860/Deep-Learning-with-TensorFlow-book)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0lOf7OSWHOvv",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "import os\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import gym\n",
        "from collections import namedtuple\n",
        "import random\n",
        "import tensorflow as tf\n",
        "from tensorflow import keras"
      ],
      "execution_count": 9,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "xtBE-yWBHZZh",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "env = gym.make('CartPole-v1')\n",
        "env.seed(1234)\n",
        "tf.random.set_seed(1234)\n",
        "np.random.seed(1234)"
      ],
      "execution_count": 10,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6aJ__wslg1Vm",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Hyperparameters\n",
        "episodes = 500\n",
        "gamma = 0.98\n",
        "epsilon = 0.2\n",
        "batch_size = 32"
      ],
      "execution_count": 11,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jbLHauFqg2C5",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "Transition = namedtuple('Transition', ['state', 'action', 'a_log_prob', 'reward', 'next_state'])\n",
        "\n",
        "class Actor(keras.Model):\n",
        "    def __init__(self):\n",
        "        super(Actor, self).__init__()\n",
        "        self.fc1 = keras.layers.Dense(100, kernel_initializer='he_normal')\n",
        "        self.fc2 = keras.layers.Dense(2, kernel_initializer='he_normal')\n",
        "\n",
        "    def call(self, inputs):\n",
        "        \"\"\"\n",
        "        input: s, output: pi(a|s)\n",
        "        \"\"\"\n",
        "        x = tf.nn.relu(self.fc1(inputs))\n",
        "        x = self.fc2(x)\n",
        "        x = tf.nn.softmax(x, axis=1)\n",
        "        return x\n",
        "\n",
        "class Critic(keras.Model):\n",
        "    def __init__(self):\n",
        "        super(Critic, self).__init__()\n",
        "        self.fc1 = keras.layers.Dense(100, kernel_initializer='he_normal')\n",
        "        self.fc2 = keras.layers.Dense(1, kernel_initializer='he_normal')\n",
        "\n",
        "    def call(self, inputs):\n",
        "        \"\"\"\n",
        "        input: s, output: v(s)\n",
        "        \"\"\"\n",
        "        x = tf.nn.relu(self.fc1(inputs))\n",
        "        x = self.fc2(x)\n",
        "        return x"
      ],
      "execution_count": 12,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2xfhtwoNhHEu",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "class PPO2():\n",
        "    def __init__(self):\n",
        "        super(PPO2, self).__init__()\n",
        "        self.actor = Actor()\n",
        "        self.critic = Critic()\n",
        "        self.buffer = []              # experience buffer\n",
        "        self.actor_optimizer = keras.optimizers.Adam(1e-3)\n",
        "        self.critic_optimizer = keras.optimizers.Adam(3e-3)\n",
        "\n",
        "    def select_action(self, s):\n",
        "        s = tf.constant(s, dtype=tf.float32)         # [4]\n",
        "        s = tf.expand_dims(s, axis=0)       # [4] => [1,4]\n",
        "        prob = self.actor(s)           # [1, 2]\n",
        "        a = tf.random.categorical(tf.math.log(prob), 1)[0]    # [1, 1] ==> [1]\n",
        "        a = int(a)  # Tensor ==> int\n",
        "        return a, float(prob[0][a])           # a, prob(a|s)\n",
        "\n",
        "    def get_value(self, s):\n",
        "        s = tf.constant(s, dtype=tf.float32)       # [4]\n",
        "        s = tf.expand_dims(s, axis=0)      # [4] => [1,4]\n",
        "        v = self.critic(s)[0]           # [1, 1] ==> [1]\n",
        "        return float(v)                 # v(s)\n",
        "\n",
        "    def store_transition(self, transition):\n",
        "        # store experience\n",
        "        self.buffer.append(transition)\n",
        "\n",
        "    def optimize(self):\n",
        "        # convert sample to tensor\n",
        "        state = tf.constant([t.state for t in self.buffer], dtype=tf.float32)         # (buff, 4)\n",
        "        action = tf.constant([t.action for t in self.buffer], dtype=tf.int32)         # (buff,)\n",
        "        action = tf.reshape(action,[-1,1])                             # (buff,) ==> (buff, 1)\n",
        "        reward = [t.reward for t in self.buffer]\n",
        "        old_action_log_prob = tf.constant([t.a_log_prob for t in self.buffer], dtype=tf.float32)    # (buff,)\n",
        "        old_action_log_prob = tf.reshape(old_action_log_prob, [-1,1])                 # (buff, ) ==> (buff, 1)\n",
        "        # rewards to go, assuming one MC\n",
        "        R = 0\n",
        "        Rs = []\n",
        "        for r in reward[::-1]:\n",
        "            R = r + gamma * R\n",
        "            Rs.insert(0, R)\n",
        "        Rs = tf.constant(Rs, dtype=tf.float32)                # (buff,)\n",
        "        # Reuse buffer roughly 10 times\n",
        "        for _ in range(round(10*len(self.buffer)/batch_size)):\n",
        "            # sample batch size\n",
        "            index = np.random.choice(np.arange(len(self.buffer)), batch_size, replace=False)   # permutation\n",
        "            with tf.GradientTape() as tape1, tf.GradientTape() as tape2:\n",
        "                v_target = tf.expand_dims(tf.gather(Rs, index, axis=0), axis=1)      # [b,1]\n",
        "                v = self.critic(tf.gather(state, index, axis=0))              # baseline, [b,4] ==> [b, 1]\n",
        "                delta = v_target - v                   # advantage, [b, 1]\n",
        "                advantage = tf.stop_gradient(delta)     # freeze advantage in optimizing PPO2 loss, [b, 1]\n",
        "                # retrieve pi(at|st)\n",
        "                a = tf.gather(action, index, axis=0)         # [b, 1]\n",
        "                pi = self.actor(tf.gather(state, index, axis=0))        # [b, 4] ==> [b, 2]\n",
        "                indices = tf.expand_dims(tf.range(a.shape[0]), axis=1)        # [b, 1]\n",
        "                indices = tf.concat([indices, a], axis=1)               # [32, 2]\n",
        "                pi_a = tf.gather_nd(pi, indices)    # pi(at|st), [b]\n",
        "                pi_a = tf.expand_dims(pi_a, axis=1)  # [b]=> [b,1] \n",
        "                # importance sampling\n",
        "                ratio = (pi_a / tf.gather(old_action_log_prob, index, axis=0))   # [b, 1]\n",
        "                surr1 = ratio * advantage                  # [b, 1]\n",
        "                surr2 = tf.clip_by_value(ratio, 1 - epsilon, 1 + epsilon) * advantage   # clipped, [b, 1]\n",
        "                policy_loss = -tf.reduce_mean(tf.minimum(surr1, surr2))   # PPO2 loss function\n",
        "                # v(s) close to MC rewards\n",
        "                value_loss = keras.losses.MSE(v_target, v)             # (b,)\n",
        "            # optimize actor\n",
        "            grads = tape1.gradient(policy_loss, self.actor.trainable_variables)\n",
        "            self.actor_optimizer.apply_gradients(zip(grads, self.actor.trainable_variables))\n",
        "            # optimize critic\n",
        "            grads = tape2.gradient(value_loss, self.critic.trainable_variables)\n",
        "            self.critic_optimizer.apply_gradients(zip(grads, self.critic.trainable_variables))\n",
        "\n",
        "        self.buffer = []  # clear buffer"
      ],
      "execution_count": 13,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zhGLIwiAYsQ-",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 442
        },
        "outputId": "9fc347d1-e6fa-4e23-8e7e-4d3b6dabdcde"
      },
      "source": [
        "agent = PPO2()\n",
        "returns = [] # average returns over 20 episodes\n",
        "total = 0 # total returns\n",
        "for i_epoch in range(episodes):\n",
        "    state = env.reset() # reset to intial state\n",
        "    for t in range(500): # maximum 500 steps\n",
        "        action, action_prob = agent.select_action(state)\n",
        "        next_state, reward, done, _ = env.step(action)\n",
        "\n",
        "        trans = Transition(state, action, action_prob, reward, next_state)\n",
        "        agent.store_transition(trans)\n",
        "\n",
        "        state = next_state\n",
        "        total += reward\n",
        "        if done:\n",
        "            if len(agent.buffer) >= batch_size:\n",
        "                agent.optimize() # train agent\n",
        "            break\n",
        "\n",
        "    if i_epoch % 20 == 0:\n",
        "        returns.append(total/20)\n",
        "        total = 0\n",
        "        print(i_epoch, returns[-1])"
      ],
      "execution_count": 14,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "0 0.6\n",
            "20 38.25\n",
            "40 56.3\n",
            "60 141.85\n",
            "80 264.65\n",
            "100 415.75\n",
            "120 426.4\n",
            "140 442.7\n",
            "160 446.25\n",
            "180 448.15\n",
            "200 472.2\n",
            "220 410.1\n",
            "240 423.1\n",
            "260 473.75\n",
            "280 478.4\n",
            "300 472.65\n",
            "320 412.25\n",
            "340 490.5\n",
            "360 462.05\n",
            "380 390.15\n",
            "400 368.2\n",
            "420 357.0\n",
            "440 469.85\n",
            "460 500.0\n",
            "480 474.75\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Z3Peln-Ru0rl",
        "colab_type": "code",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 279
        },
        "outputId": "e6020309-4a1e-4335-88df-6697edeb16cb"
      },
      "source": [
        "plt.plot(np.arange(len(returns))*20, returns)\n",
        "plt.plot(np.arange(len(returns))*20, returns, 's')\n",
        "plt.xlabel('epoch')\n",
        "plt.ylabel('total return')\n",
        "plt.show()"
      ],
      "execution_count": 16,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXhU5dnH8e+dhQRCSAQStoBhF1RARQRxq4pb3apoXaqotNhqW61Lq63d3rdvW2tba1trtdXWrWrdrUtdqFBFFkHCKktYE7YskBACAZI87x9zogGSzCSZmTOZ+X2ua66c88yZM/cJYe45z2rOOURERFqS5HcAIiIS+5QsREQkKCULEREJSslCRESCUrIQEZGgUvwOoD169uzp8vPz/Q5DRKRDWbBgQZlzLqc1r+nQySI/P5/58+f7HYaISIdiZhta+xpVQ4mISFBKFiIiEpSShYiIBKVkISIiQSlZiIhIUBHtDWVm64EqoA6odc6NNbPuwHNAPrAeuNw5t8PMDHgAOA/YDVznnPskkvGJiPjmvqFQXXJoeUYu3Lk6+vEEEY07iy8458Y458Z6+3cB051zQ4Hp3j7AucBQ7zENeCgKsYmI+KOpRNFSuc/8qIa6CHjc234cuLhR+RMuYA6QbWZ9fIhPREQOEulk4YB3zGyBmU3zyno557Z421uBXt52P6Co0WuLvbIDmNk0M5tvZvNLS0sjFbeIiDQS6RHcJznnNplZLvCuma1o/KRzzplZq1Zfcs49AjwCMHbsWK3cJBIrOlgdvLRORO8snHObvJ8lwMvAOGBbQ/WS97Phr2sT0L/Ry/O8MhHpCDpYHbxfnHP86t8rWjxm+eadUYomdBFLFmaWYWaZDdvAWcBS4DVginfYFOBVb/s14FoLGA9UNqquEhHp8OrrHT9+bRl/mrGGqpTuTR5TTjaXPDSLlxcWRzm6lkWyGqoX8HKgRywpwD+cc/82s4+Bf5rZVGADcLl3/JsEus0WEug6e30EYxMRiaraunq++8JiXlq4iWmnDKLruWsh8Pl4gPqqvYz+xyd857lFLCqq5AdfHEFqsv9D4iKWLJxza4HRTZSXA2c0Ue6AmyMVj4iIX/bW1vGtfyzkneXbuH3SML55+hCsiUQBkJOZxlNfPYFfvrWCRz9cx7LNlTx41bHkdkuPctQH6tBTlIskpA7YkFyxex/ZXTr5HYYvdu+rZdoTC/iwsIwfXzCS6ycODPqa1OQkfnj+SEblZXHXi0s4/w8f8tBXjuW4w5uuuooGJQuRjiZGG5Jr0nqQvrf8kPJSl8WVf57N4zeMo192Zx8i80/lnv1c/7d5FBRV8OvLRjP5uLxWvf6iMf0Y3juTG59cwJcfnsOSrt+k875Df8fR+KKgZCGJIVa/jYcQ1669tawvq2ZdWTXry6r5Vgun21i+m7zDOpOUdFAVRxSu/+rsp9hRvY/pt596QBVL4Zpytj05ny89OIu/Xz+OkX27heX9Yl3Zrr1c8+g8Ckuq+NPVx3LOUW0bY3xE72689s2T+M5zBXRe30SigKh8UVCykMQQo9/GW4rr8j/PZl15NaVVew946lstVF2fct/7ZHRKZnjvTI7o040R3s/jI3z9n27ZyYINO7jniyMOqYufMLgHL3z9RKY8No/LH57Nw9ccx8QhPcPyvjGjiWTcE3jCZfHplAWcMqxVK5geIqtzKn+9diz8T7tO0y5KFiIxyuE4bVgO+T0zGNgzg/weGeT37AI/b/41v7zkaFZsreLTLTt5Y/EW/jF3IwDrI9w2+tScDaSlJDVbzTK8dyYv33wi1z32Mdf9bR73TR7NxcccMkFDx9VM0s2xSnLamSgaHHK3GGVKFiIx6vmvn9jq11wxbsBn2845tu6sYcWWKng2nJEdaNfeWl5ZuIkLRvdtsRG7T1Zn/vn1Cdz45Hxufa6ArTtruPGUQc32CpLYomQhca9s117iqtIjI7f59odGzIw+WZ3pkxXZRuWXF26iel8dV58wIOixWZ1TefyGcdz+z0X88q0VXP3BmWTWbj/0QL/bkuQQShYS1+asLefbzyxknt+BhFMMfYg653h6zgaO7NuNMf2zQ3pNWkoyv7/iGPpkpZP5cROJAvxvS4pVIX5RiAQlC4lL9fWOB98v5P73VpHfI4PalBxS9jQxS3EU/pM1Z+mmSnq5LHKs8tAnwx1XhD5kPtm4gxVbq/jFJUe3qjopKcn4wRdHwsftevvE4+MXBSULiTtlu/bynecK+GB1GReO7svPLzmalLTCA475xZuf8sgHa3nnG6cw1IcYa+vqufulJWzt9Bjv3XYqWZ1TI/uGjT5kdtbs55Rfvc+ovGyeuGFcu0771JyNZKalcNGYvu2NsGPz8Rt/tPg/4YhIGM1ZW855D3zA3HXb+cUlR/PAFWPomnbod6IbTx1Ml9RkfveeP9/U/v7RepZsquQnFxwZ+URxkG7pqdx82hD+u6qU2Wua6bcfgu3V+3hj8RYuObYfXTol+PfOO1dz7YB3OT3zVfhJ5eePGKoybC8lC4kL9fWOP/5nNVf9ZQ5d01J45aaJXDluQLNVI90zOjH1pIG8sWQLyzY3UQ0UQUXbd/Obd1ZxxhG5nHd076i+d4NrJhxO727p/OrtFQSmZWu95+cXsa+unqvHHx7m6DqePfvqmLO2nNOGxc+dxMES/OuAdEhNDIBKAr7sslg16lV+fsnRTd5NHGzqyYP4+0fruf/dVfx1yvERCvZAzjl++OpSzOB/Lj7Kt26j6anJ3HrmUO56aQnvLt/GWUe2LmnV1zv+MW8j4/K7M6xXZtsDiZPqmznrytlXW8+pw8MzpiIWKVmI/1o7FUULA6AeuGJMyB/AWZ1TufHUwdz39koWbtzBMQMOa03UbfL64i3MWFnKj84f6fs8SZOPy+OR/67lvrdXcsaIXiS3YtDXh4VlbCjfzW2ThrUviEb/vi8sKOaO5xfx5NRxnDy0Y33ozlxZSnpqEicM9G+iv0hTspDwasscRC1MRXHf2yuoqqn1HvvZWVPLP1t4+9Z+U7/uxHwe/XAdv313FU9OPaFVr22tyt37+em/ljEqL4spJ+ZH9L1CkZKcxB1nD+empz/h5YWbWjXJ3VNzNtAjoxPnHBW+arQLRvfh3n+v4C8frOtwyWLGyhImDOpBemqy36FEjJKFhFcLH/xvLdnC5soatlbuYUtlDVsqa9haWcOsFk730Iw1ZKan0q1zCplpqWSmh/dPNiMthZtOG8zP3viUOWvLGT+oR1jP39gv3vqUHbv38/gN41r1LT6Szj2qN0f3y+L+d1dxweg+pKUE/7DbUrmH9z7dxo2nDg7p+FClpSQzZcLh/PqdVazcWsXw3u2o3oqi9WXVrC/fzXUx8AUgktTALVHzjac/4X9fX84TszewdFMlnZKD37av+fl5LPrxWXzw3dN585aTee7GCWGP6yvjDyc3M43fvrOqzY29wcxZW86zHxfx1ZMHcmTfrIi8R1uYGd89ZzibKvZ8No9UMM/MK8IBV40LPmK7ta464XDSU5N47MN1YT93pMxcFRi/c9rwjtXO0lq6s5CwKdq+m/4tPP/Gt0+ib1ZnsrukHlhd9JPmXxONBuD01GS+dfoQfvjqMj5YXdbuGUIPVrO/ju+/vIT+3Ttz6xntrOOPgJOG9OTEwT34438KuWxs/xY7B+yvq+fZeRs5dVgO/bt3CXss3TM6cemxeTy/oJg7zh5OTmZa2N8j3GasLCG/Rxfye2b4HUpE6c5C2m1TxR7ufmkJX/j1jBaPO7JvFodldGp/Amiup0w7etBcfnx/+mV35jfvrAz73cWfZqxhbWk1/3fx0XTuFHt12oG7iyMor97Hox+0/I3+veXbKKnay1dOiFx32RtOGsi+2nqemrMhYu8RLjX765i9tpxTw/wFIxbpzkLabGtlDQ++X8izH2/EMK46YQAsbMOJWtt9MgIDndJSkrnljKF898XFvPdpCZNG9grLeVdvq+KhGYV86Zh+Yb9jCacx/bM558je/OWDtVwz4XC6ZzQ9e+zTczfSL7szXzgiclUug3O6csYRuTw1ZwPfOG1wTDcaz1u3nZr99XFfBQVKFtKSZno21XXJ4WcjXuXpuRupr3dcfnx/bv7CkEBX0FVt6DcfI6NcLzm2H3+aUchv3lnJGUfktnv9gPp6x90vLSEjLYV7vjgiTFFGzh1nD+Od5Vv50/uF3HP+yEOeX1dWzYeFZdxx1rCIN9BPPXkgV/1lLq8s3HTAtOuxZsbKUjqlJEW0Y0SsULKQ5jXTsyl5dylPzN7A5GPz+ObpQw6su46RD/62SElO4juThnHLswW8sWQLF4xu33xHz3y8kfkbdnDf5FH06Br7de9DcjO59Ng8npizgetPGnjIOJCn52wgJcm4/PiWWqbCY8KgHhzZtxt//XAdXz6+f8yueTFzVQknDOwek9WL4aY2C2mT6bedyr2TR0WkkdNP54/qy7BeXbn/vVXU1tW37sX3DYWfZH32uPqtUaxPv4rJ758emWAj4NZJw8DBA++tOqC8Zn8dzy8o5uwje5ObGeFl9wi0o3z15IEUluz6rLdRrCnavps1pdUJUQUFShbSjKLtu1t8Pl57fiQnGbdNGsba0mpeKdjcuhc3cydmHWhthn7ZnblmwuG8sKCYwpKqz8pfX7yFyj37uXp89KqEvnh0X3p1S+PRGO1GO+OzLrOx2xYVTkoW8pmN5bt5aMYaLvjDh5z8q/f9Dsc3Zx/Zm6P6deOB6avYV9vKu4s4cNNpg+nSKYVfv/353cVTczYwOCeDCVGsm++UksSUE/P5YHUZK7bujNr7hmrmylLyDuvMoDj94nQwtVkkkmYarHen9uDL3Z5gyabA7Kuj87K4+9wjIEHzhZlx+1nDuf5vH/P8giKuDtJNtHjHbp6dV8QdUYov0np0TeNrJw/i/vdWsaioguQko6Cogh+dPzLqbQdXjRvAH6YX8tcP1vHry0ZH9b1bsre2jo/WlHHJsf1itj0l3JQsEkkz1SFd9peTlGR8/7wjOPeoPp+3Q8yLjxlB2+K0YTkcd/hh/GF6IZcem3dI9826esf7K0p4eu4GZqwqxYA7Yr8NO2TfWnget6SXwKOB/fXpwHvA7OiujZ3dpROXjc3jmXkb+e7Zw8ntFvn2klDMX7+D3fvq4npK8oMpWQgAr9488dDCDtyzqb3MjGcqr6HTvjL4vwOfq07tzqSkv7K5sobczDS+9YUhfHncAPidP7FGQlILc3xF2/UTB/LknA08OWcDt581POrv35QZK0volJzEhMHx32W2gZKFSDM61ZQ1WZ6xfzuD+3flRxeM5IwRvUhN9pr+4mRthlgzsGcGZ47oxVNzNnDTaUNiopvqzFWlHD/wMDJCWDclXiTOlYqEUZPTmSfwnVikffWkgby7fBsvLSwO2oYUaZsr9rBq2y4uOy7y401iiXpDiUjMGzewO0f3y+LRD9dRXx+ZmYFDNWNlYnWZbRDxZGFmyWa20Mxe9/YHmtlcMys0s+fMrJNXnubtF3rP50c6toQTgQn4RKKhYZDe2tJqZqzyd9zKzFUl9M1KZ0huV1/jiLZo3FncAnzaaP9e4H7n3BBgBzDVK58K7PDK7/eOkzDaOm0J+TX/4O+TCuAnlZ8/VH0iB4vBLxbnHd2HPlnp/DXIzLiRtK+2nlmF5Zw6PDdhusw2iGibhZnlAV8k0J/kNgv8dk8HrvIOeZzAagYPARfx+coGLwB/NDNzkVqNJgEVFO0AYHT/bJ8j6SASucE6Br9ApCYncd2J+fzirRUs21zpyyJSCzbsYNfe2oSrgoLIN3D/Dvgu0LA+Yg+gwjlX6+0XA/287X5AEYBzrtbMKr3jD+iSYmbTgGkAAwbE7myUsaigqJLUZGNEn25+h9IxxOAHZqK7YtwAHpi+OrBu+uVjov7+M1eVkpJknJhAXWYbRCxZmNn5QIlzboGZnRau8zrnHgEeARg7dqzuOlphUVEFI/p0i+n1AURakvXHkSxPKoHlHLjCYkZ0BgvOWFnC2PzDyExPjfh7xZpItllMBC40s/XAswSqnx4Ass2sIUnlAZu87U0QWJXTez4LKI9gfAmlrt6xuLiCMaqCko7Mx8GCWytrWLG1KmFmmT1YxJKFc+5u51yecy4fuAL4j3PuagIzDk32DpsCvOptv+bt4z3/H7VXhM+a0l1U76tjdJ6ShUhb/NebZTYRllBtih/jLL5HoLG7kECbhDf7DI8CPbzy24C7fIgtbhUUVQBq3BZpqxmrSujdLZ0jemcGPzgORWUEt3NuBjDD214LjGvimBrgsmjEk4gKiirITE9JmOmUJfHMXFUasW/9tXX1fLC6jPOO6pNwXWYbaAR3glhUVMHovOx2rystEqumPDaPe15Zwu59tcEPbqWFRRVU1dRyagJ2mW2gZJEAavbXsWJrFaP7R79fukhYNTPGxWXk8rWTB/L03I2c+8AHzF+/PaxvO2NlCclJxsQhPcN63o5EEwkmgGWbK6mrd2rclo6vme6xBvwAOHNEL25/fhGXPzybaacM5juThpKW0v6u4jNWlnLcgMPI6px4XWYbKFkkgIUbA43b6jYr8e6EQT34962n8LPXl/PnmWuYsbKE1/feQMqe0kMPDnFsRklVDcs27+TOs2NjLQ2/qBoqASwqrqRvVnrMrDImEkld01L45aWjeOy6sZRX72s6UUDIYzP+uyowiUSidpltoGSRABYVVajLrCSc04/oxTu3ntLu88xYWUJOZhpH9k3saXKULOJc+a69bNy+W1VQkpAOy+jUrtc3dJk9ZWhOwnaZbaBkEecWF1cCGown0pSLHpzFox+uo2RnTZPPLyqupHLP/oScZfZgauCOcwVFFSQZHN1P3WZFDlZbV8//vr6c/3tjOeMH9eCiMX0558g+ZD04EqpLOA5Ynw687D2iNGFhLFKyiHOLiisYmpuZUAvLixyghXVJ3vj2yRSW7OK1RZt5rWAT33txCfe8spTVqf5NWBir9AkSx5xzLCqq4KyRvf0ORcQ/Qe4EhuR25bZJw/jOmUNZsqmSVws2w/woxdaBKFnEsY3bd7Nj9361V4iEwMwYlZfNqLxsJYsmqIE7jn0+06zaK0SkfZQs4lhBUQXpqUkM75WYUyqLSPgoWcSxRUUVHN0vi5Rk/TOLtEozExY2W54A1GYRp/bX1bN0806uHX+436GIdDwJ2j22JfrKGadWbq1iX229GrdFJCyULOLUwiLNNCsi4aNkEacWFVXQI6MTeYd19jsUEYkDShZxqmGm2USf/ExEwkPJIg5V1eynsHSXqqBEJGyULOLQkuJKnNNMsyISPkoWcaig2Bu5naeR2yISHkoWcahgYwX5PbqQ3aV9C7+IiDRQsohDi4or1F4hImGlZBFntlbWsG3nXrVXiEhYKVnEmc9nmlWyEJHwCTo3lJnlAF8D8hsf75y7IXJhSVsVFFWQmmyM7NPN71BEJI6EMpHgq8AHwHtAXWTDkfZaVFTBiD7dSE9N9jsUEYkjoSSLLs6570U8Emm3unrHkk2VfOmYfn6HIiJxJpQ2i9fN7LzWntjM0s1snpktMrNlZvZTr3ygmc01s0Ize87MOnnlad5+ofd8fmvfM9GtKd3Frr21aq8QkbALJVncQiBh7DGznWZWZWY7Q3jdXuB059xoYAxwjpmNB+4F7nfODQF2AFO946cCO7zy+73jpBUKNNOsiERIi8nCzJKAc5xzSc65zs65bs65TOdc0NZTF7DL2031Hg44HXjBK38cuNjbvsjbx3v+DNMseK2yqKiCzLQUBvXM8DsUEYkzLSYL51w98Me2ntzMks2sACgB3gXWABXOuVrvkGKgoYK9H1DkvW8tUAn0aOKc08xsvpnNLy0tbWtocWlRcQWj+meRlKQcKyLhFUo11HQzu7Qt3/Kdc3XOuTFAHjAOOKK152jinI8458Y658bm5OS093Rxo2Z/HSu2VKkKSkQiIpRkcSPwPLC3lW0Wn3HOVQDvAxOAbDNr6IWVB2zytjcB/QG857OA8ta8TyJbtrmS2nrH6DwlCxEJv6DJwmujSHLOdWpNm4WZ5ZhZtrfdGZgEfEogaUz2DptCYBwHwGvePt7z/3HOudZdTuIqKKoE1LgtIpERygjuU5oqd879N8hL+wCPm1kygaT0T+fc62a2HHjWzH4GLAQe9Y5/FHjSzAqB7cAVIV6DEOgJ1Scrndxu6X6HIiJxKJRBeXc22k4n0PawgECvpmY55xYDxzRRvtY7x8HlNcBlIcQjTVhUpJlmRSRygiYL59wFjffNrD/wu4hFJK22vXofG7fv5qoTBvgdiojEqbbMOlsMjAh3INJ2iz5bGU93FiISGaG0WfyBwGA6CCSXMcAnkQxKWqdgYwVJBqO0jKqIREgobRbzG23XAs8452ZFKB5pg0XFFQzNzSQjLZR/ThGR1gvl0yXbOfdA4wIzu+XgMvGHc45FRRVMGtnL71BEJI6FkiymAAcnhuuaKJNoum8oVJdgBPofs9R7ZOTCnat9DU1E4k+zycLMrgSuAgaa2WuNnsokMA5C/FRd0rpyEZF2aOnO4iNgC9AT+E2j8ipgcSSDEhGR2NJssnDObQA2ABPM7HBgqHPuPW/qjs4EkoaIiCSAoOMszOxrBNaXeNgrygNeiWRQIiISW0IZlHczMBHYCeCcWw3kRjIoERGJLaEki73OuX0NO9704ZoN1m8ZzeTr5spFRNohlK6zM83s+0BnM5sE3AT8K7JhSVB3ruas+2fSq1s6T049we9oRCTOhXJn8T2gFFhCYCGkN4F7IhmUBFdSVcOqbbs4cXBPv0MRkQTQ4p2FtxbFMufcEcBfohOShGL2msAighOHHLJMuYhI2LV4Z+GcqwNWmpnmvo4xswrL6JaewpF9NXmgiEReKG0WhwHLzGweUN1Q6Jy7MGJRSVAfrSlnwuAeJCeZ36GISAIIJVn8MOJRSKtsLN9N8Y49fO3kQX6HIiIJIpSV8mZGIxAJ3aw1ZYDaK0QketqyUp74bFZhGbmZaQzO6ep3KCKSIJQsOpj6esfsNeVMHNITM7VXiEh0KFl0MKtKqiiv3seJg1UFJSLR09J6FktoeloPA5xzblTEopJmzSoMjK84cYgG44lI9LTUwH1+1KKQkH1UWEZ+jy70y+7sdygikkCCrWchMaS2rp6567Zz4Zi+fociIgkmlPUsxpvZx2a2y8z2mVmdme2MRnByoEXFlezaW8tEzQclIlEWSgP3H4ErgdUEVsj7KvBgJIOSps32xldMUOO2iERZSL2hnHOFQLJzrs459zfgnMiGJU2ZVVjOiD7d6J7Rye9QRCTBhDLdx24z6wQUmNmvgC2oy23U1eyvY8HGHVw7/nC/QxGRBBTKh/413nHfJDCRYH/gkkgGJYeav34H+2rrmagusyLig1CSxcXOuRrn3E7n3E+dc7ehbrVRN2tNGSlJxriB3f0ORUQSUCjJYkoTZdcFe5GZ9Tez981suZktM7NbvPLuZvauma32fh7mlZuZ/d7MCs1ssZkd26oriXMfrSlnTP9sMtJCqTkUEQmvZpOFmV1pZv8CBprZa40eM4DtIZy7FrjdOTcSGA/cbGYjgbuA6c65ocB0bx/gXGCo95gGPNTWi4o3lXv2s6S4QlN8iIhvWvqa+hGBxuyewG8alVcBi4Od2Dm3xXs9zrkqM/sU6AdcBJzmHfY4MIPAOt8XAU845xwwx8yyzayPd56ENndtOfVOU3yIiH+avbNwzm1wzs1wzk0AVgCZ3qPYOVfbmjcxs3zgGGAu0KtRAtgK9PK2+wFFjV5W7JUdfK5pZjbfzOaXlpa2JowO66M15aSnJnHMgGy/QxGRBBXKCO7LgHnAZcDlwFwzmxzqG5hZV+BF4Fbn3AEjv727iKYmK2yWc+4R59xY59zYnJyc1ry0w5pVWMbx+d1JS0n2OxQRSVChtJbeAxzvnCsBMLMc4D3ghWAvNLNUAoniaefcS17xtobqJTPrA5R45ZsIdMttkOeVJbSSqhpWl+zi0uPy/A5FRBJYKL2hkhoShac8lNdZYGWeR4FPnXO/bfTUa3zew2oK8Gqj8mu9XlHjgUq1V8DsNd6U5GrcFhEfhXJn8W8zext4xtv/MvBWCK+bSGBA3xIzK/DKvg/8EvinmU0FNhCo2gJ4EzgPKAR2A9eHdAVxblZhGd3SUziyb5bfoYhIAguaLJxzd5rZJcBJXtEjzrmXQ3jdhwQWSmrKGU0c74Cbg503kTjnmFVYzoTBPUhO0hKqIuKfUKqT7nXOveScu817vGxm90YjuERXtH0Pmyr2aIoPEfFdKG0Wk5ooOzfcgcihZnlTkp+o9StExGctrcH9DeAmYJCZNR6ElwnMinRgEmivyM1MY3BOht+hiEiCa6nN4h8EGrJ/wedTcgBUOedCme5D2qG+3jF7TTmnDMsh0LFMRMQ/La3BXQlUElglT6Js5bYqyqv3qcusiMQELWIUoz7yxleocVtEYoGSRYz6qLCMgT0z6Jvd2e9QRESULGJRbV09c9dtZ4KqoEQkRihZxKBFxZXs2lvLRHWZFZEYoWQRgz4qDIyv0J2FiMQKJYsY9NGackb26Ub3jE5+hyIiAihZxJya/XUs2LiDiUN0VyEisUPJIsbMX7+DfbX1muJDRGKKkkWMmbWmjJQkY9zA7n6HIiLymVDWs5BIu28oVAfWl/oe8L1OBCZZyciFO1f7GZmICKA7i9hQXdK6chGRKFOyEBGRoJQsREQkKCULEREJSslCRESCUrKIBRm5rSsXEYkydZ2NAWuuW8gZv5nJ9887gmmnDPY7HBGRQ+jOIga8uKCY5CTj4jH9/A5FRKRJShY+q6t3vPTJJk4dlkNut3S/wxERaZKShc8+WlPG1p01TD4uz+9QRESapWThsxcWFJPVOZUzRqgxW0Ril5KFj3bW7OffS7dy4ei+pKUk+x2OiEizlCx89MbiLeytrVcVlIjEPCULH72woJihuV0ZlZfldygiIi1SsvDJ2tJdLNiwg8nH5WFmfocjItKiiCULM3vMzErMbGmjsu5m9q6ZrfZ+HuaVm5n93swKzWyxmR0bqbhixYufFJNk8KVjNLZCRGJfJO8s/lxj3KUAAAp0SURBVA6cc1DZXcB059xQYLq3D3AuMNR7TAMeimBcvtPYChHpaCKWLJxz/wW2H1R8EfC4t/04cHGj8idcwBwg28z6RCo2v81eU86WyhouVcO2iHQQ0W6z6OWc2+JtbwV6edv9gKJGxxV7ZXHphQVFdEtP4cwRvYIfLCISA3xr4HbOOcC19nVmNs3M5pvZ/NLS0ghEFlk7a/bz72VbuXBMX9JTNbZCRDqGaCeLbQ3VS97PhkWmNwH9Gx2X55Udwjn3iHNurHNubE5OTkSDjYQ3F2+hZn89k4/rH/xgEZEYEe1k8RowxdueArzaqPxar1fUeKCyUXVVXHlhQTFDcrsyWmMrRKQDiWTX2WeA2cBwMys2s6nAL4FJZrYaONPbB3gTWAsUAn8BbopUXH5aV1bNfI2tEJEOKGKLHznnrmzmqTOaONYBN0cqlljxksZWiEgHpRHcUVJf73hxQTEnD82hl8ZWiEgHo2QRJbPXlrO5UutWiEjHpGQRJS8sKCYzPYVJIzW2QkQ6HiWLKKiq2c9bS7dw4WiNrRCRjknJIgreXNIwtkJVUCLSMSlZRMELC4oZnJPBmP7ZfociItImShYRtr6smo/X7+BSja0QkQ5MySLCGsZWXHKMqqBEpONSsoig+nrHi59s4qShOfTO0tgKEem4IjaCO6HdNxSqS0gCZgFsBH4CZOTCnav9jExEpE10ZxEJ1SWtKxcRiXFKFiIiEpSSRZh9snGH3yGIiISd2izCpGj7bn719kr+tWgz69WWLSJxRsminapq9vOnGWt49MN1JBl8+4yhXqu2iEj8ULJoo9q6ev45v5jfvruSsl37uOSYftxx9nD6ZneGgtymG7MzcqMfqIhIGChZBON1gz3YTsvm+3v+xPH5h/HYdcczKq/RVB7qHisicUbJIphmurt2dxU8dPWxnHNUb03jISJxT8miHc49uo/fIYiIRIW6zoqISFBKFiIiEpSSRQuWb97pdwgiIjFByaIZJTtrmPr4x5TTzIJF6gYrIglEDdxN2L2vlqmPz6dyz3623LiYHv2y/A5JRMRXShYHqa93fOe5ApZuruQv14zlKCUKERFVQx3s3rdX8PaybdzzxZGcObKX3+GIiMQEJYtGnp23kYdnruXqEwZww8R8v8MREYkZShaejwrLuOeVpZw8tCc/ufBIjcoWEWlEyQIoLNnF159awMCeGTx49bGkJuvXIiLSWMJ/Km6v3scNf/+YTilJPHbd8XRLT/U7JBGRmJPQvaH21tYx7Yn5bN1Zw7PTxtO/exe/QxIRiUkxlSzM7BzgASAZ+Ktz7pdhfYODphtPA14AajJ6kj5gTVjfSkQknsRMNZSZJQMPAucCI4ErzWxkWN+kmenG0/eWhfVtRETiTcwkC2AcUOicW+uc2wc8C1zkc0wiIkJsJYt+QFGj/WKv7ABmNs3M5pvZ/NLS0qgFJyKSyGIpWYTEOfeIc26sc25sTk6O3+GIiCSEWEoWm4D+jfbzvDIREfFZLCWLj4GhZjbQzDoBVwCvhfUdmptWXNONi4i0KGa6zjrnas3sm8DbBLrOPuacWxbWN7lzdVhPJyKSKGImWQA4594E3vQ7DhEROVAsVUOJiEiMUrIQEZGglCxERCQoJQsREQnKnHN+x9BmZlYKbGjjy3sCiTwpVCJffyJfOyT29evaAw53zrVqVHOHThbtYWbznXNj/Y7DL4l8/Yl87ZDY169rb/u1qxpKRESCUrIQEZGgEjlZPOJ3AD5L5OtP5GuHxL5+XXsbJWybhYiIhC6R7yxERCREShYiIhJUQiYLMzvHzFaaWaGZ3eV3POFmZo+ZWYmZLW1U1t3M3jWz1d7Pw7xyM7Pfe7+LxWZ2rH+Rt5+Z9Tez981suZktM7NbvPJEuf50M5tnZou86/+pVz7QzOZ61/mctwwAZpbm7Rd6z+f7GX84mFmymS00s9e9/US69vVmtsTMCsxsvlcWlr/9hEsWZpYMPAicC4wErjSzkf5GFXZ/B845qOwuYLpzbigw3duHwO9hqPeYBjwUpRgjpRa43Tk3EhgP3Oz9+ybK9e8FTnfOjQbGAOeY2XjgXuB+59wQYAcw1Tt+KrDDK7/fO66juwX4tNF+Il07wBecc2MajakIz9++cy6hHsAE4O1G+3cDd/sdVwSuMx9Y2mh/JdDH2+4DrPS2HwaubOq4eHgArwKTEvH6gS7AJ8AJBEbupnjln/0fILB+zARvO8U7zvyOvR3XnOd9IJ4OvA5Yoly7dx3rgZ4HlYXlbz/h7iyAfkBRo/1iryze9XLObfG2twK9vO24/X141QrHAHNJoOv3qmEKgBLgXWANUOGcq/UOaXyNn12/93wl0CO6EYfV74DvAvXefg8S59oBHPCOmS0ws2leWVj+9mNq8SOJDuecM7O47jNtZl2BF4FbnXM7zeyz5+L9+p1zdcAYM8sGXgaO8DmkqDCz84ES59wCMzvN73h8cpJzbpOZ5QLvmtmKxk+2528/Ee8sNgH9G+3neWXxbpuZ9QHwfpZ45XH3+zCzVAKJ4mnn3EteccJcfwPnXAXwPoGql2wza/hy2PgaP7t+7/ksoDzKoYbLROBCM1sPPEugKuoBEuPaAXDObfJ+lhD4ojCOMP3tJ2Ky+BgY6vWQ6ARcAbzmc0zR8BowxdueQqAuv6H8Wq9nxHigstEta4djgVuIR4FPnXO/bfRUolx/jndHgZl1JtBe8ymBpDHZO+zg62/4vUwG/uO8CuyOxjl3t3MuzzmXT+D/9X+cc1eTANcOYGYZZpbZsA2cBSwlXH/7fjfI+NQIdB6wikBd7g/8jicC1/cMsAXYT6AeciqButjpwGrgPaC7d6wR6B22BlgCjPU7/nZe+0kE6m0XAwXe47wEuv5RwELv+pcCP/LKBwHzgELgeSDNK0/39gu95wf5fQ1h+j2cBryeSNfuXeci77Gs4bMtXH/7mu5DRESCSsRqKBERaSUlCxERCUrJQkREglKyEBGRoJQsREQkKCULEZ+Y2WkNM6OKxDolCxERCUrJQiQIM/uKt0ZEgZk97E3Ut8vM7vfWjJhuZjnesWPMbI63PsDLjdYOGGJm73nrTHxiZoO903c1sxfMbIWZPW2NJ7ESiSFKFiItMLMRwJeBic65MUAdcDWQAcx3zh0JzAR+7L3kCeB7zrlRBEbFNpQ/DTzoAutMnEhghD0EZsW9lcDaKoMIzG8kEnM066xIy84AjgM+9r70dyYwEVs98Jx3zFPAS2aWBWQ752Z65Y8Dz3vz9fRzzr0M4JyrAfDON885V+ztFxBYh+TDyF+WSOsoWYi0zIDHnXN3H1Bo9sODjmvrvDl7G23Xof+TEqNUDSXSsunAZG99gIb1jA8n8H+nYSbTq4APnXOVwA4zO9krvwaY6ZyrAorN7GLvHGlm1iWqVyHSTvoWI9IC59xyM7uHwOpjSQRm8r0ZqAbGec+VEGjXgMAU0H/2ksFa4Hqv/BrgYTP7H+8cl0XxMkTaTbPOirSBme1yznX1Ow6RaFE1lIiIBKU7CxERCUp3FiIiEpSShYiIBKVkISIiQSlZiIhIUEoWIiIS1P8D3f83QDFf2k0AAAAASUVORK5CYII=\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Wt-Iln8jKQ29",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}