PINNs for harmonic oscillators.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0e16d484",
   "metadata": {},
   "source": [
    "## PINNs - Physics-Informed Neural Networks for harmonic oscillators\n",
    "\n",
    "In this notebook we will implement two harmonic oscillation PINNs, for all cases: (1) underdamped, (2) critically damped and (3) overdamped.\n",
    "\n",
    "---\n",
    "\n",
    "The following blog posts/articles were particularly useful:\n",
    "- Lagaris, I. E., et al., \"Artificial neural networks for solving ordinary and partial differential equations\", 1998.\n",
    "- Raissi, M., et al., \"Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations\", 2019.\n",
    "- Moseley B., \"So what is a Physics-Informed Neural Network?\", 2021.\n",
    "- beltoforion.de, \"Damped Harmonic Oscillator\".\n",
    "\n",
    "***\n",
    "\n",
    "In harmonic oscillators, the constants $\\delta$ and $\\omega_0$ are widely used for simplification and represent a combination of other constants that describe the system, namely: $ \\delta = \\dfrac{\\mu}{2m} $ and $ \\omega_0 = \\sqrt{\\dfrac{k}{m}} $, where $\\mu$ is the coefficient of friction, $m$ is the mass of the oscillator and $k$ denotes the spring constant.\n",
    "\n",
    "As the inertial force is opposed by the friction and the spring's restoring force (Hooke's law), the equation of motion for such a dynamic system can be expressed as `inertial force = - friction - restoring force`, or equivalently:\n",
    "\n",
    "$$ m\\ddot{x} = -\\mu\\dot{x} - kx $$\n",
    "\n",
    "$$ or $$\n",
    "\n",
    "$$ \\ddot{x} + \\dfrac{\\mu}{m}\\dot{x} + \\dfrac{k}{m}x = 0 \\ \\ \\ \\ \\ \\ (1)$$\n",
    "\n",
    "The assumption is that the solution of the differential equation that describes such system, is exponential. Therefore, the exponential ansatz of a solution to the differential equation of this system, can be formulated as:\n",
    "\n",
    "$$ x(t) = Ce^{\\lambda t} $$\n",
    "\n",
    "Calculating the first and second derivatives:\n",
    "\n",
    "$$ \\dot{x}(t) = \\lambda C e^{\\lambda t} \\ \\ and \\ \\  \\ddot{x}(t) = \\lambda^2 C e^{\\lambda t} $$\n",
    "\n",
    "Plugging these into (1), we can obtain the characteristic equation:\n",
    "\n",
    "$$ \\lambda^2Ce^{\\lambda t} + \\dfrac{\\mu}{m}\\lambda Ce^{\\lambda t} + \\dfrac{k}{m} Ce^{\\lambda t} = 0 $$\n",
    "\n",
    "which can be simplified to:\n",
    "\n",
    "$$ \\lambda^2 + \\dfrac{\\mu}{m}\\lambda + \\dfrac{k}{m} = 0 $$\n",
    "\n",
    "Solving for $ \\lambda_{1,2} $ and using the $\\delta$ and $\\omega_0$ constants:\n",
    "\n",
    "$$ \\lambda_{1,2} = -\\delta \\pm \\sqrt{\\delta^2 - \\omega_{0}^2} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eab3df23",
   "metadata": {},
   "source": [
    "### Underdamped case\n",
    "\n",
    "The underdamped case describes a harmonic oscillation with $ \\delta < \\omega_0 $.\n",
    "The solution for the differential equation (DE) in the underdamped case is:\n",
    "\n",
    "$ x(t) = e^{-\\delta t} (2A\\cos(\\phi+\\omega t)) \\ \\ $,  with $ \\omega = \\sqrt{\\omega_{0}^2 - \\delta^2} $.\n",
    "\n",
    "Where $A$ describes the amplitude and $\\phi$ the phase angle (see [here](https://beltoforion.de/en/harmonic_oscillator/) for a detailed solution of the DE).\n",
    "\n",
    "For this case, we consider the material/geometry parameters $\\delta = 2$ and $\\omega_0=20$, and try to learn the solution for the range $x \\in [0,1]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "bd23f52c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torch.nn as nn\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from PIL import Image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "69202ae5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def UD_oscillator(d, w0, x):\n",
    "    \"\"\"Defines the analytical solution to the 1D underdamped harmonic oscillator problem. \n",
    "    Equations from: https://beltoforion.de/en/harmonic_oscillator/\"\"\"\n",
    "    assert d < w0\n",
    "    w = np.sqrt(w0**2-d**2)  #omega\n",
    "    phi = np.arctan(-d/w)    #phi - phase angle\n",
    "    A = 1/(2*np.cos(phi))    #A - amplitude\n",
    "    cos = torch.cos(phi+w*x) #cos(phi+omega*t)\n",
    "    sin = torch.sin(phi+w*x) #sin(phi+omega*t)\n",
    "    exp = torch.exp(-d*x)    #exp(-delta*t)\n",
    "    y  = exp*2*A*cos         #y - exact solution\n",
    "    \n",
    "    return y\n",
    "\n",
    "def plot_result(epoch,x,y,x_data,y_data,yh,xp=None):\n",
    "    \"Pretty plot for training results\"\n",
    "    plt.figure(figsize=(8,4))\n",
    "    plt.plot(x,y, color=\"grey\", linewidth=2, alpha=0.8, label=\"Exact solution\")\n",
    "    plt.plot(x,yh, color=\"tab:blue\", linewidth=4, alpha=0.8, label=\"Neural network prediction\")\n",
    "    plt.scatter(x_data, y_data, s=60, color=\"tab:orange\", alpha=0.4, label='Training data')\n",
    "    if xp is not None:\n",
    "        plt.scatter(xp, -0*torch.ones_like(xp), s=60, color=\"tab:green\", alpha=0.4, \n",
    "                    label='Physics loss training locations')\n",
    "    l = plt.legend(loc=(1.01,0.34), frameon=False, fontsize=\"large\")\n",
    "    plt.setp(l.get_texts(), color=\"k\")\n",
    "    plt.xlim(-0.05, 1.05)\n",
    "    plt.ylim(-1.1, 1.1)\n",
    "    plt.text(1.065,0.7,\"Training step: %i\"%(epoch),fontsize=\"xx-large\",color=\"k\")\n",
    "    plt.axis(\"off\")\n",
    "    \n",
    "class FCN(nn.Module):\n",
    "    \"Defines a fully connected network\"\n",
    "    \n",
    "    def __init__(self, n_in, n_out, n_hidden, n_layers):\n",
    "        super().__init__()\n",
    "        activation = nn.Tanh\n",
    "        self.fcs = nn.Sequential(*[\n",
    "                        nn.Linear(n_in, n_hidden),\n",
    "                        activation()])\n",
    "        self.fch = nn.Sequential(*[\n",
    "                        nn.Sequential(*[\n",
    "                            nn.Linear(n_hidden, n_hidden),\n",
    "                            activation()]) for _ in range(n_layers-1)])\n",
    "        self.fce = nn.Linear(n_hidden, n_out)\n",
    "        \n",
    "    def forward(self, x):\n",
    "        x = self.fcs(x)\n",
    "        x = self.fch(x)\n",
    "        x = self.fce(x)\n",
    "        return x\n",
    "\n",
    "def save_gif_PIL(outfile, files, fps=5, loop=0):\n",
    "    \"Helper function for saving GIFs\"\n",
    "    imgs = [Image.open(file) for file in files]\n",
    "    imgs[0].save(fp=outfile, format='GIF', append_images=imgs[1:], save_all=True, duration=int(1000/fps), loop=loop)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51afe369",
   "metadata": {},
   "source": [
    "Let's start by generating the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "cbd89dc7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEWCAYAAABIVsEJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABFiklEQVR4nO3deXxU1d348c83+54A2RMCIYSExRBCBAVFEBfArVpFcalb3arF+uti7fPUWp/a2j51KXWr+rRqSxXqjuIG7soOIewQAiQhLCFkIyFkO78/ZoIhJCHLzNyZyff9es0rM/eee+/3TpL5zj3n3HPEGINSSinVUz5WB6CUUsozaQJRSinVK5pAlFJK9YomEKWUUr2iCUQppVSvaAJRSinVK5pAlHICEblJRL528jE2ichUR5ftCxEZKiJGRPycfSxlPU0gymVE5BoRWSEitSJy0P78RyIiVsfWnoh8LiI/dNK+HfIha4wZbYz53NFlXUVEpopIidVxqN7TBKJcQkR+CvwF+F8gHogD7gQmAwEujsXtvx17QoxKaQJRTicikcDDwI+MMa8bY2qMzTpjzHXGmGP2coEi8mcRKRKRAyLynIgE29dNFZESEfmp/epln4jc3OYY3dn2fhHZD/xDRAaIyHsiUiYiFfbnyfbyjwBnA0+JyBERecq+PFNEPhGRwyKyTURmtzn+IBF5V0SqRWQlkNbFW/Kl/Welff9n2qu8vhGRJ0TkMPCQiKSJyKciUi4ih0RkvohEtTnmbhE5z/78IRFZKCKviEiNvcoqt5dlc0RknX3df0RkgYj8rpPfra/9fT8kIoXARe3W3ywiW+z7KhSRO+zLQ4EPgET7e3BERBJFZIKILBORSvvv+CkRcekXDNV9mkCUK5wJBALvnKLcH4ERQDYwHEgCHmyzPh6ItC+/FXhaRAb0YNuBwBDgdmx/+/+wv04BjgJPARhj/gv4CrjHGBNmjLnH/oH3CfBvIBaYAzwjIqPt+38aqAcSgFvsj85Msf+Msu9/mf31RKDQvv9HAAH+ACQCI4HBwENd7PdS4DUgCni39Xx6Utb+Yf0W8BK29+tV4PIu9nMbcDEwDsgFrmy3/qB9fQRwM/CEiOQYY2qBmUCp/T0IM8aUAs3AfUA0tr+b6cCPuji+spIxRh/6cOoDuB7Y327Zt0Altg/uKdg+LGuBtDZlzgR22Z9PtZf1a7P+IHBGN7dtAIK6iDEbqGjz+nPgh21eXw181W6bvwG/AXyBRiCzzbrfA193cqyhgGl3LjcBRad4H78HrGvzejdwnv35Q8CSNutGAUd7Wtb+u9gLSJv1XwO/6ySmT4E727y+oP25tSv/NnBvm99LySnO+SfAW1b/Deuj44fWsypXKAeiRcTPGNMEYIyZBGBvRPUBYoAQYE2bNnXB9uF8fD+t29vVAWHd3LbMGFN/fKVICPAEMANovYoJFxFfY0xzB+cwBJgoIpVtlvkB/7Qf3w8obrNuT4fvRNfabo+IxALzsFWnhWN7nyq62H5/m+d1QFDb97w7ZbFd7ew19k/vjuJqJ5EuzltEZmJLsiPs8YcAGzrbmYiMAB7HdjUTgu19XdPF8ZWFtApLucIy4BhwWRdlDmG7whhtjImyPyKNMWHd2H93tm0/7PRPgQxgojEmgu+qlaST8sXAF232H2Vs1S53AWVAE7YqplYpXcTb2RDY7Zf/wb4syx7j9W3ic5Z9QFK7nnGDOytsL9/heYtIIPAG8GcgzhgTBSym8/cY4FlgK5BuP+df4fxzVr2kCUQ5nTGmEvgttjaDK0UkTER8RCQbCLWXaQFewFZHHgsgIkkicmE39t+bbcOxJZ1KERmI7VtyWweAYW1evweMEJEbRMTf/jhdREbar1jexNbwHSIio4Abuzh2GdDSbv+dxXjEHmMS8PNTlHeEZdjaIe4RET8RuQyY0EX5hcBcEUm2t0f9ss26AGxtX2VAk/1q5II26w8Ag8TWyaJVOFANHBGRTOCuPp+RchpNIMoljDF/Av4f8AtsbRcHsLUh3I+tPQT78wJguYhUA0uwXSV0R0+3fRIIxnb1shz4sN36vwBX2ntozTPG1GD78LsGKMVWBfRHbB+QAPdgq07bj60B+h+dHdgYU4etkfwbe2+jMzop+lsgB6gC3seWpJzKGNMAXIGtk0Iltque97BdQXbkBeAjYD2wtm2M9vdsLrYkUwFci63BvnX9VmyN9IX29yER+Jm9XI193wscd3bK0eTEqk6llDqRiKwAnjPGdJoUVf+kVyBKqROIyDkiEm+vwroRyOLkKzSltBeWUuokGdiqncKAncCVxph91oak3JFWYSmllOoVrcJSSinVK/2qCis6OtoMHTrU6jCUUsqjrFmz5pAxJqb98n6VQIYOHcrq1autDkMppTyKiHQ4soJWYSmllOoVTSBKKaV6RROIUkqpXulXbSBKKWs0NjZSUlJCfX39qQsrywQFBZGcnIy/v3+3ymsCUUo5XUlJCeHh4QwdOpQTB/pV7sIYQ3l5OSUlJaSmpnZrG0ursETk72KbnnRjJ+tFROaJSIGI5ItITpt1M8Q2rWiBiPyyo+0dIn8hPDEGHoqy/cxf6LRDKeWt6uvrGTRokCYPNyYiDBo0qEdXiVa3gbyEbUKfzswE0u2P27HNFYCI+GKbQnQmttnU5tiH0Has/IWwaC5UFQPG9nPRXE0iSvWCJg/319PfkaUJxBjzJXC4iyKXAa8Ym+VAlIgkYJufoMAYU2gffvo1up6sqHeWPgyNR1naPI5nmi5he0sSNB61LVdKqX7O6iuQU0nixOkyS+zLOlt+EhG5XURWi8jqsrKynh29qgSAL1rG8qemOVzY8Eeebbrk+HKllOfw9fUlOzv7+OPRRx912L7z8vJYvHixQ/YVFnbqSTiffPJJ6urqjr+eNWsWlZWVDjl+T7h7I3pH11Omi+UnLzTmeeB5gNzc3J6NHBmZDFXFPOz/Evf6vcFvGm/ij01zGBTix+we7UgpZbXg4GDy8vKcsu+8vDxWr17NrFmznLL/9p588kmuv/56QkJCAByWvHrK3a9ASjhxvuVkbLPBdbbcsaY/CP7BAAySGv7i/xRn+G7l4borKT5cd4qNlVLurqqqioyMDLZt2wbAnDlzeOGFFwC46667yM3NZfTo0fzmN9/NeLxq1SomTZrE2LFjmTBhAlVVVTz44IMsWLCA7OxsFiw4cRLFTZs2MWHCBLKzs8nKymLHjh0APP7444wZM4YxY8bw5JNPnhTb559/zsUXX3z89T333MNLL73EvHnzKC0tZdq0aUybNg2wDdN06NChTve7e/duRo4cyW233cbo0aO54IILOHr0aJ/fP3e/AnkX29zMrwETgSpjzD4RKQPSRSQV2IttmtFrHX70LPt1xtKHoaoE36hk/ndiCjM+9OP3i7fw7PXjHX5IpbzdbxdtYnNptUP3OSoxgt9cMrrLMkePHiU7O/v46wceeICrr76ap556iptuuol7772XiooKbrvtNgAeeeQRBg4cSHNzM9OnTyc/P5/MzEyuvvpqFixYwOmnn051dTUhISE8/PDDrF69mqeeeuqk4z733HPce++9XHfddTQ0NNDc3MyaNWv4xz/+wYoVKzDGMHHiRM455xzGjRt3ynOdO3cujz/+OJ999hnR0dEnrOtsvwMGDGDHjh28+uqrvPDCC8yePZs33niD66+/vhvvbucsTSAi8iowFYgWkRLgN4A/gDHmOWAxMAvbXNd1wM32dU0icg+2uZh9gb8bYzY5Jcis2d8lEmyXPbce2ca8TwvYur+azPgIpxxWKeVYnVVhnX/++fznP//h7rvvZv369ceXL1y4kOeff56mpib27dvH5s2bERESEhI4/fTTAYiIOPX//5lnnskjjzxCSUkJV1xxBenp6Xz99ddcfvnlhIaGAnDFFVfw1VdfdSuBdKWz/V566aWkpqYeT6Djx49n9+7dfToWWJxAjDFzTrHeAHd3sm4xtgTjcreclcrfv9nNU58W8NS1OafeQCl13KmuFFytpaWFLVu2EBwczOHDh0lOTmbXrl38+c9/ZtWqVQwYMICbbrqJ+vp6jDE97up67bXXMnHiRN5//30uvPBCXnzxRbozkZ+fnx8tLS3HX3fn/oyu9hsYGHj8ua+vr0OqsNy9DcQtRYUEMGfCYD7cuJ+ymmNWh6OU6oMnnniCkSNH8uqrr3LLLbfQ2NhIdXU1oaGhREZGcuDAAT744AMAMjMzKS0tZdWqVQDU1NTQ1NREeHg4NTU1He6/sLCQYcOGMXfuXC699FLy8/OZMmUKb7/9NnV1ddTW1vLWW29x9tlnn7DdkCFD2Lx5M8eOHaOqqoqlS5ceX9fZ8bqzX0dy9zYQt3X16Sm88NUu3lxbwh3npFkdjlLqFNq3gcyYMYNbbrmFF198kZUrVxIeHs6UKVP43e9+x29/+1vGjRvH6NGjGTZsGJMnTwYgICCABQsW8OMf/5ijR48SHBzMkiVLmDZtGo8++ijZ2dnH21ZaLViwgH/961/4+/sTHx/Pgw8+yMCBA7npppuYMGECAD/84Q9Pqr4aPHgws2fPJisri/T09BPW33777cycOZOEhAQ+++yz48tzcnI63K8jqqs60q/mRM/NzTWOnFDqque+pfxIA0t/eo7eZatUF7Zs2cLIkSOtDkN1Q0e/KxFZY4zJbV9Wq7D64MrxyRQeqmWTg3uUKKWUJ9AE0gcXjIrH10dYvGGf1aEopZTLaQLpgwGhAZw5bBAfbNzfrV4VSinlTTSB9NGMMfHsOlTLtgMd98BQSilvpQmkj84fFQfA59t6OFCjUkp5OE0gfRQXEURmfDhfaAJRSvUzmkAc4JyMGFbvOcyRY01Wh6KU6kB5efnxYdzj4+NJSko6/rqhoaHLbVevXs3cuXNPeYxJkyY5KtwTTJ06lVPdftB+eHdX0QTiAOeMiKGx2fBtwSGrQ1FKdWDQoEHk5eWRl5fHnXfeyX333Xf8dUBAAE1NnX/5y83NZd68eac8xrfffuvIkHtEE4gHyx0ykCB/H77dWW51KEp5h/yF8MQYeCjK9tMJ00jfdNNN/L//9/+YNm0a999/PytXrmTSpEmMGzeOSZMmHR/ive2w6g899BC33HILU6dOZdiwYSckltaJoD7//HOmTp3KlVdeSWZmJtddd93xXpqLFy8mMzOTs846i7lz554wXHuro0ePcs0115CVlcXVV199wphVHQ0x39Hw7p0NRe9oOpSJAwT4+ZCTMoCVu7qanVcp1S35C2HRXNv00QBVxbbXcMLI2I6wfft2lixZgq+vL9XV1Xz55Zf4+fmxZMkSfvWrX/HGG2+ctM3WrVv57LPPqKmpISMjg7vuugt/f/8Tyqxbt45NmzaRmJjI5MmT+eabb8jNzeWOO+7gyy+/JDU1lTlzOh5L9tlnnyUkJIT8/Hzy8/PJyfluwNaOhpjvaHj3jsplZWU58J2z0SsQB5mYOogt+6upqmu0OhSlPNvSh79LHq0aj9qWO9hVV12Fr68vYJtc6qqrrmLMmDHcd999bNrU8QwRF110EYGBgURHRxMbG8uBAwdOKjNhwgSSk5Px8fEhOzub3bt3s3XrVoYNG0ZqaipApwnkyy+/PD5PR1ZW1gkf/AsXLiQnJ4dx48axadMmNm/e3OE+uluurzSBOMiE1IEYA6v36FWIUn1SVdKz5X3QOm8GwK9//WumTZvGxo0bWbRoUafDp7cfFr2j9pOOyvTkZuOOxtZrHWJ+6dKl5Ofnc9FFF3UYY3fLOYImEAcZlxJFgK+PVmMp1VeRyT1b7iBVVVUkJSUB8NJLLzl8/5mZmRQWFh4fGbf91LetpkyZwvz58wHYuHEj+fn5AJ0OMQ8nDu/eVTlH0zYQBwny92Xs4EiWawJRqm+mP3hiGwiAf7BtuRP94he/4MYbb+Txxx/n3HPPdfj+g4ODeeaZZ5gxYwbR0dHHh1xv76677uLmm28mKyuL7Ozs4+XGjh3b4RDzcPLw7p2VczQdzt2B/vejrTz3RSH5v7mA0EDNzUq16vFw7vkLbW0eVSW2K4/pDzq8Ad0KR44cISwsDGMMd999N+np6dx3331Wh3WCngznrp9yDjQhdRBPf7aTtUUVnJ0eY3U4SnmurNlekTDae+GFF3j55ZdpaGhg3Lhx3HHHHVaH1CeaQBwoe3AUAHlFlZpAlFInue+++9zuiqMvLG1EF5EZIrJNRApE5JcdrP+5iOTZHxtFpFlEBtrX7RaRDfZ1zquX6oHIYH+GxYSyvqTS6lCUcjv9qbrcU/X0d2RZAhERX+BpYCYwCpgjIqPaljHG/K8xJtsYkw08AHxhjGnbSj3Nvv6kujmrZA+OIq+4Sv9ZlGojKCiI8vJy/b9wY8YYysvLCQoK6vY2VlZhTQAKjDGFACLyGnAZ0NkdL3OAV10UW69lD47izbV7Ka2qJykq2OpwlHILycnJlJSUUFamo1a7s6CgIJKTu99d2soEkgQUt3ldAkzsqKCIhAAzgHvaLDbAxyJigL8ZY57vZNvbgdsBUlJSHBB219q2g2gCUcrG39//+B3YyntY2QZy8q2WtqTQkUuAb9pVX002xuRgqwK7W0SmdLShMeZ5Y0yuMSY3Jsb5DduZZR8TQCPrFzzstEHglFLKHViZQEqAwW1eJwOlnZS9hnbVV8aYUvvPg8Bb2KrErJW/kIDFcxktu8lrSftuEDhNIkopL2RlAlkFpItIqogEYEsS77YvJCKRwDnAO22WhYpIeOtz4AJgo0ui7op9ELixPjvZYFJpMj5OGwROKaWsZlkCMcY0YWvT+AjYAiw0xmwSkTtF5M42RS8HPjbG1LZZFgd8LSLrgZXA+8aYD10Ve6fsg72N9dnJUYLYaRJPWK6UUt7E0hsJjTGLgcXtlj3X7vVLwEvtlhUCY50cXs9FJkNVMWNkNwCbzFAyKHH6IHBKKWUFHY3XkaY/CP7BDJNSgjjGppahLhkETimlrKBDmTiSfewe36UPk3GwmE0+GXDJPK8c00cppfQKxNGyZsN9Gxk94Vw2+2ZgTrvK6oiUUsopNIE4yejECKrrmyipOHrqwkop5YE0gTjJ6MRIADaVVlkciVJKOYcmECfJjA/H10fYXFptdShKKeUUmkCcJMjfl7SYUDZpAlFKeSlNIE40OjFSE4hSymtpAnGikQnh7K+up7KuwepQlFLK4TSBOFFGfAQA2/bXWByJUko5niYQJ8qICwdg2wFNIEop76MJxIniIgKJDPZnq16BKKW8kCYQJxIRMuLD2a4JRCnlhTSBOFlGXDjbDtRgTGeTLSqllGfSBOJkGfHh1NQ3UVpVb3UoSinlUJpAnCwz3t6Qvl/vB1FKeRdNIE6W3toTa/8RiyNRSinH0gTiZJHB/iRGBukViFLK6+iEUi6QER+uXXl7oaSijqc+LeDL7WUcOdZEanQoV+UO5prTB+Pnq999lLKapf+FIjJDRLaJSIGI/LKD9VNFpEpE8uyPB7u7rTsZER9OYVktjc0tVofiMRatL2X6Y1/wdt5exg0ZwOXjkmhqMfz32xu56m/L2Fel86woZTXLrkBExBd4GjgfKAFWici7xpjN7Yp+ZYy5uJfbuoXM+HAamlvYfaj2eJuI6tzra0r4+evryR0ygL9cM47EqODj695dX8qv3tzANc8v59XbzjhhnVLKtay8ApkAFBhjCo0xDcBrwGUu2NblMuJsY2JpNdaprdlTwQNv5jMpbRD/vHXiSQni0rGJ/PPWCRw+0sAPX15NfWOzRZEqpaxMIElAcZvXJfZl7Z0pIutF5AMRGd3DbRGR20VktYisLisrc0TcPTYsJhQfgR0HtSdWV2qPNTH31XUkRgXzzLXjCfL37bDcuJQB/GVONpv3VfPQu5tcHKVSqpWVCUQ6WNb+du21wBBjzFjgr8DbPdjWttCY540xucaY3JiYmN7G2idB/r6kDAxhpyaQLv35422UVh3lsavGEhni32XZczPjuPOcNF5bVcxXO6z5YqBUf2dlAikBBrd5nQyUti1gjKk2xhyxP18M+ItIdHe2dTfDY8Mo0ATSqYKDR3j5291cNzGF3KEDu7XNT85LZ1h0KL96a4NWZSllASsTyCogXURSRSQAuAZ4t20BEYkXEbE/n4At3vLubOtuhseGU3joCE3aE6tDf/5oGyEBftx33ohubxPk78vvvjeG4sNHefnb3c4LTinVIcsSiDGmCbgH+AjYAiw0xmwSkTtF5E57sSuBjSKyHpgHXGNsOtzW9WfRfcNjw2hsNhQdrrM6FLezqbSKDzft57azhzEoLLBH204aHs20jBie+qxAZ35UysUsvQ/EGLPYGDPCGJNmjHnEvuw5Y8xz9udPGWNGG2PGGmPOMMZ829W27iw9NgzQhvSO/N/XuwgJ8OWmyUN7tf0vZmRSU9/ES3oVopRL6e28LpJmTyDaDnKig9X1LFpfyuzcwUQGd91w3pmRCRGcNzKWl77dTV1Dk4MjVEp1RhOIi4QF+pEQGaQJpJ1Xlu2hqcVwcy+vPlrdNXU4lXWNvLqy+NSFlVIOoQnEhbQn1omONTUzf8UezhsZx5BBoX3a1/ghA5iYOpAXvyrUIWOUchFNIC40PDaMnWVHaGnR2QkBPtt6kIq6Rq6bmOKQ/d1xzjD2VdXzyeYDDtmfUqprmkBcKD02nLqGZkp1IEAAXl+zl9jwQM5Od8wNnueMiCUpKpj5K/Y4ZH9Kqa5pAnGh4dqQflz5kWN8vu0gl49Lwteno4EFes7XR7h2YgrfFJRTWKbvsVLOpgnEhTSBfOedvFKaWgxX5CQ7dL9X5Sbj5yPMX1Hk0P0qpU6mCcSFBoYGMCg0QBMI8HbeXsYkRZAR79jh7WPDg7hwdDxvrC2hoUkb05VyJk0gLpamPbEoqagjv6SKi7MSnbL/K3OTqaxr5LNtB52yf6WUjSYQFxseG8aOg0cwpv/2xPpok62X1IzR8U7Z/9nDo4kOC+TNtSVO2b9SykYTiIulx4ZRdbSRQ0f677hNH23cT2Z8OEOj+3bvR2f8fH24LDuRT7ce1PGxlHIiTSAuNvz4mFj9c3bCsppjrNpzmBljnHP10erycUk0NhsW5e9z6nGU6s80gbhYeqyt0bi/Ti71yeYDGIPTE8joxAhGxIXx9rq9Tj2OUv2ZJhAXi4sIJCzQr982pC/ZcoCUgSFkxDm291V7IsLFWYms2VPBPr1xUymn0ATiYiJCWkwoO8tqrQ7F5eobm1m2s5xpGTHY5wlzqlmnJQDwwYb9Tj+WUv2RJhALpMX0z668q3dXcLSxmXMyXDM3/fDYMDLjw1m8QdtBlHIGTSAWSIsNY391PUeO9a+5K77YfpAAXx/OGDbIZcecdVoCq/dUsL+q3mXHVKq/0ARigbQYW0+s/taQ/sX2MiakDiQkwM9lxzxejbVRr0KUcjTX/Ser41q78u4sO8LYwVHWBuMipZVH2X7gCFeNH+zS4w6PDSMjzlaNdfPkVJce29PVNzazvLCcbftrqGtoZmBoAKclR5KdHIWPgwbAVJ5NE4gFhgwKwc9H2NmPRoz9cnsZgMvaP9qadVoCTy7dzoHqeuIiglx+fE9TWdfAs5/v5N8riqjpoJo1ITKIW89K5fozhhDk72tBhMpdWJpARGQG8BfAF3jRGPNou/XXAffbXx4B7jLGrLev2w3UAM1AkzEm11Vx95W/rw9DBoX0q4b0L7aXkRAZRLr96suVZoyJ54kl21m65SDXOmjyKm/12baD/Pw/6ymvbeCSrES+Pz6Z7MFRRAT5cbDmGMsLy3l1ZRG/e38L/15ZxBOzs/vNVbQ6mWVtICLiCzwNzARGAXNEZFS7YruAc4wxWcD/AM+3Wz/NGJPtScmjVVpMWL/pytvU3MLXBYc4Z4Rruu+2NyIujMEDg/lks3bn7cqLXxVyy0uriA4L5P0fn828OeM4Z0QMkcH+iAhxEUFclp3Ea7efyUs3n86xxhau+tsyvVmzH7OyEX0CUGCMKTTGNACvAZe1LWCM+dYYU2F/uRxw7OQRFhoeG8buQ7X9Yv7ujaXV1NQ3cVZ6tCXHFxHOHxnPNzvLqe1nPd+66/kvd/K797cwY3Q8b989mVGJEV2Wn5oRy6Ifn0VOShQ/WZDH377Y6aJIlTuxMoEkAcVtXpfYl3XmVuCDNq8N8LGIrBGR2zvbSERuF5HVIrK6rKysTwE7UlpMGE0thqLDdVaH4nTLdpYDuLT7bnvnjYqloamFr3a4z9+Au1i4qpjfL97KRVkJPHVtTrfbNQaGBvDPWydycVYCf/hgK//4ZpeTI1XuxsoE0lFdRodjnIvINGwJ5P42iycbY3KwVYHdLSJTOtrWGPO8MSbXGJMbE+P6BtzO9KfZCZcXlpMeG0Z0WKBlMZw+dCCRwf58slnnCGlr9e7DPPDWBs5Oj+aJ2dk9nl7Y39eHJ67O5sLRcfx20WbeXV/qpEiVO7IygZQAbft0JgMn/fWJSBbwInCZMaa8dbkxptT+8yDwFrYqMY8xLMY2lLm398RqbG5h9e7Dll59gO2DblpGDJ9uPUBTP6g27I7yI8e459/rSB4QzNPX5RDg17uPA39fH+bNGcfpQwfw8/+sZ0NJlYMjVe7KygSyCkgXkVQRCQCuAd5tW0BEUoA3gRuMMdvbLA8VkfDW58AFwEaXRe4A4UH+xEcEef0VyMa9VdQ2NFueQADOHxVPRV0ja4sqrQ7FcsYYfvaf9Ryua+Dpa3OICPLv0/4C/Xx59vrxRIcFctsrq6mo1XlY+gPLEogxpgm4B/gI2AIsNMZsEpE7ReROe7EHgUHAMyKSJyKr7cvjgK9FZD2wEnjfGPOhi0+hz9JivX9QxeWFhwGYOGygxZHAlBHR+PuK9sYC3lq3l8+2lfHAzEzGJEU6ZJ/RYYH87YbxlNce4/438vv1rJv9haVDmRhjFhtjRhhj0owxj9iXPWeMec7+/IfGmAH2rrrHu+vae26NtT9Gt27raYbHhLHTy6e3XeYG7R+twoP8OTMt2j4nife+56dy6MgxHn5vMzkpUfzgzKEO3feYpEh+cWEmH28+wGurik+9gfJoOhaWhdJiwzhyrImDNcesDsUpWts/zkyzvvqq1fkjY9ldXuf1bU9d+f3iLdQea+KP38/qcaN5d9x6VipnDY/m4UWb+/X73B9oArHQ8Bjv7om1YW8VdW7S/tHq3JFxACzd0j97Y+UVV/Lm2r388OxhpDtpUi8fH+Gx2WMJ9Pfhl2/k09LSf6/2vJ0mEAultRlU0RstL7R1mpuQan37R6ukqGBGJkSwdGv/SyDGGH733maiwwL40dQ0px4rLiKI/5o1klW7K7Qqy4tpArFQbHgg4V48ve3aPZUMiw51i/aPtqZnxrJmTwWVdf2rp9DiDftZvaeCn12QQXgfe111x5Xjkzlz2CD+8MEWDlbrfCzeSBOIhUSEYbFhXnkFYowhr7iS7JQoq0M5ybkjY2luMXyxvf/cld7cYnj8k22kx4ZxVa5rhtQXEX5/xWkca2rht4s2u+SY3qSqrpG1RRUsWl/KW+tKeC+/lHVFFVQdbbQ6tON0OHeLDY8J4+sC7/sgK6k4yqEjxxiXMsDqUE6SnRzFoNAAlm45yGXZXY2e4z3eyy9lZ1ktT1+b45SG886kRofy42nDeeyT7VxXcIhJw60ZD81T7K+qZ8GqYpZsOcCGvR3fkCkCoxIimD4yjitzkkkZFOLiKL9zygQiIvcA89sMaqgcKC02lDfWllBT3+iSagVXWVdcCcA4Nxzq28dHmJYZy8eb9tPU3IKfr3dfiDe3GOYt3UFGXDgzx8S7/Pi3TRnGwjXF/HbRZt6fe5bXv9+9UXy4jsc+3sai/H20GMP4lAH89PwRjEqMIHlACAF+PtQ3NrO34iibSqv5puAQf/10B/OW7uDczFjmTk8n24L/te5cgcQDq0RkLfB34CPTnzvRO1hrT6ydZbWW/AE4y7qiCoL8fciMd05Pn76anhnL62tKWLOngolu1EvMGRatt119PHNdjiUzCQb5+/Jfs0Zx57/W8O+VRQ6/98STNTa38NSnBTz7xU4EuHnSUH5w5tBOrypGJkRw3qg47j0vndLKoyxcXcxL3+7me09/wyVjE/n1RSOJdeGkaaf8KmCM+W8gHfg/4CZgh4j8XkSc242jnzjeE8vLGtLXFVWSlRTltt82z0q33ZXu7b2xWloM8z7dQWZ8ODNGu/7qo9WFo+OYPHwQj328XYc5sdtZdoTLn/mGvyzdwYzR8Xz+86n898Wjul0llRgVzE/OG8HX95/LvdPT+WjTfqY/9gX/XL7HZTfKduu/237Fsd/+aAIGAK+LyJ+cGFu/kDIwBH9focCLGtKPNTWzubSacW7YgN4qPMifM4YNYumWA1aH4lRLtx6ksKyWH00bbuk85iLCgxeP5sixJh7/ZPupN/Byn209yPee+oa9FUd57vrxzJszjoTI4F7tKyzQj/vOH8HHP5lCdkoUv357I7f/c41LEvUpE4iIzBWRNcCfgG+A04wxdwHjge87OT6vZ5veNtSrrkA2l1bT0Nzi1gkE4NzMWHaW1bL7kPeOR/bCl4UkRQUzy4K2j/Yy4sO5fmIK81fsYfuBGqvDscy/lu/hlpdXkTIohPfmns0MB/1uhkaH8sotE/j1xaP4fNtBZs37ivX2tkhn6c4VSDRwhTHmQmPMf4wxjQDGmBbgYqdG108MjwnzqiuQdfbRbt2xB1Zb0zNtd6V/6qXVWOuLK1m5+zA3Tx7qNlWJPzlvBGGBfjzy/harQ7HEM58X8N9vb+TcjFhev3MSSVG9u+rojIhw61mpvPWjyfj6CLP/toz38p03R0t32kAeNMbs6WRd//wrcLA0iigqq6bxN4PgiTGQv9DqkPpkXXEliZFBxLmwMa83UgaFkB4bxtKt3lmN9cJXhYQH+nH16a6576M7BoQGMHd6Ol9sL+Pzbd6ZuDvz+Mfb+NOH27h0bCLP3TCe4IDuzfzYG2OSInnn7smclhTJPf9ex5NLtjulXcQ9vpb0Z/kLGb7zZZrwZY+JhapiWDTXo5PIuqIKt7/6aHXuyFhWFB6mpt59bs5yhOLDdSzesI9rJ6a4XffwH5w5lKGDQnjk/S39ZnKvv32xk3mfFnB17mCeuDobfxdcEQ4KC2T+bRO5IieJJ5fsYFH+PocfQxOI1ZY+TFrLbgAKTKJtWeNRWPqwdTH1wcGaekoqjnpMl+TpmXE0tRi+2nHI6lAc6h/f7MZHhJsmD7U6lJME+PnwwKyR7Dh4hFf7wThZr64s4g8fbOXirAR+f8VpLr2RM9DPl8euGsuz1+Vw8WkJDt+/JhCrVZWQJrY6yp0m6YTlnijvePtHlKVxdFdOShSRwf4s8aLeWEeONbFgVREXZyX0umePs10wKo6JqQN54pPtVHvZ1V9bS7cc4FdvbWBqRgyP92LOeUcQEWaeluCUXniaQKwWmUyoHCOBcna2JJ6w3BPlFVfi5yMOm+XO2fzsc6V/vq2MZi8ZdvydvL3UNjRz46ShVofSKRHh1xePoqKugac/LbA6HKfYur+aua+uY3RiBM/0Yc55d+Z9Z+Rppj8I/sEM99n7XRWWf7BtuQdaV1TJqMQIgvyd10DoaOeOjONwbQN5Tu7y6ArGGP61vIhRCRFuX404JimS7+ck849vdlNUXmd1OA516Mgxbn1pNaGBfrz4g9MJCfDOYQc1gVgtazZcMo+0wBp2mkRMxGC4ZJ5tuYdpbjGsL6l0y/GvunJOegy+PsKnXtAbK6+4ki37qrnujBRErLtxsLt+fmEGvj7Cox96T4fOY03N3PnPNZTXHuPFG3OJj3Tv3oh9YWkCEZEZIrJNRApE5JcdrBcRmWdfny8iOd3d1qNkzSbtwjuoJZj9t67yyOQBsP1ADXUNzR7TA6tVZIg/pw8d4BWzFM5fUURogK/HjDIcFxHEneeksXjDflbtPmx1OA7xyPtbWL2ngseuyiYrOcrqcJzKsgQiIr7A08BMYBQwR0RGtSs2E9s4XOnA7cCzPdjWo6TFhAKw86Dn3hXdegOhu1eddGR6Zhxb99dQUuG5VSlVdY0sWl/K98YlERboOVUmt01JJT4iiP95b7PHT3/7Tt5eXlm2hzumDOOiLMf3enI3Vl6BTAAKjDGFxpgG4DXgsnZlLgNeMTbLgSgRSejmth5leGzr/OieO8TDuqIKBoT4M8TC+Ql669yRsYBtjCJP9cbaEo41tXDtxBSrQ+mRkAA/fjEjg/ySKt5Zv9fqcHptx4EaHnhzAxOGDuRnF2ZYHY5LWJlAkoC2ncBL7Mu6U6Y723qUmLBAwoP82FnmuVcgecWVjEsZ4BF17+0Niw5l6KAQjx2d1xjD/BV7yB4cxehEz+gB19b3spPISo7kTx9u42hDs9Xh9FjtsSbumr+WkABf/nrtOJfcKOgOrDzLjj5l2l+/dlamO9vadiByu4isFpHVZWXuO/OfiDA8Nsxj50evOtrIjoNHPK4BvZWIMH1kHN/uLKeuocnqcHpsxa7D7Cyr5ToPu/po5eMj/PdFo9hXVc8LXxVaHU6PGGN44M0NFJYdYd6ccW4/hI8jWZlASoC2g/QkA+1H/eqsTHe2BcAY87wxJtcYkxsTE9PnoJ0pLcZz50fPL6kE3H8Axa5Mz4yloamFrz3wrvR/rygiIsiPi7MST13YTU1IHcjMMfE8+/lODlTXWx1Ot/1rRRHvri/lpxdkMCmtf03Za2UCWQWki0iqiAQA1wDvtivzLvADe2+sM4AqY8y+bm7rcYbHhnGw5phH3pm7rqgSEcga7HnVJ61yhw4kPNDP40bnPXTkGB9s3Mf3xyc7dYA+V/jlzEyaWwx//mib1aF0y/riSv5n0WamZcRw1zn9b449yxKIMaYJuAf4CNgCLDTGbBKRO0XkTnuxxUAhUAC8APyoq21dfAoOlxbjubMTriuqYHhMGBFuNnBfTwT4+TAlI4alWw96VG+g19eU0NhsPLb6qq0hg0K5afJQXl9bwsa9VVaH06XKugZ+NH8tMeGBPD4729IJu6xiaUuPMWaxMWaEMSbNGPOIfdlzxpjn7M+NMeZu+/rTjDGru9rW033XE8uzEogxxt6AHmV1KH02PTOWsppjbCx17w+vVi0thn+vKGJi6kCGx7rn/PM9dfe04QwICeB372922dSsPdXSYvjpwvUcrKnn6etyGBAaYHVIlugfXQU8xOABwQT4+nhcT6w95XVU1DV6dPtHq6kZsYjgMTcVfl1wiKLDdVx3xhCrQ3GYyGB/7jsvneWFh/l4s3uODvD8V4Us3XqQ/5o10iPve3IUTSBuxM/Xh6HRIR53BbKuuALwnBF4uzIwNICclAEe0w4yf8UeBoUGcOHoOKtDcag5E1IYHhvGHxZvoaHJveYMWVFYzv9+tI2LTktw6wErXUETiJtJiwmj0MN6Yq0rqiQ0wJd0L6lCmT4ylg17q9y+J9D+qnqWbDnIVbmDCfTz7Mbz9vx8ffivWSPZXV7HS9/usjqc48pqjvHjV9eRMjCER79/mkfe8+RImkDczPDYMPYcrnO7b11dWVdUydjBUZbMdeAMnjJX+oJVxTS3GK6d4PmN5x2ZmhHD9MxYHv9kO7sPWV+t29xi+MmCdVQdbeTpa3PcbqZHK2gCcTNpMWE0txj2lFv/D9MdRxua2bKv2qvqgUfEhZEUFezW7SBNzS28tqqIKSNiSPHAoWO6Q0R45PLT8Pfx4f438i3vGffEJ9v5pqCchy8bzajECEtjcReaQNyMp/XE2lhaRVOL8YoG9FYiwnkjY/mm4BD1je45rMZn28rYV1XvFV13uxIfGcR/XzySFbsOM3/FHsvieC+/lKc+K2B2bjKzcwefeoN+QhOImxnWOiqvh7SD5HnwCLxdOXdkHEcbm1m2s9zqUDo0f8Ue4iICmZ4Za3UoTjc7dzBnp0fzhw+2WvJ/sXFvFT/7z3rGDxnA/3xvTL9v92hLE4ibCQnwIykq2GO68q4rrmDwwGBiwgOtDsWhJqYOJCTA1y27kRYfruOL7WVcfXoKfv1g0D4R4U9XZhHo58M9/17n0qvCsppj3P7KagaEBPDc9eO9rrNCX3n/X58HGhYT6jFVWOuKKhk32Huqr1oF+fsyfWQcH23aT1Oze3VomL+iCB8R5kzoP1UpCZHBPDZ7LFv2VfP7xa6ZvbD2WBM/fGU1h+saeOEHuV73JckRNIG4oeGxtkEV3fUu3Fb7qo6yr6reK+7/6MhFpyVwuLaBZYXuU411rKmZhauLOW9kLAmRwVaH41LnZsbxw7NSeWXZHt5d3+HYqQ7T0NTCXfPXsqGkknnXjGNMkueO8eZMmkDcUFpMGHUNzeyrcu/7ELy1/aPV1IwYQgN8eT9/n9WhHPfBhv0crm3ghjOGWh2KJX4xI5PcIQP42X/Ws7aowinHaG4x/Pz19Xy5vYw/XHEaF4yOd8pxvIEmEDfU2hPL3RvS1xVXEuDr47VdGoP8fTl/VBwfbtpPo5tUY/1r+R5So0OZlDbI6lAsEeDnw99uGE98RBC3v7La4feHNDa3cN+CPN7JK+UXMzK4+nTv7uXWV5pA3FDrqLzu3g6SV1TJ6KQIr25YvCgrkcq6Rr4psH6OkC37qlm9p4LrJqb0y5FfWw0KC+TvN51Oc4thzgvLKSp3zDz2x5qauXv+Wt5dX8r9MzL50dThDtmvN9ME4oaiwwKIDPZ36yuQxuYW8vd6ZwN6W1NGRBMe6Md7blCN9a/lewj08+HK8clWh2K54bFhzP/hGRxtbOaa55exbX9Nn/ZXVnOM619cwcebD/DbS0dz19T+N7dHb2gCcUMiQpqb98Tatr+G+sYWr21AbxXo58v5o229say8qbCmvpG31u3lkrGJRIX0z6HD2xuVGMH8H06kqcVwxTPfsHRL77pcrygs59KnvmbD3ir+Omdcvx8gsSc0gbgpW08s970XZF2R94zAeyrfy06ipr7J0qFN3lq3l7qGZm7womHbHWF0YiTv3DOZ1JhQbn15Nb98I5+qo92b0bOitoHfvLORq59fjr+vD6/fOYlLxnrulMBW8LM6ANWxtJgwFq4uoepoI5HB7jdo27qiSqLDAkmK8v6upJOHR5MQGcTra4q5KCvB5cdvaTG89M1uxiZHMtZLe7z1RUJkMK/fOYknlmznhS8L+WDjfm44YwhXnz6YwQNPHids2/4a3lq3l38u201dYzM3TRrKL2ZkEBKgH4c9pe+Ym0qPa21Ir2H8kIEWR3OydfYZCPvDsA6+PsIVOUk8+/lODlbXExsR5NLjf7btIIWHavnLNdkuPa4nCfL35YGZI7kkK5G/frqDpz8v4KnPChg6KIQRceGEBPhyuK6RggM1lFbV4yMw87QE7p2ezog475iGwAqaQNxU6x/11v3ul0AqahvYdaiWq3L7T2Pu93OSefqznby5bi93nuPaBtYXv9pFQmQQs05z/dWPpxmTFMnfbshlT3ktn249yDcF5ewur6W+sYXwID9OTx1I7pABzBiToHeWO4AmEDeVFBVMeKAfW/f1rXeJM+SVVAJ4fQ+stobFhDF+yABeX1PCHVOGuezKa1NpFcsKy3lgZib+/WDcK0cZMiiUmyencvPkVKtD8WqW/EWKyEAR+UREdth/nvRJJCKDReQzEdkiIptE5N426x4Skb0ikmd/zHLtGTifiJARH97n7onOsK6oEh+BrOT+NbzDleOTKTh4hLziSpcd8/++3kVIgC/XeOmkUcqzWfWV5pfAUmNMOrDU/rq9JuCnxpiRwBnA3SIyqs36J4wx2fbHYueH7HqZCeFs2V/tdmNirSuqICM+gtDA/nUBe1FWAsH+vry6ssglxyutPMqi9aVcNT7ZLTtSKGVVArkMeNn+/GXge+0LGGP2GWPW2p/XAFuAJFcF6A4y4iOoqW9yqzGxWloMecWVXjv+VVcigvy5PCeJd/JKqahtcPrxnvtiJwC3u7jNRanusiqBxBlj9oEtUQBdzoojIkOBccCKNovvEZF8Efl7R1Vgbba9XURWi8jqsrIyB4TuOiPjWxvSqy2O5DuFh45QU9/UL+7/6MiNZw7lWFMLr60qdupxDlTX89qqYr6fk9wvukorz+S0BCIiS0RkYwePy3q4nzDgDeAnxpjWT9JngTQgG9gHPNbZ9saY540xucaY3JiYmN6djEVGxH/XE8tdrLOPwJvTTxNIRnw4Zw4bxL+W73HqPCF/+6KQ5haj4zEpt+a0BGKMOc8YM6aDxzvAARFJALD/7PAWXxHxx5Y85htj3myz7wPGmGZjTAvwAjDBWedhpYggf5Kigt2qIX1tUSXhQX4Miw6zOhTL3DhpKHsrj/Lhpv1O2X9ZzTHmr9jD97KTSBl08o1wSrkLq6qw3gVutD+/EXinfQGx9ZP8P2CLMebxduvadoi/HNjopDgtlxkf7lZdedfuqSAnZUC/Hg32/FFxDIsJ5alPC2hpcXwHh6c/K6CxuYW7p2nbh3JvViWQR4HzRWQHcL79NSKSKCKtPaomAzcA53bQXfdPIrJBRPKBacB9Lo7fZTITwtlZdoSGJuvno6iub2T7wRpyUvrP/R8d8fUR7p46nK37a1jSywH8OrPrUC3/Wr6Hq09PYVhM/73KU57Bkn6YxphyYHoHy0uBWfbnXwMdfs01xtzg1ADdSEZ8BE0thp1lRxiZYO3ETXlFlRgDOUOiLI3DHVyWnchflu7gr58WcP6oOIfdWPinD7cS4OfDfeenO2R/SjmT3trq5jLtDenu0A6ytqgCEe+dwrYn/Hx9uGfacDbsreKDjY5pC1leWM4HG/dz5zlpxIa7drwtpXpDE4ibS40OJcDXhy1u0JV3bVElGXHhhAfpTW0A3x+fTGZ8OI+8v6XPc4XUNzbzqzc3MHhgMLedPcxBESrlXJpA3Jy/rw/pcWFsLrU2gbS0GNYVVZAzpH+3f7Tl6yM8eMko9lYe5YUvC/u0r6c+LaDwUC2/v/w0ggO8d4pg5V00gXiAMYmRbCq1dkiTgjLbDYT9vQG9vUlp0cwcE8/TnxdQ2MspiNcWVfDcFzu5IieJs9M9614l1b9pAvEAY5IjOVzbQKmFQ5qs2WObgXC8XoGc5KFLRxPo58t9C9fT2MObC8uPHOPu+WtJiAriNxePdlKESjmHJhAPMCbR1vtq494qy2JYu6eCgaEBDNUb204SFxHEI5ePYX1xJb97b3O3t2toauHe1/Ior23g2evGExmibUvKs2gC8QAjEyLw9RFLE8iaogpy+skMhL1xcVYit52dysvL9nSrPaSxuYX7FuTxdcEhHvneGMYk9a+h8ZV36F/jcXuoIH9f0mPDLEsgFbUNFJbV8v2c/jMDYW/8cuZISivreWTxFqrrG/nJeSPw7eCO/YraBua+to6vdhziv2aN5KrcwRZEq1TfaQLxEKMTI/liexnGGJdfBawrtrV/aAN613x9hL9ck01IgC9//bSAbwoOcd/5Izhz2CD8fH2orm/knXV7+cvSHVTWNfKnK7OYrclDeTBNIB7itKQI3lhbwsGaY8RFuPYmszV7KvD1EcYO1mqWU/Hz9eFPV2ZxxrBB/PHDrdzwfysJ9vclLMiPQ0eO2e7kT4nilVtOY1SitSMLKNVXmkA8RGsd+YaSKuJGuTaBrNpVwZikSEIC9M+lO0SE749P5qKsBD7depDVuyuoPdZE0oBgzk6PJnuwtiUp76CfCB5iVGIEIrCxtIrzRsW57Lj1jc3kFVdy0+ShLjumtwjy92XWaQnMOi3h1IWV8kDaC8tDhAT4kRbj+ob0vOJKGppbmDB0oEuPq5Ryf5pAPMhpSZFs3OvaIU1W7jqMCJyuCUQp1Y4mEA8yOjGC/dX1HKx23R3pK3cdJjM+Qm9yU0qdRBOIBxln70a71j4vubM1NrewZk8FE1P16kMpdTJNIB5kTFIEAb4+rCuqcMnxNuyt4mhjMxM0gSilOqAJxIME+vkyOimCtS5KICt3HQbQBKKU6pAmEA+TkzKA/JIql8yRvnLXYdJiQokOC3T6sZRSnkcTiIfJSRnAsaYWtuxzbm+spuYWVu0+zITUQU49jlLKc1mSQERkoIh8IiI77D87HGRJRHaLyAYRyROR1T3d3hvlDIkCcHo11oa9VdTUNzEpTROIUqpjVl2B/BJYaoxJB5baX3dmmjEm2xiT28vtvUpCZDAJkUFO74n19Y5DiMDk4dFOPY5SynNZlUAuA162P38Z+J6Lt/doOZFHWLthEzwUBU+MgfyFDj/GVwWHGJ0YwcDQAIfvWynlHaxKIHHGmH0A9p+xnZQzwMciskZEbu/F9ojI7SKyWkRWl5WVOSh8C+Uv5PQDC9jbMoDilmioKoZFcx2aRGqPNbGuqIKzhuv83EqpzjktgYjIEhHZ2MHjsh7sZrIxJgeYCdwtIlN6Gocx5nljTK4xJjcmxgs+EJc+zJlmPQDLWkbZljUehaUPO+wQK3cdprHZcJZWXymluuC00XiNMed1tk5EDohIgjFmn4gkAAc72Uep/edBEXkLmAB8CXRre69UVcIIMQyiiuUto5jNF8eXO8pXOw4R6OdD7tB+0zdBKdULVlVhvQvcaH9+I/BO+wIiEioi4a3PgQuAjd3d3mtFJiMCZ/hsYVnLKIz5brmjfLmjjAmpAwny93XYPpVS3seqBPIocL6I7ADOt79GRBJFZLG9TBzwtYisB1YC7xtjPuxq+35h+oPgH8wZPpvZxyD2mDjwD7Ytd4A95bUUHDzCtIxOm5WUUgqwaEIpY0w5ML2D5aXALPvzQmBsT7bvF7JmA3DmR89COSwLPIuhl8w4vryvlmyx1QaeN9J1k1YppTyT3onuibJmk/azT4kJD2TZsB87LHkALN1ygPTYMFIGhThsn0op76QJxEOJCGcPj+bLHWU0t5hTb9AN1fWNrNx1mHNHavWVUurUNIF4sOkj46isa3TY8O5fbCujqcVo9ZVSqls0gXiws0dE4+cjx9st+uqTzQcYEOJPTop231VKnZomEA8WEeTPhNSBfLr1QJ/3dbShmSVbDjBjTDy+PuKA6JRS3k4TiIebPjKO7QeOUHy4rk/7WbLlAHUNzVwyNtFBkSmlvJ0mEA83PdPW4L1kS9+uQt5dX0pseCATdf4PpVQ3aQLxcEOjQxkRF8biDft6vY+qo418sa2MS8YmavWVUqrbNIF4gcuyk1i1u4KSit5VY324cR8NzS1cqtVXSqke0ATiBVo/+N/JK+3V9q+uLGZ4bBhZyZGODEsp5eU0gXiBwQNDyB0ygDfWlmBMz24q3Li3irziSq6bmIKIVl8ppbpPE4iXmDMhhcKyWpbtLO/RdvNX7CHI34crchw3mq9Sqn/QBOIlLspKYGBoAK8s29PtbQ5U1/PGmr1ckZNMZLC/E6NTSnkjTSBeIsjfl6tPH8zHm/ez61Btt7Z58atCmo3hzilpTo5OKeWNNIF4kVvPSiXQz5e/LNl+yrKllUf55/I9XDo2UUfeVUr1iiYQLxIdFsiNk4byzvpSNpVWdVn2Dx9sxRj46QUjXBSdUsrbaALxMneeM4xBoYH8/D/5NDa3dFjm4037WbS+lDvOSSN5gF59KKV6RxOIl4kKCeCRy8eweV81//PeZlu33vyF8MQYeCiKgv+dxk9fW01WciR3T9O2D6VU71kypa1yrgtHx3P7lGE8/2UhjWUF/Gr/LwhvquDb5lHcU347gVLN0zlHCfTztTpUpZQHsySBiMhAYAEwFNgNzDbGVLQrk2Ev02oY8KAx5kkReQi4DSizr/uVMWaxk8P2KA/MzESA578s4B0eI4ojlBLNUNnPS/5/ZPCKAJh0ldVhKqU8mFVXIL8ElhpjHhWRX9pf39+2gDFmG5ANICK+wF7grTZFnjDG/Nk14XoeEeGBWSO5aPm1vNF8NtUmhHE+BVzl+wXB0gBVete5UqpvrEoglwFT7c9fBj6nXQJpZzqw0xjT/bvkFABZAxrJqnr55BWReue5UqpvrGpEjzPG7AOw/4w9RflrgFfbLbtHRPJF5O8i0ukcrCJyu4isFpHVZWVlnRXzXtMfBP/gE5f5B9uWK6VUHzgtgYjIEhHZ2MHjsh7uJwC4FPhPm8XPAmnYqrj2AY91tr0x5nljTK4xJjcmJqbnJ+LpsmbDJfMgcjAgtp+XzLMtV0qpPnBaFZYx5rzO1onIARFJMMbsE5EE4GAXu5oJrDXGHJ9yr+1zEXkBeM8RMXutrNmaMJRSDmdVFda7wI325zcC73RRdg7tqq/sSafV5cBGh0anlFLqlKxKII8C54vIDuB8+2tEJFFEjnfHFZEQ+/o3223/JxHZICL5wDTgPteErZRSqpUlvbCMMeXYela1X14KzGrzug4Y1EG5G5waoFJKqVPSoUyUUkr1iiYQpZRSvSI9nUPbk4lIGdDbmxGjgUMODMcT6Dn3D3rO/UNfznmIMeak+yD6VQLpCxFZbYzJtToOV9Jz7h/0nPsHZ5yzVmEppZTqFU0gSimlekUTSPc9b3UAFtBz7h/0nPsHh5+ztoEopZTqFb0CUUop1SuaQJRSSvWKJpB2RGSGiGwTkQL7bInt14uIzLOvzxeRHCvidKRunPN19nPNF5FvRWSsFXE60qnOuU2500WkWUSudGV8jtad8xWRqSKSJyKbROQLV8foaN34u44UkUUist5+zjdbEacj2edHOigiHQ4w6/DPL2OMPuwPwBfYiW3+9QBgPTCqXZlZwAeAAGcAK6yO2wXnPAkYYH8+sz+cc5tynwKLgSutjtvJv+MoYDOQYn8da3XcLjjnXwF/tD+PAQ4DAVbH3sfzngLkABs7We/Qzy+9AjnRBKDAGFNojGkAXsM2/W5blwGvGJvlQFS74eU9zSnP2RjzrTGmwv5yOeDp8+F25/cM8GPgDbqer8YTdOd8rwXeNMYUARhj+sM5GyBcRAQIw5ZAmlwbpmMZY77Edh6dcejnlyaQEyUBxW1el9iX9bSMJ+np+dyK7RuMJzvlOYtIEra5Zp5zYVzO0p3f8QhggIh8LiJrROQHLovOObpzzk8BI4FSYANwrzGmxTXhWcahn1+WDOfuxqSDZe37OXenjCfp9vmIyDRsCeQsp0bkfN055yeB+40xzbYvqB6tO+frB4zHNs1CMLBMRJYbY7Y7Ozgn6c45XwjkAedimyL7ExH5yhhT7eTYrOTQzy9NICcqAQa3eZ2M7dtJT8t4km6dj4hkAS8CM41tPhdP1p1zzgVesyePaGCWiDQZY952SYSO1d2/60PGmFqgVkS+BMYCnppAunPONwOPGlvjQIGI7AIygZWuCdESDv380iqsE60C0kUkVUQCgGuwTb/b1rvAD+y9Gc4Aqowx+1wdqAOd8pxFJAXbrJA3ePA30rZOec7GmFRjzFBjzFDgdeBHHpo8oHt/1+8AZ4uIn30m0InAFhfH6UjdOeci7BPbiUgckAEUujRK13Po55degbRhjGkSkXuAj7D14vi7MWaTiNxpX/8cth45s4ACoA7btxiP1c1zfhDbzJDP2L+RNxkPHsm0m+fsNbpzvsaYLSLyIZAPtAAvGmM67ArqCbr5O/4f4CUR2YCtaud+Y4xHD/EuIq8CU4FoESkBfgP4g3M+v3QoE6WUUr2iVVhKKaV6RROIUkqpXtEEopRSqlc0gSillOoVTSBKKaV6RROIUkqpXtEEopRSqlc0gShlIft8I/kiEiQiofZ5KcZYHZdS3aE3EiplMRH5HRCEbRDDEmPMHywOSalu0QSilMXsYzWtAuqBScaYZotDUqpbtApLKesNxDahUTi2KxGlPIJegShlMRF5F9uMealAgjHmHotDUqpbdDRepSxkn/mvyRjzbxHxBb4VkXONMZ9aHZtSp6JXIEoppXpF20CUUkr1iiYQpZRSvaIJRCmlVK9oAlFKKdUrmkCUUkr1iiYQpZRSvaIJRCmlVK/8fyH5E7R/AbO7AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# initialize material/geometry parameters\n",
    "d, w0 = 2, 20\n",
    "\n",
    "#analytical solution over the full domain\n",
    "x = torch.linspace(0,1,500).view(-1,1)\n",
    "y = UD_oscillator(d, w0, x).view(-1,1)\n",
    "\n",
    "# slice out a small number of points \n",
    "# that will serve as training data\n",
    "x_data = x[0:100:20]\n",
    "y_data = y[0:100:20]\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(x, y, label=\"Exact solution\")\n",
    "plt.scatter(x_data, y_data, color=\"tab:orange\", label=\"Training data\")\n",
    "plt.title(\"Generated training data\")\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"y\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e5f956b",
   "metadata": {},
   "source": [
    "#### Vanilla FCN\n",
    "First let's try fitting the data using a vanilla fully-connected network:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "d6ef3e84",
   "metadata": {},
   "outputs": [],
   "source": [
    "torch.manual_seed(123)\n",
    "model = FCN(1,1,32,3) # (n_in, n_out, n_hidden, n_layers)\n",
    "optimizer = torch.optim.Adam(model.parameters(),lr=1e-3)\n",
    "epochs = range(1, 1001)\n",
    "losses = []\n",
    "files = []\n",
    "\n",
    "for epoch in epochs:\n",
    "    optimizer.zero_grad()\n",
    "    yh = model(x_data)\n",
    "    loss = torch.mean((yh-y_data)**2)# MSE\n",
    "    loss.backward()\n",
    "    optimizer.step()\n",
    "    losses.append(loss) # keep track of loss in each epoch\n",
    "        \n",
    "    # plot results during training\n",
    "    if (epoch) % 10 == 0: # every 10 epochs\n",
    "        \n",
    "        yh = model(x).detach()\n",
    "        \n",
    "        plot_result(epoch,x,y,x_data,y_data,yh)\n",
    "        \n",
    "        file = \"./PINNplots/\"+str(epoch)+\".png\"\n",
    "        plt.savefig(file, bbox_inches='tight', pad_inches=0.1, dpi=100, facecolor=\"white\")\n",
    "        files.append(file)\n",
    "        \n",
    "        #show plot every 500 epochs:\n",
    "        #if (epoch) % 500 == 0: plt.show() \n",
    "        #else: plt.close(\"all\")\n",
    "        \n",
    "        plt.close(\"all\")    \n",
    "            \n",
    "save_gif_PIL(\"ud_nn.gif\", files, fps=20, loop=0) # save gif"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42b375fe",
   "metadata": {},
   "source": [
    "![SegmentLocal](ud_nn.gif \"gif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "213c5888",
   "metadata": {},
   "source": [
    "Due to the fact that the training sample represents only a limited region of the full domain, the network is unable to capture the behavior of the system outside that limited region."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4c09359",
   "metadata": {},
   "source": [
    "#### PINN\n",
    "\n",
    "With a PINN we can incorporate our prior knowledge on the harmonic behavior of this dynamic system. This can be achieved by adding the underlying DE to loss function of the network. This way, the data loss aims to minimize the error in relation to the training sample available, while the 'physics loss' ensures that the learned solution is consistent with the underlying DE, by penalising the residual of the DE over a set of locations sampled from the full domain. In this case, we evaluate the physics loss over a uniform 30-point grid over the full domain ([0,1])."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8a151c3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# sample 30 uniformly spaced points\n",
    "x_physics = torch.linspace(0,1,30).view(-1,1).requires_grad_(True)\n",
    "\n",
    "mu, k = 2*d, w0**2\n",
    "\n",
    "torch.manual_seed(123)\n",
    "model = FCN(1,1,32,3)\n",
    "epochs = range(1, 20001)\n",
    "losses = []\n",
    "losses_data = []\n",
    "losses_physics = []\n",
    "optimizer = torch.optim.Adam(model.parameters(),lr=1e-4)\n",
    "files = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "41a57554",
   "metadata": {},
   "outputs": [],
   "source": [
    "for epoch in epochs:\n",
    "    optimizer.zero_grad()\n",
    "    \n",
    "    # Data loss\n",
    "    yh = model(x_data)\n",
    "    loss_data = torch.mean((yh-y_data)**2) # MSE\n",
    "    \n",
    "    # Physics loss\n",
    "    yh_physics = model(x_physics)\n",
    "    \n",
    "    # calculate first and second derivatives\n",
    "    # dy/dx\n",
    "    dx  = torch.autograd.grad(yh_physics, x_physics, torch.ones_like(yh_physics), create_graph=True)[0]\n",
    "    \n",
    "    # d^2y/dx^2\n",
    "    dx2 = torch.autograd.grad(dx, x_physics, torch.ones_like(dx),  create_graph=True)[0]\n",
    "    \n",
    "    # computes the residual of the harmonic oscillator differential equation\n",
    "    physics = dx2 + mu*dx + k*yh_physics\n",
    "    loss_physics = (1e-4)*torch.mean(physics**2)\n",
    "    \n",
    "    # backpropagate joint loss\n",
    "    loss = loss_data + loss_physics# add two loss terms together\n",
    "    losses.append(loss)\n",
    "    losses_data.append(loss_data)\n",
    "    losses_physics.append(loss_physics)\n",
    "    loss.backward()\n",
    "    optimizer.step()\n",
    "    \n",
    "    \n",
    "    # plot the result as training progresses\n",
    "    if epoch % 150 == 0: \n",
    "        \n",
    "        yh = model(x).detach()\n",
    "        xp = x_physics.detach()\n",
    "        \n",
    "        plot_result(epoch,x,y,x_data,y_data,yh,xp)\n",
    "        \n",
    "        file = \"./PINNplots/\"+str(epoch)+\".png\"\n",
    "        plt.savefig(file, bbox_inches='tight', pad_inches=0.1, dpi=100, facecolor=\"white\")\n",
    "        files.append(file)\n",
    "        \n",
    "        #if epoch % 6000 == 0: plt.show()\n",
    "        #else: plt.close(\"all\")\n",
    "        plt.close(\"all\")\n",
    "            \n",
    "save_gif_PIL(\"ud_pinn.gif\", files, fps=20, loop=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "226eda0e",
   "metadata": {},
   "source": [
    "![SegmentLocal](ud_pinn.gif \"gif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a138382",
   "metadata": {},
   "source": [
    "We can see that the PINN aims to ensure that our solution fits well over the full domain, and not only the limited region that the available data might represent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "59274792",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA9K0lEQVR4nO3deXxU1d348c93luwLZINAAknYlDWyKgoIbrhUSrVWpeJaxRZ9rE/9SWufirWtS31qa/Up1VasO2KV0roiWlBR2QoIsoUYIOwEyL7OnN8f9yZOQhImMJNJMt/363Wduedu39zB+c45595zxRiDUkqp8OUIdQBKKaVCSxOBUkqFOU0ESikV5jQRKKVUmNNEoJRSYU4TgVJKhTlNBEqFORGZKyIvhjoOFTqaCFS7EJECETk/1HEopY6niUCpEBERV6hjUAo0EagQE5FIEfm9iOy1p9+LSKS9LEVE/iUix0TkiIh8LCIOe9m9IrJHREpFZKuInGeXO0RkjojsEJEiEXlNRJLsZVEi8qJdfkxEVolIjxbiOl1E/m2vt0lELrfLzxSR/SLi9Fl3uohs8OP4WSJiRORmEdkFfNjCsS8TkXX2sVeIyHCfZQUi8lMR+UpEjorIfBGJ8ln+AxHJs8/XYhHp5bNsiIgssZcdEJGf+Rw2QkSet8/nJhEZ7bNds+dadR2aCFSo3QecCeQCI4CxwM/tZf8NFAKpQA/gZ4ARkUHAbGCMMSYeuAgosLe5E/g2MAnoBRwFnrKXXQ8kAplAMjALqGwakIi4gX8C7wNpwB3ASyIyyBjzOVAOTPHZ5FrgZT+OX28ScLodd9NjjwSeBW6zY/wzsLg+Odpm2Nv2AwbWny8RmQI8BFwFpAM7gVftZfHAB8C7dlz9gaU++7zcXrcbsBh40t6utXOtugpjjE46BX3C+vI4v5nyHcAlPvMXAQX2+18C/wD6N9mmP3AQOB9wN1m2GTjPZz4dqAVcwE3ACmD4CWKdAOwHHD5lrwBz7fe/Ap6138djJYa+fhw/CzBATivH/hPwYJOyrcAkn/M4y2fZJcAO+/1fgUd9lsXZx84CrgH+08Ix5wIf+MwPBipPdK516jqT1ghUqPXC+uVab6ddBvBbIA94X0TyRWQOgDEmD7gL6wvsoIi86tME0hd4025WOYb1xezBqlG8ALwHvGo3Qz1q//pvLqbdxhhvk7h62+9fBr5j/0r/DrDWGFP/N7R2/Hq7WzkffYH/rt/e3kemzzlpur3v+Wp0Lo0xZUCRHXcmVtJtyX6f9xVAlIi4TnCuVRehiUCF2l6sL796fewyjDGlxpj/NsbkAN8C7q5vnzbGvGyMOcfe1gCP2NvvBi42xnTzmaKMMXuMMbXGmAeMMYOB8cBlwMwWYsqs74/wiWuPfeyvsL5wL6Zxs1Crx/dZp7Uhf3cDv26yfYwx5hWfdTKbO180OZciEovVvLTH3m+/Vo7bolbOteoiNBGo9uS2O2zrJxdWk8vPRSRVRFKAXwAvQkOnaX8REaAE65e1R0QGicgU+xd5FVY7v8c+xjzg1yLS195HqohMs99PFpFhdkdvCVaziYfjfYHV3PP/RMQtIudiJaJXfdZ5Gas/YCKw0Ke8xeP76RlgloiME0usiFxqt/HX+5GIZNid0D8DFvjEdKOI5Nrn5jfAF8aYAuBfQE8RuUusDvp4ERl3omBOcK5VVxHqtimdwmPCats2TaZfAVHAE8A+e3oCiLK3+bG9XTlWp/H/2OXDgZVAKXAE60uul73MAdyN1a5eitUc8ht72TV2eTlwwD6Wq4V4hwDLgGLgK2B6k+V9AC/wVpPy1o6fZf/dzR7TZx9TgVXAMfucLATifc7jT+2YjgF/A2J8tp1lH7P+vGT4LBuK1UF8FKspaI5dPhd40We9hjhbO9c6dZ1J7A9eKdUJiEgBcIsx5oNQx6K6Dm0aUkqpMKeJQCmlwpw2DSmlVJjTGoFSSoW5TjfoVUpKisnKygp1GEop1amsWbPmsDEmtbllnS4RZGVlsXr16lCHoZRSnYqI7GxpmTYNKaVUmNNEoJRSYU4TgVJKhblO10eglGpftbW1FBYWUlVVFepQlB+ioqLIyMjA7W5uYN3maSJQSrWqsLCQ+Ph4srKysMb/Ux2VMYaioiIKCwvJzs72ezttGlJKtaqqqork5GRNAp2AiJCcnNzm2psmAqXUCWkS6DxO5rMKm0SwdX8pv3t/K4fLqkMdilJKdShhkwjyDpbxxId5HCmvCXUoSqk22r9/P1dffTX9+vVj8ODBXHLJJWzbtu2U93vDDTfw+uuvH1e+evVq7rzzzlPeP8Bzzz3H7Nmz/S4PhbDpLHbYtSWvDrKnVKdijGH69Olcf/31vPqq9ZC4devWceDAAQYOHBiUY44ePZrRo0cHZd8dUdjUCOrbzbzeE6yolOpQPvroI9xuN7NmzWooy83NZcKECRhjuOeeexg6dCjDhg1jwQLrqZ3//ve/mTRpEldddRUDBw5kzpw5vPTSS4wdO5Zhw4axY8eOhn198MEHTJgwgYEDB/Kvf/2rYfvLLrsMgLlz53LTTTdx7rnnkpOTwxNPPNGw7YsvvsjYsWPJzc3ltttuw+OxnuI5f/58Bg4cyKRJk/j0009P+Dfu3LmT8847j+HDh3Peeeexa9cuABYuXMjQoUMZMWIEEydOBGDTpk0Nxxw+fDjbt28/ldMLaI1AKdUGD/xzE1/tLQnoPgf3SuD+bw1pcfnGjRsZNWpUs8veeOMN1q1bx/r16zl8+DBjxoxp+MJcv349mzdvJikpiZycHG655RZWrlzJH/7wB/74xz/y+9//HoCCggKWLVvGjh07mDx5Mnl5eccdZ8uWLXz00UeUlpYyaNAgbr/9dvLy8liwYAGffvopbrebH/7wh7z00ktccMEF3H///axZs4bExEQmT57MGWec0eo5mD17NjNnzuT666/n2Wef5c4772TRokX88pe/5L333qN3794cO3YMgHnz5vFf//VfzJgxg5qamobkcyrCpkbgqK8RaCJQqsv45JNPuOaaa3A6nfTo0YNJkyaxatUqAMaMGUN6ejqRkZH069ePCy+8EIBhw4ZRUFDQsI+rrroKh8PBgAEDyMnJYcuWLccd59JLLyUyMpKUlBTS0tI4cOAAS5cuZc2aNYwZM4bc3FyWLl1Kfn4+X3zxBeeeey6pqalERETwve9974R/x2effca1114LwHXXXccnn3wCwNlnn80NN9zAM8880/CFf9ZZZ/Gb3/yGRx55hJ07dxIdHX1K5xDCqUZgpzyv5gGlTlprv9yDZciQIc126ILVf9CSyMjIhvcOh6Nh3uFwUFdX17Cs6eWWzV1+6bsvp9NJXV0dxhiuv/56HnrooUbrLlq06JQvt63fft68eXzxxRe89dZb5Obmsm7dOq699lrGjRvHW2+9xUUXXcRf/vIXpkyZckrH0xqBUqpDmzJlCtXV1TzzzDMNZatWrWLZsmVMnDiRBQsW4PF4OHToEMuXL2fs2LFt2v/ChQvxer3s2LGD/Px8Bg0a5Nd25513Hq+//joHDx4E4MiRI+zcuZNx48bx73//m6KiImpra1m4cOEJ9zV+/PiGjvCXXnqJc845B4AdO3Ywbtw4fvnLX5KSksLu3bvJz88nJyeHO++8k8svv5wNGza06e9tTlATgYhMFZGtIpInInOaWX6uiBSLyDp7+kWwYqlPBPpoTqU6FxHhzTffZMmSJfTr148hQ4Ywd+5cevXqxfTp0xk+fDgjRoxgypQpPProo/Ts2bNN+x80aBCTJk3i4osvZt68eURFRfm13eDBg/nVr37FhRdeyPDhw7ngggvYt28f6enpzJ07l7POOovzzz+fkSNHnnBfTzzxBPPnz2f48OG88MIL/OEPfwDgnnvuYdiwYQwdOpSJEycyYsQIFixYwNChQ8nNzWXLli3MnDmzTX9vc4L2zGIRcQLbgAuAQmAVcI0x5iufdc4FfmKMuczf/Y4ePdqczINpPtl+mO//9QsWzjqLMVlJbd5eqXC1efNmTj/99FCHodqguc9MRNYYY5q9JjaYNYKxQJ4xJt8YUwO8CkwL4vFa1XDVkHYSKKVUI8FMBL2B3T7zhXZZU2eJyHoReUdEmu2JEpFbRWS1iKw+dOjQSQVT3/ni0aYhpZRqJJiJoLlu86bfwmuBvsaYEcAfgUXN7cgY87QxZrQxZnRqarPPXj4hp6O+j+CkNldKqS4rmImgEMj0mc8A9vquYIwpMcaU2e/fBtwikhKMYPSGMqWUal4wE8EqYICIZItIBHA1sNh3BRHpKXabjYiMteMpCkYwDUNMaB5QSqlGgnZDmTGmTkRmA+8BTuBZY8wmEZllL58HXAncLiJ1QCVwtQnSZUxaI1BKqeYF9c5iu7nn7SZl83zePwk8GcwY6jXcUKZVAqU6HafTybBhw6itrcXlcnH99ddz11134XC03KhRUFDAihUrGoZu8FdcXBxlZWWnGnKnEoZ3Foc4EKVUm0VHR7Nu3To2bdrEkiVLePvtt3nggQda3aagoICXX365nSLs3MInETSMNaSZQKnOLC0tjaeffponn3wSYwwFBQVMmDCBkSNHMnLkSFasWAHAnDlz+Pjjj8nNzeXxxx9vcb2WtDTE9b59+5g4cSK5ubkMHTqUjz/+GI/Hww033NCw7uOPPx708xBI4TPonA4xodSpe2cO7P8ysPvsOQwufrhNm+Tk5OD1ejl48CBpaWksWbKEqKgotm/fzjXXXMPq1at5+OGHeeyxxxqeMVBRUdHsei1paYjrl19+mYsuuoj77rsPj8dDRUUF69atY8+ePWzcuBGgYcjoziLsEoE2DSnVNdT/qKutrWX27NmsW7cOp9PZ4iMs/V2vXktDXI8ZM4abbrqJ2tpavv3tb5Obm0tOTg75+fnccccdXHrppQ1DXncWYZQIrFePZgKlTl4bf7kHS35+Pk6nk7S0NB544AF69OjB+vXr8Xq9LQ4a9/jjj/u1Xr2WWg8mTpzI8uXLeeutt7juuuu45557mDlzJuvXr+e9997jqaee4rXXXuPZZ5895b+zvYRNH4HoMNRKdQmHDh1i1qxZzJ49GxGhuLiY9PR0HA4HL7zwQsMDXOLj4yktLW3YrqX1WtLSENc7d+4kLS2NH/zgB9x8882sXbuWw4cP4/V6ueKKK3jwwQdZu3ZtUM9BoIVNjUCHmFCq86qsrCQ3N7fh8tHrrruOu+++G4Af/vCHXHHFFSxcuJDJkycTGxsLwPDhw3G5XIwYMYIbbrihxfVaMn36dD777DNGjBiBiDQMcf23v/2N3/72t7jdbuLi4nj++efZs2cPN954I177oehNH1bT0QVtGOpgOdlhqHcWlTPpt//mse+O4MpRGUGITKmuSYeh7nw60jDUHUpynPWouYOlVSGORCmlOpawSQRxkS6SYyPYVVQR6lCUUqpDCZtEANAnOYbVO4+y+4gmA6WUqhdWieA7IzPIO1jGhEc/YuKjH/E/izayuuCI3mSmlAprYXPVEMB1A718N+M3HKYbSxnLn9eczguf7yQ7JZZrxmZyzdg+xEe5Qx2mUkq1q7BKBGx8najDG8mIT+f60k+Y6XZRmHkB86ov5Ddvl/PHD/P4/pl9ueWc7IbOZaWU6urCp2nowCb48FeQ0Bvu3gyzPkHGzSLzyAp+XfRj1p/2PN/OqmXesh1MfPQjHl+yjdKq2lBHrZTCGoa6fpC37373u1RUVFBQUMDQoUNPed/z5s3j+eefP6V9nHvuua2OW9TRhU8iOPK19RrVDUSsga4u+rWVFCb/nMQ9H/Pg3tv47OKDTBiQyh+Wbmfiox/xl4/zqapt/Q5EpVRw1Q9DvXHjRiIiIpg3b96JN/LTrFmzmDlzZsD21xmFTyLodQY4XHDm7Y3LI2Jh0j3woy8gfQQ9P7yLecmv8o/bz2Ro70R+9dZmzvvfZXyadzg0cSulGpkwYQJ5eXkAeDwefvCDHzBkyBAuvPBCKisr2bFjByNHjmxYf/v27YwaNQqwhqYePHgww4cP5yc/+QkAc+fO5bHHHgMgLy+P888/nxEjRjBy5Eh27NjR7LDTrXnllVcYNmwYQ4cO5d57722Is7lhqp944omGeK6++urAnqg2CJ8+gsTe8LN94Ipofnm3TLj+n7DkF/DZk4wo2csL1z/LioJSfr5oIzP+8gWzJ/fnvy8c2DBukVLh5pGVj7DlyJaA7vO0pNO4d+y9fq1bV1fHO++8w9SpUwHrS/6VV17hmWee4aqrruLvf/873//+90lMTGTdunXk5uYyf/58brjhBo4cOcKbb77Jli1bEJFmh4qeMWMGc+bMYfr06VRVVeH1evnTn/503LDTLdm7dy/33nsva9asoXv37lx44YUsWrSIzMzMZoepfvjhh/n666+JjIwM6dDV4VMjgJaTQD2H02oumvowbPkXvHY947MSeOvOCXxvdCZPfpTHzxdt1MtNlWpn9WMNjR49mj59+nDzzTcDkJ2dTW5uLgCjRo2ioKAAgFtuuYX58+fj8XhYsGAB1157LQkJCURFRXHLLbfwxhtvEBMT0+gYpaWl7Nmzh+nTpwMQFRVFTEwMY8aMYf78+cydO5cvv/yS+Pj4FuNctWoV5557LqmpqbhcLmbMmMHy5csbDVP97rvvkpCQAFjjIc2YMYMXX3wRlyt0v8vDp0bQFmfeDk43vPXf8PqNRH/3OR6+YhjdYt38eVk+2Smx3DIhJ9RRKtXu/P3lHmj1fQRNRUZ+c3Wf0+mksrISgCuuuIIHHniAKVOmMGrUKJKTkwFYuXIlS5cu5dVXX+XJJ5/kww8/bNi+rcNON6elfXTv3r3ZYarfeustli9fzuLFi3nwwQfZtGlTSBJCeNUI2mLMLXDxo1bN4K27EeDei07jwsE9ePTdreQdLD3hLpRSoREVFcVFF13E7bffzo033ghAWVkZxcXFXHLJJfz+978/LrEkJCSQkZHBokWLAKiurqaioqLZYadbMm7cOJYtW8bhw4fxeDy88sorTJo0qdlhqr1eL7t372by5Mk8+uijHDt2jLKysmCdklZpjaA1426D8kOw/LeQlIPjnB/zm+8M4/zfLePnizby6q1nhTpCpVQLZsyYwRtvvNHwtLDS0lKmTZtGVVUVxphmnyv8wgsvcNttt/GLX/wCt9vNwoUL+fjjj48bdrol6enpPPTQQ0yePBljDJdccgnTpk1j/fr1xw1T7fF4+P73v09xcTHGGH784x/TrVu3oJyLEwmbYahPmjHw91tg4+tw1fMweBrzP/2aB/75FS//YBzj+6W0XyxKhUBnHYb6scceo7i4mAcffDDUobS7tg5DrTWCExGBaU/BsZ3wj9mQPoJrxvZh3rIdPLF0uyYCpTqg6dOns2PHjkZ9AKpl2kfgD3cUXPFXQODvtxDlMNx0djaf5x9hy/6SUEenlGrizTffZMOGDaSk6A81f2gi8Ff3vnDZ76BwFax9jqtGZxLpcvD8ZztDHZlSQdfZmpDD2cl8VpoI2mLoFZA1AT78Nd2dlUzL7cWba/dQXKljEqmuKyoqiqKiIk0GnYAxhqKiIqKiotq0nfYRtIUIXPggPH0urHmOmWfdwGurC3lt1W5+MFHvK1BdU0ZGBoWFhRw6dCjUoSg/REVFkZHRtueyayJoq15nQPZE+HweQ8/8IWOzk3huRQE3np2Fy6kVLNX1uN1usrOzQx2GCqKgfnOJyFQR2SoieSIyp5X1xoiIR0SuDGY8ATNuFpTuhR0fcfM52ew5Vsn7Xx0IdVRKKXVSgpYIRMQJPAVcDAwGrhGRwS2s9wjwXrBiCbj+F1jDWX/5Guef3oM+STE8+8nXoY5KKaVOSjBrBGOBPGNMvjGmBngVmNbMencAfwcOBjGWwHJFwOBpsOUtnJ5qbhifxeqdR9lQeCzUkSmlVJsFMxH0Bnb7zBfaZQ1EpDcwHWj1KRMicquIrBaR1R2mw+q0S6G2Anat4LujM4iNcPLcioJQR6WUUm0WzETQ3KD9Ta8/+z1wrzGm1UeAGWOeNsaMNsaMTk1NDVR8pybrHHBGwvYPiI9yc+WoDP61fh+HSqtDHZlSSrVJMBNBIZDpM58B7G2yzmjgVREpAK4E/k9Evh3EmAInIhb6joe8JQDMHJ9FjcfLKyt3hTgwpZRqm2AmglXAABHJFpEI4Gpgse8KxphsY0yWMSYLeB34oTFmURBjCqx+U+DwNig9QL/UOCYOTOXFz3dS5/GGOjKllPJb0BKBMaYOmI11NdBm4DVjzCYRmSUis4J13HbVd7z1umsFADPG9eFgaTUfb9fnGyulOo+g3lBmjHkbeLtJWbMdw8aYG4IZS1CkjwB3DOz8DIZMZ/KgNJJiI3h9TSGTT0sLdXRKKeUXvRX2VDjdkDG6oUYQ4XIwLbcXS746wLGKmhAHp5RS/tFEcKr6jIf9G6GqGIArR2VQ4/Hyz/VN+8WVUqpj0kRwqvqeBRjYvRKAIb0SGdgjjre+3BfauJRSyk+aCE5VxhhwuGDnpw1FFw3pycqvj3CkXJuHlFIdnyaCUxURC71GQsEnDUUXDemJ18AHm3UgOqVUx6eJIBCyJ8CetVBdCsCQXgn07hbN+5s0ESilOj5NBIGQNQGMB3Z9DoCIMPm0VD7bcZhavblMKdXBaSIIhMxx4HDD18sbis7ul0J5jUdHJFVKdXiaCAIhIsbqNC74uKHozJxkRGBFXlEIA1NKqRPTRBAo2RNg33qoPAZA99gIBqcn8OkOHW5CKdWxaSIIlKwJYLwN/QQA4/sls3bnMarrWh1lWymlQkoTQaBkjLGeT+DTPDSyT3dqPF427ysNYWBKKdU6TQSB4o6ykoHPjWVn9OkOwH92HQ1VVEopdUKaCAIp62yrn8Aed6hnYhTpiVGs230stHEppVQrNBEEUtY5x/UT5GZ24z+7joUuJqWUOgFNBIGUMQacEY36Cc7o041dRyooKtNnGSulOiZNBIHkjob0XNi9qqFoWO9uAGzcWxKamJRS6gQ0EQRaxhjY+x+os0YePT09HoAt+zQRKKU6Jk0EgZY5BjzVsP9LALrFRJCeGMVmTQRKqQ5KE0GgZYy1Xgu/aR46PT2BLfv1XgKlVMekiSDQEntDQm8oXNlQdFrPePIOlukdxkqpDumEiUBEeojIX0XkHXt+sIjcHPzQOrGM0Y06jE9PT6DOa8g7WBbCoJRSqnn+1AieA94Detnz24C7ghRP15AxFop3Qel+wKoRAGzV5iGlVAfkTyJIMca8BngBjDF1gLZxtCbT7iewH2jfNzkWp0P4+nB5CINSSqnm+ZMIykUkGTAAInImUBzUqDq79BHWg2r2rAEgwuUgs3s0+ZoIlFIdkMuPde4GFgP9RORTIBW4MqhRdXauSEg7DfZvaCjKTokl/5AmAqVUx3PCRGCMWSsik4BBgABbjTG1QY+ss+s5Ara9C8aACNkpcXyefwSv1+BwSKijU0qpBidMBCIys0nRSBHBGPN8kGLqGtKHw7oXoXQfJPQiOzWWyloPB0qrSE+MDnV0SinVwJ8+gjE+0wRgLnC5PzsXkakislVE8kRkTjPLp4nIBhFZJyKrReScNsTesfUcbr3us5qHclJiAfham4eUUh2MP01Dd/jOi0gi8MKJthMRJ/AUcAFQCKwSkcXGmK98VlsKLDbGGBEZDrwGnNaG+DuunkMBsfoJBk0lJ9VKBDsOlzO+f0poY1NKKR8nc2dxBTDAj/XGAnnGmHxjTA3wKjDNdwVjTJkxxtizsdhXJnUJkfGQlGM9qAboER9FhMtB4ZGKEAemlFKN+dNH8E+++YJ2AIOxfrmfSG9gt898ITCumf1PBx4C0oBLW4jhVuBWgD59+vhx6A4ifTjsWQuAwyFkdIum8GhliINSSqnG/Ll89DGf93XATmNMoR/bNXdpzHG/+I0xbwJvishE4EHg/GbWeRp4GmD06NGdp9aQejpsWgS1leCOpnf3aHYf1RqBUqpj8aePYNlJ7rsQyPSZzwD2tnKc5SLST0RSjDGHT/KYHUvqQMDA4e2QPpyM7jFs2rs/1FEppVQjLfYRiEipiJQ0M5WKiD+D668CBohItohEAFdj3Zjme4z+IiL2+5FABFB08n9OB5MyyHo9tBWAjO7RHCmvoby6LoRBKaVUYy3WCIwx8aeyY2NMnYjMxhqwzgk8a4zZJCKz7OXzgCuAmSJSC1QC3/PpPO78kvuBOOCwlQgyk2IA2HOskoE9Tun0KqVUwPjTRwCAiKQBUfXzxphdJ9rGGPM28HaTsnk+7x8BHvE3hk7HFQndsxvVCAB2H6nQRKCU6jD8eR7B5SKyHfgaWAYUAO8EOa6uI3UQHN4GfJMI9MohpVRH4s99BA8CZwLbjDHZwHnAp0GNqitJHQRFeeCpJTUukkiXg0K9ckgp1YH4kwhqjTFFgENEHMaYj4Dc4IbVhaQMBG8dHC1ARMjorvcSKKU6Fn/6CI6JSBywHHhJRA5i3U+g/JGUY70e+RpSBpDRPUbvJVBKdSj+1AimYQ0r8WPgXWAH8K1gBtWl1CeCo18DkJkUze4jWiNQSnUc/tQIbgUW2ncT/y3I8XQ9sangjrVqBEBG9xiKK2spqaolIcod4uCUUsq/GkEC8J6IfCwiPxKRHsEOqksRgaRsOJIPQGZ3616CQq0VKKU6iBMmAmPMA8aYIcCPgF7AMhH5IOiRdSXdsxo1DQHaT6CU6jDaMgz1QWA/1hAQacEJp4tKyoGjO8Hr/aZGoFcOKaU6CH9uKLtdRP6N9RCZFOAHxpjhwQ6sS0nKBk81lO6lW4yb2Agnu/W5BEqpDsKfzuK+wF3GmHVBjqXr6p5tvR7JRxIzyEyK0ZvKlFIdhj99BHM0CZyipPpE8M2VQ3oJqVKqoziZR1WqtkrIAIerUYdx4dEKutJAq0qpzksTQXtwuiChNxyzntyZ0T2G8hoPRytqQxyYUkr511kcKyIO+/1AezRSvROqrbr1gWLrCZ+ZPsNRK6VUqPlTI1gORIlIb6wrh24EngtmUF1SYgYUWzWC+gfU6L0ESqmOwJ9EIMaYCuA7wB+NMdOBwcENqwtKzITSfeCp1ecSKKU6FL8SgYicBcwA3rLL/H6ymbJ1ywTjhZK9xEe56Rbj1qYhpVSH4E8iuAv4KfCm/czhHOCjoEbVFSVmWK9281Df5FgKispDGJBSSllO+MveGLMM6xGV2J3Gh40xdwY7sC4nsY/1al851D81jk/yDoUwIKWUsvhz1dDLIpIgIrHAV8BWEbkn+KF1MYm9rVf7yqH+aXEcKKmmpEovIVVKhZY/TUODjTElwLeBt4E+wHXBDKpLckdbzyYo3gVYiQAg72BZKKNSSim/EoHbvm/g28A/jDG1gN4SezISM79pGtJEoJTqIPxJBH8GCoBYYLmI9AVKghlUl5WY0eimsgingx2aCJRSIebPoHNPGGN6G2MuMZadwOR2iK3rqb+72BhcTgfZKbFaI1BKhZw/ncWJIvI7EVltT/+LVTtQbZWYCXWVUH4YsJqH8g5pIlBKhZY/TUPPAqXAVfZUAswPZlBdVv29BCXfXDm0+0gFVbWeEAallAp3/iSCfsaY+40x+fb0AJAT7MC6pIabyr5JBF4D+Yf0xjKlVOj4kwgqReSc+hkRORvQQXJORmKm9WonggE9rCuHth8sDVVESinlVyKYBTwlIgUiUgA8Cdzmz85FZKqIbBWRPBGZ08zyGSKywZ5WiMiINkXf2cQkgSu6IRFkp8TiEPTKIaVUSPkzxMR6YISIJNjzJSJyF7Chte1ExAk8BVwAFAKrRGSxMeYrn9W+BiYZY46KyMXA08C4k/pLOgORRsNRR7qc9E2OZbsmAqVUCPn9hDJjTIl9hzHA3X5sMhbIs/sVaoBXgWlN9rnCGHPUnv0cyPA3nk7L514CsPoJNBEopULpZB9VKX6s0xvY7TNfaJe15GbgnWYPJnJr/eWrhw518oHamiSCAWlxFBwup9bjDWFQSqlwdrKJwJ8hJppLFs1uJyKTsRLBvc0ezJinjTGjjTGjU1NT/Y+yI0rMhLIDUFcNWB3GdV7DTh2SWikVIi0mAhEpFZGSZqZSoJcf+y4EMn3mM4C9zRxnOPAXYJoxpqiN8Xc+DfcS7AGgf2o8ANsPaPOQUio0WuwsNsbEn+K+VwEDRCQb2ANcDVzru4KI9AHeAK4zxmw7xeN1Dr73EiTl0C/Nuklbh5pQSoVK0B45aYypE5HZwHuAE3jWfsLZLHv5POAXQDLwfyICUGeMGR2smDqEJjeVxUS4yOgerR3GSqmQCeqzh40xb2M9w8C3bJ7P+1uAW4IZQ4eT0PgBNWBdObTtgN5UppQKjZPtLFYnyx0FsWkN9xIA5KTEsbOoAmP0MQ9KqfaniSAUmlxC2jc5hspaD4fKqkMYlFIqXGkiCIXE3o0SQZ/kGAB2FVWEKiKlVBjTRBAKiZkND6gB6JNkJYKdmgiUUiGgiSAUEjOgtgIqrdE1MrpHIwK7jmgiUEq1P00EodDkEtJIl5NeidGaCJRSIaGJIBSaJAKwmod0mAmlVChoIgiFJg+oASsR7Dqiz/tRSrU/TQShEJMCzshG9xL0SY7hcFk15dV1IQxMKRWONBGEgsNx3CWkfesvIdV+AqVUO9NEECpNbiqrv4RUE4FSqr1pIgiV+nsJbH2TrFFItcNYKdXeNBGESmIGlO4DT601G+Ome4xbbypTSrU7TQShkpgBGCj55lk9fZNjNREopdqdJoJQaeZegqzkGAq0aUgp1c40EYRKM/cS9E2OZe+xSqrrPCEKSikVjjQRhEp9jeBoQUNRVkoMXgOFR/XGMqVU+9FEECruaKtWULS9oahvsnXlUMFhbR5SSrUfTQShlNwfivIaZrPqE4F2GCul2pEmglBKGQCH8xqeS9A9xk18lEvvJVBKtStNBKGU3B9qSqHsAAAiQlZyrNYIlFLtShNBKCX3t14P+/YT6HDUSqn2pYkglFIGWK9N+gkKj1ZS6/GGKCilVLjRRBBKCRngimqUCPomx+DxGvboJaRKqXaiiSCUHA6rVnDwq4ai7BTryqH8w2WhikopFWY0EYRaj2FwYFPDbP+0OADyDmoiUEq1D00EodZzqHXVUNkhALrFRJASF6GJQCnVbjQRhFqPodbrgS8bivqnxWkiUEq1m6AmAhGZKiJbRSRPROY0s/w0EflMRKpF5CfBjKXDqk8E+zc2FPVPi2P7wTKMfaOZUkoFU9ASgYg4gaeAi4HBwDUiMrjJakeAO4HHghVHhxebDPHpjfsJUuMorarjUGl1CANTSoWLYNYIxgJ5xph8Y0wN8CowzXcFY8xBY8wqoDaIcXR8PYbCgW9qBAN6xAPaYayUah/BTAS9gd0+84V2WZuJyK0islpEVh86dCggwXUoPYfBoS1QWwV8c+XQdk0ESql2EMxEIM2UnVSjtzHmaWPMaGPM6NTU1FMMqwPqdQZ46xpqBWnxkcRHurRGoJRqF8FMBIVAps98BrC3hXXDW++R1uuetYA1+Fz/HnrlkFKqfQQzEawCBohItohEAFcDi4N4vM4roTfEpsHetQ1F/VPjtGlIKdUugpYIjDF1wGzgPWAz8JoxZpOIzBKRWQAi0lNECoG7gZ+LSKGIJAQrpg5LxKoV7PFJBGlxHC6rprgivPvRlVLB5wrmzo0xbwNvNymb5/N+P1aTUVjKO5rHw6seJtYVy/09BpO07T2oLoXIeAb0qO8wLmV0VlKII1VKdWVhc2fxJ3s+Yfo/prO7ZPeJV24H+8v3c9uS29hctJmP93zMz8s2AQb2rgPg9HSrYrRpb0noglRKhYWwSQTlteXkHcujylMV6lCoqqviro/uoryunPlT53PHGXfw8bHNrIuMaOgn6JkQRXJsBBv3FIc4WqVUVxc2icApTgC8JrQPfCmqLOK2JbexqWgTD53zEAO7D+R7g75HnDuOhck9G105NLR3Ihu1RqCUCrKwSQQi1m0N5uRuZQiILw99yZX/vJKNhzfy24m/ZXKfyQDEuGO4OPti3o8UynyuHBraO4HtB0qpqvWEKmSlVBgIm0TgsP9UjwnNl2pBcQGzPphFpDOSly99manZUxstn95/OlUY3vUUQXkRAEN7JVLnNWzdXxqKkJVSYSJsEoHTYTUNhWJET4/Xw32f3gfAXy78C4OSBh23ztCUofSPSefNuLiGfoKhvRMB2LhX+wmUUsETNolA7BEvQtFH8MqWV9hwaANzxs4hI775q2VFhOmDrmJDVCR5O94DIKN7NInRbr4s1ESglAqesEkEDrH+1PZOBIcrD/Pkuic5u/fZXJZzWavrXjbwO7gMvLH/U8BKDiP7dGNlwZH2CFUpFabCJhHUdxa3dyJ4fM3j1Hhq+OnYnzbE0JKkqCQmR6XzL28JNdVWv8BZ/ZLJP1TOwZLQX/aqlOqawiYRhOLy0XUH17F4x2JmDp5J34S+fm1zZdZUjjodvLv+rwCclZMCwGf5RUGLUykV3sImEdQ3DbXX5aM1nhruX3E/PWN7cuvwW/3e7qwh3yenppYX8xdhjGFwrwQSolx8vP1wEKNVSoWzsEkE9Z3F7XX56J83/Jn84nzuP+t+Ytwxfm8n8WnM8ESxubqItQfX4nQI55/egw82H6DWE9qb4ZRSXVPYJIL6y0fbo2lo65GtPPvls1ze73LO6X1Om7f/VvrZJHi9vPTVCwBcPCydYxW1rNihzUNKqcALm0RQXyMI9n0Edd46/ufT/yExMpH/N+b/ndQ+ovufz5UlpSzd9SF7y/YyYUAKidFuXvp8Z4CjVUqpMEoE7XX56PNfPc/mI5u578z7SIxMPLmdZE/impJyBHhp80tEuZ3cMD6L9786wLYDepexUiqwNBEE0M6Snfzfuv/jvD7ncUHfC05+RzFJ9OwxnEu8USzYuoD95fu5YXwWsRFOHvzXV3i9oRsvSSnV9WgiCBBjDA989gARjgh+Nu5np77DflOYvedrjDE8vPJhusW4+dmlp/Px9sM88u6WkAyVoZTqmsIvERCcRPDG9jdYtX8Vd4++m7SYtFPfYb/J9Kqr4Ue9p7B011IWblvItWP7MGNcH/68PJ9b/raanUXlp34cpVTYC+qjKjuSYI41VFxdzO/W/I7RPUZzxYArArPTjLHgjuX6Cg8re5/Ng58/yJoDa/jlt35JTmocj767hcmP/ZsLB/fkWyN6MeW0NKIjnIE5tlIqrIRNIgjmncXz1s+jrLaMn4478TASfnNFQPYEnNuX8PjsLxj78jje/vptlhcuZ0TaCG66bBD5exNYue0I727aQ4TTzfCMREZnJTGqb3fO6NONlLjIwMSilOrSwiYRNNxZHOC29f3l+1mwdQHT+k1jYPeBAd03p10G294l+tA2NszcwGf7PmPJziWsO7iOz/Z+ZiW1dOjWy0Wcoxd7qnuw8asknlnTA291TzLiezGqTxIj+3ZnZJ/unNYzHpczbFoDlVJ+CptEUP9LPdB3Fj+z4RkMhlkjZgV0vwAMugTEAV/9A+mVy/he4xnfazwA1Z5qCooL2H5sO9uPbifvWB55R/Moc63CbW9eQiQflvbgnZVpeJf3xO3JYEjyMCYP6sWFg3vSPy0u8DErpTqdsEkEwbhqaPvR7byR9wbf6f8desX1Cth+G8QmQ865sOE1mPJzcHzTBxDpjGRQ0qDjHnJTVlNmJQV72n50O1uPbKe4ZjUAW42bjZuyeHxlf9IjhjN14BlcNKQXZ2R2w+EIULOWUqpTCZtEEOGIAKzB4AKhuLqYez++l3h3PLPPmB2QfTbrjO/D6zdB/kfQ//wTrh4XEUduWi65abmNyosqi9h4eCOf7/uc5YWfsqv0HY7yDi/vi+XlAzW4j17JeX2ncMmQbMb3SyHKrR3PSoWLsEkESdFJgPWgmFNRUVvBgq0LmL9xPqW1pTx13lN0j+oeiBCbd9plENcTPnoIcqaA4+Ta+JOjk5mUOYlJmZO4dywcrDjI5/s+58Ody1i6+33qkl/h3dIFvP1RBvJufwYlnsGkvqM4O6c3wzISiXRpYlCqqwqbRBDpjCTeHU9R1ckN3FZVV8VrW1/jrxv/ypGqI5zd62zuOOMOhqQMCXCkTbgi4bxfwD9+CBv/DsO/G5DdpsWkcXm/y7m83+XUemvZcGgDy3d/wkc7P6OgbBnb+ZBtu4Sn85LxVveimyuTXnHp9EvKYGBSJoNS+tA3KYGeiVFEuLQDWqnOTDrbHaqjR482q1evPqltp/59KrlpuTw84WG/t6n2VPP6ttf565d/5VDlIcalj+NHuT/ijLQzTiqGk+L1wjPnQnkRzF4FEf4Pa30yymvLWXtgLav3bWDl3i/5umQb5d5DjdYxRjB1cZi6RNzEE+mIJ8aVQLw7kcSIRJKiupMc04202O4kxyaSHJ1Ackw83aKiiY10ExvpJNrtDNzltkqpVonIGmPM6OaWhU2NACAxMpHiav8eBF/jqWFR3iKe3vA0ByoOMKrHKB6Z+Ahjeo4JcpTNcDjgoofguUvgg7lw8SMQxC/QWHcsEzImMCFjQkNZjaeG/eX72Vm8h62Hd7HjaCGFpXs5VHmA0toSqr2HOGJKKfJWQxXWdOz4fRsjYNwYrxuMC4eJwIEbh0TgIAKnROCSCNwSgdsRidsZSYQjkkhnBJHOaCJdkUQ7o4h2RRHtjiI+Mpr4iGjiIqKJjYgmITKG+Mgo4iNj6BYVS0JkNFFup3aEK9WK8EoEEYmU1JS0uk5BcQHvfP0OC7YuoKiqiNzUXH59zq8Z23NsaH+9Zp0N426HL/4Epftg+jyIiG23w0c4I+iT0Ic+CX2YkHlWi+vVeGo4Vn2Mosqj7CspYl/ZEYoqSimuLqWkupTy2koqa6uo9FRTXVdFZV0Vtd5qar011Jka6kwp1aaGSmrwemsw3lqM1IKc3GW/VuJxgXEjxoXgthIPETjF/U3icUQS4Ygm0hlNlCOaKFc00a5oYlwxxLhjiHXHEB8ZS0JkLN2j40iOiSctNoGU2DgSo6I00ahOLaiJQESmAn8AnMBfjDEPN1ku9vJLgArgBmPM2mDF0yuuF4t3LOaTPZ9QXltOaU0pJTUlHKs+RkFxAduObmNP2R4AJvSewHWDr+PM9DM7TvPFRb+BuipYMx/yl8GwK6HfZEg9HRLS2zUxtCTCGUFaTBppMWmcnhy4/Xq8Hqo91VTUVlJWU0lJdQUlVZUcqyynuLqS8tpKymsqqaitsqa6KiprK6j21FBVV0WNp4ZqbzW1nmpqvLXUmWrqvDV4TDWVppQybw1eqjFSjZEaRPy/zNh4XYiJxGGicRJDhMQS6YwlyhlHrCuOOHc88RHxJEYm0D06kaToRFJjEukR250e8d1IjokjLtLVcf6dqbATtD4CEXEC24ALgEJgFXCNMeYrn3UuAe7ASgTjgD8YY8a1tt9T6SNYsXcFty257bhyl7jITMikf7f+jOk5hnMzziU9Lv2kjtEudq+ClU/D5sVWYqgXmQgxSRARZyWFiFirP8EZaXU6OyOsyWW/OiPB6W6yzC5zRvqsG9lkuwifbdzWTW+I9Spy/Dxil3eOLzpjDLXeWoqryzhaWcqxqjKKKysoqSnjWFUZxyrLKakup6S6jNKacsprK6ioq6CirowqTxnV3nJqTQUeKjBSCY7aExzPCd5IHCYSp0ThIgqXRBHhiCbCEUmMK5EYVzzRrkhi3G4iXBE4xBAfkUhiRDfcTgdRLjfRbhcuhxOHOHA4HDhw4HI6cYoDhwhOceF0CE5x4nQ4cDmciL2OQxy4xInLITgc1jYu+9XpcOJyOKx1HD7bijWGl0PAIWJ/xJ3jMw5HoeojGAvkGWPy7SBeBaYBX/msMw143ljZ6HMR6SYi6caYfcEIaHyv8Sz81kJKa0pJjEwkISKBhIgEol3RnesfcOYYa6p9Ag58BUXboWSv1WRUeQxqyqGmDKqKrXJPNdTVgKem8ft2en5zI02ThMNFQ6LAN1lIo5fjZpqu16jMd/tmyk5QLkAEkIqQ2uzfIE0DayEe66UKKBUoBkp8XkuAEjGUSh0VUkc5ZVSIUOGASrGm3W5w1kFNB2x6ivYaIrxg7NCMPQHWGL/SuMz3fYJPheu4daTJfnyWOwC3AZc57hNoM9/tTTNlgRDo/Y2PGM79M18O8F6Dmwh6A7t95guxfvWfaJ3eQKNEICK3ArcC9OnT55SCOi3ptFPavkNxR0PGKGs6GV6PlRDqqsFTaycJn/eeWnu+xme9Gp95+xUDxgvGfm2Yp8m8OX5db50Vi/H5OvB9X69RzdW0XNaovLmykymvXyZ2/DSzrOUYo4AoIPUE6x3vmzKPMVR7PVR6vdR4vZR56zhq6vAYgwdDrddQ5/XiFWtdL+DFfrXnPRgMBg9Yr8Z+lW/WadimxVdrXY9ABR5qMTRNgWKsOd9ysf+cauNlHzV0F3fjZY3eS8P4+NKk3IuhVgx1x38Ifqup8+J0CI4mP/7qE1rTL2/T8J8TCNC3/nG7sf/ZeY0hJaZnYA7SRDATQXOnpenp9GcdjDFPA0+D1TR06qEpwBqywhFtJRTVoTmBGHtSKtCCeSdQIZDpM58B7D2JdZRSSgVRMBPBKmCAiGSLSARwNbC4yTqLgZliORMoDlb/gFJKqeYFrWnIGFMnIrOB97Bqts8aYzaJyCx7+TzgbawrhvKwLh+9MVjxKKWUal5Q7yMwxryN9WXvWzbP570BfhTMGJRSSrVORwtTSqkwp4lAKaXCnCYCpZQKc5oIlFIqzHW65xGIyCFg50lungKc2iPKgqOjxgUdNzaNq200rrbpinH1NcY0P2pKZ0sEp0JEVrc06FIoddS4oOPGpnG1jcbVNuEWlzYNKaVUmNNEoJRSYS7cEsHToQ6gBR01Lui4sWlcbaNxtU1YxRVWfQRKKaWOF241AqWUUk1oIlBKqTAXNolARKaKyFYRyROROUE+VqaIfCQim0Vkk4j8l10+V0T2iMg6e7rEZ5uf2rFtFZGLfMpHiciX9rInJADP1BSRAnuf60RktV2WJCJLRGS7/dq9PWMTkUE+52WdiJSIyF2hOGci8qyIHBSRjT5lATs/IhIpIgvs8i9EJOsU4vqtiGwRkQ0i8qaIdLPLs0Sk0ue8zfPZpj3iCtjnFuC4FvjEVCAi60Jwvlr6fgjdvzFjTJefsIbB3gHkYD2Odj0wOIjHSwdG2u/jgW3AYGAu8JNm1h9sxxQJZNuxOu1lK4GzsJ7m9g5wcQDiKwBSmpQ9Csyx388BHglFbD6f136gbyjOGTARGAlsDMb5AX4IzLPfXw0sOIW4LgRc9vtHfOLK8l2vyX7aI66AfW6BjKvJ8v8FfhGC89XS90PI/o2FS41gLJBnjMk3xtQArwLTgnUwY8w+Y8xa+30psBnrWcwtmQa8aoypNsZ8jfV8hrEikg4kGGM+M9Yn+jzw7SCFPQ34m/3+bz7HCUVs5wE7jDGt3UEetLiMMcuBI80cL1Dnx3dfrwPn+VNraS4uY8z7xhj7wc98jvWUvxa1V1ytCOn5qmdvfxXwSmv7CFJcLX0/hOzfWLgkgt7Abp/5Qlr/Yg4Yu0p2BvCFXTTbrsY/61P1aym+3vb7puWnygDvi8gaEbnVLuth7KfD2a9pIYoNrF8wvv+DdoRzFsjz07CN/SVeDCQHIMabsH4V1ssWkf+IyDIRmeBz7PaKK1CfWzDO1wTggDFmu09Zu5+vJt8PIfs3Fi6JoLlMGPTrZkUkDvg7cJcxpgT4E9APyAX2YVVNW4svWHGfbYwZCVwM/EhEJraybrvGJtZjTS8HFtpFHeWcteRk4gh4jCJyH1AHvGQX7QP6GGPOAO4GXhaRhHaMK5CfWzA+02to/GOj3c9XM98PLa7awnECFlu4JIJCINNnPgPYG8wDiogb60N+yRjzBoAx5oAxxmOM8QLPYDVZtRZfIY2r+gGJ2xiz1349CLxpx3HArmrWV4cPhiI2rOS01hhzwI6xQ5wzAnt+GrYREReQiP9NK8cRkeuBy4AZdhMBdjNCkf1+DVa78sD2iivAn1ugz5cL+A6wwCfedj1fzX0/EMJ/Y+GSCFYBA0Qk2/7FeTWwOFgHs9vi/gpsNsb8zqc83We16UD91QyLgavtnv5sYACw0q4elorImfY+ZwL/OMXYYkUkvv49VmfjRjuG6+3Vrvc5TrvFZmv0S60jnDOf4wXq/Pju60rgw/ov8LYSkanAvcDlxpgKn/JUEXHa73PsuPLbMa5Afm4Bi8t2PrDFGNPQrNKe56ul7wdC+W+stZ7krjQBl2D1zu8A7gvysc7BqoZtANbZ0yXAC8CXdvliIN1nm/vs2Lbic5ULMBrrf6IdwJPYd4OfQmw5WFcgrAc21Z8LrPbDpcB2+zUpBLHFAEVAok9Zu58zrES0D6jF+mV1cyDPDxCF1fSVh3XVR84pxJWH1RZc/++s/kqRK+zPdz2wFvhWO8cVsM8tkHHZ5c8Bs5qs257nq6Xvh5D9G9MhJpRSKsyFS9OQUkqpFmgiUEqpMKeJQCmlwpwmAqWUCnOaCJRSKsxpIlCqCRHxSOORUAM2Wq1Yo1xuPPGaSrUfV6gDUKoDqjTG5IY6CKXai9YIlPKTWOPXPyIiK+2pv13eV0SW2gOsLRWRPnZ5D7GeEbDensbbu3KKyDNijUX/vohEh+yPUgpNBEo1J7pJ09D3fJaVGGPGYt3F+Xu77EngeWPMcKxB356wy58AlhljRmCNi7/JLh8APGWMGQIcw7qrVamQ0TuLlWpCRMqMMXHNlBcAU4wx+fagYfuNMckichhrCIVau3yfMSZFRA4BGcaYap99ZAFLjDED7Pl7Abcx5lft8Kcp1SytESjVNqaF9y2t05xqn/cetK9OhZgmAqXa5ns+r5/Z71dgjWgLMAP4xH6/FLgdQESc9vj2SnU4+ktEqeNFi/1Qc9u7xpj6S0gjReQLrB9R19hldwLPisg9wCHgRrv8v4CnReRmrF/+t2ONhqlUh6J9BEr5ye4jGG2MORzqWJQKJG0aUkqpMKc1AqWUCnNaI1BKqTCniUAppcKcJgKllApzmgiUUirMaSJQSqkw9/8BL9zpuKqKjkUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the losses\n",
    "\n",
    "losses_np = []\n",
    "losses_data_np = []\n",
    "losses_physics_np = []\n",
    "\n",
    "for i in range(len(losses)):\n",
    "    losses_np.append(losses[i].item())\n",
    "    losses_data_np.append(losses_data[i].item())\n",
    "    losses_physics_np.append(losses_physics[i].item())\n",
    "\n",
    "plt.plot(losses_np, label='Combined loss')\n",
    "plt.plot(losses_data_np, label='Data loss')\n",
    "plt.plot(losses_physics_np, label='Physics loss')\n",
    "plt.title(\"Losses over epochs\")\n",
    "plt.xlabel(\"Epoch\")\n",
    "plt.ylabel(\"Loss value\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc762fe4",
   "metadata": {},
   "source": [
    "### Critically damped case\n",
    "\n",
    "Now let's extend this solution to the critically damped case. Contrary to the underdamped setting, in the critically damped case we have $ \\delta = \\omega_0 $. The analytical solution for the critically damped case DE is:\n",
    "\n",
    "$$ x(t) = e^{-\\delta t} (C_1e +t \\ C_2) $$\n",
    "\n",
    "where $C_{1,2}$ are integration constants that can be determined from the given initial conditions (see [here](https://beltoforion.de/en/harmonic_oscillator/) for a detailed solution of the DE).\n",
    "\n",
    "For this case, we consider the material/geometry parameters $\\delta = 20$ and $\\omega_0=20$, and try to learn the solution for the range $x \\in [0,1]$. In addition, we assume the initial conditions $x_0 = 0$ and $v_0 = 2$ in order to calculate c1 and c2.\n",
    "\n",
    "In the critically damped and overdamped cases, c1 and c2 are calculated by solving the equation system of $x(t)$ and $v(t) = \\dot{x}(t)$ for the initial conditions for starting position $x_0$ and starting velocity $v_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3e6931ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "def CD_oscillator(d, w0, x, x0, v0):\n",
    "    \"\"\"Defines the analytical solution to the 1D critically damped harmonic oscillator problem. \n",
    "    Equations from: https://beltoforion.de/en/harmonic_oscillator/\"\"\"\n",
    "    assert d == w0\n",
    "    \n",
    "    c1 = x0               # c1,2 are computed using the predefined\n",
    "    c2 = v0 + (d*x0)      # initial conditions x0 and v0\n",
    "    summ = (c1 + x*c2)    # c1 + t*c2\n",
    "    exp = torch.exp(-d*x) # exp(-d*t)\n",
    "    y = exp*summ          # y - exact solution\n",
    "    \n",
    "    return y\n",
    "\n",
    "# plot_result redefined with different x,y limits and text position\n",
    "def plot_result(epoch,x,y,x_data,y_data,yh,xp=None):\n",
    "    \"Pretty plot for training results\"\n",
    "    plt.figure(figsize=(8,4))\n",
    "    plt.plot(x,y, color=\"grey\", linewidth=2, alpha=0.8, label=\"Exact solution\")\n",
    "    plt.plot(x,yh, color=\"tab:blue\", linewidth=4, alpha=0.8, label=\"Neural network prediction\")\n",
    "    plt.scatter(x_data, y_data, s=60, color=\"tab:orange\", alpha=0.4, label='Training data')\n",
    "    if xp is not None:\n",
    "        plt.scatter(xp, -0*torch.ones_like(xp), s=60, color=\"tab:green\", alpha=0.4, \n",
    "                    label='Physics loss training locations')\n",
    "    l = plt.legend(loc=(1.01,0.34), frameon=False, fontsize=\"large\")\n",
    "    plt.setp(l.get_texts(), color=\"k\")\n",
    "    plt.xlim(-0.05, 1.05)\n",
    "    plt.ylim(-0.05, 0.05)\n",
    "    plt.text(1.065, 0.1,\"Training step: %i\"%(epoch),fontsize=\"xx-large\",color=\"k\")\n",
    "    plt.axis(\"off\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71ffcbd4",
   "metadata": {},
   "source": [
    "Let's start by generating the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9dda5529",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEWCAYAAABxMXBSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4CklEQVR4nO3de3xU1bnw8d+TyT2BJIRwkURuRa6GgBEQL4WiFrSW2ldRa+utLaJSPLbnbW3fc6z17cV66qW8tlrxWPWcHoUetcUW61Gq4g0FFCLILY0IAeROCLlP5nn/2DtxGCfJDJmdmSTP9/OZz8xee+291x7CfmattfdaoqoYY4wxkUqKdwGMMcZ0LxY4jDHGRMUChzHGmKhY4DDGGBMVCxzGGGOiYoHDGGNMVCxwGBNDInKdiLzh8TE2iciMWOftDBEZJiIqIsleH8vEnwUO4zkRuVJE3hGRGhHZ736+WUQk3mULJSKvisi3PNp3TC6uqjpeVV+Ndd6uIiIzRKQy3uUwJ88Ch/GUiHwP+DXwb8AgYCCwADgbSO3isiT8r+HuUEZjLHAYz4hIDnAXcLOq/reqVqvjfVW9WlUb3HxpIvIrEdkpIvtE5GERyXDXzRCRShH5nltb2Ssi1wcdI5JtfyAinwC/F5E8EfmLiBwQkSPu50I3/8+Ac4EHReS4iDzopo8RkZdE5LCIbBWReUHHzxeR5SJyTETeBUa285Wsct+Puvs/y23aelNE7heRw8CdIjJSRP4uIodE5KCI/EFEcoOOuUNEznc/3ykiy0TkSRGpdpumSk8y72QRed9d90cRWSoiP23j39bnfu8HRaQCuDhk/fUistndV4WI3OimZwEvAKe438FxETlFRKaIyNsictT9N35QRLr0h4WJnAUO46WzgDTgzx3k+yVwGlACfA4YAtwRtH4QkOOmfxP4jYjkRbFtP2AoMB/nb/737vKpQB3wIICq/h/gdWChqmar6kL3QvcS8F/AAOAq4LciMt7d/2+AemAwcIP7ast57nuuu/+33eWpQIW7/58BAvwCOAUYCxQBd7az3y8DTwO5wPKW84kmr3uRfg54HOf7egq4tJ39fBv4EjAJKAUuC1m/313fF7geuF9EJqtqDTAH2ON+B9mqugdoBm4D+uP83cwCbm7n+CaeVNVe9vLkBXwd+CQk7S3gKM4F+zyci2QNMDIoz1nAR+7nGW7e5KD1+4FpEW7bCKS3U8YS4EjQ8qvAt4KWrwBeD9nmd8CPAR/QBIwJWvdz4I02jjUM0JBzuQ7Y2cH3+BXg/aDlHcD57uc7gZeD1o0D6qLN6/5b7AYkaP0bwE/bKNPfgQVByxeGnltI/j8Btwb9u1R2cM7/BDwX779he4V/WXuq8dIhoL+IJKuqH0BVpwO4naNJQAGQCawL6isXnIty635atnfVAtkRbntAVetbV4pkAvcDs4GWWksfEfGpanOYcxgKTBWRo0FpycB/uMdPBnYFrfs47DfRvuDtEZEBwGKcZrM+ON/TkXa2/yTocy2QHvydR5IXp3azW92rdrhyhTiFds5bRObgBNfT3PJnAh+0tTMROQ24D6f2konzva5r5/gmjqypynjpbaABmNtOnoM4NYrxqprrvnJUNTuC/Ueybejwz98DRgNTVbUvnzYfSRv5dwGvBe0/V53mlZuAA4AfpympxantlLetoahD03/hphW7Zfx6UPm8shcYEnKnW1Fbmd38Yc9bRNKAZ4BfAQNVNRdYQdvfMcBDwBZglHvOP8L7czYnyQKH8YyqHgV+gtMncJmIZItIkoiUAFlungCwBKcNfACAiAwRkS9GsP+T2bYPTrA5KiL9cH4VB9sHjAha/gtwmoh8Q0RS3NeZIjLWraE8i9OhnSki44Br2zn2ASAQsv+2ynjcLeMQ4H93kD8W3sbpZ1goIskiMheY0k7+ZcAiESl0+5tuD1qXitO3dQDwu7WPC4PW7wPyxbl5okUf4BhwXETGADd1+oyMZyxwGE+p6j3Ad4Hv4/RN7MPpI/gBTn8H7udyYLWIHANexqkVRCLabR8AMnBqK6uBv4Ws/zVwmXvH1WJVrca56F0J7MFp6vklzoURYCFOs9knOB3Lv2/rwKpai9P5/aZ799C0NrL+BJgMVAF/xQlOnlLVRuCrODcfHMWp5fwFp8YYzhLgRWAD8F5wGd3vbBFOcDkCfA2nI75l/RaczvcK93s4BfhnN1+1u++lsTs7E2tyYpOmMcY4ROQd4GFVbTMYmt7JahzGGABE5PMiMshtqroWKOazNTJj7K4qY0yr0TjNS9nAP4DLVHVvfItkEpE1VRljjImKNVUZY4yJSq9oqurfv78OGzYs3sUwxphuZd26dQdVtSA0vVcEjmHDhrF27dp4F8MYY7oVEQk7EoI1VRljjImKBQ5jjDFRscBhjDEmKr2ij8MYEx9NTU1UVlZSX1/fcWYTN+np6RQWFpKSkhJRfgscxhjPVFZW0qdPH4YNG4Yk3hTzBmdOpkOHDlFZWcnw4cMj2saaqowxnqmvryc/P9+CRgITEfLz86OqFVrgMMZ4yoJG4ov238gCh0dWVxzid6/9g+37quNdFGOMiSkLHB54/M2PuPKR1fzihS1cvPgNXvpwX7yLZEyv5fP5KCkpaX3dfffdMdv3+vXrWbFiRUz2lZ3d8aSXDzzwALW1ta3LF110EUePHo3J8aNhgSPGPjpYw89WbOb8sQN4/fszGTu4D//09PvsPloX76IZ0ytlZGSwfv361tftt9/e8UYRimXgiERo4FixYgW5ublddvwWFjhipWwZ3D+Bh+6/E19zA78Y8xFF/TJ58GuTCSj829+2xLuExhhXVVUVo0ePZuvWrQBcddVVLFmyBICbbrqJ0tJSxo8fz49//OnMwmvWrGH69OlMnDiRKVOmUFVVxR133MHSpUspKSlh6dITJy3ctGkTU6ZMoaSkhOLiYrZv3w7Afffdx4QJE5gwYQIPPPDAZ8r26quv8qUvfal1eeHChTz++OMsXryYPXv2MHPmTGbOnAk4wykdPHiwzf3u2LGDsWPH8u1vf5vx48dz4YUXUlfX+R+xdjtuLJQtg+cXUdMY4PnmaVzqe52Cl/8LMqCoeB5fn3Yqj725g+9dOJqifpnxLq0xcfGT5zfx4Z5jMd3nuFP68uNLxrebp66ujpKSktblH/7wh1xxxRU8+OCDXHfdddx6660cOXKEb3/72wD87Gc/o1+/fjQ3NzNr1izKysoYM2YMV1xxBUuXLuXMM8/k2LFjZGZmctddd7F27VoefPDBzxz34Ycf5tZbb+Xqq6+msbGR5uZm1q1bx+9//3veeecdVJWpU6fy+c9/nkmTJnV4rosWLeK+++7jlVdeoX///iesa2u/eXl5bN++naeeeoolS5Ywb948nnnmGb7+9a9H8O22zWocsbDyLmiq46VAKXWk81XfG9BU56QD3zxnBEkCj6yqiHNBjel9QpuqrrjiCgAuuOACTj/9dG655RYeffTR1vzLli1j8uTJTJo0iU2bNvHhhx+ydetWBg8ezJlnnglA3759SU5u/3f3WWedxc9//nN++ctf8vHHH5ORkcEbb7zBpZdeSlZWFtnZ2Xz1q1/l9ddf7/Q5trff4cOHtwbOM844gx07dnT6eFbjiIWqSgBeaS6hP1WcIdtOSB+Uk86XJw7h2fcq+dFFY8lI9cWrpMbETUc1g64WCATYvHkzGRkZHD58mMLCQj766CN+9atfsWbNGvLy8rjuuuuor69HVaO+ZfVrX/saU6dO5a9//Stf/OIXefTRR4lk4rzk5GQCgUDrciTPV7S337S0tNbPPp8vJk1VVuOIhZxCAiq8Hjid85LKSBJtTW9x2RmF1DQ289Jmu8PKmERw//33M3bsWJ566iluuOEGmpqaOHbsGFlZWeTk5LBv3z5eeOEFAMaMGcOePXtYs2YNANXV1fj9fvr06UN1dfhb7isqKhgxYgSLFi3iy1/+MmVlZZx33nn86U9/ora2lpqaGp577jnOPffcE7YbOnQoH374IQ0NDVRVVbFy5crWdW0dL5L9xpKngUNEZovIVhEpF5HP3MogjsXu+jIRmeymp4vIuyKyQUQ2ichPgra5U0R2i8h693WRl+cQkVl3UO4bwWH6clbSJictJQNm3dGaZerwfpySk85z71XGqZDG9E4tfRwtr9tvv51t27bx6KOPcu+993Luuedy3nnn8dOf/pSJEycyadIkxo8fzw033MDZZ58NQGpqKkuXLuU73/kOEydO5IILLqC+vp6ZM2fy4Ycfhu0cX7p0KRMmTKCkpIQtW7ZwzTXXMHnyZK677jqmTJnC1KlT+da3vvWZ/o2ioiLmzZtHcXExV1999Qnr58+fz5w5c1o7x1tEst9Y8mzOcRHxAduAC4BKYA1wlap+GJTnIuA7wEXAVODXqjpVnDphlqoeF5EU4A3gVlVdLSJ3AsdV9VeRlqW0tFS9nsjpv559lh+9m8Yrqd9leF6KEzSK552Q55d/28IjqypY+3/OJy8r1dPyGJMINm/ezNixY+NdDBOBcP9WIrJOVUtD83pZ45gClKtqhao2Ak8Dc0PyzAWeVMdqIFdEBrvLx908Ke7LmwgXI+83j6RfVirDfrIFbtv4maABMGfCIJoDyitb98ehhMYYExteBo4hwK6g5Uo3LaI8IuITkfXAfuAlVX0nKN9Ct2nrMRHJC3dwEZkvImtFZO2BAwc6eSod27qvmrGD+7TbgTbhlBwG9k3jZevnMMZ0Y14GjnBX0NBaQ5t5VLVZVUuAQmCKiExw1z8EjARKgL3AveEOrqqPqGqpqpYWFHxmrvWYag4o2/ZVM3pg33bzJSUJs8YO5LWtB2jwN3taJmOM8YqXgaMSKApaLgT2RJtHVY8CrwKz3eV9blAJAEtwmsTiaufhWuqbAowZ1KfDvOePHUBNYzPvVBzugpIZY0zseRk41gCjRGS4iKQCVwLLQ/IsB65x766aBlSp6l4RKRCRXAARyQDOB7a4y4ODtr8U2OjhOURk6yfO07CjIwgc00bkk+IT3iw/6HWxjDHGE549AKiqfhFZCLwI+IDHVHWTiCxw1z8MrMC5o6ocqAWudzcfDDzh3pmVBCxT1b+46+4RkRKcJq0dwI1enUOktnxSjQicNrDjwJGZmszkU/N4wwKHMaab8vQ5DlVdoaqnqepIVf2Zm/awGzRw7566xV1/uqquddPLVHWSqhar6gRVvSton99w8xar6pdVda+X5xCJbfuqGdovM+Inws/5XH827TnG4ZpGj0tmTO926NCh1uc3Bg0axJAhQ1qXGxvb//+3du1aFi1a1OExpk+fHqvinmDGjBl09BhB6Gi5XcWeHI+B7fuOR1TbaHH2KGeAsrf+YbUOY7yUn5/fOkbVggULuO2221qXU1NT8fv9bW5bWlrK4sWLOzzGW2+9FcsiR8UCRzelquw8XMvQ/MhHvS0ekkOftGTeLD/kYcmM6Ybc6Qm4M9d5L1sW80Ncd911fPe732XmzJn84Ac/4N1332X69OlMmjSJ6dOntw61Hjy8+Z133skNN9zAjBkzGDFixAkBpWUCpldffZUZM2Zw2WWXMWbMGK6++urWMaRWrFjBmDFjOOecc1i0aNEJw6a3qKur48orr6S4uJgrrrjihDGlwg31Hm6Y9baGhI81G+SwM8qWceB/7qfB/y8Urb8fii4I++BfqGRfElNH5FsHuTHB3OkJaHIvmFW7nGWI6P9VNLZt28bLL7+Mz+fj2LFjrFq1iuTkZF5++WV+9KMf8cwzz3xmmy1btvDKK69QXV3N6NGjuemmm0hJSTkhz/vvv8+mTZs45ZRTOPvss3nzzTcpLS3lxhtvZNWqVQwfPpyrrroqbJkeeughMjMzKSsro6ysjMmTJ7euCzfUe7hh1sPlKy4ujuE357Aax8ly/8h3HWsCoKh+m/NHHuEvpLM/l8/Ow7XsOtz11UxjEpI7PcEJgqYniKXLL78cn8/pk6yqquLyyy9nwoQJ3HbbbWzatCnsNhdffDFpaWn079+fAQMGsG/fZx/knTJlCoWFhSQlJVFSUsKOHTvYsmULI0aMYPjw4QBtBo5Vq1a1zpNRXFx8wgU/3FDv4USar7MscJws9498lw4AoEj2R/VHPmV4PwDWfXzEsyIa061UtTEAaFvpnZCVldX6+V//9V+ZOXMmGzdu5Pnnn29zGPPQ4cnD9Y+EyxPNeIDhRp5oGep95cqVlJWVcfHFF4ctY6T5YsECx8ly/5h3qfNUeqEcOCG9I2MG9SU7LZk1O+xBQGOAE6YhiCg9RqqqqhgyxBkN6fHHH4/5/seMGUNFRUXrBEqho+i2OO+88/jDH/4AwMaNGykrKwNoc6h3OHGY9fbyxZoFjpPl/jHv0gIKOEK6NJ2Q3hFfkjB5aB5rd1iNwxjAGVE6JePEtJDpCbzw/e9/nx/+8IecffbZNDfHfiigjIwMfvvb3zJ79mzOOeccBg4cSE5Ozmfy3XTTTRw/fpzi4mLuuecepkxxBsVoa6h3OHGY9fbyxZpnw6onEk+GVXf7OK6s+S6NmsKzaXc6f+SXLI64I+//rdzOvS9tY8MdF5KTmdLxBsZ0M1EPq162zGnurap0foSFmZ6gOzp+/DjZ2dmoKrfccgujRo3itttui3exThDNsOp2V9XJcv+Ydz3dQKlshpyiqP/IS4e5/Rw7D/OFMQM9KaYx3UrxvB4RKEItWbKEJ554gsbGRiZNmsSNN8Z9wItOscDRCU3jL2OvvkDRjGvhiz+PevuSolySk4Q1O45Y4DCmB7vtttsSrobRGdbH0QkHqhsIKJySm9Fx5jAyUn1MGJLDWusgNz1Yb2gO7+6i/TeywNEJ+6sbABjQJ62DnG07c1geGyqrbH4O0yOlp6dz6NAhCx4JTFU5dOgQ6enpEW9jTVWdsP+Yc4/0wL6Rf+GhSof1Y8nrH7FxdxVnDO0Xq6IZkxAKCwuprKykK2bhNCcvPT2dwsLIb3u2wNEJ+1pqHH1PvsZROtSZ+fbdj45Y4DA9TkpKSusT06bnsKaqTjhwrJ4kgfys1JPeR352GsPyM3l/pz3PYYzpHixwdMK+Yw3kZ6eR7Ovc11hSlMv6XUetHdgY0y1Y4OiE/dX1DOxEM1WLkqJc9lc3sLfKm3FljDEmlixwdML+6gYG9Dn5jvEWJac6/Rzrdx3t9L6MMcZrngYOEZktIltFpFxEbg+zXkRksbu+TEQmu+npIvKuiGwQkU0i8pOgbfqJyEsist19z/PyHNqz71hDp27FbTF2cB9SfUkWOIwx3YJngUNEfMBvgDnAOOAqERkXkm0OMMp9zQcectMbgC+o6kSgBJgtItPcdbcDK1V1FLDSXe5y/uYAh2oaGNCJW3FbpCX7GHdKX9bvPNr5ghljjMe8rHFMAcpVtUJVG4GngbkheeYCT6pjNZArIoPd5eNunhT3pUHbPOF+fgL4iofn0KaDxxtR7dzDf8FKinL5YHcV/uZATPZnjDFe8TJwDAF2BS1XumkR5RERn4isB/YDL6nqO26egaq6F8B9HxDu4CIyX0TWishaLx4+2l/d+Yf/gpUU5VLX1MzWfdUx2Z8xxnjFy8Dx2amsPq01dJhHVZtVtQQoBKaIyIRoDq6qj6hqqaqWFhQURLNpRA4edx7+65998s9wBCspygVgw66qmOzPGGO84mXgqASKgpYLgT3R5lHVo8CrwGw3aZ+IDAZw3/fHrMRROFzjTNyUnxWbpqqh+ZnkZaawfpc9CGiMSWxeBo41wCgRGS4iqcCVwPKQPMuBa9y7q6YBVaq6V0QKRCQXQEQygPOBLUHbXOt+vhb4s4fn0KYjNY0A5GXFZgImEWGi+yCgMcYkMs/GqlJVv4gsBF4EfMBjqrpJRBa46x8GVgAXAeVALXC9u/lg4An3zqwkYJmq/sVddzewTES+CewELvfqHNpzuLaRFJ+QnRa7r7CkKJfXth2gur6JPuk2I6AxJjF5Osihqq7ACQ7BaQ8HfVbgljDblQGT2tjnIWBWbEsavcPHG+mXlYpIuG6ak1NSlIsqfFBZxfTP9Y/Zfo0xJpbsyfGTdLi2kbzM2HSMt2jpIH/fmquMMQnMAsdJOlLj1DhiKTczleH9s6yfwxiT0CxwnKTDtY3kxThwgI2Ua4xJfBY4TtKRmkb6xbipCpzAcaC6gT02Uq4xJkFZ4DgJzQHlaF2TZzUOgA3WXGWMSVAWOE5CVV0TqtAvM/a3zI6xkXKNMQnOAsdJOFzjDDfiRY2jdaRcCxzGmARlgeMkxHq4kVAlRbl8UGkj5RpjEpMFjmiVLePw0zcBkPfcVVC2LOaHaBkpd/v+4x1nNsaYLmaBIxply+D5RRxxm6r61ZTD84tiHjwmWge5MSaBWeCIxsq7oKmOw/QBII9qaKpz0mNoWH4mORkp1s9hjElIFjiiUVUJwDHNIo1G0qXphPRYsZFyjTGJzAJHNHIKAThGFn2p+Ux6LJUU5rBtXzU1Df6Y79sYYzrDAkc0Zt0BKRkc00z6Sq2TlpLhpMdYyam5BBQ27rYZAY0xicUCRzSK58Eli6lKzieHGsgpgksWO+kxNrEwF4ANlUdjvm9jjOkMT+fj6JGK53Fs1RvOyLjX3+bZYfKz0yjql2H9HMaYhGM1jpNwrK6Jvl0wQ9/Ewlw27LKmKmNMYvE0cIjIbBHZKiLlInJ7mPUiIovd9WUiMtlNLxKRV0Rks4hsEpFbg7a5U0R2i8h693WRl+cQzrF6P30zvK+slRTlsvtoHfurbaRcY0zi8CxwuPOF/waYA4wDrhKRcSHZ5gCj3Nd84CE33Q98T1XHAtOAW0K2vV9VS9zXCVPTek1VqaprIifD+xrHpyPlWq3DGJM4vKxxTAHKVbVCVRuBp4G5IXnmAk+qYzWQKyKDVXWvqr4HoKrVwGZgiIdljVhtYzPNAe2Spqrxp+TgSxJ7gtwYk1C8DBxDgF1By5V89uLfYR4RGQZMAt4JSl7oNm09JiJ5MStxBI7VOw/99e2CGkdGqo8xg/pYB7kxJqF4GTgkTFrofKjt5hGRbOAZ4J9U9Zib/BAwEigB9gL3hj24yHwRWSsiaw8cOBBl0dt2rM55IK8rahzgjFu1ofIogYBNJWuMSQxeBo5KoChouRDYE2keEUnBCRp/UNVnWzKo6j5VbVbVALAEp0nsM1T1EVUtVdXSgoKCTp9Mi5YaR1f0cYDTz1Fd76fiYE3HmY0xpgt4GTjWAKNEZLiIpAJXAstD8iwHrnHvrpoGVKnqXhER4N+Bzap6X/AGIjI4aPFSYKN3p/BZVbUtTVVd8wiMTSVrjEk0ngUOVfUDC4EXcTq3l6nqJhFZICIL3GwrgAqgHKf2cLObfjbwDeALYW67vUdEPhCRMmAm4N1TeGG09nF0UVPVyIJsslJ91s9hjEkYnv5sdm+VXRGS9nDQZwVuCbPdG4Tv/0BVvxHjYkblWF3XdY4D+JKE4sJcG3rEGJMw7MnxKB2rb+kc77rRWiYW5bJ57zHqm5q77JjGGNMWCxxRqqprIivVR7Kv6766kqJcmpqVD/ce6zizMcZ4zAJHlI7VNXVZM1UL6yA3xiQSCxxROlbfNQMcBhuUk86gvunWQW6MSQgWOKJUXe+nTxf2b7SYWJRjNQ5jTEKwwBGlmgY/WWnxCBy57DhUy9Haxi4/tjHGBLPAEaXjDX6y4xA4Wvo5rLnKGBNvFjiiFK/AcfqQHERsiHVjTPxZ4IhSTUNzXJqq+qSnMGpANut3HenyYxtjTDALHFEIBJSaRj/Zab64HH9iYS7rdx3FeeDeGGPiwwJHFGqbmlGF7DjcVQVwxtA8jtQ28ZGNlGuMiSMLHFGoaXCGG4lHUxVA6TBnzqq1H1tzlTEmfixwROG4Gzji0TkOMKJ/NjkZKazbYYHDGBM/FjiicLw+voEjKUk4Y2ge63Za4DDGxI8FjijEu6kKnH6O8v3H7UFAY0zcWOCIQrybqsAJHADrrJ/DGBMnFjiikAiBY2JhLslJYoHDGBM3FjiikAhNVRmpPsYPybE7q4wxceNp4BCR2SKyVUTKReT2MOtFRBa768tEZLKbXiQir4jIZhHZJCK3Bm3TT0ReEpHt7nuel+cQrDoBahwAZ5yax4ZdR2n0B+JaDmNM7+RZ4BARH/AbYA4wDrhKRMaFZJsDjHJf84GH3HQ/8D1VHQtMA24J2vZ2YKWqjgJWustdoqbBjy9JSE+Jb0WtdFgeDf6AzQhojIkLL6+AU4ByVa1Q1UbgaWBuSJ65wJPqWA3kishgVd2rqu8BqGo1sBkYErTNE+7nJ4CveHgOJ6hpaCYr1YeIdNUhw2rpIF+743Bcy2GM6Z28DBxDgF1By5V8evGPOI+IDAMmAe+4SQNVdS+A+z4g3MFFZL6IrBWRtQcOHDjZczhBdX18RsYNNbBvOoV5Gbxnz3MYY+LAy8AR7md56Oh87eYRkWzgGeCfVDWqdhlVfURVS1W1tKCgIJpN2xSvSZzCKR2ax9odR2zAQ2NMl/MycFQCRUHLhcCeSPOISApO0PiDqj4blGefiAx28wwG9se43G2qafTHbYDDUGcMzWN/dQOVR+riXRRjTC/jZeBYA4wSkeEikgpcCSwPybMcuMa9u2oaUKWqe8XpRPh3YLOq3hdmm2vdz9cCf/buFE6UKE1VAGcM7QfYg4DGmK7nWeBQVT+wEHgRp3N7mapuEpEFIrLAzbYCqADKgSXAzW762cA3gC+IyHr3dZG77m7gAhHZDlzgLneJmgY/WamJEThGD+pDn7Rk1lgHuTGmi3l6FVTVFTjBITjt4aDPCtwSZrs3CN//gaoeAmbFtqSRqWlInKYqX5JQOiyPdz6ywGGM6Vr25HgUjjf4yUqNz+x/4UwdkU/5/uMcqG6Id1GMMb2IBY4o1DU1k5kgfRwA00bkA/DOR4fiXBJjTG/SYeAQkYVdOaxHomr0B2hqVjJTEqfGMeGUvmSl+ninwpqrjDFdJ5IaxyBgjYgsc8eeiu9j03FS19gMOIMMJopkXxKlw/qxusJqHMaYrtNh4FDVf8EZS+rfgeuA7SLycxEZ6XHZEkptkzPAYWaC3FXVYtqIfLbvP87B49bPYYzpGhH1cbh3P33ivvxAHvDfInKPh2VLKLVujSMzgWocAFNHOM9zvGt3VxljukgkfRyLRGQdcA/wJnC6qt4EnAH8L4/LlzASsakK4PQhOWSm+qy5yhjTZSJpd+kPfFVVPw5OVNWAiHzJm2IlnkStcaT4kjhjaJ4FDmNMl4mkj+OO0KARtG5z7IuUmGobW/o4EitwgNPPsW3fcQ5ZP4cxpgvYcxwRqmutcSRW5zjANOvnMMZ0IQscEUrUpiqA04fkkpHi421rrjLGdAELHJEoW0bt3+4EIOP3s6BsWXzLEyI1OYkpw/vxRvnBeBfFGNMLWODoSNkyeH4RdXW1AGRWfwTPL0q44HHuqP5UHKhh91Gbn8MY4y0LHB1ZeRc01VFLGgAZNEBTnZOeQM4d5cxy+Mb22EyTa4wxbbHA0ZGqSgDqNI00GvGJnpCeKE4bmM2APmm8vt2aq4wx3rLA0ZGcQgBqSSOThs+kJwoR4ZxR/Xmj/CDNAZuH3BjjHQscHZl1B6RkUKtpZFLvpKVkOOkJ5rxRBRytbWLTnqp4F8UY04NZ4OhI8Ty4ZDF1KXlkSCPkFMEli530BHP25/oDWHOVMcZTngYOdxj2rSJSLiK3h1kvIrLYXV8mIpOD1j0mIvtFZGPINneKyO4wc5F7p3getcNmkTlkHNy2MSGDBkBBnzTGDu7L69ZBbozxkGeBQ0R8wG+AOcA44CoRGReSbQ7OkO2jgPnAQ0HrHgdmt7H7+1W1xH2taCNPTNU2NpORQJM4teXcUf1Z9/GR1iFSjDEm1ryscUwBylW1QlUbgaeBuSF55gJPqmM1kCsigwFUdRWQMGNo1Db6yUqgaWPbcu6o/jQ1q80KaIzxjJeBYwiwK2i50k2LNk84C92mrcfamtZWROaLyFoRWXvgQOebbmobmxNuSPVwzhzWj9TkJFZZc5UxxiNeBo5wU8yG3icaSZ5QDwEjgRJgL3BvuEyq+oiqlqpqaUFBQQe77FhdY3NCzTfelvQUH9NG5PPaVgscxhhveBk4KoGioOVCYM9J5DmBqu5T1WZVDQBLcJrEPFfb2JyQAxyGM2vMACoO1lBx4Hi8i2KM6YG8DBxrgFEiMlxEUoErgeUheZYD17h3V00DqlR1b3s7bekDcV0KbGwrbyzVNTaTkYBDqofzhTEDAPj7lv1xLokxpifyLHCoqh9YCLwIbAaWqeomEVkgIgvcbCuACqAcp/Zwc8v2IvIU8DYwWkQqReSb7qp7ROQDESkDZgK3eXUOLfzNARqbA92mxlHUL5PRA/uwcrMFDmNM7Hn6E9q9VXZFSNrDQZ8VuKWNba9qI/0bsSxjJGqbEncujrZ8YewAlqyqoKquiZyMlHgXxxjTg9iT4xFomf2vO9xV1WLWmAH4A2oPAxpjYs4CRwTqEnj2v7ZMOjWP3MwU/m7NVcaYGLPAEYE6t6mqOzw53sKXJMwcPYBXtu630XKNMTFlgSMC9W7gSOtGgQNg1tgBHKltYu0Oe4rcGBM7Fjgi0FLjSE/uXoFjxugBpCYn8cLGT+JdFGNMD2KBIwINTQGge3WOA2SnJXPeqAJe3PQJAWuuMsbEiAWOCLQ0VaWndL+va86EQeytqmdD5dF4F8UY00N0vythHHTXpiqA88cOJDlJ+Js1VxljYsQCRwTqu2lTFUBOZgrTP9efFzZ+gvO8pTHGdI4FjgjUd+MaBzjNVTsP1/Lh3mPxLooxpgewwBGButbbcbvn13XhuIEkCbzwgTVXGWM6r3teCbtYQ1MzIpCW3D2/rvzsNKaP7M/yDXusucoY02nd80rYxer9AdKTfYiEm3eqe5hbcgo7D9fy3s6j8S6KMaabs8ARgbrG5m55K26w2RMGkZacxJ/X7453UYwx3Vz3vhp2kfqmZtK72XAjofqkp3D+2IH8tWwvTc2BeBfHGNONWeCIQL0/0O0DBzjNVYdqGnmj/GC8i2KM6cYscESgJ9Q4wBm7KicjhT+/b81VxpiTZ4EjAk7g6P5fVWpyEhedPpgXN+3jeIM/3sUxxnRTnl4NRWS2iGwVkXIRuT3MehGRxe76MhGZHLTuMRHZLyIbQ7bpJyIvich29z3Py3MAN3B004f/Ql1eWkhdUzN/2bAn3kUxxnRTngUOEfEBvwHmAOOAq0RkXEi2OcAo9zUfeCho3ePA7DC7vh1YqaqjgJXusqfqmwLdcriRcCYV5XLawGyeXrMr3kUxxnRTXtY4pgDlqlqhqo3A08DckDxzgSfVsRrIFZHBAKq6Cgg3A9Fc4An38xPAV7wofLC6HtJUBSAiXHHmqazfdZQtn9gQJMaY6Hl5NRwCBP+srXTTos0TaqCq7gVw3weEyyQi80VkrYisPXDgQFQFD9WTmqoALp00hFRfEkut1mGMOQleBo5wj1mHjncRSZ6ToqqPqGqpqpYWFBR0al/1TYFuN21se/plpXLh+IE89/7u1gEcjTEmUl4GjkqgKGi5EAjtkY0kT6h9Lc1Z7vv+TpazQ/VNzWT0oMABcOWZp3K0tokXN9nAh8aY6HgZONYAo0RkuIikAlcCy0PyLAeuce+umgZUtTRDtWM5cK37+Vrgz7EsdDg95XbcYNNH5jMsP5Mn3toR76IYY7oZz66GquoHFgIvApuBZaq6SUQWiMgCN9sKoAIoB5YAN7dsLyJPAW8Do0WkUkS+6a66G7hARLYDF7jLnvE3B/AHtEc8ABgsKUm4dvow3tt5lA27jsa7OMaYbiTZy52r6gqc4BCc9nDQZwVuaWPbq9pIPwTMimEx21Xvd2f/62GBA+CyMwq593+28fhbO7j/ipJ4F8cY0030rPYXD9Q1urP/9bCmKnAGPrzsjEL+UraH/cfq410cY0w30fOuhjFW3zr7X8+rcQBcO30Y/oDyh3d2xrsoxphuwgJHBxr8LTWOnhk4hvfPYuboAfzn6o9ba1fGGNMeCxwdqGvsuX0cLW48bwSHahpZusZqHcaYjlng6EC9v+f2cbSYOiKfM4fl8ciqChr9NsmTMaZ9PfdqGCMtfRw9tamqxS0zP8eeqnr+ZHN1GGM6YIGjA613VfWgsarC+fxpBUwY0peHXvsHzYGYjPpijOmhLHB0oPU5jtSe/VWJCLfM+BwfHazheZurwxjTjp59NYyBhpbbcXt4jQPgi+MHMXZwX+57aZv1dRhj2mSBowMN7gU0rQd3jrdIShK+P3s0Ow/X8rTdYWWMaUPPvxp2Umvg6AU1DoAZpxUwdXg/Fq8sp8bmJTfGhGGBowMtDwCmJfeOr0pE+P7sMRw83sBjb3wU7+IYYxJQ77gadkJDk1PjSPX1nq/qjKF5XDhuIA+99g8+qbIxrIwxJ+o9V8OT1OAPkOpLIikp3GSFPde/XDwOf0D5+YrN8S6KMSbBWODoQIO/udc0UwU7NT+TBZ8fyfINe1hdcSjexTHGJJDed0WMUqM/0CvuqArnps+PZEhuBncu34S/2W7PNcY4eucVMQoN/kCvuaMqVEaqjzsuGceWT6p55PWKeBfHGJMgLHB0wAkcvfdr+uL4QVx0+iAeeGk72/dVx7s4xpgE4OkVUURmi8hWESkXkdvDrBcRWeyuLxORyR1tKyJ3ishuEVnvvi7y8hwamppJ7cWBA+CuuRPISvPxz/9dZk1WxhjvAoeI+IDfAHOAccBVIjIuJNscYJT7mg88FOG296tqiftagYd6e40DoH92GnfNncCGXUf53SprsjKmt/PyijgFKFfVClVtBJ4G5obkmQs8qY7VQK6IDI5w2y7h3FXVO/s4gn2peDAXnz6Y+1/axrqPj8S7OMaYOPIycAwBdgUtV7ppkeTpaNuFbtPWYyKSF+7gIjJfRNaKyNoDBw6c7Dn06ruqgokIP//q6QzOTWfRU+9ztLYx3kUyxsSJl1fEcE/MhU700Fae9rZ9CBgJlAB7gXvDHVxVH1HVUlUtLSgoiKjA4VhT1adyMlJ48KrJ7K+u55//WIaqzdthTG/k5RWxEigKWi4EQid6aCtPm9uq6j5VbVbVALAEp1nLM735dtxwJhblcvucsby8eR+/Xrk93sUxxsSBl4FjDTBKRIaLSCpwJbA8JM9y4Br37qppQJWq7m1vW7cPpMWlwEYPz6HXPjnenhvOHsb/mlzIAy9v5y9lNumTMb1Nslc7VlW/iCwEXgR8wGOquklEFrjrHwZWABcB5UAtcH1727q7vkdESnCarnYAN3p1DuAMctjbb8cN5fR3TODjQzV8b9kGCvMyKSnKjXexjDFdRHpDO3VpaamuXbv2pLad+JP/4Sslp/CTuRNiXKru79DxBr7y2zc5Xu9n2Y1nMWpgn3gXyRgTQyKyTlVLQ9Ptp3QHGvzNpKVYH0c4+dlp/Oc3p5LiS+LqR99h56HaeBfJGNMFLHC0Q1Wd23GtqapNQ/Oz+I9vTqWxOcDXHl1twcOYXsCuiO3wB5SA9p7Z/07W6EF9ePKGKRxv8HP5796yMa2M6eHsitiO3jbfeGcUF+aydP5ZBBTm/e5t1u86Gu8iGWM8YoGjHQ1NznzjdldVZEYP6sMfbzyL7PRkrvjd2yzfYLfqGtMT2RWxHZ/WOOxritSw/ln86eazKS7MYdFT7/OrF7faiLrG9DB2RWxHa+Cwsaqikp+dxn9+ayrzSgt58JVyvrbkHfYcrYt3sYwxMWJXxHY0+J2mKuvjiF5aso97LpvIffMmsmlPFXN+/Tp/27g33sUyxsSABY52NFpTVad9dXIhf1l0Lqf2y2TBf77Hjf+xlr1VVvswpjuzK2I77K6q2BjeP4tnb57OD2aP4bVtBzj/3td49PWK1hqdMaZ7scDRjoYm6+OIlRRfEjfNGMlLt32eM4f346d/3cyse1/jufcrCQR6/rA3xvQkdkVsR8sv4lSffU2xUtQvk99fdyZP3jCFnIwUblu6gdm/XsUz6ypbmwaNMYnNrojtsLuqvCEinHdaAc8vPIfFV00C4Ht/3MB597zCw6/9gwPVDXEuoTGmPZ4Nq94T2F1V3kpKEr488RQuKR7Ma9sO8MiqCu5+YQu/enErM8cM4PIzCpk5ZgApVuMzJqFY4GiH3VXVNUSEGaMHMGP0AMr3V/PHdZU8s243L324j5yMFGaNGcCF4wdy3mkFZKban6wx8Wb/C9thT453vc8N6MMP54zlny8czWtbD/DCxk9YuWUfz76/m7TkJEqH5XHWiHzOGplPcWGu1UaMiQMLHO349K4qa6rqaim+JM4fN5Dzxw3E3xzg3R2HWbl5P2/94xC/+p9tAGSm+phYmMvphTmcPsR5Dc3PRETiXHpjejYLHO34tI/DftXGU7Iviekj+zN9ZH8ADtc08k7FIVZXHGJ9ZRWPv7mDRnc8rKxUHyMHZDOyIJuRBVmMLMhmaH4WQ3Iz6JuRbEHFmBjwNHCIyGzg1zjzhj+qqneHrBd3/UU4c45fp6rvtbetiPQDlgLDcOYcn6eqR2Je+LJlNLy5GmEmyb8uhvPvgOJ5MT+MiV6/rFTmnD6YOacPBpy+qG37qvlgdxVbP6nmHweO807FIZ57f/cJ22Wl+jglN8N9pdM/O438rFT6ZafRPyuVftmp5GelkZeZQrI1gRnTJs8Ch4j4gN8AFwCVwBoRWa6qHwZlmwOMcl9TgYeAqR1sezuwUlXvFpHb3eUfxLTwZcvg+UU01F9KGk3IsV3w/CJnnQWPhJOanMSEITlMGJJzQnpNg5+KAzXsOlLL7iN17D5ax56jdeypqmPj7ioO1zaibTx7mJHiIzs9mT5pyfRJTyY7PZnstGSy01LISvORlpxEesqJ72nJPtJSnPf0lCRSk5NI8SWRnCQkJyXhSxKSfYIvSUhJSsLnE5KTnOWW9xSfk88ngghWQzIJycsaxxSgXFUrAETkaWAuEBw45gJPqqoCq0UkV0QG49Qm2tp2LjDD3f4J4FViHThW3gVNdTSQQhpNTlpTnZNugaPbyEpLdvo/CnPCrm8OKEdrGzlU08ih440crmnkUE0DR2qaON7QxPEGP9X1/tb3g9W1HG9wlhv9Aer9zW0GnlgSAQGSgoLJCcvuZ+TTtCQ3T7j8SR0Eo7ZWt5lO+BVt529r/23sp438ba2Idv893c8vPZ0pw/vFdJ9eBo4hwK6g5UqcWkVHeYZ0sO1AVd0LoKp7RWRAuIOLyHxgPsCpp54aXcmrKgEYKx9T50v9TLrpGXxJQn52GvnZaTAw+u1VlaZmpd7fTENTgAZ/M/Xue4M/QH1TM80BxR9Qmpvd94DiDwTwN+un6wIB/AHF33zisiqoexznszOVsaqbBgQC6uaBgBvFAm7+gH66ToPS2hvhxdki7IpoktE2Imrb+b3df9srer6stNjf3ONl4AgX3kP/+drKE8m27VLVR4BHAEpLS6P7s8kphKpdXJn8Klfy6onpxrhEhNRkcWaITI93aYzpOl72AFYCRUHLhUDoXKJt5Wlv231ucxbu+/4Yltkx6w5IyTgxLSXDSTfGmF7Oy8CxBhglIsNFJBW4Elgekmc5cI04pgFVbjNUe9suB651P18L/DnmJS+eB5cshpwiQJz3SxZb/4YxxuBhU5Wq+kVkIfAizi21j6nqJhFZ4K5/GFiBcytuOc7tuNe3t62767uBZSLyTWAncLknJ1A8zwKFMcaEIW11MvUkpaWlunbt2ngXwxhjuhURWaeqpaHp9pSTMcaYqFjgMMYYExULHMYYY6JigcMYY0xUekXnuIgcAD4+yc37AwdjWJzuwM65d7Bz7h06c85DVbUgNLFXBI7OEJG14e4q6MnsnHsHO+fewYtztqYqY4wxUbHAYYwxJioWODr2SLwLEAd2zr2DnXPvEPNztj4OY4wxUbEahzHGmKhY4DDGGBMVCxwuEZktIltFpNydyzx0vYjIYnd9mYhMjkc5YymCc77aPdcyEXlLRCbGo5yx1NE5B+U7U0SaReSyrixfrEVyviIyQ0TWi8gmEXmtq8sYaxH8XeeIyPMissE95+vjUc5YEpHHRGS/iGxsY31sr1/OlJK9+4UzdPs/gBFAKrABGBeS5yLgBZzZCacB78S73F1wztOBPPfznN5wzkH5/o4z7P9l8S63x//GucCHwKnu8oB4l7sLzvlHwC/dzwXAYSA13mXv5HmfB0wGNraxPqbXL6txOKYA5apaoaqNwNPA3JA8c4En1bEayG2ZibCb6vCcVfUtVT3iLq7GmYmxO4vk3xngO8AzeDG7ZNeK5Hy/BjyrqjsBVLU3nLMCfUREgGycwOHv2mLGlqquwjmPtsT0+mWBwzEE2BW0XOmmRZunO4n2fL6J84ulO+vwnEVkCHAp8HAXlssrkfwbnwbkicirIrJORK7pstJ5I5JzfhAYizMd9QfAraoa6JrixU1Mr1+ezQDYzUiYtND7lCPJ051EfD4iMhMncJzjaYm8F8k5PwD8QFWbnR+k3Vok55sMnAHMAjKAt0Vktapu87pwHonknL8IrAe+AIwEXhKR11X1mMdli6eYXr8scDgqgaKg5UKcXyPR5ulOIjofESkGHgXmqOqhLiqbVyI551LgaTdo9AcuEhG/qv6pS0oYW5H+XR9U1RqgRkRWAROB7ho4Ijnn64G71Wn8LxeRj4AxwLtdU8S4iOn1y5qqHGuAUSIyXERSgSuB5SF5lgPXuHcnTAOqVHVvVxc0hjo8ZxE5FXgW+EY3/gUarMNzVtXhqjpMVYcB/w3c3E2DBkT2d/1n4FwRSRaRTGAqsLmLyxlLkZzzTpwaFiIyEBgNVHRpKbteTK9fVuMAVNUvIguBF3HuynhMVTeJyAJ3/cM4d9hcBJQDtTi/WrqtCM/5DiAf+K37C9yv3Xhk0QjPuceI5HxVdbOI/A0oAwLAo6oa9pbO7iDCf+P/CzwuIh/gNOH8QFW79VDrIvIUMAPoLyKVwI+BFPDm+mVDjhhjjImKNVUZY4yJigUOY4wxUbHAYYwxJioWOIwxxkTFAocxxpioWOAwxhgTFQscxhhjomKBw5g4cOf7KBORdBHJcueFmBDvchkTCXsA0Jg4EZGfAuk4gwtWquov4lwkYyJigcOYOHHHUloD1APTVbU5zkUyJiLWVGVM/PTDmUioD07Nw5huwWocxsSJiCzHmaFuODBYVRfGuUjGRMRGxzUmDtyZ9vyq+l8i4gPeEpEvqOrf4102YzpiNQ5jjDFRsT4OY4wxUbHAYYwxJioWOIwxxkTFAocxxpioWOAwxhgTFQscxhhjomKBwxhjTFT+P9f/NvPIsdBOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# initial conditions\n",
    "x0, v0 = 0, 2 #starting position/velocity\n",
    "\n",
    "# initialize material/geometry parameters\n",
    "d = w0 = 20\n",
    "\n",
    "#analytical solution over the full domain\n",
    "x = torch.linspace(0,1,500).view(-1,1)\n",
    "y = CD_oscillator(d, w0, x, x0, v0).view(-1,1)\n",
    "\n",
    "# slice out a small number of points \n",
    "#that will serve as training data\n",
    "x_data = x[0:25:5]\n",
    "y_data = y[0:25:5]\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(x, y, label=\"Exact solution\")\n",
    "plt.scatter(x_data, y_data, color=\"tab:orange\", label=\"Training data\")\n",
    "plt.title(\"Generated training data\")\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"y\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3c72f7e",
   "metadata": {},
   "source": [
    "#### FCN\n",
    "Fitting the FCN for the critically damped case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f5d5afb0",
   "metadata": {},
   "outputs": [],
   "source": [
    "torch.manual_seed(123)\n",
    "model = FCN(1,1,18,3) # (n_in, n_out, n_hidden, n_layers)\n",
    "optimizer = torch.optim.Adam(model.parameters(),lr=5e-3)\n",
    "epochs = range(1, 801)\n",
    "losses = []\n",
    "files = []\n",
    "\n",
    "for epoch in epochs:\n",
    "    optimizer.zero_grad()\n",
    "    yh = model(x_data)\n",
    "    loss = torch.mean((yh-y_data)**2)# MSE\n",
    "    loss.backward()\n",
    "    optimizer.step()\n",
    "    losses.append(loss)\n",
    "        \n",
    "    # plot results during training\n",
    "    if (epoch) % 10 == 0: # every 10 epochs\n",
    "        \n",
    "        yh = model(x).detach()\n",
    "        \n",
    "        plot_result(epoch,x,y,x_data,y_data,yh)\n",
    "        \n",
    "        file = \"./PINNplots/\"+str(epoch)+\".png\"\n",
    "        plt.savefig(file, bbox_inches='tight', pad_inches=0.1, dpi=100, facecolor=\"white\")\n",
    "        files.append(file)\n",
    "        \n",
    "        #show plot every 500 epochs:\n",
    "        #if (epoch) % 500 == 0: plt.show() \n",
    "        #else: plt.close(\"all\")\n",
    "        plt.close(\"all\")    \n",
    "            \n",
    "save_gif_PIL(\"cd_nn.gif\", files, fps=20, loop=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fa84736",
   "metadata": {},
   "source": [
    "![SegmentLocal](cd_nn.gif \"gif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1af8047",
   "metadata": {},
   "source": [
    "#### PINN\n",
    "Fitting a PINN for the critically damped case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "97e8c1a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# sample 30 uniformly spaced points\n",
    "x_physics = torch.linspace(0,1,30).view(-1,1).requires_grad_(True)\n",
    "\n",
    "mu, k = 2*d, w0**2\n",
    "\n",
    "torch.manual_seed(123)\n",
    "model = FCN(1,1,15,3)\n",
    "epochs = range(1, 20001)\n",
    "losses = []\n",
    "losses_data = []\n",
    "losses_physics = []\n",
    "optimizer = torch.optim.Adam(model.parameters(),lr=5e-4)\n",
    "files = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1bdf5eff",
   "metadata": {},
   "outputs": [],
   "source": [
    "for epoch in epochs:\n",
    "    optimizer.zero_grad()\n",
    "    \n",
    "    # Data loss\n",
    "    yh = model(x_data)\n",
    "    loss_data = torch.mean((yh-y_data)**2) # MSE\n",
    "    \n",
    "    # Physics loss\n",
    "    yh_physics = model(x_physics)\n",
    "    \n",
    "    # calculate first and second derivatives\n",
    "    # dy/dx\n",
    "    dx  = torch.autograd.grad(yh_physics, x_physics, torch.ones_like(yh_physics), create_graph=True)[0]\n",
    "    \n",
    "    # d^2y/dx^2\n",
    "    dx2 = torch.autograd.grad(dx, x_physics, torch.ones_like(dx), create_graph=True)[0]\n",
    "    \n",
    "    # computes the residual of the harmonic oscillator differential equation\n",
    "    physics = dx2 + mu*dx + k*yh_physics # from (1) (beltoforion.de)\n",
    "    loss_physics = (1e-5)*torch.mean(physics**2)\n",
    "    \n",
    "    # backpropagate joint loss\n",
    "    loss = loss_data + loss_physics # add two loss terms together\n",
    "    loss.backward()\n",
    "    optimizer.step()\n",
    "    losses.append(loss)\n",
    "    losses_data.append(loss_data)\n",
    "    losses_physics.append(loss_physics)\n",
    "    \n",
    "    \n",
    "    # plot the result as training progresses\n",
    "    if epoch % 150 == 0: \n",
    "        \n",
    "        yh = model(x).detach()\n",
    "        xp = x_physics.detach()\n",
    "        \n",
    "        plot_result(epoch,x,y,x_data,y_data,yh,xp)\n",
    "        \n",
    "        file = \"./PINNplots/\"+str(epoch)+\".png\"\n",
    "        plt.savefig(file, bbox_inches='tight', pad_inches=0.1, dpi=100, facecolor=\"white\")\n",
    "        files.append(file)\n",
    "        \n",
    "        #if epoch % 6000 == 0: plt.show()\n",
    "        #else: plt.close(\"all\")\n",
    "        plt.close(\"all\")\n",
    "            \n",
    "save_gif_PIL(\"cd_pinn.gif\", files, fps=20, loop=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e719df8d",
   "metadata": {},
   "source": [
    "![SegmentLocal](cd_pinn.gif \"gif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3498b0f1",
   "metadata": {},
   "source": [
    "### Overdamped case\n",
    "\n",
    "Finally, in the overdamped case we have $ \\delta > \\omega_0 $. The analytical solution for the overdamped case DE is:\n",
    "\n",
    "$ x(t) = e^{-\\delta t} (C_1 \\ e^{\\alpha t} \\ + \\ C_2 \\ e^{-\\alpha t}) $, with $ \\alpha = \\sqrt{\\delta^2 - \\omega_{0}^2} $\n",
    "\n",
    "where $C_{1,2}$ are integration constants that can be determined from the given initial conditions (see [here](https://beltoforion.de/en/harmonic_oscillator/) for a detailed solution of the DE).\n",
    "\n",
    "For this case, we consider the material/geometry parameters $\\delta = 20$ and $\\omega_0=1$, and try to learn the solution for the range $x \\in [0,1]$. In addition, we assume the initial conditions $x_0 = 0$ and $v_0 = 2$ in order to calculate c1 and c2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "7b311a91",
   "metadata": {},
   "outputs": [],
   "source": [
    "def OD_oscillator(d, w0, x, x0, v0):\n",
    "    \"\"\"Defines the analytical solution to the 1D overdamped harmonic oscillator problem. \n",
    "    Equations from: https://beltoforion.de/en/harmonic_oscillator/\"\"\"\n",
    "    assert d > w0\n",
    "    \n",
    "    a = np.sqrt(d**2-w0**2)    # a\n",
    "    c1 = (x0*(a+d)+v0)/(2*a)   # c1,2 computed using the predefined\n",
    "    c2 = -c1                   # initial conditions x0 and v0\n",
    "    exp1 = torch.exp(a*x)      # exp(a*t)\n",
    "    exp2 = torch.exp(-a*x)     # exp(-a*t)\n",
    "    summ = (c1*exp1)+(c2*exp2) # c1*exp(a*t)+c2*exp(-a*t)\n",
    "    exp3 = torch.exp(-d*x)     # exp(-d*t) \n",
    "    y = exp3*summ              # y - exact solution\n",
    "    \n",
    "    return y\n",
    "\n",
    "# plot_result redefined with different x,y limits and text position\n",
    "def plot_result(epoch,x,y,x_data,y_data,yh,xp=None):\n",
    "    \"Pretty plot for training results\"\n",
    "    plt.figure(figsize=(8,4))\n",
    "    plt.plot(x,y, color=\"grey\", linewidth=2, alpha=0.8, label=\"Exact solution\")\n",
    "    plt.plot(x,yh, color=\"tab:blue\", linewidth=4, alpha=0.8, label=\"Neural network prediction\")\n",
    "    plt.scatter(x_data, y_data, s=60, color=\"tab:orange\", alpha=0.4, label='Training data')\n",
    "    if xp is not None:\n",
    "        plt.scatter(xp, -0*torch.ones_like(xp), s=60, color=\"tab:green\", alpha=0.4, \n",
    "                    label='Physics loss training locations')\n",
    "    l = plt.legend(loc=(1.01,0.34), frameon=False, fontsize=\"large\")\n",
    "    plt.setp(l.get_texts(), color=\"k\")\n",
    "    plt.xlim(-0.05, 1.05)\n",
    "    plt.ylim(-0.05, 0.05)\n",
    "    plt.text(1.065, 0.1,\"Training step: %i\"%(epoch),fontsize=\"xx-large\",color=\"k\")\n",
    "    plt.axis(\"off\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7ecc601",
   "metadata": {},
   "source": [
    "Let's start by generating the data - Note that the critically damped and overdamped cases can be easily fitted with non-physics-informed NNs, therefore for this example we select a particularly limited range to be represented with the training data, so that we see the effect of the physics loss in action."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "b225e28e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAApLUlEQVR4nO3de3xU9Z3/8dcnk0AS7jcVCRRwEQQXgUaw2iJouwW18qsPr9W24tZbpVDbbrXd31L1V7du14e6PGxltWv789euSktt0aXrekOrVgVbiSJiKdKSghhuSciFXObz++OchMkwSSZkTsZk3s/HYx6Zc873nPmcCcw73+85c465OyIiIsnysl2AiIh8OCkgREQkJQWEiIikpIAQEZGUFBAiIpKSAkJERFJSQIgcBTO70sxejPg1NpnZvEy37Q4zG29mbmb5Ub+WZJ8CQjLGzC41s1fNrMbMPgiff9nMLNu1JTOzdWb2pYi2nZEPUXef5u7rMt22p5jZPDMrz3YdcvQUEJIRZvZ14N+AfwWOA44FrgPOAPr1cC0f+r9ue0ONIgoI6TYzGwLcBnzZ3X/h7tUe+IO7X+7uh8J2/c3sTjP7i5ntNrOVZlYULptnZuVm9vWw97HLzBYnvEY6695kZu8DPzazYWb2hJlVmNn+8HlJ2P524BPAvWZ20MzuDedPMbOnzGyfmW0xs4sTXn+Ema0xsyozew04oYO35IXw54Fw+x8Lh6ReMrO7zWwfcIuZnWBmz5rZXjPbY2Y/M7OhCa+53cw+GT6/xcxWmdlDZlYdDimVHmXbWWb2h3DZz83sUTP7bju/21j4vu8xs23AuUnLF5vZ5nBb28zs2nD+AOA3wPHhe3DQzI43s9lm9jszOxD+ju81sx79A0LSp4CQTPgY0B/4dSft/gU4EZgB/A0wBliesPw4YEg4/++BH5jZsC6sOxz4CHANwb/tH4fT44A64F4Ad/9H4LfAEncf6O5Lwg+0p4D/BI4BLgN+aGbTwu3/AKgHRgNXhY/2zA1/Dg23/7tweg6wLdz+7YAB3wOOB04CxgK3dLDd84FHgKHAmpb96Urb8MP4MeAnBO/Xw8BnO9jO1cB5wEygFLgwafkH4fLBwGLgbjOb5e41wEJgZ/geDHT3nUAzcCMwkuDfzdnAlzt4fckmd9dDj249gCuA95PmvQwcIPhgnkvwYVgDnJDQ5mPAe+HzeWHb/ITlHwCnpbluA1DYQY0zgP0J0+uALyVMXwL8Nmmdfwe+A8SARmBKwrJ/Bl5s57XGA560L1cCf+nkffxfwB8SprcDnwyf3wI8nbBsKlDX1bbh7+KvgCUsfxH4bjs1PQtclzD9d8n7ltT+V8CyhN9LeSf7/FXgsWz/G9Yj9UPjoJIJe4GRZpbv7k0A7n46QHiQMg8YBRQDryccszaCD9/W7bSsH6oFBqa5boW717cuNCsG7gYWAC29kEFmFnP35hT78BFgjpkdSJiXD/y/8PXzgR0Jy/6c8p3oWOL6mNkxwAqC4a5BBO/T/g7Wfz/heS1QmPiep9OWoLfyVw8/nVPVleR4OthvM1tIEKInhvUXA2+2tzEzOxG4i6A3Ukzwvr7ewetLFmmISTLhd8AhYFEHbfYQ9BCmufvQ8DHE3Qemsf101k2+LPHXgcnAHHcfzOFhH2un/Q7g+YTtD/VgWOR6oAJoIhgCajGug3rbu0Ry8vzvhfOmhzVekVBfVHYBY5LOLBvbXuOwfcr9NrP+wGrgTuBYdx8KrKX99xjgPuAdYFK4z98m+n2Wo6SAkG5z9wPArQRj9hea2UAzyzOzGcCAsE0ceIBgjPoYADMbY2afTmP7R7PuIIJQOWBmwwn+yk20G5iYMP0EcKKZfd7MCsLHqWZ2Utjj+CXBgeViM5sKfLGD164A4knbb6/Gg2GNY4B/6KR9JvyO4DjAEjPLN7NFwOwO2q8ClppZSXg86OaEZf0Ijj1VAE1hb+LvEpbvBkZYcBJDi0FAFXDQzKYA13d7jyQyCgjJCHf/PvA14JsExw52E4zh30RwPILw+VbgFTOrAp4m+Cs/HV1d9x6giKD38Qrw30nL/w24MDzDaYW7VxN8uF0K7CQYovkXgg9AgCUEw13vExzg/XF7L+zutQQHoV8Kz9Y5rZ2mtwKzgErgvwhCKFLu3gBcQHASwAGCXssTBD3AVB4AngQ2Ar9PrDF8z5YShMh+4HMEB8Rblr9DcBB8W/g+HA98I2xXHW770cztnWSatR2KFJFcY2avAivdvd3Qk9ykHoRIjjGzM83suHCI6YvAdI7sYYnoLCaRHDSZYFhoIPAn4EJ335XdkuTDSENMIiKSkoaYREQkpT41xDRy5EgfP358tssQEek1Xn/99T3uPirVsj4VEOPHj2fDhg3ZLkNEpNcws3avCqAhJhERSUkBISIiKSkgREQkJQWEiIikFGlAmNmC8M5cW83s5hTLzcxWhMvLzGxWwrLtZvammb1hZjryLCLSwyI7i8nMYgR34foUUA6sN7M17v52QrOFwKTwMYfgUsBzEpbPd/c9UdUoIiLti7IHMRvY6u7bwitIPsKR9wtYBDzkgVeAoWY2OsKaREQkTVF+D2IMbe9EVU7b3kF7bcYQ3KTEgf8xMwf+3d3vT/UiZnYNwT2IGTeuo3u49H71jc0cqG3kQF0DB2obqaprpL4pTkPro5mG5jiNzU483vklVPLyjDwzYnkQy8sjZhDLs+B5y7w8yDMjP/F5zFrn5eVBLJwXy8sjP8+I5VnCzzzyY22nY22mjbb3rhGRD4soAyLV//rkT62O2pzh7jvDG8Q8ZWbvuPsLRzQOguN+gNLS0l5/YSl3p3x/HRvLD/DWX6v4y74aduyrY8f+Wg7UNma7vEjkGeTHksMlYToW/CzIy2sz3SZ08oyCWNvp5PXzW7YZS/EaCT/b1JIi+I6sNeE1E2tNCsLk2vLyFIzy4RZlQJTT9laFJQQ3Ykmrjbu3/PzAzB4jGLI6IiD6gtqGJp7fUsHTmz/g+Xc/YM/BBgAKYsbYYcWUDC9meskQjh9axNDiAoYW9WNocQGDCwso6pdHv1iMfvl5rY+CmBFL46/yuENz3Gl2D34mPjzohTSF03F3mprDn+3MO9w+3tqmsbntdFNz+DNs15g03dou7jQ3Bz+bWuY3p25X19hMU3PSuuF6h7fhNDXHk5Zn9++JPKNtmMWO7HmlDMyE8Mo/IkiTwyxFYLZsL2H7weu3H3xHhG9SSBfEkl/jyP1qqVXB2HtEGRDrgUlmNgH4K8Gduj6X1GYNwa0PHyEYfqp0911mNgDIc/fq8PnfAbdFWGtW/HlvDf/x4ns89vu/Un2oiSFFBcybPIrS8cOZUTKUyccNol++zkSOirsTd4IgaQ0zb51uG2ZO4xEBE08ZZqmCMDGUUoZZQpA2tROEycFa3xinKd4ctGs+HO6d1dvYnN1gNCN1kKQImuTgy89LHYSJwVeQHKxH9DhT9QKP7LG2hFpycB8Z0keunzqUe99wamQB4e5NZraE4HaFMeBBd99kZteFy1cS3OD8HIJbSdYCi8PVjwUeC9/MfOA/3b3P3NBkX00Dd/7PFlat30GeGedOH81FpSXMHj+c/JgCoaeYWXjcJZbtUnpcPO40pgi45OALwiqeIszCnl076x3ebry1B5cchB2HmbcJvsaE6UOtwdg2SJuSAz5pvxrjcbJ9d4NYWqHWfk8tMXgSg2hocQG3LTo54/VGerE+d19LEAKJ81YmPHfghhTrbQNOibK2bPnvt3bxj4+9RWVdI5+bM44l8/+GYwYXZrssyTF5eUb/HA3GVD2q1qHT5sNDmikDJ0V4dRpyRwRwGr3FLg6nDikqiOT96lNXc/0wi8ed7z+5hZXP/4mTxwzmZ1fPYcpxg7NdlkhOycsz+rUeA8m9gOwqBUQPaGqOs+yRN/ivN3dx+Zxx3HL+NAo0lCQiH3IKiIjF485Nq9/kv97cxc0Lp3Dt3Im97kCViOQmBUTEvv/kFlb/vpwbP3ki1515QrbLERFJm8Y5IvTkpvdZ+fyfuHzOOJae/TfZLkdEpEsUEBH5oKqef/j5RqaXDGH5Z6ZqWElEeh0FRERueXwT9U1x7rlkBv3zdbaEiPQ+CogIPP9uBWvffJ9lZ09i4qiB2S5HROSoKCAyrDnufG/tZsYNL+bqT0zMdjkiIkdNAZFhj2/cyTvvV/ONT0/WdZREpFfTJ1gGuTv3rfsTk48dxHl/q/seiUjvpoDIoHVbKtiyu5prz5yoSxqLSK+ngMigB196j9FDCvnMKcdnuxQRkW5TQGTIjn21/PaPe7j01HG6zpKI9An6JMuQR9fvIM/gotKSbJciIpIRCogMiMed1b8vZ+6Jozh+aFG2yxERyQgFRAb8/i/72VVZz2dnjsl2KSIiGaOAyIAnynbRLz+Ps086NtuliIhkjAKim+Jx5zdv7WL+5FEM7K+rp4tI36GA6KaN5QfYXXWIhSfri3Ei0rcoILrp+XcrMIO5J47KdikiIhmlgOimdVsqOKVkKMMH9Mt2KSIiGaWA6IZ9NQ1sLD/AvMnqPYhI36OA6IaXtu7BHc7U8JKI9EEKiG549b29DOyfz9+OGZLtUkREMk4B0Q2vvbePj35kGPm69pKI9EH6ZDtK+2oaeHf3QWZPGJ7tUkREIqGAOErrt+8DYI4CQkT6KAXEUdqwfR/98vP42xIdfxCRvkkBcZQ27qhk2vGD6Z8fy3YpIiKRUEAchea489bOSk4pGZrtUkREIhNpQJjZAjPbYmZbzezmFMvNzFaEy8vMbFbS8piZ/cHMnoiyzq7a+sFBahuama7hJRHpwyILCDOLAT8AFgJTgcvMbGpSs4XApPBxDXBf0vJlwOaoajxaG8sPACggRKRPi7IHMRvY6u7b3L0BeARYlNRmEfCQB14BhprZaAAzKwHOBX4UYY1Hpaz8AAP75zNx5MBslyIiEpkoA2IMsCNhujycl26be4BvAvGOXsTMrjGzDWa2oaKiolsFp+vNv1Yx7fjB5OVZj7yeiEg2RBkQqT49PZ02ZnYe8IG7v97Zi7j7/e5e6u6lo0ZFf02keNz54+5qTho9OPLXEhHJpigDohwYmzBdAuxMs80ZwPlmtp1gaOosM/tpdKWmb8f+Wmobmply3KBslyIiEqkoA2I9MMnMJphZP+BSYE1SmzXAF8KzmU4DKt19l7t/y91L3H18uN6z7n5FhLWm7Z33qwGYrIAQkT4uspsou3uTmS0BngRiwIPuvsnMrguXrwTWAucAW4FaYHFU9WTKljAgTjxWASEifVtkAQHg7msJQiBx3sqE5w7c0Mk21gHrIijvqGx5v5pxw4sZ0D/St05EJOv0Teoueuf9Kg0viUhOUEB0QUNTnO17a5ms4SURyQEKiC74y74amuPOCccMyHYpIiKRU0B0wbaKGgAm6BvUIpIDFBBd8N6eloBQD0JE+j4FRBdsq6hh5MB+DCkqyHYpIiKRU0B0wXt7atR7EJGcoYDogm0KCBHJIQqINFXVN7Ln4CEmjtIBahHJDQqING0PD1CPH6EehIjkBgVEmnbsqwNg3PDiLFciItIzFBBp2rG/FoCS4UVZrkREpGcoINK0Y18tQ4oKGFyoU1xFJDcoINJUvr+Oseo9iEgOUUCkacf+WsYO0/EHEckdCog0xONO+f46SoapByEiuUMBkYY9Bw/R0BRnrM5gEpEcooBIQ8sZTBpiEpFcooBIQ8t3IDTEJCK5RAGRhr8eCAJijAJCRHKIAiINu6vqGVyYT3G//GyXIiLSYxQQadhVWc/oIeo9iEhuUUCk4f3Keo4bUpjtMkREepQCIg1BD0IBISK5RQHRiYamOHtrDqkHISI5RwHRiQ+q63FHPQgRyTkKiE68X1kPwLGDFRAiklsUEJ3YFQaEzmISkVyjgOjE7qogIHQMQkRyjQKiE7sq6ynuF2Nwob4kJyK5RQHRifer6jl2cCFmlu1SRER6VKQBYWYLzGyLmW01s5tTLDczWxEuLzOzWeH8QjN7zcw2mtkmM7s1yjo7UlF9iFGD+mfr5UVEsiaygDCzGPADYCEwFbjMzKYmNVsITAof1wD3hfMPAWe5+ynADGCBmZ0WVa0d2XPwEKMGKiBEJPdE2YOYDWx1923u3gA8AixKarMIeMgDrwBDzWx0OH0wbFMQPjzCWtulHoSI5KooA2IMsCNhujycl1YbM4uZ2RvAB8BT7v5qqhcxs2vMbIOZbaioqMhU7QDUNzZTXd+kgBCRnBRlQKQ6qpvcC2i3jbs3u/sMoASYbWYnp3oRd7/f3UvdvXTUqFHdqfcIFdWHADTEJCI5KcqAKAfGJkyXADu72sbdDwDrgAUZr7ATew4GATFyUL+efmkRkayLMiDWA5PMbIKZ9QMuBdYktVkDfCE8m+k0oNLdd5nZKDMbCmBmRcAngXcirDWlwz0IfUlORHJPZN/+cvcmM1sCPAnEgAfdfZOZXRcuXwmsBc4BtgK1wOJw9dHA/w3PhMoDVrn7E1HV2p6KsAehYxAikosi/Xqwu68lCIHEeSsTnjtwQ4r1yoCZUdaWjj3VDQCMGKghJhHJPfomdQcqDtYzrLiAgpjeJhHJPfrk64C+AyEiuUwB0YGK6kOM1CmuIpKjFBAd2HOwQQEhIjlLAdGB/TUNDB+gA9QikpsUEO1oaIpTfaiJEQoIEclRCoh2HKgNTnEdpoAQkRylgGjH3pogIDTEJCK5SgHRjv1hQAwrVkCISG5SQLRjX616ECKS2xQQ7WjtQQwoyHIlIiLZoYBox76aRkBDTCKSuxQQ7dhf28Dgwnxdh0lEclann35mtsTMhvVEMR8m+/QlORHJcen8eXwcsN7MVpnZAjNLdZvQPmdfTYO+AyEiOa3TgHD3/w1MAv4DuBL4o5n9s5mdEHFtWbWvpoHhOv4gIjksrQH28MY+74ePJmAY8Asz+36EtWXV/lr1IEQkt3V6RzkzWwp8EdgD/Aj4B3dvNLM84I/AN6Mtsee5O/tqGnQdJhHJaenccnQkcIG7/zlxprvHzey8aMrKrrrGZg41xdWDEJGc1mlAuPvyDpZtzmw5Hw77a4PvQAwt0pfkRCR36ST/FCrDgBiigBCRHKaASKGyLgyIYgWEiOQuBUQKlXXBdZjUgxCRXKaASKG1B6GAEJEcpoBIQQEhIqKASKmyrpFYnjGwfzpnAYuI9E0KiBQq6xoZXJhPjlx2SkQkJQVECgdqGxmq6zCJSI5TQKRQWdfIYB1/EJEcp4BIoaquUQeoRSTnRRoQ4f0jtpjZVjO7OcVyM7MV4fIyM5sVzh9rZs+Z2WYz22Rmy6KsM1mlAkJEJLqAMLMY8ANgITAVuMzMpiY1W0hwr4lJwDXAfeH8JuDr7n4ScBpwQ4p1I3OgrpEhRTqDSURyW5Q9iNnAVnff5u4NwCPAoqQ2i4CHPPAKMNTMRrv7Lnf/PYC7VwObgTER1toqHneq6hoZWqSD1CKS26IMiDHAjoTpco78kO+0jZmNB2YCr2a+xCMdbGgi7vqSnIhIlAGR6ksE3pU2ZjYQWA181d2rUr6I2TVmtsHMNlRUVBx1sS10JVcRkUCUAVEOjE2YLgF2ptvGzAoIwuFn7v7L9l7E3e9391J3Lx01alS3i265zIZOcxWRXBdlQKwHJpnZBDPrB1wKrElqswb4Qng202lApbvvsuArzP8BbHb3uyKs8Qi6DpOISCCyU3XcvcnMlgBPAjHgQXffZGbXhctXAmuBc4CtQC2wOFz9DODzwJtm9kY479vuvjaqeltUtfYgdBaTiOS2SD8Fww/0tUnzViY8d+CGFOu9SOrjE5GrPtQEwOBC9SBEJLfpm9RJquuDgBhUqB6EiOQ2BUSS6vpgiEmX+haRXKeASHKwvomighj5Mb01IpLb9CmYpLq+ScNLIiIoII5w8JACQkQEFBBtla2iasvzDNyzEe4+GcpWZbsiEZGsUUC0KFsFjy+lutEYbLVQuQMeX6qQEJGcpYBo8cxt0FjHQYoYRF0wr7EumC8ikoMUEC0qywGo9mIGWe0R80VEco0CosWQEgCqKWZgSw8iYb6ISK5RQLQ4eznN+cXUUni4B1FQBGcvz25dIiJZovM5W0y/mIMNwC8IjkEMGRuEw/SLs12ZiEhWKCASVJ9wHvAcgxb9K5w6ttP2IiJ9mYaYEuhCfSIihykgEhwOCF3qW0REAZHg4KHwSq7qQYiIKCASaYhJROQwBUSCKgWEiEgrBUSCgy0B0V/HIEREFBAJDh5qJD/PKCzQ2yIiok/CBDWHmhnQPx8zy3YpIiJZp4BIUHOoiQH9YtkuQ0TkQ0EBkaCmoYni/jpALSICCog2WoaYREREAdGGhphERA5TQCSoaVAPQkSkhQIigXoQIiKHKSAS1DY0qQchIhJSQCQ4eEgBISLSQgERao479Y1xBvRTQIiIgAKiVU1DcB2mAf11DEJEBCIOCDNbYGZbzGyrmd2cYrmZ2YpweZmZzUpY9qCZfWBmb0VZY4vaQ80AGmISEQlFFhBmFgN+ACwEpgKXmdnUpGYLgUnh4xrgvoRlPwEWRFVfspYeRLHOYhIRAaLtQcwGtrr7NndvAB4BFiW1WQQ85IFXgKFmNhrA3V8A9kVYXxs1h4KAGKgehIgIEG1AjAF2JEyXh/O62qZDZnaNmW0wsw0VFRVHVSgEl9kAKNZBahERINqASHXNbD+KNh1y9/vdvdTdS0eNGtWVVdto6UHoILWISCDKgCgHxiZMlwA7j6JNjzh8FpN6ECIiEG1ArAcmmdkEM+sHXAqsSWqzBvhCeDbTaUClu++KsKZ2tQwx6XsQIiKByALC3ZuAJcCTwGZglbtvMrPrzOy6sNlaYBuwFXgA+HLL+mb2MPA7YLKZlZvZ30dVKwSX2QANMYmItIj0z2V3X0sQAonzViY8d+CGdta9LMrakh081HKaq3oQIiKgb1K3qm1opqggRixP96MWEQEFRKvgQn0aXhIRaaGACNXqSq4iIm0oIEI1Dc06/iAikkABEdLd5ERE2lJAhGobmilSQIiItNKYSqi+sZljB/fPdhkivVJjYyPl5eXU19dnuxRpR2FhISUlJRQUFKS9jgIi1HKaq4h0XXl5OYMGDWL8+PGY6VTxDxt3Z+/evZSXlzNhwoS019MQUygYYlJeihyN+vp6RowYoXD4kDIzRowY0eUengIiVN/YrJsFiXSDwuHD7Wh+PwoIgu5XbUOThphERBIoIIBDTXHijs5iEunFYrEYM2bMaH3ccccdGdv2G2+8wdq1aztvmIaBAwd22uaee+6htra2dfqcc87hwIEDGXn9rtCgO8HwEuh+1CK9WVFREW+88UYk237jjTfYsGED55xzTiTbT3bPPfdwxRVXUFxcDJCxcOoqBQTBAWpAQ0wiGXDr45t4e2dVRrc59fjBfOcz07q8XmVlJbNnz2bNmjVMnjyZyy67jLPOOourr76a66+/nvXr11NXV8eFF17IrbfeCsD69etZtmwZNTU19O/fn6eeeorly5dTV1fHiy++yLe+9S0uueSS1tfYtGkTixcvpqGhgXg8zurVq5k0aRJ33XUXDz74IABf+tKX+OpXv9qmtnXr1nHnnXfyxBNPALBkyRJKS0upqqpi586dzJ8/n5EjR/Lcc88xfvx4NmzYwMiRI1Nud/v27SxcuJCPf/zjvPzyy4wZM4Zf//rXFBUVHc3b3UoBAdSFPQgNMYn0XnV1dcyYMaN1uuWD/N577+XKK69k2bJl7N+/n6uvvhqA22+/neHDh9Pc3MzZZ59NWVkZU6ZM4ZJLLuHRRx/l1FNPpaqqiuLiYm677TY2bNjAvffee8Trrly5kmXLlnH55ZfT0NBAc3Mzr7/+Oj/+8Y959dVXcXfmzJnDmWeeycyZMzvdj6VLl3LXXXfx3HPPMXLkyDbL2tvusGHD+OMf/8jDDz/MAw88wMUXX8zq1au54ooruvWeKiCAOvUgRDLmaP7Sz4T2hpg+9alP8fOf/5wbbriBjRs3ts5ftWoV999/P01NTezatYu3334bM2P06NGceuqpAAwePLjT1/3Yxz7G7bffTnl5ORdccAGTJk3ixRdf5LOf/SwDBgwA4IILLuC3v/1tWgHRkfa2e/755zNhwoTWgPzoRz/K9u3bu/VaoIPUwOEhJl2sT6TvicfjbN68maKiIvbt2wfAe++9x5133skzzzxDWVkZ5557LvX19bh7l08H/dznPseaNWsoKiri05/+NM8++yzBvdA6lp+fTzweb51O5zsKHW23f//DV4KIxWI0NTV1ur3OKCDQEJNIX3b33Xdz0kkn8fDDD3PVVVfR2NhIVVUVAwYMYMiQIezevZvf/OY3AEyZMoWdO3eyfv16AKqrq2lqamLQoEFUV1en3P62bduYOHEiS5cu5fzzz6esrIy5c+fyq1/9itraWmpqanjsscf4xCc+0Wa9j3zkI7z99tscOnSIyspKnnnmmdZl7b1eOtvNJP3JDNSF96PWEJNI75V8DGLBggVcddVV/OhHP+K1115j0KBBzJ07l+9+97vceuutzJw5k2nTpjFx4kTOOOMMAPr168ejjz7KV77yFerq6igqKuLpp59m/vz53HHHHcyYMeOIg9SPPvooP/3pTykoKOC4445j+fLlDB8+nCuvvJLZs2cDwcHk5OGlsWPHcvHFFzN9+nQmTZrUZvk111zDwoULGT16NM8991zr/FmzZqXcbiaGk1KxdLpCvUVpaalv2LChy+v98vflfG3VRtZ9Yx7jRw6IoDKRvm3z5s2cdNJJ2S5DOpHq92Rmr7t7aar2GmLi8BCTvgchInKYAoLDZzEVKiBERFopINAX5UREUlFAEAwxFcSMgpjeDhGRFvpEJBhiUu9BRKQtBQRQ29CkL8mJiCRRQAB1jXF9SU6kF9u7d2/rZb6PO+44xowZ0zrd0NDQ4bobNmxg6dKlnb7G6aefnqly25g3bx6dnZ6ffPnvnqI/mwm+KKchJpEeVLYKnrkNKsthSAmcvRymX3zUmxsxYkTrdZhuueUWBg4cyDe+8Y3W5U1NTeTnp/64Ky0tpbQ05dcA2nj55ZePur7uSr78d09RD4KW+1ErIER6RNkqeHwpVO4APPj5+NJgfgZdeeWVfO1rX2P+/PncdNNNvPbaa5x++unMnDmT008/nS1btgDBZbfPO+88IAiXq666innz5jFx4kRWrFjRur2WG/2sW7eOefPmceGFFzJlyhQuv/zy1mskrV27lilTpvDxj3+cpUuXtm43UV1dHZdeeinTp0/nkksuoa6urnXZ9ddfT2lpKdOmTeM73/kOACtWrGi9/Pf8+fPbbRcF9SAIzmIa2F9vhUiPeOY2aKxrO6+xLpjfjV5EKu+++y5PP/00sViMqqoqXnjhBfLz83n66af59re/zerVq49Y55133uG5556jurqayZMnc/3111NQUNCmzR/+8Ac2bdrE8ccfzxlnnMFLL71EaWkp1157LS+88AITJkzgsssuS1nTfffdR3FxMWVlZZSVlTFr1qzWZakuQZ7q8t+p2k2fPj2D71xAPQh0FpNIj6os79r8brjooouIxYL/25WVlVx00UWcfPLJ3HjjjWzatCnlOueeey79+/dn5MiRHHPMMezevfuINrNnz6akpIS8vDxmzJjB9u3beeedd5g4cSITJkwAaDcgXnjhhdb7NEyfPr3NB/uqVauYNWsWM2fOZNOmTbz99tspt5Fuu+6KNCDMbIGZbTGzrWZ2c4rlZmYrwuVlZjYr3XUzpmwVdRXvUfTOL+HukzPezRWRJENKuja/G1rumwDwT//0T8yfP5+33nqLxx9/vN3La6dz2exUbbpyXbtUlxRv7xLkR9suEyILCDOLAT8AFgJTgcvMbGpSs4XApPBxDXBfF9btvnAstLY5RrEdimwsVEQSnL0cCpJuhVlQFMyPUGVlJWPGjAHgJz/5Sca3P2XKFLZt29Z6ZdVHH300Zbu5c+fys5/9DIC33nqLsrIygHYvQQ5tL//dUbtMi3LgfTaw1d23AZjZI8AiILEvtAh4yIPofcXMhprZaGB8Gut2XzgWWkd/iggTOKKxUBEJtfzfyuBZTOn45je/yRe/+EXuuusuzjrrrIxvv6ioiB/+8IcsWLCAkSNHtl6SO9n111/P4sWLmT59OjNmzGhtd8opp6S8BDkcefnv9tplWmSX+zazC4EF7v6lcPrzwBx3X5LQ5gngDnd/MZx+BriJICA6XDdhG9cQ9D4YN27cR//85z+nX+QtQwHnqw1f5szYRj4be6llq3DLgS7tr0gu0+W+AwcPHmTgwIG4OzfccAOTJk3ixhtvzHZZrT5Ml/tOdd++5DRqr0066wYz3e9391J3Lx01alTXKgzHPO/p98OEcCCSsVAR6fseeOABZsyYwbRp06isrOTaa6/NdkndEuUQUzkwNmG6BNiZZpt+aazbfWcvD445JJ5y1wNjoSLSN914440fqh5Dd0XZg1gPTDKzCWbWD7gUWJPUZg3whfBsptOASnfflea63Tf9YvjMChgyFrDg52dW6PiDyFHoS3en7IuO5vcTWQ/C3ZvMbAnwJBADHnT3TWZ2Xbh8JbAWOAfYCtQCiztaN5JCp1+sQBDppsLCQvbu3cuIESNSnsIp2eXu7N27l8LCwi6tp3tSi0i3NTY2Ul5eHtn5+NJ9hYWFlJSUHPGt8I4OUuv6EiLSbQUFBa3fIJa+Q5faEBGRlBQQIiKSkgJCRERS6lMHqc2sAujCV6nbGAnsyWA5vYH2ue/Ltf0F7XNXfcTdU37LuE8FRHeY2Yb2juT3Vdrnvi/X9he0z5mkISYREUlJASEiIikpIA67P9sFZIH2ue/Ltf0F7XPG6BiEiIikpB6EiIikpIAQEZGUciogzGyBmW0xs61mdnOK5WZmK8LlZWY2Kxt1ZlIa+3x5uK9lZvaymZ2SjTozqbN9Tmh3qpk1h3c/7NXS2Wczm2dmb5jZJjN7vqdrzLQ0/m0PMbPHzWxjuM+Ls1FnppjZg2b2gZm91c7yzH9+uXtOPAguG/4nYCLBDYk2AlOT2pwD/IbgjnanAa9mu+4e2OfTgWHh84W5sM8J7Z4luOT8hdmuuwd+z0MJ7uk+Lpw+Jtt198A+fxv4l/D5KGAf0C/btXdjn+cCs4C32lme8c+vXOpBzAa2uvs2d28AHgEWJbVZBDzkgVeAoWY2uqcLzaBO99ndX3b3/eHkKwR37+vN0vk9A3wFWA180JPFRSSdff4c8Et3/wuAu/f2/U5nnx0YZMENKgYSBERTz5aZOe7+AsE+tCfjn1+5FBBjgB0J0+XhvK626U26uj9/T/AXSG/W6T6b2Rjgs8DKHqwrSun8nk8EhpnZOjN73cy+0GPVRSOdfb4XOIngdsVvAsvcPd4z5WVFxj+/cul+EKluc5V8jm86bXqTtPfHzOYTBMTHI60oeuns8z3ATe7e3EfufpbOPucDHwXOBoqA35nZK+7+btTFRSSdff408AZwFnAC8JSZ/dbdqyKuLVsy/vmVSwFRDoxNmC4h+Muiq216k7T2x8ymAz8CFrr73h6qLSrp7HMp8EgYDiOBc8ysyd1/1SMVZl66/7b3uHsNUGNmLwCnAL01INLZ58XAHR4M0G81s/eAKcBrPVNij8v451cuDTGtByaZ2QQz6wdcCqxJarMG+EJ4NsBpQKW77+rpQjOo0302s3HAL4HP9+K/JhN1us/uPsHdx7v7eOAXwJd7cThAev+2fw18wszyzawYmANs7uE6Mymdff4LQY8JMzsWmAxs69Eqe1bGP79ypgfh7k1mtgR4kuAMiAfdfZOZXRcuX0lwRss5wFagluAvkF4rzX1eDowAfhj+Rd3kvfhKmGnuc5+Szj67+2Yz+2+gDIgDP3L3lKdL9gZp/p7/D/ATM3uTYPjlJnfvtZcBN7OHgXnASDMrB74DFEB0n1+61IaIiKSUS0NMIiLSBQoIERFJSQEhIiIpKSBERCQlBYSIiKSkgBARkZQUECIikpICQiQi4f0mysys0MwGhPckODnbdYmkS1+UE4mQmX0XKCS4QF65u38vyyWJpE0BIRKh8DpB64F64HR3b85ySSJp0xCTSLSGE9ysZhBBT0Kk11APQiRCZraG4G5nE4DR7r4kyyWJpC1nruYq0tPCu7Y1uft/mlkMeNnMznL3Z7Ndm0g61IMQEZGUdAxCRERSUkCIiEhKCggREUlJASEiIikpIEREJCUFhIiIpKSAEBGRlP4/a90vjPm6YjUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# initial conditions\n",
    "x0, v0 = 0, 2 #starting position/velocity\n",
    "\n",
    "# initialize material/geometry parameters\n",
    "d, w0 = 20, 1\n",
    "\n",
    "#analytical solution over the full domain\n",
    "x = torch.linspace(0,1,500).view(-1,1)\n",
    "y = OD_oscillator(d, w0, x, x0, v0).view(-1,1)\n",
    "\n",
    "# slice out a small number of points \n",
    "#that will serve as training data\n",
    "x_data = x[0:10:5]\n",
    "y_data = y[0:10:5]\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(x, y, label=\"Exact solution\")\n",
    "plt.scatter(x_data, y_data, color=\"tab:orange\", label=\"Training data\")\n",
    "plt.title(\"Generated training data\")\n",
    "plt.xlabel(\"x\")\n",
    "plt.ylabel(\"y\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d68f4e47",
   "metadata": {},
   "source": [
    "#### FCN\n",
    "Fitting the FCN in the overdamped case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "b4db5a1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "torch.manual_seed(123)\n",
    "model = FCN(1,1,15,3) # (n_in, n_out, n_hidden, n_layers)\n",
    "optimizer = torch.optim.Adam(model.parameters(),lr=1e-2)\n",
    "epochs = range(1, 501)\n",
    "losses = []\n",
    "files = []\n",
    "\n",
    "for epoch in epochs:\n",
    "    optimizer.zero_grad()\n",
    "    yh = model(x_data)\n",
    "    loss = torch.mean((yh-y_data)**2)# MSE\n",
    "    loss.backward()\n",
    "    optimizer.step()\n",
    "    losses.append(loss)\n",
    "        \n",
    "    # plot results during training\n",
    "    if (epoch) % 10 == 0: # every 10 epochs\n",
    "        \n",
    "        yh = model(x).detach()\n",
    "        \n",
    "        plot_result(epoch,x,y,x_data,y_data,yh)\n",
    "        \n",
    "        file = \"./PINNplots/\"+str(epoch)+\".png\"\n",
    "        plt.savefig(file, bbox_inches='tight', pad_inches=0.1, dpi=100, facecolor=\"white\")\n",
    "        files.append(file)\n",
    "        \n",
    "        #show plot every 500 epochs:\n",
    "        #if (epoch) % 500 == 0: plt.show() \n",
    "        #else: plt.close(\"all\")\n",
    "        plt.close(\"all\")    \n",
    "            \n",
    "save_gif_PIL(\"od_nn.gif\", files, fps=20, loop=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df780a19",
   "metadata": {},
   "source": [
    "![SegmentLocal](od_nn.gif \"gif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6b38cfa",
   "metadata": {},
   "source": [
    "#### PINN\n",
    "Fitting the overdamped case data with a PINN:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "30490ce8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# sample 30 uniformly spaced points\n",
    "x_physics = torch.linspace(0,1,30).view(-1,1).requires_grad_(True)\n",
    "\n",
    "mu, k = 2*d, w0**2\n",
    "\n",
    "torch.manual_seed(123)\n",
    "model = FCN(1,1,9,3)\n",
    "epochs = range(1, 20001)\n",
    "losses = []\n",
    "losses_data = []\n",
    "losses_physics = []\n",
    "optimizer = torch.optim.Adam(model.parameters(),lr=1e-3)\n",
    "files = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "c36dfcbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "for epoch in epochs:\n",
    "    optimizer.zero_grad()\n",
    "    \n",
    "    # Data loss\n",
    "    yh = model(x_data)\n",
    "    loss_data = torch.mean((yh-y_data)**2) # MSE\n",
    "    \n",
    "    # Physics loss\n",
    "    yh_physics = model(x_physics)\n",
    "    \n",
    "    # calculate first and second derivatives\n",
    "    # dy/dx\n",
    "    dx  = torch.autograd.grad(yh_physics, x_physics, torch.ones_like(yh_physics), create_graph=True)[0]\n",
    "    \n",
    "    # d^2y/dx^2\n",
    "    dx2 = torch.autograd.grad(dx, x_physics, torch.ones_like(dx), create_graph=True)[0]\n",
    "    \n",
    "    # computes the residual of the harmonic oscillator differential equation\n",
    "    physics = dx2 + mu*dx + k*yh_physics # from (1) (beltoforion.de)\n",
    "    loss_physics = (1e-5)*torch.mean(physics**2)\n",
    "    \n",
    "    # backpropagate joint loss\n",
    "    loss = loss_data + loss_physics # add two loss terms together\n",
    "    loss.backward()\n",
    "    optimizer.step()\n",
    "    losses.append(loss)\n",
    "    losses_data.append(loss_data)\n",
    "    losses_physics.append(loss_physics)\n",
    "    \n",
    "    \n",
    "    # plot the result as training progresses\n",
    "    if epoch % 150 == 0: \n",
    "        \n",
    "        yh = model(x).detach()\n",
    "        xp = x_physics.detach()\n",
    "        \n",
    "        plot_result(epoch,x,y,x_data,y_data,yh,xp)\n",
    "        \n",
    "        file = \"./PINNplots/\"+str(epoch)+\".png\"\n",
    "        plt.savefig(file, bbox_inches='tight', pad_inches=0.1, dpi=100, facecolor=\"white\")\n",
    "        files.append(file)\n",
    "        \n",
    "        #if epoch % 6000 == 0: plt.show()\n",
    "        #else: plt.close(\"all\")\n",
    "        plt.close(\"all\")\n",
    "            \n",
    "save_gif_PIL(\"od_pinn.gif\", files, fps=20, loop=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c119e662",
   "metadata": {},
   "source": [
    "![SegmentLocal](od_pinn.gif \"gif\")"
   ]
  }
 ],
 "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}