figure1_nonnormal.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>.container { width:100% !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "import matplotlib.pyplot as pl\n",
    "pl.rc(\"figure\", facecolor=\"gray\")\n",
    "pl.rc('text', usetex=True)\n",
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "sns.set()\n",
    "sns.set_style(\"whitegrid\",{'axes.grid' : True,'axes.edgecolor' : 'k'})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "build lower triangular matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0., 0.],\n",
       "       [2., 1., 0., 0.],\n",
       "       [3., 2., 1., 0.],\n",
       "       [3., 3., 2., 1.]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def make_theta(diag_value,delay_value,triag_value,N):\n",
    "    return np.diag([diag_value]*N,k=0)+np.diag([delay_value]*(N-1),k=-1)+np.tril(triag_value*np.ones((N,N)),k=-2)\n",
    "make_theta(1,2,3,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "N=100 #network size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create initial condition ensemble approximately uniform on unit hypersphere:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARMAAAEBCAYAAABbt0EoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xd8VFXawPHfnZI+6aRTkgATaoDAUtRdpNhQ11Usr65YVnHtXWxrr+z6btF9Xd1dV1nXgrHsquAiKjZAIPSSSwklIY30npnM3PePmWAIk8xMMjN3yvl+PnwId+bmPvcyeXLuOeeeR1IUBUEQhMHSqB2AIAjBQSQTQRA8QiQTQRA8QiQTQRA8QiQTQRA8QiQTQRA8QiQTQRA8QiQTQRA8QiQTQRA8QiQTQRA8QiQTQRA8Qqd2AO7KzMwMB6YBFYBF5XAEIRhpgXRg49GjRztd3Sngkgm2RPKt2kEIQgg4DfjO1TcHYjKpAPjggw/IyMjo9407duxgwoQJPgnKW8Q5+I9gOA9XzqG8vJwLL7wQ7D9rrgrEZGIByMjIYOjQof2+sbq62ul7/J04B/8RDOfh5jm41Y0gOmAFQfAIkUwEQfAIj93mGI3GHKBOluUGB68tBBqAHGC5LMsNjrZ5KhZBEHzPIy0To9E4D3gPW2Lo/VoOkCPL8mpZll8Fnne0zRNxCIKgHo8kE1mWVwOb+ni5uwXSbV4f2wRBCGC+6DNJAup6/Duxj21CEFIUBbFoeWgIxKFhQSVWs5nOmho6q6rpqD5G57FjdDU3Y25qoqvJ9relrQ2ryYy1y4xi7sJqMgGwVqdD0miQdDoknQ5dVCTa6Bh0MdHoYqLRGwyEJSURPiSZ8ORkwockE5aUhEavV/msBVf5IpnUcmLLo66PbW7ZsWMH1dXVTt9XVFTk7rf2O74+B8ViQTlWg1J9DGt1te3vY8egsenkN0dFIkVFQaT977g4JL0OtFoknQ6tVmt7n9V6/I9isWDq6ISODpSaWigrQ2lrh/b2E7+3JCElJSIlJyMNSUYzJBkpLc22TZK8fyEcCIXPU1VV1YC+r9eSidFojLeP0BQCN3RvAzb3sc0tEyZMcDr5pqioiIKCAne/tV/xxTmYGhpo3iPTLMs0y3tp2X/geItC0umIzMokKj+fyIx0wlOGEJGaSnjKEMKTkpC6k4UHzsHS2YmpppbOY8forKmho6KSttIy2krL6Ni7D4vVCoA2OhrDqJHEjBqJIc9I7Ngx6KKiBncRXBAqn6fS0tIBfW+PJBP7MO9UoMFoNCLL8magCMiVZbnEaDQesI/4TAGutw8Nn7DNE3EIrrF0dtK0ew8NW7fRsGUrbYePALbEEZObQ9pZZxAzahTRI4YTkZGORuebu2FteDiRmRlEZp78mITVbKa97CgtBw7QvHcfLXv3U/b+h7bWjkaDYdQo4vInEJ8/EYNxtLg9UoFHPiWyLBdia2303Jbb4+tX7V+u7m+b4D3mpmbqNmygdu16GnfsxGoyIel0xI4dw/BFvyR27BhicnPQhIWpHapDGr2e6OwRRGePIHXeXAAsHR00791H47btNGzfQVnhB5QtL0QbHUVCQQFJM6aTMGUS2shIVWMPFaIDNoiZm5up/X4dNWvX0bhjJ1ithKekkHrmfBImTyJ2/Di04eFqhzlg2ogI4idOIH7iBIYDXS2tNO7cSd3GTdRt2ETNN98i6fXET8onZfZPSfzJNL9NlsFAJJMgo1gsNGzbTtXqL6n7YQNKVxcR6WlkXXgBSTNnEJ2bo1rnpbfpYqJJmjGdpBnTUSwWmvYUU7t+A7Vr1yFv3IQ2OorkU2aRcvpsDGPygvY6qEUkkyBhamig8rNVVK1ajam2Fp0hhrQzzyBl7ulE52SH3A+OpNUSN34ccePHkX3NIhp37qL6q6859s13VK1aTWRWFukLzmbI7J+hixK3QZ4gkkmAaykpoeLjTzn2zXcoXV3ET55E9q+utjXpRSckYEss8fkTic+fiOWG66hZu47KFZ9R8spfObzsTVLmnE76grMddvwKrhPJJAApikLj9h2ULi+kaecuNBERpJ4xj4xzF4gfCCe0kZGkzp1D6tw5NMt7qVixksr/rqJixUqSZs0ga+FFxORkqx1mQBLJJIAoikLD1m2UvvsezXuKCUtMZMTVi0idPw9dTLTa4QUcg3E0BuNoRlxzFRUff0rFpyup/X4dCVMLGHrJQgzG0WqHGFBEMgkQjbt2cfiNN2mW9xKWlETODdeTOm+OGJ3wgLD4eIZfeQWZv7iAihUrKf/PJ2y/7wESp09j+C+vIGpYYK+u5isimfg5a00te555jrofNhKWlEjujTeQMvd00R/iBbqYaIZespCM8xZQ/skKjn7wEVtuv4uUOacz7PJL1Q7P74lk4qfMTc0cefsdTCv/S2NEBMOvvIL08xYE9LyQQKGNjGToxReRduZ8Spe/T+XKz6j55ls0s2ZgnTBBtAb7IJKJn1EUhWNfreHgP5bR1dKCtmAyU269hbD4OLVDCzn62FhyrruGjPPO4dDr/6R2zTdskfeSs/g6EqZMVjs8vyOSiR9pKy3jwF9epWnnLgxGI7k3Laa4tlYkEpVFpKaSt+QeNix/D776ht2PP0XijOnkLL6O8CSxFE83kUz8gGKxUPbBR5S+sxxtRAS5N/+a1HlzkTQaqK1VOzzBTpubw+RfXMDRj/5D2fJCttx6BznXX8uQ2T8LuUmBjohkorL28nL2/eFFmuW9JJ0yi5zF14mWiB/T6PUMvfgikmbOYP+f/sy+P7xI7dr15N50A2EJCWqHpypR6kIliqJQ8elKtt5+N+1Hyxl9953k3Xe3SCQBIiorkwnPPsmIa66iYes2ttx6B7U/bFQ7LFWJlokKulpa2fenl6j7YQPxkycx8tabCE9KUjsswU2SVkvmBeeTMLWAvS/8nuJnniP9vAWMuOrKkBy6F8nEx1r2H6B46e8w1dQy4tqryTj/XHG/HeCisjKZuPRZDr2+jIqPP6VpdzHGe+8kMj1d7dB8Stzm+IiiKFSs/C/blzyI0mVh/DNPkvnz80QiCRIavZ6c639F3gP30VlVxba776N+8xa1w/IpTy7b6LA6n32N13qgxL4pESgAltj/TIXjtXeCkrWri4N/e43Klf8lfspkRt95G/rYWLXDErwgacZ0orOzKX72eXY/+QwjFv2SjAvOD4lfGoNumbhQnW+qLMuSLMu59qUcl8iyXIKt8NZBYH4wJ5KulhZ2P/E0lSv/S+aFFzD24QdEIglyEakpTHjuaZJmTufQ68vY94c/HV+gO5h5omXSb3W+nonCaDQu7rH26xL72rFBq72igt1PPENndTUjb7uZ1Llz1A5J8BFtRATGe++mbMT7HPnX23RWH2PMQ/eji4lROzSv8UQySeLHWxjoozqf/XanZ32cnB6r0xfaWysu8/e6OdaKSkxvvQNWhbArLqMsPo6yAcYSCrVaAoXb55Gbjf6iC2j66GM23HE3YZdfhhSnbss04OrmOLAY6G6VIMvyUgCj0VgHvALMd+eb+XPdnMYdO9nz5tuEx0Qz9rFHiMrKHPD3CpVaLYFgwOdRUEDDpEkUP7sU5c23GfvYw0QNG+b5AF3gzbo5nhjNcbU636XdHbNGo3Gh0Wi8z769BFvHbVCo/WEDux5/irDkJCY898ygEokQPOInTmDCM0+iWK3seOA3tOw/oHZIHueJZFII5MKJ1fnsX/fUM2F0V/rr3u52RT9/VLtuPcXP/Zbo7BFMeOYpwpPFRDThR9HZI5j4/NNooyLZ+cjjNO/br3ZIHjXoZGLv6+iuzreYH6vz9b4xK+mxz2pgnn1I+VKCoKJf7fofkH/7vxhGj2Lc44+ijzWoHZLghyJSUxn/9BPoYqLZ9ejjNMt71Q7JYzxV0c9Rxb7cXu8p6GOfgB/Rqf1hI/LSF4gZmcvYRx8WpROEfkWkpDDh6SfY+fCj7HrsScY98SiGUSPVDmvQxAzYQarfshV56e+IzsmxJxLvF9AWAl/4kCGMf/pJdIYYdj/+FG1lZWqHNGgimQxCy/4DFD/3WyKzMhn32G/QRYsV4gXXhScnMe7xR5C0WnY98gSdx2rUDmlQRDIZoPaKSnY/8TT6WANjH3lYlJoQBiQyPZ2xjz6Mpb2dXY89gbmpSe2QBkwkkwEwNTSy+zHbMN/YR38jlu4TBiUmJ5sxD99PZ/Uxip9ditVsVjukARHJxE1Ws5ni55Ziqqtj7G8eFPNIBI+IGzeOkbfdQtPuPRz4v1dQFEXtkNwm1jNxg6IolLz6N5r3FDP6nrtExTfBo4acdgrtZWWUvrOcyKFZZF14gdohuUUkEzdUrviMqlWryVp4IUNOO0XtcIQgNPSyS2gvO8rhZW8SNWwoiVMD5zEEcZvjosZduyj522skTJvKsCv+R+1whCAlSRIjb7uZ6Jxs9v3+T3RUOX+Y1V+IZOICc2Mje3/3ByLS0hh91+22EhSC4CXa8HDy7rsbRbEiL/1dwHTIip8KJxSrlX1/fAlzczPGe+8Sk9IEn4hIS2PUrbfQsv8AB197Xe1wXCKSiRPl//mE+qLNZF9zFTE52WqHI4SQpJnTyTj/XCpXfEbtuh/UDscpkUz60VJSwuFlb5I4Yzpp55yldjhCCBq+6JdE52Rz4OW/YGpoVDucfolk0ger2cy+P7yILjaWkbfcGBILAgv+R6PXM/rO2+hqa+fAn1/26/knIpn0oXR5IW2HjzDy5l+jN4jlBAT1RA0bxvArL6duw0aqv/hS7XD6JJKJAy37D1BW+AEpc2aTOG2q2uEIAhnnnUvs+HEcfO11TPX1aofjkEgmvSgWC/tf+j/C4uPJ/tU1aocjCABIGg0jb/o1VpOZg3//h9rhOOT1Ilz211+hV8EtZ/uopWLFZ7QePIRxyT1BXZZACDyRmRlkLbyQ0rffJWXuHBImT1I7pBP4oggX9Cq45eI+Pmeqq+fIW+8QP3kSSTNnqB2OIJwk66JfEJmZQclfXsXS2al2OCfwxG1Ov0W47JbIspwgy/ISN/bxuUNvLMNqMpGz+Fdi9EbwSxq9ntwbb6CjsoqjH/5b7XBO4KsiXCcU3HJxn355ugiXtewopjXfoD11FrsrKqCiwt2QvCYYClgFwzmA/5yHZkwepYUfUJ2WimRw73Y8oItwOSi4NejSFp4swqUoCjsLP0SJj6fglpvQRvrPgtDBUMCq5zmsKSpl2co91NS3k5wQyaKzxzC7oP//R3/hT/8X7RmZbLnlduJ27GLUrTe5vF9AF+Hqo+CWq4W7fKJ+4yaadu9h6GWX+FUiCTZrikp56b1tHKtvRwGO1bfz0nvbWFM0sA9vKItMTyP9nLOo/uJLWg8dUjscwDdFuBwV3HK4jxoUi4VDy94kIiOD1Plz1QojJCxbuYdOs+WEbZ1mC8tW7lEposA29NKL0UVHc/jNt9UOBfBBES5HBbf62cfnjn3zLe2lZQy/8nI0OrFWlDfV1Le7tV3ony4mhoyfn0f9xk20HChxvoO34/HEN3FWhMtRwS1H+/iaYrFQVvgB0dkjxFCwDyQnRHLMQeJIThC3lgOVvuBsyv/9MaXvLGfMQ/erGktIz4CtXf8D7WVHybr4IjEU7AOLzh5DuF57wrZwvZZFZ49RKaLAp4uOJuP8c6nbsJGWEnVbJyGbTBRFoXR5IZFZmSTNmK52OCFhdsFQbrk4nyEJkUjAkIRIbrk4P2BGc/xV+rnnoI2OpvRddSvthmwnQX3RZtoOHWbU7bciabXOdxA8YnbBUJE8PEwXHU36OWdRVvgBHZWVRKSlqRJHyLZMKj5ZQVhiIsk/PVXtUARh0NLOPgtJq6X8kxWqxRCSyaSt7CgNW7aSdtYZYgRHCArhSYkknzqLqs+/oKu1VZUYQjKZVK5YiaTTkXrmfLVDEQSPyTjvXKwdHVStVmcBpZBLJl1t7VR98RXJp51KWHy88x0EIUDEjMwlZvQoqj5frcryjiGXTGrXrsPa0UGaaJUIQSh1/jzaS8to2bvP58cOuQ6D6i+/IiIjHUOeUe1QQlYgP+zn75JPPYWDf/8HlatW+7wWdki1TNorKmnatZuUOaeLSWoqEQ/7eZcuKpLkU2dR8933WDo6fHrskEomx9Z8DZJEyuyfqR1KyBIP+3lfyumzsXZ0UL/Jt2uvhNRtTs133xM7bizhQ5LVDsWvefM2RDzs532xY/LQx8dT8/06kk89xWfHDZmWSduRUtrLjpI8a6baofg1b9+G9PVQn3jYz3MkrZakWTOo31SEpd13STpkkkntuvUAJIrncPrlzduQNUWldHR2nbQ9XK9lWl4K1z61ivPv/jfXPrVK9KEMUvKps7CaTNRt8t1SQSFzm1O7bj2GPCPhSW4vNxtS3LkNced2qLvF0ztRGaL0nJafwRebyo6/1t0aAsQozwDF5uWhT4indu06hpzmm1sdX9XNWWz/sgB4XpblEke1dDwRiyOdx2poPXiIEVcv8tYhgoara468XLiVFesOH/+3swTgqMUDEBGuY2NxtcPW0O/f2cL/vrVZDB8PgKTVklBQQO26dSgWi08eZvV63Rz7amqb7K99jm1BaehVS2ewcfSnYetWAOKnTPbmYYKCK2uOrCkqPSGRdOtOAI5uVfpr8ThKXgBWqyKGjwchYcpkLK1tNMt7fXI8T7RMnNXAiQfmY1vndTO21gvYaun4ZAGG+i3bCEtMJGqY+M3mTPdv/963LwDXPrWKmvp2JE3fc3SsVts07p4tFQN9t3hiovQ0t5mdxtXdbyNaJ66Lz58IGg31RZuJHev9Bai8XjdHluVCo9HY3fKYx4/LNJ5QS8e+LqzHKRYLjdu3kzhtmpio5qLea4707u9QrK4999FptvDCW7YOQEkCCei5Z+8WkDNi+Ng9uphoYvOM1G/ZyvArr/D+8bx+BKBHH8p8WZYvtm/rXUvHrYdlXC3CtWnFSrqaW6iPj/WbAkruUjvuv31U4bC/wx29nzuTJMhI0HCw2nmrpFtslFb1a6H28d1lTk7C8u33bFq7Fik8HPDvIlwu1cCx18653v71Qmz9LEv5sZaOW1wtwpWJxEEgf8ECwpOT3D2M6vyh8FPTW54vQ6kouJVIAE7Jz6KgQL1i3f7wf+GuBq2OXd98R25EJAlTJvt9ES5ndXO6k8ersiw32G9tHNXS8YrmYpnwIckBmUh8ZU1Rab9zPPxlQtnGYuctUeFEBuNo0Gho2lPs9WN5vW6O0WicAvwVKDIajfXADY5q6Qw2jr407SkWTwj3w5UZr9PyUtQLsAfRZ+I+bWQkMTnZNO32/rNPXq+bI8vyZiChn328NqKjNDZhqq3FkJfnrUMEvP5mvHZ3wvpLiyAmSg+IJQzcFTN6FNVfrkGxWr16nKCeTm+trAQgJtftLpmQ4cqMV1daBL4YKDOZLWIJgwGIycnG2tFBh/3nwVuCOpkoVbbfqFHDh6scif9y5cE7V/pMfLFKYKfZygtvbRZLGLgpOsf2y7S15KBXjxPUycRaXU1EWiq6KP/oQPRHfc147fngnaOH8/yN6E/pW9SwoUhaLS1eTiZB/aCfUlVN1MiRaofh1xzNeJ2Wl8KqDUfostiaG67MUFWbv4w4+SONXk9kViZtR0phnPdmwgZtMrGazSh19WIKvQt6z3i9/DcrjieSQCHqFfcvMiPdlky8KGhvczprakBRiEhXp1RiIAuElkhvYjSnfxEZGXRUVnl1RCdok0lHpW1KcERaqsqRBI7uyWuBSIzm9C8yIwPFYkFpaHD+5gEK2tuc48kkVbRM+tJzvkZMlJ62ji4sLj7E529e/WiHaJ30IyLd9ktVqW/02jGCNpl0VlWBVktY4knz5QROfhI4EG9tegr0+L0tPMn2OInS3Oy1YwRtMjHV1yMZYpA0QXsnNyh9rXwWyNYUlfLqRzuOJxZDlJ7FF0wQLRYgLNH+LG5Tk9eOEbTJxNzQCFHRaofhd7pvbfpa3SxQhes1/PHdLSeMQjW3mfnDO1sA0UGrCQtDFxvr1ZZJ0P7aNjc2IcWIZNJTz6nowSZMr3U4nG2xKmJ2rF1YfBxKa5vXvn8QJ5NGpOgotcPwK8F4a9Otvz4TMTvWRmcwoHixjk7wJpPmZogUsyJ7CtUfKjE71kYXEwNerD8clMlEURQUsxlJr1c7FL8Sij9UWo0kZsfaacLDoMt7z1kFZTKxmky2L3RB2788II4e6gtmhig9d1w2OeQ7X7tp9GEoZu8lE18V4TrpdWf7DIZitt8/60LnB8cVPR/qC8ZO2J4MUXreevIctcPwK5owvX+3TFwownXS6872GSyrqTuZiJZJb7MLhvLaw2cQ7EU/2ju7xBT7XjR6P08mOC/C5eh1Z/sMitVsu82RRDLpU7D3n3RZxJBwb5JeD2bvzRT2RDJJ4sTyFr0rgzt63dk+g9JlH0u3Nnlvgk6gC4X+k1AdvepLR0UlKApWLz05HLC/uvsrwmWtrwdAMhgCrmiSI944BwOwYGosX2xrorHNgiT5ZulFX/JG0a5A/jyZFStoJLZs2dLv+/y5CJej110q3NWf/opwddbWsQnAYgm4okm9ebPwU0EBXLPw5GdagkG4Xst1F+RT4MGRnEAswtXTwW07KNfpA7oIl6PXHe7jKZow+/wSL3Y2BYvuKfbBkEgMUXokYEhCJLdcnC+GhHuxmsxeHeEcdMtEluUSo9F4oEcR8p5FuHIdvW4fGna0j0dowsJsX4hk4lQwTbEXQ8H9s5pNXh3h9HoRrn5eP2mbp2i6L5hIJk4FSyelRhPsg92DZzWZvTrCGZQzYCWtFkmnQxHJxKlgGSK2BugKcb5kNXm3ZRKUyQRAFx0NXnxCMlj4Sx3hwdJoJDFJzQlLayuEh3vt+wdtMtF7ee2GYLCmqJQvNpWdtN0XpT49zWpVRJlQJ8zNzUhR3luWI3iTSZxIJs701fkaqPNNOs0WXnhrM9c+tUokFQfMjY1IXqxuGbzJJD4OWlvVDsOvBUvna2+imPnJrF1dtqVMDQavHSN4k0lcPEpLC0qg/pr1gWDpfHVEFDM/kbm+HhQFKVYkE7dFpKaA2Yy50XurcQe6YH8+J1hbXgPRWWubZC4ZYr12jOBNJvZKfp0DfM4gFMwuGMotF+cH7RyNYG55uau7KJ0UL5KJ2yJSbcmk+yIKjs0uGMqdl00O+BZK73QYrteK5Rp76CgvB0lCSvBeUbqgTSbhqbb5Ex2VlSpH4v+6WyhDEiIDctGkcL2Ws2cOPx6/eDbnZO3l5YSnDPHqDNiAXYLAGW14OFJ8HG1HRI++K2YXDD3+w9e7dKi/6zRb2FhczWsPn6F2KH6rvbyCyIwMvLc2fRC3TACklBRaDx1WO4yA07OlEihEZ2vfFIuF9tIyIrOyvHqcoG2ZAEipKbR/vw6ryfTjk8SCS3q2VABeLtzKinX+m5hFZ2vf2o+WYzWZiMnJpt6LxwnqlokmNQWsVtrKTp4yLrjnxoWTMET5Zx0i0dnav5aSgwBE52R79ThBnUwk+4hOy/4DKkcSHE7Lz1A7hJNoNJLobHWi9eBBJL2eyKxMrx4nuG9zEhPQGQw0F+8l7Yz5aocTsNYUlfptrR3FqohE4kTL3n1Ejxjx4zo/XhLcLRNJwpBnpLm4WO1QAlb3yI4vE4k7Hb+ir6R/VrOZ5n37iR2b5/Vj+aqi32L7lwXA8/alHF8BlgBTAWRZ9viKawCxY/Ko37gJc1MT+ljvzf4LVr5e1lGntdUGdqUlJPpKnGvZfwDFbCZ2jPevky8q+s0DNtlf+xx4xf7SPOAgMN9biQTAkGcEoGmPaJ0MhDeHXA1R+hM6dQ1Rem6/1FYb2NFzQzqtJBaNdlPTbtvDjoHSMnFWnS8emI9tBfrN2FovAEtkWS4c6EH7q5vT076WZtDr2ff5ag4FaO1hNWu1xEZpaWzzfMtEr5WYnx/DxOzoXq9UU1RUfVJdn7goLXPzY3u93/ZeXwq0ujmm79ciJSexff/+49ucnYOadXOSgJIe/z6hOp8sy4VGo7G75TGPHxeQzumxOn2hLMs9v4dT/dXN6VZUVMTU6dPZnT+R9qPlAVnzRO1aLddRygtvebQSCUMSIll09hinrYpmSgkr3oPU1k5YeBjZ2dkerYPjLrX/L9xl6exkQ2kZaWfOJ8cetyvnMNC6OT4ZzenRhzJfluWL7duWAhiNxjpstz5eG26Jn5RP/aYiOqqqjj8AKLhmdsFQ9hysdWnCmiRBTKSeljYzMVF62ju76LL8uJ5MuF7r8q1J7yn93QsedcckONe0azdWk4mEKZN9cjynycTeeRrv4KXVsixvxsXqfEaj8T7s9XHsHbY59oRSwo+3Pl4RPzkfgPrNW0k/+0xvHioo3bhwEmOyk1i2cg819e0kJ0QyLS+Fb7eVHy/eZYjSs/iCCSf8oPccUna1NdLNUcdv94JHIpm4pr5oM5qwMGLHjfXJ8Zwmkx71bfpSCNwAJ1f0626R2JPHq/biW/Ow9bF095fk4OGKfr1FZmYSkZZK3Q8bRDIZoN7T68GWZFzZZyC3B311/IpncFyjKAr1mzYTN2EcWi+uSN/ToEdz7H0d3dX5FnNiRT+MRuMU4K9AkdForAdusI/ezLMnmUvxcEW/3iRJImnmDBq376CrpcWbhxI8pK/5I2JeiWtaDx6ko7KSxBnTfXZMr1f0s98KnbQiS499Bjyi446kmTM4+uG/qduwiZQ5s31xSGEQFp095qRlEMS8EtfVfr8ONBqSfJhMgnoGbE8xo0YSlpRE7br1aociuKD3gk1iXonrFEWh5ru1xE+c4NOJmkH9bE5PkkZD0qyZVK78DHNTM3ovrtIteIajfhrBudYS2y1O1sJf+PS4IdMyAUiZMxulq4uab79TOxRB8JpjX3+DpNORON13tzgQYskkJieb6OwRVH/5ldqhCIJXWM1mqr/6msRpU33e+g6pZAKQMud0WvYfoPXwEbVDEQSPq9uwia6mJlLP6P1Ui/eFXDJJ/ulpSDodVas+VzsUQfC4qs9XE5acTHz+RJ8fO+SSSVh8HEmzZlL9xVd0tYkJUELwaK+opGHrNlLnno6k9f1DrSGXTAAyzj0HS3s7x75ao3ZfgeAvAAAQ20lEQVQoguAxFZ98iqTVknaWOrO8QzKZGIyjiRk1iopPV6BYrWqHIwiD1tXaStXqL0k+9RTCEr1Xta8/IZlMANLPPYf2o+XUb96idiiCMGhVq7/A2tFBxvnnqhZDyCaT5FNnET4kmbLl76MoivMdBMFPWc1myv/zKbHjxhKT69UH8PsVsslEo9ORedEvaJZlGnfsVDscQRiwqtVfYqqpYeglC1WNI2STCUDq3DnoExIoe+99tUMRhAGxms2UFX6AIc9InArDwT2FdDLRhIWR+Yuf07h9B427dqsdjiC47Xir5LJLkCRJ1VhCOpkApJ11BmGJiRx+45+i70QIKFaTydYqMRqJn5Svdjg+q5tzUo0cZ/v4ijY8nKH/cykH/vwydet/IGnmDDXCEAS3lX/8KaaaGkbdfovqrRLwQd0cuxNq5Li4j8+kzj2dyKwsDi37F9auLjVDEQSXmBsbKSv8gIRpBcRPnKB2OIBnbnOc1c0BW42cBFmWl7ixj89IWi3DF/2SjvJyqv67Ss1QBMElR95ZjqWjgxFXL1I7lOO8XjfH7oQaOS7u0y9Xi3C5WjRJ0UposkdQsuxfHDXEIEX3Lg6lnkAr/ORIMJwD+Md5WKuPYVr5X7QFk9lTVQVuFs3y5yJcTjmokTPo1ehdLcLlzqrobWlpbL39buK2bGfU7bcMNkSPCLTCT44EwzmAf5yHYrWy48HfYI2JYcrtt7q9LKOqRbgGWzenjxo5n/e3j1qisrLI+Pl5HH3/Q1LPmEfsGO/XZxUEd1St/pLmPcWMvO1mn67v6gpf1M1xVCPH4T7+YOglCzn29bcc+L+/kP+/v0Wj1zvfSRB8wNTQyOE3/kns+HGkzDld7XBO4vW6OY5q5PSzj+q0ERHk3riYtiOllL6zXO1wBOG4g39/DUtHB7k3LvaLoeDevF43p9frhf3t4y8SpxaQMncOZR98ROKM6RhGjVQ7JCHE1Xy/lppvvmPY5ZcRlZWldjgOhfwM2L5kX3s1YQnx7Pvji1jNZrXDEUKYqa6eAy+/QsyoUWQtvFDtcPokkkkfdDHRjLzlJtpLyzj0xj/VDkcIUYqisP+lP2PtNDHqjltVWY7RVSKZ9CNhymTSF5xDxcefUvvDRrXDEUJQxSefUl+0hRFXX0lUVqba4fRLJBMnRlyziOjcHPa/+BKdx2rUDkcIIc3yXg79YxkJ06aSdvZZaofjlEgmTmj0eoz33InV3IX8wu/FszuCT5ibmpF/+wJhSUmMvuNWJI3//6j6f4R+IDIjg5E3/5rmPcUceu11tcMRgpxitbLvjy9iqm/AeN/d6GJi1A7JJSFTuHywhvz0NFoOlFD+0X+IGjGCNBUqpgmh4fCbb1G/qYicX18fUNMSRMvEDSMW/ZL4SfmUvPJXmvYUqx2OEISqv1zD0fc/JO2sM1SrfzNQIpm4QdJqMd57F+EpQyh+dikdlZVqhyQEkaY9xez/88vETZxA9vW/8stZrv0RycRNupgYxjx0P4rFwq7Hn8Lc2Kh2SEIQaD9aTvGzzxOeMgTjfXej0QVeD4RIJgMQlZXFmIcfwFRTy+4nn8HSLmoWCwPXWVvLrkcfB2Dsww+iNxhUjmhgRDIZoNgxeYy+5y5aDpQg//YFMWQsDIi5uZldjz5BV0srYx95mMjMDLVDGjCRTAYhafo0cm9cTH3RFva+8AeRUAS3dLW1s+fJZ+ioqCTvwSXEjMx1vpMfC7wbMz+TdsZ8LO0dHHrtdfZKEsa77/Dr5ycE/9DV1sbuJ56med9+8u67x28WhR4MkUw8IPPn54GicOgfb7BXq2H0HbeJhCL0qau1ld2PP03zvn0Y77mTpJnT1Q7JI0Qy8ZDMC85HsVptxbwsVkbfeZtYpU04SVdrK7see5LWAyXk3Xd3UNVp8noRLvuyjPX8uBp9IlCArSjXCYW5PBGLmrIuvABJo+HQP95gT2srefffizYyUu2wBD9hqqtn9xNP01ZaivG+u0maERwtkm6+KMI1VZZlSZblXPvqa0vsyzaeUJhrsHH4i8wLzmfkrTfTsH0HOx95HHNTs9ohCX6grewo25c8SHtFBWMeuj/oEgn4oAhXz0RhNBoX91iusXdhrqCROm8OefffS+vBQ+x44CExUzbENct72XH/Q1g7Oxj/1OMkTJmsdkhe4asiXN23Oz1LWpxQmMveWnGZp4tweZxOi/7yS2l/t5CiO+4h7JKL0AwfNqBv5Q+FnwYrGM4B3D8Py+49mD/6GMkQg/6Ky9jb1AgqX4uALsJltxg4XjbDQWGu+e58M28U4fK4ggLaZ0xnz1PP0vGvd8i98QZS581x61uofg4eEAznAO6dh2K1cuTtdykr/BBDnpG8B+4jLN5R+SnfCugiXD1c2iOBOCrMFZQiMzKYuPRZipe+wP4X/0zroUOMuOpKMdITxLra2tj3+z9Rt2EjKfPmkvvr60Pi/9sXRbi69UwYjgpzBS1dTAxjH3mIQ68vo+LjT2mW95J3392EDxmidmiCh7UePIT82xdor6gk+/pfkb7g7IB7+negvF6Eq4eSHvucVJhrsHH4O41OR85112K87x7aS8vYeuc91G0Kjn4EwbaKfMXKz9h27/10tbUz/olHyTj3nJBJJOCjIlz2fxf0sU8hIST5lJlEZw9HXvoCe558hvTzzmX4lZejDQ9XOzRhgLpaWtj/0svUrltP/JTJjLr9VsLi49QOy+fEDFgVRGZkMOH5Zzj8xj+p+PgTGrZsYdQdtwXUEn2CTd2mIg78+S+YGxsZftWVZF5wfkAs/uwNoXnWfkAbHk7O4usY9/gjWNo72H7fAxx56x1RPTBAdLW0su+PL7LnyWfQGWKYuPTZ4zOgQ5VomagsflI+k//0e0r++ndK332Pmu/XkXvTYuLGjVM7NMEBRVGoW7+Bklf/hqmhgayLL2LopReHxGiNMyKZ+AFdTDSj77yN5NNOoeSVv7HzwUdImTuHEVdfqXZoQg/W2jp2P/E0DZu3EDV8GHkPLhG3pj2IZOJHEqcWEDdhPGXLCzn64b+p27ABTpmJdeJE8ZtPRZaODsre/xDT+x9iCQsj+1fXkHbOWQG5Tqs3iavhZ7Th4Qy/8gqG/Ow0Sv76Go2ffc6W7TsZsehKEmf8JKSGGtWmWCxUff4FR955F3N9A5oJ45hy152EJSaoHZpfEsnET0UNG8a4Jx5l47vLkb5bS/FzS4kdO4ZhV/wPceNFf4o32fpFfuDQsn/RUV6OYUweeUvuZV9bq0gk/RDJxI9JkoR21EgmX7yQqtVfcOTtd9n50CPEjh/H0EsvJm7CeNFS8SDFaqV23Q+UFb5Pa8lBIrOyyHvwfhJ/MtV2nYPkYUVvEckkAEhaLWlnnsGQ2T+jatVqyt7/kF2/eYzYsWPIvOgXJEyZHNJDkoNl7eqi5pvvKHv/A9rLjhKRkc7I224mZfbPxPKbbhDJJIBow8PJOG8BaWfOp+rzLyh7/0P2PPkMERnpZJx7DilzThcru7nB1NBI1arPqfxsFabaWqJGDGf0PXeRPGuGSCIDIJJJANKEhZG+4GxSz5xP7dr1lP/nE0pe/TuH//U2KXPmkDpvDtEjhqsdpl9SFIWWvfuoWPEZNd99j9LVRfykfHJvXEzC1AJx2zgIIpkEMI1Ox5CfnsqQn55Ks7yX8v98QuXKz6j4+BOic3NJnTeHIT89FV1MjNqhqq7zWA3Hvv6G6q/W0F52FG1kJGlnnkHaOWcSlZWldnhBQSSTIGEwjsZ4712Ym5o49vU3VK3+kpJX/srBv/+D+MmTSJ41k8SfTA2pxGKqq6duw0Zqvl9L446doCjEjh1Dxs3nk3zqLHRRUWqHGFREMgky+thYMs47l/RzF9BacpDqr76mdu066jduQtJqicufSNL0nxA/OZ+I1FS1w/UoRVFoP1pO3YaN1K3fQPPevaAoRGSkM/SyS0iZ/VMi0tLUDjNoiWQSpCRJIiY3h5jcHLKvvYqWffupXbee2nXrOfDyKwBEZKQTPymf+En5xI4Zgz428Apmd9bW0bh9O43bd9CwbQem2loAonNzGHb5ZSTN+AmRQ4eKvhAf8FgysZe8qOu1ulr3ayfV1emv1o7gWZJGg8E4GoNxNMOvupL2o0dp2LKNhq3bqP7iKypXfAZAZGYGBqMRQ95oDKNHE5mV6VfT+C2dnbSWHKR57z5a9u2jee8+Oqtsi4rrDAbiJo4nfuJE4qdMIiIlReVoQ4+ninDNw1Yv53p6LcHYo65O9/qvrxiNxud7b8O+9KPgXZIkEZWVRVRWFhnnLcBqNtMs76W5WKapWKZuUxHVX35le7NGQ2RGBlHDhxE9YjgRaWlEpKYQnpqCPi7OK7/tFYsFU30DnceO0VFVRduRUtpLy2grLaWjqhqsVgDCkpMxjBpJ+tlnEZc/gegRI8RcG5V5aqW11UajcVMfLzuqq3PAwTZBBRq9nrjx445P0VcUhY7KSlr2l9B25Ahthw/TeuAAtd+vPXG/sDDCkpMJi49DZzCgMxjQxxrQRUejCQtD0uvQ6MPQ6PVYDh6kurEJxWLB2mVB6erC0tZGV0sLXS2tdLW2YG5qxlRTQ2dt3fGEAbYJexEZ6URnZ5N82qnEjMzFMGqUmNbuh3zRZ+Koro5LtXb64/d1czxIlXOIioC80bY/QLjJhNLQaPvT2IBS34ipqRFTWxtKTQ1KWzu0tZ2QCHra52ijXg+REUgRERAZiZSWhtY4Gik2FikuFik+HikxAUWrpRVoBWoADpbYakGqIBQ+T4FQN8ejAqJujgcE0jkoioLVZEIxm7F2/zGZ2bVrFxPyJyJptUhaHZJOizYy0q/6Y1wRSP8XffH3ujnOOKqr42qtHSGASJJkWxi71+LYmqpKMSQbAjxRN8ehHnVzHNXVcVhrRxCEwOWR7m/7MO9U4FKj0TjFvrkIHNfV6afWjiAIAcpTozmF9Kp/07NuTh91dU7aJghC4BID84IgeIRIJoIgeIRIJoIgeEQgzjPRApSXlzt9Y1VV1YDHzP2FOAf/EQzn4co59PjZcmu5uUBMJukAF154odpxCEKwS8f26ItLAjGZbAROAyoAi8qxCEIw0mJLJBvd2UlSFMU74QiCEFJEB6wgCB4hkokgCB4hkokgCB4hkokgCB4hkokgCB4hkokgCB4RiPNMBEFwgbsVIwZ7vKBNJr6+kILz6xoo192F83gFWIJtDR9kWfa7ZTTcrRiBB6pDBOVtjv1Cvoftw9D7te4Ludq+psrzvo7PVUajcaHRaJxnNBoX21ek6/36K0ajMd7+HlVX+Hd2XQPlursY5zxsS1rP98dEAscTnDsVIwYtKJOJGhfS0wLwQ+3sugbEdce1OJfIspwgy/ISH8XkaUmcuO6y29UhHAna25x+DLrMho+4+qEudLBdDc6ua6Bcd1fizLG3BKcAhfZlSENeKCaTQCE+1H6qR19DHfAKMF/diNzmleoQAZdMvFR6IyD52Yfa2XUNlOveb5z2ztnuzssSHPTL+SsnFSMGLeCSiZdKb/icC0kx0D7UDq+rv113Fzg7j+5zAds198vz6FExosFoNGL/TBUBubIslxiNxgM9WrUeqQ4RcMnEFWpcSHe5kBQD6kPdz3X1q+vujAvnsdo+ulYHTMN/z8PtihGDJdYz8WP21ksJtg/1q7IsNxiNxgPdHwr7690f6mf9dd6GEBpEMhEEwSOCcp6JIAi+J5KJIAgeIZKJIAgeIZKJIAgeIZKJIAgeIZKJIAgeIZKJIAge8f+eehYi7TPo1AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "num_points=1000\n",
    "np.random.seed(4)\n",
    "x0_data=np.random.randn(num_points*N).reshape((N,num_points))\n",
    "x0_data=x0_data/(la.norm(x0_data,axis=0)[np.newaxis,:])\n",
    "# print(\"init |data\")\n",
    "# print(x0_data)\n",
    "fig,ax=pl.subplots(1,1)\n",
    "ax.scatter(x0_data[0,:],x0_data[1,:]) #plot 2D projection\n",
    "r=1.0\n",
    "xvec=np.linspace(-r,r,1000)\n",
    "ax.plot(xvec,np.sqrt(r**2-xvec**2),'r')\n",
    "ax.plot(xvec,-np.sqrt(r**2-xvec**2),'r')\n",
    "ax.set_ylim(ax.get_xlim())\n",
    "ax.set_aspect(aspect=1.0/ax.get_data_ratio())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot transients over various matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>iter</th>\n",
       "      <th>x</th>\n",
       "      <th>$\\alpha$</th>\n",
       "      <th>$\\beta$</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>0.095557</td>\n",
       "      <td>0.95</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>0.099261</td>\n",
       "      <td>0.95</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>0.097620</td>\n",
       "      <td>0.95</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>0.100146</td>\n",
       "      <td>0.95</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>0.098375</td>\n",
       "      <td>0.95</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   iter         x  $\\alpha$  $\\beta$\n",
       "0     0  0.095557      0.95      0.0\n",
       "1     0  0.099261      0.95      0.0\n",
       "2     0  0.097620      0.95      0.0\n",
       "3     0  0.100146      0.95      0.0\n",
       "4     0  0.098375      0.95      0.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d=0\n",
    "avec=[0.95,1.0,1.05]\n",
    "bvec=[0,0.005]\n",
    "num_steps=150\n",
    "df_std_all=pd.DataFrame()\n",
    "df_norm_all=pd.DataFrame()\n",
    "for b in bvec:\n",
    "    for a in avec:\n",
    "        df_std=pd.DataFrame(np.array([np.std(np.linalg.matrix_power(make_theta(d,a,b,N),k) @ x0_data,axis=1) for k in np.arange(num_steps)]).T,columns=(np.arange(num_steps)).tolist()).melt(var_name='iter',value_name='x')              \n",
    "        df_norm=pd.DataFrame(np.array([la.norm(np.linalg.matrix_power(make_theta(d,a,b,N),k) @ x0_data,axis=0) for k in np.arange(num_steps)]).T,columns=(np.arange(num_steps)).tolist()).melt(var_name='iter',value_name='x')              \n",
    "        df_std[r'$\\alpha$']=a\n",
    "        df_std[r'$\\beta$']=b\n",
    "        df_norm[r'$\\alpha$']=a\n",
    "        df_norm[r'$\\beta$']=b\n",
    "        df_std_all=df_std_all.append(df_std).reset_index(drop=True)\n",
    "        df_norm_all=df_norm_all.append(df_norm).reset_index(drop=True)\n",
    "df_std_all.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1e72d2fe400>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPMAAAC9CAYAAABmpEsjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzsvXmcHFd16P+9tfW+zL5pRjOjpSXZI8navBsbC2yz2sYsgUB4YecHCbwESELCI+QlBEIILwSwjX8kIYEEvMEDzGJjG7zgRZJlyZbcWkajkWbfenrv2u77o3q0LyNpxjOS+/v59Kerq6pv3ama0+fcc885V0gpqVChwvmPMtcdqFChwsxQEeYKFS4QKsJcocIFQkWYK1S4QKgIc4UKFwgVYa5Q4QKhIswVKlwgVIS5QoULBG2uO3A6EolEHNgIVAPdyWTyoTnuUoUK85LzQTO/DdiSTCbvBD4z152pUGG+MueaOZFIdALjyWQyVf58G5ACOoEfloV4SkN3T7fdlpYWH7AeGACcme53hVlFBZqAZ/v6+kpz3ZnzhTkV5kQisRH4EvABYEtZsDuTyeSXy8fvAD5UPv2DnJlmXg88NoPdrfDyczXw+Fx34nxhToU5mUw+lEgkNh2xa0orT7ERDmnrO/G09ZZpNj8AcN9999Hc3DwDvT3M9u3b6erqmtE2Z7Pd2Wx7Ntrt7+/n1ltvhfIzrDA95tzMPoYajjalq4/Q3ilgE4c19elwAJqbm2ltbZ3RTg4PD894m7PZ7my2PZt9pjI8OiPmmzAfR9l7vWiu+1GhwnxnvgnzGN4U1BTj59rg9u3bGR4ePtdmjmPz5s0z3uZstjubbc90u0NDQzPa3iuF+SbM91A2o8ve6+mOj09KV1fXjJuBmzdvZu3atTPa5my2O5ttz0a7Bw4cmNH2XinMtTf7NmAdkEokEiSTyS2JRGJveZy8Bs/LXaFChWkw197se/C08ZH77ixvViK9KlQ4A86HCLAKFSpMg4owV6hwgTDfHGAVjkG6DtI1kXYJXAvpWkinhHRMpGOBtMB1j/6SooLQEKrhvTQ/wprENbPetlJ57Bcilac6T5DSBaeAkx/BtbK4xQlcM+0JsZg6CRAChIIQKgjV+4wovwNSeidKF6SLdG3AQc3sodhvASBUH4q/CsVfjaJHUIxwRcAvACpPcI6QjuUJbSmFWxjFKU6gTe6iNJQDRUWoPoQWRDGiM3M9fQw1UOdtuzZucRInPwzSARRUfzVKqNF714Mzcs0KLy8VYX6ZkFIirSxOYRwnN4hrpjwtKlSEHkTxVyP1KtRg3az3RSgawggD4UN9c+0CztgLWIAw4mjRNrRAHULVZ70/FWaGigNsFpHSxSlOYI7vonjgUYr9T2CldiGljeKvRQ3WowZqULQAQijl70hcMw2AObaDfPdPAcjv/Qmpp/8WgOxL32f04Y952zu/x+hDHwEgt/fHTDz5OQBKg5vI7bn/cHuufdJ+CiFQ9CBqoB41UA/SwRzdTuHAI5jjSVwrNzs3qMKMcl5o5kQi8UG80M7uZDJ5zlFhs4mUEmlmsHODONmDuG4JIQwUI4KixI4617VyKHqIYt/jFPufhPAbSW/9F/J77qfxtl9T6P01uV13E+h4PU4phZXaDYDQAihGxNvWQyiBGq+9Utl0Bgq9v6LY9wShxbcw+fztRPb9CrniQYoHH8EtpgguvgUxNc4+BkULoGgBpOtgZw9iT+5DjbSixztRtMBs3boK58i8F+ayIG8qR4cdmd88r5COiZ0fxk7vR5ppUHQUI4qmxA+d4+SHKfY/SXDRm8glf0B667/QeNtD2JmDFPseg8Wvwd9yJWrQ047hFe8mvPxdAERWvJvIincDEFp8C6HFt5S3bya0+ObyOe8hsuI9AMTWfYbISk/DB9quZ6QQwnRcsvsexMkepNj8ekrb/hlXupjLPoZluzhSIqcc4wI0VcFQDfyqH//4AfypXkJ1K9CirYcsiQrzhzkX5tNVGsHLmPph+fTqEzYyh7hmBjtzEDtzAKT0NHCw4dBxK7WH7M7vEV39MUpDm5l89u8xai7CqFtFpOsD4NqEEm8nvOwdTOzcia9+Ob76NQCovvjJLnsctutiOi6WLbFcl5zpULR1iuMpClYH+4XKyMFJaP0ThFNEjuSoLjkIAZmcSc22T2NVrabQ7v14SAlF28WVEtuVuFJFSJfY6DMsaBqkrnVtxQN+BOUQZPCKa9x5ypNnibmOzZ5OpZEUnhCnTtrQy42UOMVx7FQ3TmEEFAPFX40QClK62JleUk9/kciqDyMUg9LA77AXvxl/y9XUv+GHqOFWhBAYNSvO6LKO61JyXEzbpWS7ZE2bvOmStxycY+aaNVVBFQJNFUT8GhFDoSo45czyeX/Giv8PCYRcBxFuRQnUYgiH4Lb/RanjXTjxi4/rQ8HSeXH/flpLDouWXH42d++Co/x/vCmZTKYSicQHE4nEmrkYDs51bPZ0Ko18CZhSUeecEnkuSClxi2Oo6e0UB0ZQtBBqsKE8Tk4z+ugnCHa+gUD7jUingLTyGE0rabjlgUNmqWKET3sNy3UxbUnJdsiaDjnTJmc6mLZbnk4WCAG6ItBVQdhQUc5FSyoqheWf9DbzB1EKfSAdRGkMJd+HU7Xy0KkBXcUXr6dvoIdYTfvZX/PCIj5lWQJrX5Ga+QQcV2kkmUzeWf61qwbumItOTQmxOZ70xsPSRSsL8fhjn0ENNRG95I9RQ83lQIwQdTf82xEtCKSUIB2ktMF1QLpYjk3JcjBtT2B7x8bJ7N2PWw4OkUJFUzU0VSOg6YR9p35cUkqEEBzMuLgS2qIKW4Y8L3YV8Mt9FpoC1y/UebDHQhHe9jMDNroClzRoTKgtWBu+g6pp+Pbcha/3XjJX/gfSV3voOooQGL4wg0P7AP9M3+7ziqlCk4lE4tN4/793z1Vf5pswn5Bz+aU71+IEwkqjFHoRVhqpBkH1Y4w/Qc8v/oPCwg/gL/lwXZf+l16CmndBxkJs3wTSRDhFcIu4TgnHMbEdScmVFG1JyZFIygvdS4EiQBdQ6N2PEAKBRCKwkeQclbyjUWO4PJaqp4TBNXUm/36gDleq/H5nkb95IYoA/rIrw19sjWIoks+vzPCPW6NoCnx+JXx3axpNkXS6Gf51axS9vP3V56IEVMnnVmb47NYougKfX5nm0aHrWR5soaEvTdP4vzIeuZKS0QSAbeXQRyZpbV7xSi9OsK5cDWcLHBoazknG33wT5hmvNHK2xQlcM4uV2o2dLaIYnTi5QfJ77ie65hPsG/opkQC0LWmHRX+Ka+ZwzRRuKYNp29iupGTZ5G2NnBPFdAWgQtk0DioKuioQQuBImDQF1T7J/S/mUMM1vKmtyD+9EGasqPK/103yocerKNiC775qhEd7q5k0Be9eth/R5yKkZHFsghtbNRTNT+eCGv7UbyAUjY6GWv6u2kFTBOZomjvfUI0Agnotd7ZIhICQXstX61wcF1qjtbxHetq7fWE1927JsrbhWr7QnEHreQY32k6w4woACvk04XA1ZpFKcYKj6ZyrC883YZ7xSiNnirRLWJPd2OkeQMEc3oKv6VLszAHy+x7A334DZvwywi1VZAc2Y9qSgi3J2Qp5R8Ut31IhdDRFQdcFPkWwP6syVFDZUGfyg+4Ae9I6n105zl88G2V3xs99G8d5eLSa/X1h3rSwiF+FoOY5td69OIcjBQiVL65PE9AkUonxl2u9ktJShnlHxyTCGUFMOmwwhGdFFKppD0WRaoB9oxDSD88rh43D283hw9NMN3Yejvi6/5YwOVMy7Pp47+BXubWmjg/3/wI33AFa0yw+hfOKNZQ1cXkmZk6GgnCGwpxIJKLJZDJ95HYikWhPJpM9Z3Px+VRpRDoWdrYPa2I3TnEcoQUQ0iH11F8TXvlRZMOVqGs/x+DYKINjA+QUiRQGCIGieD8CETlC79gkT5jrecfiIs89/zDR/Et0RVIsKOTplHkC+y1i+sfZn70Uffg3fFO5DxkTsFnjrpCKiAWw+6/iI8tvRJiTGAce5jqjGjdQj1uqJ2pUwbFzvEKA8CEV35ThDq6JyPcjOABo+MwSwqpDaqHDSRmnQVcEcb937r/d3EJMt3Cf+E/2KiuoX/uxGbv35ytlhZMq/7/GKft45qo/pxTmKUFNJBKfwsvZIZFI3IunPTsSicQP8Mziu87m4vOh0oh0bezcAOb4LmRhFBFewNjD70LUrqK06P0UF3+cwVI18uA+LBFg1PLjUzO8mI3DwV/wqvAuCul+otIbEawAPjx4Cbd1Qp2eI6CbuL5qXGMhaRFEDSi8ribMjYEJnNwSSi1vAmkhXItcJkU4oOH6vJGGMMfRhx5FyMOhmFJo2LGLKC75IABK7gBuoBGUY2KoFQOpGOUvOWj2AEp6BwgDN9CENOKg+qZ9n+qCCuDjdv/f8ny6ii9YJSJnfdcvGNbNpfAeyymF+QiN+1AymXwukUh0JJPJfcCfJRKJjvKx6Uc2zCNcu4CdHSA3spOiI8lv/Rpu5gADic+j1d3M3lIDS8b2sym/nN/1ZvizjufoG9iHnh/FbP09dqd1ripsBz1HLriE5+02VjfHKGi13LOqiKZoLFjxRgCKgF5+mUf2IdSGGWo79Hm4vx/tiIL9briD7NqvIaxJlOIISmkYURxB6l4mlbCyhHb8PVKoOKGFONFl3ivU7uU0TyFUXDWE1OMgLZRcL+R6kf4aXH8DaKFp37f3XdqGPvo0/q1fInXZd1GUirk9X5iWmZ1MJp8rv+8Dz8QGYslkcuss9m3Gsc0imcwImcFtTBZK5Hp+Tqz/Rxy86J9IWSt4dnwJr8se5JfFK/jfz1exZcUX2Zjdxeu1LByADi1OXyRBVpr8XmcWteMjFIVCGMnq8jUCCA7lEx+ZZ3y2CIE04jhGHIclRx2SikFh0ftRcvvRMrsw+n+O6H8AJ9BM/uLPeie5Nhw5By10pBEDKRHmJGpxFKnHkMFmpH56XasIQY9Zy670OozeNFe0v3KFeb6tSHqmY+ZP4YVXpvDM7nkvzNnsJC+80MOebT/CF19EcWwnq8fuYE/r35DK1dKXvpQ1Bx+i2h7jMrWf+p6fcWnHn/PVdTZ63kLElpIPLMDxN4EaolZYqCO96JNZ7wKCcsLC1DQTIASH9ggVKXRQDVB8SNUAYSAVDYR2bsKuGtjVl0D1JZ7Gt3No6V3lHGVvWi207a+x410ElATIpsPXE8IbPwM4BZTJHUgjjhtsBe3U+cwtzYv45fDHeXO9e8rzKry8TEuYE4nEq5PJ5MPA5mQy+Q/lfbHTfG1esPX7byQUDbFK38X2/mWU1Bq+W7yVutEhXlX4V4LhHEyARJD11ZKhDmfoGZYKh151OQBuQcMu2UiRRyJIFwzGJxU47G5CoCAEKAJUpRyhJbx3QQmFPAouCi4CT8OhKKD4kXoYoQYpmiY500ERIBDl80AoR2yfSvi1kCfcU7gOVs069PHNtDnP4mZ+iVV3BWbdVUcLrBpAqgGw86ip7Uh/I26gyfsBOgGKELzLuIe65+5n8spHzvrZVJhZpquZX5tIJF4LrEkkEm8FHuQ8KYW7UBugSXexpWCxspcRJ0V9fQzpT3HvyKsYcmIkrYXstVpwUPELk9s37CemlfjxPpXBrE2VkiOujhNVigRVk4X6BI22yoSlM276EIqGUq63pagamqaDGqQo/chDsqdhqFNmt5edJG0XIUuIfBaki5PNkB5KY6thbCWEo/rxHpE84i8SqAqoQvHeFYGmKCiKQBOgKN5+BYGqxbDb3oG64FayPY9SZ+/Ad/DHiNIYpfbfO/5maUGkDCBKo6ilMdzwQqRRfULr4ReppUykbuUPZvyJVThbpjtm/rMjPycSievxFkE/Ky/2y8nmps+iFHZSa+5iqdZNmzYImQcgA50nC5Pu9d7+B8ARviFXCkx0dNtFSUkaJDRIiSpObG66qPQ5dYy6NQzTwJBsZEKpx/DHual5ElSFbZPVhDSHqGYhFTB0A7/MIdxJcMEVfmw9jqsGcMtTTxJvSO5IsGwXiYPrlst/HSP4ZbuflLucycZrCNUMohpBRMkhlHoWX74bs+X1YJQNLSG8sbO0UNN7PE98qO04z3ei82LGF156Zg+jwqxyuqmpKfP6KJLJ5K9nr0szy803vZeGxkZGxgd5bvcefrJvH0NjA4xOZsApEFYKhHRBTUBQ55c0B0yCioXP8OEPRNH81ei+AKh+QKBm93JgPEtdfQNZW2G0ZJCzXIqmRckysSyLDv8YncYwmXyekVSJmDrOMvEkAVH0OmWD2RejaDTzbP/VbCotZ7e9AIlCWLWoNizuXLUFXZH8YrCWnCOo0YpUGzkiQT/VAR2/MTUVVS7op57kBuCJc6YcbZZSG3AsiTQLNKRHiKR+hzb2LKPxa8nUXovfF8SnKeiqBr4qhJVBTb2IG1nkOc7KLBr9HpeO/IxU/W9n6clVOFNOp5nfChwnzC8n5Yn5jXjz2d1n40E0dJ2Whlaa6xdw9ZrLGBgdpn+oj/2jw2wfcdk0orF5WFCwvXHp0jhcFLdYHhxmma8foaoYPpWITyOqVCOcUQDCmktYKx5xpakJqCA2rQSqIbHAO+JKSd6aRCmNQnEEke8nWOjhb6p2AlAkwG5nES+wii1WF7riadifDrewM3t8Ub+vX7yZpdU6j47Vk0wHaAhYNPht7xWw8atHjue9v0tTyp9Uz2zO1G+kGF9F3egDNEz8iqrMJnprbmUksBSJxKepRHw+gpqLb3InItTmjaWF4IH0KiYmg694M7tcPOM15Y8fOCJ7alrHZ5LTCfPby9lKVcDeZDL5kdnqyCl4G948d3cikTinsboQglgwQKxtIYtbWhnJllgxPMLrJgdRCv10T9o8M+rnqUGF+3vgXtlCQGthXXWOtVVplsUs4oqOnbXJjufwaxohn0pAV9CU01TeOGKKichiAPKWw3hqiKg9hN/sI5Hvpku+wO+pUBqooxDo5PZl40wo9YybOuOWwZhpMFLy0erL4Svl6BkP89Ohekz36Ov/UWKUG1sy9OZ0dkz6iJpF/HGFiH70kMAy6uhv/gMC+b00Dt1D68T/pSf6JyBUbBfGciZjSAzVoNbch98uQqS9YmZzqLDGa5LJ5FsTicQavHTdD033+ExzOmF+KJlMvr3csY5EIvH+ZDJ5xuPk01UTOdWv1VSEzVSq2Zle+2RoqkJTLEBTrI1sqZnB9DIYGmBRzT7etyRFSYuyacTg6QGbJ/vCPDbsDZ5bQzZr/T6uDhVYFppkIm8yaEtUAQFDI+xTCeoqqnL6KSdDV4lE6xnNxwlEVjAhJqgP2wQLPQSKPcTSTxNPP029XksuuJRcKIEdmzJ1BQ5h3t9+gA+27mbc0ui3auhzaukrhkjEPIth01iAu/bUAHXQDY1+i0S0xHWNWTbUFg71pRBcRM/C/4lmT4JQ0axxfE4Bzd8CgO1CX8FPsHCAWtth0civuHR0bs1s5U9/8h7gD2ep+e+4X3njd09zzkbgWWAqFHnjGR6fUU4nzIeylpLJ5L5EInHGwjTNaiLT+bX6IPCZM73+dAj7NBbXRemoiTBRaGd/fx+F4RdYH8/yqpZalEv9HEgV+d3ePp4YjfLz4XZ+NKaiC5e1VRmuqEmxLj5BWGQYy5oMuC6aIoj4NEKGhk8TqCfR3EFDZbKk4LiAEFhGHZNGHZOx9ahOjmB+N6FckqrJJ6mafJKCr5VMpIt8YFG5CL6CqwWJa1DtG6VL9uPGA5hGHS4hbm5Nc0Vdnu39BSbUWnZnfOyY9LMoYrKhtsBLkz7+bW8Vl1QXWF1dYElERwHqR35CKLeT4bqbmYxfhqaA5tMwnRgDIwd4Znwp6YqZvYiysJ7l8RnldML8mUQi8UPgV8AmvOSHMxpDT7OayJS2PpLU1Pi4fOxOPE0+a5lUqiKoDfmoXdJJrrWZwYPbGRzci6lWURc2eMeKEO8s7Ke3f5h+J8wTqUaeGI3x1d0LgYXU+0yurE1xeXWKrsg4BavIRMFCFQK/JggYGgFNxdCUQ5pbCEHYpzFZsI7rj6OGyERWk4msRrXThHMvEclup370ARwlSCZ8EZlwF47mjaldxQf4UFwTf+EAruLH9NXTGAiiRiepqzscv+2Uh9R5WyFrK/x7dzX/3g21Pptr6nPcWP8OLnO/S+PwPejWGKO1N4FQMFSB7a+iITRJbfvbZutRTIuy5jyd9nzFcLrY7EngbYlE4gN4QvftGbjmcdVEyte650QnH6HZU3g/KC9Ldc6Q30/nonW0VFcxPriDgZLBmBMlVLLRhcsV1eNc3mDxp4kDDBYNnhiN8cRojAeHqrm/rx5VSFbGslxalWJNbJQGJ0WpYJFRbBRAVRWCmopPV9CEmJpXOimOFmUytoHJ6DoCxf1EstuJpTcRS28iF0wwGV2LZXgF9F3F8OqSuSb+4n5cJYwqj7YMyj4w1tQUWFNTIGUqbBkP8NhwiP97MErccFjQ9ofUDv2ImolHUJ0cQw23gfD8A13mwyzY+yUmmx+bjdt/vrCXU+cvn+74jDLdeeaZEOKzoqydF53t98+10ghSouQLaOZBAkTIpzOUbJv+4VFUvYAQAhW4xtfHNS1gNwt25uJsStfxbKaO2/e1Aq3U6EXWhofoCgyyMjxGUHFIIT3zWkoylktIF9jjKe+hnGLIPUEVaNegh9ZQb+6gJr+LcP4lJtUFDPpWkVePXhVDlcNEsMkOZigqcaQ48WPvUqGrCXL1KpqQjIy6fGPknVzq1vIufsKA08KYkQDgudRK7FI7fwCv5EojD+HlL3+5rHSOdc6e7viMMhfFCWa8msipONtKI0cinVUU+h5D0ULY2TB7XniUWLiBUcuH64JfP3pM3AS8mglggrGSxuOjcX4zEufxsQX8amIhhnBYFU1xZc04l1WNU+czIVMim80Q9gdwAENV8GkKhqaglcM5jydOngUUnWuIZLcRzTxHIv8zcoFFpGKXYxmH63ZNjKeojypAmpLRgKOdPBr3yJ+CZbbOV/b8Po/kuljsW8jNzRmEgA7pY1L1SgK/UiuNTM2wlGdZUpTz7xOJxN5kMrnoZMdni7kQ5jmvJnKmCFVHjy3CGn8JxRdFwyHqDxCLhsgWLUbzFkXLxq8rxzm6anw2b24Z5c0to5iuYPNEhEeH4zwyHOdr3d5KFItDWVaFh1mh9dKluiiAiyRbspHlRSB9uopfFeiqinqML81V/UzGNpCOrCaaeY5YejPBwl5yoeVMxK7E0cIgwFGCCNfGV+zHVdOYvgbPJD8F1zdlWVeT5+vJZXx7b4ja0otsXGDRVXqeBYVnmOQVbWZTduR++Zh9i051fLaYdWGeT9VEzgUt3IQ1vrNcMrfsvEIQ8RuE/To502Yka5aF+sRTU4YiubwmzeU1aT6T2M9zYxq/Gali02QN9w11cC+dxIdKrI2NcWlshK7IOIYiQYLluBQtCdioiiCge5pbE8rh1VwVg8nYpWTCq4ilnyWa2Uowv4dU7FJSsqN8joajRFDcIv7CPkyjAVs/dUp6zHD5i4uH+f93x7neupton8299lvI5ha+0r3Z8wohT+N4OV9paWlpB/Y99dRT52xmT1HofxJZStP93M9oqg0dt9yqRJIr2YzkTCzXJaCpKKeYb3ZdyWC2hAL0ZCTPjoR40WzluXQNBVcjqNisiY1yWXyE1ZEx/KoX8CEl2I5LOWMan6rg15XjtLZmpahO/ZZgoZuiiJKq20jRf7gYAtJBdfI4WgTT13jSsfSh0yVs7z3IW0tf40X/DWwN3cpNK1bOipl92WWXAXT09fX1zGjjFzDzraDfvEYN1GPmTu5MEwjCPp2QT2OyaDGSNRFAQFdOmHmkKIKqgM5ItkRUs7k8cpAbomksV7AtU83TqTqenazl8YlGDOFwSXSMS+MjrI2NEpp6chIsKSkWbcBGVwV+TcHQVKQeZ7juTQQKPcRGH6Zx+D4y4S7G41chFR8IFUeLoLp5/IUeSr4WXPXkC8MJASsXLiB3YDGJwsMsK/76FW9mzycqwnwGqP44yJMvjTqFQBD3G4QNldGcRbponXA8DRDQPI3qOIfb1RXJ2tgYa2NjOFKwIxvn6VSd95qsRxMuXZFxLouPcFlshJBmo5UtAFdCtuQgSw6qAL+mYvnaGAi+iQ5lJ9HMFgKFHkarN1IMLAQOj6X9xf2YRuNpze5k4DrWFL7Nvblrue5MbmCFWaUizGeAOIMCeACaotIYUYn6NQbSRUzHIaAfk94kBCGfylDmxMMdVUi6IhN0RSb4wwW72J2P8nSqjqdS9Xyrt5a7xFLWRMe4unqINdExDMXF0Mo/GhKKtkPOcsiZAqXqcqr8nTRMPETjyP1MRtYyEb/Cq4aiaDgyiFEaQJFFTKOBk82Piegi8qMGYV9lIfb5REWYz4izK/ET1DUWVgUYypjkSjZBQz2qYohfVZiO50IRkAilSYTSvLt5L3vyER6faOSJiXqenqwnqNhcGh/m6uohLgpPoJaXZdUAUxGUHIc+t5bB8G10ln5HTWYzvlI/I7Wvw9EiIBQcNYJmpVAck5K/BSmOz60MGhq/s9dznf4I2bO6IxVmg4owv0xoikpzzM9ozmQ8bxI6QqB19WTzyCdHCFgSyrAklOE9LXt4IRPnsYlGnkrV88h4M1Vaiauqh7i+pp8F/ny5D4qXBikNupWrGVca6cz9hqaB7zFUcxNWcKE3haWGy+PoXkr+VtwTLEr3YOlKniu0zVqWw/nC2aRAHruv/P5rDsdc3H02JXwrwvwyIhDUhXzoqmAwUyJcFmhRrhV2tuXxVCFZFZ1gVXSCD7Qm2TJZw2MTjTwwvICfDLexLJTiqlA314XT+BT30ELq2WCCpNFAe/rnNI/+iAOhq8hHV+PTVBwliOIW8RV7KPnbjpuPXuQbpc4+eO435TzmbFIgE4nEl07wnc/g5eq/9Vz6c5ok3AqzQdxvUB8yyJvOjLftU1wurxrh053buePiJ/j95j1M2gZ3Da/hAy9cyV0HltJXPFzMr6TF2V31FtLGQtpyjxEbf5SRTIGcaWPhVVfxFXtRXPOo62zUHuHNgVe8J/uoFMfy59MdP913zpoyvUTzAAAgAElEQVSKZp4jqoMGJcclV3IIGKeo+XMOxHWLmxt6eXN9L5uGDZ4sLObXY038YnQBa6KjvKHuAF2RCVxh0B29iZbckzQVnscvC+zm1WSEQ8inElYdfMXeozR0iioc9s9Kv6fLvq8Zs5rP3PEJ83QZWWeTAnmy72xMJBJ34yVmvDWZTJ5xunFFM88ZgrqQJxiuO7uBO0LAssAYf9y+g9svepK3N3azNx/lC3sv4U9e2sCjY404qPSFr6IvdDlV5h4S+V/hV1zyJYeRPBQsF1+x99BSORF/4NRlfyucKV8sm9lfxDO9z5iKZp5DNEWlNmwwki29bNeM6RZvberh5ob9PD7RwE+G2/iX3hXcO9TObY09XFW1FlfotGZ/S2f6Abpjr0OiMWFCyTKJcAAr0IZhHkSRMz9MOBPKmnMu85nPJgXyuH1lp9lU/PZDnKUwVzTzHBM2NM52yutc0BXJdTWD/OOyZ/h0xzZ8isPX96/gkzs38AvzWvZHXk3UOkB7+leAi09TKGGQyuYQhQG6zUacyr/PQ5S90qdIgTz2+HH7EonExnLSEXhj6LNKPqpo5jlGUxQifg17lk3tkyEEbIiPsi42yqbJWr4/sIi/717FQ7Fm/lcdrCg+TFv2EXrDr0ZTFRwRJJsd47fuq/jxxGX81Zz0en5wNimQ5ampE6VFfjuRSEw1fVbJRxVhngeEDW3Wx82nQykL9ZrYGD8bXsDdgx3clvkwdzVLrig+QlGtZjh4CaoicAiyQXmaV1c9QZZ3zGm/55qzSYE8SVrkOU1LQcXMnhf41LkwtE+MJiRvbjjA15Y/zdLgJO85+FGettbQnHuSaMnzXquKxl6zkUcKa+a4txWOpCLM8wBdVTm8lMz8oNYo8VeLt/LO5m7eN/wn7LMX0Jp5GM31osk69UGuDcz7uhKvKCrCPC/wosDmkSwDnul9S0Mvf9Sxi0+OfwJck+b0b0BK1vt2YojTZ5BVePmoCHOF07I+NsqtrSn+Mf1OaqxuwsW99Nr1OLLy7zOfuJAdYCpAf3//jDXoWjkKg2MMjuVwpUSZwQzA1EQO2yzgOLPzSHLZPKZ99sLXQo4BpYZHRxfSkHqGIXch/YpKbmhoxgvwHfHMZic07gLlQi4bdBVUymCc51zd19f3+Fx34nzhQtbMzwJXAwPA3IYqVThTVLyKxS/b0i4XAhesZq5Q4ZVGxYNRocIFQkWYK1S4QLiQx8wVKsw65WoiG/EWQ+xMJpN3JhKJ+KnWHJ8tKsJcocJZUs50uiOZTE5lQd2dSCQewktxnNVF4k5ExcyuUOHseRvw4BGfx4EPTa0r/nJTEeYKFc6NYwPUHzzhWS8DlampChXOkrKZ/eccFuA4nol9z9nU8DpX5r0wl2/YRrw1nbvnyoSpUGG+cz6Y2W8DtpSLgn9mrjtTocJ8Zc692WXX/viUK7+8nnMKz1z54VRl/7KGnrbp0tLS4gPWUwnnPB85FM7Z19f38lU7PM+ZU2EuFzT7El7Noy1lwe4sl1UhkUjcweEVAj7ImWnm9VQSLc53rgYqiRbTZE6FOZlMPpRIJDYdsWtKK0+xEQ5p6zvxtPV0y1sMANx33300NzfPQG8Ps337drq6uma0zal2L1q+FOkUca0C0soizSzSziGPK2srgal1nwUIgUAgkd6q6NL13vEO9+zrob29HRSB0LyF4hVfDEUPILQgQpzdiGs27kV/fz+33norlJ9hhekx52b2MdRwtCldfYT2TgGbOGItn9PgADQ3N9Pa2jqjnRweHp6RNqXrIK0sjpnGLYzRYvRShwKaAE0gQjqoNeBWYaWSaOEFoKhkd3yXwMIbUPwxUr/7POEV70WLdTD+208RXflh9KqljD/xF0Qu+kP06mVktt1OU2wlbQs7KA5txojHERpIZ6RcrUhB9dWihBpR/VUomv9lvxcnoTI8OgPmmzAfR9l7vei0J56E7du3Mzw8PIM98ti8efOZf0lKcAsIK4OwxlGsSUCCEEhhgNDofv5XuP5mpGIQ2vsPlBrejB1aQnTnp8i3vhc7cjGR3T9ipBjDDi8nWCgxfrAfdxQCbphU3zDuqCRQdJg42I8cNgntfwSxoJVd235LJPk58gs/jB3sJLj/DorNb8MJtKEUn0bqMRACV4shfQ1IIwbi9P8iZ3UvTsHQ0NCMtvdKYb4J8xjeFNQU4yc7cbp0dXXNuObYvHkza9eunda5Urq4pTROfhgn1490SoAKogloRmh+Uk9+Dn/rtezPNRLq+RciKz9EePm7mchcRP2ilRj1qyk1/B8aqpag+qth1aOHL7DysiO2Lz+8veqqw9urX8XOnTtJLOnAar2dhmg7bmGUybEYTYuW4ZRSjD/6BWqu+2e0aDvFgafw1UUQhkCLLESLLEDRDy82d7b3YrrMdOWS2WQaS7p+Gs9/AzO4fOuJmG/CfA9lM7rsvT4vyz9KKZFmGjs3hJM9iOuaCKHjltLgltBiixj60RsILNxIdM3/xDXTnpCrAepu/C5qpA2haFRf9XeH2vQ3XXrO/VI0P7661QCovji113/T218YI7b+z9BrLibf/RPSW/6J+jfeh50+SKH3YfzNV6DFOtBjHSh66Jz7caEwjSVd1wDry8c34gWYfJEZWL71RMzpPHPZsbUOeHsikVhTjprZW/7DP8hZVvafK6Rdwkr3Uux7nOLA77AzB7BzQ9hjL6H6q5h44s+ZfO6fEYpKpOt9+BdcixCC2o23E1p8CwB61dIzGrMe1wcpyy+X6QYEqYEaQotvBi2Ab/FbiG/8Dravnuz+X5PZ/m0sLU52/yNkd9+POZ5EOtZZ9+8C43TLs76dcnRYebh422x2Zq692ffgaeMj902ZHOdFpJeUErc06Qluts9bNL04iTn2IqHFNzO56SuYo9toaLma2PpPe2YyEFoy/ecqHQvpmuDaSNdCIo8pmi+88bjgsHf7cAcBCeYEhcwQtitJmw7bR6HR71KwHe7dq3BFs4IQgrteFNy2tJ74SIpv7n4Hv9++kabeIo1b/xm/4cNe8wXqIvtpaL8M3R899xt4Dqxbt25Wl3TdtGnTuS7pOoYX4jnF1BDynJdvPRHzzcw+b5DSxS2MYqb24pZSICWloS0E228gt+fHZF/8NwKt1xFd/dHy1I84ZOKevFEH18oh7SLglveB0AMoehihhxB6CEU1QNERig6KhhAqCAVXCgq2y57RHD5NIV20+MaTPVzVXsWBCYNv/XuR929opimk8ckH9vHJyxvojKt8Y/sg0aCPFXGH5LhFJjdJBJe+tB9bDTNasPjYwU/zhUstGvOSiac/S3/XR1h9+e/N+n0+z5la0fHL5bH1FF9MJpNfLlumX2IGlqaBijCfMdJ1EKVRin2PI60cTmEMLdJGafBpJp/+AlqwnvDStxFafKs3j+uLnbgdKZFOEWnlAdcTWiePYkRQogtR9BBC8yNUP0I5uuKs5bhMFix+/tIwzVE/miL41E93cOOyehZWBfjA3dt434ZWLm2r4gdb+6kLGTQYGpcurGJRXRUXN0X4xq0h1jTHaIz6uGn1CkK6iqoKbrgCXNfFtYtctbpALp9icryXr1wXJR4QvHjgIBtKB+kdGSG2f+fLcMdPTllzztslXZPJ5JbyInF345nbm2Zq+dYTURHmaSJdBzs/hD2xCzW3C2QX0rUZ/dX7iKz6COHEO6jZeAdG7cqTLkIu7SKulQNcpATFF0GLLvQEWA9hj0Yxqpd7prJQKNmQL9j899YDNER0aoM+3vuDrdyQqOO6xTX8wX9v5Z2XNPO21c1kSjY+TWFDW5y/f/1yru2sYUldiLevasZvqDz/3BY+fcu6Q325fGH1Cft4mABQBTTDohUUizkGBnvwKy79jd+iKWDT27OF6upLZuoWn488BNyBp3lPtKQrR0Qzfhq4o3zelFCf9fKtJ6IizKdBSomTH8aaSCLtPEKPEuj9DpP5hVRd9jnil34WX/OVCNXAV7fq0HekayLtIrKYwrEySKeIEDpCD6OoOqgG0sxgmRks28F24L+3DrF1ZJiuuMltvzS4usnmfyx3+dQDPl7brvKJ9RGWVwkuqnK5ukXlofdfwrLGauojft50UdOhPnc1HW8NnOwHZrr4/SE62i+iuqqBbdt+zdCz36DF3Q2v+uk5tXs+c7olXcvbU8eeTSaT95Rnac55+dYTURHmU+AUJzDHX0KWUpgj2yj2PU788s/jhJdh1HhxLIGO1yPdEm5xAtfK4RTHcQqj4JRAKChGHCVYixZqxFH82JbJPXsUpBCsDXTzD48P0BIS3LZE0N2fR+b3s2L9RWzs8HNFfYmLaiwef4fOwniYgOZyz+sBOQS5YS7zSUgpWPkIrr8a1RdHGOFDY/TZIBarpW1hF3sO3MiEvJHTeAEueKaxpOtrjjmWYobGyMdSEeYT4Fo5rIk9mBN7KPU/RjjxdlwzjTW5D2mmKdXfQKB1AWZqD7KYwnWKSDvvhUtrftRIG9JXw7Zn74bcEzSIYcYnBqliDENY3JX/d1w9ylWN/8kXA094vq4kfCUGSKgy/o5v3fRqsrt+SPqhrxIH0qqfbKAWNVBLoG0joaW34dpF7PR+NKHh2P1Y6R4EEoSGGmxADdZ543YtMKP3p7lxES1VMQL+EweSVJgbzgthLnsCx/Em22ctkEQ6FlZ6P8WBp1G0AG5xnMzz30SLtuNvvwGj6QqsdA9q5iVKY2mkbWIXhtnT34eb7aeRA4zkXT5a+jpfuw70vp+jySIy1kAmsJRRo441rbX8V0sDDVUxZOZjuOa7ESggFHr299Le3oEW9hJDfI0biF/6l7ilFE5hFLcwhlMcPdRfa/wlxn79YQC0yEL0moswaleg13QhFB0714+QEmFEEGXLYSaCPnTDx2D3o2wwHyB3zW/Pub0KM8O8F+ayIG8qewaPTImcMaR0sbP9mKPbkHaJ8Yc/SnDRmwmv/BDV1/0fhOrDGtnGjnHYNZqnU3G442mDDxf/ipAo0AgMO9WYtYvpF+1cXWtwUUOEsZu+R13YT1BXaD/RhWMdR310RlSMmuWHPuvRdvToCb8JgBZrp+rKv8We3Ic58RKlgd9R6HkAf8s1VF/zZYTjI7/vZ/jq16AW+ikefAzhi6JF29ECdQj17Fe+Cy68gd+MrWDd6U+t8DIx58J8uuIEeBPzPyyffjoX7BkhpYudG8ZO7Sb11F8j7SLV1/4TLPsQT+WbWb3/KX6YdOgeSPLxpicJT3bTZgc50PYZaus7eGboZm7oWkwu0kVLtJm6kEYn8Ppy+2Hf7N5e1Rcn0Hb9EX+PxMkPgWMCYI69yOSz3sxHMNBOXn0dvqbLcc0MFgIt0ubFXRvhM7724pZWqvRKJNh84nwoTpDCE+IZKyruOiZuYZzxzV8jtetHRK/7JtudFTx9MMOrN/2Cew+28c09TTy3+DPcWjyACEpkPorWuJ7musuIiyhvuagO+GMAas/w+tJ1gKmcY/fwGutuyQsYEYJDucpCmXausRACLdR46LNRt5r61/+AwsHfMrH752Se/yaZ579JaOlbiV7yCezsQez0PtRQI3qs86Rz4ifipU3fZ3X2noqZPY+Y63DO6RQn+BKHQ+LOKotKujbDB7YSiDSwY9MPCO74e4YTn+PgoEl+pJ6OZ+6llBuhVQ5RM/A071r/FV67yk90YAVG9DUYDevQYx0I18F1Lfbu2oE54QNpg3TKwulpeuG6SFymUhsVYSA1HUXREYqGVHQU1Y/QfKCoCKEihBdKidBAqN73XQukF76JlEeHcEq8YgSKN8UlFOO4wBIoC3d0IZEV7+agWMeStiqKvb9Gi7YjFBUn209+3wMEF70RJzfoCXV8ybQ0dTxxC4/3LuEVPcs8z5hzM/sYjitOUF7u44OJRKIab4L+jPjR9z9Bfdhmg/VzfmvchFSDjBRX05ArsLb0M3yhAfThp7yTw9VoVWtoidosdUeRkZs9TekUscZ3esoUiSgN4+QML/yyLMDggGqgaCGEHkSoPoSqAYonhIoKZcGdqhKi6DEUfwxFD6FoIZyogr/5+HRCKd1yXLaNdExwLVynhGtmwMp6KZbS9jzZCC9qTAsglKMfrxZqJLz8XYc+W5Pd5Lt/Qn7P/fhbryW46Gac3CBarNP78VKNk97XtpoqAnbjSY9XePmZb8J8Qs4l33Nt8Wc0+bw453Xmzxlyaoi2vhPbLvHIxAImnGVsMxezpbSUETdObU+ab43/BxElxyNjVeRMhxp1kmolTZWSIa5kaNRSjO0rAQ4KR2cmuVLgCB1FNbCVIKNOnIISo6TGsdQ4Ra0B6a+nszqOqgj6sqApLgHVxa/rbH50L8JXhVSDoE4ne0oB4uBaCNf74RHWEIqVAmlSDu4Gxc/OnceGXy5BLP8HjNGHkP0PUzzwCFaki2LTLUijFifYiTSqy2b/0Wx74htc7vyM3DW/fcUXJyg7aU9YBvp0+c4zyXwT5hkvTvCsdiML7R46lH1ElQIN6hgNI18HIBHQcRWDd8pHUKWJmEpuyHtvv+8D16eQFxGyxMgSZYKFuL4ltNTVsSejsHVEoeiqlFwVIS00abI8brO6xiafT3NgeIxatZ8F6g5iSs5rOA1yWOD46ngh3cZzxcU8by5mu7WYSddPQCvww9eOEw0EuH1nkJcmBCG/j4hPJ+LXiPpU3tLVzIrGMPsnCvSMF6iNGNSEDOpCBtVBHZ+morglhOPFj+987jGWtDeCFJ7lcFRgyWW41ifJ7b6P/J77aL34SoQWwM4PoUdV9Orlx6VljuTeyeO9K7kEXunFCe7AGw4eN8tyunznmeakwpxIJD4ATMcjIoCJZDJ51wz0Z8aLE9z6nq/T2tpKKl/iqRc283zyOcYGk5QKEwREiaAoEfZp1EWCLKyOsKixkXC0FtVfgxKoQTFix41Hd+7cSefy5XQCrz3FtR1X0m5KcqZL3nIZLmSxswcxCgdoEX04mV7WyyTX+g67DSbUZrrVlSx01mGrF2G4OaySpD+XJWur5G1BwZY0RvzkbYd7nx/ge8/1HXVdAbxlZRPvXd9K/2SRu54eRrPbWDZcRVPIpcWXY3GojxU1CrruR/eF0bUAkRXvJrzsnQhFxbVyjD/8MQLtNxFccgu++rVoofpD15gvZva6desePdH+TZs2XVs+/jU4YaDaJzZt2rR13bp17wXee7Lvn45kMvmhskCfiKPyncsO31njpMKcTCa/PZsXhqOKE6QSicTUHzxVnGANMxi3Gg/6uHHDFdy44QoAsiWL3+7q48HkAD/unuCFbhtrj1ffckkcLm82uK5D59ULJXVnGWehKoIqv6DKP+WNri6/Vh46px5wzSzW+EscfOlRmpQhqocfQ774CwC+ULUU3/or8C94FWqo2Yv3RkWGbAjoXNywlPesW8BozmQkV2I0ZzGaM1neECbm1+iZkIwXTEbSJk8NDuKWRwUr6sN88fp6DGuYjz84RGPApTWq0FHlZ1mNzvqaIoH4MrI7/o3iwUeJrPwIwc43oFctQgiFbU9+p+LNPj2ny3eeUeb98jRnS0tLSzuw76mnnppWDbC8afHEnj4eTA7zm+5xtg1blMq1ITuiLpc3wrWtKte1amQHulmSWIZAgFDLU0hl55ZQT+hZng47d+5k+fLlSNfGGn+J0uAzlAafxRx9HqSLGmrCv+Ba/AuuRY22gWuh6CG0WAdasP6kDivHlTz17GaWX7yS/nSRnvE8RdulLR5gIFPkSw/vpn8yz0DGwiqPNMK65P43BwilniHe/S0MaxTf4rdRvfqj+BvX8MsnHuBg7wtcsvymWTGzL7vsMoCOvr6+nhltfBYoa+a7jx0zJxKJL1FOsCh/3ntk3PZMM10z+0RR+/KI/TNlZs8ZQUPnNSvaec2KdgBM2+XJnjEeSg7xyN5R7t+X4/u7XMCiMbCAaw4INrb72dgRoCnolquBWOAWcMvTSYfuzjSnkqYQioZRezFG7cVELv5DnOIExb7HKB58lNzuu8kl/ws10kaw8w0E2q7HGn0RS+xEi5aDQI4J2VQVgV8TVIcMqkMGFzcdXSHkpmUNFCyHbMlmx+AEz/X0MzI+hG6NYUUSXHvwH/mr6Le5cfcP2CSWs7RYorUqOi/M7HnOKfOdZ5o5NbPnM4amcO3iOq5dXAeA7bg8cyDFg8kRfrx1Hz/bU+CHO7IAtFcFuLKjmtcmmrlhaT11YeNQiR9cC+mYJ5lK8vCmsU6eDKH6qwgtehOhRW/CNbMUDz5KvvunXhDIttvxNV5KcNHNSAH2ZDdqqNkrvuebXlkfQ1MwNIVYQKclHuA1y5opWg7jkxPs79/Lh1ce4Ad9H+ebQ2+mZ7Cdj4yPsWriP1lv/7JiZp+a0+Y7zyRn5M1OJBLRZDKZTiQSq5PJ5NZEIhHFmwvumZ3uzR80VeGK9mquaK/mDbVZVq2+hGd6U/wyOcxDu0e5d9sA39vShwA6a4Jc01nDaxN1vGZJHdWhGMfqYemYSLuAY+VwSxO4hTGENYFTGAGhomhBxAkK+ylGmGDnGwh2vgE700u++2fk9z3AxOOfQQ0vILT0bfgXXIOT7UcNN6LHFk1bqI/Er6s019bSXFtLV+dy3n5wB8/sKfCNnQJr7/3URLt5JPpHbDi723nBcIQ3uzORSHSW4yL2JpPJRSfLd54tzmjMnEgkPoU3qO/AK4OyBS/8ct6Z2Gc6Zj4TTlQr2rRdnugZ41fJUX69e4TnB9JYjkQRsLQuxDWdNdyQqOP6JXVE/SdOcNj87FOs7lqKW0xh5wdxS+lyRKfPq/+lnPi3V7o2xQOPkE3+F9bYDoQeIZR4O8H2mxCqhhpewAv70qzdcOU5/d2jowfZsfNJdm7/Ba81/5O9NW+nasUnXvFj5vnCGWnmZDL5DwCJROISvEitjcC+WejXeYehKVy3uI7rFtcByylaDr/pHuOXLw3zyJ4xvvPMAe58qhdVCJY3hLl2kae5X72oluBUQoaio/qrUf3V6PFOzzw3/1979x4bVZUHcPx759V26PTdabHFQqm9BVqo5bEsu+ACZQXMikbETTYxmxhwE7OJughG3Oi6ieiauIluWCGrMetqXHTVaFZQUbQCK09tqNIptIAFCoVO37TTee0fc0sKtDAd5nn7+/zD9La5PUzvb865557z+3XhudiKt7cFv9eNYjBhMNsu2/GkGEykFC0hpWgJAxfq6Pnhn/TU/YNex9uMK11FysSlmDpPMdCehzmtKOTdUjk5hVSU3waei7hqtzLpwlY6eDgM754Ih1AXjXQ4HI5O4D/hbIyeJJuN3K7auV0NPJvtcbn54ugFPm04z85jbWzac4K/7T6B2ahQnm9jUUkOE5V+KjxeLKbAoFwxWjBqCQn8WWX4B7rxXDyPt7sZr8uFYjBjsNgu67EtOeVkLfgL7vYGuutepef71+g9+i7mnOW42yfh7W7GnD0Vo9UeUjaSzMw88gqmUVdbSQXhXfklbkzQwaz1xouBHdr98iKHw/FF5JqmL6lJZu4sH8+d5YFcXe29A3x29DyfNVygpqmNF79qwg+s372dmYUZLC2zs6zMzvTxNgyGwM4pJSkdS1I6/ozJWsWMs3h6msHrCaTgHTKLbc4sJWv+87jbG+j89mVSzrxN2xe7sE1/EK+nD/O4fMzZU0LKQlJYqLLXb70qe7eIrdH0zLMIDKmfUFV1EoGZuYgHs7YSrJrAaoth178mosxxFlZVFrCqsgCA1u5+Nm3fT31/CjVNTr7eVs+GbfXkpVqYNzGL5VMCwX1TegqKolwKbHNGCb5+J+7OE3j7WkExXbZqzZxZSvbClzi2byvpFz6gY88fScqfS+qMB/H2tWHOmYpp3E2j6qWtKTayvWcwGz0ReW9EaK4ZzKqqLiKQ5aMLaNR64mgPrVcRGA00abOCugjmK9ltyfyq2MrT2mTS92e7+O+RVj6pb+XThvO8X3cWhcBk2sKSHJaV2VlcErjfNlpzMVpz8bl78XSfxtN9Enw+DEkZKEYziqLgSZtO7ux76D32Ht21r9C243fYpv6WFG8fPtt5LFlTAtsyg6AoCluVB2hpc7Ixgu+JGJ3r9cyPA35t+2GTqqoZDofjvdH+kutlE7nWTpLBHVNaDx2WMh6JYFp+GtPy01i3sIQBt5ea4218fCTwGGzLNyd55X8nSTYZmFmYzvIyO3eVj6csLxVLVinm9Il4elrwdDbhc7kuJR1QDCZSS1eRUvgLOg++SPfhzfQ1f05a1SP4+tux2CsxJmcG1T5LchpnvTLMjifXC+Z7tYkuVFVdDCxRVfUJArub3gnmkVSQ2USC2UmyBlgfxM/pjsVspLrUTnVpYDKttaef7fXn2Xakla+a2tiw3cGG7Q4mZllZXJLNimn5LCopIMVWeClxv+J24ve6UIxJGK12suY/R9+pGjr3bcT55cPYKlbj81zEkl2OKe3m6w67l/X/ncfth+lgTzTegrg2mi2Q2r/RL+k6GMja68+1RqCqajoEl8styGwig731UB2Db472vS0EevKELPMaTvbUZO6fNYH7Z03A6/Vx8FQnH9SdZZujldf3n+LVfc2kWozMLcpk+ZQ8VkytwmPtx+914XN1XRp+pxQuwJJTTse+Z+mu3YSrZS+2yodIclcGht3XWHKalpaHr6cuiv/r+BTCFsj1RKika0iPprQg/zzE33lVNhHtnO8O98NDevYO4AAR3A+aiIxGA3OKMplTlMmzd0yhpauPD+vO8dEP59hzwsmOoxdY+xEUp5m4tyqXe1UTJa4fMSg+DEmBZ9pZ81/gYtOHdB38K+01a0mf+Rh+Tx9JuTNG3LyRm2rB0CvD7ITYAhkvtN455J0mhw8fprW1NYwtCgh3do1wnndWEsy61YSrIoeDrQPsbO5j79kBNu5sYuNOKLIZWGB3cUduI5PTwGhOA6UUQ8mTWE+8jHP3Blz5KxjIqcaTNnXYjCfOlt0mjxAAAAWBSURBVHpsijdsbR5qNJlGtP3Mrx84cOD1cL4Odj/zdYy0BTK6JV0jmJwg7NlErqWioiIqyznj9bzzgN8De/cfwJdbzFuHTvNJQytvNPp4o3ECpZkGVhQN8JuKNIpLVYzlc+nY+2eUU++TYe4gffKjJBfMuCrJ35u7JlHkOwGM7UwjIYpuSdcI7poKezYRcX0mg8LMiVn8dGIWLo+X2tOd/OvQabbVt/LCdz5erO1lbl4Xd5dauafyaVIz36bn8Ga8F8+T8ZMnsRZVX7ZhQzGa8ftiP8we2oNG4vUNumoLZCRLugaXkPkGDMkmcp+qqlXakGIwm8gaIryTRFwtyWRkTlEWL91dQe0fbuPjB+awcvpN1HcmsbZmgOmvtvHw8eU4y57E3d6As2YtPY0f4HN1XTqHs8+vZQMV17ADbSZ7cAukqqrVWicG0SrpGmRygsHvjTjM1ia23r3i2OBUvC4XgCQSq8XE0il5/FK1c6F3gHe+PcnWQ438u97NVkclqyc8xaP9z+GsWQd+H6m3rMRgSeW7lDv48NyMMb9oJMQtkNEt6TrWkxOMNQaDgt2WxEMLSrl/TjG7Go7z2u4jvHZcpcb0DG/Z/4Tny3V4PH4ypv2aUsMxypJrgbti3fSYcjgcVz1duaKk61UlX5GSriJabMkmlk2/hXklE9hTV8umb9zc0/IMb+Y+hbtmPQOGFKoGtlFsbaA91o0Vl0T8nlkkrnRrMktnz+GVVbN5/GfZPNb7FC63m+YvN5BszcIXZA0sER3y1xDXpCgKBfk3c9/CJTy9qJAtPEKOtxnDxR8ZfipFxIoMs0VQUsdlMPfWhZgUL9v23M5S4yexbpK4gvTMImgpVhsVU+eRc8sy/IBPeua4IsEsRiUjI4+pJRV8PzAZlz+0XGIiMiSYxahNKJyCX1FIUQZi3RQxhASzGDWbLROnLwOvXD5xRf4aYtQMBiOFplaMgyVwRVyQYBYhyTB0Xf+HRFRJMIuQ9Cm2WDdBXEGCWYSk3ytLFOKNBLMISboyYkJVESMSzCIkLd48mf6KMxLMIiQn/TfjldXAcUWCWYSk2NCIGSlPE0/0/NFqBDhz5kzYT3zu3LmIJJ2L1Hkjce7uLhctRgOdEWjzkL/ZyIm7xVVGVWw9kRQUFPwc+DrW7RA3ZP7p06d3xboRiULPPfN+YD7QAnhj3BYxOkZgPMPnnBYj0G3PLMRYIxNgQuiEBLMQOiHBLIROSDALoRMSzELohASzEDohwSyETuh50YgII1VViwGnVpJ0sLpnB4GSpVsdDkfHcMdi1uAxSJfBnCgXlVZBcD2Bkrc4HI4d8dh2rRzp8wQqFh7SArtYK4qGqqqbVVV9/spjaHW4RXTobpg95ELboZWODVsx6wioBo4DS7RAjsu2OxyOHcCBIYcGP3AGVY9wTESRHnvmRLqo1mv1qwclStuzgaYhX2eNcExEkR6DOZEuqmJtCFtFoCB9IrVdxBk9BnPCGHJ/6QQ2A4di26KgtXH5B41zhGMiinR3z0yCXFSqqq5UVXWd9mUTgQmvhGg7gVHEZABVVTMIfAgNd0xEkR6DOVEuqg4CbYVAIMdtQGgz7LOA+1RVrXI4HE1Ao3aLsAZYPdyx2LV4bNLlfmZVVdcQ6O2qgC3x8HhnOFo7ncBsYKP2rDYh2i7ijy6DWYixSI/DbCHGJAlmIXRCglkInZBgFkInJJiF0AkJZiF0QoI5gamq+k6s2yDihwRzgtK2SwpxiQRzAlJVtYrAXucm7bUQsgIsUWlD7NWy3FMMkp45cRVLIIuhJJgTkDa0PjDktRASzAlqcLvhSofDERfbJEXsyT2zEDohPbMQOiHBLIROSDALoRMSzELohASzEDohwSyETkgwC6ET/wdXzpu/YFXSaAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 144x180 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax=pl.subplots(2,1,figsize=(2.0,2.5))#ensemble statistics\n",
    "dfs=(df_std_all,df_norm_all)\n",
    "leg_str=(False,'full')\n",
    "tmax=125 if d==0 else num_steps\n",
    "xticklabels=[0,25,50,75,100,125]\n",
    "yticklabels=[1e-2,1e0,1e2]\n",
    "\n",
    "for pit,df in enumerate(dfs):\n",
    "    ax[pit].set(yscale=\"log\")\n",
    "    sns.lineplot(ax=ax[pit],x='iter',y='x',data=df,style=r'$\\alpha$',hue=r'$\\beta$',ci='sd',legend=leg_str[pit],palette=sns.color_palette(\"colorblind\", 2))\n",
    "    ax[pit].set_ylim(2.3e-5,2.3e2)\n",
    "    ax[pit].set_xlim(0,tmax)\n",
    "    ax[pit].set_xticks(xticklabels)\n",
    "    ax[pit].set_yticks(yticklabels)\n",
    "ax[0].set_xticklabels([])\n",
    "ax[0].set_xlabel('')\n",
    "ax[0].set_ylabel(r'$\\sigma_{h^i_t}$')\n",
    "ax[pit].set_ylabel(r'$\\left|\\left|h_t\\right|\\right|$')\n",
    "ax[pit].set_xticklabels([r'$0$' ,'',r'$50$','',r'$100$',''])\n",
    "ax[pit].set_xlabel(r'$t$');\n",
    "ax[pit].legend(bbox_to_anchor=(1.0,2.1),frameon=False)\n",
    "# fig.savefig('fig1_theory_panelb_d0.pdf',format='pdf',dpi=300,bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot Fisher memory curves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_Jk(alpha,beta,rho,n_steps,num_steps,N):\n",
    "    M = rho * np.eye(N)\n",
    "    for i in range(N-1):\n",
    "        M[i+1,i] = alpha\n",
    "    for i in range(N-2):\n",
    "        for j in range(i+1):\n",
    "            M[i+2,j] = beta\n",
    "    M_k = M.copy()\n",
    "    e_k = np.zeros((4*N+1,N))\n",
    "    e_k[0,0] = 1\n",
    "    C_n = np.eye(N)\n",
    "    for t in range(4*N): \n",
    "        e_k[t+1,:] = M.dot(e_k[t,:])\n",
    "    M_k = M.copy()\n",
    "    for t in range(n_steps):\n",
    "        C_n = M.dot(C_n).dot(M.T) + np.eye(N)\n",
    "    J_tot=0\n",
    "    J = np.zeros(num_steps)\n",
    "    for i in range(num_steps):\n",
    "        J[i] = np.sum(e_k[i,:] * np.linalg.solve(C_n,e_k[i,:]))\n",
    "        J_tot += J[i]\n",
    "    return J"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "J_tot: 20.32295558140331, alpha: 1.05, beta: 0.005, d: 0.2\n",
      "J_tot: 12.086536304234857, alpha: 1.05, beta: 0.005, d: 0\n",
      "J_tot: 20.328096176553625, alpha: 1.05, beta: 0, d: 0.2\n",
      "J_tot: 12.08727212122481, alpha: 1.05, beta: 0, d: 0\n",
      "J_tot: 16.1637115854401, alpha: 1.0, beta: 0.005, d: 0.2\n",
      "J_tot: 5.301425802978977, alpha: 1.0, beta: 0.005, d: 0\n",
      "J_tot: 17.132327989761464, alpha: 1.0, beta: 0, d: 0.2\n",
      "J_tot: 5.18737751763962, alpha: 1.0, beta: 0, d: 0\n",
      "J_tot: 12.094319320360942, alpha: 0.95, beta: 0.005, d: 0.2\n",
      "J_tot: 3.175767594604019, alpha: 0.95, beta: 0.005, d: 0\n",
      "J_tot: 11.939270436426913, alpha: 0.95, beta: 0, d: 0.2\n",
      "J_tot: 3.0327020748310467, alpha: 0.95, beta: 0, d: 0\n"
     ]
    }
   ],
   "source": [
    "dvec=[0,0.2]\n",
    "n_steps=int(1e4)\n",
    "J_store=np.empty((len(avec),len(bvec),len(dvec)),dtype=object)\n",
    "for ait,alpha in enumerate(avec[::-1]):\n",
    "    for bit,beta in enumerate(bvec[::-1]):\n",
    "        for dit,d in enumerate(dvec[::-1]):\n",
    "            J_store[ait,bit,dit] = compute_Jk(alpha,beta,d,n_steps,num_steps,N)\n",
    "            print(f\"J_tot: {np.sum(J_store[ait,bit,dit])}, alpha: {alpha}, beta: {beta}, d: {d}\" )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANEAAAC9CAYAAAAtDJjOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXd4VNXWh98zPb13EgghHFpAegARlKIoNhTRa+/t2vVib3yI2LuCiliwgGABRaVJkRqKhHYoCRBSSe/Tzvn+mMkQIEDIDJkknPd5eMjpa8pv9t5rr72WoCgKKioqTUfjbQNUVFo7qohUVNxEFZGKipuoIlJRcRNVRCoqbqLztgGNQRTFu4BiIEOSpE3etkdFpT4tviVyCihNkqQfgbu9bY+KyrF4vSUSRbEjUCxJUqlz+2qgFOgIzAaSnP8DhHrFSBWVk+DVlkgUxZHAHByCqRNUR0mSFkuSNB2Y6jxVFY9Ki8WrIpIkaTGQVm9XXStUx0hgHxDs3C5uJtNUVBqN17tzxxAGZNTbDpUkabooineJohgKTGvsjeLi4oxAfyAXsHvWTJU2hhaIATZkZ2ebT/filiaiBnF27U6X/sBKT9ui0qYZCqw63YtamoiKOHr84073LRdg3rx5xMbGumVUc5Oenk5KSoq3zWgSrdH2nJwcxo0bB87vzOnS0kTkcmOLohgMuDMnZAeIjY0lPj7eA6Y1HwUFBa3O5jpas+00sdvvVRE53dn9gFJRFJEkaZMoivucXrs+wJ3etE9FpTF4VUTOCdQfj9lXN/5Z3PwWqaicPi0+YkFFpaWjikhFxU1UEamouEmbF1FxaZG3TVBp47R5Ee3ct8PbJqi0cdq8iCqqK7xtgkobp82LKK9EjVlVObO0eRFl5DYpkkNFpdG0eRGVmW3eNkGljdPmRRQVGultE1TaOG1eRLLWx9smqLRx2ryIducWeNsElTZOmxfR5tqO5JXXetsMlTZMmxdRleLLpMW7vW2GShumzYsIYOaGLLU1UjljnBUissuK2hqpnDHOChFZ7IraGqmcMc4KEYHaGqmcOc4aEVnsCmv2l3jbDJU2SEvL9uNxBBQ+GJjPfePv8LYpKm2UNt8SKQj8uz/L22aotGHavIgAwgODvG2CShvmrBBRXlm5t01QacOcFSI6VHbaOcpVVBpNm3csAFjRk/G2AUFwbGt8o2h/lzpOUvEMZ0VLFKipcgkIQK7O954xKm2Os0JExXIwduWseKkqXqDNf7MEIMsWyQslt3rbFJU2SpsXkR4r5YofpbK/t01RaaO0eRH5aGqpln0YatrmbVNU2ihtXkSBmios6NlsTsKiOJyRGt8oL1ul0pZo8yKK8rGjICBZEhhd/BU+t5Wr7m0Vj9LmRRTvHAr5a2vJqjGoyyFUPE6bF1HH0LqUWTpAUBfnqXicNi+iHvFxAFjRADI2u6y2Rioepc2LKCoyAYBSOYA+hl1YZTVxiYpnafMi+mZrBRrsFMmB1MgmQF0qruJZ2nwA6pbcCgIEC2WyPzWyAVCw2FGXiqt4jDbfEi24YyAxxirs6PATzLzcPZP3BpczzG8nsix72zyVNkCbb4kAkvzN7KqFKG0hf+2tokprJNfix+OHc4iLaudt81RaOW2+JQLoG+Nwc0doK+gRHci65+9g46MjiItqh8ViZs8B6bhrbFW55MwZga0qr7nNVWllnBUi6hHvCPOp1kXw0QMPodPpiHDWLXr3p+9JeX8b6bvTj7qmdO0rmLP/oXTd5Ga3V6V1cVaIKDkuHoAd1WGs2riKGT9/SYdnv0XK2M41Q8/jqV5ldO/UHYBnZ37KIx+/R+WOLwGZom2zsFaefsnK3PJahn/0T6Nd6e62fKf7PBXPcVaIqEuHLmiQqVR8eXDuBu5aFUQ/vwNY7Xbaxybywo23odE43oq8ChtSQQWK4nA63Jj3JBe9/qXrXo9O+5jXv5vp2v74p1nMX7HQtf3r8t/YsG0dkxbtZlVmMfd/v4K99bqLUsZ28g/nuLYP5h2itKLM1fLtWz6FWqsdAEVRqCzLx1xR4Nq22WXssuK63i4rWGyy63kv/bWbGuf1J0LtqnqWs0JERqORWH2p429q8BdqKVZC6ZbU/bhzP7ntCj4KnQqyBYDLfVdwqWmp6wu3sQB2F9W4zn9lnY1vN2S4tu9eWMzUPzcxc0MWsgI/7a7luXnLXMf7f7ydl+b84tru/MY6Xvj6S1fL13XpSGatWQqA1VJL4KT1TP75TwAOrZ6KYeJvvL0ig5qs5eyeNQLjxAV88E8mM9cfQFZg2toDrMwowlq2n5l//cE5by7HapdRZDvvrszglb+2UfDbfzBn/8OHc75gR96RTEhyPXGqNJ6zQkQAPQIr0WLHip7/JJayqjyByV99wt0ffMS+A7tc55WufcXVCgH8x/8vLvdb7RobLX/uXj79772u49uevZKP7rzetf37LT0xhvdAVhxfSJ1GQA7s7Dr+/gVabhjS27X95qBqxgTtcD3zueAvGObjWPuk1B7m+W77GT1ohON47j882XkPg9uHIJuLsdnMPD4kin8yi5Gd12sEmLc1l/LN71O06UM6hPpgK89k/weBmOQKXlsqUZW9BpBZuK+KR752CFpRFJKe/47x730DwNaDuVzwznzmbjnoOG4zY7cfaeHU7uMRzgoXN0DXcD1/FgnstcQxmCIuisjm5W1xDDJsZfnWLSS17wKAOW+tqxVyIVsw565t8L5B/oFHbUfHivy0cwkWu0NENhkW7K0hr7yW6EATN4+56qjz775oJIdm3I/ifOYt/gtQLEZsVXkYgxJ48bYHXOfGXz2fV1xbV9Kt05U8VF5L0itLsMiO30NZga83HeLpe2/jge52Hg3vgb0yG1O7YdwQt58rzl9I9S4NKDIvh3zGQWNf4EYUayU3BS2lZ1/HMvrAouWk51jpH1SAosQz5r2f0JTt4tf/3c3mIh/Gz1xDVrmNSYt38+G4nqf5abQtzh4RRYegSBrKCWBn7h7mPvkYY976lvv7i0y4aILrvLjrN7j1nEmLdrtaoTrqwowa+rId2/IBoMiUrptM+AXvN/l5U9PsrufpAuKJGfc7tqpcDi+8ERQbABGaYiKVVdiq8tD5RfPSE5+57hHf9UIOPixhDEumxmqnV1woKQlB6PxiCKku52C5lbqo+CFRVnI2fMg9Yy4EOlFltuFnPPLVslXlUvD7DURePAudX/QpX1Nro1V050RRHCmK4lR37tG9fXsADFiIDg4g0M+fFU/f4RLQYx+9ye2TH6cgd89R1+WW1zLs/WVsmXVpowbiaw+UuFqhOk5WkaKhlk9QrCds+dx5XkOCVWR7g258rSkEn9hUNMYgfA06pk4YzQ1XPwLAW8szMGgcOcjsssJjiwr5ovg8dIGJZJRZCXhmIa/MmuG4v6JQsqZtTxe0lpYoDRjvzg16dOoG5JNkKmZBfgzVVaX4+gUDUFFdzfcHw8m1dWbJuyu4vN3P7C60srk6lgI5DAGFNwsTSH77Tj62PUDui6NP+JxNjw47LbuObflyy2sZ+9FiPmcaUc4W4mSczvNO1FWtPbSSnDkjTtpSlNZYCfbRk1tey8y0bCxOLVrsCmW1Vn5/4Ar+b0seYdZCbk3MI7XyW+y2G5nw5Wrm7RzN8uh5ZKRtIa1oHe+PH9hom1sDzSIiURQ7AsWSJJU6t68GSoGOwOy6/WeSAF9/2umLCTVY2Vnrz0Pvv86nT052HvNl36TrmbV0MV9uhI/3J2NDy0BDOgWWUBQEvq8axX0Bc+htW8qS5aUEBwSSnJhCYEgcimxDka0IWhNC/SyRTWDSot1sKZR5syaBtxvZpWssDXVVc8trufr9H3jbNBFDvedZbDIVZhtdX1vGbQMSeG3ZXn66pT+frTuAxX50a1ZrUxj24WoqLQ7Hw+wbL2Hk19Gsza6Aoo1AMtssSdxf/AQUFHCv6X7kflPYeKiMm/vHe+z1eYszLiJRFEcCU4E7gU1OQXWUJOk15/FpwN1n2g6AZP8q8iptJOmyWFTSnt1v+aDXOD54jW8Ut9+Vxe0XwoJ1q7lv3nb+tSSjxY4dHXY0vF9xLQB/zAewMvfCVVw5agJzV67gnt9y8BHMGAUbRo0No8bOtHEi/ftcyB9rlvHx0jT0GgWdRkErgFYQePnqUXRI7M3idcuYt3Y9tbKer7OTUNDwQ9VIbt3yP84bmMfqnRIrtm5CwHmtBrTAPVdej29gFKs3LmOrtAmNABoUNIKMFpnrr7gPnSmQtLTfyNi3EbsisNsWT7SulC1Vcbx/y1U8M38Ta0pCuFwzhYtXrqdj0TwmjhnIlbP2sbfERmGVhZUZRQC8uXwfqzKLG3xv6wQEsOaAoys5e9NeXjJOYR4z2GuLcx3fpj+PjRsP8dbyvfTY+zg9Rj6NMaIX4Ogezt+Rx9iuUei0rWK0ceZFJEnSYlEU0+rtqmuF6hgJrtapPqWSJC32pC1iqIY1pbE84P8dr1fcwtzq4VzrvwQ4OrXw2IGDuWxOEaDgSP8IIGDEwgdhUwkfOoUKq5YB3Ryu6naR0VzWIYuy/N2YZQGLoses6ClfcjeZK3LIsI5kR+21WGUNdkVw/ZtotgKwNa+SLw91wqxokZ3Ps6Hjw9LL6bluMn/mDGHKjk7HvZ6bLyrDNzCK77cc4oOdnY87ft0l1WAK5KP1hczM6F3vSBjj/ZZRUDmWLzcXAgL5cjhfVFzMe3veJCt7NZHCWywsbE+3CB9qCrcDsXwcO50DcSJjV6XwcGoon24q45b+8UQFGHn+D8eE8ryxkczMrALgg9W5rDU4xlFrzCkA3Oa/gDfThnD7CF/sisDzhy7kttwIeilllOekkUUi477ZzoNDE3nn8h5UmW3kVZhJCvdrwifePHhjTBQGZNTbDgWQJOnHk1xzDdBRFMU+kiRtauqDu0UHU7vPyIU+61lqHsib5dcz0ieNcG1Zg+fXtUJ1yAgsqR3AbTumstfvSuLOHw5AatdupHbtRuY7hgbvM0a/mPsm/X5Cu17fbKRGgaNFCwurU8nc+ihP3/oUE6+OwC4r2BUFuwwHv+5D6bd7KRcUbrH7c2W0P4oioJgiibp6CbICev8AAJ6bMJ47K2pRFGhnqqTvh+msrulGp0kL0KLDhg4dNvoZdnCJfxohqZP5v/AB/Pl1KVd00vLKmlj+478I3wPf0RXwF77inbVwo99CInYV0tu4G5gEQKfCGfy6fSwAFkXDaqd49lrb8XjgN+ywdiCtMpgNcx2xiqWmzoz/aiMAcdoCHhwRDkBiqC8ALy3azRt/7+OPOwcyWow84XvoTU5LRKIoBgL9gGAcrUmaJElnvPiPJEnTgelNuTY9PZ2CAkfYTIjToyTZEnkpeDpXFbzKY8UP8nn4ZHSCzMaNG4+61n7M22NFz7+WzugqfmNBrpnhzvMfW/APiT5VPO/U0B/VqcToCull2AtAjWxgw4a1aDT6Bm3Mr3SUfjlWtHY0vFdyGU//8SjVnZ446ppwyx6X3kK0lYRQ6bwol+Ich4dxU72Qvzp5HwaKamUgCh02bM7n2dCxxSJSYAukNnMj1doL+OWSUGwyHCxRuK/sB9e9Hg78gQ8rrubZ4BnoBJnyqGv43PA9fWyLKNpvBMYyMsHE1sMWgrW19AqqZEdVEOcPe5S3fy2gzpeYFHik66cTFIqUcFbtOUyYSUNq5UI2bkzlr22HAfgzbSdhlY5UZ1VWmawKG11CG/7ROl3y890rcNAoEYmi2BuYgOOnMg2HgEKAp0VRVIAfJEna0shnFuFsfZw03Mn2ECkpKcTHOwav3bp359aVC9lm7chFvmt5IeRzni65j2dL7mFyyCf07dvXdV3U/CLXl7s+xcYuDH8sk8E2Owad1nHfDVtICPB3fEOB50vvZJTPepeIzsv7hCtWr+Wzhx522PTcZ1zTycpzNzsjH75zxNI1JNpNls6U5/1OWXwRowePptpi5fmvZjLI3Jk+xt1UySZmVF7KcNNGUgwZVMg+zN+RztjeXenXYyDF5WV8suAXxvbpQc8ufThcWgLk0BB2NHxQdgWPGBcTNDiITvGdyCoooGfOOuyljvFJji0cjSDzW9Sj6ASZQ7YINgbdwYTRPYgODWPThtUsV6roc875+FDNwc86oFgrwR9eWvsSKN2pU78/ZYCjiuFN/dszc8NBDln0BOgquGW5wqIHgyi1lXBd70jevL4PAHnltUz8bSe/bi+k8OWL0GpO7sjJLa/lum828v0NfYkONDV4TlaWe3kITykiURRHAEWSJD3ZwOG5znN6i6J4gSRJSxvxzB9xOhJEUQwGmtw9O118TCaS9Vlss3QEYILfEnJt4bxfcQ0lcgAz9q2iXdK5ACd1YwMuAQF8dM/tAGS+cw8AP0dORCscGWjfHfATg3rf5toWg+3Ehji+PIpy8ni1tS/cS+RzcTwXso/Rg6HGbOGdXRH4BXVyiEgx8U75tYRpykgxZFAu+/Pyv2FE+OygX4+B5JeW8OymECJ9ttGzSx9yCg+77m1rQLTrDJfSL3MM83b/S6f4TuzOyuB/q+Dr8BhidUVk2GJ5ufQOuusziNKWsMcaz6NLK+gXsYWYsJGgNTG07xDnHQOJv30f9qocMrb9wQ8LErHX667uqvDhivY1LC8IZFhSKDPWH2RcSgzYQnhhiZZJK0s4WFrDOTlp7MpN5M6524kOMDE33dHEbssrp1dsELKsoDmBmCYt2s3KjGIe+3U7s27o2+A57tKYlihNkqSGBw1OJEnaLIpig4VRnQ6DfkCpKIpIkrRJFMV9Tq9dHxxeu2ajuz6TpbV9URQQBHg46AeCNRW8WnYTfadlMKnrQm667CZMYWKTnxGrKzxq+46A+SSeN9e1/eNjR5yRp3KJB5j01Lx+JFQo1N8X6+uXsf/dawCI0JSyO248grOTFKs9TMWLQzDoHQsRk+PiKXzaiK+vIzype4dE4PhFiHVse3wYJYf3EhzqaL3P7d6D7P/uouKnnQAMMO4gLeZm/DWOINwhpq3sv01DZELDoT9aUwhaUwjvFdpRNAeh3sSwrAisyqqhs18RCRaZcB8NA+KDGSVG8n16MbsKLZwfa+PXnBh+fXMVAuCnk0lNCGHtwTJWZBTTKzaIZ//YxVvLMxjTJYJnRnamX7xj/i+3vJbP1h1EAb7fksMbl3YjJsinQTvd4ZQ+xPoCEkXxnHp/9z7Recfs/1GSpL6SJE2scwpIkjRdkqTFkiS91hxzRPVJTRlAsRwE4zYT999qBuTOZJthGGsfGESkv5F7tqeS8vpSJn/wBHvXTsdauq/R9y40JHNdwcsctgcftV/jG+Wx5QeCICAImnrboBVkNILi2vbzD0VvdHxZdFotoaExmEwO75ZO1/C4rA6D0Zeodj0x+oYAYDT5E9OhH3pny2oQbIRoK4/aTuh2CSM+WMJdn8054X0biqywosdfY+ZB3y8JXzOBVFbQM+MJFEUhOcKPomoLfz50OTOvdvygKUClTUOCLp8IPwP/5jiG4+sOlBBk0vHL9nw2Zx/5Gk5adCSjkwI89Mv2k772ptLYMdFspx0ZoihmSJJU7mx9xkmSNO+MWHaG6JecBOvy2bBrBx3bd+fqjrV8vqcDb1LIlmevYfaGHby2uIbn9rfjhf12uurnkxpwiNRogW7xUXRJSMI3uD26gHg0PuFHtSRfRH1HWuYBvkxcfFycXOGSB1yhL8dOoEb5Gxscf0X5G0/4OjS+UeRVmHmo6FHeC3uLCG2paz+cPF6tKc8rNCTzwKEbj3pW/eedp1vNuJo/ECwfNnj9pkeHUbjkASq2f3F01ITGgH+X6wjo8TZfHlyCLqArgiAQY6rl5/xKZEVhbfbRkeJ7zZEkhWvZX1yNXVbYcKiUG/q044sNWUiHHQ6W3PJaZm7IwlZvecdP2/JcgcCepFEikiTpGlEUE3FMmn7m/DsNaPzPdAuhr9gTDX+SdqCACcAjl4xi+jub+OD3xbxxbwrXDezOdQO7sz2vnNnrtvP3biuzDrfn82Id7AANduJ1aSRofyNcV0GkyUaEyY7RYODz/UOQFS2fr83kctMyEoO0BPgYMVFLxfYZgEzF9i/w6TAGnX8MgkYHGj0HH+4Igg5Bq0cQdPybvp1zevcDjQ7FbgWN7rhuX/u7spg6dytpeQ2Ltv7y9mNFe6rxXkOc7AcC4LHeUJGeRe3BGTCo4fufKOzIcngLpthBmGIHuXaLugNAEHPnvMzMf/sfdcn2/ApmXRFDQmwHpIJKKs12UtuHsCqzmN0FjjmqSYt2Yz9mvCmfJBDYHRrt4pYkKVMUxSmSJG0GV3cu4xSXtTj8fEwkmwrZUuDojnRuF89FUX/yWUY0T+bvIzwqCYDu0YG8dLnjQzXb7KTnVrArt5CdWYfYla/hYFkoadUKBSV6zLL2qGdYZA0XLe3o2haQ0fMVBsGKQbBhmFaAgWwMghW9YMOAY79esKFFRiPIaBYsd/yN7OiuoaBxdt20KJgVPQurByCjZfrqDIrSv8FfY0YQFMf5tio0TECzsgzdunsdUQ51EQ0ah0v5yD0d/+sExdk9dERUaDUCGo2GCtmHzzOHIStaPluTQUrJd0T72gn3NxEV5E+oUXD9SBjz5mMu2II+tAsa3dG/+KcTIX/v5RNI9XmZiD1TeT5gFG+VTaBQdnQxZVnh/35ZxkNDYvGLHw5Aj0gfxAh/0p2LDNceKMF6TPdR4czkGxRO5R1qrcTFxXUAMteuXetycddx7dvTWJbnS96r1yMIGjbs3s2g6Tv5YlgZN15602k9R1EUdh+u5Jy3VmC2HYkp02sFXrggAdlcSv6/M7HawYIei6LHigEharBj2y5jtjkSSlrsCnYFamrN6PQGx8Sq4giFkRXHmFxWBOwKFFt0VNu1ONzFCgZBxqitf76AjIAd7Qlt9xQCMiGaCkI15YRqygnWVBCsraB3cBXnhMmIET6EJ52PIaw7Wr8YtKaQRt1XURTef+c6LhXmcdAWxf+K72eDxbEa2SAonN/BwIPDU/hu/udMuiCKbysu4Ped+ay4fwiCIDDykzVUWWyseXDoSZ+TlZVFamoqQGJ2dvb+0339jXFxN8p1fRoubq/TLz6Y2dm+ZOxPJymxF/07d2bvoyY6xCac9r0EQeDdlZnHuaoFILsKXgqaQ0XQ7OPGAQFdrCcMLt24ceNRc1bHkutciAd1ohXQ6vRIT40gXFvCoRkiiv3IOELR+hJz8w40vlGOqAcZ7IqCrCjObcW5zTHbCnnlZsZ8ug5zvaBTo1bD1//pjaxAbmEuGWvep8jmT7EcSJE9iEx7HEVyNLOz9ZDtEFnC8ny66b8mRlvMNZ0sDOiSgj4wAVPCiBOKShAEHnzke6zlB4gqlljbfpSrW7tkz2FGTVvLoyP1bLck8PauIN65KZmnRiS7rt9bWMV5HR1TkgUVZp5euIub+rbjvKSwE763TeGUIpIkaakoik8A+xpyIoiieBWQKEnSGx617AwytHs3WLufVelbSUp0BD7WCSjrkER8u9Nzb59sTY856vRWyjaGky38eyno0+PWDAmKjaq0KU2KCH9vZSYKx3aLFJbtK+TDcT0pXDKFiqAfj3cWdL+V0p5TSM+rYFtOKRszfdh5OIqlFTBji4bQrWW00x5kTMDDPDA4Dv+wJHwSLkAflHicDfrA9ugDHevBqvb9SsXW6SQOcSSP+Wd/CSP6DCAlJuCoa2qtdrLKalwxdz56LTPWH6RDqE/ziwhAkqTXnROqr+II+QnBEWlQiiNaYe5Jb9DC6J3cDR9BYvX+Qm6ut/+DBb/zv+XVbLrpAF1SGj/4PvmaHvdWyjZEc4r2VIv+TugsyFtL8gh/kiP8HROodAWgpNrCd5uz2ZxVyM7cMKZmJ/Pxb5Uk6PIYYvo/Hh8aR4dBD5Fv8eU/szYdF2kgaHTUHFiEUDyEcP1LvLxoN3emJnDDOdEcOrCJcT9V8lBqGHHSiyjKTSQ7RRRg0hEfbGLv4aomvQ8n43QcC5uBzR63wAsYdFp6B5WwvuDo2KvLBw1h4vLFPP5TGvO7DUfQeiY2y9M0p2hPteivzllwzqTvSAkXeHhA8km7oiG+Bu4bkgg4Wpy/9xbyddp+Mgsj+XR/Er/9WcAdq29lubkfq2rOYdIiiQ+v6uW63jfxYiLHzqZyx1d8UPsN12beRapxB3l/f8GSzZvYmP8EW7ZsIqV6MV/1TWF40ijXte1DfDlQUnOsSW5zyslWZ9BpmyM1zsj22ljKSw659sWHBfFYfxO/l6ewYOEnXrSu9TEmMJ1BAQdO+7rhncL5/Np+LP3vCP55YCgBwXG8WHony2p6IyMwY20Gh3IzUeQjYVR+nS4n6rK5jL9nHpkXreC6/sl8bxnHrXlPEOGrIzM3iwBNFecWvkGk7sicVvsQH/aXVHvk9danMauePhVFsdepT2tdnNs5ERs61mxdf9T+p668iPamCh5aZaIs918vWdf6mPLQK9x3w0S37jGwfQjbnryI//SOReeMhTMrWq774BsWv9ebopXPHOXAEbR62o98HVN0P+44fzA/3NiX9vrDSJYE7il8guXV3Y/K65Ac7o+vXuvx/HqNEVEI8LooikWiKP4giuIdoih28KgVXuDcFEdU8D97jo7g9dFr+fS6gRyWQ9gg7WroUpUzSF6FmXnpefUiDQT+MZ/DhYcmM3phCAdXvAwcncXVWpaJbdMkrugoE2XZSaE9iMW1A3i7dDyVO75yhVo9P7oz2/93/gmDVZtKY8ZEUyVJWgKuCdaRwI/OJRAjmmM90ZkgPNCPJFMxa3Ntxx0b2T2J/c9GERbo7wXLWic9XvqOlHB4PPX4FbanQ0OeR71GYEBCMOsOJHPTX1uZH/k9lpx/XBEZPvHDKV0/hdrs1fQz+mNXYGltX/oad6HYLY1OP9ZUGhOAuqTe35slSXpdkqR+ONYX3XXGLGsGBkbJpFXEYK0+fhY7LNAfRVH4dMFccnb95gXrWhfjQrdxfsihU594ChoMVJUVqix23r+iByvN5zBx9l+UbZ2OI4xqJoaYQRhjh1Cbs4ob/X7j0cDvMCtGRN1BdpnbUev0TBZUmBk9bQ3zt3s2B3mTl4dLkpQhiuKSU5/ZchnsVXOSAAAWvUlEQVQuxvPtgRo2pq8kdeBlxx3fU1DO/X9r+WnDen5+pAuG4CQvWNk6ePkBx9jj2NXBp8vJvIGKorA+q4xpaVcys+ISuhkyeDroS/zXvEjEqOnkz78Ka/Eudlo7ALC5/SRe3V1E0X8uBBxu7sV7ChmWFMal3T2XRLLJ6VScyUeOnxlrRYzu4whsXLJ9T4PHO0cFMWlEDH9U9eONGVOw16p1Xr2JIAh8Mjaad8Pf53r/P8m3h3F30ZPs3boAwRBA3I3/0u5WiRURTwFwz+COfHxViquKho9eS6ivntzy4yPY3cGdlqifJw3xBgnhQSSaylh+SOaZE5zzxIUDeGHZAp7NuZxnn1111LEof2OTIqLbIp1f+J5e4TJPDm76YsbGULZuCmN9VzPW9DfX+i3i0vzXebdsPO+tmEjkmC/RByXy5S0J/F9JNQn2HdTk/07FBh0hqc8BEB1gJK/Cs0n4W0dirzPI0FiFdRXxmCsaLuSl0QiubJ/H0tCanLOVm6LTuTii4dwNnqR+hESSPpsJ/ouZVz2cPXvWuc7x0WuJzv2W3B+Gkpl7iI2HjswVRQeYyKtoIS1RW+GCbh34KqOUdVtWcN7QCae+QKVBnr3XM2OiU3HscopJpTV8N/kvphWNoNfuOfgmXoJG74sh0rEI+9vcJFYr5zmSGwJ92gWRXaa2RB5l5DmO+aKlO/d715BWzrEV/JqL+GAfbu4by5zqkWz+9WEOfBRK7twL0QclIej9uC3hIB9c2cN1/mtjuzHLmTnIU5z1IooN9qOTbxkrsk+deach2up6rNOl60s/MOGDb73y7KdGdcdoMHJL1Yd85PsOC/dUULTiMYwxgwgs+YeYdROozVl9xp5/1osIYHiCgfXViVQXpJ/65GPIX3QPinz8hO3Zxh1x27k8xr0kiE0lMcyX3+5IxeQTxNu747ij8Ck+TStCH9gBa1kGiwsiyCl1xMwt2n2Ybq8tY2+h56K5VREBI1K6Ua348M/mFQ0eP1kCj/tWBZAz/1pk29lddnHi3f/HjeMf9drzhySGsv1/51P5ysVcmOTPpNLbKCEU69BvuH3fVayqdngNNQLsKqgku8xz0dxnvWMBYESPZJiTybJd2Yy+6PjjJ3Jjv7JkD88uhAsOa7hfc3a/ldUWGxo3y8p4ApNeyyuX9aHv2yv4aWcZj9x3OcuCS+kR7Vi0Fx3gWJvkSQ/d2f3JOwn3M9I9sIrlhwORzWVojA3moTyOp0ckM6h9CMM6XoKg0WCrOIS9KhdjdP9TX9zG6PV/s+kVZuepc7t52xTOiQ0k2b+Gn/LiuH//b3Te9wu16cUEXTaP6ABHr8KTIlK7c05GJ4ey0SxyeN/pRTKd3ykcjUbDnsOVXPDu76R9N57yrdPPOofDfQm7uKpdkbfNAByRDZclG0kzdyFr0UPkZW9nYbkjTVaIjx6NAIVVllPcpfGoInJycd9eWNGzaFPTUoMfKqtla007xhVMZfHCdyn47dqzKkzokTte5rqrHva2GS5SOydiR8supRvl1VU8m+HID67RCFzZI4b2IZ5LJ6x255yc2zEcX62NvzKruV62ORIrngbndwpn9QNDuWzGeq4veoVX09/nyvz+KNYq5Nrjf6E1vlG0v8u9agQtibIa6ykrNDQnqcnJwAH+rY7jJt+NrH3wXNexOTd7NmJNbYmcGHVahrXTsLyyGzWHVjbpHt2iA1j70LkMaB/OI0UP8Yc8tkEBwdGV+doC/abM4ZZPvDNP1BBxwT7EGCr5tzISzEXIP/WmJmv5GXmWKqJ6XNyrC1n2KLZtXdTke4T7GVl09yAmj+nCDTc7soidDcOjhxIlJrRv1toEp+ScEEcXG2Ch+Xyyqh1OhfvmbmXQe037oWwIVUT1GNPd8Yb/sSPbLceAQafhqRHJ+Bp0VMomxhW8ysLqVE+Z2SL5760vMf6KB71txlF0jTBy0BZJxXnzeHDvKDbXODLhygpkFnsuYYkqonp0DPMjKdDO3+WJWPI9E0hZpfggoPDf4id4tPhBSuUjS87rZ7Bp7RRWmSmtsXrbjKNIjgrBhg673YLU9z2uTHJ83cP9DBRVWz2WsEQV0TGcKyxnrbkHmd8NI/Mdg+vfgenxp764AaK0JfwQ+SwPB37Pb9VDuCjvHZbWOPOyCW3n7R/06o/cPm2Wt804ii7OHOzpa79ByV9N/o/nAxDhb8AuK5TWekb0bedT9BBDdWuoVYxsMHc9an9THQEa3yj0gp0HAucwL/JJQrXlTKu4AsEUgSAI2KryyJ03hppDDYcctRYe67SXG5I8n13UHbomOBZe7y11tDjm8ixsVXmE+zmScnpqrkh1cR9DqnEbBiwsq+3Luaatbt+vvhs7ERhls1NaYyMqYCI5ZbXMXrWdywt3kvfjSHzajyZk8MtuP9Mb3HPzC8CZX090OkQE+BKgqWa/LQYAWVEoXTeZrp0ncXO/dhi0nmlD1JboGHw0Fgab0llS0/+MeNWMOi1RztCTL9OyeHRZNVdXzmRf13cw56eR810q/tsfb3WR4TlltRRWtayVvvbqPBK1OWTaYgHQY6Nyx1f0DKnli2t70yHU1yPPUUXUACNMG8iyR7HH1rRxUGN58oJO/HBjX/IrrYxeFM+z/r9Q2WsKsjHKNdlrzttwXJWHlsiwN37k7uktZ54IHNUC2+kKyLZHuPYpsp3SdZNRFMVjjgW1O9cAF/hs5LlSWFrTj876MxdVIAgC43vFcpEYyStL9/Duigx8jCOZl96F4g3z653pyHsX5a8n98UGwsxbAE+LmUSGhnvbjKMw560lRtuLJTVHqsUjWyg6tIl2T/7O5DFdeGy4+2nQ1JboGDS+UURri+mh38vi2n5H7T9TBJh0TLm4K9KTFzDpoi4U1zbc8uRXWin6+zHs1QVnzJamcusNz3HJxfd624yjiLt+A8ndRmHGyLIuP5Oc/SOam0pJvmEVNlmmxEMuebUlOoY6R8BVf0m8/JdEsb4jfe5JR9CevHS9J4gPPnVQZOWOrwhOdST4spYfQOcX0yJKwBworsaoa3m/yQnhwYAFMcJI7asXoRUUBI1AsI/eY/NaqohOwGXdo3npr90sLW1P14NL8E1sGd2o+Dv3o9E7Clcd/v16rGUZ+IvX4d/9JowRvTgwPb5Bd/yZDngd8fZceoVaefq8llVAJCEyGjjI3i0/Ebf6MkLPnUxQn4cJ8aCIWt5PRwvhnNhA2gWZWGoZTOWuljNg1uj9sMsKOWU1BA98Bp92wyhPn0bOrP5kzxpwwvmsMx3w+lL3g9zRtWVFLAB0iHGEch3Iy+If7RjWWRwiD/bRe6w7p4roBAiCwCXdolhZ25OSPQuRLZXeNsnFT+m5JL6yhP+ui6GgzzQS7jhA2PB3nCNn73D9tc8w5qK7vfb8ExEdEooeG7n2SPaUyGytddTmvaFvO8Z288w4VxXRSbisWxTVdh2rqzpStXtOsz031NTwx1KXMKVvu2DuSm3P7H9z6PnmcsZ+s4cVhquIufbkdVnP5NzTnsOVZJV6vpSju2g0AtGGcnLlCO7r7cOjQxwTrw8N7ci9gzt45BnqmOgknN8pHD+Dlr8Zxej0TwnocWuzPPfPK6NPWvc0McyX969M4cXRIp+s2c9Hq/fz+PztXNwlEgCzoscoHN9VyfoskXY3b2t0DonTYcy78+gVauHpYb09fm93iTXVkm8JoTZrBnlzLyT22pXYZYUKs41gH/cdRmpLdBJMei2jO0ewpKYv5rw0zAUtq+5zmJ+BZ0Z2Zv8zI1lw+0A0GoEa2cB5uR/zcNHDrDN3Oyrqwr/7LS4BFS1/nJK1k7CW7vOILZN75XBfz5azsrU+UT4yhbZAiuUgZvBfAJ7+fSfRL/7lkVwYqohOwdhuUeRU69hu70JF+ufeNqdB9FoNnZyl5q0+cYz1/Ye/a/vwn8OTuCDvQ94vH0+hIZnQIZMAUBQZa7FE6dr/49DMruT/Oo7a7H/c+kJNGP8kI0be7pHX42mifDUctgcx1e8b/i51jIlCfPVY7DI1VveXo6jduZMQ8+JfrsoPl+dNhoXAwvktuqRKr3slZgDvm238tC2PL9OyeHdvNNfe8gEA+RVm/A1aoq+cj63iEBXbPqf830+ozjif0PNeI6hP05KN7MirwNeg9eAr8RxRgUYqFD8+HdcBX50GRZEJcXbjSmvcHyeqIjoJJyqd0hpKqvgZddzQtx039G3HodIa4oIcSQuf/WMXP2zJ5rJu0Uw4J5bR/Z+l25+p5Fda4Fvg2yPhRqfzY3HZBz/RK9TM08M8myzeE8QEOgJNt82eQHTFKuLvPEiIj2OCuqTGSqCb92/xIhJFMRhHseVQIEOSpMVeNqnV0a5eJMRt/R1BtfO25vLt5myCTDrKahv+Nc6vNCPbatDoTh1J8UbfPEKCIk55njeICQkGzOQUFvCD/xRethsI8XV040pqLAS6qYLWMCa6BtgkSdJ0YKK3jWntDOoQyqfje5H7wmgW3D6Ay09RuzT76z7UZq866TkAV1z5P4Zd0Dzey9MlNswRGHtYDiajTKBc9qFzuB8vjO5MXKD7+eeapSUSRbEjUCxJUqlz+2qgFOgIzK7b3xBO8dS1SBnNYO5ZgUGn4eKuUVzcNYqvNp646veCshTOnT2WpHFz8G0/6oTn/ZtThp+hZXZsYiOigWwK7SHMuiIYP18bGkMAL4x2JLnPqnYvc+sZb4lEURwJzMEhmDpBdZQkabFTIFMbeau7aEEtkadcwy2d+w/dyI4OL+HTbhj5FWaySxvOknPVRz8z8cuWlWOhjqjQSARkDtuDKVx8D5U7v0VRFHLKaimudn+J+Bn/6ZAkabGz0ngdda1QHSPB1TrVp7Ru/OM8Nh2HEJuW57cJRPkbG3QihGtKKd/yAWHD324uU7zG8vsG0z/+YgStlg+Wb2by34foFQ4X9+rEmC6RpCaEoNNqeHdAIcEtdEyk1+kI01ZSIIexQPMfOln7MEpWaDdpES+M7sxtXd1b4eqN9jeMo7tloQCSJP3Y0MnOlmwqDuGlAc0WoHWsZ6r3W8vxN2iZl/wzFdvnEZz6PFpTSHOZc8Y40Y9FlL+RoR3DXNvXdvfDJi1haVEMry1VmLJkL3FBJvY/M5JLLnsMaFk5FuoTpq+h1JjEO5XdeEmIQafVEGDUeSQItWV2YuvhbI3cX37oAcalxPDiXxLVl9yPsvNrKtI/Jbj//7xtlts01o3dPbEzkx98haKlD5CzfQobAm6istNdaDUC323KJsxPT9ipb+MVwgxWii06pCe6g+LwzIX46Cmtbp0iKsLZ+jgpPpMPS09Pp6DAMytBRa0VRYFp6yq5JXggRevfYp8yBLQmj9y/Pi31Fx2AsHsxdo7hgn1vw9YZZP5r4dm89+iuz+TdsLfJXAmyPpTS1N+8bamLAI2Z3dW+HJjZE1tgTyq7TcWIlf15haSn57p1b2+I6EecXTKnx+2MjnFSUlKIj/dMwpE+isLzG5axsVzP82Onkjt7OJ206wnq+4hH7l/Hxo0bTxqA2jLohznvCnK+H9zgUY21uEW9htjV61lX7k+BHM4ivyd4vm9f4tavxi4rpKQkuHXv5vDOXQ30AyaIothHkqQMYJ9zrHMXcOeZtsFTCILAlSkxLNtbRHVwP0wJIylNewPZ2rKSFjYXxmjPlig5k0T46SmVA8BeS5EmDoAHzk3kwaGJbt+7ObxzP+Jofervm+78s9VFH1yVEsPUpXuZvz2fCanPkzv7PMr//Zjgfo972zSVkxDu51iLpdcqvD7UEW1+ZYpjbVFWlnvL5ltDxEKLom+7IOKDTcxLz8UUm4pP+9GUbXyrRa18VTmeyEBHlHtRrYbs7xxV8wqrzGw65H45GFVEp4kgCIxLieFP6TDltVaCU59HrimkfMv73jbNqzwd9CU3+v/ubTNOSESQYx1VrtCeiUxDURQ+WXOAfu+sxGp3LzmmKqImcE2vWCx2mV+352OKGYBvx0sp3fA6tir3vDytkbp8fBf4bKSfUTpuf0shMtgxn5cefAOHNfEIgkCwybEcotzN5RCqiJrAwIQQ4oNNzPk3B4DQoa+i2M2UrH7By5Y1P+3vyiLxYQurQh9ja9JrFA9dQ+LDlhZXj7YuO2ucP/wxzoBsrSLE1yGiMrN7c0WqiJqARiNwdc9Y/pQOU1pjRR+STFDvB6jc/mWLW0LeXNyxdQD/XXvyiHBvEhniaBmzM9aQ8/1gLIXbXAvzTrQUpLGoImoiR7p0eQAED3gKjU84Rcsf88i6/dZGtJ9CSpTnJ509hdFoJEBTTZFFy4fyI2Ta2xHk7M5VmlUReYUBCcEkBPu4unQaYxAhg17EnL2KKukHL1vX/Oht5Zgsed4246QEa2sol/0J0MvIukDESD++vO4cOkf4n/rik6CKqInUVXT4a/dhSpzh9AE9bsMQ1Y+i5Y9hr3FvjUprQ5btmK2eqTx3pgjWWSiT/bivaxWibwnhfkZu7BvvqhfVVFQRucE1vWKx2hV+2e5I0StotESM/ATZXELxiie8bF3zYpcVyqs8V5H7TBBksFEqB1C57TMqdn6DLCusyiwiq8S9pJOqiNygX3wQHUKOdOkADBE9Ce73BJU7v6E6Y4EXrWteJkd9y+Ptm22pV5MIMciUKwF8Z7+Rn81jUIDzPlzNXDcDUFURuUFdl27R7sNHrZAMHvA0hoheHP7rTmyVOSe5Q9theFg+g+K8X+LlZAQbBcpkP/72vwPfwEi0GoEAo44K1cXtXa45JxabrPDztiODakFnJGLM1yi2ag7/dTuK7H6CwJbOL8p4FpsHeduMkxJi0lBu9+WXyw1cEe1ofQJNOios7n0+LX5RnhtoAXJyzmxLEKEodNBXMevvLYyKqZ9G14+aLs+T888zFPIIwf0bP0bKz893OyiyuXl8bTvCDBamX9Rybdebq7BXlrHms8sICk8i8pJZ+JtLKMp3jYmalH1SaKtzGnFxcecCK71th0qrYmh2diPygx1DW26JNgBDgVyg7fenVNxBC8Tg+M6cNm22JVJRaS5Ux4KKipuoIlJRcRNVRCoqbqKKSEXFTVQRqai4iSoiFRU3UUWkouImbXmyVeUM0piaU6dTh6q5aMDuaThK9vQDVxWT07K7TYqoJX54J8MTH2RzUq9Sx53Apno1p15zHp8miuLUY/fRjBU9GuJYu527RwKZwHRJkiY29Fo4hd1trjvnRhExb1L3QY5yCqhFvwZnpY5T1ZxqsA6VN2nAboCJkiSFSJJUV0DutO1uiy1Ri/vwGsHEY+oztbbX0FDNqQbrULVAOjpbqD440l2ftt1tUUSt5cOrj9sfpErTqNdtKwam0YQqJW2uO9cakSTpNWdXYzGOD7K10VDNqWatQ9UURFG8WhTFuiptGTjGn6dtd1sUUYv/8OrjqQ/Sy/yIs5phvZpTDe1raZRypGJJXT3g07a7LYqoNXx49fHIB9mcNKbmVEusQ9WA3YuBkc79E2ii3W1yPZEoinfh+FXvg8N12WLcww3htLcY6A9Mcc6xtKrXcDbTJkWkotKctMXunIpKs6KKSEXFTVQRqai4iSoiFRU3UUWkouImqohUVNxEFdFZhiiKc7xtQ1tDFdFZhHOJhYqHUUV0liCKYh8c65IynH+reAg1YuEswtmVu1MNIfIsakt0dtFRFZDnUUV0luDswqXV+1vFQ6giOnuoC/G/WpKkFrW0orWjjolUVNxEbYlUVNxEFZGKipuoIlJRcRNVRCoqbqKKSEXFTVQRqai4iSoiFRU3+X9OEULXs9xpZgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 180x180 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mkr=['^','s']\n",
    "sty=[':','--','-']\n",
    "col=sns.color_palette(\"colorblind\", 2)[::-1]\n",
    "\n",
    "fig,ax=pl.subplots(1,1,figsize=(2.5,2.5))\n",
    "for ait,alpha in enumerate(avec[::-1]):\n",
    "    for bit,beta in enumerate(bvec[::-1]):\n",
    "        for dit,d in enumerate(dvec[::-1]):\n",
    "            ax.semilogy(J_store[ait,bit,dit],linestyle=sty[ait],\n",
    "                         marker=mkr[dit],color=col[bit],\n",
    "                         markevery=(int(29*np.random.rand(1)),30))\n",
    "ax.set_ylim([5*1e-3,1])\n",
    "# ax.set_ylim([5*1e-6,1])\n",
    "ax.set_xlim([0,150])\n",
    "ax.set_xlabel(r\"$t$\")\n",
    "ax.set_ylabel(r\"$J(t)$\");\n",
    "# fig.savefig(\"J_tt.pdf\",bbox_inches='tight')"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}