es_brachistochrone.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "based on https://github.com/hardmaru/estool/blob/master/simple_es_example.ipynb\n",
    "\n",
    "# The Brachistochrone\n",
    "\n",
    "The Brachistochrone is a physics problem posed by Johann Bernoulli in 1696.\n",
    "\n",
    "What path gives us the shortest distance between two points? As most people know, in flat space the answer is a line. \n",
    "\n",
    "What path gives us the *shortest time* between two points? Assuming we are in a uniform gravitation field, the Brachistochrone curve is the path that minimizes the amount of time it takes to get from point A to point B. What does this Brachistochrone curve look like?\n",
    "\n",
    "To set up the problem, we assume that point $a$ is located at $(0,0)$, and that point $b$ is at point $(x_f,y_f)$ and all points on the path have $y \\leq 0$. Starting with a particle at rest at point $a$, and using conservation of energy $E = T + U$, we know that at point $a$ the total energy is 0, and at any other point $(x,y)$ on the curve, the energy is $\\frac{1}{2}mv^2 - mgy = 0$.  From this we get the following equation for the velocity of the particle:\n",
    "\n",
    "$$ v = \\sqrt{2 g y}$$\n",
    "\n",
    "$s$ is the path the particle takes, so we know $\\frac{ds}{dt} = v$. We can then rearrange the equation to get $\\frac{ds}{v} = dt$ and integrate to get\n",
    "\n",
    "$$ \\Delta t = \\int_a^b \\frac{ds}{\\sqrt{2 g y}}$$\n",
    "\n",
    "We can express our goal as trying to find the path $s$ such that $\\Delta t$ is minimum. Rearranging the variables we have\n",
    "\n",
    "$$ \\Delta t = \\int_a^b \\frac{\\sqrt{dx^2 + dy^2}}{\\sqrt{2 g y}} = \\int_0^x \\sqrt{\\frac{1+ y'(x)}{2 g y}} dx = \\int_0^x f\\{y,y';x\\} dx $$\n",
    "\n",
    "So $\\Delta t$ is a functional (its arguments are functions). We need to find the function $y$ and $y'$ that will minimize its value. Classicaly, this is a problem in the Calculus of Variations. Solving it requires plugging $f$ into the Euler-Lagrange equation and solving the resulting differential equation to find $y$.\n",
    "\n",
    "$$ \\frac{df}{dy} - \\frac{d}{dx}\\frac{df}{dy'} = 0 \\qquad \\text{Euler-Lagrange} $$\n",
    "\n",
    "Plugging in $f$ into the above equation eventually leads to the differential equation of a cycloid.\n",
    "\n",
    "$$ \\left(\\frac{dy}{dx}\\right)^2 = \\frac{2k-y}{y} $$\n",
    "\n",
    "So in other words, the shape of the Brachistochrone curve is a cycloid. Cycloids cannot be expressed in closed form, but can be expressed in parametric form. \n",
    "\n",
    "$$ x = k(\\theta-\\sin\\theta) \\qquad y = k(1-\\cos\\theta) $$\n",
    "\n",
    "k is the radius of the circle the cycloid is being drawn by. k is found by using the condition that the cycloid must pass through the point $(x_f,y_f)$. \n",
    "\n",
    "## Approximating the cycloid \n",
    "\n",
    "We can expand the cycloid function in terms of any set of basis functions we want in order to get an approximation of the function and its derivatives.\n",
    "\n",
    "### Fourier basis \n",
    "\n",
    "One possible basis is the fourier basis:\n",
    "\n",
    "$$ y(x) \\approx a_0 + \\sum_{i=1}^n a_n \\cos{\\frac{xn}{2k}} + b_n \\sin{\\frac{xn}{2k}} \\qquad y'(x) \\approx \\sum_{i=1}^n \\frac{n b_n}{2k} \\cos{\\frac{xn}{2k}} -\\frac{n a_n}{2k} \\sin{\\frac{xn}{2k}} $$\n",
    "\n",
    "And we can get the fourier coeffcients with the following equation:\n",
    "\n",
    "$$ a_n = \\frac{1}{2\\pi k} \\int_{-2\\pi k}^{2\\pi k} y(x) \\cos(\\frac{nx}{2 k}) dx \\qquad b_n = \\frac{1}{2\\pi k} \\int_{-2\\pi k}^{2\\pi k} y(x) \\sin(\\frac{nx}{2 k}) dx $$\n",
    "\n",
    "Since we can only express the cycloid parameterically, we change all the variables to $\\theta$:\n",
    "\n",
    "$$ \\text{d}x = k(1 -\\cos(\\theta)) \\text{d}\\theta \\qquad x$$\n",
    "\n",
    "$$ a_n = \\frac{k}{2\\pi} \\int_{-2\\pi}^{2\\pi}(1-\\cos{\\theta})^2 \\cos(\\frac{n}{2} (\\theta-\\sin{\\theta})) d\\theta \\qquad b_n = \\frac{k}{2\\pi} \\int_{-2\\pi}^{2\\pi}(1-\\cos{\\theta})^2 \\sin(\\frac{n}{2} (\\theta-\\sin{\\theta})) d\\theta = 0 $$ \n",
    "\n",
    "All the $b_n$ will be zero because the cycloid is an even function. We can then plug in this $y(x)$ and $y'(x)$ in our $\\Delta t$ equation to actually see what the minimal time is. Of course, to get a good approximation of the cycloid, we will need quite a few terms in the fourier expansion.\n",
    "\n",
    "### Polynomial basis (Chebyshev)\n",
    "\n",
    "Another basis is a polynominal basis defined on the interval $[-c,c]$. The Chebyshev expansion of a function can be expressed as \n",
    "\n",
    "$$ y(x) \\approx c*a_0 + \\sum_{i=1}^n c*a_i T_i(x/c) $$\n",
    "$$ T_1(x/c) = 1 \\qquad T_2(x/c) = \\frac{x}{c} \\qquad T_{n+1}(x/c) = 2\\frac{x}{c}T_n(x) - T_{n-1}(x/c)  $$\n",
    "\n",
    "Where $T_1(x)$ are the Chebyshev polynomials of the first kind. Taking the derivative of the polynomials of the first kind produces the polynomials of the second kind.\n",
    "\n",
    "$$ \\frac{\\text{d}T_n(x)}{\\text{d}x} = nU_{n-1}(x) \\qquad U_1(x) = 1 \\quad U_2(x) = 2x \\quad U_{n+1}(cx) = 2xU_n(x) - U_{n-1}(x)  $$\n",
    "\n",
    "So the derivative of $y$ will be \n",
    "\n",
    "$$ y'(x) \\approx \\sum_{i=1}^n a_n i\\,U_{i-1}(x/c) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sys, math\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.integrate as integrate\n",
    "import cma\n",
    "from es import SimpleGA, CMAES, PEPG, OpenES"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#constants\n",
    "XBEGIN = 0.0001 #don't start integration from zero-- messes up integrator\n",
    "YBEGIN = -0.0001 #don't start integration from zero-- messes up integrator\n",
    "XFINAL = 10.0\n",
    "YFINAL = -5.\n",
    "TERMS = 20 #how many terms to include in the expansion\n",
    "g = -9.8 #m/s^2\n",
    "\n",
    "# Below are the optimal parameters gotten by solving the parametric equation with Mathematica.\n",
    "# Can we do the same with ES?\n",
    "# For reference, Mathematica equation to numerically solve the parametric equation of a cycloid is:\n",
    "# FindRoot[Evaluate[{K*(t - Sin[TFINAL]) == XFINAL, K*(1 - Cos[TFINAL]) == YFINAL}], {K, 5}, {TFINAL, 3.2}, MaxIterations -> 1000]\n",
    "K = 2.585999608432 #radius of the cycloid that draws the brachistochrone curve. \n",
    "TFINAL = 3.50837   #angle of the cycloid at (XFINAL,YFINAL). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.special as special_funcs\n",
    "C = XFINAL\n",
    "\n",
    "def get_cheby_path_func(theta):\n",
    "    order = theta.size\n",
    "    def f(x):\n",
    "        cheby = np.array([special_funcs.eval_chebyt(i,x/C) for i in range(order)])\n",
    "        cheby = np.transpose(cheby)\n",
    "        return np.dot(cheby,C*theta)\n",
    "    return f\n",
    "\n",
    "#use chebyu for the derivative\n",
    "def get_cheby_path_func_deriv(theta): \n",
    "    theta = theta[1:]\n",
    "    order = theta.size\n",
    "    def f(x):\n",
    "        cheby = np.array([(i+1)*special_funcs.eval_chebyu(i,x/C) for i in range(order)])\n",
    "        cheby = np.transpose(cheby)\n",
    "        return np.dot(cheby,theta)\n",
    "    return f\n",
    "\n",
    "def delta_time(path, path_prime, xb=XBEGIN, xf=XFINAL):\n",
    "    #equation gotten from conservation of energy\n",
    "    #due to conservation of energy, values of f(x) > 0 are impossible.\n",
    "    #they will create a negative inside the sqrt and return an NaN\n",
    "    def functional(x):\n",
    "        p = path(x)\n",
    "        if p > 0:\n",
    "            p = -0.00000001\n",
    "        result = np.sqrt((1+path_prime(x)**2)/(2.0*g*p))\n",
    "        return result\n",
    "    time = integrate.quad(functional,xb,xf,limit=200)[0] \n",
    "    return time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 234,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What are the chebyshev coefficients for the ideal cycloid path? \n",
    "#n is the number of coefficients to approx with\n",
    "# x = k(deg-sin(deg)); y = k(1-cos(deg))\n",
    "# dx = k(1-cos(deg))d(deg)\n",
    "# limits are 1 - deg = sin(deg)\n",
    "lim = 3.50837\n",
    "\n",
    "def cheby_cycloid_coefficients(n):\n",
    "    def f(deg,i):\n",
    "        x = (K/C)*(deg-np.sin(deg))\n",
    "        y = (K/C)*(1-np.cos(deg))\n",
    "        coef = 1./(C*np.pi) if i == 0 else 2./(C*np.pi)\n",
    "        cheby = special_funcs.eval_chebyt(i,x)\n",
    "        results = (coef*cheby*(y**2))/np.sqrt(1-x**2)\n",
    "        return results\n",
    "    thetas = [integrate.quad(lambda x: f(x,j),lim,-lim,limit=200)[0] for j in range(n)]\n",
    "    return thetas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 235,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.044442696873506306, -0.0, -0.011447188632346877, -0.0, 0.008658537392531045, -0.0, -0.0044668242738176555, -0.0, 0.002855195474102763, -0.0, -0.0020064688185464684, -0.0, 0.0014978925121412474, -0.0, -0.001173656592045836, -0.0, 0.0009110581607894403, -0.0, -0.0007810377769037651, -0.0]\n",
      "-0.1064483732330774\n",
      "-0.5039518942760242\n",
      "-0.12013772246576382\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/anaconda3/envs/handsonml/lib/python3.6/site-packages/ipykernel_launcher.py:14: RuntimeWarning: invalid value encountered in sqrt\n",
      "  \n"
     ]
    }
   ],
   "source": [
    "a = cheby_cycloid_coefficients(20)\n",
    "b = get_cheby_path_func(np.array(a))\n",
    "c = get_cheby_path_func_deriv(np.array(a))\n",
    "print(a)\n",
    "print(b(.0001))\n",
    "print(b(10))\n",
    "print(c(10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 228,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"stuff for fourier series implementation\\n\\ndef cycloid_coefficients(n):\\n    # n is number of sine coefficients and n+1 the number of cosine coeffcieints. \\n    # change of variables to theta is necessary since cannot write in terms of y(x).\\n    # sine terms are all 0. cosine terms are non-zero.\\n    coefs = [0]*(2*n+1)\\n    def cos_coef(c,theta):\\n        return ((1-np.cos(theta))**2)*(np.cos((c/2.)*(theta-np.sin(theta))))\\n    def sin_coef(c,theta):\\n        return ((1-np.cos(theta))**2)*(np.sin((c/2.)*(theta-np.sin(theta))))\\n    for i in range(n):\\n        coefs[2*i] = -1.0*(K/(2.*np.pi))*integrate.quad(lambda x: cos_coef(i,x),-2*np.pi,2*np.pi,limit=200)[0]\\n        coefs[2*i+1] = 0\\n    coefs[0] = coefs[0]/2.0\\n    coefs[2*n] = -1.0*(K/(2.*np.pi))*integrate.quad(lambda x: cos_coef(n,x),-2*np.pi,2*np.pi,limit=200)[0]\\n    return np.array(coefs)\\n\\ndef fourier_y(c_list):\\n    # Given a list of coefficients c_list, f(x) is the fourier expansion. \\n    # f(x) = a0 + sum(an*cos(nx) + bn*sin(nx))\\n    # even entries in c_list are the a_n, odd entries are the b_n. \\n    #c_list should always have odd len(c_list).\\n    a = c_list[::2]\\n    b = c_list[1::2]\\n    def f(x):\\n        cosines = np.array([np.cos(x*n/(2.*K)) for n in range(len(a))])\\n        sines = np.array([np.sin(x*n/(2.*K)) for n in range(1,len(b)+1)])\\n        result = (np.dot(a,cosines)+np.dot(b,sines))\\n        # the result cannot have f(x) > 0 due to conservation of energy\\n        if np.isscalar(result):\\n            return result if result < 0 else -0.00001\\n        return np.array([i if i < 0 else -0.00001 for i in result])\\n    return f\\n\\ndef fourier_yp(c_list,x):\\n    #derivative of the function f(x)\\n    #f' = sum(n*bn*cos(nx)-n*an*sin(nx))\\n    a = c_list[::2]\\n    b = c_list[1::2]\\n    cosines = np.array([(n/(2.*K))*np.cos(x*n/(2.*K)) for n in range(1,len(b)+1)])\\n    sines = np.array([-(n/(2.*K))*np.sin(x*n/(2.*K)) for n in range(1,len(a))])\\n    result = np.dot(a[1:],sines)+np.dot(b,cosines)\\n    return result\\n\""
      ]
     },
     "execution_count": 228,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'''stuff for fourier series implementation\n",
    "\n",
    "def cycloid_coefficients(n):\n",
    "    # n is number of sine coefficients and n+1 the number of cosine coeffcieints. \n",
    "    # change of variables to theta is necessary since cannot write in terms of y(x).\n",
    "    # sine terms are all 0. cosine terms are non-zero.\n",
    "    coefs = [0]*(2*n+1)\n",
    "    def cos_coef(c,theta):\n",
    "        return ((1-np.cos(theta))**2)*(np.cos((c/2.)*(theta-np.sin(theta))))\n",
    "    def sin_coef(c,theta):\n",
    "        return ((1-np.cos(theta))**2)*(np.sin((c/2.)*(theta-np.sin(theta))))\n",
    "    for i in range(n):\n",
    "        coefs[2*i] = -1.0*(K/(2.*np.pi))*integrate.quad(lambda x: cos_coef(i,x),-2*np.pi,2*np.pi,limit=200)[0]\n",
    "        coefs[2*i+1] = 0\n",
    "    coefs[0] = coefs[0]/2.0\n",
    "    coefs[2*n] = -1.0*(K/(2.*np.pi))*integrate.quad(lambda x: cos_coef(n,x),-2*np.pi,2*np.pi,limit=200)[0]\n",
    "    return np.array(coefs)\n",
    "\n",
    "def fourier_y(c_list):\n",
    "    # Given a list of coefficients c_list, f(x) is the fourier expansion. \n",
    "    # f(x) = a0 + sum(an*cos(nx) + bn*sin(nx))\n",
    "    # even entries in c_list are the a_n, odd entries are the b_n. \n",
    "    #c_list should always have odd len(c_list).\n",
    "    a = c_list[::2]\n",
    "    b = c_list[1::2]\n",
    "    def f(x):\n",
    "        cosines = np.array([np.cos(x*n/(2.*K)) for n in range(len(a))])\n",
    "        sines = np.array([np.sin(x*n/(2.*K)) for n in range(1,len(b)+1)])\n",
    "        result = (np.dot(a,cosines)+np.dot(b,sines))\n",
    "        # the result cannot have f(x) > 0 due to conservation of energy\n",
    "        if np.isscalar(result):\n",
    "            return result if result < 0 else -0.00001\n",
    "        return np.array([i if i < 0 else -0.00001 for i in result])\n",
    "    return f\n",
    "\n",
    "def fourier_yp(c_list,x):\n",
    "    #derivative of the function f(x)\n",
    "    #f' = sum(n*bn*cos(nx)-n*an*sin(nx))\n",
    "    a = c_list[::2]\n",
    "    b = c_list[1::2]\n",
    "    cosines = np.array([(n/(2.*K))*np.cos(x*n/(2.*K)) for n in range(1,len(b)+1)])\n",
    "    sines = np.array([-(n/(2.*K))*np.sin(x*n/(2.*K)) for n in range(1,len(a))])\n",
    "    result = np.dot(a[1:],sines)+np.dot(b,cosines)\n",
    "    return result\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 229,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/anaconda3/envs/handsonml/lib/python3.6/site-packages/ipykernel_launcher.py:14: RuntimeWarning: invalid value encountered in sqrt\n",
      "  \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr4AAAFwCAYAAABejXgsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xd4lFX6xvHvSSdAQu8lVKVJkVUUBQVEUcGCq4ii6FrW3ta+q+66ll1/KoruurIqInYUFEVAl2LDLkqTjvROEkoIJDm/P56EhEhJYJI3M3N/ruu9ZvJOeZ8AxjtnznmO894jIiIiIhLpYoIuQERERESkPCj4ioiIiEhUUPAVERERkaig4CsiIiIiUUHBV0RERESigoKviIiIiEQFBV8RERERiQoKviIiIiISFRR8RURERCQqKPiKiIiISFRQ8BURERGRqBAXdAEVjXPOAQ2ArUHXIiIiIiL7VRVY7b33JX2Bgu9vNQBWBl2EiIiIiBxUI2BVSZ+s4PtbWwFWrFhBSkpK0LWIiIiISDGZmZk0btwYSvkJvYLvfqSkpCj4ioiIiEQQLW4TERERkaig4CsiIiIiUUHBV0RERESigoKviIiIiEQFBV8RERERiQoKviIiIiISFRR8RURERCQqKPiKiIiISFRQ8BURERGRqKCd20RERESiVeavsHAsbFsFlevDERdA1YZBV1VmFHxFREREos3uHfD5PfDjM+BzC89/dhf0fgaOuiq42sqQgq+IiIhINNm6Esb2hw0z7etGPaDu72D1l7BmBnx8NeTlQKdrg62zDCj4ioiIiESLravgrZMgfTFUqg2nvwJpp9pj3sOX98FXf4dpt0D9Y6Hu0YGWG2pa3CYiIiISDbI2w9snW+hNSYOLvikMvQDOwfF/g5ZnQ+4umHqzheEIouArIiIiEul8Hnw0BLYshKpN4IJpkJr22+c5B72egbgkWPU5LJtY3pWWKQVfERERkUj39cOwdIIF2rPfg5Sm+39u1YbQ6Xq7/80/yqe+cqLgKyIiIhLJln0MX9xn93v/C+p0OvhrOt8IOFg53UaJI4SCr4iIiEikylwOH14IeOhwBbS/rGSvS2kMzU6z+7NfLLPyypuCr4iIiEgkysmG8b+HnZugThfoNbx0r2//B7v95fWIWeSm4CsiIiISiabfBmu/gaTqMGCMze8tjWb97DWZv8LG2WVTYzlT8BURERGJNPNehZnP2v1+oyG1WenfIz4ZmvS2+0s+CF1tAVLwFREREYkkG2bB5Cvtfre/QPPTD/29mp9pt4vHH35dFYCCr4iIiEikyM6E8QMhJwuangLH3X9479csPzSv/Rp2bT38+gKmLYtFREQkMmRthk1zYdtK2LYadqyDnJ2Qt8s2cIivCokptlVvteaQ2sI2cYiJkDjkPUy6LH+TisZw+msQE3t475nSxHr+Zv4Kq2dAWt/Q1BqQCPmbFhERkaiSuwvWfQ8rpsOqz2DDTAu7pRWXDHWPhvrdoPFJ0PhkiK8U8nLLxdcPw8J3ISYe+r8NybVC874NT7Tgu+ozBV8RERGRcrFjgy2yWvQe/Pox5Oz47XNSmkJKGlRpAMl1bYFWTIJtxbtrK+zKhO1rIX0xZCy291j1mR3fPWZdDBr3gta/h9YDIaFquX+bh2TR+/DFn+1+72eh/rGhe+9GJ8K80bDys9C9Z0AUfEVERKTi2r0dFrwDc0fBiqk2ZaFApVrQqIcd9Y6Bmu1sKkNJ+TzYPB/WfA2rv4RlE2HrCtvad+kE+N+10PJsaHcZNO1j4bki2jQXPrrY7ne8Fo66MrTv3/BEu137tY20xyaE9v3LkfMR0pA4VJxzKUBGRkYGKSml+I9HREREQsN7G4GdPRIWvA27txU+VqezhdEWA6B2x9CGUe9h0xxYNA7mjoYt8wsfq9EGutwIbYdAfOXQXfNwbVsDb3SHjKX2C8B5n0BsfGiv4T08WwOy02HIjyXb8riMZWZmkpqaCpDqvc8s6esUfItR8BUREQnIrm0w9xWY+YyNYhao1gLaXgptLz60frSHwntY+y3MedlGmwvCd1J1OOpq6HIzVK5bPrXsz84t8GZP2DgLUpvD4K8guXbZXOutXjbi3vcF6HB52VyjFA41+EbkVAfn3LXA7UB9YA5ws/c+/CemiIiIRKL0xbbZwuwXITvDzsVXgSMHQbuh0OD48p9m4BzUP8aOEx+BOS/BD09DxhL45lH4YRh0uBK63g4pjcu3NrC2ZWP7W+itXA/O+7jsQi/YSPuKqbD+x7K7RjmIuODrnLsAGAZcC3wBXA185Jxr671fHmhxIiIiYryH5VMsQC75EMj/BLp6K+h0vQXe0szXLUuJKdDlJqtryQfwzSM2L/jH4fDTc9DuUvjdnVC9ZfnUs2MDvNvPulokpsLASdaerSzV7WK3638o2+uUsYib6uCc+xr4wXt/TZFz84Bx3vu7S/D68p3qMP8tqPe78vvoRkREJEg52fDL6/DDk7Dh58LzzfpB5xsg7VRwFXx/Le9h+f/g64dgxTQ752LgiEFw7N1Qq33ZXTt9CYw9Azb/Yv2IB04sDKVladNcGNnO5jdfn3H4/YEPk6Y6AM65BOBo4NFiD00Gjt/PaxKBxCKnyq9vyTf/hM/uhAbd4YJpkdNAW0REpLgdG+Hn5+DHZ2xjCbAQ1e4yC7w1WgdbX2k4Z10emvaBVV9aAF46AX55zY5m/eDo26BJr9BO0Vg8Hj66xBaZVW1s0xtqHBG69z+Q6kdYz+Pd2yF9UfldN8Qq+K9UpVYLiAXWFTu/Dqi3n9fcDWQUOVaWWXXFHfF76w+4+gv4+pFyu6yIiEi52fQLfPxHGNEYvviLhd4qDeHEf8BVK6D38PAKvcU1PB7O/RAu/h5aDQQcLP0IxvSBV7rAnFGwO+vwrrFjI0wcCuMGWOit3w0unFG+4TMmFmq2sftFFx6GmYia6uCcawCsAo733s8ocv5eYIj3/sh9vGZfI74ry22qw9zR8NEQcLEw6HNo0K3srykiEi12bbPNCravtcC1KxN277BNC3KyIDfbfv66WPvULSbORiGTqkNidUisBpVqWlBLqBL0dxM+fJ5tMPHj8Pz5u/nqdIGut9nmEKFuuVVRbFkI3w+zxXA5+YE3IcUW6rW5GBocV/JPeLevhZn/tmkhu7YCDo6+GU58NJheuhOG2EYWJzwEx95T/tcvQlMdzEYgl9+O7tbht6PAAHjvs4Hsgq9dea8abXNR/scjr1vz6SE/hs8uMSIiFcGubbDhJ5vzuGUhbFkA6Qutr+nu7aG7TmI1+3i5aiOo2gSqt7YRsBpH2tcBz3msELI2weyXbEpD+uL8k8567na91TZCqKibQIRK9VbQ51no/jdb+DZrhG33+/PzdiRVh6Z9bcONOp3t31Ol2oU7y2UshXXf2S8OyyaDz7X3rd0J+vzLgnNQ9oz4zguuhsMUUSO+sGdx2/fe+2uLnJsLvFchF7cB7EyHUR1h63Kb63Tai+VzXRGRcJO721ayr5lht+u+t523OMD/y+KSrd1Tcl0LHfHJdq5gK1ufBz4H8vKP3dusP2p2ut1mbbSR4gOJS7IgXKOthZmCI7lWSL/9CsnnwaovYNZ/Yf6bNooONsrZ7lKbv1u9VbA1BsnnwcpPYc5Im6O7c3PpXl+/Gxx9q22fHPSiv4Vj4f1zoe7RcPF3gZaiDSzy5bczewX4IzADuAq4Emjnvf+1BK8PZgOLlZ/CmycBHvq/Da3PK79ri4hUVLm7Ye03sHI6rJhuayL2NYpbpQHUbG8Bq3pru01tYedDMUUhOxO2rbTtbDNXQOYyG1ne/Ivd5mbv+3VVGhWG4Lpd7KP+qo0iY9Rz83z72HvuaPvzKFCnM3S8BtoMrlg7nFUEeTnWBq2gH+6Gn206Q8HmGC4GUtLsU4TGJ0PzMwpHWSuCTb/AyDb293pDZqBBXMG3iPwNLO7ANrCYDdzivf+0hK8Nbue2z++Frx+2EYlLfrYfjiIi0Wb7OlsctPRDWDYpf25jEUk1oOEJULerjTzV7WIjukHJy7WPsjfPgw2zLNCs/8FWvu9LpdqFIbjgNrVZxQ/D3tuipsXv2cjfuiIjfglVodV50PFq+wi/on8vFU3OTsDZvN2K/GeXuxueTrYAf+WvkNIksFIUfEMk0OCbuxteP95+mDTpZW1Kgv5YQ0SkPGxeYB+TL/nAtoktOnWhUi1o1NOOxj2tR2o4/GzMzrS5x+t/LAzDm+ZaaCgusVphCK7TxQJ99ZbBf58702HVZ7bRxJLxRebtYgsC006FtkNsDm98cnB1Svl5qa39kjdwov39B0SL2yJBbDyc/iq80tl+yHz3BPzuT0FXJSJSNnash1/etI/L136z92N1utjHvM3PhHpdgw+AhyIxBRqdaEeBnJ2wcbbNTV7/A6z7ATb+bPOJl0+xo0B8lcIpEjXa5E/haG3TN8piVDB3twWagqC+8rP87WmL/BISmwhNekPLs6DFWVC5bujrkIqteiv7d1L0l6AwouBb0dRoDScPg4+vgs/vsR8wdTsHXZWISGjs3gGL3rOwu2xS4Yp1FwtNT7E+qM1Pt3AXieKSLMjX61p4LneXjQSv+yE/DH9vI8W7t9lo66rP9n6P+MpQraV1kqjSwI7KDSC5joXthKoQX9WuhbcpCngL3dnp+Yv20mH7apuikbEsf87yfKuluOqtofFJ9veTdpraukW7ai3sNn1JsHUcIgXfiqjDFdbibNE4mDDYmmLrIyQRCVd5uTaSOW80LHy3cCEP2JbtbS6GIy6I3tHD2ASo08kOLrdzeTm2eGz9D7B+poXSLQssbOzebsF4w0+hryUhxeqo3QnqHwONToKqDUN/HQlfqc3tNkPBV0LFOThlhK383PwLTL/degKKiIQL7y2YzR1tW7huX1P4WGozC7ttLgrbbU/LXEwc1GpnR9shhedzd1uf1/SFsG0VbFudf6zKb7u21Vqv7dqavzmHA/KPuESbS1xwJNexDgKpaXZbvXV4LLKTYCn4SplIrgWnvQzv9IWf/mX7frc4M+iqREQOLHM5zHvNRnc3zSk8n1TDRnULdq5SuDo0sfE2JS6ct/iV8FYQfNMX2y+4YfbfsoJvRZZ2Chx9C3z/JEy6HC6dFb0fBYpIxbUzHRa+Y2F3xbTC87GJ0KK/hd1m/YLZYlVEQis1DXA2ZSlrIyTXDrqiUlHwrehOeBh+/QQ2zoJJl8E5H4bdb1ciEoFyd1mv3XmjbTeqohs4ND7Jwm6rgZBULbASRaQMxCVBlYa2oUvGEgVfCbG4JDjjNRjd1f4nM/NZ6Hx90FWJSDTyHlZ/aWF3/lt7b71aM38u6pEXBtrUXkTKQbXmFnzTl0D9Y4OuplQUfMNBrfbQ4zGYeiN8erttY1irXdBViexfXo61SUpfZLc71tuxfR1kb4GcrMIjL9cW8sTEQUy83Sak2A6GSdVtEU6l2raTYZVG+bcN7HlSPjbOsQVqv7xuC6sKVGkARw620d3aR+nTKJFokdoMVn4KmUsP/twKRv/nCBedr7cWZ8smwoeDYPA3EF8p6KpELMzu1Yx/lvUE3dfuVKHiYqxvaY38hv7Vj7DuANVb2+r0mNiyu3a0yFhmQfeX1+3vtEB8FWg9ENoMsSkN+rMWiT5VG9vt1hXB1nEIFHzDhXNw2kgY1dF2/Zl2M5zyn6Crkmi0bQ2smGrH8in7b2kTl2RN9lPSoHI9a52UXAcSq1tf6rhKdrhY28QgL8eO3F3Wjil7C+zMP3ast4/Vtq60tk15u+3rbSv33ukKIC7ZPiWp3dGOOp1sNDKhapn/0YS9zF+tf/gvb8KaGYXnYxMgrR+0GWw7qamvuEh0U/CVclG5Lpw+Gsb0hZ+ftykPRw4KuiqJdAX9WBe9Z6Fow8xiT3BQ40jbVrXu0RY292yrWgbbzPo8C8IZS62h/5YF1uh/ywLrbZqzw7a/Lb4Fbmrzwsb8dbvYVrCV60f3x/PeW8uxhWPt73b9D0UedNCkl83ZbXWuTTsREYGwDr7Oe3/wZ0UR51wKkJGRkUFKSkrQ5ezbF3+Br/5uI1gX/wDVWwZdkUSiTXNhziiY/6ZNXdjDWWhsfLIFo4Yn2DapFUFers0rLtjVasNPtuvVtlX7fn5yHfteih7VWpRNYK8odm2DldNh2WSbPpW+qPAxF2N/ny3PgSPOj9xtg0Xk8GycDS93sF+Ir9t88OeXgczMTFJTUwFSvfeZJX2dgm8xYRF883LgrV62f3udLnDhl7Yjj8jh2rnFdtqa+7LN2y0QVwnSToUWZ9lH3cm1gqvxUGRtKgzB63+0Y/M8Gz0uLqFq/hSJLoVhuGZb2zggHOVk2/e7Yhr8OhlWfW5TRQrEJkLTU6Dl2dBiQNi1JhKRAGRnwDP5rQpv3Abxlcu9BAXfEAmL4As213FUJ9i5CbrcBCcPC7oiCWcbZ8OPw2HuK9ZpAaxrQrPToe0ltvlApM3r3J1li7bW/1AYhjfOgpydv31ubALUbF9sdLhjID/sD8h727523few+gtrPbb227177ILNu047FdL6WujV/GcRKa3hKbY19tB5UPPIcr+8gm+IhE3wBVjyIYzN38Z4wFhodXaw9Uh48Xm28cAPT9lCtQK1j4L2V9j88Wgb/cvLgc2/FAbhgiM7Yx9PdtbLslpLqNYKqrfKv9/S5r+VZdcVnwdbV9nCwvTFFtg3/gzrf7JfhourVAsadLeQm9bXaozmuc0icvhGtrMpcQMn206z5exQg68Wt4Wz5mfA0bfB94/D5MttwY4ax8vB5OXCgrfh64dspBess0Krc6DzDdDwxOgNRTFx1hGiVnvbjAFsFDVj6W/D8PY1FjrTFwOTfvteSdWt5VqV/COxOiSm2pGQfxubYH/2MbH53S18kR7HO2D3DtskoqAPctYGG9HNXGrdL/bFxdhiw/rHQcPuFnirt4rev1MRKRtVG1vwDbMFbgq+4e7Eh22u79pv4IML4ILp9j9TkeLycmDeq/D1w9YBAWyjiI7XQKfrIKVxsPVVVK5gZLe59a8tsH2djQ5vWWgLxNLzb7csstBa0Ipt05yyqSsmzqYspDazOci1O9pofY226vEtImUvTDs7KPiGu9gEOPMNeKULrPkKpt0GvYcHXZVUJN5bq6rP7oYt8+1cUg3ocrON8CZVC7a+cFW5rh2Ne+593nubGrF9tY3Obltlt9npdn5Xpt1mZ9giM59ro/A+115ftMdxXCX7uyrogVypjvVErtYcqjTU7nUiEpwqjex228pg6ygl/dSMBKnNoN8rMK4/zHwGGnSDNhcFXZVUBKu+gE/vsEVOYHM9u94Ona7Rgqay4pz9MpFUzUZiRUQiUZX6drt9bbB1lFIEN6uMMi3OhG5/sfuTr4QNPwdbjwRr8wJ47xx44wQLvXGVoNuf4Q+L4Zg7FHpFROTwVC4IvmuCraOUFHwjyXH3Q9O+tjDm/YGwMz3oiqS87dpmUxpebm/TG1wMdLgS/rAIuj9YcTaaEBGR8KbgK4GLiYUzXoOqTWyRzcSh+27QL5HHe5j/FrzUBr551OaONusHl86Gvs9rBy4REQmtyvXsdvu6sMoaCr6RplJNGDDGFr0tfg+++UfQFUlZ2zQPxvSxrh7bVtqc77Pfh3M+hJptgq5OREQiUXJdwNnC3KyNQVdTYgq+kaje76DXM3b/83th8QfB1iNlI3cXzPgbvNIJlk+BuCQ4/q9w6Rxo0V99W0VEpOzExtuCaYBt4TPdQcE3UnW4Ao66CvAwYbA1mZbIseZrGH00fHm/BeDmZ8DQuXDcferhKiIi5aNK+M3zVfCNVM5Br+HQqIftpT1uAGTtYytTCS+7t8PUW+C142zXtUq14YzX4ezxNsVBRESkvFQOv5ZmCr6RLDYB+r9juzulL4YPzofc3UFXJYdq2ccwsj38MAzwtqXu0Llw5CBNaxARkfIXhp0dFHwjXXItW+gUX8XmgU67JeiKpLR2b4dProV3+kLmMuvace5H0G+U/f2KiIgEYU9nBwVfqUhqd4DTRwMOZj4LPz4TdEVSUqu+hFEd4ad/29edroOhs6HZacHWJSIiohFfqbBangUnPGz3p9wIC8cFW48cWE42fHoXvHmiTVOp0gjO+xh6P6Nd10REpGJIrmu3O9YHW0cpKPhGk2PuLNLp4UJY/VXQFcm+rP8JXv0dfPsPawre9hK4dBY07RN0ZSIiIoWSa9vtjg3B1lEKCr7RxDno/ay1vsrZCeP6w5ZFQVclBfJy4OtHLPRunGUdGwaMhX4vQ1K1oKsTERHZW3Idu9WIr1RYMXFwxhtQt6vttPLuaWH1m1rE2rIQ3ugBn99j2w23PNvm8rY6O+jKRERE9q1S/ojvzs2QlxtsLSWk4BuNEqrAOR9Y39f0xTD2DOv1K+XPe/jpORjVCdbMgIQUOO1lGPBu4W/SIiIiFVGlmvl3POwMj70CFHyjVeW61hIrqSas/dY2uNidFXRV0WXbGvul45NrIGcHNOllc3nbXaK+vCIiUvHFxFmOgLCZ7qDgG81qHAHnTbIuASumwfiBtv2tlL2F78LLHWDpRxCbCCcPs64NKU2CrkxERKTkwmyBW9gEX+fcvc65L51zO5xz6ft5ThPn3Hjn3Hbn3Ebn3NPOuYTyrjWs1D0azvkQ4ipZCJtwkS2ykrKRnQETh8L7A+1joTqdYcgP0OUmcGHzn6OIiIgpmOebpeAbagnA28C/9/Wgcy4W+BCoDJwADAIGAo+XV4Fhq9GJcNY42+J4wRiYdHnYTFIPKyum22YUc162kHvsPTD4K6jZNujKREREDk2YdXaIC7qAkvLe3w/gnBu6n6f0BdoCjb33q/Ofexsw0jl3r/c+s1wKDVdpfeGMN2H8eTD3FRv17TfK5u/I4cnJhi/+At/9H+Ahtbn92TbsHnRlIiIih0dTHQJzHDC7IPTmmwQkAkfv70XOuUTnXErBAUTvtlitzob+b0FMPPzyOnwwSHN+D9eGn60v73ePAR46XAGXzFToFRGRyKCpDoGpB6wresJ7vwXYlf/Y/twNZBQ5VpZVgWGh1bnWSis2ARa+Y3NRc3YGXVX4ycuBb/6592YUZ70HfUdoy2EREYkcBcE3TKY6BBp8nXMPOOf8QY6upXhLv6/L7Od8gUeA1CJHo1JcLzK1OBPOfh/ikmDJBzC2P2RrpkiJbZoLr3eHz+60EfPm/a1NWcsBQVcmIiISWgVzfMNkxDfoCZzPAG8c5DnLSvhea4Fji55wzlUH4ik2ElyU9z4byC7ymhJeLgSmT4d77oHkZKhUCSpXhnr1oGFDaNIE2rWD1q0hPr78aiqQdiqcM8G2NV7+CbzZE879EKo0KP9awkVeDnz7fzDjfgu8ialw0pPQbqj68oqISGQKszm+gQZf7/1GYGOI3m4GcK9zrr73fk3+ub5YqP0+RNcIrTVr4MsvD/ychATo0AF69ICePe2oVq186mtyMpw/zTZZ2DATXjsOBk6Emm3K5/rhZONsmHgZrPvOvm52OpzyPFRtGGxdIiIiZalgA4sw2bnNeX+gWQAVh3OuCVADGADcDpyY/9Ai7/22/HZmM7HR3dvznzsSGOe9v6EU10kBMjIyMkhJSQnhd7APK1fCN99AVpYd27ZZGF65EpYuhVmz7FxRcXFw0klwzjkwcCDUrVu2NQKkL4F3T4MtCyGpOpz5FjTtU/bXDQe5u+Dbf8JXD+aP8laDk5+CtkM0yisiIpFv60p4vrF1gbp5V7n9vy8zM5PU1FSA1NJ07gqn4DsSuHQfD53svZ+W/5wmwL+AXkAW8Brwp/zpDCW9TvkF34PJy4Nff4Wvv7ZpEVOnwvz5hY/HxUH//nDlldC3L8TGll0tOzbatIc1X1kP2p7/B11uju5wt2K6bTe8eZ593fxMOOU/mg4iIiLRY3cWPJ1s969Pt2l+5SDig295qVDBd18WLYKxY2HMGBstLtCkCdx6K1xxhc0VLgs5O+GTP9oGDGCjmn3+A/GVyuZ6FdWOjfDpnwr/HJLrwElPwJGDo/sXARERiU5PVYacHfCHxVCteblc8lCDbyS1M4sOLVvC7bfbKPCsWXDTTVCjBixfDjffDE2bwt/+Bun73NX58MQlwakvwcnDwMXaRhevHWtdDKJBXi78/F946Yj80Oug4x/hsl+gzUUKvSIiEp0qhc88XwXfcNa+PQwbBqtWwXPPQfPmsGkT3H8/tGhhj2WXeJZHyTgHXW6C8ybbSOfGWTD6aPjpPxDJnx4snwqju8LHV8LOzVD7KLjwS+jzb5v3LCIiEq0q1bLbrFD1Kyg7Cr6RICkJrr7a5v++/jq0aQObN8Mtt9j9t98OfSht0gsu+dnanhVMgXjvHNi2+uCvDSdbFsK4s+HtXtbZIjEVej4OF30HDboFXZ2IiEjwCjo7ZGnEV8pTXBwMGgQ//wwjRkD9+tYd4vzz4fTTYcmS0F6vcl04d4ItdIuJh8Xvwci28PPz4PNCe63ytnUVfHIdjGxn35eLhU7XweWLoOutEBtAb2UREZGKqFItSEix7kYVnBa3FVPhF7eVxvbt8Nhj8OijNuUhKQn+8hebIxzqTTE2/AyTr4C139rXjXpYW686nUJ7nbK2fS188w/46d+Qmz9NJO00OOlxqNk22NpEREQqIu/LfZ2LujqESEQF3wILFsA118CUKfZ1167wyitw5JGhvU5eLvw4HL74M+zeDjjr/ND975DSOLTXCrX0xfD9kzD7RcjJsnMNT4DuD0LjkwItTURERPam4BsiERl8wX4bGz3aukBs2WKjv//8J1x3HcSEeMZL5q/w2d3wy+v2dVwSdLwGjr4VqjYK7bUOh/ew+kv47nFYNA7I/2+h/rFw/IO2SYc6NYiIiFQ4Cr4hErHBt8CqVXD55TB5sn192mk2+lurVuivtfZbmP4nWPmpfR0TbyPAnW8IdgrEjg0wbzTMegE2zSk836wfHH2bLdxT4BUREamwFHxDJOKDL9hI57//DX/6k22V3LixdX449tiyudaySfDNI4UBGKDe76DDFdBqYGH/v7KUtRmWjIcF78CyiZC3287HJcGRF9mCNc3hFRERCQsKviESFcG3wKxZMHAgLFxoi90efxyuv77Isq7iAAAgAElEQVTsRjtXz7B5tIvGFQZPFwuNe0KLs6DxyVCrnW2JfLhyd9mI84qpsHyKhW6fW/h43a7Q/nI48kJIqnb41xMREZFyo+AbIlEVfAEyM+EPf7AtkAEGD4YXXrA5wGVlx3qYMwrmvWq9cYtKqmlzbGt1gFrtIbU5VGkAletDXOLez83LhewM2LbKji0LYePP1mFi46zCRWoFah8FLc+F1gPtvUVERCQsKfiGSNQFX7DpCMOHw223QU6OTXkYNw7q1Sv7a6cvgYXvwq+TYdUXttf3/sTEQWyizRXO3WkbZxxIpVrWkaHxydC0L1RvGdLSRUREJBgKviESlcG3wNSpNvVhyxab9zt+PHTsWH7Xz90N676H9T/CxtmwaTZsXWG7weUeYOvlpJpQtSGkpNmobq0OUOsoqHGEFqmJiIhEIAXfEInq4As23/fMM633b+XKtgVy//7B1uQ9ZKfD7h020pu7G+IrQVxlSKhiC9REREQkahxq8NWWxbK3Vq3gq6+gTx/b+e2cc+DFF4OtyTlIqm6jutVaQM0jIaUpJNdS6BUREZESU/CV36peHSZMgKFDITfXFr89+qiNvIqIiIiEKQVf2bf4eBvpvfNO+/ruu+HWWyEvL9i6RERERA6Rgq/sn3M20vvEE/b1sGEwZAjs3h1sXSIiIiKHQMFXDu6WW2D0aIiLg9deg/PPh127gq5KREREpFQUfKVkLrrIevsmJtrtuefCzoP00RURERGpQBR8peTOOAPef992dfvwQzj7bMjKOvjrRERERCoABV8pnb59reNDcjJMmmQ9f7dvD7oqERERkYNS8JXSO/lkmDgRqlSBKVMs/GrkV0RERCo4BV85NCeeCJMnQ9WqMG2azfnNPsC2wiIiIiIBU/CVQ3fccYXTHiZOhEGD1OpMREREKiwFXzk8J5wA771X2O3hkktstzcRERGRCkbBVw5fnz7wzju229sbb8AVV2iHNxEREalwFHwlNM44A15/HWJjYeRIuPFG8D7oqkRERET2UPCV0Bk4EEaNsq2On30WHnww6IpERERE9lDwldAaPBieecbu338/PPdcsPWIiIiI5FPwldC79lq4777C+2PGBFuPiIiICAq+UlYeeAD++Eeb53vRRbbRhYiIiEiAFHylbDhnUx7OOw927YKzzoLvvw+6KhEREYliCr5SdmJjYfRo2+J42zbo1w8WLQq6KhEREYlSCr5Stgo2tujcGTZssPC7cWPQVYmIiEgUUvCVspeSYlsbp6XZiO9ZZ0FWVtBViYiISJRR8JXyUa+ehd9q1eDLL2HIEO3uJiIiIuVKwVfKT5s2Nu0hIcG2OL7jjqArEhERkSgSFsHXOZfmnHvBObfUOZflnFvsnPurcy6h2PM6OOem5z9nlXPuPuecC6pu2YeePeGll+z+448XbnYhIiIiUsbCIvgCR2K1Xg20A24B/gg8XPAE51wK8DGwGvgdcAPwJ+DW8i5WDmLwYHjoIbt/003w/vvB1iMiIiJRwXnvg67hkDjnbgeu8d43z//6GuARoK73Pjv/3F1YAG7kS/iN5gfojIyMDFJSUsqmeLGNLa6+GkaMgEqVYPp0+N3vgq5KREREwkBmZiapqakAqd77zJK+LlxGfPclFdhc5OvjgOkFoTffJKABkLa/N3HOJTrnUgoOoGpZFCvFOAf/+hecdpp1eDjzTFi2LOiqREREJIKFZfB1zrXARnKfK3K6HrCu2FPXFXlsf+4GMoocK0NUphxMXBy89RZ06gTr11v4zSzxL20iIiIipRJo8HXOPeCc8wc5uhZ7TQNgIvC29/6/xd6y+HQGt5/zRT2CjR4XHI0O/TuSUqtaFcaPh/r1Yc4cGDQIcnKCrkpEREQiUNAjvs8AbQ5yzC54cn7onQrMAK4q9l5r+e3Ibp382+IjwXt477O995kFB7D1kL8bOTSNGtkCt0qV4KOP4Lbbgq5IREREIlBckBf33m8ESrR/rXOuIRZ6vwcu894X3/1gBvCwcy7Be78r/1xfrMvDstBULGWma1d45RU47zx4+mk44gi49tqgqxIREZEIEvSIb4nkj/ROA1ZgLcpqO+fqOeeKjvC+BmQDI51z7Z1z5wD3AE+UtKODBGzgQHg4v0PdjTfC5MnB1iMiIiIRJSyCLzZy2xLohS0+W1PkAMB7nwGcgs3R/Q74F/BE/iHh4q674JJLIDcXfv97mDs36IpEREQkQoRtH9+yoj6+FUB2NvTpA59/Ds2awddfQ+3aQVclIiIiFUQ09vGVSJWYCGPHQvPmsHQpnHuuhWERERGRw6DgKxVTrVrwwQeQmmojv1deabu9iYiIiBwiBV+puNq0sQ0uYmOt48MjjwRdkYiIiIQxBV+p2Pr2heHD7f6998KYMcHWIyIiImFLwVcqvmuusfZmAEOGwLffBluPiIiIhCUFXwkPTzwBp58OO3fCgAGwYkXQFYmIiEiYUfCV8BAbC6+/Du3bw9q1Fn63bQu6KhEREQkjCr4SPlJSYPx4qFMHZs6Eiy6yjS5ERERESkDBV8JLWhqMG2e9ft9/33Z6ExERESmBUgdf59xI51yPsihGpESOOw5GjrT7//d/8N//BlqOiIiIhIdDGfGtCkx2zi10zt3jnGsY6qJEDmrQIHjgAbt/zTUwZUqg5YiIiEjFV+rg670fCDQEngF+Dyxzzn3knDvPORcf6gJF9uu+++DCCyEnBwYOhAULgq5IREREKrBDmuPrvd/kvX/Ke98ZOAZYBLwCrHbOPemcaxXKIkX2yTl48UWb+pCeDmecAZs2BV2ViIiIVFCHtbjNOVcf6Jt/5AITgHbAXOfcLYdfnshBJCXB2LHQtCksWmQjv7t2BV2ViIiIVECHsrgt3jk30Dn3AfArNt3hSaC+9/5S731fYAhwX2hLFdmPunXhgw+galWYPt3m/HofdFUiIiJSwRzKiO8aYAQWeo/x3nf13j/nvd9a5DmTgPRQFChSIu3bw5tvQkyMTX947LGgKxIREZEK5lCC7y1AA+/9dd77mft6gvd+i/e+2eGVJlJK/frBsGF2/667bAqEiIiISL5D6erwivd+Z1kUI3LYbrgBrrvOpjpcfDH88EPQFYmIiEgFoZ3bJPIMGwZ9+8KOHdC/P6xaFXRFIiIiUgEo+ErkiYuDt96Ctm1h9WoYMAC2bw+6KhEREQmYgq9EptRU6/RQq5ZNdxgyBPLygq5KREREAqTgK5GrWTMYNw4SEmyh2513Bl2RiIiIBEjBVyJb9+7W3gzg//4Phg8Pth4REREJjIKvRL6LLoKHH7b7N92kNmciIiJRSsFXosNdd8HVV1ubs8GDYcaMoCsSERGRcqbgK9HBOXjmGTjzTNi509qcLVgQdFUiIiJSjhR8JXrExcEbb0DXrrBpk+30tn590FWJiIhIOVHwlehSubK1OWvWDJYssZFf9fgVERGJCgq+En3q1oWPPoIaNeCbb+DCCyEnJ+iqREREpIwp+Ep0OuIIeP99SEyE8ePh+utt4ZuIiIhELAVfiV7du8Orr9rCt//8B+6/P+iKREREpAwp+Ep0GzgQnn3W7j/4IDz1VLD1iIiISJlR8BW55hoLvQA33wyjRwdbj4iIiJQJBV8RgHvvtV3dAIYOhQ8/DLQcERERCT0FXxGweb5PPAFDhkBuLpx3Hnz2WdBViYiISAgp+IoUiImBF17Ye3e3n34KuioREREJEQVfkaLi4+Gtt+DEEyEjA/r2hV9+CboqERERCYGwCb7Oufedc8udczudc2ucc6845xoUe04H59x051yWc26Vc+4+55wLqmYJU5UqWY/fTp1sS+PevWHRoqCrEhERkcMUNsEXmAqcDxwBDARaAGMKHnTOpQAfA6uB3wE3AH8Cbi33SiX8VasGH38M7drB6tUWfn/9NeiqRERE5DA4H6a7VTnnBgDjgETv/W7n3DXAI0Bd7312/nPuwgJwI1/CbzQ/QGdkZGSQkpJSRtVL2Fi3Dnr2hPnzoVkz+PRTaNQo6KpERESiWmZmJqmpqQCp3vvMkr4unEZ893DO1QAuAr703u/OP30cML0g9OabBDQA0g7wXonOuZSCA6haRmVLOKpbF/73P2jRApYuhV69YM2aoKsSERGRQxBWwdc59w/n3HZgE9AEOKvIw/WAdcVesq7IY/tzN5BR5FgZmmolYjRsCFOmQNOmsHAh9Oljc39FREQkrAQafJ1zDzjn/EGOrkVe8hjQGegL5AKjii1eKz6dwe3nfFGPAKlFDn2OLb/VpImF34YNYe5cG/lduzboqkRERKQU4gK+/jPAGwd5zrKCO977jcBGYIFzbh6wAugGzADW8tuR3Tr5t8VHgvfInxqxZ3qEmkDIfjVvbuH35JNhzhw46SSbBtGwYdCViYiISAkEGnyLBNlDUZBQE/NvZwAPO+cSvPe78s/1xbo8LDvkIkWKat3aFrj16mUL3nr2tDDcpEnQlYmIiMhBhMUcX+fcMc65651znZxzTZ1zJwOvAYuxwEv+19nASOdce+fcOcA9wBMl7eggUiItWlj4bdYMFi+GHj1gyZKgqxIREZGDCIvgC2QB5wL/A+YDLwKzgZ4FXRy89xnAKdgc3e+AfwFP5B8iodW0qYXfVq2sv2/PnrBgQdBViYiIyAGEbR/fsqI+vlIqa9ZYl4e5c6FOHfjoI+jSJeiqREREIlpU9fEVqTDq14epU6FzZ2txdtJJ9rWIiIhUOAq+IoerTh2YNs26PWzdCqedBu++G3RVIiIiUoyCr0gopKTAhAlw7rmwaxf8/vcwYkTQVYmIiEgRCr4ioZKUBG+9BVddBXl5dvvgg6B59CIiIhWCgq9IKMXGwnPPwZ//bF/fdx8MHQrZ2Qd8mYiIiJQ9BV+RUHPORnqfe86C8KhR0LcvbNoUdGUiIiJRTcFXpKxcfbXN+01JsZ6/xx0HCxcGXZWIiEjUUvAVKUt9+8IXX9iGFwsXQrduMH160FWJiIhEJQVfkbLWvj18/TUccwxs3mwbXgwfrkVvIiIi5UzBV6Q81K1rvX4vvBBycuDGG+HSS2HHjqArExERiRoKviLlpVIlePVVeOIJW/T2yivQvTssXRp0ZSIiIlFBwVekPDkHt9wCn3wCtWvDzJnQtStMmhR0ZSIiIhFPwVckCCedBD/8UDjvt18/uPtu2L076MpEREQiloKvSFAaNbIOD3/8oy10e/RR6NFDUx9ERETKiIKvSJCSkuDf/4YxY6BaNfjqK+jUybY+FhERkZBS8BWpCAYOtPm+xx8PmZlwwQXW9SE9PejKREREIoaCr0hF0bSpTX24916IibGtjtu3h4kTg65MREQkIij4ilQkcXHw97/D559Dq1awapUtfLvqKti6NejqREREwpqCr0hFdNxxNvXhxhvt6xEjoF07GD8+2LpERETCmIKvSEWVnAxPPWU7vjVrBitWwIABcO65sHJl0NWJiIiEHQVfkYquZ0+YPRvuvNOmQowdC23awLBhtv2xiIiIlIiCr0g4SE62Pr8//GDTILZtsx3gOnWCyZODrk5ERCQsKPiKhJMOHWzh23/+AzVqwJw5cOqpcMYZ8MsvQVcnIiJSoSn4ioSbmBjr8rBoEdx8s01/mDDBWp/deCNs3Bh0hSIiIhWSgq9IuKpeHZ580kZ9BwyA3FwYPtwWwv35z7BlS9AVioiIVCgKviLhrnVreO89+OQT6NLF5v8+9JAF4AcftJ3gRERERMFXJGL07g3ffQfvvmvTHjIy4L77oHlzC8IaARYRkSin4CsSSZyDc86Bn36C11+30eBNm2zqQ5MmcNtt6gEsIiJRS8FXJBLFxMCgQTb/d/Ro6waxbRs88YSNAF92GcyaFXSVIiIi5UrBVySSxcXBRRfZCPCECbYZxu7dMHIkHHWUff3WW3ZOREQkwin4ikQD56BfP9v++Kuv4LzzIDYWPv0ULrgAmjaFv/4V1qwJulIREZEy47z3QddQoTjnUoCMjIwMUlJSgi5HpOysXAnPP2/HunV2LjYWTj/dpkKccQYkJARbo4iIyD5kZmaSmpoKkOq9L3H7IgXfYhR8Jers2mWdIJ55Br74ovB8zZo2TeKyy2xrZBERkQpCwTdEFHwlqs2bZ/N/X3ll72kPHTvC4MFw/vmQlhZUdSIiIoCCb8go+IoAOTnw8ccWgseNs1HhAsccY/OCf/97aNw4sBJFRCR6KfiGiIKvSDGbN8OYMfDmm7Y4Li+v8LHu3W0U+OyzrU+wiIhIOVDwDREFX5EDWLsW3nnHQvDnn0PRnx+dO8NZZ9nRsaN1khARESkDCr4houArUkKrVtlI8Dvv2KK4oiPBTZsWhuATT4T4+ODqFBGRiBM1wdc5lwh8DXQEOnvvZxZ5rAPwDHAMsBn4D/CgL8U3qeArcgg2bIAPP7T5wJMnQ1ZW4WPVqlmLtAED4LTTwH5QiYiIHLJoCr5PAa2AfhQJvvmBdQEwFXgIaA2MBP7qvX+8FO+v4CtyOHbsgE8+gffeg/ffh40bCx+Li4MePaB/fztatAiuThERCVtREXydc/2AJ4CBwBz2Dr7XAI8Adb332fnn7gJuABrtb9Q3fwQ5scipqsBKBV+REMjNhRkzYPx4O+bN2/vxtm0LQ3C3braBhoiIyEFEfPB1ztUFvgfOBjYCS9k7+I7CvvmzirymM/AD0Nx7v3Q/7/sAcH/x8wq+ImVg0aLCEPzppxaMC9SqZbvF9e8PfftC1arB1SkiIhVaRAdf55wDJgBfeO//7pxL47fBdzKwzHt/VZHXNQBWAcd772fs57014isShPR0mDjRQvCECfZ1gYQEOOkkmxfcv79apYmIyF4ONfjGlF1JB+ece8A55w9ydMWmK6RgUxkOpHiKd/s5X/gC77O995kFB7D1kL8hESm5atVg0CB49VVYvx6mToVbb4WWLW3DjMmT4frrrUNEx47w5z/D11/v3T1CRESkFAId8XXO1QJqHeRpy4A3gP7sHWBjgVzgVe/9pYc61WEfNWlxm0iQvIf58wunRBRvlVa3Lpx5po0E9+kDlSsHV6uIiAQi0qc6NMFGfAs0ACYB5wFfe+9X5i9uexhb3LYr/3V3AjdygMVt+7iWgq9IRbJpE3z0kXWImDgRthb5UCYpCXr3thB85pnQsGFwdYqISLmJ6OBb3H7m+KYC84EpWABuhbUz+5vamYlEiF27bFFcwWjw0mIf5HTpYiF4wADbSU67x4mIRKSoD7755zsAz2IbWGwBnsOCrzawEIk03sOcOYUh+Kuv9t5CuXFjOP98m0d89NEKwSIiESSqgm9ZUvAVCVPr19vucePH28K47dsLH2vRwgLwoEHQvn1wNYqISEgo+IaIgq9IBNi5EyZNgjfesLnBO3YUPtaunQXgCy6AVq2Cq1FERA6Zgm+IKPiKRJjt2+GDDywET5hg84QLHHssXHqpBeHq1YOrUURESkXBN0QUfEUiWHo6vPeeheCPPy7cOS4hAc46C4YOtV3j4uICLVNERA5MwTdEFHxFosS6dfDaazByJPz8c+H5evXg4ottJFjzgUVEKiQF3xBR8BWJQjNnWgB+9VXYuLHwfLducNVVNh84OTmw8kREZG8KviGi4CsSxXbtss0yXn7ZukPk5Nj51FQbBb76aujQIdgaRUREwTdUFHxFBLCpECNHwvPPw5Ilhee7dbMAfP75GgUWEQmIgm+IKPiKyF7y8uB//7MAPG7c3qPAl1wC118PrVsHW6OISJRR8A0RBV8R2a+1a+Gll2DEiL23Sz7tNLjhBruNiQmuPhGRKKHgGyIKviJyUHl51g7tmWdst7iCn6MtW8J118Fll9mIsIiIlAkF3xBR8BWRUlm8GP71L3jhBcjIsHOVKxdOg2jbNtj6REQikIJviJQk+Obm5rJ79+7yLUzCUkJCAjH66Ds6bN8Oo0fD8OEwZ07h+d694cYb4YwzIDY2uPpERCKIgm+IHCj4eu9Zu3Yt6enpwRQnYScmJoZmzZqRkJAQdClSXryHqVMtAL//vk2LAGje3OYBaxqEiMhhU/ANkQMF3zVr1pCenk6dOnVITk7GORdMkRIW8vLyWL16NfHx8TRp0kT/XqLRr7/aNIgRI2DLFjtXpYqF3xtugFatgq1PRCRMKfiGyP6Cb25uLgsWLKBOnTrUrFkzuAIlrGRkZLB69WpatmxJfHx80OVIUAqmQTz9NMydW3j+9NPhppvglFNAvxiJiJTYoQZfTT4soYI5vclqWC+lUDDFITc3N+BKJFCVK9umF7NnWzeIM8+0oDthApx6KrRrB889ZwFZRETKjIJvKenjaikN/XuRvTgHffrYdsgLFtiit6pVYd48uOYaaNQI7rjDpkiIiEjIKfiKiAShZUt46ilYuRKGDYMWLSA9HR57zBbCnXcefPppYY9gERE5bAq+climTZuGc65UnS7S0tIYNmzYAZ/jnGPcuHGHWx49evTgtddeK/HzzzvvPJ544onDvq5IiaWk2Dzf+fOtC0Tv3tYJ4p13oGdPOPpoGDkSdu4MulIRkbCn4Cvl7ttvv+Wqq64q8+t88MEHrF27lkGDBpX4Nffddx8PPfQQmZklnicvEhqxsdC/P3zyCcyaBVdeCUlJ8OOP1gWiaVO47z5YsyboSkVEwpaCr5S72rVrl8siwaeffprLLrusVBtIHHXUUaSlpfHqq6+WYWUiB9G+PTz/vE2DePRRm/u7fj08+KAF4Isvhm+/DbpKEZGwo+B7OLy3VdhBHKWY95eXl8c//vEPWrZsSWJiIk2aNOGhhx4CoFevXlx//fV7PX/Tpk0kJiYyZcoUALKzs7njjjto3LgxiYmJtGrVihdeeGG/13vnnXdo164diYmJpKWl8fjjj+/1ePGpDgsXLqRHjx4kJSXRtm1bPv744wN+P6NGjaJmzZpkZ2fvdX7gwIFccsklAGzcuJFPPvmEAQMG7Hl82rRpJCQk8Nlnn+059/jjj1OrVi3WFBlFGzBgAK+//voBaxApFzVrwp13wtKl8NZbcPzxsHs3vPoqHHOMff3mm3ZOREQOznuvo8gBpAA+IyPDF5WVleXnzp3rs7KyCk9u2+a9RdDyP7Zt8yV1xx13+OrVq/uRI0f6RYsW+c8++8yPGDHCe+/9q6++6qtXr+537ty55/lPPfWUT0tL83l5ed57788//3zfuHFj/+677/rFixf7Tz75xL/xxhvee++nTp3qAb9lyxbvvfffffedj4mJ8X/729/8/Pnz/UsvveQrVarkX3rppT3v37RpU//kk096773Pzc317du39yeddJL/8ccf/fTp033nzp094MeOHbvP72fHjh0+NTXVv/XWW3vObdiwwSckJPgpU6Z4770fO3asr1y5ss/Nzd3rtbfffrtv2rSpT09P9zNnzvSJiYn+3Xff3es5EyZM8ImJiXv9mRyqff67ETkc337r/ZAh3sfHF/48aNjQ+4cf9n7DhqCrExEpFxkZGR7wQIovTc4rzZOj4Yi04JuZmekTExP3BN3idu7c6WvUqOHffPPNPec6derkH3jgAe+99/Pnz/eA//jjj/f5+uLBd/Dgwf6UU07Z6zm33367b9u27Z6viwbfSZMm+djYWL9ixYo9j3/00UcHDL7ee3/NNdf4fv367fl62LBhvnnz5nvC+pNPPumbN2/+m9dlZ2f7zp07+/PPP9+3a9fOX3HFFb95zk8//eQBv2zZsv1ev6QUfKXMrFnj/f33e1+nTuHPhaQk7//wB+9//jno6kREytShBl9NdTgcycmwbVswRwnnyM6bN4/s7Gx69+69z8cTExO5+OKLefHFFwGYOXMmP/30E0OHDt3zdWxsLD179izx9bp3777Xue7du7Nw4cJ9buIwb948mjRpQqNGjfacO+644w56nSuvvJLJkyezatUqAF566SWGDh26p29uVlYWSUlJv3ldQkICo0eP5p133iErK2uf3SUqVaoEwI4dOw5ah0hg6tWDBx6A5cvh5ZehSxfr/PDCC3DUUdCrF4wdq2kQIiJFxAVdQFhzznZkqsAKQtyBXHHFFXTq1ImVK1fy4osv0rt3b5o2bVri1xflvf/Npg3e738+8r4eK8mmD507d6Zjx46MGjWKU089lVmzZjF+/Pg9j9eqVYstW7bs87VffvklAJs3b2bz5s1ULvZ3uHnzZsAW4YlUeImJcMklMGQIfPGF9QZ+912YOtWOunXt8csvhyOPDLpaEZFAacQ3wrVq1YpKlSrxv//9b7/P6dChA127dmXEiBG89tprXH755Xs9lpeXx/Tp00t0vbZt2/L555/vde7LL7+kdevWxMbG7vP5y5cvZ/Xq1XvOzZgxo0TXuuKKK3jppZd48cUX6dOnD40bN97zWOfOnVm7du1vwu/ixYu55ZZbGDFiBN26deOSSy4hLy9vr+fMnj2bRo0aUatWrRLVIVIhOAcnnABvv22L4e68E+rUgXXrbFOMNm2ge3d48UX71EhEJBqVZl5ENByUZo5vmHjggQd89erV/csvv+wXLVrkZ8yY4f/73//u9Zznn3/eJyQk+GrVqv3mexw6dKhv3LixHzt2rF+yZImfOnXqnjnBxef4fv/993stbhs5cuRBF7e1bdvW9+7d28+cOdN/+umn/uijjz7oHF/vbX5PcnKyT0hI2LPYrkBOTo6vU6eOHz9+/F7njjvuOH/uued6771fs2aNr1Wrlv/nP/+512svvfRSf/nllx/sj7VEwvnfjUSAXbu8HzvW+zPP9D4mpnAucJUqNhd46lTviy0AFREJB1rcpuC7X7m5uf7vf/+7b9q0qY+Pj/dNmjTxDz/88F7P2bp1q09OTvbXXnvtb16flZXlb7nlFl+/fn2fkJDgW7Zs6V988UXv/W+Dr/fejxkzxrdt23bPtR577LG93q9o8PXeFtCdcMIJPiEhwbdu3dpPnDixRMHXe++HDBnia9Sosc8ODHfddTXgqDoAABljSURBVJcfNGjQnq//+te/+vr16/uNGzfuOTdu3DifkJDgf/zxxz3fa0pKip8xY8ZBr10S4fzvRiLMqlXeP/KI9y1b+r0WyjZs6P2tt3r/3Xfe5y8OFRGp6A41+DrvtQ98Uc65FCAjIyODlJSUPed37tzJ0qVLadas2T4XTYW7FStWkJaWxrfffkuXLl2CLqfETjnlFNq0acPTTz/9m8fWrVtHu3bt+P777/fMWT6YZ599lvfee4/JkyeHpL5I/3cjYch7+OwzGDUKxoyBjIzCx1q3hsGDYdAgOOKI4GoUETmIzMxMUlNTAVK99yXeblVzfKPc7t27Wb58OXfeeSfdunULm9C7efNm3njjDaZMmcJ11123z+fUrVuXF154geXLl5f4fePj4xk+fHioyhSpeJyDHj3gv/+1+b9jx8L559v2yAsWWKeII4+Etm3hnnvgm2+g2Dx4EZFwpRHfYqJtxHfatGmcfPLJtG7dmjFjxtChQ4egSyqRtLQ0tmzZwl/+8hf+9Kc/BV3OfkXqvxuJQFu3wrhx8Npr8MknkJNT+FiDBnDWWXb06AGl7PYiIhJqhzriq+BbTLQFXylb+ncjYSk9HT76yILwhAl7d4FISrLw27evHe3b2yiyiEg5OtTgqz6+IiKyt2rV4MIL7cjOhin/396dh0lVXvse/y666YYGBASZaVqZZFAJilFEAYEYozEON2BiQAPEGCPK0YsgXBxijkeDAxqN16BwkGg86g2i5pqDGEMOgwniLBDUMKiNTGIPMjRNr/PHW9Vd3XTT1Rxgd1G/z/O8z67ae1fVKnY1terda7/vn0MS/PLLkJ8PCxeGBtC+PYwYAYMHh+HSevRQIiwi9ZYSXxERqVl2Npx/fmjusHp1ReL7l7/Apk3hQrknnwz7H3dcSIAHDYIzz4R+/ZKeaVJE5HBTqUMVKnWQQ0mfGzmq7d4Ny5aFmuClS+Fvfws9xIkaNAgXy/XvX9H69YNwilJE5KCo1EFERI6sRo3g3HNDg5D0vvUWLFkShkz7+9/DyBGrVoX2u99VPLZDh5AQV20dO4ZkWUTkMEiZxNfM1gNVB2O9x92nJOyTCzwCnAvsAp4G/re7lxypOEVE0lZ2dihvOPNMmDQprNu0KSTDiW3jxlArnJ8f6ocTNWwInTtDly6QmxuW8dtt24bWqhVUMwW6iEhtUibxjbkVmJVwv/xSYzPLAP4IbAUGAa2AuYABE45gjCIiEte+PVxwQWhxO3bAP/4Ba9ZUbh9/DHv3wj//GVpNGjSA1q1DEtymTVi2bBkuyqupNW8eWlbW4X/PIkej4mLYujX8jbZtC8cck5IXsqZa4lvk7l/UsO1bQG+gs7vnA5jZTcC/m9m0utR/SPLi4wDv2LGDFi1aJPWYvLw8Jk6cyMSJE2vcx8yYP38+F198cdKxDBkyhH79+jFz5sykH5OMg4lFRA6gZUs444zQEpWWhl7gDRtC27ix4vZnn4Wyie3bw4QaW7aEVleNGoUv7ObN67ZMvN2smXqc5ehXUhKGNfyP/wi1/Bs2VN7epQtcfDHccAMcf3w0MR6EVEt8J5vZdOBT4DlgRkIZw5nAB/GkN+Y/gWzgVOD16p7QzLJj+8Q1O+RRSyUrVqygSZMmUYdR7vbbb2fNmjU888wzUYcikt4yM0NJQ24unH129fuUlsK2bSHp3by5YvnVVwduRUXh8bt3h3YwSXOipk33T45btAg90O3ahR6xdu0qWps26m2W+s891ObPmwfPPBN+aCbKyQk/+oqKQiL84IPw6KNhxsfJk1OiPj+VEt8HgbeAHcDpwL8BxwPjY9vbAZsTH+DuO8ysJLatJrcAtx3yaKVGxx13XNQhVPLiiy8yKV6PKCL1W2ZmRTJZF6Wl4cu6oAAKC2tf1rStJNbXUlwc2uefJx/DsceGi/fidct5eZVvH3dcSp46lqPA+vXh4tN588LU5XHt2sEVV8B3vhNGZImf2S0qCsMZzpwZ6vSnToV334W5c0Otf33m7pE14HbAa2mn1fDYy2LbW8Xu/xb4z2r2KwEuP0AM2cAxCa0j4AUFBZ5o165dvmrVKt+1a1fFyrIy95LiaFpZmSdr3759fvfdd3vXrl09KyvLO3fu7L/85S/d3X3o0KH+85//vNL+27Zt86ysLH/ttdfc3X337t0+adIk79Spk2dlZXm3bt388ccfd3f3119/3QHfsWNH+eOff/557927t2dlZXmXLl383nvvrfT8Xbp08QceeKD8/tq1a/3ss8/27Oxs79Wrly9cuNABnz9/fo3vqbi42EePHu1NmjTxdu3a+b333uuDBw/2G264wd3d77jjDu/bt+9+j+vfv79Pnz69/P7GjRu9YcOG5fHXFsvcuXO9SZMmvnbt2vLnuO6667x79+5eXFy83+tV+7kRkdS1e7f75s3uH33k/uab7n/+s/v8+e5z57rPnOk+dar72LHuF1zgfuqp7h07umdmuoe+tAO3xo3d+/Rxv+QS9ylT3GfPdl+yxH3r1jr9ny+SlB073GfNcj/nnP0/h1dc4f6nP7nv3Xvg5ygrc3/8cfeGDcNjY9/BR0JBQUE8TzzG65B7Rt3j+zBQ2/nl9TWsfyO27AZsB74Avpm4g5m1BBpSpSc4kbvvAfYkPKaWcBKU7oSHmia//6F0fTE0TK5c4JZbbmHWrFk88MADDBo0iE2bNrFmzRoAxo8fz3XXXcd9991HduxX2lNPPUWHDh0YOnQoAGPGjGH58uU89NBDnHLKKaxbt45t27ZV+1orV65k5MiR3H777YwaNYply5Zx7bXX0qpVK6666qr99i8rK+PSSy+ldevWvPHGGxQWFh6w9jdu0qRJvP7668yfP5927doxdepUVq5cSb9+/QAYO3Ysd9xxBytWrGDAgAEAvPfee7z99ts899xz5c/z4osvcs4559CiRYukYhkzZgwvv/wyV1xxBcuWLWPRokU89thjLF26tF6Vb4jIYZKdHcoW2rRJ/jFlZeGCvk2bQq1yvG55w4bQ07ZhQ6ht3rULPvwwtKpatoSePUPr0wdOOilMF92xo3qJJXlFRfDSS6Fu909/qjiDYRaGJRw9Gi69NNSxJ8MMxo0LZy0mTw49v/VcpImvu28Dqs+gaveN2HJTbLkcmGZm7d09vu5bhKR25cFHmdqKiop48MEHefjhh7nyyisB6Nq1K4MGDQLgsssuY8KECSxYsICRI0cCMGfOHK666irMjLVr1/Lss8/y6quvMnz4cABOOOGEGl/v/vvvZ9iwYUyfPh2AHj16sGrVKmbMmFFt4rto0SJWr17N+vXr6dSpEwB33XUX559/fo2vUVxczBNPPMGTTz7JiBEjAJg7d2754wE6derEeeedx5w5c8oT3zlz5jB48OBK8S9YsIDvfe97dYrlscce4+STT+b666/nD3/4A7fddlv5a4iI7KdBgzAEW6tWIVmtTkkJfPppGNli7dow6kV8uXFjSJzfeCO0RC1ahOfs27ciGe7bN5RVyOFVWlpR8lJSUtFnWlYWlpmZ0LhxqItt3DjUeB/pHymlpaEE4dVXQ1uypCLZBejdG8aMCeUMCd+hdTZ8OKxYoRrfQ8XMzgTOIFygVgAMAB4AXnT3jbHdFgKrgHlmNgk4FrgXmOWHa0SHzJzQ8xqFzOSmAF29ejV79uxh2LBh1W7Pzs7mRz/6EbNnz2bkyJG88847vPvuu7zwwgsAvPPOO2RkZDB48OCkXy+eSMadddZZzJw5k3379pFR5Uro1atXk5ubWylpPfPMMw/4Gp988gklJSWV9jv22GPp2bNnpf1+8pOfMHbsWO6//34yMjJ46qmnuO+++8q3FxYWsnjxYmbNmlWnWFq2bMkTTzzBeeedx8CBA5kyZcp++4iI1ElWFnTtGtp551XetmtXSIj/8Y8wZfSHH8L774f7X30VkpklSyo/pkOH/RPi3r01fXRtSkpCL/wnn4QfHJs3wxdfhGW8FRSEntPdu+v23Gbh379Fi/DDJN5atqx8Pz78XtVRRZo12z+x3LcvxFJYGC7YjJ9F+PhjePvtkPRWjbN7dxg1KrSafogdjBRIeiFFEl9Cr+0owkVo2cAGwni+v4rv4O77zOwC4DfAUhImsDhsUZklXW4QlcaNG9e6z/jx4+nXrx+fffYZs2fPZtiwYXTp0iXpxydy9/3KRfwA02JXt622cpMDPV+i7373u2RnZzN//nyys7PZs2cPl112Wfn2V155hV69epW/17rE8te//pWMjAzy8/P5+uuvK01vLSJySDVuHJLXk06qvH7PnpD8fvBBaO+/H5br11dMELJwYcX+ZiGxjifE8aS4e/cwcUi6KC4O40R//HFIcD/5pOL2xo2hx7YuMjNDCYxZSP7MQisthZ07K57PHb7+OrS6XBSZKH6c4q9RdYrw6hxzDAwZAiNGhNajR1qXx6RE4uvubxF6fGvbbyNw4eGPKHV0796dxo0b89prrzF+/Phq9znppJM47bTTmDVrFk8//TS//vWvK20rKytj8eLF5aUOB9K7d2+WVOl5WLZsGT169Nivtze+/8aNG8nPz6dDhw4ALF++/ICv0a1bNxo2bMgbb7xBbm4uADt27GDt2rWVeqYzMzO58sormTNnDtnZ2Vx++eXkJPR2LFiwgIsuuqjOsSxbtoxf/epXvPTSS0yZMoUJEyYwd+7c2v5pREQOrexsOPnk0BIVFYVe4aoJ8ZYtIcH7+GOIndUDQjJ14omVe4f79AnDymWmRJpQmXsYhqtqUhu/vbnGy36CnJzwAyEvr/LQdPGZA1u2DMPZNWsWlgcaxcA9TPiwa1dIgnfuDL30X365f9uxo2IZH1kkPprI3r3h+eLLqrKzQ29xXl5FO/nkMBJDt24p0xt7JKTgJ1rqolGjRkyePJmbb76ZrKwszjrrLLZu3cqHH37IuHHjyveLX+SWk5PDJZdcUr4+Ly+PK6+8krFjx5Zf3LZhwwa2bNlSXhOc6KabbmLAgAHceeedjBo1iuXLl/Pwww/zm9/8ptr4hg8fTs+ePRkzZgz33XcfhYWFTJs27YDvqWnTpowbN45JkybRqlUr2rZty7Rp02hQzR/2+PHj6dWrFwBLly4tX19aWsorr7zCokWL6hRLUVERo0ePZsKECZx//vnk5uZy2mmnceGFF/L973//gHGLiBwRzZpVP0HIli0VyXBiQlxcHG6//z78/vcV+zdsGBKobt0qt65dQz1olBf07twZek3js/xVbYW1VDgee2x4H/H3E1927RqS3EPVI2oWyliyskLZwsFwDz278QQ4cQyGnJzQo1vfhxCrR5T4poHp06eTmZnJrbfeSn5+Pu3bt+eaa66ptM8PfvADJk6cyA9/+EMaNWpUadujjz7K1KlTufbaa9m+fTu5ublMreHKzf79+/Pss89y6623cuedd9K+fXt+8YtfVHthG0CDBg2YP38+48aN4/TTTycvL4+HHnqIb3/72wd8TzNmzKC4uJiLLrqIZs2acdNNN1FQULDfft27d2fgwIFs376db36zYtCPxYsX07RpU0499dQ6xXLDDTfQpEkT7rrrLgD69OnDPffcwzXXXMPAgQPp2LHjAeMWEYlMmzbhyv1zz61Y5x5O78eT4HhbsyYkWx99FFp1mjcPtcTx1r59mEq6pqmiGzasaJmZoRRg795QVxtfFhaGHtGCgorJR7ZvDyNixMs38vPD+tp06FB9Ytu1a+i1TRVmYcbBKt/NcnAs2XrJdGFmxwAFBQUFleo2d+/ezbp16zj++OP3SwyPBp9++il5eXmsWLGC/v37Rx3OIePunHjiifz0pz/lxhtvLF9//fXXU1paWmNP9KFytH9uROQoVVYWelTjpRHxMoH47eKILuxOlJMTpso94YSK1rVrWOblhdpoOWoVFhbSPPSiN6/LIAbq8U1ze/fuZdOmTUyZMoUzzjjjqEp6t2zZwrx58/j888/58Y9/XGlb3759ax09QkQkbTVoAJ07hxYb072ce6gjTuyBjbcvv6zoqY332sZP0cdbVYm9wfGpn+O9xS1ahN7ZxJ7leDvmmLS+SEsOjhLfNLd06VKGDh1Kjx49eP7556MO55Bq27YtrVu35re//S0tq5zWuvrqqyOKSkQkxZlVDLF14ol1e6x7GIKrtDSUO2RkKHmVI0qJb5obMmRI0sODpZqj9X2JiKQss5DwpuJoEXJU0PgWIiIiIpIWlPjWkXoRpS70eREREak/lPgmqWFstpSdO3dGHImkkpLYnOjVTd4hIiIiR5aKbJKUkZFBixYt2LJlCwA5OTm1Tq0r6a2srIytW7eSk5NDpurZREREIqdv4zpo164dQHnyK1KbBg0akJubqx9JIiIi9YAS3zowM9q3b0+bNm3YW9N82SIJsrKyqp1KWURERI48Jb4HISMjQzWbIiIiIilGXVEiIiIikhaU+IqIiIhIWlDiKyIiIiJpQTW+NSgsLIw6BBERERGpxsHmaaaZpSozs47AZ1HHISIiIiK16uTunye7sxLfKiwMuNoBKDpCL9mMkGh3OoKvKYeOjl/q0zFMfTqGqU/HMLVFdfyaAfleh2RWpQ5VxP7xkv7l8D+VMLFBkburviLF6PilPh3D1KdjmPp0DFNbhMevzq+li9tEREREJC0o8RURERGRtKDEN3p7gDtiS0k9On6pT8cw9ekYpj4dw9SWMsdPF7eJiIiISFpQj6+IiIiIpAUlviIiIiKSFpT4ioiIiEhaUOIrIiIiImlBiW+EzOxaM1tnZrvNbKWZnR11TJIcM7vFzFaYWZGZbTGzF8ysZ9RxycGJHU83s5lRxyLJM7OOZvY7M9tuZjvN7B0zOzXquCQ5ZpZpZr+MfQ/uMrN/mtmtZqbcpJ4ys3PM7CUzy4/9n3lxle1mZrfHtu8ys7+YWZ+o4q2OPlwRMbNRwEzgX4FvAP8FvGJmuZEGJskaDDwCnAGMIMyCuNDMmkQaldSZmQ0ArgbeizoWSZ6ZtQSWAnuB84HewE3AV1HGJXUyGbgGuA7oBdwMTAImRBmUHFAT4F3CMavOzcCNse0DgC+AV82s2ZEJr3YaziwiZvY34C13/1nCutXAC+5+S3SRycEws+OALcBgd/9r1PFIcsysKfAWcC3wf4B33H1itFFJMszsbuAsd9eZshRlZi8Dm919XMK6/wfsdPfR0UUmyTAzBy5x9xdi9w3IB2a6+z2xddnAZmCyuz8WWbAJ1OMbATPLAk4FFlbZtBAYeOQjkkOgeWz5ZaRRSF09AvzR3RdFHYjU2UXAm2b2XKzc6G0z+0nUQUmdLAGGmVkPADM7BRgE/P9Io5KDdTzQjoTcxt33AIupR7lNZtQBpKnWQAbhV1CizYQPjaSQ2K/c+4El7v5B1PFIcszscqA/4XScpJ4TgJ8R/vbuAk4HHjKzPe7+ZKSRSbLuIXQarDGzfYTvxWnu/vtow5KDFM9fqsttuhzhWGqkxDdaVetMrJp1Uv89DJxM6KmQFGBmnYEHgW+5++6o45GD0gB4092nxu6/HbuI5meAEt/UMAr4EfBD4EOgHzDTzPLdfW6kkcn/RL3ObZT4RmMbsI/9e3fbsP8vJanHzOzXhFOu57j7Z1HHI0k7lfD3tjJ02AOht+kcM7sOyHb3fVEFJ0nZBKyqsm41cFkEscjBmQHc7e7PxO6/b2ZdgFsAJb6p54vYsh3h7zOuXuU2qvGNgLuXACsJowEkGgEsO/IRSV3Fhmx5GLgUONfd10Udk9TJa8BJhB6meHsTeArop6Q3JSwFqg4h2APYEEEscnBygLIq6/ah3CRVrSMkv+W5TeyapsHUo9xGPb7RuR+YZ2ZvAssJwynlAv830qgkWY8QTs99Dygys3jvfYG774ouLEmGuxcBleqxzexrYLvqtFPGA8AyM5sKPEuo8b061iQ1vARMM7ONhFKHbxCGwpodaVRSo9hION0SVh1vZv2AL919Y2ws9Klm9hHwETAV2Ak8feSjrZ6GM4uQmV1LGPOuPeFL+F80FFZqiA3jUp0fu/u/H8lY5NAws7+g4cxSipldCPwb0J3Q23S/u8+KNipJVmxs1zuBSwinw/OB3wO/iJ0ZlXrGzIYAr1ezaa67XxW72Ps24KdAS+BvwM/rU4eCEl8RERERSQuqoxERERGRtKDEV0RERETSghJfEREREUkLSnxFREREJC0o8RURERGRtKDEV0RERETSghJfEREREUkLSnxFREREJC0o8RURERGRtKDEV0RERETSghJfEREREUkLSnxFRI5CZnacmX1hZlMT1n3TzErM7FtRxiYiEhVz96hjEBGRw8DMvgO8AAwE1gBvA39094mRBiYiEhElviIiRzEzewQYDqwATgEGuPvuaKMSEYmGEl8RkaOYmTUGPgA6A6e5+3sRhyQiEhnV+IqIHN1OADoQ/r/vEnEsIiKRUo+viMhRysyygL8D7xBqfG8ETnL3zZEGJiISESW+IiJHKTObAfwvQm1vMfA6UOTuF0YamIhIRFTqICJyFDKzIcBEYLS7F7p7GTAaGGRmP4s0OBGRiKjHV0RERETSgnp8RURERCQtKPEVERERkbSgxFdERERE0oISXxERERFJC0p8RURERCQtKPEVERERkbSgxFdERERE0oISXxERERFJC0p8RURERCQtKPEVERERkbSgxFdERERE0sJ/A68UhjCkiNrqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x119fb05c0>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Delta t: 1.9743860485041633 s\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "xaxis = np.linspace(XBEGIN, XFINAL, 1000)\n",
    "thetas = cheby_cycloid_coefficients(TERMS)\n",
    "thetas = np.array(thetas).reshape(-1,1)\n",
    "func = get_cheby_path_func(thetas)\n",
    "deriv_func = get_cheby_path_func_deriv(thetas)\n",
    "plt.figure(figsize=(8,4), dpi=100)\n",
    "y_points, = plt.plot(xaxis, func(xaxis), color=\"red\",label='cycloid y(x)')\n",
    "yp_points, = plt.plot(xaxis, deriv_func(xaxis), color=\"darkorange\",label='cycloid dy/dx')\n",
    "plt.legend(handles=[y_points, yp_points], loc=3)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.show()\n",
    "print('Delta t:',delta_time(func,deriv_func),'s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Alternative functions\n",
    "\n",
    "To see how other functions compare to the optimal cycloid, we can decompose them into their fourier coeffcients and plug them into $\\Delta t$ just the same. One such example is the straight line that connects $a = (0,0)$ to $b  = (x_f,y_f)$. The coeffcients for the line are:\n",
    "\n",
    "\n",
    "$$ a_n = \\frac{1}{2\\pi k} \\int_{-2\\pi k}^{2\\pi k} \\frac{y_f}{x_f}x \\cos(\\frac{nx}{2 k}) dx = 0 \\qquad b_n = \\frac{1}{2\\pi k} \\int_{-2\\pi k}^{2\\pi k} \\frac{y_f}{x_f}x \\sin(\\frac{nx}{2 k}) dx $$\n",
    "\n",
    "All the $a_n$ will be zero because the line is an odd function. If we use few terms in the fourier expansion, the result will not look like a straight line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 224,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What are the chebyshev coefficients for a line? \n",
    "def cheby_line_coefficients(n):\n",
    "    def f(x,i):\n",
    "        u = x/C\n",
    "        y = -0.5*u\n",
    "        coef = 1./(C*np.pi) if i == 0 else 2./(C*np.pi)\n",
    "        cheby = special_funcs.eval_chebyt(i,u)\n",
    "        results = (coef*cheby*y)/np.sqrt(1-(u**2))\n",
    "        return results\n",
    "    thetas = [integrate.quad(lambda x: f(x,j),-C,C,limit=200)[0] for j in range(n)]\n",
    "    return thetas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 225,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAArUAAAFwCAYAAACmbYPbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XmYjmUbx/HvNWOsY0ZlT2bGviWEqNdaKEWWIlmyVwghS6ko+15CZY0skSyppIVIRUiLfZtRESEzMrZm7vePS0ITY8zM/dzz/D7HcR+Ze57lnHne4+3X5byu0ziOg4iIiIiIlwW4XYCIiIiIyPVSqBURERERz1OoFRERERHPU6gVEREREc9TqBURERERz1OoFRERERHPU6gVEREREc9TqBURERERz1OoFRERERHPU6gVEREREc9TqBURERERz0vndgGpyRhjgLzACbdrEREREZH/lBU44DiOk9gn+FWoxQbaX9wuQkRERESuKh/wa2If7G+h9gTAzz//TEhIiNu1iIiIiMhlYmJiuOWWW+Aa/2bd30ItACEhIQq1IiIiImmINoqJiIiIiOcp1IqIiIiI5ynUioiIiIjnKdSKiIiIiOcp1IqIiIiI5ynUioiIiIjnKdSKiIiIiOd5LtQaYzoZY/YZY04bYzYaY6q4XZOIiIiIuMtTodYY0xQYBwwGygJrgI+MMfldLUxEREREXOWpUAv0AKY6jjPFcZxtjuN0B34GnkzowcaYDMaYkL8vIGtqFtu8OYweDSdPpua7ioiIiPgfz4RaY0x64HZgxWXfWgHc+R9P6wdEX3T9kmIFXuabb2DOHOjVC8LDYehQiIlJrXcXERER8S+eCbVAdiAQOHTZ/UNA7v94zlAg9KIrX4pVd5ly5WDKFChYEI4cgWefhbAwGDAA/vgjtaoQERER8Q9eCrV/cy772iRwzz7Qcc44jhPz9wWcSPHqzkufHtq1g+3bYdYsKFYMjh+HgQNtuH32Wfj999SqRkRERCRt81KoPQLE8e9V2Zz8e/XWZ6RLBy1awE8/wfz5ULo0nDhh2xHCw6FnTzh40O0qRURERLzNM6HWcZyzwEag1mXfqgV8lfoVXZvAQHj4YfjuO1i8GMqXh9hYGDMGIiKgSxf4+We3qxQRERHxJs+E2vPGAO2NMW2NMcWNMWOB/MDrLteVaAEB8OCDsH49fPQR3HknnDkDEybY/tuOHWHvXrerFBEREfEWT4Vax3HeAboDLwCbgapAXcdxolwtLAmMgXvvhS+/hM8+gxo14Nw5mDwZihSBxx6DHTvcrlJERETEG4zjJLjHKk06f1ZtdHR0NCEhIW6X8y9r18LLL8PHH9uvjYGmTeG556BUKXdrExEREUkNMTExhIaGAoSe3+ifKJ5aqU3r7roLli+3rQn164PjwLx5cOut0KgRbNrkdoUiIiIivkmh1gdVqABLlsDmzXZzmTGwaBHcfjs88IAd7CAiIiIi/1Co9WG33WaPAfvpJztyNyAAPvgAKleGWrVg9Wq3KxQRERHxDQq1HlCiBLz9tt041ratPfv200+hWjWoWhU++cS2KoiIiIj4K4VaDylUCKZOhV274Ikn7NSyNWugdm27ertsmcKtiIiI+CeFWg8KD4dJk+x5tt26QcaMsG4d1Ktn+27few/i492uUkRERCT1KNR62M03w7hxEBkJvXtDlix2YlnjxnYc79y5EBfndpUiIiIiKU+hNg3IlQuGD4eoKOjfH0JCYMsWePRRKF4cZsywgx1ERERE0iqF2jTkppvs8IaoKPvPG2+0/bdt2kDRovDmm3Ykr4iIiEhao1CbBmXLZldsIyNhxAjImRP27YPHH7ebzcaPh1On3K5SREREJPko1KZhWbPCM8/YQDtuHOTNC7/8Al27QoECMHo0/Pmn21WKiIiIXD+FWj+QObM9JWHvXntqQlgY/PYb9OplT1IYMgRiEj1ZWURERMT3KNT6kQwZ7Pm2u3bZ824LFoSjR+G552zQffFFOHbM7SpFRERErp1CrR8KCrKTybZvt5PKiheH48fhpZfsym2/fnD4sNtVioiIiCSeQq0fS5cOmjeHn36CBQvgttvgxAkYNsyG2x494OBBt6sUERERuTqFWiEgAB56yA5uWLIEype3pyOMHQsREdClC+zf73aVIiIiIv9NoVYuMAbq14f162H5crjrLnuu7YQJ9iiwDh3sZjMRERERX6NQK/9iDNSpA2vWwOefQ82adiLZlClQpAg89pjtxxURERHxFcZxHLdrSDXGmBAgOjo6mpCQkBR9L8dxOBefom+Rqtats722n3xivzYGGje25+CWLOlubSIiIpLyggLAGJPi7xMTE0NoaChAqOM4iT50VKE2hZyNcxjzw9EUfQ8RERGR1NKj9E2kD/TdUKv2AxERERHxvHRuF5BWBQXY/6JJ63bsgJEjYf58iD/fblG9OvTtC//7n6uliYiISDIK8vGlULUfSLLYvdv23L71Fvz1l71XpQr07w+1atkeXBEREZGrUfuBuKpQIXs6wu7d8OSTkD69PT2hTh2oVAmWLQM/+u8nERERSWUKtZKswsJg4kTYtw+6d4dMmey5t/XqQblysHDhP20KIiIiIslFoVZSRN68diJZZCT06QPBwbB5s51cduutMGcOxMW5XaWIiIikFQq1kqJy5rS9tpGR8PzzEBoKW7dC8+ZQvDjMmGEHO4iIiIhcD4VaSRU33QQvvQRRUTBoENx4I+zaBW3a2Cllb7xhR/KKiIiIJIVCraSq0FB47jkbbkeMsCu5kZHwxBNQsCC8+iqcOuV2lSIiIuI1CrXiiuBgO2I3MhJeeQVuvhl+/RW6dYOICBg1Cv780+0qRURExCsUasVVmTJB166wZw+8/ro9PeHQIRt4w8Nh8GCIjna7ShEREfF1CrXiEzJkgMcft32206bZc2+PHrXDG8LD4cUX4dgxt6sUERERX6VQKz4lKMhuHtu2DWbPhhIl4Phxu8ksLMyO3z182O0qRURExNco1IpPSpcOHn0UfvwRFiyA226zPbbDh9uV2x494MABt6sUERERX6FQKz4tIMAObPjuO1i6FCpUsKcjjB0LBQpA586wf7/bVYqIiIjbFGrFE4yxo3bXrYOPP4b//c+eaztxoj0KrH17u9lMRERE/JNCrXiKMVC7NqxeDStXQs2a8NdfMHUqFC0KrVrB9u1uVykiIiKpTaFWPMkYqF4dPvsM1q6F++6DuDiYNctuLmva1PbjioiIiH9QqBXPu/NO+PBD2LABGjQAx4H586F0afv1xo1uVygiIiIpTaFW0ozbb4dFi+D776FJE7uau2QJlC8PdevC11+7XaGIiIikFIVaSXNKl4Z33oGtW6FlSwgMhI8+siu6d98Nq1bZ1VwRERFJOxRqJc0qVgxmzoQdO+zpCOnSweefQ40aULWqPUVB4VZERCRtUKiVNK9gQZg82R751amTHcn75Zdw771QqRK8/77CrYiIiNcp1IrfyJ8fJkyAvXvh6achUyZYvx7q14eyZeHddyE+3u0qRUREJCk8E2qNMc8ZY74yxsQaY467XY94V968MGYMREZCnz4QHGw3lz38MJQqBbNn27NvRURExDs8E2qB9MACYJLbhUjakDMnDBsGUVHwwgsQGgrbtkGLFlC8OEyfDufOuV2liIiIJIZxPNZMaIxpDYxzHCdbIh6bAchw0a2swC/R0dGEhISkUIXiVdHRtj1hzBg4etTeCwuDvn2hTRvbiysiIiIpKyYmhtDQUIBQx3FiEvs8L63UJkU/IPqi6xd3yxFfFhoKzz5r2xJGjoRcuewq7pNP2s1mr74KsbFuVykiIiIJSeuhdigQetGVz91yxAuCg6FXL9i3zwbZm2+GX3+Fbt0gIsIG3j//dLtKERERuZirodYYM8AY41zlKp/U13cc54zjODF/X8CJZCxf0rhMmeCpp+xRYG+8AeHhcPgw9O5t/zx4sG1ZEBEREfe52lNrjMkOZL/KwyIdxzl90XNak8ie2gTeLwSIVk+tJMW5c/ZkhCFDYNcuey80FLp2he7d4cYb3a1PREQkLUhqT22a3iiWwHMVauW6xcXB/PkwaJAdxQu2ZaFTJ+jZ056qICIiIkmT5jeKGWPyG2PKAPmBQGNMmfNXsNu1iX8JDIRmzeDHH+3AhjJlbI/tiBG2LeHpp20ProiIiKQez4Ra4CXgO2AgEHz+z98BSe65FbkeAQHQuDFs2mRH7VasCKdOwbhxUKCAXbmNinK7ShEREf/gufaD66H2A0lJjgOffgovvwxr1th76dJBq1bQrx8UKuRufSIiIl6Q5tsPRHydMVCrFqxeDatWwd1323G706ZB0aLQsqWdWCYiIiLJT6FWJAVUq2ZXbb/6CurWhfh4ePttKFkSmjSBH35wu0IREZG0RaFWJAVVrgwffAAbNkDDhrZFYcECuO02aNDA3hcREZHrp1Arkgpuvx3ee8+u0DZtalsVliyBChXgvvvsiq6IiIgknUKtSCq69VaYN8/21rZqZY8HW74c7roLataElSvtaq6IiIhcG4VaERcULQpvvQU7dkD79hAUZANtzZpQpYoNugq3IiIiiadQK+KiggVh8mTYvRs6d4YMGWDtWtuScMcdsHSpwq2IiEhiKNSK+ID8+eG112DfPujRAzJnhm+/hQcftBPLFiywJyiIiIhIwhRqRXxInjwwejRERkLfvhAcbDeXNWkCpUrB7Nn27FsRERG5lEKtiA/KkQOGDrVjdl98EbJls5vLWrSA4sXtQIdz59yuUkRExHco1Ir4sBtvhAED7MrtkCGQPbvtv23XDgoXhtdfhzNn3K5SRETEfQq1Ih4QGgr9+tlwO2oU5MplV3GffBIKFIBXXoHYWLerFBERcY9CrYiHZMkCPXvaDWXjx0O+fHDgAHTvDhERMGIEnDjhdpUiIiKpT6FWxIMyZYIuXWwrwhtvQHg4HD4MffrYPw8aBMePu12liIhI6lGoFfGwDBmgY0fYuRNmzIAiReDYMXj+eQgLs/88etTtKkVERFKeQq1IGhAUBI89Blu3wty5ULIkxMTYFduwMOjdGw4dcrtKERGRlKNQK5KGBAbCI4/Ys20XLoSyZeHkSRg50vbcdu8Ov/7qdpUiIiLJT6FWJA0KCIBGjWDjRli2zI7cPXXKnpJQoIA9NSEqyu0qRUREko9CrUgaZgzcfz98/TWsWAFVqsDZs/Z820KF7Hm3u3e7XaWIiMj1U6gV8QPGQK1asHo1fPEF3HOPHbc7bRoULWonlW3b5naVIiIiSadQK+JnqlaFTz6xq7f33w/x8TB7tt1c9vDD8P33blcoIiJy7RRqRfxUpUq233bjRmjYEBwH3n0XypSBBx+Eb791u0IREZHEU6gV8XPlysF778GPP9qTE4yBpUuhYkW4915Yu9btCkVERK5OoVZEAChVyp5xu22bPfM2MBA+/hj+9z+oUQM+/9yu5oqIiPgihVoRuUTRonY62c6d0KGDHeywahXcfbcNuMuXK9yKiIjvUagVkQQVKABvvgl79kCXLnYk71dfwX332daEJUsUbkVExHco1IrIFd1yC4wfD/v2Qc+ekDkzbNgADRrYTWULFkBcnNtVioiIv1OoFZFEyZMHRo2CyEjo1w+yZrXjeJs0sf24b79tz74VERFxg0KtiFyTHDlgyBAbbgcMgGzZYPt2aNkSihWDqVPt1DIREZHUpFArIkly443w4osQFQVDh0L27Lb/tn17KFwYJk2C06fdrlJERPyFQq2IXJeQEOjb167cjh4NuXPD/v3QqRMULAjjxkFsrNtViohIWqdQKyLJIksW6NED9u6F116zG8wOHICnn4bwcBg+HE6ccLtKERFJqxRqRSRZZcoEnTvD7t32SLCICPj9d7uaGx4OL78Mx4+7XaWIiKQ1CrUikiLSp7fDG3bsgLfegiJF4NgxeOEFCAuD/v3hyBG3qxQRkbRCoVZEUlRQELRqBVu32jG8JUtCTAwMHmxXbnv3hkOH3K5SRES8TqFWRFJFYCA88og92/a996BsWTh5EkaOtOG2Wzf49Ve3qxQREa9SqBWRVBUQAA0bwsaN8MEHUKmSPfrr1VftaN4nnrAnKYiIiFwLhVoRcYUxULcufPUVfPIJVK1qhza88YY957ZtW9i1y+0qRUTEKxRqRcRVxsA998AXX9irVi07bnf6dDuhrHlz248rIiJyJQq1IuIzqlaFFSvgm2/ggQcgPh7mzIFSpeChh2DzZrcrFBERX6VQKyI+54474P33YdMmaNQIHAcWLrSby+rXh/Xr3a5QRER8jUKtiPissmVtmP3pJ2jWzG4ye/99G3rr1IEvv3S7QhER8RUKtSLi80qWtG0I27ZB69b2eLAVK6BKFahRAz7/3K7mioiI/1KoFRHPKFLEbiDbtQs6drSDHVatgrvvhrvugo8+UrgVEfFXngi1xphwY8xUY8w+Y8wpY8weY8xAY0x6t2sTkdQXEWGP/tq7F556CjJmhK+/tkeEVagAixfbTWYiIuI/PBFqgWLYWh8HSgJPA08AQ9wsSkTclS+fHdqwbx/06gWZM9uhDg0bQpkyMH8+xMW5XaWIiKQG43j07+qMMc8ATzqOU+AKj8kAZLjoVlbgl+joaEJCQlK6RBFJZUeOwNixMH48nDhh7xUrBs8+azeapUvnbn0iInJ1MTExhIaGAoQ6jhOT2Od5ZaU2IaHAsas8ph8QfdH1S0oXJSLuyZ4dBg+GqCgYOBBuuAG2b4dWraBoUZgyxU4tExGRtMeTK7XGmILAJqCn4zhTrvA4rdSK+LGYGJg4EUaPtqu4APnzQ58+dgxvxozu1iciIv/myZVaY8wAY4xzlav8Zc/JCywHFlwp0AI4jnPGcZyYvy/gRAr+OCLiY0JCoG9fiIyEMWMgTx7Yvx86d4YCBWyrQmys21WKiEhycHWl1hiTHch+lYdFOo5z+vzj8wIrgXVAa8dxrml/szEmBIjWSq2Ifzp9GqZOheHD4eef7b0cOaBnT+jUCbJmdbc+ERFJ+kqtZ9oPjDE3YwPtRqCF4zjXvKdZoVZEwPbVzpwJQ4faY8HA9t927w5du0K2bO7WJyLizzzZfpBY51doVwE/A72AHMaY3MaY3K4WJiKelD49tG8PO3bYcFu0KPzxB7z4IoSFwXPP/dODKyIi3uCJUAvUBgoBNbEnGBy86BIRSZJ06aBlS9iyBebNg1Kl7OayIUMgPByeeQZ++83tKkVEJDE8EWodx5nhOI5J6HK7NhHxvsBAaNoUvv8eFi2CcuXg5EkYNcpOL+vaFX7RgYAiIj7NE6FWRCQ1BARAgwawYQN88AFUqmQ3l40fb09LePxxO71MRER8j0KtiMhljIG6deGrr+DTT6FaNTh3Dt58EwoXhjZtYNcut6sUEZGLKdSKiPwHY+Duu2HVKli9GmrXhrg4mDHDjt999FHbjysiIu5TqBURSYQqVeDjj+Gbb6BePYiPh7lz7eayhx6C775zu0IREf+mUCsicg3uuAOWLrUhtnFje2/hQru5rF49WL/e3fpERPyVQq2ISBKUKQPvvgs//WTbEAICYNkyG3pr14Y1a9yuUETEvyjUiohch5IlYfZs2LbNbiBLlw4++QSqVoXq1eGzz8AjgxtFRDxNoVZEJBkUKQLTptlTER5/HIKC4Isv4J574M474cMPFW5FRFKSQq2ISDIKD4fXX4e9e+3QhowZ7eay+++H8uXtcIf4eLerFBFJexRqRURSQL588MordlhDr16QJQts2gSNGsFtt8E779jjwUREJHko1IqIpKDcuWHkSIiMhOeeg5AQu7nskUdsP+7MmfDXX25XKSLifQq1IiKpIHt2GDQIoqLgpZfghhtgxw547DEoWhSmTIGzZ92uUkTEuxRqRURSUbZs8PzzNtwOGwY5ctj+2w4doFAhmDABTp92u0oREe9RqBURcUHWrNCnj+25HTsW8uSBn3+GLl0gIgLGjIGTJ92uUkTEOxRqRURclCULdO9uV2snTIBbboHffoOePe1JCsOGQUyM21WKiPg+hVoRER+QMSN06gS7d9v+2gIF4MgR6NfPhtuBA+GPP9yuUkTEdynUioj4kPTpoV07u4ls5ky7ieyPP2DAAAgLg2eftWFXREQupVArIuKD0qWDli1hyxZ7pu2tt8KJEzB0qA23vXrZNgUREbGuOdQaY2YYY6qmRDEiInKpwEBo0gQ2b4bFi+H22yE2FkaPtm0JTz1lN5iJiPi7pKzUZgVWGGN2GWOeNcbcnNxFiYjIpQIC4MEH4dtv4cMPoXJlOHMGXnsNChaExx+3JymIiPiraw61juM0Bm4GXgMeBiKNMR8ZYx4yxgQld4EiIvIPY+C++2DtWvjsM6heHc6dgzffhMKFoXVr2LnT7SpFRFJfknpqHcc56jjOK47jlAUqAruBWcABY8xYY0zh5CxSREQuZQzUrAkrV8KaNVCnDsTFwVtvQfHi0KyZHccrIuIvrmujmDEmD1D7/BUHfAiUBLYaY56+/vJERORq/vc/WL4c1q2D+vUhPh7mzbObyxo3hu++c7tCEZGUl5SNYkHGmMbGmGVAFLYFYSyQx3GcxxzHqQ20BF5I3lJFRORKKlaEJUtsiH3oIbua+957UK4cPPCADb0iImlVUlZqDwKTsYG2ouM45R3Hed1xnBMXPeZj4HhyFCgiItemTBlYsMC2HzRvbjeZffABVKoEtWvD6tVuVygikvySEmqfBvI6jtPZcZzNCT3AcZw/HMeJuL7SRETkepQoAW+/Ddu3Q5s29uzbTz6BatXs9emn4DhuVykikjyScvrBLMdxTqdEMSIikvwKF4Zp02DXLnjiCTu1bPVqqFUL7rzTruIq3IqI12mimIiInwgPh0mTYM8e6NYNMmaEb76x/bbly8OiRXaTmYiIFynUioj4mXz5YNw4iIyEZ56BLFlg0yZo1Ahuu82enBAX53aVIiLXRqFWRMRP5coFI0bYcNu/P4SE2M1lzZrZfty33rKDHUREvEChVkTEz2XPDi+/DFFR8NJLcOONdipZ69ZQtChMngxnz7pdpYjIlSnUiogIANmywfPP25Xb4cMhZ07Ytw86doSCBeG11+DUKberFBFJmEKtiIhcImtW6N3bBtqxYyFPHvjlF3jqKShQAEaPhpMn3a5SRORSCrUiIpKgzJmhe3fYuxcmToT8+eG336BXL3uSwtChEBPjdpUiIpZCrYiIXFHGjPDkk/ac26lTbSvCkSPw7LMQFgYDBsAff7hdpYj4O4VaERFJlPTpoW1bO6Fs1iwoVgyOH4eBA224ffZZ+P13t6sUEX+lUCsiItckXTpo0cIe/zV/PpQuDSdO2HaE8HDo2RMOHnS7ShHxNwq1IiKSJIGB8PDD8N13sGSJnUoWGwtjxkBEBHTpAj//7HaVIuIvFGpFROS6BARA/fqwfj189BHceSecOQMTJtj+244d7WYzEZGUpFArIiLJwhi491748kv4/HOoUcNOJJs8GYoUgccegx073K5SRNIqhVoREUlWxthA+/nnNuDeey/ExcHMmVC8uB3D+9NPblcpImmNQq2IiKSYu+6yLQnr18ODD4LjwLx5cOut0KgRbNrkdoUiklYo1IqISIqrUAEWL4bNm+3mMmNg0SK4/XZ44AH45hu3KxQRr/NMqDXGLDXG7DfGnDbGHDTGzDLG5HW7LhERSbzbbrPHgP30kz0WLCAAPvgAKleGWrVg9Wq3KxQRr/JMqAVWAk2AokBjoCDwrqsViYhIkpQoYQc47NhhBzqkSweffgrVqkHVqvDJJ7ZVQUQksYzj0f/XMMbUBxYDGRzHOZfI54QA0dHR0YSEhKRofSIikniRkTBihB3De/asvXfHHdC/P9x/v21XEBH/EBMTQ2hoKECo4zgxiX2el1ZqLzDG3Ag0B766UqA1xmQwxoT8fQFZU61IERFJtPBwmDjRnmfbrRtkzAjr1kG9erbv9r33ID7e7SpFxJd5KtQaY4YbY04CR4H8wINXeUo/IPqi65eUrVBERK7HzTfDuHF25bZ3b8iSxU4sa9zYjuOdO9ceDyYicjlXQ60xZoAxxrnKVf6ip4wEygK1gThgpjFX/EupoUDoRVe+lPpZREQk+eTKBcOHQ1QUPP88hIbCli3w6KP2rNsZM+xgBxGRv7naU2uMyQ5kv8rDIh3HOZ3Ac/MBPwN3Oo7zdSLfTz21IiIedPw4vPYajB0Lx47ZexER0LevnVSWIYO79YlI8klqT62XN4rdAuwHajiOsyqRz1GoFRHxsD//hEmTYNQoOHzY3suXz7YqtG8PmTK5W5+IXL80HWqNMRWBisCXwB9AAeAlIA9Q0nGcM4l8HYVaEZE0IDYWJk+2JyYcOGDv5c4NvXrB449DcLC79YlI0qX10w9OAY2Az4AdwDTgJ6BaYgOtiIikHZkz21MS9u61K7dhYfDbbzbUhofDkCEQk+h/FYpIWuCJldrkopVaEZG06dw5ePttG2Z377b3smWDrl1t+L3xRnfrE5HES+srtSIiIv8pKAjatIFt22y4LV7cbi576SW7ctuv3z89uCKSNinUiohImpEuHTRvDj/9BAsWwG23wYkTMGyYDbc9esDBg25XKSIpQaFWRETSnIAAeOghO7hhyRIoXx5OnbJHgkVEQJcusH+/21WKSHJSqBURkTTLGKhfH9avh+XL4a674MwZmDABChWCDh3sZjMR8T6FWhERSfOMgTp1YM0aWLkSata0m8umTIEiRewAh+3b3a5SRK6HQq2IiPgNY6B6dfjsM1i7Fu67D+LiYOZMKFECHnkEfvzR7SpFJCkUakVExC/deSd8+CF8+y08+CA4DrzzDpQuDQ0bwsaNblcoItdCoVZERPxa+fKweDF8/z00aWJXcxcvtvfvvx++/trtCkUkMRRqRUREsCu077wDW7ZAixb2BIUPP7QruvfcA1984XaFInIlCrUiIiIXKV4cZs2CHTugXTt79u1nn9le3KpVYcUK26ogIr5FoVZERCQBhQrZ0xF274ZOnSB9ent6Qp06UKkSLFumcCviSxRqRUREriAszJ5ru28fdO8OmTLZc2/r1YNy5WDhQoiPd7tKEVGHUhAFAAAgAElEQVSoFRERSYS8ee1EsshI6NMHgoNh82Y7uezWW2HOHHs8mIi4Q6FWRETkGuTMCcOG2XD7/PMQGgpbt0Lz5rYfd8YMO9hBRFKXQq2IiEgS3HQTvPQSREXBoEH26127oE0bO6XsjTfsSF4RSR0KtSIiItchNBSee86u3I4cCbly2T8/8QQULAivvgqnTrldpUjap1ArIiKSDIKDoVcvu6HslVfg5pvh11+hWzeIiIBRo+DPP92uUiTtUqgVERFJRpkyQdeusGcPvP46hIfDoUPwzDP2z4MHQ3S021WKpD0KtSIiIikgQwZ4/HHYuROmT4fCheHoUejf34bbF1+EY8fcrlIk7VCoFRERSUFBQdC6tT0hYfZsKFECjh+3m8zCwqBvXzh82O0qRbxPoVZERCQVpEsHjz4KP/4I774Lt91me2yHD7crtz16wIEDblcp4l0KtSIiIqkoIAAaN4bvvoOlS6FCBXs6wtixUKAAdO4M+/e7XaWI9yjUioiIuMAYO2p33Tr4+GP43//subYTJ9qjwNq3t5vNRCRxFGpFRERcZAzUrg2rV8PKlXD33fDXXzB1KhQtCq1awfbtblcp4vsUakVERHyAMVC9Onz6KXz1Fdx3H8TFwaxZdnNZ06a2H1dEEqZQKyIi4mMqV4YPP4QNG6BBA3AcmD8fSpe2X2/c6HaFIr5HoVZERMRH3X47LFoE339vV2qNgSVLoHx5qFsXvv7a7QpFfIdCrYiIiI8rXRrmzbNn3bZsCYGB8NFHcOedtgd31Sq7mivizxRqRUREPKJYMZg5E3bssKcjBAXB559DjRpQtao9RUHhVvyVQq2IiIjHFCwIkyfD7t3QqZMdyfvll3DvvVCpErz/vsKt+B+FWhEREY/Knx8mTIC9e+HppyFTJli/HurXh7Jl7eSy+Hi3qxRJHQq1IiIiHpc3L4wZA5GR0LcvBAfbzWUPPwylSsHs2fbsW5G0TKFWREQkjciZE4YOhagoeOEFCA2FbdugRQsoXhymT4dz59yuUiRlKNSKiIikMTfeCAMH2nA7eDDcdJPtv23bFgoXhtdftyN5RdIShVoREZE0KjQUnn3WtiWMGgW5ctmg++STdrPZq69CbKzbVYokD4VaERGRNC44GHr2hH37bJC9+Wb49Vfo1g0iImDkSPjzT7erFLk+CrUiIiJ+IlMmeOop2LMH3ngDwsPh8GHo3dv+efBgiI52u0qRpFGoFRER8TMZMkDHjrBzp908VrgwHD0K/ftDWJjdZHbsmNtVilwbhVoRERE/FRQErVvbExLmzIESJexK7csv23Dbp49dyRXxAoVaERERPxcYCM2awY8/2oENZcrYHtsRI2xbwtNP2x5cEV+mUCsiIiIABARA48awaZMdtVuxIpw6BePGQYECdiRvVJTbVYokTKFWRERELmEMPPAAfPMNrFgBVarA2bMwaRIUKgTt2tlzb0V8iUKtiIiIJMgYqFULVq+GVavgnnvsuN1p06BoUWjZ0vbjivgChVoRERG5qmrV4JNP4KuvoG5diI+Ht9+GkiWhSRP44Qe3KxR/57lQa4zJYIzZbIxxjDFl3K5HRETEn1SuDB98ABs3QsOG4DiwYAHcdhs0aAAbNrhdofgrz4VaYARwwO0iRERE/Fm5cvDee3aFtmlT26qwZAlUqAD33WdXdEVSk6dCrTHmPqA20MvtWkRERARuvRXmzbO9ta1a2ePBli+Hu+6CmjVh5Uq7miuS0jwTao0xuYDJQEsgNpHPyWCMCfn7ArKmZI0iIiL+qmhReOstO6WsQwc72GHlShtsq1SxQVfhVlKSJ0KtMcYAM4DXHce5lm6dfkD0RdcvyV+diIiI/K1AAXjzTXvkV+fOdiTv2rW2JeGOO2DpUoVbSRmuhlpjzIDzG76udJUHngJCgKHX+BZDgdCLrnzJ+xOIiIhIQvLnh9deg337oEcPyJwZvv0WHnzQTixbsMCeoCCSXIzj4n8uGWOyA9mv8rBIYB5QD7i42EAgDpjtOM5jiXy/ECA6OjqakJCQay9YREREkuT332HsWBt0T5yw94oXh+eesxvN0qVztz7xHTExMYSGhgKEOo4Tk9jnuRpqE8sYkx+7Uvu3vMDHwEPAOsdxEtVWoFArIiLirmPH4NVX4ZVX4Phxe69QIejXzw5zCApytz5xX5oOtZczxoQD+4CyjuNsvobnKdSKiIj4gJgYmDABxoyBI0fsvbAw6NsX2rSxvbjin5Iaaj2xUUxERETSlpAQuzobGQmjRkGuXBAVBU8+aTebvfIKxCbqrCMRy5MrtUmllVoRERHfdOoUTJ0Kw4fDL+ebCnPmhJ49bdDNqkM5/YZWakVERMSzMmWCLl3sUWBvvAHh4XD4MPTpY/88aNA/PbgiCVGoFREREZ+RIQN07GiHOMyYAUWK2M1lzz9ve26ffx6OHnW7SvFFCrUiIiLic4KC4LHHYOtWmDsXSpa0m8sGDbLhtndvOHTI7SrFlyjUioiIiM8KDIRHHoEffoCFC6FsWTh5EkaOhIgI6N4dfv3V7SrFFyjUioiIiM8LCIBGjWDjRli2zI7cPXXKnpJQoIDdTBYV5XaV4iaFWhEREfEMY+D+++Hrr+GTT6BqVTh7Fl5/3Q5xaNfObjYT/6NQKyIiIp5jDNxzD3zxhb3uuQf++gumTYOiRaFFC9i2ze0qJTUp1IqIiIinVa1qV22//tqu4sbHw+zZdnPZww/D99+7XaGkBg1fSEBcXBznzp1LvcIkRQQFBREYGOh2GSIiksq++86ekvDee//cq18f+veHChXcq0sSJ6nDFxRqL+I4Dr/99hvHdbpzmpEtWzZy586NMcbtUkREJJX99BMMHgzvvAN/x506dexZt3fd5W5t8t8UahPhaqH24MGDHD9+nJw5c5I5c2YFIQ9zHIfY2FgOHz5MtmzZyJMnj9sliYiIS3bsgKFD4e23IS7O3qte3YbbGjVsf674DoXaRLhSqI2Li2Pnzp3kzJmTm266yZ0CJdkdPXqUw4cPU6RIEbUiiIj4ub17YdgwO6ns7y7DO++04bZOHYVbX5HUUKuNYuf93UObOXNmlyuR5PT356keaRERKVAA3nwT9uyBLl3sSN6vvoL77oOKFWHJkn/aFMR7FGovo5aDtEWfp4iIXO6WW2D8eNi3D3r2hMyZYcMGaNAAypSBBQv+aVMQ71CoFREREb+UJw+MGgWRkdCvH2TNasfxNmkCpUrZHty//nK7SkkshVqPq169Ot27d7/wdXh4OOPGjXOxon9MnTqV2rVrJ/rxy5Yto2zZssTHx6dgVSIiIpfKkQOGDLFjdgcMgGzZYPt2aNkSihWDqVPt1DLxbQq1acy3335Lx44d3S6DM2fO8MILL/D8888n+jkPPPAAxhjmzJmTgpWJiIgk7IYb4MUXbbgdOhSyZ7f9t+3bQ+HCMGkSnD7tdpXyXxRq05gcOXL4xGa3hQsXEhwcTJUqVa7peW3atGH8+PEpVJWIiMjVhYRA3762LWH0aMidG/bvh06doGBBGDcOYmPdrlIup1D7HxwHTp5057qenZeXtx8YY5gyZQoNGzYkc+bMFC5cmKVLl17ynK1bt1K3bl2Cg4PJlSsXLVu25MiRIwm+/smTJwkJCeHdd9+95P77779PlixZOHHiBADz5s2jfv36F75/+vRpSpYseckq8r59+wgNDWXy5MkX7tWvX5/169ezd+/epP8SREREkkGWLNCjh91Q9tprdoPZgQPw9NMQHg7Dh8P5f+2JD1Co/Q+xsRAc7M6V3P/1N3DgQJo0acIPP/xA3bp1ad68OceOHQPswIlq1apRpkwZNmzYwPLlyzl06BBNmjRJ8LWyZMnCI488wvTp0y+5P336dB566CGyZs0KwJo1ayhfvvyF72fMmJHZs2fz1ltvsXjxYuLi4mjZsiU1atSgQ4cOFx4XFhZGzpw5WbNmTfL+EkRERJIoY0bo3Bl277ZHgkVEwO+/29Xc8HB4+WXQMFL3KdT6gdatW9OsWTMKFSrEkCFDOHnyJOvXrwdg0qRJlCtXjiFDhlCsWDHKli3LtGnTWLlyJTt37kzw9dq3b8/HH3/MgQMHADhy5AjLli2jbdu2ABw/fpzjx4+TN2/eS55XpkwZBg0aRIcOHXj66afZs2cPU6ZM+dfr33zzzURGRibjb0BEROT6pU8PHTrAzp3w1ltQpAgcOwYvvABhYdC/P/zHX3RKKlCo/Q+ZM8Off7pzJXdLbOnSpS/8OUuWLGTNmpXDhw8DsHHjRlauXElwcPCFq1ixYgDs2bMnwderWLEiJUuWZObMmQDMmjWL/PnzU7VqVQBOnToF2NXZy/Xs2ZOiRYsyfvx4pk+fTvbs2f/1mEyZMhGrZiUREfFR6dJBq1awdSvMm2eP/4qJgcGD7cpt795w6JDbVfofhdr/YIztpXHjSu55AUFBQZf9bObCsVnx8fHUq1ePzZs3X3Lt2rXrQkhNSPv27S+0IEyfPp02bdpcGHRw0003YYzhjz/++NfzDh8+zI4dOwgMDGTXrl0JvvaxY8fIkSNHkn5WERGR1BIYCE2bwvffw3vvQblydm/MyJE23HbrBr/+6naV/kOh1s+VK1eOLVu2EB4eTqFChS65smTJ8p/Pa9GiBfv37+fVV19ly5YtPPbYYxe+lz59ekqUKMHWrVv/9by2bdtSqlQpZs6cSe/evf/1mNOnT7Nnzx7Kli2bfD+kiIhICgoIgIYN7VSyDz6ASpXs0V+vvmpH8z7xhD1JQVKWQq2f69y5M8eOHaNZs2YXTh1YsWIFbdu2Je4KMwJvuOEGGjVqxDPPPEPt2rXJly/fJd+vU6cOX3755SX3JkyYwNdff83MmTN59NFHeeihh2jevDlnLzrR+ptvviFDhgxUrlw5eX9QERGRFGYM1K0LX30Fn34K1arZoQ1vvGHPuW3bFv7jLyklGSjU+rm8efOydu1a4uLiqFOnDqVKlaJbt26EhoYSEHDl/3m0a9eOs2fPXtggdrEOHTrw4YcfEh0dDcD27dt55plnmDhxIrfccgtgQ+7x48cvGdAwd+5cmjdv7hNn7YqIiCSFMXD33bBqFaxeDbVq2XG706fbCWXNm9t+XElexrmeQ1E9xhgTAkRHR0cTEhJyyfdOnz7Nvn37iIiISHCDk/zb7Nmz6datGwcOHCB9+vT/+n6TJk0oW7Ys/fr1S9Tr/f777xQrVowNGzYQERGRLDXqcxUREV+wbh0MGgTLltmvjYFGjeyJCWXKuFubr4mJiSE0NBQg1HGcmMQ+Tyu1cs1iY2PZsmULQ4cO5fHHH08w0AKMHDmS4ODgRL/uvn37mDhxYrIFWhEREV9xxx3w/vuwaZMNs44DCxdC2bJQvz6cP2lTroNCrVyzESNGUKZMGXLlynXFVdiwsDCeeuqpRL9uxYoVadq0aXKUKCIi4pPKlrVh9qefoFkzu8ns/fdt6K1TBy7bjiLXQKFWrtmAAQM4d+4cn3322TWtxIqIiIhVsiTMmQPbtkHr1vZ4sBUroEoVqFEDPv/cruZK4inUioiIiLikSBG7gWzXLujYEYKC7Aazu++Gu+6Cjz5SuE0shVoRERERl0VE2KO/9u6Fp56CjBnh66/tEWEVKsDixXB+bpL8B4VaERERER+RL58d2rBvH/TqBZkzw8aNdrhDmTIwfz5c4Rh5v6ZQKyIiIuJjcue243ajouDZZyFrVvjxRzuWt1QpmDXLnn0r/1CoFREREfFR2bPD4ME23A4cCDfcANu3Q6tWULQoTJlip5aJQq2IiIiIz7vhBnjhBYiMhKFDbdjduxc6dLAjeCdOhNOn3a7SXQq1Hle9enW6d+9+4evw8HDGjRuX6nXMmDGDbNmyJfvrtm7dmgYNGiT764qIiHhRSAj07WvD7ZgxkCcP7N8PnTtDgQIwdizExrpdpTsUatOYb7/9lo4dO7pdRqJFRUWRIUMGYmISPQVPRETE72XJAk8/bVdrJ0yAW26BgwehRw8ID4fhw+HECberTF0KtWlMjhw5yJw5s9tlJNqSJUuoXr06ISEhbpciIiLiORkzQqdOsHs3TJ5sV2t//92u5oaFwUsvwfHjbleZOhRq/4PjOJyNc+dyruOU5cvbD4wxTJkyhYYNG5I5c2YKFy7M0qVLL3nO1q1bqVu3LsHBweTKlYuWLVty5MiRK77PjBkzyJ8/P5kzZ6Zhw4YcPXr0wvciIyMJCAhgw4YNlzxn/PjxhIWFXfLzLVmyhPr16wMQFxdHjx49yJYtGzfddBO9e/e+5LG///47uXPnZsiQIRfurVu3jvTp07NixYpr+C2JiIikLenTQ/v2sGMHzJxpN5H98Qe8+KINt889B1f5V7vnpXO7AF91Lh7G/HD06g9MAT1K30T6wOR7vYEDBzJixAhGjhzJ+PHjad68OVFRUdx4440cPHiQatWq0aFDB8aMGcOpU6fo06cPTZo04fPPP0/w9datW0fbtm0ZMmQIjRo1Yvny5bz44osXvh8eHs4999zD9OnTKV++/IX706dPp3Xr1hhjADh+/Dhr1qxhxowZAIwePZpp06YxdepUSpQowejRo1m0aBE1a9YE7Cr0tGnTaNCgAbVr16ZYsWK0aNGCTp06Ubt27eT7hYmIiHhUunTQsiU8+igsXAiDBtmjwIYMgVdegSefhJ497ZFhaY1Wav1A69atadasGYUKFWLIkCGcPHmS9evXAzBp0iTKlSvHkCFDKFasGGXLlmXatGmsXLmSnTt3Jvh6r7zyCnXq1KFv374UKVKErl27UqdOnUse0759e+bOncuZM2cA+P7779m8eTNt2rS58JgPP/yQW2+9lVtuuQWAcePG0a9fPxo3bkzx4sV5/fXXCQ0NveR169atS4cOHWjevDlPPPEEGTNmZNiwYcn2uxIREUkLAgOhSRPYvBkWLYJy5eDkSRg1yk4v69oVfvnF7SqTl1Zq/0NQgF0xdeu9k1Pp0qUv/DlLlixkzZqVw4cPA7Bx40ZWrlxJcHDwv563Z88eihQp8q/727Zto2HDhpfcq1y5MsuXL7/wdYMGDejSpQuLFi3ikUceYdq0adSoUYPw8PALj7m49SA6OpqDBw9SuXLlC99Ply4d5cuX/1c7xqhRoyhVqhTz589nw4YNZMyY8Rp+GyIiIv4jIAAaNIAHH4Tly+Hll+343fHj4fXXoU0b238bEeF2pddPK7X/wRhD+kB3rr//ej65BAUF/etniz8/QDo+Pp569eqxefPmS65du3ZRtWrVBF8vMT2/6dOnp2XLlkyfPp2zZ88yZ84c2rZte+H7586dY/ny5Tz44IPX/PPs3buXAwcOEB8fT1RU1DU/X0RExN8YA/fdB2vXwqefQrVqcO4cvPmmPee2TRvYtcvtKq+PZ0KtMSbSGONcdunvna9TuXLl2LJlC+Hh4RQqVOiSK0uWLAk+p0SJEnzzzTeX3Lv8a7AtCJ9++ikTJ07k3LlzNGrU6ML3Vq5cSbZs2ShTpgwAoaGh5MmT55LX+euvv9i4ceMlr3n27FmaN29O06ZNGTRoEO3atePQoUNJ/vlFRET8iTFw992wahWsXg21a0NcHMyYAcWK2V7cLVvcrjJpPBNqz3sByHPRNcjdcryvc+fOHDt2jGbNmrF+/Xr27t3LihUraNu2LXFxcQk+p2vXrixfvpwRI0awc+dOXnvttUtaD/5WvHhxKlWqRJ8+fWjWrBmZMmW68L2lS5deaD34W7du3Rg2bBiLFi1i+/btdOrUieOXnUPy3HPPER0dzauvvkrv3r0pXrw47dq1S4bfhIiIiH+pUgU+/hjWrYN69SA+HubOhVKl4KGH4Lvv3K7w2ngt1J5wHOe3i64/3S7I6/LmzcvatWuJi4ujTp06lCpVim7duhEaGkpAQML/86hUqRJTpkxh/PjxlClThhUrVtC/f/8EH9uuXTvOnj17SesB2FB7eetBz549adWqFa1bt6Zy5cpkzZr1kt7dVatWMW7cOGbNmkVISAgBAQHMmjWLL7/8kkmTJl3nb0JERMQ/VawIS5faENu4sb23cKHdXFavHpzfW+7zzPWciZqajDGRQAYgPfAzsAAY6TjO2Ss8J8P55/wtK/BLdHT0vw77P336NPv27SMiIkIbj5LR4MGDmTdvHj/++OOFe5s2baJmzZr8/vvv/+r3TW76XEVERK7Nli32CLB58+zqLcDDD8P8+anz/jExMX+ffhTqOE6iR456aaX2FeARoAbwGtAdmHiV5/QDoi+60tjhFb7rzz//5Ntvv2X8+PF07dr1ku/99ddfjB8/PsUDrYiIiFy7kiVh9mzYvt1uIEuXzvbb+jpXV2qNMQOAF6/ysAqO42y4/KYxpjHwLpDdcZwEpyRopdY9rVu3Zu7cuTRo0IA5c+YQGJiM0ySugT5XERGR6xMZCaGhcMMNqfN+SV2pdTvUZgeyX+VhkY7jnE7guTdjV14rOY6zLpHvFwJEK9T6D32uIiIi3pLUUOvq8AXHcY4ASZ1EXPb8Pw8mUzkiIiIi4lGemChmjKkMVAJWYntjKwBjgaWO4+xPzvfyysY5SRx9niIiIv7BE6EWOAM0xfbfZgCigMnAiOR6g783LcXGxl5ynqp4W2xsLPDvqWoiIiKStngi1DqOswm7UptiAgMDyZYtG4cPHwYgc+bMyT6uVlKP4zjExsZy+PBhsmXL5tpGNREREUkdngi1qSV37twAF4KteF+2bNkufK4iIiKSdinUXsQYQ548eciZMyfnzp1zuxy5TkFBQVqhFRER8RMKtQkIDAxUGBIRERHxEC9NFBMRERERSZBCrYiIiIh4nkKtiIiIiHieX/bUxsQkeuKaiIiIiKSipOY0408Tl4wxNwO/uF2HiIiIiFxVPsdxfk3sg/0t1BogL3Aild4yKzZE50vF95Tkpc/Q2/T5eZ8+Q+/TZ+h9bnyGWYEDzjUEVb9qPzj/i0l04r9eF00kO+E4jnoePEifobfp8/M+fYbep8/Q+1z6DK/5fbRRTEREREQ8T6FWRERERDxPoTZlnQEGnv+neJM+Q2/T5+d9+gy9T5+h93niM/SrjWIiIiIikjZppVZEREREPE+hVkREREQ8T6FWRERERDxPoVZEREREPE+hNoUYYzoZY/YZY04bYzYaY6q4XZMkjjGmnzHmW2PMCWPMYWPMYmNMUbfrkqQ7/5k6xphxbtciiWeMudkY87Yx5qgxJtYYs9kYc7vbdcnVGWPSGWMGnf/34CljzF5jzAvGGOUOH2WMqWqMed8Yc+D8/182uOz7xhgz4Pz3TxljVhljSrpVb0L0P64UYIxpCowDBgNlgTXAR8aY/K4WJolVDZgAVAJqYSfvrTDGZHG1KkkSY0wFoCPwg9u1SOIZY24A1gLngPuAEkBP4LibdUmi9QGeALoAxYHewDPAU24WJVeUBfge+5klpDfQ4/z3KwC/AZ8YY7KmTnlXpyO9UoAxZh2wyXGcJy+6tw1Y7DhOP/cqk6QwxuQADgPVHMdZ7XY9knjGmGBgE9AJ6A9sdhynu7tVSWIYY4YBdzmOo7/l8iBjzDLgkOM47S66txCIdRynpXuVSWIYYxygoeM4i89/bYADwDjHcYafv5cBOAT0cRzn/+3cS4iVZRzH8e+/wqhNBBUWmCjlzrTsAiImXVyEFF2giIRaVEy1sHaNi6CgC4LdFIogsiijNkHYQgoLzKhGjLISXEgSplBC2UWF8d/ieSemcYQ50zjPeeL7gcPhvGfzW8w75/f+3+d9Xq4WdhQntVMsImYAi4DNY77aDCye/kSaAmd17werptBkrAc2ZeaHtYOoZzcCQxHxbrcMaEdE3Fs7lCZsK3BtRMwDiIgFwBLgg6qpNFlzgJmM6jaZeQT4hD7qNqfVDvA/dA5wKuXqZbQDlD8INaS7Ol0LbM3MnbXzaOIi4g7gMsptMrVnLjBAOf+eBK4EXoiII5n5etVkmohnKAOBXRExTPldXJ2ZG+vG0iSN9Jfxus3sac5yQpbak2fsuo4Y55j63zrgEsqEQY2IiFnA88DyzDxcO48m5RRgKDMHu887uodSBgBLbf+7HbgLuBP4FlgIPBcR+zJzQ9Vk+i/6uttYaqfez8Awx09lz+P4Kxz1sYh4kXILdGlm/lg7j3qyiHLObS/DdqBMipZGxEPA6Zk5XCucJuQn4Lsxx74Hbq2QRb1bAzydmW93n7+JiNnAo4Cltj37u/eZlHNzRF91G9fUTrHMPApspzw1P9r1wLbpT6RedduWrANuAa7JzD21M6lnHwHzKdOhkdcQ8Caw0ELbhE+BsVvpzQN+qJBFvTsTODbm2DD2jlbtoRTbf7pN9wzR1fRRt3FSe3KsBd6IiCHgM8p2QhcCL1VNpYlaT7lldhNwKCJGpu6/ZuZf9WJpojLzEPCvNdAR8Qfwi2ujm/EssC0iBoF3KGtq7+te6n/vA6sjYi9l+cGllO2gXq2aSifU7RZz0ahDcyJiIXAwM/d2+3wPRsRuYDcwCPwJvDX9acfnll4nSUQ8QNnT7XzKj+vDbgfVhm4rk/Hck5mvTWcWTZ2I+Bi39GpKRKwAngIupkyK1mbmK3VTaSK6vUufAG6m3KLeB2wEHu/uaKrPRMQyYMs4X23IzLu7B6cfA+4HzgY+Bx7sp0GBpVaSJEnNc22LJEmSmmeplSRJUvMstZIkSWqepVaSJEnNs9RKkiSpeZZaSZIkNc9SK0mSpOZZaiVJktQ8S60kSZKaZ6mVJElS8yy1kiRJap6lVpIaExHnRsT+iBgcdeyqiDgaEctrZpOkWiIza2eQJPUoIm4A3gMWA7uAHcCmzFxVNZgkVWKplaRGRcR64DrgS2ABcEVmHq6bSpLqsNRKUqMi4gxgJzALuDwzv64cSZKqcU2tJLVrLnAB5X/57MpZJKkqJ7WS1KCImAF8AXxFWVP7CDA/Mw9UDSZJlVhqJalBEbEGuJcnBncAAACCSURBVI2ylvZ3YAtwKDNXVA0mSZW4/ECSGhMRy4BVwMrM/C0zjwErgSURMVA1nCRV4qRWkiRJzXNSK0mSpOZZaiVJktQ8S60kSZKaZ6mVJElS8yy1kiRJap6lVpIkSc2z1EqSJKl5llpJkiQ1z1IrSZKk5llqJUmS1DxLrSRJkpr3NxO+YL8FJ5PBAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11912e4e0>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Delta t: 2.2516269000305043 s\n"
     ]
    }
   ],
   "source": [
    "xaxis = np.linspace(XBEGIN, XFINAL, 1000)\n",
    "thetas = cheby_line_coefficients(TERMS)\n",
    "thetas = np.array(thetas).reshape(-1,1)\n",
    "func = get_cheby_path_func(thetas)\n",
    "deriv_func = get_cheby_path_func_deriv(thetas)\n",
    "plt.figure(figsize=(8,4), dpi=100)\n",
    "y_points, = plt.plot(xaxis, func(xaxis), color=\"blue\",label='line y(x)')\n",
    "yp_points, = plt.plot(xaxis, deriv_func(xaxis), color=\"skyblue\",label='line dy/dx')\n",
    "plt.legend(handles=[y_points, yp_points], loc=3)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.show()\n",
    "print('Delta t:',delta_time(func,deriv_func),'s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Evolution Strategies and the Brachistochrone\n",
    "\n",
    "Rather than using the Euler-Lagrange equation to solve the problem, we can use ML to find the optimal curve. Essentially we will try to minimize the $\\Delta t$ function by crawing around in fourier space until we find the coefficients that return the minimum function.\n",
    "\n",
    "A random set of fourier coefficients will create a random curve. This curve may go above the $y=0$ line, which violates conservation of energy. A random curve will also likely not respect the boundry conditions of $a = (0,0)$ and $b  = (x_f,y_f)$. If we want to use an ML algorithm to find the brachistochrone curve, we need to create a reward function that uses the $\\Delta t$ function as a basis, and then has extra terms that heavily penalizes curves that do not follow the boundry conditions and conservation of energy. \n",
    "\n",
    "We can then use this new reward function as the fit function for the CMA-ES or OpenAI-ES algorithms. Fourier coefficients that return a lower $\\Delta t$ and respect the laws of nature will recieve a higher reward and pass on their coefficients to the next generation. \n",
    "\n",
    "For the case where we are using few fourier coefficients, the approximation of the cycloid may not be optimal, and the ES alogorithms may very well do better than our benchmark."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cost_function(c_list):\n",
    "    c_list = np.copy(c_list)\n",
    "    def cost(x):\n",
    "        #equation gotten from conservation of energy\n",
    "        #due to conservation of energy, values of f(x) > 0 are impossible.\n",
    "        #they will create a negative inside the sqrt and return an NaN\n",
    "        return np.sqrt((1+fourier_yp(c_list,x)**2)/(2.0*g*fourier_y(c_list,x)))\n",
    "    time = integrate.quad(cost,0.0001,XFINAL,limit=200)[0] #don't start integration from zero-- messes up integrator.\n",
    "    #we also want to constrain the path to begin at (XBEGIN,YBEGIN) and end at (XFINAL,YFINAL)\n",
    "    #i.e. F(XBEGIN;theta) - YBEGIN = 0 and F(XFINAL;theta) - YFINAL = 0\n",
    "    begin_penalty = math.fabs(fourier_y(c_list,XBEGIN)-YBEGIN)**3\n",
    "    end_penalty = math.fabs(fourier_y(c_list,XFINAL)-YFINAL)**3\n",
    "    return -1.*time-begin_penalty-end_penalty\n",
    "\n",
    "fit_func = cost_function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NPARAMS = (FOURIER_TERMS*2)+1 #Our parameters are all our fourier coefficients\n",
    "NPOPULATION = 50    # use population size of 50.\n",
    "MAX_ITERATION = 200 # run each solver for 300 generations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "rand_coefs = np.random.randn(NPARAMS)\n",
    "yaxis_rand = fourier_y(rand_coefs,xaxis)\n",
    "plt.figure(figsize=(8,4), dpi=100)\n",
    "rand_y, = plt.plot(xaxis, yaxis_rand, color=\"green\",label='rand y(x)')\n",
    "line_y, = plt.plot(xaxis, yaxis_line, color=\"blue\",label='line y(x)')\n",
    "optimal_y, = plt.plot(xaxis, yaxis_cycloid, color=\"red\",label='cycloid y(x)')\n",
    "plt.legend(handles=[optimal_y,line_y, rand_y], loc=3)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.show()\n",
    "\n",
    "print('CYCLOID:')\n",
    "print('  fitness:',cost_function(cycloid_coefs))\n",
    "print('LINE:')\n",
    "print('  fitness:',cost_function(line_coefs))\n",
    "print('RANDOM:')\n",
    "print('  fitness:',cost_function(rand_coefs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# defines a function to use solver to solve fit_func\n",
    "def test_solver(solver):\n",
    "    history = []\n",
    "    for j in range(MAX_ITERATION):\n",
    "        solutions = solver.ask()\n",
    "        fitness_list = np.zeros(solver.popsize)\n",
    "        for i in range(solver.popsize):\n",
    "            fitness_list[i] = fit_func(solutions[i])\n",
    "        solver.tell(fitness_list)\n",
    "        result = solver.result() # first element is the best solution, second element is the best fitness\n",
    "        history.append((result[0],result[1]))\n",
    "        if (j) % 2 == 0:\n",
    "            print(\"fitness at iteration\", (j), result[1])\n",
    "    print(\"fitness score at this local optimum:\", result[1])\n",
    "    return history"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# defines CMA-ES algorithm solver\n",
    "cmaes = CMAES(NPARAMS,\n",
    "              popsize=NPOPULATION,\n",
    "              weight_decay=0.0,\n",
    "              sigma_init = 0.2\n",
    "          )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cma_history = test_solver(cmaes)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cma_coefs = cma_history[-1][0]\n",
    "plt.figure(figsize=(16,8), dpi=150)\n",
    "for i in range(len(cma_coefs)-1):\n",
    "    yaxis_temp = fourier_y(cma_history[i][0],xaxis)\n",
    "    plt.plot(xaxis, yaxis_temp, color=(0.0, 1.0-(i/len(cma_coefs)), (i/len(cma_coefs))),alpha=0.05+0.15*(i/len(cma_coefs)))\n",
    "yaxis_cma = fourier_y(cma_coefs,xaxis)\n",
    "cycloid_y, = plt.plot(xaxis, yaxis_cycloid, color=\"red\",label='cycloid y(x)')\n",
    "cma_y, = plt.plot(xaxis, yaxis_cma, color=\"orange\",label='CMA-ES y(x)')\n",
    "plt.legend(handles=[cycloid_y,cma_y], loc=3)\n",
    "plt.show()\n",
    "\n",
    "print('CYCLOID:')\n",
    "print('  fitness:',cost_function(cycloid_coefs))\n",
    "print('CMA-ES:')\n",
    "print('  fitness:',cost_function(cma_coefs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# defines OpenAI's ES algorithm solver. Note that we needed to anneal the sigma parameter\n",
    "oes = OpenES(NPARAMS,                  # number of model parameters\n",
    "            sigma_init=0.5,            # initial standard deviation\n",
    "            sigma_decay=0.999,         # don't anneal standard deviation\n",
    "            learning_rate=0.1,         # learning rate for standard deviation\n",
    "            learning_rate_decay = 1.0, # annealing the learning rate\n",
    "            popsize=NPOPULATION,       # population size\n",
    "            antithetic=False,          # whether to use antithetic sampling\n",
    "            weight_decay=0.00,         # weight decay coefficient\n",
    "            rank_fitness=False,        # use rank rather than fitness numbers\n",
    "            forget_best=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "oes_history = test_solver(oes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "oes_coefs = oes_history[-1][0]\n",
    "plt.figure(figsize=(16,8), dpi=150)\n",
    "for i in range(len(oes_coefs)-1):\n",
    "    yaxis_temp = fourier_y(oes_history[i][0],xaxis)\n",
    "    plt.plot(xaxis, yaxis_temp, color=(0.0, 1.0-(i/len(cma_coefs)), (i/len(cma_coefs))),alpha=0.05+0.15*(i/len(cma_coefs)))\n",
    "yaxis_oes = fourier_y(oes_coefs,xaxis)\n",
    "cycloid_y, = plt.plot(xaxis, yaxis_cycloid, color=\"red\",label='cycloid y(x)')\n",
    "oes_y, = plt.plot(xaxis, yaxis_oes, color=\"orange\",label='OpenAI-ES y(x)')\n",
    "plt.legend(handles=[cycloid_y,oes_y], loc=3)\n",
    "plt.show()\n",
    "\n",
    "print('CYCLOID:')\n",
    "print('  fitness:',cost_function(cycloid_coefs))\n",
    "print('OpenAI-ES:')\n",
    "print('  fitness:',cost_function(oes_coefs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "cycloid_fitness = [-cost_function(cycloid_coefs)] * MAX_ITERATION\n",
    "plt.figure(figsize=(16,8), dpi=150)\n",
    "\n",
    "oes_fitness = [-i[1] for i in oes_history ]\n",
    "cma_fitness = [-i[1] for i in cma_history ]\n",
    "\n",
    "cycloid_line, = plt.semilogy(cycloid_fitness, color=\"black\", linewidth=0.5, linestyle=\"-.\", label='Cycloid Fit')\n",
    "oes_line, = plt.semilogy(oes_fitness, color=\"green\", linewidth=1.0, linestyle=\"-\", label='OpenAI-ES')\n",
    "cma_line, = plt.semilogy(cma_fitness, color=\"blue\", linewidth=1.0, linestyle=\"-\", label='CMA-ES')\n",
    "\n",
    "plt.legend(handles=[cycloid_line, cma_line, oes_line], loc=1)\n",
    "\n",
    "# Set x limits\n",
    "#plt.xlim(0,2500)\n",
    "\n",
    "plt.xlabel('generation')\n",
    "plt.ylabel('-1*fitness')\n",
    "\n",
    "# plt.savefig(\"./rastrigin_10d.svg\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}