{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Tutorial: Learning Multiscale PDEs Using Fourier Feature Networks\n", "\n", "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/mathLab/PINA/blob/master/tutorials/tutorial13/tutorial.ipynb)\n", "\n", "This tutorial demonstrates how to solve a PDE with multiscale behavior using Physics-Informed Neural Networks (PINNs), as discussed in [*On the Eigenvector Bias of Fourier Feature Networks: From Regression to Solving Multi-Scale PDEs with Physics-Informed Neural Networks*](https://doi.org/10.1016/j.cma.2021.113938).\n", "\n", "Let’s begin by importing the necessary libraries.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## routine needed to run the notebook on Google Colab\n", "try:\n", " import google.colab\n", "\n", " IN_COLAB = True\n", "except:\n", " IN_COLAB = False\n", "if IN_COLAB:\n", " !pip install \"pina-mathlab[tutorial]\"\n", "\n", "import torch\n", "import matplotlib.pyplot as plt\n", "import warnings\n", "\n", "from pina import Condition, Trainer\n", "from pina.problem import SpatialProblem\n", "from pina.operator import laplacian\n", "from pina.solver import PINN, SelfAdaptivePINN as SAPINN\n", "from pina.loss import LpLoss\n", "from pina.domain import CartesianDomain\n", "from pina.equation import Equation, FixedValue\n", "from pina.model import FeedForward\n", "from pina.model.block import FourierFeatureEmbedding\n", "\n", "warnings.filterwarnings(\"ignore\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Multiscale Problem\n", "\n", "We begin by presenting the problem, which is also discussed in Section 2 of [*On the Eigenvector Bias of Fourier Feature Networks: From Regression to Solving Multi-Scale PDEs with Physics-Informed Neural Networks*](https://doi.org/10.1016/j.cma.2021.113938). The one-dimensional Poisson problem we aim to solve is mathematically defined as:\n", "\n", "\\begin{equation}\n", "\\begin{cases}\n", "\\Delta u(x) + f(x) = 0 \\quad x \\in [0,1], \\\\\n", "u(x) = 0 \\quad x \\in \\partial[0,1],\n", "\\end{cases}\n", "\\end{equation}\n", "\n", "We define the solution as:\n", "\n", "$$\n", "u(x) = \\sin(2\\pi x) + 0.1 \\sin(50\\pi x),\n", "$$\n", "\n", "which leads to the corresponding force term:\n", "\n", "$$\n", "f(x) = (2\\pi)^2 \\sin(2\\pi x) + 0.1 (50 \\pi)^2 \\sin(50\\pi x).\n", "$$\n", "\n", "While this example is simple and pedagogical, it's important to note that the solution exhibits low-frequency behavior in the macro-scale and high-frequency behavior in the micro-scale. This characteristic is common in many practical scenarios.\n", "\n", "Below is the implementation of the `Poisson` problem as described mathematically above.\n", "> **👉 We have a dedicated [tutorial](https://mathlab.github.io/PINA/tutorial16/tutorial.html) to teach how to build a Problem from scratch — have a look if you're interested!**" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "class Poisson(SpatialProblem):\n", " output_variables = [\"u\"]\n", " spatial_domain = CartesianDomain({\"x\": [0, 1]})\n", "\n", " def poisson_equation(input_, output_):\n", " x = input_.extract(\"x\")\n", " u_xx = laplacian(output_, input_, components=[\"u\"], d=[\"x\"])\n", " f = ((2 * torch.pi) ** 2) * torch.sin(2 * torch.pi * x) + 0.1 * (\n", " (50 * torch.pi) ** 2\n", " ) * torch.sin(50 * torch.pi * x)\n", " return u_xx + f\n", "\n", " domains = {\n", " \"bound_cond0\": CartesianDomain({\"x\": 0.0}),\n", " \"bound_cond1\": CartesianDomain({\"x\": 1.0}),\n", " \"phys_cond\": spatial_domain,\n", " }\n", " # here we write the problem conditions\n", " conditions = {\n", " \"bound_cond0\": Condition(\n", " domain=\"bound_cond0\", equation=FixedValue(0.0)\n", " ),\n", " \"bound_cond1\": Condition(\n", " domain=\"bound_cond1\", equation=FixedValue(0.0)\n", " ),\n", " \"phys_cond\": Condition(\n", " domain=\"phys_cond\", equation=Equation(poisson_equation)\n", " ),\n", " }\n", "\n", " def solution(self, x):\n", " return torch.sin(2 * torch.pi * x) + 0.1 * torch.sin(50 * torch.pi * x)\n", "\n", "\n", "problem = Poisson()\n", "\n", "# let's discretise the domain\n", "problem.discretise_domain(128, \"grid\", domains=[\"phys_cond\"])\n", "problem.discretise_domain(1, \"grid\", domains=[\"bound_cond0\", \"bound_cond1\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A standard PINN approach would involve fitting the model using a Feed Forward (fully connected) Neural Network. For a conventional fully-connected neural network, it is relatively easy to approximate a function $u$, given sufficient data inside the computational domain. \n", "\n", "However, solving high-frequency or multi-scale problems presents significant challenges to PINNs, especially when the number of data points is insufficient to capture the different scales effectively.\n", "\n", "Below, we run a simulation using both the `PINN` solver and the self-adaptive `SAPINN` solver, employing a [`FeedForward`](https://mathlab.github.io/PINA/_modules/pina/model/feed_forward.html#FeedForward) model.\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "GPU available: True (mps), used: False\n", "TPU available: False, using: 0 TPU cores\n", "HPU available: False, using: 0 HPUs\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "b9db39fafd844f7fa5e3cb0b39307330", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Training: | | 0/? [00:00" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGzCAYAAAABsTylAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAABh20lEQVR4nO3dB3xT5RoG8Kd70UHpYpW996YsFREUF04UZImAMpThABXBiaJXEWXLVGQpggoyZClQ9h4tG8roYHXSnft7v5Calra0penJeP73HpumafLlJPQ8eb9x7HQ6nQ5EREREZshe6wYQERER5YVBhYiIiMwWgwoRERGZLQYVIiIiMlsMKkRERGS2GFSIiIjIbDGoEBERkdliUCEiIiKzxaBCREREZotBhagI+vbti8qVK2e7LiEhAa+88gqCgoJgZ2eH4cOHm7QN8vjSDi2cO3dOPcd58+Zp8vi2TPb7+PHji/U+5XWU+5XXlcjcMKiQTTh8+DCeffZZVKpUCa6urihfvjweeughfPfdd8X2GJ999pn6g//aa6/hxx9/RK9eve76OxkZGShXrpw6SPz1118wNz///DMmTZoEc2EISIbNwcEBwcHBeOqpp3DgwIFst5WfDx06NNff/fXXX++4bzn4y8+uXr2adZ0EQbmuYcOGyO1sIzkfw9zJe3TFihVaN4OoUBhUyOpt374dzZs3x8GDBzFgwAB8//33qvJhb2+Pb7/9ttgeZ+PGjWjdujXGjRuHl156Cc2aNSvQ71y5ckVVRxYuXAhLCSoS+G7dulWgMGYKL774ogqDc+bMQY8ePbL2fc6wkpePPvoo1+CRX9Bdvnw5LF1eQUVeR3k95XUlMjeOWjeAyNQ+/fRTeHt7Y/fu3fDx8cn2s+jo6GJ7HLmvunXrFup3fvrpJzRt2hR9+vTBu+++i8TERHh4eMDcSSVBKlNakX0mYdCgbdu2eOKJJzBt2jTMmDEj399t3LixCjS//fYbnn766bs+lpubGypWrKjCjdxenru1kcqUbETmiBUVsnqnT59GvXr17ggpIiAgINfwINUQOUD5+vrihRdeQERERJ73v3nzZnXwOnv2LFatWpXVvXC3/n75BCsHS7n/559/Xn2/cuXKO24nn/w/+eQTVKhQAe7u7njggQdw9OjRO253/fp1vPnmm2jQoAFKlSoFLy8vPPLII6qSlFt7lyxZosKRjKmRcCQHeuPnef/996vnc/78+aznZBiXk3OMyldffaW+l9vmNGbMGDg7O+PGjRtZ1+3cuRMPP/ywCpDynO677z5s27YNRdWxY0f1VV6Du5H9XbNmzQJXVaTy9v777+PQoUPq9SqK9evXo127duo9KK9NrVq11L7PGXT79++PwMBAFQIbNWqE+fPnF2m8lHFXloFcliAs92l4PQ1jnPIaozJ16lT1b8fFxUV1UQ4ZMgQ3b97Mdht5n9SvXx/Hjh1T7015PaVrdeLEiYXeT0S5YVAhqyfl7L179+LIkSMFqr707t0bNWrUwNdff60GxG7YsAEdOnS44w+0QZ06dVQ3hJ+fn/q0Lpdl8/f3z/exfv/9dzUAVw6cEhbkD35u3T8ffPABxo4dqw5cX375JapWrYrOnTurg46xM2fOqLL+Y489ptr+1ltvqS4LCQGXL1/O9blKEHnnnXfw+uuvq4Npp06dVGAS7733nno+8rwMzymv8SoStORAt3Tp0jt+JtdJe0uXLq2+l24a2Z9xcXGqm0y6I2TfStjYtWsXihpGRZkyZe56W6kcSPCQAFfQ4CHdS/KeKGyXkZBQKa9JSkqK+v3//e9/KhQaBzPZ5/L6yz7u2bOnep0lxEmQKK7uSblvCRzt27fPej0HDRqU5+0l6EgwkYAibX7mmWdUtUpey7S0tGy3lRAqwVPeo3Lb2rVrq/eVOY67IgukI7Jy69at0zk4OKgtJCRE9/bbb+vWrl2rS01NzXa7c+fOqdt8+umn2a4/fPiwztHRMdv1ffr00VWqVCnb7eT7Rx99tMDteuyxx3Rt27bN+n7mzJnqcaKjo7Ouk8vOzs7qfjMzM7Ouf/fdd+VoqdphkJycrMvIyMj2GGfPntW5uLjoPvroo6zrNm3apH63fPnyuri4uKzrly5dqq7/9ttvs66Tx835PA33K7edO3du1nWyb5s1a5btdrt27VK3W7BggfpenkONGjV0Xbp0yfZ8kpKSdFWqVNE99NBD+e4zw+N++OGHupiYGF1kZKRu8+bNuiZNmqjrf/3116zbyvdDhgy543e//PJLXXp6umpHo0aNstoxbtw49XO5XwPZvx4eHury/Pnz1c+XL1+e52Pk5ptvvrnjfnOaNGmSus1PP/2UdZ28P2WflipVKtvrJLeTthq3MbfXyPB8jMlzMX7PGMjrKLeVfWT8vuvcuXO299T333+vbjdnzpys6+67775sr7FISUnRBQUF6Z555pl89w1RQbCiQlZPZveEhoaqT7HyKVpK0l26dFHlaalqGMhgyczMTFUdkJkfhk2qHfJpetOmTcXWpmvXrmHt2rVqUKiBfGLNWZX4+++/kZqaimHDhmUr4+c29Vk+LUs3hWE2kTyGoZth3759d9xeKkeenp5Z38usqLJly2L16tVFek7du3dXlStDdUNI95K068knn1Tfy9iQkydPqgqFtM+wj6U69OCDD+Kff/5Rr8HdSCVGKlaGSpQ85hdffFGgMSc5qyoFnQUjlY6iVFUMXY7SrZfXc5N9Ls/F+P3g5OSkKl1SdduyZQtKkuF9J+8zw3tKyGB06VKUSpwxeZ8ZjxmSrr6WLVuqKh/RvWJQIZvQokULFUSkRC3dCzJuIj4+Xh2cpW9dyAFUDkByMJKDoPF2/PjxYh14KwdwKZ83adIEp06dUpuMMWnVqlW27h/DmA9pkzFpk6ErxUAOgt988426rYQD6bKR28nYitjY2DvakPM+JQhVr169yGtpPPfcc+qgJs9NyL5ctmyZGicjBzfDPhYyeDjnPv7hhx9U90hubc1p4MCBqqtKuuUkHMlr8/bbbxeqvRI85PkWNHgYwo2ErcJM8ZUAJ4N9ZaaZjD+Rrj4Jo8ahRV5neT2MQ4GhW9Hw85JkeDwJucYkgEjXY872yPipnIOM5f1pPC6JqKg464dsivyhldAimwyo7NevnzqYyid0OXAY1jPJbQaEfGosLoYwIgew3MgnUTkgFIaM9ZCxLC+//DI+/vhjNRBYDnzyqbggVYp7JWMZZPyDHIRloOiOHTtw4cIFVekwMLRDxmDI+JfcFGQ/y0FdxtPcC0PwkHEguQ1izivcyL6VcNOtW7cC/Y4MypZKkVTkpBKxZs0aFeZkTM66devuebZNXrOQpKpWUvJ6DoUdz0OUGwYVslmytoqQdUxEtWrV1B/WKlWqqBBjKjIzRdZ2kYXCZKCrMTmQy5oWsn6JHEQN61pIJcI4uMTExNzxafWXX35Rsy5mz56d7XoZqCrVlZwM1Q0Dee5S2ZHFzQwKOxVXqgeDBw9GeHi4OhjLDJDHH3886+eyj4VUWO41aBQH6a6QGVUffvih6ho0RbgREhila0s2GegsoVIGK0t4kf0gr7NUvuT1N66qhIWFqa/5rW8ilYvcBnrnVoUp6OtpeDx5HY3fd9IdJO9fc3jtyHaw64esnhwMcvtkZxiLYShvy/gGORDJQSvn7eV7GVNRnNUU6aqQrifjTcbHSHgx3EYOCDJWQVbQNW5TbrNvpO052y3VokuXLuXajgULFqjuL+OgI6FNumoMZNpyQbpijMfZSDsWLVqkHltmuxivCyPTviWsyHRmGXuRkwSwkmTcnWM8Xulu4Ua6jOR9UhDSpZeToZokXV2ia9euiIyMzOo2E+np6ep1lwpTzkBrTPanvEYSdAzkdcxtRpO8FnnNXjMm7zupPk6ePDnbe0pCsDzWo48+etf7ICourKiQ1ZOBqElJSWqZdZk2KZ8KpaIhBwVZf0K6fwx/8OXTtYxfkXEaUtqXwabyCVL+6Mu4CFmn5F5JCJEDlSwilhv5ZC9tlgGwsrCZPOaECRPUQV8OaPv371fdUzmrJPJz6ZKQ59OmTRs1NVkeK68uJOkakrU95PZRUVEq/MgBWAZMGgcL2U8jR45U3WVy0DSukOS2Lo1UdaRqICFIKizGpFogY1EkDMn6HPLYMqhZwpQESqm0/PHHHyhJhu6cgq5qK+FGqiGG983dyGsiXT9ycJdKhYynkfVJZFyH7H8h7y2Z+iuVGhlzI+9LCY4yhVleF+NBzznJmBeZCizvbxl8K+91WfhOqoI5B1HL6ykDZeX1ka46qR7KuKicZMyQ/DuQMCbTjuU9KdUVabe8D4wHzhKZXIHmBhFZsL/++kv38ssv62rXrq2mesq0y+rVq+uGDRumi4qKuuP2MsW1Xbt2aiqnbPJ7MgU1PDz8nqcn7927V03lHDt2bJ63kWnScpsRI0ao72V6qEzHLVu2rM7NzU13//33644cOaIeL+f05FGjRmXdTqY+h4aGqumjsuWcnrxo0SLdmDFjdAEBAer20vbz589na0tCQoKuR48eOh8fH/U7huec2/Rkg1mzZqmfeXp66m7dupXrc9y/f7/u6aef1pUpU0ZNn5b7ff7553UbNmzId/8ZTzG+m/ymJ+c1PTe/6cnG0tLSdNWqVSvQ9GR5Tk8++aSuXLly6r0nX1988UXdiRMnst1O3ov9+vXT+fn5qds1aNAg1/2bc3qyYQp+/fr11e/VqlVLTXPObXpyWFiYrkOHDur1Np7ennN6svF0ZHn/Ozk56QIDA3Wvvfaa7saNG9luI++tevXq3dHOvKZNExWWnfzH9HGIiMyFrEwrVQ/pmpHuJiIic8YxKkRERGS2GFSIiIjIbDGoEBERkdniGBUiIiIyW6yoEBERkdliUCEiIiKzZfELvsmS05cvX1YLIhV2uW8iIiLShow8kYUhZfHBnCfktKqgIiElrxU+iYiIyLxFRESolZqtNqgYlpaWJ2o4lTwRERGZt7i4OFVoyO8UEVYRVAzdPRJSGFSIiIgsy92GbXAwLREREZktBhUiIiIyWwwqREREZLYsfowKEdG9TpFMT09HRkaG1k0hsioODg5wdHS856VDGFSIyGalpqbiypUrSEpK0ropRFbJ3d0dZcuWhbOzc5Hvg0GFiGySLBZ59uxZ9alPFpySP6RcNJKo+CqV8kEgJiZG/TurUaNGvou65YdBhYhskvwRlbAi6zjIpz4iKl5ubm5wcnLC+fPn1b83V1fXIt0PB9MSkU0r6qc8IiqZf18m/Rc6bdo0NGzYMGsxtpCQEPz1119ZP09OTsaQIUNQpkwZlCpVCs888wyioqJM2SQiIiKyICYNKrJ2/+eff469e/diz5496NixI5588kkcPXpU/XzEiBH4448/sGzZMmzZskWdt+fpp582ZZOIiIjIgpg0qDz++OPo2rWrGkRTs2ZNfPrpp6pysmPHDsTGxmL27Nn4+uuvVYBp1qwZ5s6di+3bt6uf5yUlJUWdH8B4IyIi83P//fdj+PDhsHbjx49H48aNS+zx5s2bBx8fn3u+n82bN6sB5Ddv3oQ5K7HOWVmjYPHixUhMTFRdQFJlSUtLQ6dOnbJuU7t2bQQHByM0NDTP+5kwYQK8vb2zNp45mYhsTd++fdUBRirWxlasWGFRM5fkgCvtffjhh7NdLwdOuV4OpIXZJ926dTNBK603NLZp00ZNz5djqTkz+ayfw4cPq2Ai41GkmvLbb7+hbt26OHDggJoOmDMVBgYGIjIyMs/7GzNmDEaOHHnH2ReL2/ErcTh2OQ5XE1Jw81Ya0tIz4ehgD79Szqji54GGFXzg7+lS7I9LRFQQMoPiiy++wKBBg1C6dOkSfWz5kCmzOYqDLAj2999/Y9OmTXjggQdgaVNwLXmhQGdnZwQFBcHcmbyiUqtWLRVKdu7ciddeew19+vTBsWPHinx/Li4uWYNzTXnG5OX7LmLUsoOY8FcYpm0+jR+2nsX0Lafxyarj6D9/D1p8+jcenvQPvl4XjsMXY9Ublogsm/w7TkpN12Qr7N8QqUbLQUaqzPnZunUr2rdvr6aKyoe6119/XVW2DaRyIZUYY/IBUqod4ty5c+o2S5YswX333acC0sKFC3Ht2jW8+OKLKF++vJre3aBBAyxatAiF5eHhgZdffhmjR4/O93YRERF4/vnnVdt8fX3VeEdpm6HrZf78+Vi5cqVqq6Ea8+yzz2Lo0KFZ9yEVBflZWFiY+l6mzMrjS1AyDC2Q/RMQEKCeZ7t27bB79+47ukpkUogMV5DjkezfnE6fPo2qVauqx87tdZXrpM3SgyD3Iev4yOMa3LhxA71791YB1N3dHY888ghOnjxZqGqSPFepohh+LuNAv/3226z9I/sut66fX3/9FfXq1VPtqly5Mv73v/9lu1+57rPPPlOvmaenp3oOM2fOhEVXVCSxVa9eXV2WF1ZedNlZ3bt3V28S2UHGVRWZ9WMOCa92kBfa1/CDfykX+Lg7w9nRHinpGYiJT0F4ZDxOxSQgLDJebZM3nkJZb1d0qhOIttX9ULesF8r6uMLJgdMeiSzJrbQM1P1grSaPfeyjLnB3LvifZFmoTg4YPXr0UAc5mbyQ2wFTulU++eQTzJkzRy2+JQdP2WRMYGFIkJCDVpMmTdRBXKrk8jf9nXfeUR8YV61ahV69eqFatWpo2bJloe5bDtpynPjll19UuMitgtOlSxdVnf/3339VFUaekzy3Q4cO4c0338Tx48dVhd3wvCTMSEV/xowZWfcjB2s/Pz91gJahBnI8kvuWLhDx9ttvqwO1hJ5KlSph4sSJ6nFPnTql7s94X3z11VcqjEiYMO6ikvbI7/Tv31+1MTfyGN98840aDiGhQHoRDh48mPVzCRYSTH7//Xe1b2Ufy3hP+ZBflEqWHHNPnDiB+vXr46OPPlLX+fv7ZwU9AxmSIWFQXg85RsuY0cGDB6uZudImA3kffPzxx3j33XfVayZFCAmxUpgwhRJf8E0WWJLUKm9w2eEbNmxQ05JFeHg4Lly4oN6MWnumWQW15eVGYio2hUdj/bEobDkRgyuxyfhxx3m1GXi5OqJMKRd4uTnBwQ5wtLeHTCnXf7WDo70d7FV3csn3KWvRjW13+3Htbj9fddnwvf7/+rRv1MZs1+X83dtf7e3s4ObsAE8XR3i4OKK0uzPKl3ZD9YBSCPQq2gJDRJbgqaeeUoM4x40bpyYn5CTVlp49e2aNTZCJDZMnT1YHFVk+ojALcMl95JyVKQHBYNiwYVi7di2WLl1a6KAiFYU33ngD7733Xq7jTKSaI8eOH374IWsMjgQS+ZArIaFz586qYiTHFuMPulJRkPuVgCbhRg70Y8eOVb/z6quvqq8tWrRQVQupMsk+kUqSVDDErFmzsH79erVv33rrraz7lYP9Qw89dEc75cD+2GOPqecxatSoPJ+vHOeknVIVk+OgVCUM+8wQULZt25YVoBYuXKiqYVL5eu6551BYMgZFigbyPPMrBMjklgcffFDtIyGTYGSfffnll9mCioQmCTBCQpSELum6s8igIuNJ5AWXFyE+Ph4///yzemPIm1l2nCROGW8iSVVSo7zRJaS0bt0a5q60hzOeblpBbclpGdh++irWH4vGwYibOBkdj7QMHeKS09VG2gn0ckH7Gv7o2iAI99UMgIM+GRLlys3JQVU2tHrsopBxKjJz0jg0GMindPmELwc6424Hw+kD6tSpU+DHad68ebbvZWyGVHQkmFy6dElVyCUoFHWVXzngSfVDKj/yqT7n85CqhnQ1GJOqjlSN8iIVBDm+SCVFDtRSDZIgMWXKFPVzud7QPSL3I9WVtm3bZv2+hAgJEFKtyW9fGMKHhBeZ3Xq3mU4SNiZNmqQqMlIVkgO/zJKVMCWPJV9btWqVdfsyZcqoEJCzHcVN7l+61IzJ/pC2yustVTwh66MZSHCU8BMdHW2ydpk0qEjDpZ/NMKpYnpyEFEMSlRQmq9ZJRUXe4FIumzp1KiyNq5MDOtYOVJvIyNQh9lYariem4FpCKuKT05Ehg64ydUjPlD8S/32V60uaFsNpdPI/nXz9rxFyWV1nfDnrx/pL+uv++11D2w3XCdmPSWkZSEhOR2JKOq4mpuLijSScu5qIqLgU/LL3otqq+nlg8APV8WTjcuyWo1zJH93CdL+Ygw4dOqi/nfLB0PhTr0hISFCDbY3HPxjIB0jDc845jkIO2DnJWA5j8ilbuhTkICbjU+TncoCWwFIUUh2R5/Dhhx+qMJHzeUgV3jhwGUgXRl7kucn+kQ/IMuZCQokch+R4c+TIEVUByS3g3U3OfWFoh1SGZJyOjN/Ib/ykVEekB0HGxkjFRqoTsj8lOBWFvb19gV7D4pKz+0n2s4RfUzHpv8jcSpHGpOwoydaQbq2FfGr39XBWW/UArVtju2SA4oELN7HuWBRWHLiEM1cT8eayg/h2wwkMvr86nmlaQY09IrJ0Mk1ZuoBylt6bNm2qSveGcYK5kQOsfJg0kK6HgpxNWrom5NP3Sy+9pL6XA5WMg5BZnUUlVXXpmpIAlPN5SPePDHLNKwBIxSS3GTjSzSVdOBJUpNohB3UJLxIMJLAYKigytkbuQ56XjE8xHOxlHEtB1oKRrqc///xTVUckOK5bt+6OClDO20sVRTZZoV3GzMiYGqlypaenqwkohq6fa9euqWCT176V11CClzGZxGIcKPLaP8bkseX5G5PvpQvIUE3RAv9Kk9WST8Ztqvth/BP1sPWdjhj9SG01vTzi+i2MWX4Y93+5CT/8ewbXE4v2CZDIXEhFQ8aiyEE+Z3eKVA1k8KwcuCSEyMwY45kw0m30/fffY//+/WoFcRm7UZABmzLeRaoBcv/SZSCVm3s9BYp8eJWKSs7nIc9NBsFKMJLBtNJtJVUSqRRdvHgxazaKdHPJAf3q1atZFQWpokhYkxXRZRaP4TqpzkgXjqE6Il9lUKiMRVmzZo36nQEDBqjQJsMUCkLuQwYVS9eNDHuQSlBuZByMfJCXcHHmzBn89NNPKrhIQJL9Ks9THltmFB08eFCFQZldlbNbxvg1lNduwYIF6jWWMUs5g4vsHwk/MoBW9k9uFRAZVyPjRmWgrIROGVQs742iVJ2KE4MK2YRSLo549b5q+Pftjhj7WF0EeLrgcmyymm7e6rO/0X1GKCZvOIk9564jLcN0JUwiU5EBnjkPPtLNId0JctCRKcoyRuODDz5QXRTGMzikK0J+LjOI5KBUkHEm77//vqp0SPVADvwyTqE4FlyTJSxk7IYxac8///yjuqtkQK988pfwIGNUDBUWObBLRUnCh1QYDJUBCXHSrSQVJ1nLS0h7pbpgGJ9iXJmSoQgye0mem4yLkeEKhVmnRh5Dpi9LV8yjjz6abSq4gbRHqjxSzZHXSLqA5HQyMhbFMFBYurqkCywkJETd1+rVq/MMkPIayABYmbUkg4NlTKgMuzAmr6tURaQqI/tHxtTkJM9ZxhzJbCQZ3yPvFXlf5exSLGl2OgtfAESmo8n4F1mS31RrqpD1kQHQy/ddwqJdF3D4Umy2n7k7O6BZpdJoVcUXrauWQeOKPmqxP7IucpCTT+ZVqlQp8unniajo/84Kevy2rFFjRMU4ALpHq2C1nYlJwPbT1xAq25lrqivo35NX1Sak+vJUk/Lo27Yyynq7ad10IiKbwqBCNq+qfym1vdS6kppBdCI6HjvPXMfOs9ew7dQ1RMenYMY/ZzB32zm82LIiRj5UC97uxbN8OBER5Y9BhciILMQnqxLL1qdNZaSmZ6qF/WZvPYtdZ69jfuh5rDp8BeMer4fHG/3Xz09ERKbBjneifMj05S71grB0UAgWvtIK1fw9cDUhFcMW7VdbbJLp1iogIiIGFaICk/M4rX6jPd54sIZaK+ePg5fRZdI/2Hp7LAsRERU/BhWiQnBxdMCIh2ri19faoIqfByLjkvHS7J0Y/eshdf4nIiIqXgwqREUgU5ZXvd4OvVrrV7BcvDsCD/xvM6ZsOoWEFJ7fiYiouDCoEN3Dyrcfd6uPX14NQe0gT9xMSsOXa8PRZsIGfLrqGCKu330ZciIiyh+DCtE9al7ZF38Oa4dvujdSJz6UM2bP+vcs7vtyEwYu2KPOrG3h6yoSEWmGQYWoGMjKtU81qYC/R96HOX2bo30NP2TqoE6I2GPWTjzy7b9YvOuCWhGXiKDO1SNn3b158+Y93Y+cu0buR85lRNaJQYWomNdh6Vg7ED/2b4X1IzrgpdbBcHNyQFhkPEYvP4z2Ezep2UKssFBRyUE5v238+PGwVnLOmZznE5LzFMnZn+XcNGSduOAbkYnUCPTEJ90a4K3OtbF0TwTmbT+HSzdvqfVXfj94GV8+2xA+7s5aN5MsjByUDZYsWaJOHCdnDDYwnHhPSCCWk+/J2XytlZxoT06ISNaLFRUiE5Pl9gd0qIqNb96H4Z1qwNnBHuuPReHRyVtx6OK9lb2pmEmlKzVRm62AVTY5KBs2OaGbVFEM34eFhcHT01OdvVfOvuvi4oKtW7fmWokYPnx4trMHy5mXJ0yYoE4e5+bmhkaNGuGXX37Jty1Tp05FjRo11MnmAgMD8eyzz2b9LCUlBa+//joCAgLUz9u1a4fdu3fneV9SCZIzHBubNGkSKleunPXz+fPnY+XKlVnVI+k+yq3rR84Y3bJlS/X8y5Yti9GjRyM9/b/ZePK8pW1ytmFfX1+176y5EmXprDdmE5nhGizDO9VEpzqBGPLzPpy/loTnpodi4rMN8WTj8lo3j0RaEvCZRqdGePcy4OxRLHclB+avvvoKVatWRenSpQv0OxJSfvrpJ0yfPl2Fj3/++QcvvfQS/P39cd99991x+z179qiD/Y8//og2bdrg+vXr+Pfff7N+LiHg119/VeGiUqVKmDhxIrp06YJTp06pcFBYb775Jo4fP67OuDt37lx1ndzP5cuXs93u0qVL6Nq1qwpnCxYsUOFtwIABKiwZhxFp18iRI7Fz506Ehoaq27dt2xYPPfRQodtGpsWgQlTC6pf3xh/D2mHE4gPYEBaNNxYfwLErcXi7S2214i3Rvfroo48KdcCV6sdnn32Gv//+GyEhIeo6CTlSjZkxY0auQeXChQvw8PDAY489pqo4EkaaNGmifpaYmIhp06Zh3rx5eOSRR9R1s2bNwvr16zF79my89dZbhX5O0qUllR5pa35dPVLlkXEr33//vaq01K5dW4WZd955R3WT2dvrOxIaNmyIcePGqcsSzOT2GzZsYFAxQwwqRBrwcnXCzN7N8b914Zi6+TRmbDmDE5Hx+PbFJupnpBEnd31lQ6vHLibNmzcv1O2lypGUlHTHQTo1NTUrfOQkt5VwIoHm4YcfVttTTz0Fd3d3nD59GmlpaapCYeDk5KS6Y6QqYkpy/xK2JKQYSDsSEhJw8eJFBAcHZwUVY9JFFB0dbdK2UdEwqBBpRKonbz9cG7WCPPH2L4ewKTwGT03Zhlm9m6Oq/38DIqkEycGtmLpftCSVDmNSRcg500yChIEcxMWqVatQvnz2bkgZ55EbqaLs27dPjRNZt26dqlZI10p+41Dyc7c2FjcJTsYk2Mg4HTI/HExLpDEZn/LLq21Q1tsVp2MS8dh3W/Fj6DlkykIsRMVAxpkYzxYSxoNP69atqwKJdOdUr1492ybdKHmR2USdOnVS408OHTqkBrZu3LgR1apVg7OzM7Zt25YtdEiIkcfKq42RkZHZwkrOtVHkPmUWU37q1KmjxpwY34+0Q4JVhQoV8v1dMk8MKkRmoEEFb6wc2hatqvgiKTUDY1cexVPTtmPnmWtaN42sQMeOHdXgVxlcevLkSTU248iRI1k/l4O4DFYdMWKEGmQqXTdSLfnuu+/U97n5888/MXnyZBUmzp8/r+5bKhK1atVSFZ3XXntNjUVZs2YNjh07pga0SvdS//79c70/mYkTExOjQo88/pQpU9TsJWMyA0gCkUzHvnr1aq4Vl8GDByMiIgLDhg1TA2lllpA8Xxk4axifQpaFrxqRmQjwdMWiAa0x/vG6cHd2wMGIm+g+cwdemb8bJ6PitW4eWTCZbTN27Fg1E6dFixaIj49H7969s93m448/VreR2T9SlZAxJ9IVJNOVc+Pj44Ply5erECS3l9lCixYtQr169dTPP//8czzzzDPo1asXmjZtqsbBrF27Ns9ZSHIfMhBWAopMjd61a5cKT8Yk7EgQkjE4UoExrtgYSNfV6tWr1e/L/bz66qsqHL3//vv3sAdJS3Y6C18iU6aqyVoCsbGx8PLy0ro5RMUiOj4Z3/59Up2VOSNTB5kM9GyzChjxUE2U9XbTunlWITk5GWfPnlUHYpm6SkQl+++soMdvVlSIzLS68ulTDbBuRAd0qReozhu0dM9FdPxqC6ZvOY20DA76IyLbwKBCZMaq+ZfCjF7NsXxwG7SoXBq30jLw+V9h6DZlG87E6GdqEBFZMwYVIgvQNLg0lg4Kwf+ea4TS7k44ejlOzQ5afTj7TA4iImvDoEJkIWSdh2eaVcBfb3RA66r62UGDF+7DtM2neTZmIrJaDCpEFibI2xULX2mNvm30J2v7Yk0YPv7zOMNKEXG/EZn3vy8GFSILXdV2/BP1MPYx/eJZc7adxTu/HlIzhKhwK5PK2h5EZBqGf185VwIuDC6hT2TB+rerAi9XRxVSZFZQYkoGvuneGM6O/AxyNw4ODmotEMP5XeQcNcbnhyGie6ukSEiRf1/y70z+vRUVgwqRhXuueUWUcnHE64v3Y9XhK0hMTce0ns3g5lz0Pwy2wnAWXp6Mjsg0JKTkd7brguCCb0RWYsuJGAz6cQ+S0zLRsoovZvdpDk+eiblA5PwxpjwBHpEtcnJyyreSUtDjN4MKkRXZdfY6+s/bjfiUdNQIKKXOxFzZz/LPBkxE1ocr0xLZIKmkLBrYGgGeLjgZnYAnvt+KlQcucWYLEVksBhUiK1O/vDf+HNYOTYN9EJecjjcWH8BrP+1DTHyK1k0jIio0BhUiKxTg5Yolg0IwolNNONrbYc3RSHT+Zgt+P3iZ1RUisigMKkRWysnBHm90qoHfh7ZD3bJeuJGUhtcX7Wd1hYgsCoMKkZWrW84LK4e2xfBONbKqKw99s4VjV4jIIjCoENlIdWV4p5pZ1ZWbSWlq7IosFJeSnqF184iI8sSgQmSj1RV7O6jVbF+YuQPXE1O1bhoRUa4YVIhstLoyr19Ltfz+/gs30X1GKCJjk7VuGhHRHRhUiGxUh5r+WD64LYK8XNWaK88zrBCRGWJQIbJh1QNKYdmrIQj2dceF60no+cMOXE3gjCAiMh8MKkQ2rqKvO34e0ArlvF1xOiYRL/2wEzeTOGaFiMwDgwoRoUJpdywc0Br+ni4Ii4xHnzm7EJ/Mk/QRkfYYVIhIqeLngYWvtEJpdyccvBiL/vP2ICk1XetmEZGNM2lQmTBhAlq0aAFPT08EBASgW7duCA8Pz3ab5ORkDBkyBGXKlEGpUqXwzDPPICoqypTNIqI81Az0xI/9W8HT1RG7zl3HwAV7kZzGdVaIyEqDypYtW1QI2bFjB9avX4+0tDR07twZiYmJWbcZMWIE/vjjDyxbtkzd/vLly3j66adN2SwiustJDWXqsruzA7aeuorec3Yhjt1ARKQRO10JrqEdExOjKisSSDp06IDY2Fj4+/vj559/xrPPPqtuExYWhjp16iA0NBStW7e+633GxcXB29tb3ZeXl1cJPAsi27Dr7HX0n7cb8SnpqB3kiZm9miO4jLvWzSIiK1HQ43eJjlGRxghfX1/1de/evarK0qlTp6zb1K5dG8HBwSqo5CYlJUU9OeONiIpfyyq+WDSwNfxKOasBto999y82hrFblohKVokFlczMTAwfPhxt27ZF/fr11XWRkZFwdnaGj49PttsGBgaqn+U17kUSmGGrWLFiibSfyFa7geT8QI0r+iAuOR0vz9uDr9eFIyOTJzMkIisLKjJW5ciRI1i8ePE93c+YMWNUZcawRUREFFsbiehO5XzcsGRQa/RqXUl9P3njKfSesxPR8VzFloisJKgMHToUf/75JzZt2oQKFSpkXR8UFITU1FTcvHkz2+1l1o/8LDcuLi6qL8t4IyLTcnF0wMfd6uOb7o3g5uSAbaeuoeu3W7H91FWtm0ZEVs6kQUXG6UpI+e2337Bx40ZUqVIl28+bNWsGJycnbNiwIes6mb584cIFhISEmLJpRFQETzWpgN+HtkXNwFJqqf2es3fim/UnkMmuICKyxFk/gwcPVjN6Vq5ciVq1amVdL2NL3Nzc1OXXXnsNq1evxrx581R1ZNiwYer67du3F+gxOOuHqOTdSs3A+N+PYskefdfrYw3L4qvnGsHVyUHrphGRhSjo8dukQcXOzi7X6+fOnYu+fftmLfg2atQoLFq0SM3o6dKlC6ZOnZpn109ODCpE2lm2JwJjlh9GeqYOrar4Yk7fFvBwcdS6WURkAcwiqJQEBhUibW07dRWv/rhXrbciU5rn9WsBd2eGFSKywHVUiMj6tK3uhwX9W8LTxVEtEjdgwR6kpHPZfSIqHgwqRHTPmgSXVmHFw1k/I2jk0oMcYEtExYJBhYiKLaxM79UMTg52WHXoCj7846ia+UdEdC8YVIio2LSv4Y//Pd8YMo5+fuh5fL/xlNZNIiILx6BCRMXqiUblMO6xuury/9afwM87L2jdJCKyYAwqRFTs+ratgiEPVFOX319xGGuOXNG6SURkoRhUiMgk3uxcC92bV4SMqX198QHsOHNN6yYRkQViUCEiky34+OlT9fFQ3UCkpmdiwPw9OHo5VutmEZGFYVAhIpNxdLDHdy82QcvKvmpBuB6zdmLv+etaN4uILAiDChGZlJz/Z1af5mhc0Qext9LQ84edWHc0UutmEZGFYFAhIpPzdnPCzwNa4f5a/khOy8Sgn/Zi2ubTXGeFiO6KQYWISoSc/2dW7+Z4qXUwJJ98sSYMby47xOX2iShfDCpEVGKcHOzxSbcG+PCJerC3A37ddxE9Z+3EtYQUrZtGRGaKQYWISlyfNpUxr19LeLo6Ys/5G3hyyjaER8Zr3SwiMkMMKkSkiQ41/fHb4LaoVMYdF2/cwtNTt2FjWJTWzSIiM8OgQkSaqR5QCisGt0Xrqr5ITM1A//l7uOQ+EWXDoEJEmirt4YwFL7fCiy0rqkG27/52GL/svah1s4jITDCoEJHmnB3t8dlTDdC3TWX1/du/HORaK0SkMKgQkdksuT/u8bp4oYX+/EBvLD6AI5e45D6RrWNQISKzCiufdKuP9jX8cCstA6/M34PI2GStm0VEGmJQISKzOz/QlJ5NUSOgFCLjkvHKgt1ISk3XullEpBEGFSIyO16uTpjdpwV8PZxx5FIchi8+gAzpDyIim8OgQkRmKbiMO2b2agZnB3usOxaFT1Yd07pJRKQBBhUiMlvNK/vif883UpfnbjuHOVvPat0kIiphDCpEZNYeb1QOox+prS5/vOoY1hzhtGUiW8KgQkRmb1CHqujZSn/W5TcW78eus9e1bhIRlRAGFSKyiGnLcsbljrUDkJKeib5zd2HnmWtaN4uISgCDChFZzLTlqT2bqjVWklIz0Hfubp7EkMgGMKgQkcVwdXLArN7NcX8t/6wF4XgSQyLrxqBCRBYZVp5rVkEttS8nMfxqbTh0MoCFiKwOgwoRWRwnB3tMfLYh3niwhvr++02nMGrZQaSmZ2rdNCIqZgwqRGSxA2xHPFQTXzzTAA72dli+7xJenrcb8clpWjeNiIoRgwoRWbTuLYLxQ5/mcHd2wNZTV/Hc9FCeyJDIijCoEJHFe6BWAJYMDIFfKReERcbj6anbcP5aotbNIqJiwKBCRFahQQVv/Da4Dar6e+BybDJ6/rATV2Jvad0sIrpHDCpEZDUq+rpj8cDWqFzGHRdv3EKv2bsQxzErRBaNQYWIrEqApyt+eqUVgrxccSo6AcN+3o8MmcdMRBaJQYWIrE6F0u5qgK2rkz22nIjBF2vCtG4SERURgwoRWaX65b3x5bON1OWZ/5zBL3svat0kIioCBhUislqPNyqHoQ9UV5ffXX4Ye8/f0LpJRFRIDCpEZNVGPlQTD9UNRGpGJgb9uBeXb3ImEJElYVAhIqtmb2+HSd0bo3aQJ64mpGDAgj1ISk3XullEVEAMKkRk9TxcHNWJDH09nHH0chzeXHYQmZwJRGQRGFSIyGbWWJn+UjM4Odhh9eFIfLvhpNZNIqICYFAhIpvRsoovPulWX12WoPLzzgtaN4mI7oJBhYhs7iSGr91fTV1+b8VhrNh/SesmEVE+GFSIyOa83aUWerWuBJ0OGLH0AH7acV7rJhFRHhhUiMjm2NnZ4cMn6qFnq2AVVt5fcQSTN5yETr4hItsJKv/88w8ef/xxlCtXTv1hWLFiRbafyx+FDz74AGXLloWbmxs6deqEkyc5wI2ISmbasoxXef3BGur7r9efwId/HONsICJbCiqJiYlo1KgRpkyZkuvPJ06ciMmTJ2P69OnYuXMnPDw80KVLFyQnJ5uyWUREinyAkgXhxj9eV30/b/s5jFx6AGkZmVo3jYhus9OVUK1T/iD89ttv6Natm/peHlYqLaNGjcKbb76prouNjUVgYCDmzZuHF154oUD3GxcXB29vb/W7Xl5eJn0ORGS9Vh64hFFLDyI9U4eOtQMwtWdTuDo5aN0sIqtV0OO3ZmNUzp49i8jISNXdYyANbtWqFUJDQ/P8vZSUFPXkjDcionv1ZOPyalE4OePyxrBo9J69C3HJaVo3i8jmaRZUJKQIqaAYk+8NP8vNhAkTVKAxbBUrVjR5W4nINjxQOwA/9m8FT1dH7Dp3HT1m7WBYIdKYxc36GTNmjCoTGbaIiAitm0REVqRFZV8sGRiCMh7OOHIpDq/M24PktAytm0VkszQLKkFBQeprVFRUtuvle8PPcuPi4qL6sow3IqLiVLecFxb0bwlPF31l5Y3F+zkbiMjWgkqVKlVUINmwYUPWdTLeRGb/hISEaNUsIiKlXjlvzO7bAs4O9lh7NAqTN3LpBCKrCyoJCQk4cOCA2gwDaOXyhQsX1Cyg4cOH45NPPsHvv/+Ow4cPo3fv3momkGFmEBGR5ucGekp/bqBJf5/E2qN5j58jIgsMKnv27EGTJk3UJkaOHKkuyyJv4u2338awYcMwcOBAtGjRQgWbNWvWwNXV1ZTNIiIqsOebV0TfNpXV5ZFLDuBEVLzWTSKyKSW2joqpcB0VIjI1WQBOpiuHnrmGSmXc8fuQdvB2d9K6WUQWzezXUSEishRODvaY0rMpyvu44fy1JLyxZD8yOLiWqEQwqBARFYCvhzNm9m6mFoTbHB6Db9af0LpJRDaBQYWIqBAzgb54pqG6/P2mU1hz5IrWTSKyegwqRESFXGq/f7sq6rKcG+gkB9cSmRSDChFRIY15pDZaV/VFYmoGXp6/G9HxPOM7kakwqBARFZKjDK7t0RTBvu6IuH4LfefsRjzPCURkEgwqRERFUKaUCxa83BJ+pZxx7Eoc+szZhdhbDCtExY1BhYioiCr7eWBev5bwdnPCvgs30fOHHbiRmKp1s4isCoMKEdE9qF/eG4sGtFbTl+Vsyy/O2oGY+BStm0VkNRhUiIiK4WzLSwa2RoCnC8Ii49F9Zigu37yldbOIrAKDChFRMagR6Imlg0LU6rVnYhLx3PRQnLuaqHWziCwegwoRUTGOWVn6agiq+Hng0s1beG5GKMIjuc4K0b1gUCEiKkZSUZHKSu0gTzVWRbqBDl28qXWziCwWgwoRUTHz93TBkoEhaBLsg5tJaeg9ZxdXsCUqIgYVIiIT8HZ3wk/9W6FxRX1YeWn2TtUdRESFw6BCRGQiHi6OmNu3BWoElEJUXAoG/bgHyWkZWjeLyKIwqBARmVBpD2fMe7ll1jorY5Yfhk6n07pZRBaDQYWIqAQG2Mq5gRzs7fDb/kuYs+2c1k0ishgMKkREJSCkWhm817WOuvzZ6uPYfuqq1k0isggMKkREJaRf28p4uml5ZGTqMHTRfq5eS1QADCpERCXEzs4Onz3VAPXLe+F6YioGL9yHlHQOriXKD4MKEVEJcnVywLSezeDl6ogDETfx6arjWjeJyKwxqBARlbCKvu6Y9EJjdXlB6Hms2H9J6yYRmS0GFSIiDXSsHYhhHauryzJlOSwyTusmEZklBhUiIo0M71QT7ar74VZaBvrP24Po+GStm0RkdhhUiIg0IuuqfPdik6yzLb8yfw8SUtK1bhaRWWFQISLSeOVaWWa/tLsTDl2MRb+5uxhWiIwwqBARaayynwfmv9wSnq6O2H3uBsMKkREGFSIiM9Cwgg8WvtKKYYUoBwYVIiIzDSt95uxCXHKa1s0i0hSDChGRGYYVWRBu7/kb6DFrh1rFlshWMagQEZlhWFk8MARlPJxx5FIcXpgZiug4Tl0m28SgQkRkhuqW88KSQSEI9HLBiagEPD8jFJGxDCtkexhUiIjMVPWAUlg2qA0qlHbDuWtJeGn2TnYDkc1hUCEiMmPBZdyxeGBrBHm54lR0AvrO3YWkVM4GItvBoEJEZOYqlHbHT6+0zFoUbtTSg8jM1GndLKISwaBCRGQBqgd4Ymbv5nB2sMdfRyLxzd8ntG4SUYlgUCEishAtKvtiwtMN1OXvNp7CuqORWjeJyOQYVIiILMgzzSqgX9vK6rJ0AZ29mqh1k4hMikGFiMjCvNu1DppXKo34lHS8+uNeDq4lq8agQkRkYZwc7DG1Z1P4e7ogPCoeY5Yfhk7HwbVknRhUiIgsUICXK6b0aAoHezusPHAZ87ef07pJRCbBoEJEZKFaVvFV3UDik1XHsefcda2bRFTsGFTIOkkZPD0VSIkHMjO1bg2RybzctjIea1gW6Zk6DF64D9HxXGafrIuj1g0gK5CWDNw4B9y8AKTfAmAHePgD3uUBrwqAfT55ODMDiAkDLu4BrhwArp0G4iOB9GTA0RUoFQCUrgT41QL8agJ+NQCfSoCD0Vv31k0g+jgQfRSIur3J9ylx+p/bOQClAgH/WkBgPSCogf6r3J+jy3/3k56if/yr4UDMCf3Xa6eAlARAlwG4lwG8ygEBdfWb3IdvVcDeIf/9kxwLxF7Ubxmp+n1Tpjrg4Xeve54IdnZ2+OKZhgiPjMfJ6AQM/Xm/OvuyjGMhsgZ2OgsfgRUXFwdvb2/ExsbCy8tL6+ZYlow0/UH98n7g1g0pQ9wOGBUBn2DAu0L2A7mQt8v1M8ClvfpwcXE3EHkIyMxj1oGjG+BXHShTA/AsC7h66QNB0lXg6in976YmFK7dDs6ARwBgZ69vd2p80Z6/vaM+fDi56cPIreuArpDVFwlTEoAC6un3meyHpGv/BRPZ8mpf6cpA+eZAhRb6TQKUo3P222SkAwmR+hAom9x3apK+zRJ01GPXvfN1Ipsjy+t3m7INCSnpeLltFXzweF2tm0RULMdvBhVLJF0ZN87qqxhJt/ukDQcuCRpy8HX1lo9auQSMffqQIZuEBKlc5EfChdyvBA4JFHGXgeSbd97OxUtf+XAupa+SJMboD9KZaXd/Pk4eQPmm+k0qJxKQ5Pmk3dJXV6TdV6XCcRK4djL3NkvlRiocgVLpqK8/eHuV1YcaVdG4BEQfA6KO3K66HNFff8fz8Ab8pXJTS/9VApabjz4Uyb6Wfa4qN8f0laC0JBSIm6++wiTBJiEauHn+zttI5Ueeu7x+8trJ7WR/SzUnP/IcJfBUbgdUaQ+Uawq4lMq98qXCU4T+NZTvXTz1VauAOvrLZNHWHLmCV3/apy6/27U2BnaopnWTiPLEoGLu5CBhOMjk1zUi5MAfeRg4vw04t03/NbewYMzeSR8w3P30B7qbEbl/spdAU66J/kAv1Kf3CP3BLK+DsIMLULbh7WrA7YqAVBOMg5GhGiAHZEPAkPCSHKc/WMvjSveHHCBlu1v3Sda+yARib1cW5K0r9yNhKrcDc37kdyUESBVFKhTy+7Kv5KCd83ncLTCqAHQMiL+ifz0l2EjgkE32qwQUZ4/svyvdVZf3/VeVkk1VtfKo/HiV11dgpH0S4qTNCVH6wJXb73kH3w6YrvrXPSFG/9rmR0JZpTZApbb6rz4VC7APMvTdWfI4Bd1vZFLTNp/GF2vC1GVZxfbFlsFaN4nI8oPKlClT8OWXXyIyMhKNGjXCd999h5YtW1pPUJFP4ue26rcL24EbF4CU25/m5ZO6VEFkDIUccD0D9d0a0gUhYUTGSFw++N/tDaTCoQ5c/voxIRIqJAgkXsu7qyErYDT7bytdJfegJG8LCQM3zgPJN/TVDamWSFvV2I4cXRR0b2R/S/VIgo+EOQmXqhuugv69kVeQM1TKDO8v2eIv51+9kgAiFR7pLpIqmVRZJGTlJOFTwqiMy3FyBzKky+4aEB+lDz3yVbrw5L0q72NPGb9TG6jQEqjSQf/+4vtEExP+Oo4ZW86oy5891QA9WjGskPmxmKCyZMkS9O7dG9OnT0erVq0wadIkLFu2DOHh4QgICNA+qMinZul2kHEcUvaXMCB/mOUTpJToJShIsJBPunJgkU+VEkxiwvW/IwcO+dQr4z/uhbMnUClE/2lXSvxlGwEOTnlXa+QAooLLVf1BRA54Ekp44LB+ElalgiVVGxncLO8dd1/9IGT5mlvlQ37n4i59te78duDygbt3Od2NhJuKrf7rjpLgI9UgGccj700VeGSLBhKj//t3JbeTQdMSdNxK31sbbJT8Wf/4z+OYs+2s+v7jJ+uhV4h+2X0ic2ExQUXCSYsWLfD999+r7zMzM1GxYkUMGzYMo0ePvuP2KSkpajN+onL7Yg8q4X8Bh38BzmzSf4q8V/61gcrtgcpt9QMvVReDvb4SIn+s5dO0bHJZ/ohLuV+6C2RWiXSNBDbIPtOFyJRkWnfELv14HnlPygBoqcBIt1apIMAz6HYVMOi/wcgS5CWUS9g5+68+LN8r/zpAcCugYmsguLW+iphfF5N0RcnYIwlD0pV3t25VKyZ/2j9bfRyz/tWHlQ+fqIc+bRhWbJZOp/93Lf8+5EOIfICRf8/5VWvl3/SZzcDpDUDIUP2/QVsLKqmpqXB3d8cvv/yCbt26ZV3fp08f3Lx5EytXrrzjd8aPH48PP/zwjuuLPaisex/Y/t1/nwzLNtbPsJCqiYQICRgSYAyfBmUMgHyV8CHjJnyrAUH1geAQfQVEggmRrZA/KzJF/Ny/+k3GKcnYJxnPoro7b48HktCjvgbox1WlGgLPUeD66TvvV8KRGizsqa/AyPgc400NkL79J03GC8m/WfUBob1+3I38YbYh8uf98zVhWd1AYx+ri/7tqmjdLCoJ6an6cXBntug/PMiHiJwfumUAvww5kHF0avybdPGm6T8sy9g749u3GwF0Gm97QeXy5csoX748tm/fjpCQkKzr3377bWzZsgU7d+7UrqIigxxPrAGqPagfMJpXNwsRmYZ8CIjYCVzYof8q3VEFmUWmSNUlx582CUhBDfUfHGSGmHQxybir1MTbHzaib1c2r+i7byVUyR9u+QNu+NAhFR4Lq9LIn/gv14Zj6mZ98Huvax0M6FBV62bR3ciAeRk+cHEXcOXg7fWcMvUflmUso6GyKZt0mapxblf0AUOGHEhFVK1rlYMEeAkouf0sJ/lQLmPOqnUEanfV/7vRIKhYXF+Ci4uL2kxOzWZpbvrHIaLcSUCo87h+EzKgW/5wS7VFpqhL8JAxLFmb7+2vt6smcZf00/DP3q7qyMB0WVRQtqKSaql0Q8lgYQku0i3r7J59TJsMNJZB6LLujZpVlqAPPNI2mfYuA45zzgIz8YJwb3WpBUd7O0zeeAqfrj6ODJ0Or97HqcsmobpMzupnAspAdam+y3tVxjSqRTBvb8YzFQ2zGeUDsgQMCScy0zOv9akKSh5TqolqcLss/1Dzv/eezMqUgC6zH6WdUkVR/64cbi9KKcMO6gFOrtCapkHFz88PDg4OiIqKyna9fB8UFKRZu4jIDMlYGDV9uk3Bbi/jWWSr/4z+e/mDbPikKWvgyKdPqaYYqiZquz37TrqmZJZcWqJ+DZ5LcgDZre9aOrlWvyl2+gAiBx0ZxC7dT3er+siBQCo0atzN7fE3cgDLSUrwhnVvpNtMyvDyCVfCkjwvCUkyOLqAYWVk51qwt7fDpL9P4vO/wlTNaRDDSvaAIStTRx3WT0IwvDdkHxsWwVQz8HJU1OSAL90qUvmTWZ3yVcZ13Y2s2SQLYEqVRKp5ub1v5L1YQRaEbK4PD8K48qcGo8sYstTbtw/UvzcMQw5kbGReY7pkzKPM6JPNzD+Um8VgWpmKLFOSDYNpg4ODMXTo0FwH01rk9GQisnxyQJJFEtVg4X/01R35FJrXwn1yYJODi3yClU+qcoCRcTvyKTYntTSBDEx2vz3AXg5GV+6+UrKsRWQ4kMlXWejQ0E1tWGJADr4yhf32dv7SJVy8Goc4uKNOjZqoXLeFfi0ltcLxXWYFyn3KgVGqU3Iwl7aqxSb99b9fwOCkSFeG3I8ESDndhfok76d/TlJxuFsXm7RFQpysQRQdph/TJEFSBn5LmJM2ySKUMtvRt4r+qwSDrOcRrQ8lWWsZ7bn7+lTSbSJtk9dK7kMqZlI9k6n7OW8n+0PeA7IYps6wCOal2883l8UmJYRK12TFlvrXUr5KQLKz3vWJLGKMimF6sgyenTFjhgosMj156dKlCAsLQ2Bg4F1/n0GFiDQjg+hldpMcdOWAbZgVld8MPamSqLE3O4GIHUDkkbyngktVR9a9kQOWVHwkuKjwcUrftXSH2xUemckhFYGCrpysHstZPwZBQosc1KXbQNolVSLDua9k2QXDObRyI2vpqPNpGVaHLqc/UEt1KO6ivjsk+vYmASqvZRtkzIVMSJCqkUxVl+cu10n3nzxvWTJCgkVuoS8/at845r1v5DHkVBYSHCVgqvWpruorWtKVmNfrJM9RpuLLEhJSzZAB3/l1mcjsGwks0g6hxp3IStoWNxrDNoKKkKnJhgXfGjdujMmTJ6tKS0EwqBCRRZOQI5UWCTxyYJQBvmo8Q8XbswzzqCzI2jfqnFu79V1TF/fm8kndTl/dkYqCLHUgm4c/Mu0c8fOWQ0iIPoumThfQ3Pk87HP7lJ8bGW8ha/JIN4g6mN/Sh5Bcg9NdGBY1dPXRhwB1yoyzBR80LaFDTpkh60oZTvQps8fkuagFFM/pN7nPnNPl5XlIIDNUo+Sr3FdeEyekoiaLKcrzNKxPJbPPpKtFKicFXV2bLDOo3AsGFSKi2wMy1UKPt1cLlkG+UuHI45N9fHIanvx+G85cTUS7amUw/+lAOFzZr59hItUDNSbGSX8wlqqGDMSUMQ9lquV+EkxZUVkCl3SnqHNhhesHFkv1QG4v4UsGZxpOmyGVi9yWbZBAIANLpYoTc1wfMgxjOKTCJFUaaYNUfiSgSCWrIKR9Us2SKo78jpzeggtgaopBhYiI8hUeGa/OuHwrLQPDOlbHqM61tG4S2ZC4Ah6/LWtBACIiKja1gjzViQvFdxtPYWNYAWarEJUwBhUiIhvWrUl59A6ppC4PX3wAEdcLMQCXqAQwqBAR2bj3Hq2DxhV9EJecjld/2ovktHs8ISVRMWJQISKycS6ODpjasylKuzvh6OU4jFt5VOsmEWVhUCEiIpTzccPkF5uo9cWW7InAkt1FmG5MZAIMKkREpLSv4Y+RnWqqy2NXHsWRSwVcW4XIhBhUiIgoy5AHqqNj7QCkpmei37zdHFxLmmNQISKiLHLiwm+6N0atQE/ExKegz5xduJqQ41w2RCWIQYWIiLLxdnPC/JdboryPm1q59oWZOxAdl6x1s8hGMagQEdEdgrxdsfCVVijr7YpT0QnoPnMHrsTe0rpZZIMYVIiIKFeV/TywdFCIqqycvZqI7jN24OINjlmhksWgQkREearo644lg1oj2NcdF64nqbDCAbZUkhhUiIgoXxVK68NKFT8PXLp5Cy/OYmWFSg6DChER3VVZbzcsGtAalcu44+KNW2qALcMKlQQGFSIiKvAA28UDQ7LCilRWLt/kAFsyLQYVIiIqVFhZNFBfWYm4fgv95u5GXHKa1s0iK8agQkREhe4GWjigNfw9XRAeFY8hC/chPSNT62aRlWJQISKiQpMpy3P6tICbkwP+PXkVX64L17pJZKUYVIiIqEgaVPDG/55vpC7P2HIG645Gat0kskIMKkREVGRdG5TFy22rqMujlh3E+WuJWjeJrAyDChER3ZPRj9RG02AfxCen47Wf9iE5LUPrJpEVYVAhIqJ74uxojyk9m8LXwxnHrsRh3MqjWjeJrAiDChERFctMoG9faAw7O2DJnggs3R2hdZPISjCoEBFRsWhfwx8jO9VUl8euPIIjl2K1bhJZAQYVIiIqNkMeqI6OtQOQkp6JwQv3ITaJi8HRvWFQISKiYmNvb4dvnm+MCqXd1NmWRy49gMxMndbNIgvGoEJERMXK290J019qpgbZbgiLxsS1XAyOio5BhYiIil398t74/OkG6vL0LaexIPSc1k0iC8WgQkREJvF00wp4s7N+cO24349yJhAVCYMKERGZdHBt3zaVodMBb/96CD/vvKB1k8jCMKgQEZHJ2NnZYdzjdbOW2X/3t8P4ccd5rZtFFoRBhYiITB5Wxj5WBwPa68PK2BVHMH87x6xQwTCoEBFRiYSVd7vWwaD7qmaNWZm99azWzSILwKBCREQlFlZGP1wbg++vpr7/+M9jmPXPGa2bRWaOQYWIiEo0rLzVpRZe71hdff/p6uNq+jJRXhhUiIioxMPKyM61MLxTDfX953+FcTYQ5YlBhYiINDG8U00Mu11ZeX/FYWw4HqV1k8gMMagQEZFmRj5UE881qwA5HdDri/bjVHS81k0iM8OgQkREmnYDffZ0A4RULYPE1Ay8+tM+JKaka90sMiMMKkREpCknB3tMfrEJAr1ccCo6Qa1gq5OlbIkYVIiIyBz4e7pgas+mcLS3w6pDVzBnGxeEIz0GFSIiMgvNKvni/UfrqMufrT6OXWeva90kMgMMKkREZDb6tKmMJxuXQ0amDkN/3oeY+BStm0QaY1AhIiLzGlz7VAPUCCiF6PgUNRNIQgvZLgYVIiIyKx4ujpj2UlO4Ozsg9Mw1fL0+XOsmkYYYVIiIyOxUD/DEF880VJenbDqNjWFcDM5WmSyofPrpp2jTpg3c3d3h4+OT620uXLiARx99VN0mICAAb731FtLTOX+eiIiAxxuVQ982ldXlNxYfwMkoLgZni0wWVFJTU/Hcc8/htddey/XnGRkZKqTI7bZv34758+dj3rx5+OCDD0zVJCIisjDvdq2D5pVKIz45HX3n7kZ0fLLWTaISZqcz8ao6Ej6GDx+OmzdvZrv+r7/+wmOPPYbLly8jMDBQXTd9+nS88847iImJgbOzc4HuPy4uDt7e3oiNjYWXl5dJngMREWnnemIqnp66DeeuJaFOWS8sfKUVfD0Kdowg81XQ47dmY1RCQ0PRoEGDrJAiunTpohp+9OjRPH8vJSVF3cZ4IyIi6yWhZG6/lvAr5YLjV+LQY9YOFV7INmgWVCIjI7OFFGH4Xn6WlwkTJqgEZtgqVqxo8rYSEZG2qvh5YPHAViqshEXGM6zYkEIFldGjR6s57vltYWFhpmstgDFjxqgykWGLiIgw6eMREZH5zATKGVauJXBBOGvnWJgbjxo1Cn379s33NlWrVi3QfQUFBWHXrl3ZrouKisr6WV5cXFzURkREthpWWuOFmTtUWJGvCwe0QoCnq9ZNI3MIKv7+/morDiEhIWoKc3R0tJqaLNavX68G1NStW7dYHoOIiKxP9YBSWDKotaqonIxOwAsz9GGlrLeb1k0jSxqjImukHDhwQH2VqchyWbaEhAT1886dO6tA0qtXLxw8eBBr167F+++/jyFDhrBiQkRE+armXwpLB4WgvI8bzlxNxPMzQhFxPUnrZpElTU+WLiJZGyWnTZs24f7771eXz58/r9ZZ2bx5Mzw8PNCnTx98/vnncHQseKGH05OJiGzXxRtJ6PnDTpy/loSKvm749dU2CPBiN5AlKOjx2+TrqJgagwoRkW2LjE1WFZUL1/XrrEi3kJerk9bNIktfR4WIiKg4BHm74sf+/62zMuxnnnHZmjCoEBGRxatUxgPz+rWAq5M9tpyIwTfrT2jdJComDCpERGQV6pf3xudP68+4/P2mU1h7NO/FQ8lyMKgQEZHV6NakPPq11Z9xedTSgzgVrZ9pSpaLQYWIiKzujMstq/giISUdg37co76S5WJQISIiq+LkYI8pPZoiyMsVp2MS8c6vh2DhE1xtGoMKERFZHX9PF0zp2QSO9nZYdegKFoSe17pJVEQMKkREZJWaVfLFmK511OVPVh3D/gs3tG4SFQGDChERWa2X21bGI/WDkJahw9Cf9+NGYqrWTaJCYlAhIiKrZWdnh4nPNkQVPw9cunkLQxftQ1pGptbNokJgUCEiIqvm6eqEqT2bwt3ZAdtOXcPYFUc4uNaCMKgQEZHVk3MAffdiE9jbAYt3R3DlWgvCoEJERDbhwTqBGP9EPXV58sZTmPQ3w4olYFAhIiKb0TukMt67PRNo0t8nMXnDSa2bRHfBoEJERDZlQIeqGPNIbXX56/Un8P1GhhVzxqBCREQ2Z9B91fDOw/qw8tW6E/hqbTgH2JopBhUiIrJJr91fLauyImdb/ujPYwwrZohBhYiIbLqy8vGT+gG2c7edw5jlh5GRybBiThhUiIjIpvUKqYyvnmuUNXX53eWHWVkxIwwqRERk855tVgHfvqBfZ2XJngh8tvo4w4qZYFAhIiIC8Hijcvj8mYbq8qx/z2Le9nNaN4kYVIiIiP7zfPOKWQNsP/7zGDaFR2vdJJvHoEJERGRkYIeqeL55BciY2td/3o8zMQlaN8mmMagQERHlOOPyJ90aoGVlX8SnpGPQj3uRmJKudbNsFoMKERFRDs6O9vi+ZxMEeLrgZHQC3v71EAfXaoRBhYiIKBcBnq6Y9lJTONrbYdWhK5i99azWTbJJDCpERER5aFbJF2Mfq6suT/grDKGnr2ndJJvDoEJERJSP3iGV8FST8mrF2mGL9uFK7C2tm2RTGFSIiIjuMrj2s6caoE5ZL1xNSMXghfuQkp6hdbNsBoMKERHRXbg5O2D6S03h5eqI/RduqjVWqGQwqBARERVApTIeapl9Ozvgpx0X8Mvei1o3ySYwqBARERXQA7UD8MaDNdTld387jB1nOLjW1BhUiIiICuH1jjXwUN1ApKZn4pX5e3D4YqzWTbJqDCpERESFYG9vh+9ebILWVX2RkJKOXnN24sglhhVTYVAhIiIqJFcnB8zq3RyNK/rgZlIaXpy1A/sv3NC6WVaJQYWIiKgIPF2d8GP/lmhRuTTik9PRa/Yu7Dl3XetmWR0GFSIionsIK/P6tURI1TKqG6j3nF1cvbaYMagQERHdAw8XR8zp2wLta/ghKTUD/ebtwtaTV7VultVgUCEiIiqGBeFkzErH2gFITsvEy/N3Y1NYtNbNsgoMKkRERMU0wHb6S83Q+fbU5YE/7sG6o5FaN8viMagQEREVE2dHe0zp2RSPNiiLtAwdhv68H7vOcoDtvWBQISIiKkZODvb49oXG6FIvEKkZ+srKmZgErZtlsRhUiIiIipmjgz0mdW+CRrfXWRmwYI+aFUSFx6BCRERkogG2P/RujiAvV5yOScRbyw5Cp9Np3SyLw6BCRERkIv6eLpj6UlM4OdjhryORmPHPGa2bZHEYVIiIiEyoaXBpjHu8nro8cU0Ytp3iGiuFwaBCRERkYj1bBePZZhWQqQOGLdqPSzdvad0ki8GgQkREZGJ2dnb4pFt91C/vheuJqRj8014kp2Vo3SzbDirnzp1D//79UaVKFbi5uaFatWoYN24cUlNTs93u0KFDaN++PVxdXVGxYkVMnDjRVE0iIiLSdEG4aT2bwcfdCQcvxuLDP45q3STbDiphYWHIzMzEjBkzcPToUXzzzTeYPn063n333azbxMXFoXPnzqhUqRL27t2LL7/8EuPHj8fMmTNN1SwiIiLNVPR1x+QXmsDODli0KwJLdl/Quklmz05XgnOlJIhMmzYNZ87oRz3L5ffeew+RkZFwdnZW140ePRorVqxQQacgJOx4e3sjNjYWXl5eJm0/ERFRcfh+40l8te6EWsn2l1dD0LCCD2xNXAGP3yU6RkUa4+vrm/V9aGgoOnTokBVSRJcuXRAeHo4bN27keh8pKSnqyRlvRERElmTw/dXRqY7+nECvzN+DiOtJWjfJbJVYUDl16hS+++47DBo0KOs6qaQEBgZmu53he/lZbiZMmKASmGGTcS1ERESWxN7eDl93b4RagZ6Ijk9Br9k7EROfonWzrCOoSNeMjF7Ob8vZbXPp0iU8/PDDeO655zBgwIB7avCYMWNUZcawRURE3NP9ERERacHL1QkL+rdEhdJuOHctCS/MDEVkbLLWzTI7joX9hVGjRqFv37753qZq1apZly9fvowHHngAbdq0uWOQbFBQEKKiorJdZ/hefpYbFxcXtREREVm6QC9X/NS/FXrM2qGW2X9+RigWvtJKDbqlIgYVf39/tRWEVFIkpDRr1gxz586FvX32Ak5ISIgaTJuWlgYnJyd13fr161GrVi2ULl26sE0jIiKyOJX9PLD01RD0/GEnzl9LUmHlp1daoZp/Ka2bZt1jVCSk3H///QgODsZXX32FmJgYNe7EeOxJjx491EBaWW9FpjAvWbIE3377LUaOHGmqZhEREZmdCqXdsXRQCKoHlMKV2GR0nxGK8Mh4rZtl3dOT582bh379+uX6M+OHlAXfhgwZgt27d8PPzw/Dhg3DO++8U+DH4fRkIiKyFtcSZGDtLhy7EgdfD2fVLVS3nHUe2wp6/C7RdVRMgUGFiIisyc2kVPSeswuHLsbC281JhZUGFbxhbcxyHRUiIiLKn4+7sxqj0jTYB7G30tDjhx04EHETtopBhYiIyCynLrdCy8q+iE9OR7+5u3AmJgG2iEGFiIjIDJVyccTcfi3QsII3biSloe/c3biaYHuLwjGoEBERmSkPF0fM7tMCwb7uuHA9CUN/3of0jEzYEgYVIiIiM+bv6YI5fZvDw9kBO85cx8S14bAlDCpERERmrnqAJ758rpG6PPOfM1h16ApsBYMKERGRBejaoCwGddCfouatXw7iZJRtLAjHoEJERGQh3upSCyFVyyApNQODftyLuOQ0WDsGFSIiIgvh6GCP73o0QVlvV5y5mohRSw8iM9Oi1229KwYVIiIiC+JXygXTXmoGZwd7rD8WhWlbTsOaMagQERFZmMYVffBxt3rq8lfrwrHlRAysFYMKERGRBereIhgvtgyGnLHv9UX7cf5aIqwRgwoREZGFGv9EXVVdkXMC9ZmzyypXrmVQISIislAujg6Y2asZKpR2w7lrSXh53m6rmwnEoEJERGTBArxcMf/llijt7oRDF2PR64eduJmUCmvBoEJERGThqvmXwo/9W6mwcvBiLF6YucNquoEYVIiIiKxA/fLeWDwwRE1fDouMV2ElKi4Zlo5BhYiIyErUCvLE0kGt1YJwp6IT8PyMUFy8kQRLxqBCRERkRar6l8LSQSGo6OuG89eS8Pz0UJy9arlTlxlUiIiIrExFX3csG9QG1fw9cDk2Gc9ND0V4pGWexJBBhYiIyAoFebtiyaAQ1A7yVANrX5gZirDIOFgaBhUiIiIr5VfKBYsHtkajij64kaRfFC7iumWNWWFQISIismI+7s5Y0K8lagaWQlRcigorspKtpWBQISIisnLe7k5Y8HIrlPN2xZmriRi+eD8yM3WwBAwqRERENjJmZWbv5nBxtMem8BhM+vsELAGDChERkQ0tCjfh6Qbq8uSNp7D2aCTMHYMKERGRDXm6aQX0a1tZXR619KBaGM6cMagQERHZmHe71kGrKr5ISEnHoB/3qK/mikGFiIjIxjg52GNKz6YI8nLF6ZhEvLXsIHQ68xxcy6BCRERko2usTH2pKZwc7PDXkUjM+OcMzBGDChERkY1qGlwa4x6vpy5PXBOGbaeuwtwwqBAREdmwnq2C8WyzCpBlVYYt2o8L18xr5VoGFSIiIhtmZ2eHT7rVR/3yXriemIo+c3epr+aCQYWIiMjGuTo5YE6fFijv44azVxPRb95uxCWbxzL7DCpERESEAC9XzH+5BXzcnXAw4iZe+mEnbiZpX1lhUCEiIiKleoAnfn6lNXw9nHHoYix6zNqpeTcQgwoRERFlqVvOC4sHtlbTl49dicMLM0MRE58CrTCoEBERUTY1Az2xZFBrBHq54ERUAj7+8xi0wqBCREREd6jmXwpLBoagU50AfPSkfq0VLThq9shERERk1ir7eeCHPi00bQMrKkRERGS2GFSIiIjIbDGoEBERkdliUCEiIiKzxaBCREREZotBhYiIiMwWgwoRERGZLQYVIiIiss2g8sQTTyA4OBiurq4oW7YsevXqhcuXL2e7zaFDh9C+fXt1m4oVK2LixImmbBIRERFZEJMGlQceeABLly5FeHg4fv31V5w+fRrPPvts1s/j4uLQuXNnVKpUCXv37sWXX36J8ePHY+bMmaZsFhEREVkIO51OpyupB/v999/RrVs3pKSkwMnJCdOmTcN7772HyMhIODs7q9uMHj0aK1asQFhYWK73Ib8rm3HYkUpMbGwsvLy8SuqpEBER0T2Q47e3t/ddj98lNkbl+vXrWLhwIdq0aaNCiggNDUWHDh2yQoro0qWLqsDcuHEj1/uZMGGCemKGTUIKERERWSeTB5V33nkHHh4eKFOmDC5cuICVK1dm/UwqKYGBgdlub/hefpabMWPGqPRl2CIiIkz8DIiIiMhizp4sXTNffPFFvrc5fvw4ateurS6/9dZb6N+/P86fP48PP/wQvXv3xp9//gk7O7siNdjFxUVtBoaeKykhERERkWUwHLfvNgKl0EFl1KhR6Nu3b763qVq1atZlPz8/tdWsWRN16tRRXTU7duxASEgIgoKCEBUVle13Dd/LzwoiPj5efWUXEBERkeWR47gM5Si2oOLv76+2osjMzFRfDYNhJazIYNq0tLSscSvr169HrVq1ULp06QLdZ7ly5VT3j6enZ5GrNHkxDNSV++dAXdPhfi4Z3M8lg/u5ZHA/W/5+lkqKhBQ5juen0EGloHbu3Indu3ejXbt2KnTI1OSxY8eiWrVqKqCIHj16qO4g6RqSsSxHjhzBt99+i2+++abAj2Nvb48KFSrAlOTF4T8E0+N+LhnczyWD+7lkcD9b9n7Or5Ji8sG07u7uWL58OR588EFVIZEw0rBhQ2zZsiVrjIk0cN26dTh79iyaNWumupU++OADDBw40FTNIiIiIgtisopKgwYNsHHjxrveTsLLv//+a6pmEBERkQXjuX7yIZWfcePGZZtlRMWP+7lkcD+XDO7nksH9bDv7uURXpiUiIiIqDFZUiIiIyGwxqBAREZHZYlAhIiIis8WgQkRERGaLQYWIiIjMlk0HlSlTpqBy5cpwdXVFq1atsGvXrnxvv2zZMnWyRbm9rBOzevXqEmurLe3rWbNmoX379mpFY9k6dep019eGivaeNli8eLE6BUW3bt1M3kZb3M83b97EkCFDULZsWTXNU859xr8fxb+fJ02apBYYdXNzU8u+jxgxAsnJySXWXkv0zz//4PHHH1fL2MvfgBUrVtz1dzZv3oymTZuq93L16tUxb9480zZSZ6MWL16sc3Z21s2ZM0d39OhR3YABA3Q+Pj66qKioXG+/bds2nYODg27ixIm6Y8eO6d5//32dk5OT7vDhwyXedmvf1z169NBNmTJFt3//ft3x48d1ffv21Xl7e+suXrxY4m235v1scPbsWV358uV17du31z355JMl1l5b2c8pKSm65s2b67p27arbunWr2t+bN2/WHThwoMTbbs37eeHChToXFxf1Vfbx2rVrdWXLltWNGDGixNtuSVavXq177733dMuXL5elSnS//fZbvrc/c+aMzt3dXTdy5Eh1LPzuu+/UsXHNmjUma6PNBpWWLVvqhgwZkvV9RkaGrly5croJEybkevvnn39e9+ijj2a7rlWrVrpBgwaZvK22tq9zSk9P13l6eurmz59vwlba5n6WfdumTRvdDz/8oOvTpw+Dign287Rp03RVq1bVpaamlmArbW8/y207duyY7To5mLZt29bkbbUWKEBQefvtt3X16tXLdl337t11Xbp0MVm7bLLrJzU1FXv37lVdCsYnN5TvQ0NDc/0dud749qJLly553p6Kvq9zSkpKUmfY9vX1NWFLbXM/f/TRRwgICFDn4iLT7Offf/9dnYhVun4CAwNRv359fPbZZ8jIyCjBllv/fm7Tpo36HUP30JkzZ1T3WteuXUus3bYgVINjocnO9WPOrl69qv5IyB8NY/J9WFhYrr8TGRmZ6+3leirefZ2TnFlb+k9z/uOge9vPW7duxezZs3HgwIESaqVt7mc5YMp5z3r27KkOnKdOncLgwYNV+Jalyal49nOPHj3U77Vr1056CpCeno5XX30V7777bgm12jZE5nEsjIuLw61bt9T4oOJmkxUVshyff/65Guj522+/qQF1VDzi4+PRq1cvNXDZz89P6+ZYtczMTFW1mjlzpjpLfPfu3fHee+9h+vTpWjfNqsgAT6lUTZ06Ffv27cPy5cuxatUqfPzxx1o3je6RTVZU5A+zg4MDoqKisl0v3wcFBeX6O3J9YW5PRd/XBl999ZUKKn///bc6yzYV334+ffo0zp07p0b7Gx9QhaOjI8LDw1GtWrUSaLn1v59lpo+Tk5P6PYM6deqoT6bSxeHs7GzydtvCfh47dqwK36+88or6XmZmJiYmYuDAgSoYStcR3bu8joVeXl4mqaYIm3zl5A+DfLLZsGFDtj/S8r30JedGrje+vVi/fn2et6ei72sxceJE9UlozZo1aN68eQm11nb2s0yzP3z4sOr2MWxPPPEEHnjgAXVZpnZS8byf27Ztq7p7DEFQnDhxQgUYhpTi288yli1nGDGEQ557t/hocizU2fDUN5nKNm/ePDXFauDAgWrqW2RkpPp5r169dKNHj842PdnR0VH31VdfqSmz48aN4/RkE+3rzz//XE1L/OWXX3RXrlzJ2uLj4zV8Fta3n3PirB/T7OcLFy6oWWtDhw7VhYeH6/78809dQECA7pNPPtHwWVjffpa/ybKfFy1apKbQrlu3TletWjU1Y5PyJn9XZSkI2SQSfP311+ry+fPn1c9lH8u+zjk9+a233lLHQllKgtOTTUjmfwcHB6uDokyF27FjR9bP7rvvPvWH29jSpUt1NWvWVLeX6VmrVq3SoNXWv68rVaqk/sHk3OQPERXve9oYg4rp9vP27dvVcgZy4JWpyp9++qmaGk7Ft5/T0tJ048ePV+HE1dVVV7FiRd3gwYN1N27c0Kj1lmHTpk25/r017Fv5Kvs65+80btxYvS7yfp47d65J22gn/zFdvYaIiIio6GxyjAoRERFZBgYVIiIiMlsMKkRERGS2GFSIiIjIbDGoEBERkdliUCEiIiKzxaBCREREZotBhYiIiMwWgwoRERGZLQYVIiIiMlsMKkRERARz9X+2LoAnyaHalAAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# define the function to plot the solution obtained using matplotlib\n", "def plot_solution(pinn_to_use, title):\n", " pts = pinn_to_use.problem.spatial_domain.sample(256, \"grid\", variables=\"x\")\n", " predicted_output = pinn_to_use(pts).extract(\"u\").tensor.detach()\n", " true_output = pinn_to_use.problem.solution(pts).detach()\n", " plt.plot(\n", " pts.extract([\"x\"]), predicted_output, label=\"Neural Network solution\"\n", " )\n", " plt.plot(pts.extract([\"x\"]), true_output, label=\"True solution\")\n", " plt.title(title)\n", " plt.legend()\n", "\n", "\n", "# plot the solution of the two PINNs\n", "plot_solution(pinn, \"PINN solution\")\n", "plt.figure()\n", "plot_solution(sapinn, \"Self Adaptive PINN solution\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can clearly observe that neither of the two solvers has successfully learned the solution. \n", "The issue is not with the optimization strategy (i.e., the solver), but rather with the model used to solve the problem. \n", "A simple `FeedForward` network struggles to handle multiscale problems, especially when there are not enough collocation points to capture the different scales effectively.\n", "\n", "Next, let's compute the $l_2$ relative error for both the `PINN` and `SAPINN` solutions:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Relative l2 error PINN 2593.94%\n", "Relative l2 error SAPINN 1861.15%\n" ] } ], "source": [ "# l2 loss from PINA losses\n", "l2_loss = LpLoss(p=2, relative=False)\n", "\n", "# sample new test points\n", "pts = pts = problem.spatial_domain.sample(100, \"grid\")\n", "print(\n", " f\"Relative l2 error PINN {l2_loss(pinn(pts), problem.solution(pts)).item():.2%}\"\n", ")\n", "print(\n", " f\"Relative l2 error SAPINN {l2_loss(sapinn(pts), problem.solution(pts)).item():.2%}\"\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Which is indeed very high!\n", "\n", "## Fourier Feature Embedding in PINA\n", "Fourier Feature Embedding is a technique used to transform the input features, aiding the network in learning multiscale variations in the output. It was first introduced in [*On the Eigenvector Bias of Fourier Feature Networks: From Regression to Solving Multi-Scale PDEs with Physics-Informed Neural Networks*](https://doi.org/10.1016/j.cma.2021.113938), where it demonstrated excellent results for multiscale problems.\n", "\n", "The core idea behind Fourier Feature Embedding is to map the input $\\mathbf{x}$ into an embedding $\\tilde{\\mathbf{x}}$, defined as:\n", "\n", "$$\n", "\\tilde{\\mathbf{x}} = \\left[\\cos\\left( \\mathbf{B} \\mathbf{x} \\right), \\sin\\left( \\mathbf{B} \\mathbf{x} \\right)\\right],\n", "$$\n", "\n", "where $\\mathbf{B}_{ij} \\sim \\mathcal{N}(0, \\sigma^2)$. This simple operation allows the network to learn across multiple scales!\n", "\n", "In **PINA**, we have already implemented this feature as a `layer` called [`FourierFeatureEmbedding`](https://mathlab.github.io/PINA/_rst/layers/fourier_embedding.html). Below, we will build the *Multi-scale Fourier Feature Architecture*. In this architecture, multiple Fourier feature embeddings (initialized with different $\\sigma$ values) are applied to the input coordinates. These embeddings are then passed through the same fully-connected neural network, and the outputs are concatenated with a final linear layer.\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "class MultiscaleFourierNet(torch.nn.Module):\n", " def __init__(self):\n", " super().__init__()\n", " self.embedding1 = FourierFeatureEmbedding(\n", " input_dimension=1, output_dimension=100, sigma=1\n", " )\n", " self.embedding2 = FourierFeatureEmbedding(\n", " input_dimension=1, output_dimension=100, sigma=10\n", " )\n", " self.layers = FeedForward(\n", " input_dimensions=100, output_dimensions=100, layers=[100]\n", " )\n", " self.final_layer = torch.nn.Linear(2 * 100, 1)\n", "\n", " def forward(self, x):\n", " e1 = self.layers(self.embedding1(x))\n", " e2 = self.layers(self.embedding2(x))\n", " return self.final_layer(torch.cat([e1, e2], dim=-1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will train the `MultiscaleFourierNet` using the `PINN` solver. \n", "Feel free to experiment with other PINN variants as well, such as `SAPINN`, `GPINN`, `CompetitivePINN`, and others, to see how they perform on this multiscale problem." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "GPU available: True (mps), used: False\n", "TPU available: False, using: 0 TPU cores\n", "HPU available: False, using: 0 HPUs\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "406bab5254d34c1dac9234243b2f3eb4", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Training: | | 0/? [00:00" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# plot solution obtained\n", "plot_solution(multiscale_pinn, \"Multiscale PINN solution\")\n", "\n", "# sample new test points\n", "pts = pts = problem.spatial_domain.sample(100, \"grid\")\n", "print(\n", " f\"Relative l2 error PINN with MultiscaleFourierNet: {l2_loss(multiscale_pinn(pts), problem.solution(pts)).item():.2%}\"\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is clear that the network has learned the correct solution, with a very low error. Of course, longer training and a more expressive neural network could further improve the results!\n", "\n", "## What's Next?\n", "\n", "Congratulations on completing the one-dimensional Poisson tutorial of **PINA** using `FourierFeatureEmbedding`! There are many potential next steps you can explore:\n", "\n", "1. **Train the network longer or with different layer sizes**: Experiment with different configurations to improve accuracy.\n", "\n", "2. **Understand the role of `sigma` in `FourierFeatureEmbedding`**: The original paper provides insightful details on the impact of `sigma`. It's a good next step to dive deeper into its effect.\n", "\n", "3. **Implement the *Spatio-temporal Multi-scale Fourier Feature Architecture***: Code this architecture for a more complex, time-dependent PDE (refer to Section 3 of the original paper).\n", "\n", "4. **...and many more!**: There are countless directions to further explore, from testing on different problems to refining the model architecture.\n", "\n", "For more resources and tutorials, check out the [PINA Documentation](https://mathlab.github.io/PINA/)." ] } ], "metadata": { "kernelspec": { "display_name": "pina", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.21" } }, "nbformat": 4, "nbformat_minor": 2 }