feynman-scratchpad.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "253c51d2-080c-4d65-bbac-4779e75a0ff2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29361d53-f583-47e3-a9b8-d7b2cd695b70",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pylab as pylab\n",
    "import warnings\n",
    "import matplotlib.cbook\n",
    "\n",
    "from plot_params import params\n",
    "\n",
    "warnings.filterwarnings(\"ignore\",category=matplotlib.cbook.mplDeprecation)\n",
    "\n",
    "pylab.rcParams.update(params)\n",
    "cols_default = plt.rcParams['axes.prop_cycle'].by_key()['color']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c77527c1-7515-4295-aecf-f5d71ebfb0af",
   "metadata": {},
   "source": [
    "## scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7344d3ff-8ea3-4ef1-854e-d693fc48b385",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import spence\n",
    "from scipy.optimize import minimize"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a7e4028-8249-4afd-aa6d-cd32c536e826",
   "metadata": {},
   "source": [
    "$$f(x) = \\left( \\sum_i p_i x^i\\right) \\mathrm{Li}_2\\left(\\sum_i q_i x^i\\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "004f7250-5d29-4d97-8d8f-b48a970a997d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def function_target(c, x):\n",
    "    \"\"\" Target function, f(x) = P(x) Li2(Q(x)) where P and Q and polynomials with coefficients c\n",
    "    \"\"\"\n",
    "    p, q = np.split(c, 2)\n",
    "    P = np.polyval(p, x)\n",
    "    Q = np.polyval(q, x)\n",
    "        \n",
    "    return P * spence(1 - (Q + 0.j))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "db0d7389-4a75-4bbf-bc4f-192b67969a4b",
   "metadata": {},
   "outputs": [],
   "source": [
    "degree = 3  # Degree of polynomials\n",
    "\n",
    "c0 = np.random.rand(int(2 * degree))  # Randomly initialize coefficients\n",
    "x = np.linspace(-10, 1, 20)\n",
    "\n",
    "function_true = function_target(c0, x)  # True target function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "de7e3165-b59f-417e-8416-a9ff1e2b28c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import minimize\n",
    "\n",
    "def mse(c, x):\n",
    "    \"\"\" Mean squared error loss\n",
    "    \"\"\"\n",
    "    return np.mean(np.abs(function_target(c, x) - function_true) ** 2)\n",
    "\n",
    "opt = minimize(lambda c: mse(c, x), x0=np.random.rand(int(2 * degree)), method='SLSQP', options={'ftol': 1e-8})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "f147deac-8642-4054-a58a-c0286eaaf763",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, '$q_i$')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlcAAAEdCAYAAAAhEEITAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABNfUlEQVR4nO3deVxU9f4/8Ndn2ARFh8UF3MGtcimkrGxVtH2HzMpKS7j3d5fq3pJLlpK5BLfbvd/vvffbZSo1LUtBb2qlKZVpViZiuZS5jDu4wiiKAsN8fn/MGRxxGAaY4Zwz83o+HjxylnPmzQyc3nzO53xeQkoJIiIiIvIOg9oFEBEREfkTNldEREREXsTmioiIiMiL2FwREREReRGbKyIiIiIvYnNFRERE5EXBnjxJCJECIE1KmeHh89MBlCk3o6WUpmbWR0RERKQrbpsrIUQSgDHKzWhPdqg0VmYpZaFyO1UIkc4Gi4iIiAKB8GQRUSFEKoAxUso0D567SUo5tLH7iIiIiPyRV+dcCSGMABJcPJSkPEZERETk1zyac9UEybgw18qZRXms0NVGQohKACEAjnu5HiLSro4AaqSUEWoX0hI8fhEFJLfHL283V0bYG6n6XDVczkIMBkNw586d4zx5kYqKCkRGRjaxNM9VVlYiIsI3x3vu+1K+/Dz1+p7odd9N+SyPHj0Km83mkzpaWZOOXwB/5v1p3/ws/WffXj1+SSkb/QKQCiDfw+ftcXH/HgApbrYradeunQTg8mvq1KnSmb1s3xk6dCj33Yr79uXnqdf3RK/7dvVZTp06tcHfbQAV0oNjkJa/AJTExcW1+H3yFr3+7Oh13/ws/WffTfks4+LiJIAS2cBxwdvrXFkauD8ajYxeRUZGNnjwys7O9nKZRNRasrOzXf5ex8XFAUCF2vUREXmbt5urIrhessEIwOzl1/KZ9PR07rsV9+1Len1P9LpvUp9ef3b0um9f0ut7otd9e5MvlmLYI6VMbOy+eo+XxMXFxZWUlHhSM4QQ8KRu0gd+nv6jKZ9lfHw8SktLS6WU8T4uy6eaevxStuHPvJ/gZ+k/vHn8atHIlRAiQQiRX+/uPGUhUcdz0gHktOR1iIi0rLKyEsnJyUhOTobJxPWSifyVyWRCcnIyTpw4AQANzqx3O3IlhEgAkAEgBfb1q0wANkopC5THUwDkA+gtpbQ4bZcO+2lAIzyIv+HIVWDj5+k/OHLl8Tb8mfcT/Cz9hzePX26XYpBSmgFkunm8EECUi/t9+qfb1KlTfbl7amX8PP0HP0vP8H3yH/ws/Yc3P0uP5lz5WnP+8iMibTmyeAkiEhPQ/sorPXp+II9cEZG2nNmxAxVbt6LzfffBEBra6PN9OueKiAgAdm7ejs+WrsGOV6aidOEitcshIvKYtNlw4D8mHJ43H788/ydIq7XF+/T2Cu1EFGBsNhtmF6zFvvgh2BLVE68Nu0btkoiIPHay8AtU7tkNADBedx1EcMtbIzZXRNQiaz4uxL5a+zB6n64xiOzVS92CXFAusnEsZNzoRTbKNpOUf1oA388lJaLWd7rMgt3vf4hgAKExseiS+pBX9svTgkTUbLXnq9Duk3wkndyDdrZqjHv6PrVLuoTj6mUpZYFypXOZ83IxDWyTp2yTqzRVRmW9PyLyIx/MWYr/dL4am6N7I378Uwhq08Yr+2VzRUTNdnTJEhiOH8WII9uQfXNvdIi55OJhLchQrmwGACgNVkZDTxZCGAE87FhyRmECkOWzComo1Zm37cS6YzU4HxSKXd0GIPrGG7y2bzZXRNQsVceO4UjBYgBAREICut4+SuWKLqU0SgkuHkpSHnMlGfWyUJV1/JK8WRsRqcdms2HOwi9hEwKQEhMeuBEGg/daIjZXGlBYWAghBDIyMpCZmYnMzMyLbmdkZCAqKgq5ublql+o1JpMJBQUFSEtLg8ViUbscaoZlpo9QWWtfyqV7ejpEUJDKFbl0SaOksCiPudJQyLxFWViZiHRu/Yq12FUTAgAYFmXAgOSBXt2/Zia0O+IjAHswo17CGb3BYrFg06ZNSEq68Idxbm4uMjMzkZBgP5bn5OT4TayG2WzG6tWrkZ+fD6PRCKPRqHZJ1EQbv/gOi8+0RXjfkXg0xorIKy73eFuTyQSTydRofISXGKFMSK+noQYKUspiIcRFAfTKKJcR9lEw3YTQE9Glqs5XYcHaXwARijBZiyeeetDrr6GZkauIiAgUFRWhqKgooBorB+fGyhWj0VjXaOmd2Wyu+15SUlJUroaaylpdjfdX/wgAqAkKxqBHm3Z1TXp6OoqKihAbGwsAlV4v8FLGZmyT6XS1IGAf5bKggaasoqICQgiXX9nZ2c14eSLylYI5H6Nc2K9wvrtfDGK6xDb43OzsbJe/16WlpQAQ2dB2mmmuApmnIzfunmexWC45vcbTbeQLyz74FEdhPzCN7tEeXXp2VbkitywN3B8N96NXJgCFQogUIUSKMiHeiAZGrSIjIyGldPnF5opIO44dLMXKfacAAB1Rjfsev8ft87Ozs13+XsfFxQFARUPbaea0YCBryujN0KFDMWbMGCQkJGDjxo0YNWoUoqOjMXHiRCQnJyMvLw8AkJaWhsLCQpSXl9dta7FYYDKZkJSUhOLiYqSkpDQ6Ypabm4uEhASUlZUhOTn5ouebTCZER184e5KaeuFK9YZeq7CwEHl5eTCbzcjNzfWoBtIOy7ETWPbLMcAQDKOsxpjxj6pdUmOKYG+k6jOikdN7Uspix7+V04IW54B6ItKfpfOXoUbYW59xKUMQGtZ41E1z6L65OmB6G+f27lW7jDrhvXujR/pEn+w7JSUFGRkZWL16NSZNmlR3ai0pKQlZWVlYvXp13XPz8/MRFXXxZfFpaWl185xSUlIwdOhQbNq0qcHXS0tLQ1ZWVl3zU1BQUPfvtLQ0vP3223WjaQUFBSgoKKhrsBp6LUcj6fgeSF/mzVmO8wb7YeOR6/shLCJc5Yrck1JahBCuRqjM7holIURqvaUYHgYwy9v1EVHrObNjB4Z89ylCjN1xMnEgrkm53mevpfvm6tzevajYtk3tMlqN80iR84hPY6cWCwsLL5k8npycfFFD5MxsNqO4uPii18jMzERSUlLd6UbnfaWmpiIxMRGpqalNfi3Sh52bt+O7slpACCQE1+Cme25VuyRP5Qkh0h0rrCsLiOY4HlSuAMyRUqY5bZMlhChUmjMjgDQppfbWmiAij0ibDQfz3oYBwOAzR3DFb17x6evpvrkK791b7RIu0hr1NGdie3Gx/QxHYWHdWopITExssCmr31gB9tGmhISEulOF9ZWVlcFisTT5tUj7bDYb5hSshRShEFJiwsO3eHVNGF+SUuYKIdKFEClQJrfXi7JJAJAihDA6jWZlAnhYaaxiADg3XkSkMye/+BJnd+8CAHR+6EGEdeni09fTfXPlq1Nw/ighIeGi+V1NvVKvKU1dS1+LtOXwmm9w9sw5IDwUw2OC0WfwALVLahJ3uYDKZPUoF/cRkR84XWbBvz7ZhGFh7dG1XajX8gPd0cefntRijonl9TV0RWFDz2/sMaPR2OTXIm2rPV+F8vlz8bj5a9xRtkOT+YFERA35YM5S7AiPxbzEW2BLfdRr+YHusLnyEwkJCTCbL1z85HxKDrgwcuT8nIYaJMf+HFf3OZjNZpjN5rp9OTdLJpMJOTk5Hr9W/UbLMceLtOfokiWoPnECBkjc+WCKVvMDiYguYd62C+uO1QAAugfXYPCdI1rldXV/WtDfmEymuiv4MjIyMGrUqLqr6oqLizFr1ixYLBYkJiZetNhqQkICRo0ahYKCAhiNxrqJ746r9gD7nKnMzExcffXVddu4mweVn5+PzMzMukbI0XA59uW8TAOAi+pp6LWKi4uRk5ODoqIiJCYm1n1vhYWF2LRpU91SEqQNFSWlKHXKD4wdxdO7RKQP9vzAL2ATIT7JD3RHSClb5YXcFiFESVxcXFxJSYnapRCRkxmvmnD6+EmMKN2K66dNRuQVV3ht3/Hx8SgtLS2VUsZ7bacq4PGLSJvWfboG/163BwAwrIPA81kTvLbvxo5fPC1IRC4VffE9tp4Lwv52nbB5yK1ebayIiHypLj8QsOcHjne/Eru3sbkioktYq6sxf/VmAECwrMVTT9yhckXa5gieT05O9puAdSI9Wzx3ab38wI5e2a/JZEJycnKjwfOcc0VEl7g4PzASXXp1U7kibXMEzxOR+o4dLMWKvRZABHmUH9gU6enpSE9Pd5wWbDB4niNXRHSRuvxAAEZZg4fH369uQURETXBk4UL0OWWfAzkuZbDP8gPd0UxzxWF1Im24OD+wL9p4OT/Q02F1IqKmOrNjB6q//hJ3HS7GH9uW4ZqU4arUoZnTghxWJ1Lfzh9/9nl+oKfD6kRETeHIDwQAQ3AIhk4cp1otmhm5IiJ1SSnxfv4aSCHs+YFpN+smP5CI6Kdlq3Bqj33phdbID3RHMyNXRKSu8m/W45af10B0uQIde/dEnyGXqV0SEZFHTpdZ8M/v9iOoz0jccXYfrmqF/EB32FwREWrPV+HQnDlobz2PB07txOW//ZPaJREReeyDOUtxVgQDocFoe93oVskPdIdj/kRkzw88fhwA0PWJJxAaGalyRUREnnHOD+wRVI2Uh0apXBGbK/KyzMxMZGZmql0GNUGJ+QCWfrEZVmFgfiAR6cqF/EDR6vmB7qhfAcFisSAzMxNCCAwdOhQFBQVN3ofJZEJBQQHS0tLqgpbVkJGR0az6vUEr74HezJ2/El937I+5fUYgdvwEiKAgtUsiIvLI+hVrsasmBAAwzGjAgOSBKldkxzlXGmA0GpGTk4OCggLk5OQgJaVpIwdmsxmrV69Gfn4+jEYjjEajbwr1QEJCgiqvq6X3QE+KvvgeW87Zm6mexjbodOVglSsiIvJMXX6gCFXyAx9Qu6Q6HLnSkOY2BGazua6paWpj5i/4HjRd/fzA8eOYH9hcXASZqPX5Kj/QHWYLEpFbyxbUyw/s3V3livSLiyATta7zJ0/iu52lQGg7r+cHuuPpIshsrjSsuLgYEydOREpKCkaNsl/9sHr1aowaNapudKawsBB5eXkwm83Izc1FSkoKkpKSYLFYYDKZkJSUhOLi4rr7CwsLkZmZiTFjxiAhIQEbN26s27er+1NSUhrcl0Nubm7d7ejo6Ea/r9zcXCQkJKCsrAzJyckX7ctkMl20j9TU1Lp/u/ueXL0H1DDLsZP2/EARbM8PfOo+tUsiIvJY6Xvz8Njur1EUk4hrHntIlfxAd/yiuVq7dT/WbTvg9jk9OnXAuJEX5pPsP2rB+19ubXTfk8feeNHtGR+ua/C5Nw7sgZsG9Wx0n55KSkpCVlYW8vLykJWVBaPRiOTkZIwcORKbNm0CcOEU2OrVqzFp0qS6bdPS0urmH6WkpGDo0KHYtGkTUlJSkJGRUfd8x6m0pKQkl/e725fjsaysrLpmpri42O33VP/5BQUFdf9OS0vD22+/XXd6tKCgAAUFBXUNlrvvydV7QA2bP2cZzgun/MC2jPgjIn04s2MHTn71FYIB3N4nFn1Srle7pEv4RXN14nQldhw80aRtzlbVNHkbAG63uaxHbJP31xjH5GxHw2E0GmE2m91uU1hYeMmk7uTk5LpGxXlkyHmEx9X97vblGEFy3oe7ESOz2XzJ8zMzM+tG2hzfn0NqaioSExORmpra6PdEntv548/41sf5gUREvlA/P7D70+NVrsg1v2iuYttHYEB3941Nj04dLrrdNiyk0W1ccbdNbHvf/PXvyak2Z47Ro8LCwrr7EhMTL2pMGrqqr/797vZVXFzcpKsD6zdWgH20KSEhoe5UYX1lZWWwWCwefU/UOCklzIuWoJ21A84EhzM/kIh05bMPP8PGc+1xc3Ab9HvwXlXzA93RTHPluNoGuDBhzFM3DerZ5NNxPTsbLznl54nmbKOGhISEi66aa8kVdA3tyxvrWTWlOfPm9xSoyr9Zj9itP2CCMODMyHtUyQ80mUwwmUyNXm1DROTsdJkFS7YcxlljNxyPjMU/VM4PdEczf7I6rrYpKipqUmNFl3KcrquvOQtruttXUlJSo6coPdlXY48ZjUavfk+BypEfCADh7dri1qfHqFJHeno6ioqKEBsbCwANXm1DRORsgSM/EMB91/RBSHi4yhU1TDPNFbluFDxtHpyf5xjRcW58Gpto3hB3+0pISLikwSosLERZWZnLfTme73xqz2w2w2w2172O8/dhMpmQk5Pj8fdU/71yzPEiu92LFl/IDxw3DsHt2qlcERGRZ8zbdmGtxvID3dHMacFAZrFYMGvWLJjNZmRmZiIrKwupqakoLi6uW2LAMXE7MzOzLi4nJycHxcXFyMnJQVFRERITE+uullu9ejUyMzNx9dVXA7A3No55UrNmzYLFYkFiYmLdKGFD97vbFwDk5+dfshSDxWKpu7Kvvvz8/LrvwbEvx7arV6++aJkGAB7V0dB7UFhYiE2bNiEvL6/lH5LOlZgP4LXtZ9C/61UYHXYGsaO1fWAiInKw2WyYu/AL2ESIpvID3RFSSrVrgBCiJC4uLq6kpETtUoj80sxXTXUxN5PuHISkm65RuSI4FuErlVLGq11LS/D4ReRb6z5dg3+v2wMAGNZB4PmsCSpX1PjxS9utHxG1WNGXF/IDh0TUaqKxIiLyRF1+IKDkB7bOSuwtxeaKyI9Zq6sxf5UjP9DG/EAfYbYgkW+okR/oDrMFieii/MBR3dsxP9BHmC1I5H3VZWWIXfc5ekT1xbl2HVotP9AdZgsSBbj6+YFjxjM/kIj04/B78xB9+jjSTh9H11emai4/0B2eFiTyU875gWOu68P8QCLSjTM7duDkl18CAKKGXYv4a4aqXFHTeDRyJYRIB+BYvChaStnopAJlGwAwAoCUMrc5BRJR05X+shPflVkBYUBCcA1uvneE2iUREXmk1lqLz99dhO4AgoND0P1p9a8ObKpGmyulSTJLKQuV26lCiHR3DZYQYpJzMyWEMAohcqSUmV6pmogaJKXE6fnv4fE9+7EmbhAmPHOv5teEISJyWLFoBfKDuqJzQls8PaQzwuK0mR/ojidH3AxHYwUAUsoCABmNbJPofENKaYEygkVEvmVZ/y0qtm5Fp6rT+P3AKFXyA4mImuN0mQX//ekwAKCyTVv0Sbtf3YKayW1zJYQwAnCVrJukPNaQZCGE54m8ROQVteercHD2bABAcGQk4h97VOWKiIg855wfmJbcGxGR+ozpamzkKhkX5lo5syiPNSQTwCbHvCshRKpyHxH50OLZi1Fc0wY2MD+QiPRFb/mB7jTWXBlhb6Tqc53Mq1BOI44EkCeE2AOgWDk1SEQ+UmI+gGX7z2Jl1yR8ednNzA8kIt24kB8odJMf6I4nlRubulPllGAKgCgAhQD2CCFS3G1TUVEBIYTLr+zs7KaWQBRw5s5fCauw/0qPvvtGiKAglSuyy87Odvl7XVpaCgCRatdHROpbv2ItdtaEAACuMQoMSB6ockUt01hzZWng/mi4H73KlFLmSiktUsoMAKMA5LubpxUZGQkppcsvNldE7m36Srv5gdnZ2S5/r+Pi4gCgQu36iEhd9fMDnxx/r8oVtVxjzVUR7I1UfUYAZlcbKCNU+c73KacJM2EfzSIiL7JWV2Pe58wPVBOzBYmab2vBcpyrtf9bC/mB7nglW1BKaRFCuBqhMruZQ2Vs4H6XzRgRtQzzA9XHbEGi5qkuK0PQ0nw8XWPFtoQk3Pf4OLVLcsvTbEFP5lzlOa227lhUNMfpdoIQom6kSlkHa4yL/YxyXi+LiFrOcrzMnh8IMD+QiHTn8HvzUHv+HMJra/DYuLt0lR/oTqMrtEspc4UQ6crpPqNyn/O4dwKAFCGE0Wk0K1MIkQPgJOzztowAZnmvbCICgPmzlzI/kIh06fTPv9TlBxqHXYv2V12pbkFe5FG2oLuoG2U0KqrefWZwXSsin6rctw/hP/+IsE6XoWsYmB9IRLpRa63Fa/O/QqfOV2B4uVmX+YHueNRcEZG2SClx0PQ2rizbi/4Vpeg1bbqu14QhosCyctEKHJRhOBjbBx0HX4HrdZgf6A6PxkQ65MgPBICet6Wg+8B+KldEROSZ02UWLFHyAyNlDR4ef7+6BfkAR66IdObcmUpseW8B2oL5gZ5SLsRxXPkc7W6qQ71tgAtzTXN9Ux1RYHHOD0wd2ku3+YHusLki0pmFcz7GqugrMVRGYuwDNzM/sBFKk2R2XK0shEgVQqS7a7CEEJOcmykhhFEIkSOl5FxSohbYu13JDxQCPQzVGJU6Wu2SfIKnBYl0pMR8EIWHzsImDDgY2xOdRnNdXg9kOC8DoywXk9HINonON5QroY1er4wogNhsNsz5yCk/8EF95we645/fFZGfmjt/RV1+4BN3DkVwSIjKFWmbErmV4OKhJHdxXACSlYxUIvISf8sPdIfNFZFOaDk/UMOS4ToH1aI81pBMAJsc866EEKng8jJEzVZT5X/5ge6wuSLSAeYHNpsRrgPo3QXPO9bvGwl7QsUeAMVuIr9QUVEBIYTLLwbPEwEnP/kUo/b9gNjzpzWfH+gsOzvb5e91aWkpAEQ2tJ1mmisGnxI1bNmCz/wmP9DT4FMvMjZ1A+WUYArsCyQXAtijpFS4FBkZCSmlyy82VxToqsvKUPrRQvSoPIn06j24f9w9apfksezsbJe/13FxcQBQ0dB2mrlakMGnRK7Z8wOPAiLYL/IDPQ0+9RJLA/dHw/3oVaaU0jHpPUPJT80XQvR2N4JFRJdy5AcCQK+JzyAk1D/yA93RzMgVEblWOP+/zA9sviLYG6n6jADMrjZQRqjyne9TThNmwj6aRUQe2vb9Zvzww3YA/pcf6I5mRq6I6FKV+/ah95rlSIuIwb4+Q5gf2ERSSosQwtUIldnNCJSxgftdNmNE5FqttRZzln6Pwz2vRd8zRzB5wni1S2o1HLki0ihHfqCUNvQ6V4b/98x9frsmjI/lOa227lhUNMfpdoJy2g9A3TpYY1zsZ5TzellE5N7KRStwWNpPAfbq0xNt4uNUrqj1cOSKSKOc8wM73nkHwnv1VLkifZJS5goh0pXTfUblPuerZhIApAghjE6jWZlCiBwAJ2Gft2UEMKu1aibSuwpHfqAIRqSswdgJD6tdUqtic0WkQefOVOKv/92Aq9p1Ql9xjvmBLeQu6kYZjYqqd58ZXNeKqNk+mLPs4vzA9g2uWuCX2FwRadCiOR/DHNIB5p7X4TeDOzI/kIh0w54fWO33+YHucAIHkcaUmg9i9aGzAIAuqMYNaVwwlIj0QUqJOQsv5AeOf+CGgJwrGnjfMZHGzWF+IBHp1DeffY2d1RfyAy+7epDKFamDzRWRhhSv2VCXHzg4nPmBRKQftpoamD//CkG22oDID3SHc66INMJaXY33VhYDCEWwtGHCEzwdSET6cWzZclyxdzO6heyA4YExuskP9AWOXBFpxHI/yg8MNMxGpUDnyA8EgM6dY3HzWP8ctfI0G5UjV0QacOpEGZb9fBQw+Ed+YKBhNioFOvPc+XX5gT3SJ8IQ7J/thafZqBy5ItKA44vyMfjkHhikDQ8zP5CIdGTLt8WYWRqGH2L6IPKaYQGTH+iOf7aWRDpSuW8fTq1cgRulDTd0aYuh9z6tdklERB6ptdbiveUbUB0UivWdL8M9Y5htDmho5IpzFigQ2fMD34GUNghDEAZOfMrv14TxdM4CEWmfc37grXHh6N6vt8oVaYNmRq44Z4EC0U8r1uDUzzvQFoGTH+jpnAUi0raL8gNtNRg7PrDyA93RTHNFFGjOn63EW1//ivN9UzDi1B5cyfxAItKRi/IDk3uhbYfAyg90R3fnH6qPH8ee13NQc/q02qUQtcjCOUtxSoSgKigEHYYPZ34gEenGPkd+IBCw+YHu6Kq5qjl1GjtfnoKtm39B1ox5OHn4qNolETVLqfkgVh88AwDojGrc/ehdKldEROQZKSVmMz/QLV29G0FtI1DZuy+W9LgWh0Q4XvnfJSg1H1S7LKImc84PfJL5gUSkI0Wr1jE/sBG6aq4MwcFI+vMfMDjKfmVCmQjF1LxPsP/nXSpXRuQ55gcSkV7ZrFa0WbwAdx/ciI41ZwI6P9AdXTVXABAcEoLnM5/CjbH2rvm0CMW0uV9g5+btKldG1LgL+YFgfiAR6c6xpctQVVqCAadLMPnGXgGdH+iO7porAAgKMuC3f3oct3W1L5Fz1hCCmR9+gy3fFqtcGZF7zA/0T1ynjwKBc35gm27d0enuwJsr6uk6fbpsrgDAYDBg/B/G4v4+RgDAeUMw/rq0CBu/+E7dwogaYK2owKbNvwIA8wP9jGOdvqKiIqSnp7t9bsW2bag9d66VKiPyngV5+TgkwgD4d36gO+np6SgqKkJsbCwA+G+24CPPPISxAzsDkKgRQfhw2XqUf/e92mURXaLkgwW4z7wedx7ahMeG92V+YACq2P4zdr0yFfOmvInjB4+oXQ6Rx7Z8W4yVp4LxQe+bsDN5JPMDG6H75goA7nv8bky4ugdiqypw//7vYX49ByfXrFG7LKI6lfv24fhnKyAAXNszBjfcM0LtkkgFx1eswA8deuDzoC6Y/K+PYd7Oi3FI+xz5gQAQBIkbxnISe2P8orkCgNEPjca0sTegnQGQtlrs/dvfcfyzFWqXRQSbzYYDTvmB3Sc+AyGE2mWRCno9+0dU978CAHBahODVeV9h8zebVK6KyL36+YE9mB/YKL9prgAg+ppk9J2WjaDwCAASi/JX46O8RWqXRQHu62Vf4u0z7XG0TQd0vOP2gMgPJNcMISH4feYE3NUtHABQJYLwxvJifLn0S5UrI3KtLj8QUPIDOVfUE37VXAFA5MCB6DfjNWzr0h/rOl+Bj/dWYM7/vA+bzaZ2aRSAzp+txEff7cLhiBgs7j0cHceMUbskUpnBYMC43z+KJwZ1hkFK1AoDTN+ZUTDnY7VLI7rExfmBPZkf6CHNNFfevJS5bd++GP67CWgrrQCAz0ur8NYb82CrrfVGqUQes+cH2ofT77m8C8KjjOoWpAGeXsrs7+587G784ea+CJG1AAQKfj2JvL+/j9pa/iFI2rDv59318gNvU7ki/dBMc9WUS5k90ffKy/DK+FFoL2sAAOvKavH31+fAWl3d4n0TeaJ+fuA9jwXemjCueHopcyC47s6b8dKDw9DWZj9OfXW0Cp++8R/YrFaVK6NAJ6XE7I+YH9hcfv1O9RqQgFd/ezeiYT9wbawQyJ01B9Xnq1SujALBXKf8wCduT2J+ILl02bAhePWZ2xAjqzGofD/ivlmJ3dNeQ21lQPedpLKTGzaizZFDgJS4pgPzA5vKr5srAIjr1Q3Tnn0QnYV9xGrLuWBMnzEH587wwEW+U7xmA35yyg8cesswlSsiLevWrzdm/DkN97WrhABwevNm/Jr1EqrLytQujQKQzWrFkdmzMbr0JzxV8gOeHH+P2iXpjt83VwAQG9cJ014Yi+4Ge4O1syYEC2f8C9azZ1WujPxR/fzA8cwPJA8YO8XishnTYbz6agCAZe9+zJg+G7u37lS5Mgo0x5Yuw/nSEgDAVal3Iyauk8oV6U9ANFcA0CHGiOy/jENCSA0utxzE5du+wc7JL6Pm1Cm1SyM/84lTfmBK93aIY36g3/PWBTlBbdog8aUsRI8ajU+6J+PX0Ci89v4aFK8t8mK1RA0rLzmK/QsLANjzAzvedafKFWmLpxfkCCll61XVUBFClMTFxcWVlJT4/LWqzp3Hvr/9HRUb7BmEbbp1R7/p0xAaE+Pz1yb/Z62owIbf/BFfteuOQ8Y4/GPKk4y5aUB8fDxKS0tLpZTxatfSEr44ftlsNix4ayE+OWifvmCQEhOu64WU+1O89hpErvx1+jvYUV6Fm49ux73PPwPj0CS1S9Kkxo5fATNy5RAW3gb9/vIiYm69FQBw5vBh/H362zi8c5+6hZFfKPlgAcJOn8TtJT9i6p2D2FhRsxgMBjz+u7F4coh9LSybEHjn+/1YNPu/apdGfmzLt5ux6YzA2ZA22Jc4hI1VCwRccwUAIjgYvZ57FrF33oVPuw3F5rCOyH5nJcxbf1W7NNKxyn37cHzFSgBA5KDB6HLTcJUrIr27Y+zd+OOt/RFqs6/Rt2RnGd56cz7XwiKvq7XW4r1PvgcABEkbnnpstMoV6VtANlcAIAwG9MiYiO59ewEAKgwhmD7/K/yycYu6hZEu2Ww2LHt7Eaok7PmB6cwPJO+49vYb8VLatWinrIX19bFq5MycjaoqLilD3rNy0QoctjnyA9swP7CFPGquhBDpQohU5cvjFT6FEJOcttXcZAGDwYCJzz2Ou3q2AwBUGkLwev53+HHtRpUrI735etmX+Lg2BrP7jsT5kXciolcvtUsiPzLg6sF4deLtiJX2K563nZH49rW/opZXPJMXOOcHtpM1GDv+fnUL8gONNldKM2WWUhZIKQsAlHnSYAkh8gGYpJQmAMUA8ltcrQ8IITDut2OQ2j8aAFBlCMYbn27G95+vU7ky0gtHfiAAWINDcMXYh1SuiPxR17698NqfHkZPUYU7D21C5E8/4Ne/cC0sarkFcy/kB6YNZX6gN3gycpUhpSx03FAarAx3GwghJgFYLaW0KNuYAQxtQZ0+lzr+AYy7Mg5CSlhFEP73yx1Y89/CxjekgOecH3jfZZ1h7BitckXkr6I6x2D6lAm4dkA3AEDlvr3Y8cKLqNi/X+XKSK/2/bwbXx9lfqC3uW2uhBBGAAkuHkpSHmtIFoBFzncoDZam3fXInZh4fQIM0gabMMD0/V7s/vRztcsiDWN+ILW2kPA2SMz6CzrefjsAoOT0ebzwr09QtOYHlSsjPVry0edO+YHDmR/oJcGNPJ4MwNWYs0V57JKhHaXpMgKIVuZZWQAkSSlzW1Bnqxlx3wiEh4fh319sR0rJT7D8ZzmO1laj871c/p8uZc8PtMfcMD+QWosIDkaP//dbIDoGpo3HcDqoDd5c8RMmlJ9GygOam95KGnVqYxGGb16FiOjesA0YiMuuHqx2SX6jsRbVCHtzVJ+7k/zJAMwAjMo8rUIAxUKIvGZVqILrRg9HziPXYki1/ds8+PbbKF24CFpYcJW0g/mBpCYhBHqOHYMHr7kw2v7Ohv1Y+O4StUsjHbBZrTj49jsIgsSw80cx8bepapfkVzwZ/zM2Y78JUspixw2lwXrY3anEiooKCCFcfmVnZzejhJbpetVg9J85A8HtOwAAti9cgvn/Mx82G9eXIcBaU8P8QA9lZ2e7/L0uLS0FAM6cbaHbH7kLz468DKE2KwDgv7vK8X9/m49aa63KlZGWOecHxj86FiFGo7oF+ZnGmitLA/dHo+HRq7IGtiuDfVTLpcjISEgpXX6p0VwBQERiAvq/PgtVHeOR3+t6fHakBv/KnYvaGqsq9ZB2mD9Zheoq+yRQ5ge6l52d7fL3Oi4uDgAq1K7PG7yVLdhcw0bfgMkPD69bC2vtcWUtrPNcC4sudfzgEeR99QssIRHMD2wir2QLKiNNe6WUUfXulwCiHFcDutimXEop6t1fDmCk84iW02Otli3YHId278f0t1fAIuzzaZLa2vB85pMICQ1VuTJSg7WiAtvSf4OqM2ewvedgjH09C+HtGHPTVMwW9L6S3fsx8+3PcEK5erVXcA2ynn8EHWKM6hZGmvLXGe9gU4VAkLRh2j1DkHjDNWqXpDstyhZUmidXI1RmV42V0zZmIcQlVxm6aqz0oFufnnj19/chFva/CovPGjBr+mxUnTuvcmWkhpIPFsB6pgJBkHhg7B1srEgz4vv0tK+FZbCPWJVXVmNH9jRUnzypcmWkFVu+3YxNFfaxj35tbGysfMSTOVd5zouGKv/OcbqdoCwY6iwTTmthCSFSAbT+WLkXde4eh9eeT0WcsJ8K+rk6BK/NmIPK035xVoM8dGqP+aL8QOP116lcEdHFojrH4NWXnsLVbaqRuv87BO3bjR0vvIhzBw6oXRqprH5+4HjmB/pMo82VYwkFIUSK0iRBWXXdIQFAivNkdWWh0T1K/M0k2Ce4Z3q1chVEdY7FtBfHomeQfQRrtzUUU19/H6dPlqtcGbUGm82GnNmfY2l8Ek6HtmV+IGlWm3YR+NMr6eifcjMAoPrECex4MRP7f9iscmWkppWLVjI/sJV4tFqYlNIkpSxUllYw1XusUEp5yfwrZZtcx5cXa1ZVZLQR2VlPoG+ofVL7QVsocnPmoaacDZa/+3rZlzDXhmJnh3j8PGw08wNJ00RQELr/JgNdnxgHAPgpNBpZi4uwavEqlSsjNdjzAw8BYH5ga+BSrM0Q3i4CL09+CgPDa9GmthrDzfaMr6pjx9UujXzEOT+wjbTi8fH3qVwRUeOEEIhLS0Ps759FYdwQ2IQBszcexEfvLOa6fQHGOT8wlfmBPsfmqpnCwsLwl6zx+H28DR2rKnC+5DB+zczEeQ1cMUTet8gpP/Be5gfqjhAiXQiRqnx5Ejy/SQghla9yp6/VrVGvt/W6bSSeHT0QYTb72lcf77bg//42H1Yrl5UJBM75gd0N1RjN/ECfY3PVAsGhIRj63P9D5/vsoxjVJ07gk6lvYNePP6tcGXlTqfkgVh20X7jQGdW4l/mBuqI0U2ZlWkMBgDIPGqxCAInKV2/layLsF+vo0tUp12PymOGIVNbCWneiBq/PnMOrngPA8UUL0fvMUUBKTGB+YKvgO9xCQgh0e3oC4seOxb62HbEk+jLMWLAO277jxFF/wfxA3ctQUiIA1F1wk9HQk5WLcxZKKc3Kl0WZUxqt1+VkHPoNvQLTMu5ER2lfqmFbpQFTZsyF5aRF3cLIZ05tLILYtAEPHNiA5ztXMz+wlbC58gIhBOIfHQvbyDtgEwacNwQj5+ONKPrye7VLoxYq/voH5gfqmNIoXbLmHoCkhuK4lGbqoiZKCJFe/2IevYpL7IHpLz6CXspaWPutIXj5rx/hxEFOafA3NqsVB995FwAQFNEWSRMeVbmiwMHmyosenpiKR67oCEiJGhGEv3++Fd98ukbtsqiZbFYrPlixEYA9P/CpccwP1KFkuF4I2QI3cVzOlCbM7L2S1NchNhqvvjweg9rY51xFnT6OQ6+8jHP79qtcGXlT8cKlOKvMA44f+wjzA1sRmysvu3/cvZhwdXcIKVErDPi/tbtRWPC52mVRM5xctRq37VqP3hVHkdK9HeITmB+oQ0Y0nHXqqXTn04quaC143hNhEeH4y+SncX9H4K5Dm2A9eQI7Mv+Ciq3b1C6NvOD4wSP43y0nMa/PrTje6zJ0vJtzRZujucHzbK58YHTqbfjtjX0QJG2wCYF3Nh7E8vnL1C6LmsBaUYHD77+P6OozGHvOjMeefkDtkqj5jM3dUBm1imnseVoMnvdEUEgwxvxpAno9+QQAoLbyLH6ZMhXrl+ryokhy8t68T1AtgnAyLBJR99wLQ3Cw2iXpUnOD5zXTXKmdKu9tN919C567bRBCZC0gBD7YfhyfzylQuyzyUMkHC2CtsP/edH96AkLC26hckf/wNFXeSywN3B8Nz0av0gHs8Vo1GiSEQJeHHkTvP/8JMjgYKzoPwj+/O4APTflcC0untn67GUVKfuBlYVZcO3q4yhUFHs00VxERESgqKkJRURHS0xtdhkYXrh5xLSbddzXCZC3iKssQ+fEHOPz+BzxgadzuH3/Bqu9/gQ2C+YE+kJ6ejqKiIsTGxgJApY9frgj2Rqo+IzybRzUKTTuFqFsxt9yC6BezYG4fBwBYaj6Nf70xj2th6UyttRZzlfxAg7Rh/KOjVK4oMGmmufJXg66/Ci+PvQGPWHYg1FaL0oULcfDtdyBtNrVLIxdsNhtmF3yNwi6D8H7izYh7ZgLzA3VMWULBVXNkrh/Z1YAU+NlkdncSr78aL4+9sW4trPUnrZg1Yw7Ocy0s3XDODxwR1wY9+ru6WJZ8jc1VK+h75eUY8vprCO3YEQBwdPlyLH8jD9aaGpUro/rWLvsSZqt9Hau+cVFon8ADkx/Ic140VPl3jtPtBCFEvpvtLT6sTXP6XnU5pv32LnRS1sLafs6AKdPnwnIiIAbwdI35gdrB5qqVtImPx4Cc19Emviu+6TQAC8pCkTtzDqqrqtUujRTnz1biQ6f8wHETmB/oDxzB8UKIFCFEqnKf88TOBAApDax7ZUaAnBZ0Fte7O6a/+Ah6B9kbrAO1IZj8xiKUmA+qXBm5c1F+YFIP5geqiM1VKwrt2BG9pk/Hvo69AQBbzgVhxozZOHfG19NOyBMX5QcO6MT8QD8ipTRJKQuVCBxTvccKpZRRrk4TSikTPTx96Hfax0Yj++UJGBxun3N1EiGYmvcpynf79fx+3ao4cBBFh08DUPID025XuaLAxuaqlbXrGINXX3gE3Qz2Eatfq0MwbdZ7OFN+SuXKAtul+YF3q1wRkfrCwtsg86WncXOM/X8Vw0p/xr7Jk3H6py0qV0b1HZkzF0/t/hJXn9iF8fczP1BtfPdVYIyNQnbm40gIts+52lsbiqm5C2A5dlLlygLXRfmBtyUhOJT5gUSAfS2s37zwFJ4baMSV5ftQW1mJXdnZKFu7Tu3SSHFqYxFOFW1EmM2KhwZ3xeXXMD9QbWyuVNKuQySmvPQkBoTZh9wPy1BMeXMRjh0qVbmywHNJfuCtzA8kciaEwLWPP4Tef/4zRHAwpNWKzf/4NxaZ8mHjlc+qujg/MAJdnxinckUEsLlSVZuIcLz00ngMaWs/OB1DKKb+82OUHzikcmWBQ9bWYv/HnyDCWoVgacOTzA8kalDMLTejb/ZUVLdrj8U9r8USx1pYvPJZNYtm/xfLZQwqg0IRP3Ys8wM1gs2VykLDQvFi1ngMM9o/isSTB3DglVdw7iCvymkNJz5fhe67f8SEXYV4um87dGV+IDWDvyVMuNN+yBDEZb0Ma0gYAODbslrMnDkH5yvPqVxZ4Dl+8Ag+21OOn6J7Y1m/m5kf2Ao8TZgQWlgtXAhREhcXF1eipHcHIpvNhs/fmo/YlYshAAS3b4++r76Ktn0S1S7Nb1nPnMG29AxYKyoQ2rETBr71bxjCwtQuK2DEx8ejtLS0VEoZr3YtLRGox68jew9hVt4yHIX9d6Z7UA0mP/8wjLG8yra1vDHjnbqYm+dGDGDMTStq7PjFkSuNMBgMuON3T6Kbcr7cevo0tr48BdvWb1K5Mv+1c55TfuCE8WysiJqgS+9ueO3FsUgIsl/5fLA2BJPfWIhDew6oXFlgYH6gtrG50pi4tDT0yMhAjQhCQcfByFm2CRtWrVe7LL+z+8dfMGO/wBddBiF40JUwDr9e7ZKIdKd9TBSmvjweQ8JrAQAnEYqpphX4ZdM2lSvzb8wP1D42VxrU6e67YB2XjsMRMagRQfifL37G10u/ULssv+HID6wVBvwY3RsRDz/C/ECiZgoLb4NJk5/GrbH2K27PimD874KvUL75R3UL82MX5Qd2YX6gFrG50qgb0+7AM9f2gkHaYBMG/OfbvVj50Wdql+UX1i6/kB94XXQQ+l15ucoVEelbUHAQMl54Cg8ldkBYbQ3uObAB5mnTcHLN12qX5ncqyk9dnB844X51CyKXNNNcBdLVNp4a+UAKfj9iAIKkDVIIzN1cgv/O+a/aZena+bOV+PDbC/mBT0y4V+WKAo+nV9uQ/qRNTMX02y9DZ2slpNWKvX/7G44sXgItXDjlLxbMWcr8QB3QTHMVERGBoqIiFBUVIT09vfENAsT1t92IP915JUJlLSAEFv5ahgVvfaR2Wbp1aX5gjMoVBZ709HQUFRUhNjYWABis6We6j7wZ/V59FUERbQEAaxavxD//+h6s1VwLq6XOHz6MrhvXoEtlOfMDNU4zzRU1bOjNVyPzoWvRRtpXc1+2/yzm/GM+/xpsoiP7DjE/kKgVRA4ehAE5r+N4515Y3i0Z35bVYsbMOQypb6GD78xGlzMn8OjedXjh0VuYH6hh/GR04oprBuPlx29BO1kDg7QhYsNaHHjrP5CMnvDYnPc+q8sPHMf8QCKfCu/VE0NfyUSswX4l4S/ngzBl5nsoZ4ZqszjyAwGg06gUdB54mcoVkTtsrnSkz6D+mPr0bXjozG4knjmK4ytWYN/f/wFptapdmubt+n7TRfmBycwPJPK5Tj274rVJTmth2UIx+c18HNy9X+XK9KW6qhqF8xZDgvmBesHmSme69+uNe17LRHj3HgCAk2vW4IdZf0PVufMqV6ZdsrYWNR/Mw1jzOnQ9V878QKJW1D7aiOxXJuDKCPsIVhlCkP32SvxStFXlyvSjYO5S5Lftg4W9hqNd6iPMD9QBNlc6FBodjX6zZiIisQ/KQtvhP6ci8dqMOaisOKN2aZp0YtVqVO7bi67nyvDnYV2ZH0jUykLbhOHFl57GiI72q9zOimDMzP8e33GB5EYdP3gEK83lAIDzbduj5713qlwReYLNlU6FdGiP/jNeww/9r8O54DDstoYie9Z8nC6zqF2apljPnMHh+fMBAKEdO6HLQw+qXBH5Iy4l07ig4CCk//lJpPY1AlKiRhiwbNlanPjiS7VL07T35n2CamWu6OMjBiE0LFTligKbp0vJsLnSsaC2bfGnrAnoE2K/xPmALRRTcj/EydJjKlemHQvy8vEL2kGC+YHkO1xKxnOpTz+EjGE9EXeuHHcf3Ih9//gHSvMLePWzC5fmB96gckXk6VIybK50LiKyLV55eTyuULK9jiAUU/+xGEf3HVK5MvXt/vEXrDhWi6U9rsGGwbcyP5BII259cBSyn74NEW3sf+wcnjcP+9/K41pYTpgfqG9srvxAWFgYsrKeQlKk/S+/EyIUU95ajgO/mlWuTD2O/EApBISUGJ12G/MDiTSkw+BBGJD7OkJjYiEBLNx8ENNnzEHlmbNql6YJK/Mv5AfeyvxA3WFz5SeCQ0Px58zxGB5tPzd/SoRi2uxV2PPTDpUrUwfzA4m0L7xnTwz421+xs89QbI5JwI6qIEyZOQ9lR0+oXZqqKspPYcmPF/IDH2V+oO6wufIjQcFB+N0LTyAlrg0AQFitOPTXXJz9dafKlbWu82cr8dG3uwEwP5BI60JjYnDPK88iUVkL65AtFJP/XoCDO/eqXJl6ij/6GNXKFDTmB+oTmys/YzAY8Myzj+GhxPZI27cebU+dwM6XX0HFlsBZU2bRnKWwCPuoFfMDibSvfbQRU1+ZgKS29rmj5QhB9rursP2HLSpX1vrOl5Sg3aqleHpXIW6VZcwP1Ck2V34qbWIarpr4FACB2vPnsCv7Vexf953aZfkc8wOJ9Cm0TRheeOlpjOxk/8PorAjGrMUbsH7lOpUra10H334X0mpFO2s1Hn/mfuYH6pRmPjWuE+N9HW+/DQkv/BnCEITtER3x0ifb8MWSVWqX5VNz5jE/UOs8XSeGAo8hKAgT//QE0vpFAVLCKgz411e/YvmCT9QurVWUO+UHxo5KQdu+fVWuiJpLM80V14nxjeibb0LnFydhVfxVqBUGvLPhAD79YLnaZfnE6a1bEbn3VwTbapkfqGGerhNDgeuhCQ/iN9f1QpC0QQqB/au/QumifL9eC6u6qhrZBRvwbcf+sLVtx/xAndNMc0W+0+2G6/CH0QMRrByo5m89hvx3CtQuy6tkbS0Omd7BtSd2YeL+tXj6SUZEEOnZLfen4MW7r8SwU/sx7MQuHJ4/Hwfe+o/fBtUXzF2Ko6INvu00AIdSHmR+oM6xuQoQ16RcjxfvTUKYtE8YXbz7FN775wLYbDaVK/MOR34gAPS97y507tVN5YqIqKWuvPFq/OaliQjr2BEAcHzFCmyZmYvKCv9aC8s5PzAW1bh7HK9w1js2VwFkyPCheOmRGxAh7X/5rTh8DqY35+u+wTp9shw7FiwCAIR27IjODzI/kFoX54z6Tnj37hjwxl8R3rMXqg1BmH0yFK/Mmo+Tpf6zFhbzA/XD0zmjQgvnsIUQJXFxcXElJSVqlxIQ9v68GzPfK0SFslzBsPbAH158AsEh+pz8/a/cudhwogrDju/EI+MfQOxNzN/Sg/j4eJSWlpZKKePVrqUlePxqHbVnz+KdmSZ8VdseABAlq/HShNvRvX9vlStrma3fbsaMZcUAgMtCazB12m9Urog80djxiyNXAaj35X2QnX4XoqR90b7jh49g79/+DpsO5zLs/vEXrD9pRY0hGIe69kX0DcwPJPJHQW3bYvyU3yGprX2kvVyEYursVdi24SeVK2u+S/IDHxutckXkLWyuAlTXxO6Y9ocHMLS2HPcc3IhT67/BnhkzYauqUrs0j9XPD5yQehPXhCHyY6FhoXjhpQlI6Ww/bVYpgvH6kh/wzYq1KlfWPMwP9F/8P1EA69itC56fkgGjspbKqaIi/Jr9KipPnVa5Ms8wP5Ao8BiCgvDM8+Mwpn80hLIW1r/X7MSyD/S1FlblqdMX5weOv0/lisib2FwFuODISPR77VVEDhoMG4DFZcGY8vr7sBw7qXZpblWdrcRH3+4CALSRtcwPJAowD4x/AL8d3rtuiZkFW4/ig39/qJu1sE4uXoJRBzfBWH3Wnh9obK92SeRFHjVXQoh0IUSq8tXkFT6FEHlNL41aS1BEBPpmT4F56K3YFtUTh2QYpr65CMcPHVG7tAYtmrsUFmEfTr93QEfmBxIFoJvuHYlJ9yQh3GZFkLSh7drPsf9f/9b8WljnS0pwbPky9Kk4gt8bDjM/0A812lwpzZRZSlkgpSwAUNaUBksIkQogpQU1UiswhIbiwazfYXCEfbLoUYRi6j//i5I9+1Wu7FJH9h3C5wec8wPvUrkiIlLL4BuGYuqTI/Dg6Z3oVlmGE6tWYc/MWag9f17t0hp06J13lQZQoHfGM5wr6oc8+UQzpJSFjhtKg5Xhyc6FEMZm1kUqCA0LxaSsp3BNBwEAKBOhyDatwL6fd6lc2cU+e39ZvfxArglDFMh6XdEX98zIQkRv+7IMlo0b8eXLM3Ci9JjKlV2qaNU32LZ9DwAgNmUk8wP9lNvmSmmOXF2+kORh45QCoLDRZ5FmBIeE4LnMp3BzR/tE8dMiBNPmfoGdxdtVrsyuYvt2DNqwEvcd+AHDI2qYH0hEAIDQ6Gj0nzUT7a+8EsfD2mO+6IpX/rEE+381q11aneqqasz+Yis+6n0jCnsMZX6gH2ts5CoZQJmL+y3KYw0SQiSBjZUuGQwGZDz/OG7vZl98ttIQghkL12PL+k2q1iVra3HwPyYIAP2rypD+G67ETkQXBLVtiz5Tp2DfVTehOigE5SIE2bNXY+t3m9UuDYA9P7BMmSva+8orEBIVpXJF5CuNNVdG2Bup+lw1XPVFSyldbUs6YDAY8NTvx+KBPkYAQJUIwobZH+L0jz+qVpNzfmCXhx5EWKdOqtVCRNpkCA5GRtYzGN3F3sScE8HI+bgI33z6tap1XZQfKKvxwBO8wtmfeTLnytjUnQohUp3naXmioqICQgiXX9nZ2U0tgbxkzDMP4bHBnXH98V9x1bGd2PXqNFi+39DqdZw6WY43P9+Kw+FRzA/UmezsbJe/16WlpQAQqXZ93sBsQW0xGAyY8Nw4PDIg5sJaWGt3Yen85arVdHF+4EDmB+qUV7IFhRApAPKklIn17i8HMFJKWeximwQAkFKaldtGAJvq76PeNszm0riydd9g75tvQlqtEIYgdPvjH9B55IhWe/1/5c7FN2W1AIBJN/VG0p2t99rkG8wWpNawdtmXMK3fA6uwjyXc1jUcT/7ukVa9Qm/rd5sxYynzA/1JS7MFiwBEu7jfCKChWYIpAFKFEJOEEJMAZAGIVm5zSQadir7xBvSZPBmGkFBUCgOmL/8RH8/9uFVe25EfCAC9gmtw5e23tMrrEpH+3XTvCEy6174WFgB8fvgc5v313VZbC6vWWov3ljM/MNC4ba6UOVOu5leZG5pPJaU0SSlzHV8A8gCUKbc5wV3HOiQPRd9p2fik57U43qYDPtpxEh/+Z5FPX5P5gUTUUoOHD0X2UyMQJathrD6LhO8+x+4ZM1tlLazP81fiEPMDA44n/5fKc140VPl3jtPtBCFEvi+KI+2JHDgQ48aOQhtpP0W3dF8F3v3HfNhsNp+83trlXzE/kIharOflfTH9jw9inO0QImqrcaqoCDtfmowai8Vnr1lbWYl932yAkJL5gQGm0eZKGX2CECJFWW0dUkrnGZsJAFJcrXulPD8PQIIQIkdZnoF0buC1V+LlR29EO1kDAFh9pBr/98Z7qLXWevV17PmBOwEo+YHj7/Hq/okosMR07YxrZ2aj/VX2/xWd3bUL7095E/t2+GYtrNJF+bjmwI94cs9XmDCsN/MDA4hH51eUU32FSgSOqd5jhVLKKFenCZXnj5JSCillpqsJ8KRPfYZchikTRqOD0mB9U2bD31+fA2t1tdde45L8wE6xXts3EQWmoIgI9JnyMmJGjMBWYw98HtETr85Zja3fenctrPMlJTi6dCkAoFevrrju/lFe3T9pGyevULP16J+AV//fPYiBvaEqOiPw+sw5qD5f1eJ9Mz+QvKm54fPKhTiObXlBjp8wBAej13PPQlx9HQBlLaylRVj7yRqvvcaet2fX5Qd2z5gIwbmiAYWfNrVIl55dMe3ZVHQR9gbr10pgw4w3WjxR9OTCjzCwbD+ElBg3+irmB1KzNTd8XplLalJG64sBcG6pHxFC4Kk/PoZHL4+tWwvrrXW78fG8ZS3e9zefrcUbFdHYEtUT0SOZHxiI2FxRi8XEdcS0F8ciwVCF+w9sQPiPG7BrylRYz55t1v4qtm/HuW/WIuXIVrzQwYLkEdd6uWIKME0On1eWkVntmO6grNs31JdFkjrufeI+/O7GRATLWkgh8NHPx/Hu/3zQ7It0qquqseDr7TgXHIav4wahw8NjvFwx6QGbK/KK9tFGTHtlPAb16w4AOPPLL9j50mRUN/FKHFlbi4N59ml9hpBQDHzmCW+XSgGkBeHzWQAuWmfEsTAy+Z8b7r4Vf7kvGRE25SKd0vN4c9ZsVFc1fQ6pc37gnX2iERPf2au1kj6wuSKvCQ4PR59XXkbUdfZ5DLtLy5E1630cPeD5ytVFS1agbP9BAECX1IeYH0gt1eTweaXpMsK++HGqcqX0JF8VSNow8PokZI8fhSipzCGtEFg845+oPXfO430cP8z8QLJjc0VeZQgJQcKkF2G4OQWLe16LwyIcU/+9FAd27m1021Mny/F/RYcxu28K9nQfwPxA8gYjmh4+nwx7AoVRmadVCKBYCJHX0AbMRvUPPS5LxPRnH0RXUYX+pw4jcfPX+DVrMmrKyz3a/r25nzI/0M80NxuVzRV5nQgOxpDnfoehsW0AABYRimnvfo49W3a43W7+u0txTgSjMjgMMSNTENQmrDXKJf9nbMY2Cc5LxygN1sMNnUqMjIyElNLlF5srfYmJ74zXssbhkU4SBgCVe3Zjx4uTcP7wYbfbbf1uM4rsFzjjstAaXHvbjb4vlnwuOzvb5e91XFwcAFQ0tB2bK/KJoOAg/P7FJzGii71BOiNCMP39r/HzDz+5fH79/MARD/Cqd/IKSwP3R6Ph0auyBrYrQwOnEsm/RLSPRP9XJiN25EgAwNljxzErdz5+Wr/J5fOZH0j1aaa5qqysRHJyMpKTk2EymRrfgDTPYDAg/bnHcU9P+8jpOUMwchZvQPHXP1z0PJvNhjnMDwwYJpMJycnJOHHiBABE+PjlmhM+b4br0S53DRn5GUNwMHo++0fEjXkEK+Ovwq/hHZG7rBhrln95yXO//u8q5gfSRTTzf7CIiAgUFRWhqKgI6eker/FHOvDYbx9GWn/7/9+qRBDe/OwnfLdyXd3jaz/5CnuYHxgw0tPTUVRUhNjYWACo9OVrNTN83gLALIS45P+QTJkILEIIdH38UQy94SoIKVErDPjPejOWzF1a95zaykpELf8Io0p+RCdrJfMDCYCGmivybw+NfwBPXhVXt1jfhys34sTadfb8wPWO/EAr8wPJF5oTPp8Jp7WwlJxUDqkHqHsevxe/v7kPQmQtAIFFO07g7X+8D5vNhtJF+agtL8eQ8v2YcscVzA8kAECw2gVQ4LhjzJ0ID/8S+Wu24sH932H/G2tRPORmWIT9YHTPgE7MDySvk1LmKhE2KVBO9zUUPu+0aGiBECLaeQkGKWVmK5ZNGjP8zltgjDbib0u+R6UhBF8cqULZNBNu/Wk1DADa9R+AmFtvUblK0go2V9Sqbrl3BIZ0i8Lemd/AarXh51NWIMKeH3gf8wPJR+oHztd7rBBAVFO2ocB0xbVXIju6A15/dyXKRCg2nw+BtdMVuK3kJ3RPf4b5gVSHPwnU6qKSrkK/16YhJKItHtn7DUaWbsGTtyUxP5CINK9Hv9547bmH0E3YA+q3G7uj5MY70LZfP5UrIy1hc0WqaHfZZej/+kxEDRmMu+8cjqRbh6ldEhGRR2LiOmHaS09iWHugf1gtbnoqVe2SSGN4WpBUE9G7N/pNf03tMoiImiwisi2ef+lptcsgjdLlyBVXPPYv/Dz9Bz9Lz/B98h/8LP2HNz9LIaX02s6aXYQQJXFxcXElJZ4F/AohoIW6yTv4efqPpnyW8fHxKC0tLZVSxvu4LJ8SQpR06NAhrk+fPgDs63g1tlYff+b9Bz9L/+HJZ2kymWAymbBlyxbU1NScklIaXe5LCz8UbK4CGz9P/xGozVVTjl/KNvyZ9xP8LP2HN49fujwtSERERKRVbK5c8GW2IffduvT6nuh136Q+vf7s6HXfvqTX90Sv+/YmNlcu6PUHQ6/79iW9vid63TepT68/O3rdty/p9T3R6769SStzrmoMBkNw586dPXp+aWkp4uLifFbPiRMnHKGy3Hcr7NuXn6de3xO97rspn+XRo0dhs9msUsoQnxTTSpp6/AL4M+9P++Zn6T/79ubxSyvNVSWAEADHPdwkEkCF7ypCBIBK7rvV9u3Lz1Ov74le992Uz7IjgBopZYSPamkVzTh+AfyZ96d987P0n3177filieaKiIiIyF9wzhURERGRF7G5IiIiIvIiNldEREREXsTmioiIiMiL2FwREREReVGw2gU0hRAiHUCZcjNaSqmP1cTIJSFECoA0KWWG2rVQyyi/m0YAiQAsUspMdSvSJh7D/AePX/7DF8cv3SzFoHzzZilloXI7FTw46ZIQIgnAGOVmgpQyTc16qGWEEOnOv4dCiBzwc70Ej2H+gccv/+Kr45eemqtNUsqhjd1H+qH8z2UMD076JYQwAkiRUhbUu68cQJSU0qJOZdrDY5h/4fFL/3x5/NLFnCvlm01w8VCS8hgRqSMBQL7zHU4HJFe/swGJxzAiTfLZ8UsXzRWAZFyYp+DMojxGRCqQUhYDqD8aY1T+aW71grSLxzAijfHl8UsvzZUR9oNQfa4OVkTUipQDlLMsACaeEryIETyGEWmOr45ferpa0Kh2AUTknhAiAUCSlHKU2rVokFHtAoioYd48full5MrSwP3R4F9+RFqSA4ATfC9laeB+HsOItMNrxy+9NFdFsB+E6jOC8zqINEG5hHkiTwe6xGMYkYZ5+/ili+ZK+WZd/XVn5oGcSH3KGk55jt9HIYRRGWIn8BhGpGW+OH7porlS5ClvAIC6NyNHxXqICHUrVRdJKc3KbSMA55XIyY7HMCKN8dXxSzeLiAIXVjiGfSidKxvrlPIXQQaAFNjXEjEB2Oi8kBvpg/JZ7nHxkEVKGdXa9Wgdj2H6x+OX//Dl8UtXzRURERGR1unptCARERGR5rG5IiIiIvIiNldEREREXsTmioiIiMiL2FwREREReRGbKyIiIiIvYnNFRERE5EVsroiIiIi8iM0VERERkRexuSIiIiLyIjZXRERERF7E5oqIiIjIi9hcEREREXkRmysiIiIiL2JzRURERORFbK6IiIiIvIjNFREREZEXBatdAAUWIUQqgGgAQwFkAnhYeWgogBwppVmt2oiIPCGESFf+WQbAAiABQJGUsli1okhT2FxRa4uWUpqEEBKARUqZCdQ1XfmwN1lERJokhMgHsFFKmavcTgeQJ6UU6lZGWsLTgtRqhBBJABYJIYzKXbOcHrYASHJ6bqoQYlLrVUdE5J7SSCU5GiuFEUCx03N47CKOXFHrcQyZK6NUxVJKi9PDSbA3WI7nFrRqcUREjcvBxX8UAsAoAIWOGzx2EcCRK1LHKAALXdxX6OK5RESqU0bcjbj0OJUCYHVr10PaxuaK1JACpwOUctBKgX2Cu2NYPUed0oiIXEoALozAA4AQIkW5r1C5zWMXAWBzRa1MaaQSYL9i0OFtABlSSrMQIlUZVk9Roz4iIleUpsrimDOq/DcHynwrHrvIGedcUWtzjFoZlcmhibBfaVMI2OcrCCESAHBJBiLSmjQAWUKIPU738dhFlxBSSrVroAAihMiD0xIMDTwnB8BqR8NFRKRFypIyo5xOC/LYRQB4WpBanyeTP1OllIWO+QxERFqjLC2Deo0Uj10EgM0VtSJlyDwBQFEjTy1Qlmto7HlERGpJgdP6VgoeuwgATwtSK1H+ksuBfT2rXACz6q1zRUSkecpE9iwA6bDH3+TVW1SUiM0VERERkTfxtCARERGRF7G5IiIiIvIiNldEREREXsTmioiIiMiL2FwREREReRGbKyIiIiIv+v9nWHyaOzWM3AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 2, figsize=(10, 4))\n",
    "\n",
    "ax[0].plot(np.split(c0, 2)[0][::-1], label=\"True coef.\")\n",
    "ax[0].plot(np.split(opt.x, 2)[0][::-1], ls='--', label=\"Inferred coef.\")\n",
    "ax[0].set_xlabel(\"$i$\")\n",
    "ax[0].set_xlabel(\"$p_i$\")\n",
    "ax[0].legend()\n",
    "\n",
    "ax[1].plot(np.split(c0, 2)[1][::-1])\n",
    "ax[1].plot(np.split(opt.x, 2)[1][::-1], ls='--')\n",
    "ax[1].set_xlabel(\"$i$\")\n",
    "ax[1].set_xlabel(\"$q_i$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef2c1660-b0c1-4d01-87a8-48c15dc3c878",
   "metadata": {},
   "source": [
    "## jax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "307d4429-3d1d-4031-9061-b54142351709",
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax\n",
    "import jax.numpy as jnp\n",
    "from jax import vmap, jit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9c707d6-cdfd-4ae8-8eb7-560a054f857b",
   "metadata": {},
   "source": [
    "$$\n",
    "\\operatorname{Li}_{2}(x)=\n",
    "\\begin{cases}\n",
    "\\frac{\\pi^{2}}{3}-\\frac{1}{2} \\log (x)^{2}-\\sum_{k=1}^{\\infty} \\frac{1}{k^{2} x^{k}}-i \\pi \\log (x), x \\geq 1\\\\\n",
    "\\sum_{k=1}^{\\infty} \\frac{x^{k}}{k^{2}},|x| < 1\\\\\n",
    "-\\frac{\\pi^{2}}{6}-\\frac{1}{2} \\log (-x)^{2}-\\sum_{k=1}^{\\infty} \\frac{1}{k^{2} x^{k}}, x \\leq -1\n",
    "\\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "341f0bc6-9a93-40b6-951d-a3a48aae097f",
   "metadata": {},
   "source": [
    "Approximation to $\\mathrm{Li}_2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "905f9e69-826b-485c-895a-dbd68a27d215",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Li2(x, k_max=3):\n",
    "    \"\"\" Target function, f(x) = P(x) Li2(Q(x)) where P and Q and polynomials with coefficients c\n",
    "    \"\"\"\n",
    "    k_ary = jnp.expand_dims(jnp.arange(1, k_max), 1)\n",
    "    x = jnp.expand_dims(x, 0)\n",
    "        \n",
    "    result = jnp.where(x >= 1.,\n",
    "                jnp.pi ** 2 / 3 \\\n",
    "                - 0.5 * jnp.log(x) ** 2 \\\n",
    "                - jnp.sum(1 / (k_ary ** 2 * x ** k_ary), axis=0) \\\n",
    "                - 1.j * np.pi * jnp.log(x), \n",
    "            jnp.where(x <= -1,\n",
    "                - np.pi ** 2 / 6 \\\n",
    "                - 0.5 * jnp.log(-x) ** 2 \\\n",
    "                - jnp.sum(1 / (k_ary ** 2 * x ** k_ary), axis=0),\n",
    "            jnp.sum(x ** k_ary / k_ary ** 2, axis=0)))\n",
    "    \n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "31d33f61-d5ab-4ad6-9f9c-b7b74fe9b064",
   "metadata": {},
   "outputs": [],
   "source": [
    "def function_target(c, x):\n",
    "    \"\"\" Target function, f(x) = P(x) Li2(Q(x)) where P and Q and polynomials with coefficients c\n",
    "    \"\"\"\n",
    "    p, q = jnp.split(c, 2, axis=-1)\n",
    "        \n",
    "    P = jnp.polyval(p, x)\n",
    "    Q = jnp.polyval(q, x)\n",
    "            \n",
    "    return P * Li2(1 - (Q + 0.j))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "768b617c-ee5b-44c1-92f3-b3deee7f8b43",
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 1701\n",
    "key = jax.random.PRNGKey(seed)\n",
    "k1, k2 = jax.random.split(key)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "cf7df201-0497-413a-a1f5-07247fe525c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Check batching behaviour for future\n",
    "\n",
    "# c = jax.random.uniform(key=k1, shape=[5, 6])\n",
    "# x = jnp.repeat(jnp.expand_dims(jnp.linspace(-10., 1., 20), axis=0), repeats=5, axis=0)\n",
    "\n",
    "# jax.vmap(function_target)(c, x)\n",
    "# mse_batched = vmap(mse)\n",
    "# mse_batched(c, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f2b6a4c1-3462-42b2-94d4-7510cd417cf9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DeviceArray([[-3.7451019e+02-1.2855370e+03j,\n",
       "               3.7494454e+00-4.4372296e+02j,\n",
       "               1.4077888e+01+0.0000000e+00j,\n",
       "              -8.8849468e+00-1.6449562e-07j,\n",
       "               1.0438957e+00+4.6510980e-08j,\n",
       "               2.1050968e+00+1.3794192e-07j,\n",
       "               2.2508863e+01+9.3865225e-08j,\n",
       "               1.2021977e+02+7.3201200e-08j,\n",
       "               4.2288300e+02+7.2235608e-08j,\n",
       "               1.1194384e+03+7.4742928e-08j]], dtype=complex64)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "degree = 4  # Degree of polynomials\n",
    "\n",
    "seed = np.random.randint(42)\n",
    "key = jax.random.PRNGKey(seed)\n",
    "k1, k2 = jax.random.split(key)\n",
    "\n",
    "c0 = jax.random.normal(key=k1, shape=[int(2 * degree)])\n",
    "x = jnp.linspace(-5, 5, 10)\n",
    "\n",
    "function_true = function_target(c0, x)\n",
    "function_true"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "4f0640b0-761e-4e8c-aef9-fa10b259624c",
   "metadata": {},
   "outputs": [],
   "source": [
    "@jit\n",
    "def mse(c, x):\n",
    "    \"\"\" Mean squared error loss\n",
    "    \"\"\"\n",
    "    return jnp.mean(jnp.abs(function_target(c, x) - function_true) ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "aa00a0c3-f006-4b60-8c91-5638965187f4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DeviceArray(41530.46, dtype=float32)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chat = jax.random.normal(key=k2, shape=c0.shape)\n",
    "mse(chat, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "b69894be-6d1e-4d28-b685-399173b28088",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DeviceArray([ 2.3718036e+05, -5.4135121e+04,  1.0235893e+04,\n",
       "             -2.3057173e+03, -5.8207676e+03,  3.9627759e+03,\n",
       "             -4.0384882e+02,  9.6801949e+01], dtype=float32)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jax.grad(mse)(chat, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "6f0f9c3b-4838-4e82-81c9-4f681a8ccaf1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DeviceArray(0., dtype=float32)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check that truth gives zero MSE\n",
    "mse(c0, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "ae6556fc-7ba7-4480-a20c-6ec35fd96f8a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate truth coefficients\n",
    "\n",
    "seed = np.random.randint(42)\n",
    "key = jax.random.PRNGKey(seed)\n",
    "\n",
    "c_hat = 0.5 * jax.random.normal(key=key, shape=c0.shape)\n",
    "x = jnp.linspace(-5, 5, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7061c21d-fd59-4421-97a7-b16f8c5cc77c",
   "metadata": {},
   "source": [
    "Optimize with SGD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "d0124c8b-76b6-4e29-9b62-afcdf1008e65",
   "metadata": {},
   "outputs": [],
   "source": [
    "import optax\n",
    "\n",
    "alpha = 1e-8\n",
    "tx = optax.sgd(learning_rate=alpha, momentum=0.9, nesterov=True)\n",
    "opt_state = tx.init(c_hat)\n",
    "loss_grad_fn = jax.value_and_grad(mse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "0473cc29-aca5-4e7c-90a2-d4be2a50aaa9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loss step 0:  302693.4\n",
      "Loss step 100:  9010.9375\n",
      "Loss step 200:  4167.376\n",
      "Loss step 300:  2311.2239\n",
      "Loss step 400:  1790.7037\n",
      "Loss step 500:  1687.3295\n",
      "Loss step 600:  1665.0557\n",
      "Loss step 700:  1656.0605\n",
      "Loss step 800:  1649.4174\n",
      "Loss step 900:  1643.2598\n",
      "Loss step 1000:  1637.2294\n",
      "Loss step 1100:  1631.2513\n",
      "Loss step 1200:  1625.3083\n",
      "Loss step 1300:  1619.393\n",
      "Loss step 1400:  1613.5049\n",
      "Loss step 1500:  1607.6433\n",
      "Loss step 1600:  1601.8065\n",
      "Loss step 1700:  1595.9901\n",
      "Loss step 1800:  1590.1946\n",
      "Loss step 1900:  1584.4207\n",
      "Loss step 2000:  1578.6671\n",
      "Loss step 2100:  1572.9271\n",
      "Loss step 2200:  1567.2056\n",
      "Loss step 2300:  1561.4967\n",
      "Loss step 2400:  1555.8021\n",
      "Loss step 2500:  1550.1227\n",
      "Loss step 2600:  1544.4536\n",
      "Loss step 2700:  1538.7955\n",
      "Loss step 2800:  1533.1466\n",
      "Loss step 2900:  1527.508\n",
      "Loss step 3000:  1521.8759\n",
      "Loss step 3100:  1516.2515\n",
      "Loss step 3200:  1510.6315\n",
      "Loss step 3300:  1505.018\n",
      "Loss step 3400:  1499.4062\n",
      "Loss step 3500:  1493.8002\n",
      "Loss step 3600:  1488.1957\n",
      "Loss step 3700:  1482.5917\n",
      "Loss step 3800:  1476.9889\n",
      "Loss step 3900:  1471.3885\n",
      "Loss step 4000:  1465.786\n",
      "Loss step 4100:  1460.1816\n",
      "Loss step 4200:  1454.5739\n",
      "Loss step 4300:  1448.9657\n",
      "Loss step 4400:  1443.3505\n",
      "Loss step 4500:  1437.7327\n",
      "Loss step 4600:  1432.1079\n",
      "Loss step 4700:  1426.4794\n",
      "Loss step 4800:  1420.8394\n",
      "Loss step 4900:  1415.1962\n",
      "Loss step 5000:  1409.5428\n",
      "Loss step 5100:  1403.882\n",
      "Loss step 5200:  1393.3906\n",
      "Loss step 5300:  1387.7103\n",
      "Loss step 5400:  1382.0322\n",
      "Loss step 5500:  1376.3505\n",
      "Loss step 5600:  1370.6669\n",
      "Loss step 5700:  1364.9799\n",
      "Loss step 5800:  1359.2906\n",
      "Loss step 5900:  1353.6014\n",
      "Loss step 6000:  1347.9081\n",
      "Loss step 6100:  1342.2098\n",
      "Loss step 6200:  1336.5088\n",
      "Loss step 6300:  1330.8037\n",
      "Loss step 6400:  1325.094\n",
      "Loss step 6500:  1319.3805\n",
      "Loss step 6600:  1313.6638\n",
      "Loss step 6700:  1307.9426\n",
      "Loss step 6800:  1302.2151\n",
      "Loss step 6900:  1296.4834\n",
      "Loss step 7000:  1290.7474\n",
      "Loss step 7100:  1285.0046\n",
      "Loss step 7200:  1279.2568\n",
      "Loss step 7300:  1273.5043\n",
      "Loss step 7400:  1267.7444\n",
      "Loss step 7500:  1261.9803\n",
      "Loss step 7600:  1256.2123\n",
      "Loss step 7700:  1250.4352\n",
      "Loss step 7800:  1244.6545\n",
      "Loss step 7900:  1238.8652\n",
      "Loss step 8000:  1233.0717\n",
      "Loss step 8100:  1227.2704\n",
      "Loss step 8200:  1221.4652\n",
      "Loss step 8300:  1215.6537\n",
      "Loss step 8400:  1209.8339\n",
      "Loss step 8500:  1202.5671\n",
      "Loss step 8600:  1196.7444\n",
      "Loss step 8700:  1190.9086\n",
      "Loss step 8800:  1185.0591\n",
      "Loss step 8900:  1179.2002\n",
      "Loss step 9000:  1173.328\n",
      "Loss step 9100:  1167.4451\n",
      "Loss step 9200:  1161.5521\n",
      "Loss step 9300:  1155.647\n",
      "Loss step 9400:  1149.731\n",
      "Loss step 9500:  1143.8059\n",
      "Loss step 9600:  1137.8705\n",
      "Loss step 9700:  1131.9249\n",
      "Loss step 9800:  1125.9692\n",
      "Loss step 9900:  1120.0046\n",
      "Loss step 10000:  1114.0316\n"
     ]
    }
   ],
   "source": [
    "for i in range(10001):\n",
    "    loss_val, grads = loss_grad_fn(c_hat, x)\n",
    "    updates, opt_state = tx.update(grads, opt_state)\n",
    "    c_hat = optax.apply_updates(c_hat, updates)\n",
    "    if i % 100 == 0:\n",
    "        print('Loss step {}: '.format(i), loss_val)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "737a25c1-6531-456e-bc1d-842c6615c09e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DeviceArray([-2.1205392 ,  1.1840112 ,  0.07616048, -1.1845839 ,\n",
       "             -0.36922392,  0.11557341, -0.25100276,  0.6766643 ],            dtype=float32)"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "20bb3664-13e0-4e86-9520-5ea45fb465a7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DeviceArray([-1.7128636 , -0.66194534,  1.661196  ,  0.39116493,\n",
       "             -0.60357046,  1.650306  , -2.3163958 ,  0.9106825 ],            dtype=float32)"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c_hat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be89aa88-6fab-4dd1-8b8f-aa21b3c3556d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}