cart_pole.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.random as npr\n",
    "from scipy.linalg import solve_discrete_are"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.set_printoptions(precision=3, suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "hz = 60.0\n",
    "_sat = 5.0\n",
    "def f(s, a):\n",
    "    th, x, thdot, xdot = s\n",
    "    sdot = np.array([\n",
    "                thdot, xdot,\n",
    "                9.81*np.sin(th)/1.0+_sat*a[0]*np.cos(th)/1.0-0.2*thdot,\n",
    "                _sat*a[0]-0.1*xdot])\n",
    "    return s + sdot/hz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wrap2Pi(x):\n",
    "    th = np.fmod(x+np.pi, 2.0 * np.pi)\n",
    "    if th < 0:\n",
    "        th = th + 2.0 * np.pi\n",
    "    return th - np.pi\n",
    "\n",
    "Q = np.diag([10., .2, 0., 0.05, 0., 0., 0.])\n",
    "Qf = np.diag([20., 0.2, 0., 0., 0., 0.,0.])\n",
    "R = np.diag([0.01])\n",
    "\n",
    "Rinv = np.linalg.inv(R)\n",
    "\n",
    "def l(x): #loss in z space\n",
    "    x[0] = wrap2Pi(x[0])\n",
    "    return np.dot(x, np.dot(Q, x))\n",
    "\n",
    "def ldx(x):\n",
    "    x[0] = wrap2Pi(x[0])\n",
    "    return 2.*np.dot(Q, x)\n",
    "\n",
    "def m(x): #loss in z space\n",
    "    x[0] = wrap2Pi(x[0])\n",
    "    return np.dot(x, np.dot(Qf, x))\n",
    "\n",
    "def mdx(x):\n",
    "    x[0] = wrap2Pi(x[0])\n",
    "    return 2.*np.dot(Qf, x)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "state_dim = 4\n",
    "action_dim = 1\n",
    "\n",
    "num_x_obs = 7\n",
    "num_u_obs = 2\n",
    "num_obs = num_x_obs + num_u_obs\n",
    "\n",
    "def z(s):\n",
    "    th, x, thdot, xdot = s\n",
    "    return np.array([th, x, thdot, xdot, np.sin(th), np.cos(th), 1.])\n",
    "\n",
    "def v(s, a):\n",
    "    th, x, thdot, xdot = s\n",
    "    return np.array([a[0], np.cos(th) * a[0]])\n",
    "\n",
    "def dvdu(s, a):\n",
    "    return np.array([[1.0, np.cos(s[0])]]).T\n",
    "def dvdz(s, a):\n",
    "    th, x, thdot, xdot = s\n",
    "    out = np.zeros((num_u_obs, num_x_obs))\n",
    "    out[:,0] = np.array([0., -np.sin(th) * a[0]])\n",
    "    return out "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "npr.seed(0) ## set the seed for reproducibility \n",
    "num_trials = 60 ## number of resets\n",
    "horizon    = 100 ## how long we simulate the system for\n",
    "M = num_trials * horizon ## M sized data\n",
    "\n",
    "_A = np.zeros((num_obs, num_obs))\n",
    "_G = np.zeros((num_obs, num_obs))\n",
    "cnt = 0\n",
    "s_bnds = np.array([6.24, 1.0, 2.0, 2.0])\n",
    "a_bnds = np.array([1.])\n",
    "## loop to collect data\n",
    "for k in range(num_trials):\n",
    "#     x_t = npr.uniform(low=-s_bnds, high=s_bnds)\n",
    "    x_t = npr.normal(s_bnds*0., s_bnds)\n",
    "    for t in range(horizon):\n",
    "        u_t = npr.normal(0.*a_bnds, a_bnds)\n",
    "        x_tpo = f(x_t, u_t)\n",
    "        u_tpo = u_t\n",
    "\n",
    "        z1 = np.concatenate([z(x_t), v(x_t, u_t)])\n",
    "        z2 = np.concatenate([z(x_tpo), v(x_tpo, u_t)])\n",
    "        cnt += 1\n",
    "        _A += np.outer(z2, z1)/cnt\n",
    "        _G += np.outer(z1, z1)/cnt\n",
    "        \n",
    "        # reset for next loop\n",
    "        x_t = x_tpo\n",
    "        u_t = u_tpo\n",
    "\n",
    "        \n",
    "# _A /= M\n",
    "# _G /= M\n",
    "\n",
    "K = np.dot(_A, np.linalg.pinv(_G))\n",
    "A = K[:num_x_obs,:num_x_obs]\n",
    "B = K[:num_x_obs, num_x_obs:]\n",
    "\n",
    "# A[np.abs(A)<1e-3] = 0\n",
    "# B[np.abs(B)<1e-3] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = int(1.5*hz)\n",
    "u = [np.zeros(1) for t in range(T)]\n",
    "def update(x, u, A, B):\n",
    "    zt = z(x.copy())\n",
    "    zot = []\n",
    "    for t in range(T):\n",
    "        zt = z(zt[:state_dim])\n",
    "        zot.append(zt.copy())\n",
    "        zt = np.dot(A, zt) + np.dot(B, v(zt[:state_dim], u[t]))\n",
    "    rho = np.zeros_like(zt)\n",
    "    for t in reversed(range(T)):\n",
    "        Bdz = np.dot(B, dvdz(zot[t][:state_dim], u[t]))\n",
    "        rho = ldx(zot[t].copy()) + np.dot((A+Bdz).T, rho)\n",
    "        Beff = np.dot(B, dvdu(zot[t], u[t]))\n",
    "#         u[t] = np.clip(-Rinv.dot(Beff.T.dot(rho)), -1., 1.)\n",
    "        du = Beff.T.dot(rho) + 2.0 * R.dot(u[t])\n",
    "        u[t] = np.clip(u[t] - 0.1 * du, -1, 1)\n",
    "    return u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 3.14   0.    -0.083  0.083]\n",
      "[ 3.687 -1.084  0.243 -2.484]\n",
      "[ 3.53  -4.612 -0.403 -5.955]\n",
      "[  3.075 -10.117  -1.177  -6.16 ]\n",
      "[  2.163 -13.39    0.54   -1.666]\n",
      "[  4.034 -14.145   2.269  -1.077]\n",
      "[  2.43  -13.35   -4.608   3.01 ]\n",
      "[ 1.728 -9.284  3.533  6.77 ]\n",
      "[ 5.696 -4.774  2.39   4.687]\n",
      "[ 6.422 -1.554  0.083  2.878]\n",
      "[6.464 0.192 0.097 1.326]\n",
      "[ 6.396  0.543 -0.124 -0.249]\n",
      "[ 6.368  0.09  -0.055 -0.923]\n",
      "[ 6.348 -0.888 -0.041 -1.445]\n",
      "[ 6.3   -2.223  0.06  -1.544]\n",
      "[ 6.29  -3.609 -0.117 -1.792]\n",
      "[ 6.309 -4.947 -0.064 -1.67 ]\n",
      "[ 6.229 -6.374 -0.064 -1.573]\n",
      "[ 6.226 -7.371 -0.121 -1.008]\n",
      "[ 6.128 -7.872 -0.088 -0.008]\n",
      "[ 6.162 -7.284 -0.002  1.187]\n",
      "[ 6.175 -5.938  0.035  2.034]\n",
      "[ 6.198 -3.997 -0.028  2.548]\n",
      "[ 6.209 -1.654  0.028  3.04 ]\n",
      "[6.267 1.013 0.207 3.398]\n",
      "[6.407 3.59  0.031 2.485]\n",
      "[ 6.415  5.147 -0.093  1.152]\n",
      "[ 6.384  5.698 -0.007  0.186]\n",
      "[ 6.378  5.533 -0.051 -0.626]\n",
      "[ 6.39   4.761 -0.079 -1.411]\n",
      "[ 6.371  3.346 -0.049 -2.042]\n",
      "[ 6.364  1.473 -0.029 -2.513]\n",
      "[ 6.335 -0.774  0.045 -2.749]\n",
      "[ 6.319 -3.165 -0.033 -2.918]\n",
      "[ 6.259 -5.593  0.021 -2.703]\n",
      "[ 6.319 -7.704 -0.554 -3.187]\n",
      "[ 6.107 -9.786  0.063 -1.447]\n",
      "[  6.132 -10.46   -0.127  -0.306]\n",
      "[  6.151 -10.067   0.053   1.121]\n",
      "[ 6.194 -8.787  0.131  1.974]\n",
      "[ 6.18  -7.028  0.011  2.461]\n",
      "[ 6.197 -4.723  0.051  3.064]\n",
      "[ 6.209 -2.029  0.099  3.517]\n",
      "[6.372 0.987 0.109 3.284]\n",
      "[6.417 3.201 0.019 1.957]\n",
      "[6.376 4.299 0.043 0.895]\n",
      "[6.374 4.685 0.047 0.128]\n",
      "[ 6.385  4.456  0.028 -0.649]\n"
     ]
    }
   ],
   "source": [
    "xt = np.array([3.14, 0., 0., 0.])\n",
    "\n",
    "trajectory = []\n",
    "\n",
    "for t in range(int(40*hz)):\n",
    "    u[:-1] = u[1:]\n",
    "    u[-1] = np.zeros_like(u[-1])\n",
    "    u = update(xt, u, A, B)\n",
    "\n",
    "    u1 = u[0]\n",
    "    xtpo = f(xt, u1)\n",
    "\n",
    "    z1 = np.concatenate([z(xt), v(xt, u1)])\n",
    "    z2 = np.concatenate([z(xtpo), v(xtpo, u1)])\n",
    "    cnt +=1 \n",
    "    _A += np.outer(z2, z1)/cnt\n",
    "    _G += np.outer(z1, z1)/cnt\n",
    "    \n",
    "    K = np.dot(_A, np.linalg.pinv(_G))\n",
    "    A = K[:num_x_obs,:num_x_obs]\n",
    "    B = K[:num_x_obs, num_x_obs:]\n",
    "    \n",
    "    xt = xtpo\n",
    "    \n",
    "    trajectory.append(xt.copy())\n",
    "    if t % 50 == 0:\n",
    "        print(xt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEGCAYAAACkQqisAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hc1Zn48e+ZombJsiXZlots2ZIr7gg3DBhjTDMBQgglhQQSspCEJSHJkk2yS7LJL5tGChtISEiDACGUhBI6uBtccO+SJVuukiWr15k5vz/OjC3bkiyN5s69M/N+nkfPSKPRve/VSPPOae9RWmuEEEKIcLjsDkAIIUTskiQihBAibJJEhBBChE2SiBBCiLBJEhFCCBE2j90BREJOTo7Oz8+3OwwhhIgpGzZsOK61HtSXY8RFEsnPz2f9+vV2hyGEEDFFKbW/r8eQ7iwhhBBhkyQihBAibJJEhBBChE2SiBBCiLBJEhFCCBE2SSJCCCHCJklECCFE2OJinYjoXm1zO+XVTRyra6GpzU9zu5+Wdj/NbX5afQHSkz0M7p/M4IwUhg1IYVhmKi6XsjtsIWJOuz/A4ZpmXEqRluQmLclDsscV1/9PkkTilNaaN7Yf5dGlJWw+WNurn03xuhg2IJWMZA/9gh/pyR6S3C68HoXX7SItyW2+lxR8TOjrZA/9kt30S/KQm5mC1+0iENAcb2gFBenJHlK9bpTq+p8qENA0tPlQmMcD3T4+FhxvaGVV8XG2H65j77F6SiobqW1up6HVh9Yaj9tFktuFx63wuFwkuRWe4Ndel4skj4v+qZ7gY1x43er034mGVl+AVp+fgNa4XS48LoXbpc64Dd7v7uL+0NduRZrXTWaal/RkL163ed49LnVarF63IjXJQ/8U8zfS3fPk8wcoq2pk++E6CgalM3l4ZhR+89bQWrPnWAPv7qpg77F6qpvaKK9uYn9VE77A6Xs0KQXpSR4yUjxkpHiDt6c+97pdtPsD+Pyadn8ADQS0xh/QBLSmtT1AU5sflwsCAfAFAuZvxGM+ZuVn8fmLx9jzi0CSSFw63tDKPU9+yNqyagoG9eP+y8cxdkg6uZmppCe7SfGaj1Svm2SPi/oWHxX1rVTUt1Be3cy+ygYO1zbT0OqnsdVHdWMTjW0+2n2aNn+ANl+ApjYfgXPsZ5bkcZGR7KGxzUdLe+Dk/R6XIjs9iRSvm8xUL0op6prbcSkIaDhc00yr79Tjhw9IZW5BNgWD0hk7OJ2xQ9LJTPXS7te0+vy0+QK0+gJk90ticP8Uq36tXWrzBWhu99PU5qO+xUd9Szt1LT5ONLax8UAN7++rYm9Fw8nfyZicfkzPG8DANC/pKR4UivaAeRHx+QO0BW99AfOi0h78ndc2t1Pn95kXnOALTMeX7CSPmxSvC7dS+AJ+/AGNL6DxBwLBW43Przu/P3jrP9eT2o0Ur4uc9GR8fs2IgamMz83A41LUt/rYe6yBPcfqT3ter5kylG9ePYERA9PCPmc0aK1paQ9Q1djKvspGVhYf560dxyg93gjAsMwUstKTGDs4gysn5zIqux8KaGrz09jmo7nNH/y7MH8b9S0+jje0UXq8kfoWH23+QIek7EIpcCmFWylcLkWyx0Wq1432g1uZNxn+gKamqY1WX4AxOf1s/f2oeNjZsKioSEvZE+NEYxsf+81qDtU0819LzuPjRSPwuCM/9KW1ptUXoKHVR1Or39y2+YK3fupb2impbKSx1UeK101+dhooRWOrj9rmdqob2mhu91PX0k5AQ0aKB7R515bbP4XB/ZMBqGv2setoPZsP1lBZ33rOuIZmpjA9bwBjh2TQP8VDWpKH1CQXqV4Pbf4A1Q2tVDe1m3/A9uALdIcX6jZ/gNb2U7cdk1SbP3DyhVZr8Otzv+j2S3Jzfn4Wc8ZkMa8ghynDM3E7uGtD61NJpanNb1pLLb4zktypd81t/gDNbeZ5rKxv5XhDG26Xoriigf1VjQQ0pCW5KRyczoTcDMbn9mdCbgbv7Kzg0WXFBAKwYPwglkwbxmUTBgOw+1g9+yobcSnon+LF63HR3GZaWCleF5mpSSS5XTS3+9lX2UD5iSZqmkzibm03XbRtwVZZa/C5CwQ06Ske3C6FS6mTybeh1by4t7T78WtNoMNzG2oNtPtPf469bsXcghwWTxrC4klDbHnjEilKqQ1a66I+HUOSSHz5whPreXdXBU/cOZs5Y7LtDiei6lra2XusgZKKBhrbfHjcLpLdLpK9pnvlcG0Lm8tr2FReQ/mJJrr601bBF6cUrwtvsGsm1DWQ1OF4yR73yfuTPa6T3Tlul3mH2PGdYorXTVqSm4xgt05Gipf+KR7yc/rhtSCJx4PDNc38bsU+/rX1CMfqzv0GoStulyIz1fy+U4Kt62SPm2Sv6+TnSkF9i4+ANkkiEPzj6JfsISPZQ2qS+2SCMR+cfJ6T3C5Sk0yreVR2GpOHZ9I/xRupX4OtJIkESRIx3t9XxS2Pvc83rhzPPQsK7Q7HVu1+04/c0u6nqc10NSW5XWT1S2JAWpKjWwOJJhDQrN9/gpXFx0n2uBg7OJ3Cwem4lKK2uZ12f4AUrxuPW9EcbB35/Bqvx8Xo7H6MGCgTQcIViSQiYyJx5JGlJQzOSOaOC0fbHYrtvG4XmakuMlPj4x1jPHO5FLNGZzFrdJbdoYgwSBKJE4dqmlmxt5J7F44lxeu2OxxhhbrDULkbAj5IyYQh50GSvYOqQkgSiRMvbz6M1nDjzBF2hyIiKRCAzU/Dut/D4Q9P/54nBSZeC5f8B+SMtSc+kfAkicSJ5XsqmZCbwchsZ0+XFL1Qfwz+fjscWANDJsOi78Lw803yaKyE4rdh699hx0tw6Tdh7pfBLf/SIrrkLy4ONLX5WF92gs9cmG93KCJSGirhT1ebLqzrHoHpt5lpZR1NuBou+Qa8ej+8/SAUvwM3PwGpA20JWSQmmXsYBz4orabNH+CisTl2hyIiIeCH5++A2kPwyRdgxifOTiAhGblw85Nw/aNQ/gH86VpoPB7deEVCc2wSUUqVKaW2KqU2KaVk/m43Nh2oQSmYOVLegcaFZT+C0uVwzU9h1NxzP14p01K57W9QtReevBHam62PUwgcnESCLtVaT+/rPOZ4t+1QLQWD0umXLL2TMe/wJlj2Y5h2K8z4ZO9+tmAh3PQnOLIZXvoyXa62FCKCnJ5ERA9sPVTLlBguZieCAgEzvtFvEFz5v+EdY/xVsPDbZsB97WORjU+ITjg5iWjgTaXUBqXUXWd+Uyl1l1JqvVJqfWVlpQ3hOUNFXQsV9a2SROLBtufg0Hq4/HuQOiD841x0P4y9wgy2V++LWHhCdMbJSWS+1nomcBXwRaXUxR2/qbV+TGtdpLUuGjRokD0ROsDuY/UATBrW3+ZIRJ8E/KYba/B5MPXmvh1LKVjyc3B54KV7TQtHCIs4thNda30oeFuhlHoRmAUstzcq5ymragJgtM3loEUfbXvBDIp//C/gisB7u8zhpkXzyn2ma2taHxOTsFZDBZS8B8f3mK/TB0PeLBg6veuZeQ7hyCSilOoHuLTW9cHPFwPfszksR9p/vJEUr4vBGcl2hyLCFQjA8p/A4Ekw4drIHXfm7bD+cXjv+3De9eCRvxHHqSmHd75nujJ1AFAmaehg63HYTFj4LShcZGuY3XFkEgGGAC8Gd0nzAE9prV+3NyRnKqtqYlRWv5jf+S+hlS6F47vhhsci0woJcbnMKvcnPwrr/wBz7o7csUXfaG0mPrz9XZMw5n4RptwEQ6aY563uMOx+DVY/bKZsT74RlvwCUpzXbe3IJKK13gdMszuOWHCgupH8bOnKimlrfw9pOaa1EGkFC2H0JWa8ZfptpnCjsJffZ6Zgb34KCi+HJQ/BgJGnP6b/MLjgTjPNe+XPzfNXVQyfeB7SnTUG7OSB9YRworGN//7nNh5dWkKgl1uTBgKa/VVN5Mt4SOyqKYc9r8HMT1vT3aQULHoQmqth9f9F/viidwIBePlek0AWfBM+8fezE0hHnmRY8ADc+jRU7jGtytaG6MXbA5JEbPb15zbz5zX7+dHru/jj6rJe/eyx+hZafQFGZknRxZi18QlzW3SHdecYPtNU+137W2itt+484tze/m/Y9FeTQBY80PNB83FXmLpox7bD3z9jWjMOIUnERnuO1fP2zgq+tngcFxZm8+jSEtp8PZ+OuT84M0u6s2KU1rDlbzD6YhiQZ+25LrwPWmrhw79Yex7Rtc3PwOpfQdGdpnx/b429HK75GRS/Be88GPHwwiVJxEbPf3gQj0tx66yRfHbeaI43tLK6pOfF8/ZXNQIwSsq/x6aD6+BEWd/XhfTEiCIYNR/W/Br87dafT5yu9iC8+jXzHFz1o/Cn7RZ91iSh1Q+bKcEOIEnERquKj3NBfhbZ6cnMH5tDvyQ3b2w/2uOfL6tqwutWDM1MsTBKYZktz5q9QSYsic755t8HdYdg2/PROZ8wtIaX7wPth+t/De4+btm8+PuQPRZe/ndHFNqUJNIFM2jdyM4jdTS0Rr7/sa6lnR2H607uK53idTOvMIc1JVU9PsaBqibyBqbhccvTGHMCftj+Ioy7MnrTNgsXQc44qakVbdtfNF1Qix6Egfl9P15SmqlIULMfVjzU9+P1kSOn+NopEND89YP9PLq0hMO1LQC4FFyQn8WN54/guunDSPb0fQ/zDftPENAwO5hEAGblZ/HWjmNU1LcwOOPcrYuyqkbZyTBWHVgDTcetmdbbFaXggs/Da1+HQx+aAXdhLV+rqWE2ZDJc8LnIHXf0RTDl47DqFzDtFsguiNyxe0newnbQ0u7n7r9u4Dv/3M6IrDR++NEp/N9tM7hnQSHHG1r5xnNbuPB/3+O3y0poaff36VxrS6vxuBQzOuwBUpRvPl9fduKcP691cHqvDKrHpp2vgDvZrBOIpmm3QFK62bNdWG/jk6bFcPl3wdX3N5+nWfx90x36xn9G9ri9JC2RIK01X312E29sP8a3r5nInfNHn1wFvmQq3L94HKuKq/jt8hJ++Nou/ry6jPsXj+ejM4eHtVp8bWk1U0Zkkpp06g9r8vBMktwuNh+s4eopQ7v9+erGNhpafTK9NxZpDTtfhsLLIDk9uudO6W8G8jc+CZf/D/TLju75E0nAbwbAh58PBZdF/vgZQ8yOljnjIn/sXpCWSNAjS0v419ajfPOqCXzuojFnJQalFPPH5vDEnbN56vOzGZSRzP1/38wnfv/ByVlSPdXc5mfLwRpm5Weddr/X7WJcbjo7Dted8xihwov5OZJEYs7hjVB30KzdsMOsz4O/FTY9ac/5E8WuV+BEKcy717oiihOXwCBJIrZbtqeSn725m49MG8ZdF4855+PnFeTw4j0X8v9umMLWg7Vc8YvlPPH+fnQPd5LbWH6Cdr8+Oaje0aSh/dlxuO6cxzo1vVe6s2LO3jdBucyguh0GT4S8OaY1IrsfWmfVr2DgaPveLERJQieR3UfrueZXK7jjT+sYNySDH350So+7plwuxW2zR/LWVy9h9uhsvvOPbXzzha09SiRrS6tRCoryO08iVY1tVNa3dnuM/VVNKAUjBqb2KF7hICXvwrAZkHb28x81Mz5pyo6Xr7Uvhnh2ZLPZYGzO3ZEfC3GYhE4iSR4XQ/qn8Jl5+fztrrlh7VGem5nCHz9zAXcvKOCZdeU8tfbAOX9mbWk1E3P7k5l69nzxScNMgbztR7rv0tpf1ciwzNSIzBQTUdRcAwfXm8KIdjrvBvD2g42ygt0SG580EyemftzuSCyX0ElkdE4//vCZC/jOkklkpoW/AMjlUnzjivHMGZPFz97cQ31L1yuCW9r9fHjgRKddWQAThmYAnHNcpLSqSTaiikVlK8yiM7uTSHI6TL4Btr3ouIJ+Ma+9xSwknXgtpA489+NjXEInkUhSSvGtqydR3djG71aUdvm4NSVVtLQHWDC+83LO/VO8jBiYys5uWiJaa0orGySJxKKS98wU2xEX2B0JzPg0tDfCjn/YHUl82fUKtNTAzE/ZHUlUSBKJoCkjMrlqci5/XFlKbVPnrZE3dxwlPdnD3IKup1ZOHNqfXUe7rrZa3dhGXYtPSsDHopJ3If+ivpe+iIS8WWbgd8uzdkcSXzb9FTJHQv7FdkcSFZJEIuzey8ZS3+rj8VVnt0bqW9p5ZfMRFk0c3O1YxsTcDPZVNnS5oLEsODNrjCSR2FK9z0z5tLsrK0Qps5te6XKo73nNNtGNpmrYtwym3BjZXSodLDGuMoomDu3fZWvkb+vKqW/1ccf80ec8RkDD3mOd91XvqzRJRLqzYkyo6mrBpfbG0dHUjwNaijJGyu5/mTGvSdfZHUnUSBKxQGetkcZWH79dvo85Y7KYOmJAtz8/YagpyNfVuEhZVSMel5LpvbFm33uQmQfZhXZHckrOWBg6Xbq0ImXHP81OhUOn2x1J1EgSsUDH1kh1YxsAv36vmMr6Vr5+xYRz/vyorDRSvW52Hu08iZRUNDIyS6r3xpRAAMpWmg2orFq9HK6pH4cjm+D4XrsjiW0ttaa1OfEjznuOLSSvQha5b9E4Wn0B7vrLeh5dWsKjy0r42PkjOH/Uuaf8uVyK8bkZXbZEdh2tY+LQKJUPF5FxfDc0n4BR8+yO5GyTbwQUbP273ZHEttIVEGiH8VfbHUlUSRKxyPjcDB66eRrbDtfyo9d3MXdMNt+77rwe/3xohtaZK+AbWn2UVTUxMbieRMSI/avNrROTSEauaSFteVbKoPTFvqVmAacTpm9HkVTxtdCSqcOYOyab4w1tjBuS3qtqvxOHZvD02gMcrWthaOapsY9dwdbJpGHSEokpB9ZAeq6ZUutEUz8O//wiHNpgttIVvbfvPci/EDxJdkcSVdISsVh2ejLjczN6XS4+1F2168jp60V2hJLI0MzIBCisp7VpiYya69y+8onXgssrCw/DVVMOVcUwxkEz76JEkohDjc8Nlj85Y1xkU3kNOelJDOmfbEdYIhw1B8ze5qMutDuSrqVkmqnHO/4pXVrh2LfU3I5ZYGMQ9pAk4lD9U7yMzEpjy8Ga0+5fV1ZN0aissDbCEjY5sMbcjpxrbxznMuk6k/CObLI7ktizbymkDzFl9hOMJBEHmz06i/f3VRMImHeGh2uaKa9u7rJ4o3Co/avMO/3Bk+yOpHvjrwblNq0R0XOBgEkiYxY4t7vSQpJEHOzCwhxqm9tPdmmtLD4OIEkk1uxfY1ohTi+DkZZlZmlJl1bvVGyHpuMJ2ZUFkkQcbV5BNkrBmzuOAfD6tqMMH5DKeTIzK3Y0VELVXud3ZYVMus7U+Dq23e5IYkcCj4eAJBFHG9w/hfmFOTy/4SBlxxtZvqeSJVOHynhILDng4PUhnZmwBFCw82W7I4kdpcshZxz0H2Z3JLaQJOJwn70wn0M1zVzzqxW4XIrb5+XbHZLojf1rwJMaO7WU0geZEvF7XrM7ktigNRxcB3mz7Y7ENpJEHG7hhCHce9lYhg5I5aGPT2PYACm6GFMOrDaL92JpAdq4K80e4XWH7Y7E+aqKTTmbvFl2R2IbSSIx4KuXj+Ptr17CkqmJ2VyOWS11cHRr7HRlhYy/ytzuecPeOGJB+VpzO0KSiBAi0srXgg7EzqB6yKAJMGAU7Hnd7kic7+BaM307Z5zdkdhGkogQVjmwGlye2OvqUMp0ae1bCm1NdkfjbOXrYHiR86dvWyhxr1ycra3RdL9U74P2FrujiX37V8PQaZAUgztQjr8SfC1QuszuSJyrpQ4qdsTem4QIkyq+wqxleOM/YdtzpvsFwJsGYxfDnLth5Bx744tFbU1wcD3MvcfuSMIzaj4kZZgurdAYiTjdoQ2ATrjS72eSJJLoag/Bn64xM3Fm321mEvlazAvg9hdNVdeJH4FrHjLTP0XPHFxnNigaNd/uSMLjSYLChWZwXeuELOdxTgfXASrhS+dLEklkfh8891loPA6feRXyOryjmn4bLP4fWPMILP8JlF8ItzwNI863L95YUrYSlCu2W3HjrjIlUI5sgmEz7I7GecrXmkkIKYm9LYNjx0SUUlcqpXYrpYqVUg/YHU9cWvpDKP8Arv3F6QkkJKkfXPJ1+Py74EmBPy85VeJBdK9spVlgmBLDJWoKF5nb4rftjcOJAoHgIsPE7soChyYRpZQb+DVwFTAJuFUp5fASqDHm2A5Y+RBM/yRM+Vj3j82dDJ97GwbmwzOfCPYFiy61N8Oh9WaXu1iWPsi0QPZKEjlLVTG01CT0+pAQRyYRYBZQrLXep7VuA54Brov4WY5th19Og7/eZGYkJZJ3vmsGThf/T88enz4YPvkCpGXDkx+DE2WWhhfTDq4DfxvkX2R3JH1XuMishWg+YXckznIwuMgwwWdmgXOTyHCgvMPXB4P3naSUuksptV4ptb6ysjK8s3hTzRzv8g/gj1dDY1XYAceUAx+YWTcXfcWU/+6p/kPhUy9CwA9//wz4Wi0LMabFw3hISOHlZsaedGOerjy4yDB7rN2R2M6pSeSctNaPaa2LtNZFgwaFOWsoawx87HG4/WVoqoJ/3R/ZIJ1q9a8gdSDM+kLvfza7AK5/BA5vhDe/E/nY4kHZSsidGh8DrsPPh5QB0qV1poPrzNTeBF5kGOLU38AhIK/D1yOC91lj6DS46H4zpbVspWWncYTqUtj1KhTdAUlp4R1j4hKY+yVY+1spGX6m9uD06PwYndp7JrfH7L1e/LZsVBXSUgsVO2U8JMipSWQdMFYpNVoplQTcArxk6Rnn3Qv9R8DrD5jumni17vfgcsMFn+/bcRY9aJLvq18z/1TCOLgO/K3xMR4SUng5NByFY9vsjsQZQosMZWYW4NAkorX2AV8C3gB2As9qra3dai0pDS7/rin7sfXvlp7KNn4fbHnW1EXqP7Rvx3J74dpfQmMFvPO9yMQXD8pWAio+xkNCCi8zt3vfsjcOpygPLjIcntiLDEMcmUQAtNb/0lqP01oXaK1/EJWTTr4RBk+Clb+Iz6b7vqXmRX/qzZE53rAZZlxl3ePBfyxh1odMhdQBdkcSORm5kDsFit+xOxJnOLgWBk+M7TVAEeTYJGILpeDCf4fKnbD3Tbujibwtz5hB0nFXRO6YC79lXmRefyA+E29vtLeY7qxYLXXSncJFUP6+KTqYyEKLDBO8XlZHkkTONPlGSM+FDX+yO5LIam+B3a/BpOvAkxy54yZnwMJvm8V121+M3HFj0cG1ZjxkdByNh4QUXg4Bn1T1rdprxgBlfchJkkTO5PbC1JtMS6QhzPUnTlS6DNoaTDHFSJt2Kww+zyxgTOS1I6UrzPqQWNvJsCfyZkFyfxkXkZ0MzyJJpDPTbjPvurY9Z3ckkbPzZbNC3Yp3yS63Wfl+oszM/kpUpcvNOFE8rA85k9sLYy6Rqb4H15ou4exCuyNxDEkinRkyCQZNNOsp4kHAb7qyxi2ObFdWR4WXQcFCWPbjxJzy29YYrJcVh11ZIYWXQ90hqNxldyT2KZdFhmeS30RXxl9pdqaLh5pB5R9A03GYcI2157nsv01RurWPWXseJzqwxrReR19sdyTWCVX1TdQurZZak0BlPOQ0kkS6Mu4q0P74mNa453Vwec07SSsNm27WoKz5NbTWW3supyldYX7H8bQ+5EyZw80U+EQtDX9wPbKT4dkkiXRlRJGpL1Xyrt2R9N2+peYPPxrz2i/+hmm9JdrYSOly8zcTi/up90bhZabV1dpgdyTRF9rJcLhszNaRJJGuuNww6sLYr6XVVA1HtphB0WgYcb7p9lj9sBknSAQttWb3v3geDwkpvNyUuS9dbnck0Ve+1rTEZJHhaSSJdGfUhVCzH2oP2h1J+MpWABrGLIjeOS/5D1MVef0fondOO+1fbcqlx/N4SMjIueDtl3hdWgG/aYmMnG13JI4jSaQ7oZ3pylbZG0df7FsGSenRbYLnzTJJa9Uvoa0peue1S+lycCcnRl+5Jyk41fetxJrqW7EDWutMEhWnkSTSnSGTITkT9sdwl1bpMrP4ze2N7nkv/gY0VsLGJ6J7XjuUvGt+x94UuyOJjsJFUHPAbBGbKA68b27zpCVyJkki3XG5zWDpoQ/tjiQ8dYfNP/roKI2HdJR/IYycZ1oj8byKvfagmfYZqnSbCBJxqu+B9yFjGAwYaXckjhNWEgnu8ZEYhs80G9DEYrdM6N2TXWU4Lv6aWZy28Ul7zh8NJe+Z24IESiIDR0HOONOllSgOvG+mbytldySOc84kopRaqpTK7/D1LMymUYlh2AyzXiQWN+Q5uA48qaaMtx0KFprm//KfmgKQ8ajkHcgYakqDJ5LCRWZCQXuz3ZFYr6Yc6g7G9xqgPuhJS+SHwOtKqXuUUj8AfgN81tqwHGTYDHMbi11a5WtN/NEeDwlRCi79FtQfjr+qyGBm7JS8Z5Jlor1DLVgIvhaTSOJd+QfmVpJIp86ZRLTWbwD/BvwSuAO4Wmsdg6+oYcoYCulD4PBGuyPpnfYWOLLZ/i08R19s9tdY8bPY7BLszuFNpsxLwUK7I4m+UfPAnRQfi3HP5cAaU7x08Hl2R+JIPenO+g7wMHAx8CCwVCllcREmB1EKhk43i8liyZFNEGi3v2S1UmbjqsYKWP+4vbFEWsk7gIIxl9odSfQl9TPTXUNjQvHswPtmgo3bY3ckjtST7qxsYJbWeo3W+rfAFcB91oblMEPOM7OcfG12R9JzoX0PnFAsbtQ880K78ufxVS6j+B1TL6xftt2R2KNgIVRsh/qjdkdinZZaOLZd1od0oyfdWfdprZs7fL1fa21xJT+HGTzJVGiNpXnxB9fCwHxIH2x3JMal/xlfq9hbas3EhUTsygoJXXs8t0bK1wFaxkO6IetEeiI086Zih71x9JTWwX0PHNAKCcmbZcZH1vw6PtaNlC43s/YSaWrvmYZMhn6Dgt16cerAGlDB9WKiU5JEeiJnHLg8Zr1ILKgth4ajzujK6mj+V01cm56yO5K+K37HlJNx2u84mlwu001Z8h4EAnZHY40Da2Do1PivztwHkkR6wpNktsOMlZbI4eAkgGEz7Y3jTGMWmCnHq34Jfp/d0YRPa/Pue/TF9k2fdorCy8yGZ8e22h1J5LU3my7L/Pl2R382qhoAABwsSURBVOJovU4iwfUiNyulEmuqwuCJsZNEjmw2TfAhk+yO5HRKmdbIiVLY8Q+7owlf9T5TOyqRx0NCxiwwt/E41bd8rSl7n58A1Zn7IJyWiALmAy9EOBZnGzwJTpTFxh4ZR7fAoAngTbU7krNNWGK6B1f+PHarwO55w9yOTaz5JZ3KyDVjI/GYRMpWmjdjMqjerV4nEa31r7XWX9Zaf8SKgBwrZ5y5rSqxN46eOLIZhk6zO4rOuVww/yumjMzeN+2OJjx7XoNBE83sNwEFl5q1FLHwBqs3ylaYKdyyCVW3epRElFLXKKW+oZT6r9CH1YE5TnahuXX6NN/6o9BwzLlJBGDKTZCZBysesjuS3mupNaU+xl1hdyTOUXCZ6faJ5X13ztTWZPZUl/GQc+rJivXfADcDX8Z0Zd0EjLI4LufJGmNund4SObLZ3Do5ibi9MO/LUP5+7NVeKn7HrBkaf5XdkTjHyLngSYmvLq3yD0zFh0TY8riPetISmae1/jRwQmv9XWAuMM7asBwoKQ36D3d+S+TIZkBB7mS7I+nejE9BWk7stUb2vA6pWYmxi2FPeVPMVtLxlERkPKTHepJEQqvVm5RSw4B2YKh1ITlYdkFsJJHsQkjOsDuS7iWlwZy7zZ4UR7bYHU3P+H1mHGfsYrNhmTilYCEc32026YoHZSvNdHSn/x85QE+SyCtKqQHAT4APgTLgaSuDcqzsQqja6+xZRU4eVD/TBZ8z1VFX/tzuSHrm4DpoPgHjr7Q7Euc5WQIlDlojbY1waIOMh/RQT5LIj7XWNVrr5zFjIROA71sblkNlF5qB1aZquyPpXFO1Wa0eK0kkdQBccKdZM+L0sSYws7Jc3sQuddKVwRPNtgnxkERC4yGjZTykJ3qSRNaEPtFat2qtazvel1CcPkPraHDVsNPHQzqac495YV71C7sj6Z7WsPMV8+5UpnyeTSnTGtm31GzWFctKV5jxkDwZD+mJLpOIUipXKXU+kKqUmqGUmhn8WACkRS1CJ3F6EqncZW5jafOcjCEw81Ow6WmoO2x3NF2r3AXVJTDxWrsjca6Chaa7L9b23jlT2UoYPhOS0+2OJCZ01xK5AvgpMAJ4CPhZ8OMrwH9aH5oDDRhpCjFW7bU7ks5V7DAzh5xS/r2n5t0LaFjp4NbIzpcBBRMSZz+2XhuzAFBQHMNdWq0NcPhDmdrbC10mEa31n7XWlwKf0Vpf2uHjOq11YpU8CXF7TSI5UWZ3JJ2r2GX6pmNtv++Bo2DarWYfdqe2Rna+DHmzTZkP0bl+OWY8LpZLw5e/b9YByaB6j/VkTGSVUupxpdRrAEqpSUqpOy2Oy7kGjHJmEtHalKoP7X0Say7+utmfw4nrRk6UmXpkE5fYHYnzFV5mChe21NodSXjKVprehrzZdkcSM3qSRP4IvAEMC369h0TbHrejgflwYr/dUZyt7jC01prCi7Fo4CiY/gn48M/OW2uw8xVzO0GSyDkVXGbeDJQutzuS8JSugOHny3hIL/QkieRorZ8FAgBaax9g2fQLpdSDSqlDSqlNwY+rrTpXWAaOguZqaKmzO5LTVQY3zBrssPLvvXHx10yLymmtkR3/gNwpkDXa7kicL2+WWftTHINdWs01ZjxktJR+742eJJFGpVQ2oAGUUnMAq9uqP9daTw9+/Mvic/VOqHJrjcNaI6FdF2O1OwvMeNOMT8KHfzH7dTjBiTKzyHDyjXZHEhvcXvMiXPKOsxfldqZsBeiA2a1R9FhPkshXgZeAAqXUKuAvmGKMiWlAsPak08ZFKnZB+hBIy7I7kr65+GtmYsDyn9gdibEtOIfkvI/aG0csKVxo3gTEwgLSjkreNVseS120XjlnEtFafwhcAswDvgCcp7W2utjRl5RSW5RSf1BKDbT4XL0Taok4bVykYkdst0JCMkdA0R2w8a/OeBHa9jyMmGW6MUXPnCyBEmNdWiXvmVlZniS7I4kpPd2UahYwDZgJ3KqU+nRfTqqUelspta2Tj+uAR4ECYDpwBLM2pbNj3KWUWq+UWl9ZWdmXcHondSAk93dWd1YgYBbDDYqDJAJmC113Eiz7sb1xVOwym2dN+Zi9ccSarDEwcHRsjYtUl5ptm2XL41475z7pSqknMC/qmzg1oK4x3Vph0Vov6snjlFK/A17p4hiPAY8BFBUVRa/zVSnzrtRJ3Vm1B6C9KT5aImBWsRd9FtY+Bgu/DQPy7Ilj2/OgXDDpenvOH8sKL4NNT4GvFTzJdkdzbvveM7cyHtJrPWmJFAEXaq3vCW6L+2Wt9b1WBaSU6lhm/gZgm1XnCtuAUc7qzoqHQfUzzbnHDMy+/6g959catj1nVi5nDLEnhlhWcJl5Y3Pgfbsj6ZmS98x+QTlj7Y4k5vQkiWwDorlM98dKqa1KqS3ApZgyK84yMN90Zzll9knFDnMbq2tEOjMgz8yI+vDPZupltB3eCNX7pCsrXKMvMov2YmFcJOCH0mVmr/hYq/bgAN0VYHxZKfUSkAPsUEq9oZR6KfRhVUBa609pradoradqrT+itT5i1bnCNjAffC1mL3MnqNhl9iyPt+qyF94LbQ2w/g/RP/e25011YSm4GJ7kDFMFNxbqaB3eaFbYS1dWWLobE/lp1KKINSdnaJU5o5ZSxc74aoWE5E4xA50f/AbmfjF6fet+n0kihYvMRAoRnsKF8M73oP6oM/5PurL3LUBJEglTdwUYl2mtlwFXhz7veF/0QnSgk2tFHDAu4vfB8T3xNR7S0bx7TYtvy9+id87it6H+CMz4RPTOGY/GBXeA3PO6vXGcS/FbMKII+mXbHUlM6smYyOWd3HdVpAOJKaHZQrUOWFV9ohT8rfGbRMYsgNypsPr/zFTmaNj4BPQbdOpFUIRn8CTzhmuXs4pOnKbxOBz6EAo7e5kTPdHdmMjdSqmtwPjgwr/QRylg9WJDZ/OmQlqOMwoFhgbV4zWJKGVaI8d3m3eMVqs/Brtfg+m3mRIeInxKwfirzW6HbY12R9O5kncBDWN7tOpAdKK7lshTwLWYkifXdvg4X2v9ySjE5mwD8qCm3O4ozKA6CnLG2x2Jdc67HvqPgNUPW3+uzU+bKrQz+rSeVoRMuNq0lEveszuSzu1907whHDrD7khiVndjIrVa6zKt9a1a6/0dPqqjGaBjZY5wTktkYD4kxfGOxW4vzLnbFMg79KF15wn4YcMfYdSFkFNo3XkSyci5kJIJux3YpRXwm1X1hYvA1dPiHeJM8psLV+ZIqC23f61I5a7YLv/eUzM/bcrNrPqldefY/S8z4272F6w7R6Jxe2HsFaaL0O+zO5rTHd5otnUYK+MhfSFJJFyZI8yK3OYT9sXga4WqYhgch9N7z5TSH2Z9Hnb8M9iFZ4E1j5hy9LL5VGRNXGJerMsdtnp975umrI3Uy+oTSSLhCs3QsnPfi6pisx90IrREAOZ8EbxpsNyCwoyHNsCB1TDrC+ByR/74iazgMnAnn9oh0il2vwbDi2J/+wSbSRIJV+YIc2vnuEioZlY8LjTsTL9s09W07XnTFRFJKx4yffczZUA94pLTzbv9Xa/Y3/0bcmI/HN0iFQkiQJJIuDJHmttaG2doVe4C5U6sonHzv2Jm07z5nci9IB3bbl7gZt8df6VjnGLiEvO/ctQhqwN2BVtFE6Xrsq8kiYQrLQs8qfa2RI7vNWXpY6HUdqSk9IcFD5iZWluejcwxl//E7GgnA+rWGXelGX9wSpfWzpdhyGSz94noE0ki4VIquFbEzjGREshOwKmoRXdA3mz419eh9lDfjnVoA2x/EWb/m/SNW6lfDoycd6oFYKeGClOiXrqyIkKSSF/YuVYkEDAD69kJ1JUV4nLD9Y9CoB1e+lL45VC0hje+bUqczL8vsjGKs024xqxrsnvb412vAlqSSIRIEumLzDz7xkTqD4OvGbIL7Dm/3bILYPH/mLIVKzrdQfncdr5kZmQt+KYpXS6sNeEac7vrVXvj2Pmy6cZKlFmNFpMk0heZedBYCe3N0T/38b3mNpEG1c9UdCdMvRne+75ZP9IbTdXw6tdgyBSYebs18YnTDRxlimna2aXVfMJsQDVhiWxAFSGSRPriZDXfPvbLh6Oq2Nwm4phIiFJw7S9hxAXwwhd6VxLltf8wC+CufwTc3W2rIyJqwhIoX2sKXdph58tmbdV519tz/jgkSaQvTq4VsaFLq6oYvP0gY+i5HxvPvKlwy9NmXOOZ26Du8Ll/ZsuzsPVZuPjrMHSq9TGKUyYuAbR9tbS2Pme6sobNtOf8cUiSSF9khloiNiWR7AJpkgOkD4LbnoHWenj6Vmht6PqxhzbAS182RRYvuj96MQpj8CQYONqeLq36o1C6HKbcJP83ESRJpC/6DzNz3+2YoXWiDLJGR/+8TjXkPLjxcbOY7ckbzZ7ZZ6rcDX+9CdKHwE1/lv1C7KCUGWDftwxa6qJ77u0vAhomfyy6541zkkT6wu013UnR3lckEDDnHDAyuud1uvFXwk1/Mq2NPy2BgxvM/VrDjpfg8cvNCv9PvWhaL8IeE68107P3vhnd8279uxnYHzQuuueNczKi2FeZI6LfndVYYTb6Ce31Lk6ZdB3ckgov3gW/X2i6HP1tZp/23Clw85Nm/xVhnxEXmNI1e96AKVFqFVTsMm8uFv8gOudLIJJE+iozDw6tj+45Qy0faYl0btxiuG8rfPiEKdTo9sKoeWY6sHRh2c/lhjELzFRbraMzPrHxCXB5YNot1p8rwUgS6avMEWaNQiAQvd3RavYHz50XnfPFouQMmHuP3VGIroxZANueC26qNtHac/naYPMzMP4qU35FRJSMifTVgDzTv9sQxXnvoXpdAySJiBg1ZoG53bfU+nPteR2ajsMMKfNvBUkifWVHSfjackjNklIdInYNyDPrNaKRRNb93rTaZQdDS0gS6Ss7djisOSCtEBH7xiyAspXgb7fuHJW7zdhL0R1SmcAikkT6yo5V6zK9V8SDMQugraF35Wp6a+3vwJ0kO1ZaSJJIXyVnQMqA6K0V0TrYEpHpvSLG5V8EKOu6tFrqYPPTMPlGGVC3kCSRSBgQxZLwjcdNCXiZmSViXVoWDJsOJe9Yc/zNz5iWzqzPW3N8AUgSiYzMkdErfVIbmpkl3VkiDoy70lT1baiM7HG1hrWPwfDzzYewjCSRSBiQZ7qztLb+XDWSREQcGX8VoGHvG5E97r6lULUXLpBWiNUkiURCZh601UNLjfXnkjUiIp7kToX+I2BXhEvDv/8IpGXDeTdE9rjiLJJEIuHkNN8ojIvUlENKpvkQItYpZVojJe9GbofQo9tMccfZd4M3JTLHFF2SJBIJ0dxXpObAqQWOQsSDCVebySKRmqW18ueQlA6zPheZ44luSRKJhND4RFRaIgfMXtVCxItR803LetvzfT9WdSlsfwGKPgupA/t+PHFOkkQiIS0bPKnWt0ROrhGRloiII54ks9vgjpeg+UTfjrX6YVOtd84XIxObOCdJIpGgVHT2FWmqhvZGSSIi/sz4lNkjZ8uz4R+j5gBsfBKm3Qr9h0YuNtEtSSKREprma6WasuC5JImIODNsOoyYBat+ZUq3h+Pd75s3dJd8I7KxiW7ZkkSUUjcppbYrpQJKqaIzvvdNpVSxUmq3UuoKO+ILS2YUVq3LGhERzxY8AHUHYeNfev+zhzfClr/BnHtO1bMTUWFXS2Qb8FFgecc7lVKTgFuA84ArgUeUUu7ohxeGAXnQWBm5aYqdkSQi4lnBQhg5D979gSnv01OBALzxLbPl7vyvWBef6JQtSURrvVNrvbuTb10HPKO1btValwLFwKzoRhemk/uKWFj+pOaAKfYoa0REPFIKlvwcWuvh5X/veQWIdb+H/avgsu9ASn9rYxRncdqYyHCgY5/QweB9Z1FK3aWUWq+UWl9ZGeG6O+GIxr4iMjNLxLvBE2DRg7DrFfjgt+d+fFUJvPVfULgIZt5udXSiE5YlEaXU20qpbZ18XBeJ42utH9NaF2mtiwYNGhSJQ/ZNNBYcShIRiWDuF2HcVfDmt82mVV1pa4JnbwdPMnzkYdOSEVFnWRLRWi/SWk/u5OOf3fzYIaBjUagRwfucL2MoKLd1M7RkHxGRKJSCGx412+c+fRsc3Xr2YwIBeOlLcGwb3Pg49B8W/TgF4LzurJeAW5RSyUqp0cBYYK3NMfWM22P+kK0aE2mqgvYmaYmIxJA6ED75PCT1gz9cBct+DCf2m610D2+CJz9qVrhf9l8wdpHd0SY0WzYdVkrdADwMDAJeVUpt0lpfobXerpR6FtgB+IAvaq39dsQYFiun+Z4oM7eSRESiGJAHd74Jr9wH7/3AfIQk94drHoIL7rQvPgHYlES01i8CL3bxvR8AP+jse443MN+6rT6rSsxt1hhrji+EEw3IMy2S6lJTmbelzqxGH3cV9Mu2OzqBTUkkbmUXwOanoK3RNMMjqboEUCZRCZFoskbD7C/YHYXohNPGRGJbdqG5DbUaIqmqxHSXyf4IQggHkSQSSSeTyN7IH7u6BLKlK0sI4SySRCIpu8DcRrolojVU7YOsgsgeVwgh+kiSSCR5U02X0/EIt0SaqqC19lSSEkIIh5AkEmnZBVBVHNljnpyZJUlECOEskkQiLXusSSI9LR7XE9XBJCItESGEw0gSibTsQmitM2XhI6WqBJRLSp4IIRxHkkik5QRnaFV2Vuk+TMf3mPUhnqTIHVMIISJAkkikDT7P3B7bHrljVu6CQRMjdzwhhIgQSSKRlpELadmmumgk+FpNd9bgCZE5nhBCRJAkkUhTCoZMjlwSqSoG7ZeWiBDCkSSJWCF3ClTsBL+v78eq2GluB0sSEUI4jyQRKwyZDL6WU1Nz+6Jyl9nsKmds348lhBARJknECrmTzW0kurQqdpry757kvh9LCCEiTJKIFXLGg8trdmDrq8pdMqguhHAsSSJW8CTBsOlwcF3fjtNab2Zm5U6LTFxCCBFhkkSskjcbDn0Ivrbwj3FkC6BhqCQRIYQzSRKxSt5s8LfC0S3hH+NIsDts2PTIxCSEEBEmScQqebPM7YH3wz/G4U2QMQzSB0cmJiGEiDBJIlbJyDX1rvavCv8YhzdKK0QI4WiSRKxUcBnsW2ZKl/RWc41ZrT5UkogQwrkkiVhp7GJob4T9q3v/swfXARpGzol4WEIIESmSRKw0+mJwJ8OeN3r/s/tXg8sDI4oiH5cQQkSIJBErJaXB2Mth2/Pgb+/dzx5YY6b2JvWzJjYhhIgASSJWm/4JaKyAvW/1/Gfam80ak5FzrYtLCCEiQJKI1cZeDhlDYfXDPd93vXSFWWNSsNDa2IQQoo8kiVjN7YWL7ocDq2H7Cz37mb1vgDcNRl1obWxCCNFHkkSiYebtMOIC+OeX4K3/hj1vdt0q0dp8f8wC8KZEM0ohhOg1SSLR4EmCm/9qxjhW/RKeuglef6Dzxx7aALUHYPxV0Y1RCCHC4LE7gISRMQQ+9QK0t8Bb/wUf/AZGzYNJ153+uE1PgSf17PuFEMKBpCUSbd4UuOIHkDsVXnsA2hpPfa+lFrY+BxOXQEqmfTEKIUQPSRKxg9sLV/8E6g/Dyp+fuv/930BrLcz7sn2xCSFEL0gSscvIOTDlJlj1K6jcY9aFrPgZTPyI7B8ihIgZMiZip0XfheK34bcXQcBv1pNc8zO7oxJCiB6TJGKnzOHw+XdhzSPgSTbdWLJ3iBAihkgSsVvWGLjmp3ZHIYQQYZExESGEEGGTJCKEECJstiQRpdRNSqntSqmAUqqow/35SqlmpdSm4Mdv7IhPCCFEz9g1JrIN+Cjw206+V6K1lj1hhRAiBtiSRLTWOwGUUnacXgghRIQ4cUxktFJqo1JqmVLqoq4epJS6Sym1Xim1vrKyMprxCSGECLKsJaKUehvI7eRb39Ja/7OLHzsCjNRaVymlzgf+oZQ6T2tdd+YDtdaPAY8BFBUV9XC3JyGEEJFkWRLRWi8K42dagdbg5xuUUiXAOGB9hMMTQggRAY5abKiUGgRUa639SqkxwFhg37l+bsOGDceVUvv7cOoc4Hgffj6WybUnrkS+/kS+djh1/aP6eiBbkohS6gbgYWAQ8KpSapPW+grgYuB7Sql2IAD8m9a6+lzH01oP6mM867XWRed+ZPyRa0/Ma4fEvv5EvnaI7PXbNTvrReDFTu5/Hng++hEJIYQIhxNnZwkhhIgRkkSMx+wOwEZy7Ykrka8/ka8dInj9SmuZHSuEECI80hIRQggRNkkiQgghwpbQSUQpdaVSardSqlgp9YDd8VhBKVWmlNoarIq8PnhfllLqLaXU3uDtwOD9Sin1q+DvY4tSaqa90feeUuoPSqkKpdS2Dvf1+nqVUrcHH79XKXW7HdfSW11c+4NKqUMdKmNf3eF73wxe+26l1BUd7o+5/wulVJ5S6j2l1I5ghfB/D96fKM99V9dv/fOvtU7ID8ANlABjgCRgMzDJ7rgsuM4yIOeM+34MPBD8/AHgR8HPrwZeAxQwB/jA7vjDuN6LgZnAtnCvF8jCLHLNAgYGPx9o97WFee0PAl/r5LGTgn/zycDo4P+CO1b/L4ChwMzg5xnAnuA1Jspz39X1W/78J3JLZBZQrLXep7VuA54BrrM5pmi5Dvhz8PM/A9d3uP8v2ngfGKCUGmpHgOHSWi8Hzlyg2tvrvQJ4S2tdrbU+AbwFXGl99H3TxbV35TrgGa11q9a6FCjG/E/E5P+F1vqI1vrD4Of1wE5gOInz3Hd1/V2J2POfyElkOFDe4euDdP9Lj1UaeFMptUEpdVfwviFa6yPBz48CQ4Kfx+vvpLfXG2+/hy8Fu2z+EOrOIY6vXSmVD8wAPiABn/szrh8sfv4TOYkkivla65nAVcAXlVIXd/ymNm3bhJnnnWjXCzwKFADTMVWyf2ZvONZSSqVjql7cp8+o/p0Iz30n12/585/ISeQQkNfh6xHB++KK1vpQ8LYCU2pmFnAs1E0VvK0IPjxefye9vd64+T1orY9prf1a6wDwO8zzD3F47UopL+YF9K9a6xeCdyfMc9/Z9Ufj+U/kJLIOGKuUGq2USgJuAV6yOaaIUkr1U0plhD4HFmO2Jn4JCM06uR0I7e/yEvDp4MyVOUBth66AWNbb630DWKyUGhhs/i8O3hdzzhjTugHz/IO59luUUslKqdGYitlridH/C6WUAh4HdmqtH+rwrYR47ru6/qg8/3bPKrDzAzNDYw9mNsK37I7Hgusbg5ldsRnYHrpGIBt4B9gLvA1kBe9XwK+Dv4+tQJHd1xDGNT+Naba3Y/pz7wzneoE7MIONxcBn7b6uPlz7E8Fr2xJ8MRja4fHfCl77buCqDvfH3P8FMB/TVbUF2BT8uDqBnvuurt/y51/KngghhAhbIndnCSGE6CNJIkIIIcImSUQIIUTYJIkIIYQImyQRIYQQYZMkIoRFlFIDlFL32B2HEFaSJCKEdQYAkkREXJMkIoR1/hcoCO7j8BO7gxHCCrLYUAiLBKupvqK1nmxzKEJYRloiQgghwiZJRAghRNgkiQhhnXrMVqVCxC1JIkJYRGtdBaxSSm2TgXURr2RgXQghRNikJSKEECJskkSEEEKETZKIEEKIsEkSEUIIETZJIkIIIcImSUQIIUTYJIkIIYQI2/8H/y2IrorCRQAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "trajectory = np.stack(trajectory)\n",
    "plt.plot(trajectory[:,0])\n",
    "plt.plot(trajectory[:,1])\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('theta, x')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}