Figure 1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import ode\n",
    "from matplotlib.colors import LinearSegmentedColormap\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from numpy import sin\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_cmap(color1, color2, color3):\n",
    "    colors = [color1, color2, color3]\n",
    "    cmap_name = 'list'\n",
    "    n_bins = 100\n",
    "    cmap = LinearSegmentedColormap.from_list(cmap_name, colors, N=n_bins)\n",
    "    return cmap\n",
    "\n",
    "def build_hankel(data, rows, cols):\n",
    "    X = np.empty((rows, cols))\n",
    "    for k in range(rows):\n",
    "        X[k, :] = data[k:cols + k]\n",
    "    return X \n",
    "\n",
    "def HAVOK(X, dt, r, norm, center=False, return_uv=False):\n",
    "    if (center): \n",
    "        m = X.shape[0]\n",
    "        X̄ = X - X[m//2,:]\n",
    "        U, Σ, Vh = np.linalg.svd(X̄,full_matrices=False)\n",
    "    else:\n",
    "        U, Σ, Vh = np.linalg.svd(X,full_matrices=False)\n",
    "    V = Vh.T\n",
    "    polys = true_polys(X.shape[0], dt, r, center)\n",
    "    for _i in range(r):\n",
    "        if (np.dot(U[:,_i], polys[:,_i]) < 0):\n",
    "            U[:,_i] *= -1\n",
    "            V[:,_i] *= -1\n",
    "    V1 = V[:-1,:r]\n",
    "    V2 = V[1:,:r]\n",
    "    A = (V2.T @ np.linalg.pinv(V1.T) - np.eye(r)) / (norm * dt)\n",
    "    if (return_uv): \n",
    "        return A, U, V\n",
    "    return A\n",
    "\n",
    "def true_polys(rows, dt, r, center): \n",
    "    m = rows // 2\n",
    "    Ut = np.linspace(-m*dt, m*dt, rows)\n",
    "    poly_stack = []\n",
    "    for j in range(r):\n",
    "        if (center):\n",
    "            poly_stack.append(Ut ** (j + 1))\n",
    "        else: \n",
    "            poly_stack.append(Ut ** j)\n",
    "    poly_stack = np.vstack(poly_stack).T\n",
    "    Q = np.empty((rows, r)) # Perform Gram-Schmidt\n",
    "    for j in range(r): \n",
    "        v = poly_stack[:, j]\n",
    "        for k in range(j - 1): \n",
    "            r_jk = Q[:, k].T @ poly_stack[:, j]\n",
    "            v -= (r_jk * Q[:, k])\n",
    "        r_jj = np.linalg.norm(v)\n",
    "        Q[:, j] = v / r_jj\n",
    "    return Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAADnCAYAAADl9EEgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAD5ElEQVR4nO3doY/XdRzH8c8dbESZkQ3UmdVR3NSAMxEMBnU2m5uJWQgkRyBYHImNZmNIMBBIToO6WZyYnSKbkUmhwfEPsMPgm3t+z8cj3m2vz68899suvG9nb29vAT27B/0BgCcTJ0SJE6LECVHihKij+/720cPZP+Xevzs6v547ObsP/4XdIztP/PGz/hzAvyNOiBInRIkTosQJUeKEKHFClDghSpwQJU6IEidEiROixAlR4oQocUKUOCFKnBAlTogSJ0SJE6LECVHihChxQtTOfv9l7K97D0bv1p46fmxyfq07P87ur7XWC2/Ov8Hh5m4tbIs4IUqcECVOiBInRIkTosQJUeKEKHFClDghSpwQJU6IEidEiROixAlR4oQocUKUOCFKnBAlTogSJ0SJE6LECVH73q1djx6O3q396Y9/JufXGy8dH91fa611+9rs/qsfze5z8NythW0RJ0SJE6LECVHihChxQpQ4IUqcECVOiBInRIkTosQJUeKEKHFClDghSpwQJU6IEidEiROixAlR4oQocUKUOCFKnBB1oEelp9345e/xN94/fWL2ge8uze6/fWF2n6dzVBq2RZwQJU6IEidEiROixAlR4oQocUKUOCFKnBAlTogSJ0SJE6LECVHihChxQpQ4IUqcECVOiBInRIkTosQJUeKEqEN9t/ZZuPzt76P75955eXR/ffPp7P57V2b3DwN3a2FbxAlR4oQocUKUOCFKnBAlTogSJ0SJE6LECVHihChxQpQ4IUqcECVOiBInRIkTosQJUeKEKHFClDghSpwQJU6Icrc27rOvfxvd//KDV0b311fvzu6vtdbHN+ffmORuLWyLOCFKnBAlTogSJ0SJE6LECVHihChxQpQ4IUqcECVOiBInRIkTosQJUeKEKHFClDghSpwQJU6IEidEiROixAlR7tb+z3149efR/eufvD66v9Za6/Jrs/vnfp3dd7cWtkWcECVOiBInRIkTosQJUeKEKHFClDghSpwQJU6IEidEiROixAlR4oQocUKUOCFKnBAlTogSJ0SJE6LECVHihChxQpSj0ox664vvx9/44fyZ2QcuPj+7//l9R6VhS8QJUeKEKHFClDghSpwQJU6IEidEiROixAlR4oQocUKUOCFKnBAlTogSJ0SJE6LECVHihChxQpQ4IUqcECVOiHK3ls178cKt0f0/L50d3V+7R9ythS0RJ0SJE6LECVHihChxQpQ4IUqcECVOiBInRIkTosQJUeKEKHFClDghSpwQJU6IEidEiROixAlR4oQocUKUOCFq/7u1wIHxzQlR4oQocUKUOCFKnBAlToh6DIeUYKnv/8pXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Linear dynamics\n",
    "np.random.seed(4)\n",
    "plt.figure()\n",
    "vals = np.array([1, 2, 3, 4, 5, 6, 7, 8])\n",
    "plt.imshow(np.diag(vals,1) + np.diag(-vals,-1), cmap=create_cmap('tab:blue','white','tab:orange'))\n",
    "plt.xticks([])\n",
    "plt.yticks([])\n",
    "plt.axis('off')\n",
    "plt.savefig('pdf_files/fig1a.pdf', bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply HAVOK to time series x(t) = sin(t) + sin(2t)\n",
    "ϵ = 10\n",
    "t = np.linspace(0, ϵ, 10001)\n",
    "dt = t[1] - t[0]\n",
    "rows = 41\n",
    "cols = t.size - rows + 1\n",
    "X = build_hankel(sin(t) + sin(2 * t), rows, cols)  \n",
    "A, U, V = HAVOK(X, dt, 4, 1, center=True, return_uv=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAADnCAYAAAD7GCa6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2dd7gcZfXHP3tvekiAEJCONCk6gEtRBHEZfoB0EUFRYBFQB0eF2JDFBuoqFiDq4NDElSKC9CIgDAsWqisyKgJSQ01ISO+5+/vjzCWXkDszuzuz9XyeZ5+Q7Dvve0jmfued856SqVarKIqiKN1JX6sNUBRFUdJDRV5RFKWLUZFXFEXpYlTkFUVRuhgVeUVRlC5GRV5RFKWLUZFXFEXpYlTkFUVRuhgVeUVRlC5GRV5RFKWLGdFqAxRFUXoNx/IOAbYG7gD+Y7tmavVlMlq7RlEUpbk4lnczcHDw2xcRsb8DuMt2zTeSXEtFXlEUpYk4ljcamAmMX83XA8BDwAm2az6exHrqk1cURWkue7B6gQfR5N2AV5NaTEVeURSluewf8f2DSbpsVOQVRVGaS5TI35HkYqmKvFEyckbJ2DLNNRRFUToFx/LWB3aMGHZ7kmumJvJGyTgCeSLdaZSM9dNaR1EUpYPYL+L7WcAjSS6YisgbJeMzwNXAKGAL4HajZKyVxlqKoigdxIfDvhxg4G7bNVckuWCiIm+UjIxRMgrAhavMvSNwk1Eyxia5nqIoSqfgWF4fsG/YmIc2vWUbo2SMS3LdxETeKBl9wDnAD4YZ8kHgKqNkaJatoii9SBaYHDbgqcmP7AD8ySgZaye1aJI7+QuBUyPGHApcZJSMTILrKoqitD0DDIS6amaOfZkFo+cAfAD4Y7BxbpgkRf52IE767PHA2QmuqyiK0tYYJSMzZ+z0z4eNmbbWfwf/czlwlp/3B5JYOzGR9/P+HwAr5vCvGSXja0mtrSiK0s6MX7LmjycuXneDsDEvrhT54/28f1tSayd68Orn/QuBb8Yc/mOjZHw6yfUVRVHaDaNknLLu/M2+2l/tH3bMsr6lvDLxaYBT/bx/RZLrpxFCWQR+HnPsxUbJOCwFGxRFUVqOUTKOAc7bZM62oeNemfg/VvQt/4Gf96cmbUPiIu/n/SowBYjzNOoDfm+UjA8lbYeiKEorMUrGgcClVGGT2eEiP3vM9L8B30rDjlSSoYIDg08Df4wxfDQSQ79TGrYoiqI0G6Nk7AH8ARix5uJ1mbhkndDx28x432eCDXLipFbWwM/7y4AjgftjDJ8I3GGUjG3SskdRFKUZGCXjvcAtwFiI3sVXqb4wesXYRGrHr45UC5T5eX8B0v3k3zGGrwfcbZSMLdK0SVEUJS2MkvEe4E/Am2VcNpm9Xeg1GTJ3pNn+L/VSw37en4WU1nw+xvCNAM8oGZuma5WiKEqyGCVjW+Bu4E3fTN9APxvO3Srq0kRLC69KU+rJ+3n/JaT62owYwzdDhH7DdK1SFEVJhqCk+t2IR+JNNpi3JSMHRodduiK4LjWa1jTEz/tPAgcA82IM3xJx3awXOVJRFKWFGCVjM8AD3rYxjfLHAw/Yrjk7DbsGaWpnKD/v/x04DFgaY/i2wF1GyQg/llYURWkRRsnYCNmJr9bFvHG0yKfqqoEWtP/z8/49wBFIfYYoDCTqRmvRK4rSVhgl4x2IwK+2+924pROZvHCjqGm6T+QB/Lx/C/AJxB8Vxc5IRbYJ6VqlKIoSD6NkTAbuAoYN+47hqpkF/H3wN5Vibkwixq1Cyxp5+3n/WuA44lWufD9wq1EyxqdrlaIoSjhGyZiE7MDfEzYuhqvmzt03PataKeYOqBRzdwK3JmTiW2iZyAP4ef9K4MSYwz+I7ugVRWkhwQ7+bqQByLBkqhk2nhOe27neGn9fguQQ3YZ0jDIrxVzimf8tFXkAP+9fCoTWWR7CBxEf/ZopmqQoivI2gmi/e4BIIZ68YJPq2OVrhI7ZZOK9eSTAZChRjZdqpuUiD+Dn/V8BX445fHfgTj2MVRSlWRglYwOgTISLZpBdpx14Q9j340a+xqgR81f31dGVYm79mg0MoS1EHsDP++cChZjDd0Pi6CelaJKiKMpgmGQZCK9PsJKTN529XWiOz5pjnh7uq1HE92zEom1EHsDP+z8Ezoo5PItkxq6bokmKovQwQYmVe4F3xbxkynf9T94E1d3DBq01dliRBzi5UsyNjbleJG0l8gHfBX4Uc+yOwD1BvKqiKEpiGCVjc0TgVxsHvyrrLO/7QemldTaaMHraU5AZVlv7MkuZOPqFsKkmA5+qydgQ2k7kg5rKBeB7MS95N1AOfGaKoigNY5SMrRCBf2ec8QfNG/OPc15b+xvAV+cs2Xxc2NiJo5+nLxOaInQL8K+YpkaSqVZTq3DZMEbJ+Bbx3TdPAfv4eX9aiiYpitLlDKkmGatI4nGzx7PPAsljqlbh7y9NYdnA8JHe71z7djaY8NCqf7wQ+A0wNVsoP1mH2cPS1iIPYJSM04jvvpkG7Ovn/SdSNElRlC7FKBm7Ih3tImtmZapw/Ozx5BauTFSdt2Rj/vXaCaHX7bSBw9iRMwd/+xLwC+CibKE8q167w2g7d82q+Hn/bOArMYdvAvzFKBk7p2iSoihdiFEy9kGqScYS+JNWEXiAmQvDA3DGjJjJmBEzAR4GPglsni2Uz05L4KEDRB7Az/vnAF+MOXwy4qPfO0WTFEXpIoyScQSSeRqewYQI/OfeWIM9VxH4ahVmLQrPcp0w+sUnMhn2BN6XLZR/ly2UlzVgdiza3l0zFKNkfA5wYw5fCnzCz/vXp2iSoigdjlEyPoPoSuSmt68KJ89ag90Wv70RyIKl7+CxVz8XNcX7bdd8sC5D66QjdvKD+Hn/AqTWTZwn0yjgD0bJCHeQKYrSkxglIxOc+V1IDC3sr4I9jMADzFoYWZDsJcRN01Q6SuQB/Lz/a6R6ZZwyxX3AJUbJ+Fq6VimK0kkYJSMD/ISYQR0jq/ClmRPYZRiBB5i1KFLkb7BdcyC2kQnRcSIP4Of9y4HDgcUxL/mxUTLODv5hFUXpYYySMYIqvyZmQMe4gQxfe30iOy0ZNeyYRcsmLV+4LDIn87r4ViZHR4o8gJ/3bwb2B+bGvOTrwK+NkjH8v5SiKF2Nff4u66y3vO8RMhwfZ/yaKzKcPmMi2ywdOdyQN4Cz/zP92LMjppoF3Bff0uToWJEH8PP+fUAOmB7zkuORmvRawVJReohKMbfeDT/Z86fPjlr+6vQRAzvGuWbd5X2cMWNNNl0+YnVfPwXYwCbZQvkbS1esuU/EdDfZrhmn5WnidFR0zXAYJWNr4E5ipiAD/wEO9PP+86kZpShKy6kUc+8BprwwYvkx564zb9SsEfFc4hsv6+drr09krYG37YPvAc4BbssWygMAjuVtBLwYMeWhtmveXJv1ybDaR1Sn4ef9p4ySsQci9O+Occn2wINGyTjYz/uPpGudoijNpFLM9SGu3CnAvv8avZRfTJrP4r54G9qtlozgyzMnML76psAvA64EzssWyo+u5pLDIqZcAPwp1uIp0BU7+UGC+vK3Ij1h47AQONrP+zelZ5WiKM2gUsyNA45FuittC3DfuMX8Zq0FrIgZcrHD4pF8YdYERlczAK8DvwLOzxbKrw53jWN5dwFh7pprbNc8Kp4FydNVIg8QNPu+FnmSx6EKnOrn/Z+nZ5WiKGlRKeY2QPzjFkFJgipVrp2wiJsnLoo9z+4LR3HSG2swgsx/gPOAy7OFcugEjuWtA7wG9IcMO9p2zatiG5IwXeGuGYqf9xcYJeNQ4CIknj6KDDDVKBlbAF/x836c+Pu2wLG8DDAeKeUwGnjVds05rbVKUZpDpZh7L+KS+QTwZvjLMqpcsvZ87h+3NPZcB84bw5Fzx93RR+Zc4M5soRx393sw4QK/FCmX0DK6bic/SBAT/22kCUlcbgaO8fN+3LDMVHEsrw94H7AXUnxt8mo+q2ZnzAD+t5rPU7ZrvtEcyxUlHQJ/+8GIuOdW/X523wC/nDSPp0bHC2TJVOHA+WP+ctTc8Va2UP53rfY4lncjcGjIkNts1zyo1nmTpGtFfhCjZOSRXf2wga6r8DhwmJ/3n0rPquFxLG8kIuofRRK+kmyG8jLSkOBq4N5WhXQpSq1Uirk1kBDoU4CtVjfm2ZHLmbrOPN7ojxdB019l6frL+/O3n/RoXa4Ux/LWQPz2w6fBwkm2a15Sz/xJ0fUiD2CUDBPJNlsz5iWzkeJmd6Rn1UocyxsD7IsI+6FAMxqUz0DOLq5BBL9j3FRK71Ap5jZBKtB+Bhg2v+WvY5dw6drzWRY3p73Ky2Q42M/7/6jXNsfyPob8/AzHALC+7Zoz6l0jCXpC5AGMkrE94hvbLOYlA8A3gJ8GLQkTx7G87YM1DidGidMUmc5Kwb9PBV9pNZVibjfEJXMkIT7vFVS5euJCbp8Qt8IJAI8BB/l5Pyq2PRTH8q4Ejg4Zcq/tmrlG1kiCnhF5AKNkrA/cBOxaw2VXAif5eT/+MX0EjuWtBXwH2aGEHdq0gueAIvAb2zVTr3WtKINUirkRwEcQcf9A1Pj5mQHOnzSff4+p6Ta9Aziq0XM3x/JGI2/Dw/f5g1Nt15zayDpJ0FMiD2CUjHGIcEclMAylAhzu5/3QFutRBAepnwZ+CKzbyFxN4DmkmfplKvZKmlSKuYlICfEvETNr/cUR4n+fHjODNeBi4PN+3m/4fnYs7wCio2Y2s12zIc1Igp4TeQCjZPQDtbQVBHlqH+Hn/T/Xs6Zjee9HejnuUs/1LeQZROwv14NaJUkqxdzmiLCfSPiO+C08MmYJF609n8XxK28l7np1LO8i4KSQIY/YrlmLxyA1elLkBzFKxrFI5E3Y6fhQViA3y8/i3iyO5a2P1KzO12VkNMuRE/4VwEYprQEShnkWcKX67JV6qRRzGcQVMwU5i4ot1cupck3t/vfEgygcy+sHXiH8bfwM2zWLSa3ZCD0t8vBmd/brqU0gbwaO9/P+sM13g0QlG/Fvx96lDMPTiI0+IuiDn5nAXNs1q8GaY4HNkRCzoZ9tgE0btGGQJ4Cv2K55a0LzKT1ApZgbCXwMEfead7gz+1fgTJrP06NqeplMJRzasby9gHsjhm1nu+Z/k1y3Xnpe5OHNA9lriXHYM4QXkAOct/VrdCxvFOL/O7YBs3wk7PM6wB8U8npxLO/dwFHBJ7KFTQyuBE6xXfP1BOZSupRKMbc2Ev74RWDjeuZ4dPRSLpw0nwUxC4wF3AQcm0Zio2N55yHx+sPxuO2a2ye9br2oyAcYJWM08EvC/WyrsgxpRjJ10H3jWN6ayAMjqr706ngMuAK43nbNVJKxgjeM9yChaR8H3tXAdK8jP7y/b/QhpHQXlWJua0QIj0dKb9TMCqpcO3Eht9bmngE5Q/qun/cTb7UXBE88h2SgD0fRds0zkl67XlTkhxCUQjgZmEptdX2uB06w7p+6BnLibtS49GvAaUgkS9N6QAaCbyARPxYwps6pbgI+b7vmS0nZpnQegb/9Q4hL5hCkLlRdzOpbwfmT5scuTxCwAMj7ef/aeteNwrG8vQEvYtgutmv+PS0bakVFfjUYJeNDwB+Q2jCxmDx/k5cO/9epI/urI9arYanlyAPle60uLBYcEJ+GPOTiHkQPZQ7wVeAS3dX3FpVibhTyVjgFeG+j8/1z9NIFzqR5fUv6GFvDZc8i/ne/0fXDcCzvN4QHUbwAvLOdfgZU5IfBKBmbAr8nRm36jea8i/2eOIHRK2q5J7kT8Wm3xeHMII7lbYhEEH2W+sTeAz5ju+YziRqmtB2VYm4d5A3QJoEaS0sy1ad/ts7caU+MXp6r8dLIQIgkCGrVvEq4++kntmt+PU07aqWje7ymSZD49CGk1dewvGvGrhz4uFWLwD+LZPV9uN0EHsB2zZdt1/wSsCXgIKVSa8EEHnMsr2VNEpR0qRRz21aKOReYBnyfxgX+3j+NXzTlsxvMWlKjwK8Avobs4FMV+ICPEn2+UGqCHTWhO/kYGCXjI8BvGFrgrArZl/Zjt2k1VREtAZbtmjWfJLUKx/I2QRLHwmp0DMdU4Ou2a9b6oFDajMDfvg/ikjkwgSmXA79bmBn4+ckbvrEPclgat1IswEvAx/28/9cEbImFY3l3I5uY4fi77Zptl+yoIh8To2RsjpTo3YUq7P78R9jxlb1rmeJM4Mx28tXVgmN5hyCt0GpNuLofOMp2zYaKQSmtoVLMjUEe8FOoPaBgdcwCXMDJbzRzFPBb4IM1znE7Eh7ZtPBdx/I2RaJqwg6Tv2S75i+aY1F81F0TEz/vPwvsCfz83a9+MLbAD7CC+za/era7+ykPdKrAAwSd5t+NZAjXwu5AxbG8ekJKlRZRKebWqxRz3wGeB35N4wL/BOK/3yS/0cxv5jea+X9IyHAtAj8AFJAKks3OzziWcIFfBvyuSbbUhO7ka8SxvP2qVP+YIRP5gFzat5g7t7mUF9d60/XuAKf5eX9BqkamjGN5JiL2W9Rw2QDwLeBHzQwTVWqjUsy9B2mEfQz1Hbyvyt3Iudbt2UJ5wCgZ6yE/Bx+rcZ5XgKP9vB+VaZo4QajxE8DWIcNusF3z8CaZVBMq8jXgWN52wAPAxKixC0bO4bbtLmDm+LeFjj8HfM7P+3cmb2HzcCxvPOJHPZXa4qFvAY7TVoTtQ+Bv3x9xyeyXwJRLkaS+87KF8mPwZg7KscC51N4U5zbgBD/vv5aAbTUTFBe8P2LYR23XvL4Z9tSKinxMgq7sDxFj9zpr7Cvctt0FzB8dqmO/Bb7s5/2ZCZnYEoIfgMuRaJy4PAt8xHbNx9KxSolDpZgbiwjvqcB2CUw5Azm3OT9bKL8pyEbJeCdwAbU/QBYBXwYuSKtxTxwcy/sV4moajpnAhu0aYKAiH4OgFs2dSEhlKC9PfIrbt7mEpSNi9RiZgaR+X9XKm7hRgiYolyKhoXGZB3zMds2OfqPpRCrF3PpIbLtFDQl/IfwHcclckS2U34wcC0p6fxH4ATCuxjkfBo7x8/6TCdhXN0FrzlcIaT0I/NJ2zS82yaSa0YPXCAJ/3PnEEPg5o1+fcdu2F8QVeJBSpVcCNwfJVx2J7ZqzkRjiryOxy3GYANzmWN4JqRmmvIVKMbdjpZj7DZKV+U0aF/g7gA8D78kWypesIvAG8DfEPVOLwK9ASlrv0WqBDziEcIGHNoyNH0ot9Vl6lSlIU4MoXlpzyeTdlvcv2xmpQFlLeYODgH8bJeN04Fd+3u+4eu1B5NBPHMt7CMkUfkeMy/qBSxzL2wL4VidHH7UrlWKuD4lrn0J4jHdcFgOXIf72/6z6pVEyxiARMKdTu748jezeH2jYyuSI6gPxH6Bt6tSsDnXXhOBY3sFI8a2og8VFwJ62a1YAggiCi4BD61j278CX/Lz/tzqubQscy9sAEfpawuOuBE6wXXNJOlb1FpVibjxwHOJvb6TS6CCvIVExbrZQnrHql8HB6mGI22bzOua/CDmjmt+QlQniWN47kKSrsD7Mp9mu+eMmmVQXKvLD4Fje4OvmGjGGf8x2zbdUvgtu+uORmz7qdW91XIaEW75Sx7Utx7G8kUjDlK/WcNm9wOEaeVM/lWJuQ+ALiL997QSmfAy5h6/KFsqrfQAbJWM74Dzqi8yZBlh+3o/ql9p0HMv7MvCzkCEDwCa2a77cJJPqQkV+NTiWNwlp3r1ZjOHftF3zB8N9aZSMDZDerkfUYcp8JExxqp/3O3KH61jeRxGfZZyHJcB/gQNt13w2Pau6j0oxl0VcMp8gGTfsLYg//Z5sobxakTBKxprAt5E+rbWuWUXeDAp+3p/XiKFp4VjeP4EdQobcabvm/s2yp15U5FeDY3mXIckgUVwBHBvHl2yUjMORm7qeYk5PAaf4ef+PdVzbchzL2wm4Fdgw5iXTgYNt13w4Pas6n0ox1w8cjIQZ7pXAlIuQGk1Ts4XyE8MNMkpGH+IK+hHxzl5W5XHgpHZ2STqWtyPwaMSwT9mueWUz7GkEFflVcCzvQESQongA2LuWYmNGyVgL+cH4XJ3m3QJ8zc/7bVe9MgrH8jZG/l7DdkZDWQAcartmVIOGnqNSzK2BuAJPpbb8hOF4GXnbvDBbKIdWczRKxp7AT4H31bHOMuCHQLHd30wdyzsHeTMajnnA+rZrLmySSXWjIj8Ex/ImAv8muhflNGBX2zXrysALmpJcRHia9HAMIDHpZ/p5f1o967eK4O/3D8C+MS9ZghQ3uyk9qzqHSjG3CeJv/yz1nfOsyt8Rl8w12UI5NJEnCIksIm8O9fAgsnv/V53XN43gPOlFwiPkLrFds5ZWoS1D4+Tfyo+JFvgFiCuh7hTroP7GjsiuZlmNl/chIZ1PGSXjZ0bJSCKZpSnYrjkXCRe9JOYlo4HrHMv7ZHpWtT+VYm7XSjF3JZIp/HUaE/gqcAPi3tk1WyhfESbwRsnY3CgZlwH/pD6Bn4e8cezRCQIfsD/RIdBtHRs/FN3JBziWlwPuiTH0SNs1/5DUukbJ2AaJTPhwnVPMA34CnNtO4WdhBAlmBaThRByqSA9ZNz2r2ovA334Y4m/fI4EpFyAP159nC+WnowYbJeMdwBlIlE4tdd6HcilysPpqnde3BMfyriG8gNozwFadktehIg84ljcOCRWL8m9ebbvmx5NePwi3PAgR+3p9rNOR9PEL2t3fOYhjeZ9CytiOinnJN2zXPDtFk1pOpZibCJyARKzUE2++KtOAnwMXZwvl2VGDg4iZryAPl6guSMPxMPBFP+8/WOf1LcOxvMlIbHzYPfld2zXPbJJJDaMiDziW9zPkpg5jJrC97ZrT07LDKBmjkcOeb1L/D9jLyMHYRZ2ws3cs70OI+yCuC+KHwBmdsouKS6WY2wwR9pOIUeU0Bg8i/vZrs4Xy8qjBQQLfqUhNm3rXn470By75eb8jy0k7lldANkthbNlJPYx7XuQdy3sfkvQUdT5xjO2aVzTBJIySsRFyPtCIL3oW0n7vl03qf1k3QbjancQvBeEgXXg6UkiGUinmdkce7EfQ+BnZAHAdcE62UI4qjQuAUTI2QxLWTgLG1LnuciQ650w/78+pc46WExQifJbwUN8/266ZRLhq0+hpkXcsbzSS9LR9xNBbgUOavXsMwtXOA3ZuYJr5SJnXc/y837aZeY7lvQu4C9gk5iWXIWUQInep7UalmBuBFHSbArw/gSnnIfWSfpEtlGMlkRklY1vgNCQfpJHkqduQsN631bHpNIID/qiN3Am2a17aDHuSotdF/kwkYy+MucC7W9WjNPDXH4G8QjZSg2QpkujyUz/vP5WAaYkT9NG8i/ihpVchyWgdIfSVYm5NZMf8JSCJqqPPI29rl2QL5blRg4N7aTfga8hDppZmL6vyN+B0P+/f18AcbUMQDPAgsGvIsNeRMgaxc2PagZ4VecfydkDihKN2MZ+zXfPCJpgUilEyRiAV8b5LdJhnFHcgLo/b2q3iZVAU6k7iJ01djWQetq3QV4q5LRBhP5H45R3CGCzhe0NMf/tY4ONIjH0jb4UA/0Iio27p5B4Iq+JY3u7I32sY37NdM2pT2Hb0pMg7ljcCyViNuuHvAfZpp0O+4Af288gPWq1t1FblOaSTzyXt1KHKsby1ETdAXFfGNYjQ15pzkBpBS709EJfMR2jc374CSSQ7N1sox4paCToynYy8PTR6rzyP9Oi9st02BkngWN7vgaNChiwDNrNds+MKBvaqyH8FiUAJYxFg2K4ZGVPcChIKdRtkCeL6+KWf9x9p1LYkcCxvDeBG4tdAvxY4utVCXynmRiIx1lMIf/WPyxwkO/oX2UL5hajBQV2Z/0N27QfTmEsGpHvZ9+mg0NxaCdyEzxBeUvgy2zWPa5JJidJzIh9UmHwGWDNi6Jdt1zy3CSY1hFEy1kVaCH6B6P+nOFSQQ83ftapx8iBB67XfE78u//XAJ1rRa7NSzK0NfAZpd9eoOw3kHp0KXJotlCOrNBolY2ukX+sxJBNf/xKSZHexn/cXJDBf2+JY3tlIJnEYu9iu2dbNQYajF0X+R0hUQRgPIE1AOua1NNjZW8jOvpauVMOxAvGNXw7c4Of9lhRiCuqI/BYpoRuHG5F6N00R+koxtxXykP00jb9RAfwZqd9+c7ZQDr3/jJKxDuJrP5ZkonRAKp7+CLjcz/tt2Zg6SRzLG4/UqQnL0+i4sMmh9JTIBx2LngbGhgxbCrzXds2ODAkLfPafRnYmcerhx2E+4g65DCg32ycbnKH8BvhUzEtuRspPpOJeCPzteyEumUNp3CWyHDlAPjdbKIe6y4KEuYOQUr8HUn/JgVX5J1KA7Npu9LkPh2N5JyM9nMM4wnbN65phTxr0msifjxxEhfEd2zXPaoY9aWKUjJHI7vd0YLsEp34VEdGbgLv9vB+7a3kjOJbXj9RCOTbmJbciP5yJCX2lmBuFHM5NAbIJTDkbyWH4ZbZQHjZE1ygZayOCfhhS42hCAmsP8lckPPf2boqWiYNjeX1Ij9ZtQoY9h9Sp6dgHX8+IvGN5WyJdh8JCJl9B/kHbvkZ0XIKDuP0Qn/2BNL7rHMpCxKVzIxJS93qCc7+NQOgvRmqpx+EWpDVjQ0JfKeYmIT0AvkD8xidh/A9JcitlC+XVlp4IMlEPCz57kUy3p0EGD9odP+/3bGMWx/IOQKK4wviK7ZrnNMOetOglkb+c6Nf9k7u50qFRMrZA/PYn0nhI3aoMIHHGtwAe8A8/7yceux7svi5CinjFoW4ffaWYexdSzyUPjKv1+tVQRuLbb8kWym8pyRC42T4A7I1ExeyYwHqr8jwrQ2ZTfSB3Ao7l3UF4X9r5wMa2a3ZsqQboEZEPmnL/k/Bd7NPAdq0OwWsGQ5JjbGCXlJaZhxwilpF8g0eTEv1A6C9A4r/jcB0SdRP5bxv423PIAXa9DTKGsgzZNZ+bLZT/MfiHRskYgxyW7h183kf8apy1cieS/HZrL/nbw3As791IYlcYv7Bd80vNsCdNekXkbwIOiRjWEf0ak8YoGbsh5xQfI5lszOGYi4j+fUim8T8aKZwWCP2vkC5JcbgG+ORwmbGBv/0TiLgnsYueBbiAk99o5ivARogff2fE/bO5ETMAABc5SURBVLI70hQlLV5D6rBc4Of9J1NcpyNxLO8Cwu+dKrCN7ZptWQKkFrpe5B3L+wByuBSGD+zUDVUN68UoGeMQ/++xyCtsWGJIUjwH/AOJza8AlVoaTARC/0uiD9MHeVutm0oxN5mV/vb14649HFWqTzw+evllztrznpnfX30PIuxZkglrjWIRkitwGXBXGu6ybsCxvHWQsMmwqps3264ZNz+jrelqkQ+KDt0DfChi6KG2a97cBJM6gqAr0NFIYk2jtU5qZTrwJOI++9/QX1e38w/+jX+BuJ7icCVw3O6bnrU1K/3tNZXYrVJlbl+V10asYPqIFbzWP8Czo5a/9syo5XMXZKrrkUkkKS2+OXKPX4aEP0YmTvU6juWdjoSLhrFPtzSR73aR3w8pxhXG/cAe7VSfpp0wSsb2iNgfQWNVMJPgDaTe92tIuv10YEammplx5D9PO3LSog0OiDPJWmOeemnbda/aiMwAK4BlmSrLMit/nd83wNy+KnP7BpjXN8Dc/qr82ldlTv8Ar/evYHHruyM/jJw1XNFpDd1bSZBc9yziPhsOH9ixWzSha0U+2OE9TPRONGe75r1NMKnjCfrRHoYkAH2AZMMxG6Oa4UPPHMV20z8Qa/iT6z6It+XvINMx9/9S4G4kWujmdu4N0M44lncC0Y3kT7Rd89fNsKcZJBl7224cQbTA36ECHx8/7z+BdKz6cdAu7iBE9PcjPIs4fTJV7t3iajLVPradEZ3h/64Z72NFZgX3bnF1Owv9G0hS143AHeqKaYygFtJ3I4a9jrj0uoauFPkgDf77MYYW0ralW/Hz/nQkA/XSICQzN+SzC42X1q2dTJV7t7yKDH1sM2O3yOHbTf8AA5kqf9786nZ5J1mM5BqUET/7g37e7/qQ3iZiEd157Fed1hQkiq4UeaSuR1iqMsA1tmtWmmFMtxOUNvhj8MEoGROBPVkp+jvTJNGvZqqUt7ySvmofW78enQLw7tf2oMoK/rL5ta0Q+qXImdA9iLA/6Of9rhKYdsGxvAnAGRHDliBhuV1F14l80Iz3OxHDViANEJQU8PP+XCRd/DZ4s0LmHkh99cGQwiTK8a6WaqaKt9UVUM2w9czo4KD3vLYX1UyVv77zurSF/klWhos+AjzQrNo/ClOAyRFjzu/EpiBRdJ3II0k9Uf0zf2O75hPNMEYBP+/PYYjoAwQ+/feyUvTfC2yZ1JrVzADe1pfTRx9bznxv5Hjj1Q8xkBng/s1uSELoVwBPIElfg6L+aPDwU5qMY3mTga9GDJsP/LAJ5jSdrhL5IKJmSsSwpUDHV5nsdAKf/h0MCXE1SsZ4YAtgK0Twh/66KTW6fKqZAe7e6rdkqn1sMSs6iXXHV/ammhnggU1viiP0y5BQvP+t5vN8L9Ri7yC+QXTlzp/ZrjmjGcY0m64SeSRVPMoRe77tmpFt1JTmE3Qg8oPPm1SKucy8voGDZvUPnLYkU91T4terzO0Pfg1+P69/gIWZKsszVZZnYBlVlmcGuGvrEvs++Wk2f8OItGGnl/ehf2Dk9L++89pHybwZiz8d3vLfrwEvakZp++NY3sZINnMYM5FGLV1Jt4n8qRHfL0NCAJUOoFLMjUUSsU6dMNC3/YSB+s5uB6iyZMKN/G/R2KXzF28VWQTMeHWv9YxX93rQds1v17Wg0k58m+gaQUXbNbvWldb6vL2EcCxvMyQ2PoyruvFgpduoFHPvqBRzZwIvABcC2zcyXx8ZfywDn15r1EuTCSKAYvAtx/KiDvCVNsaxvK2JLkn9ItGdoTqabtrJ20Q/tKY2wxClPirFnIGcqXyKZMru3obUb787WyhXs4BjeR9FulrtG+P67wbnPGd2S4p7j3EW0YX2zuy2uPhV6YqyBo7lrQFMo4ub8XYrlWKuD9gfKfH7fwlMuRhp/H1etlB+fHUDHMsbi7Qw3CfmnN8Hvq1C3zk4lrcTUuE0jKeA7YcrP90tdIu75jjCBR6k3ZrSJlSKubGVYu6zwL+RHXejAv8q8E1gk2yh/LnhBB7Ads1FSP2dcsy5vwl8P9jVK53BD2KM+Va3Czx0wU4+qCn+OOEVEp+jw5vxdguVYm59xLV2MrBOAlP+E3HJXJUtlGvq5epY3njER//BmJecDZyuO/r2xrG8PZEGNWE8CuzcCz0kukHkD0SKOIXR8c14O51KMbcj4m8/mmT87bcg4n5PtlCu+yYOXH23Ixm5cfgJcJoKfXsSvG3dh5TVCONA2zXjHsJ3NN1w8BoVNjmf6NKiSgoE/vYDEHGP6/8OYxHwG2BqtlBOJGPZds35juUdgAh9nDrFXwP6Hcv7qgp9W3Io0QL/Z+Tfuyfo6J18LzXj7SQqxdw4pI3gFKILxcXhFaT704XZQnlmAvO9jaCA1W1EC8QgU4EpKvTtg2N5ayJnPGENQQD2tF0zqiVo19DpB6+nRHxfRcRBaQKVYm6DSjH3fSS+3aVxgf8H8rB4Z7ZQ/mFaAg9gu+Y85K0jypc7yCnAVD2MbSt+QrTA39pLAg8dvJMPig5No0ea8bYzQ/ztnwRGNjhdFQlvPAe4rxF/ez0EPvpbiO4LPMivgc/qoX5rcSzPRDpnRfFe2zUfTduedqKTffKfJboBs4ZNpkTgbz8QEXczgSkXIk1IpmYL5acSmK8uAh/9QciDZu8Yl5wATHAs7xjbNbUoWQsIoqQuijH0sl4TeOhQd01QM96OGPYY0oxBSZBKMTeuUsxZSNjqzTQu8C8DpyPx7V9opcAPYrvmAuBg4u0MAY4EbnAsb1x6VikhfA+pXhrG68BXmmBL29GR7hrH8j4JXBEx7ATbNS9thj29QKWY2xB5sFrApCSmRFwy12QL5bbcAQeZsTcSrwQCSOjeId1c7KrdcCzv/UjLxKizkaNt17yqCSa1HR0n8sFB10OElxSeAWza7TUpmkGlmNuJlfHtSfjbb0Li25vub6+HQOivR0ovxOER4ADbNV9PzyoFwLG80cjh/HYRQ28CPtKrkVCd6JN/P9E147uuGW8zCfztByHiHscvHcVC5IByarZQ/l8C8zUN2zUXOZb3EeAPyN9JFLsA9zqWt6/tmi+na13P802iBX4u8PleFXjoTJ/8iRHfL6MLm/E2g0oxN75SzJ2M+NtvonGBfwnpyrNxtlD+YqcJ/CDBhuGjwNUxL9ke+ItjeVF+YqVOHMvbEbm3oviK7ZovpW1PO9NR7prgYOtVwlt5/dZ2zXyTTOoKAn/7FxB/+9oJTPl3VvrblyUwX1vgWF4/Ev9/UsxLXkZcN4+lZ1Xv4VjeCOBBpDdwGHcD+/byLh46byd/GNG9Gn/eDEO6gUoxl60Uc5chBdxOpzGBrwI3IPHlu2YL5Su7SeABglj4zxK/VdyGyI5+v/Ss6km+SrTALwQ+0+sCD523k7+d8AOwf9quuVOz7OlEhvjbvwzkEphyAeJv/3mnumNqJTj8/xZwZsxLVgCfs11Tayg1iGN52yGHrVEt/U61XVObBNFBB6+O5W1IdChbqRm2dCKVYm48kEcKum2dwJQvIm9NF2cL5TcSmK9jCHaHZzmWN4d4CXf9wMWO5W2O1DDvnJ1VG+FY3trI22KUwN8P/DJ9izqDTnLXHEO4vSuAK5tkS8dQKeY2qhRzP0RKQDg0LvCPIOGUW2QL5Z/0msAPJdgpngjErUl+BnB5EPqn1IBjeSOBawjvGwGwFDhRy0yspCPcNcHr8b8Ib+h8q+2aBzfJpLanUsxlEZfMx2n8ja2KxIqfC/y1E+Lbm4ljeUciyXlx8wjuAw63XXNWelZ1D8HP//lIYEAUZ9iuWUzZpI6iU9w1OxMu8KCumkF/+8GIuMctsBXGfFb6259OYL6uxHbNaxzLm4vE0q8R45K9gL85lneg7ZrPpGtdVzAY+RXFo0glSmUIneKuOS7i+9lIHZWeJIhvt4EnkDT8RgV+GtIcY5NsoXyKCnw0tmvegbQRjJsAtQ1wf9CqThkGx/I+TLxzj8XA8bZrdlVEVxK0vbsmKEb2MuH9QF3bNU9ukkltQ6WY2xjZ5XyWZOLbH0bCA6/ttvDHZuFY3sZIO8odYl6yAnmgnqcHsm8liKR5AJgYY/gnbNf8fcomdSSdsJM/kOiGz79thiHtQqWY27lSzF0OPAucRmMCPwBci3REel+2UL5KBb5+bNd8EdnR3xnzkn7kwXq1Y3lxxKwnCPpF3EI8gf+uCvzwdIJPPip79Snkad/VVIq5flb62/dKYMr5wMWIv/3ZBOZTAmzXnOtY3sFIeY2oMhyDfAzYwbG8I2zXjGpp2dUEb+/XEl0+GOD3wFnpWtTZtLW7Jniav0x41MI3bdf8QZNMajqVYm4N4Hgkvn3LBKZ8gZXx7XMSmE8ZhiAq5HSglvtzIZI4dXk6VrU3wd/ZxUgzligeBj5ku+aidK3qbNpd5L9IeJmCKvBO2zVfaJJJTSPwt38R8bevlcCUDwE/A67LFsrLE5hPiUnQ/+BSYFQNl/0KaRS+JB2r2g/H8vqAHxOvucdLwK62a76SrlWdT7uL/MOElxX2bNfcp1n2NINKMbcLUuL3KBp3pw0g8e3nAPdrfHvrcCxvL+A6os+XhvIw0uyi66ObHMsbg4RBHxVj+EJgT9s1/5GuVd1B24q8Y3nbA/+OGHa87ZodHx8f+NsPQcRd/e1dimN5myLlit9Xw2WLkLrpU7s1i9OxvElI6G/ccNIjbNe8LkWTuop2jq6JOnBdgBzOdCyVYm6NSjH3BSS+/XoaF/hpSIW+jbOF8hQV+PYicCvuRW11VcYibra/Opb37lQMayGO5b0T+CvxBf4MFfjaaMudfFC3+wWkVOtwdGzd+BT87Q8iLhn1t3cIgZ/+IqCW5t/LgO8DP7Jdsy374taCY3k7IzkF74h5yeXAcZpPUBvtKvL7A7dHDNvHdk2vGfYkRQr+9uuAc7KF8v2N2qY0n2Bnfi2S/VoLPlKE6+HkrWoOjuUdiBQci/uQ+zOwn7b1rJ12FfkrgE+GDJmGRNXErf7XMlLwt88DLkH97V2BY3kTkH/PI2u8dAB5ezvTds35iRuWIo7lfQaJHuqPecktSEbrgvSs6l7aTuSDrL9XEV/kcBRt1zyjSSbVhca3K3EJYsNPQYpr1fqG93pwndPuIuhY3pZIJdNDarjMBb5ou6a6IeukHUX+RCQyJIxtbdd8ohn21EpK/vafAderv727cSxvN6TqZz0HrDOAs4Ff2a65MFHDGsSxvPFAAQkKqCVX4HTgbPXBN0Y7ivxdQFjs+wO2a+7eLHviUinmdkZKDqi/XamboKHI6UiDkXruo9cQsXdbnQkavKEcBfwU2LiGS5cBn7Zd84pUDOsx2krkHctbB7lJw3x1J9uu6TbJpFC0noySFo7lGciuPiwZMIxXEbH/bSuakziWtwPiXqy17PUcpKHKPclb1Zu0m8gfj6R/D8dSYH3bNVvacm6Iv/0UYKsEpnwBmApcov52ZRDH8kYgZzrfA8bUOc1y4C4kCeuGNH92gsJiewOfQHpA1JqHMw04sNcLtCVNu4n8jcChIUNus13zoGbZsypD6rd/DvW3K03CsbytkDe8RpvBLAP+hAj+jbZrzk7AtnHA/sBHkQPVNeuc6jFE4F9q1CblrbSNyDuWtwYSKRDW5Pgk2zUvaZJJb6Lx7UqrCYp3fQYoApMSmHIZEnv+OPC/IZ9nV1cULfCvjwcmB59tgcOBAwiPhIvD5YBtu+bcBudRVkM7ifyRyA5jOAYQV82MZtiTYnz71Gyh/FwC8yk9SBBi/CWkUmMSb5OrUkXcJk8jP3OTh3zCNmD18BgSHnlfwvMqQ2gnkb8SODpkyL22a+bStiOF+PbnWelv152KkgiO5a2J3KNTqN9F0ireQIquXajx7+nTFiIfhI3NACaEDDvFds2w2vINkYK//QEkI1H97UpqOJa3FiL0pxKvVV4rqQIXII1+ZrbamF6hXdr/7UO4wINUaUycFOLbrwXOVX+70gyCw9PvOJY3FbmPTwHWaK1Vq+UvwJe0BnzzaZed/EXASSFDHrFdc9ek1kshvn0eK+Pbn0tgPkWpi2BnfwSyadmH+PVh0uJfyGHxVZq52hpaLvJBWeFXgHVDhhVs1/xho2ulEN+u/nalbQl6JH8EEXyT5gl+BYkgu852zcebtKYyDO3grtmDcIEHuWHqRv3tSi9iu+bryBvmxYHgH44I/t4kK/hV4G+sFPbnEpxbaZB22Mmfh+ysh+Nx2zW3r2du9bcrytsJwjC3R95mV/2E9aBdhuSyDH5eRnztN2pD7falpTv5IMHi8IhhNe3i1d+uKOEESUcPBJ+34Fje2kjo8Ea8XdTnqV+982i1uyYLbBoxJpbIV4q58ayMb1d/u6LUQVDb5pHgo3QBrRb5j0Z8/zwQGnJVKeY2YqW/fe0EbFJ/u6IoXUOrRT7KVXP9cK+HlWIui7hkPk4y/vbrgZ+pv11RlG6iZSLvWN52wHYRw97iqqkUc32s9Lc3WpEPVvrbf6H12xVF6UZauZOP2sXPQMKy0vC3a/12RVF6gnYW+Rt23/Ss9SvFs5L0tz+E1G+/Tv3tiqL0Ai0RecfyNiWirdmWk27cEniO5Pzt5wJ/yxbKGgKmKErP0KqdfOguvj+zmMnjfbPBNbRfqqIoPU9bivzaY5+iLzNQ79zTEH/7xepvVxSl12m6yDuWty7wwbAxk8bVVdPoISS+/Vr1tyuKoghNF/kJo174wrylmw7bxT2TWcZaY56OO92gv/0c4H71tyuKoryVpoh8EN9+EDClv+/ovcPGrjXmafr7lkVNqf52RVGUGKQq8pVibhyQR+Lb37ViYCRzFm8ees06Y/8b9vULwM9Rf7uiKEosUhH5SjG3IWADFjBp8M/nLtmMauiSA6w19snVfaHx7YqiKHWQqMhXirmdkKbCRwMjV/1+9qItQ6+fOPp5RvYvHvyt+tsVRVEaJDGRrxRzJyJ+8mGZvThc5Nca+z8Qf/sliL/9maTsUxRF6UWS3MnfCiwFRq3uy8XL12Tx8smhE4wZMfsc4HvZQnl2gnYpiqL0LMOGMtZKtlB+FbhyuO/nRLhqoPrKk68f+VUVeEVRlORITOQDzh3uiyhXDWTu1NZiiqIoyZKoyGcL5ceAu1b984FqX2ToJHB7krYoiqIoye/kQaJh3sK8JZtOX1EdE3ZNFfhTCrYoiqL0NGmI/B3AYPGZR4BPPj79U6FRN8AjtmvOTMEWRVGUniZxkc8WygPAKUgRst2yhfLvqvTvF3HZHUnboSiKokCmWk33rNOxvMnAdCATMuyDtmv+JVVDFEVRepA03DWrsi/hAj8XeLAJdiiKovQczRD5/SO+v9t2zciyk4qiKErtpCryjuVliBZ5DZ1UFEVJibR38jsA60eM0UNXRVGUlEhb5KN28U/Yrvl8yjYoiqL0LK0Wed3FK4qipEhqIu9Y3nhgz4hhKvKKoigpkuZOPscwZYcDlgL3pri+oihKz5OmyEe5au6zXXNBiusriqL0PGmK/IcjvldXjaIoSsqkIvKO5W0ObB0xTEVeURQlZdLayUe5al4G/pXS2oqiKEpAq0Reu0ApiqI0gcRF3rG8kcA+EcPUVaMoitIE0tjJvx+YEPK9doFSFEVpEmmIfFRUzcPaBUpRFKU5pCHyWspAURSlTUhU5B3LWxfIRgxTkVcURWkSSe/ko7pAzUG7QCmKojSNpEU+Theo5QmvqSiKogxDYiLvWF4f6o9XFEVpK5LcyW+IhEeGoSKvKIrSRDLVanKJp0FP1x2QHf3+SD35wXLD/7Vdc7vEFlMURVEiSVTkVyVoHJJDBP9F2zV/nNpiiqIoyttIVeQVRVGU1pJ2j1dFURSlhajIK4qidDEq8oqiKF2MiryiKEoXoyKvKIrSxajIK4qidDEq8oqiKF2MiryiKEoXoyKvKIrSxajIK4qidDEq8oqiKF3M/wPIYC4NuBm+WQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 468x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the Orthogonal Polynomials and the Frenet Serret Frame\n",
    "plt.figure(figsize=(6.5,4))\n",
    "\n",
    "# Orthogonal Polynomials \n",
    "#plt.subplot(2, 1, 1)\n",
    "#plt.title('Orthogonal Polynomials', fontsize=16)\n",
    "colors = ['peru','tab:green','tab:purple','darkcyan']\n",
    "labels = [r'$p_1$',r'$p_2$',r'$p_3$',r'$p_4$']\n",
    "for i in range(3):\n",
    "    plt.plot(U[:,i], colors[i], label=labels[i], linewidth=8)\n",
    "#plt.legend(bbox_to_anchor=(1.02, 1.0), loc='upper left', fontsize=14)\n",
    "plt.xticks([])\n",
    "plt.yticks([])\n",
    "plt.axis('off')\n",
    "plt.savefig('pdf_files/fig1b.pdf',bbox_inches = \"tight\")"
   ]
  },
  {
   "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.7.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}