a9b1bd2826
Tutorial update and small fixes * Tutorials update + Tutorial FNO * Create a metric tracker callback * Update PINN for logging * Update plotter for plotting * Small fix LabelTensor * Small fix FNO --------- Co-authored-by: Dario Coscia <dariocoscia@cli-10-110-13-250.WIFIeduroamSTUD.units.it> Co-authored-by: Dario Coscia <dariocoscia@dhcp-176.eduroam.sissa.it>
396 lines
31 KiB
Plaintext
Vendored
396 lines
31 KiB
Plaintext
Vendored
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Tutorial 5: Fourier Neural Operator Learning"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "8762bbe5",
|
|
"metadata": {},
|
|
"source": [
|
|
"In this tutorial we are going to solve the Darcy flow 2d problem, presented in [Fourier Neural Operator for\n",
|
|
"Parametric Partial Differential Equation](https://openreview.net/pdf?id=c8P9NQVtmnO). First of all we import the modules needed for the tutorial. Importing `scipy` is needed for input output operation, run `pip install scipy` for installing it."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 29,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"\n",
|
|
"from scipy import io\n",
|
|
"import torch\n",
|
|
"from pina.model import FNO, FeedForward # let's import some models\n",
|
|
"from pina import Condition\n",
|
|
"from pina import LabelTensor\n",
|
|
"from pina.solvers import SupervisedSolver\n",
|
|
"from pina.trainer import Trainer\n",
|
|
"from pina.problem import AbstractProblem\n",
|
|
"import matplotlib.pyplot as plt"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Data Generation\n",
|
|
"\n",
|
|
"We will focus on solving the a specfic PDE, the **Darcy Flow** equation. The Darcy PDE is a second order, elliptic PDE with the following form:\n",
|
|
"\n",
|
|
"$$\n",
|
|
"-\\nabla\\cdot(k(x, y)\\nabla u(x, y)) = f(x) \\quad (x, y) \\in D.\n",
|
|
"$$\n",
|
|
"\n",
|
|
"Specifically, $u$ is the flow pressure, $k$ is the permeability field and $f$ is the forcing function. The Darcy flow can parameterize a variety of systems including flow through porous media, elastic materials and heat conduction. Here you will define the domain as a 2D unit square Dirichlet boundary conditions. The dataset is taken from the authors original reference.\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 36,
|
|
"id": "2ffb8a4c",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# download the dataset\n",
|
|
"data = io.loadmat(\"Data_Darcy.mat\")\n",
|
|
"\n",
|
|
"# extract data\n",
|
|
"k_train = torch.tensor(data['k_train'], dtype=torch.float).unsqueeze(-1)\n",
|
|
"u_train = torch.tensor(data['u_train'], dtype=torch.float).unsqueeze(-1)\n",
|
|
"k_test = torch.tensor(data['k_test'], dtype=torch.float).unsqueeze(-1)\n",
|
|
"u_test= torch.tensor(data['u_test'], dtype=torch.float).unsqueeze(-1)\n",
|
|
"x = torch.tensor(data['x'], dtype=torch.float)[0]\n",
|
|
"y = torch.tensor(data['y'], dtype=torch.float)[0]"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "9a9defd4",
|
|
"metadata": {},
|
|
"source": [
|
|
"Let's visualize some data"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 88,
|
|
"id": "c8501b6f",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAEjCAYAAAARyVqhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA73klEQVR4nO3de3RTZbo/8O9O2qaFtilI6QW5FbkI2KIonSpXQUv1IKDDpaPSMoIzHhid1YWjdQkUbz3iURmBA+oMFAdUdAbhzBEZsQocB1AB63Xk0NrSIm1pC73T5rLf3x/+Gg295H1pQ3bD97NWFmTnyc67k+ynT5K9n1cTQggQERERGZjJ1wMgIiIi8oQFCxERERkeCxYiIiIyPBYsREREZHgsWIiIiMjwWLAQERGR4bFgISIiIsNjwUJERESGx4KFiIiIDI8FC3VbWVlZ0DQNlZWVHmMHDRqE9PR01/V9+/ZB0zTs27fPtSw9PR2DBg3q+oESGcxnn32GG2+8ET179oSmacjLy3PtTxdDdt8pKiqCpmnIycm5qMe5GDk5OdA0DUVFRV263smTJ2Py5Mlduk7qGAsWonY0NjYiKyvLragh6u7sdjvmzJmDs2fP4sUXX8Rf/vIXDBw40NfDMqRvv/0WWVlZXV7s0MUJ8PUAiC6F48ePw2TquD5/9dVXoeu663pjYyNWrVoFAPwkRX6joKAAJ0+exKuvvopFixa5lj/++ON49NFHfTgy4/n222+xatUqTJ48udU3SO+//75vBnUZY8FCHjU2NqJHjx6+HkanWCwWjzGBgYGXYCREvnXmzBkAQEREhNvygIAABATwT4KsoKAgXw/hssOfhLq5lt+dv/vuO8ydOxfh4eG44oor8NBDD6GpqcktduvWrRg7dixCQkLQu3dvzJ8/HyUlJW4xkydPxujRo3H06FFMnDgRPXr0wGOPPeb67fk///M/sX79esTFxaFHjx649dZbUVJSAiEEnnzySVx55ZUICQnBzJkzcfbs2Vbjfe+99zBhwgT07NkTYWFhuP322/HNN9+4xXz55ZdIT09HXFwcgoODER0djV//+teoqqpq8zmorKz0uO0XHsPSlp//Dl9UVITIyEgAwKpVq6BpGjRNQ1ZWFjZv3gxN0/D555+3WsczzzwDs9mMH374ocPHIvKF9PR0TJo0CQAwZ84caJrm+vawvWNYZPJGW6qrq5Geng6r1YqIiAikpaWhurpaapx2ux2rVq3C0KFDERwcjCuuuALjx4/H3r173eI+/PBDVz6JiIjAzJkz8a9//cvj+lv25Qv9PE/k5ORgzpw5AIApU6a4ckDLT8RtHcNy5swZ3HfffYiKikJwcDASEhKwZcsWt5if59JXXnkFQ4YMgcViwQ033IDPPvtM6vm5XLGc9hNz587FoEGDkJ2djcOHD+Oll17CuXPn8NprrwEAnn76aSxfvhxz587FokWLUFFRgbVr12LixIn4/PPP3T5tVVVVISUlBfPnz8c999yDqKgo123btm2DzWbD7373O5w9exarV6/G3LlzcfPNN2Pfvn145JFHkJ+fj7Vr12LZsmXYtGmT675/+ctfkJaWhuTkZDz77LNobGzEhg0bMH78eHz++eeuYmHv3r34/vvvsXDhQkRHR+Obb77BK6+8gm+++QaHDx9ulVQ9bfvFiIyMxIYNG/DAAw9g9uzZuPPOOwEA8fHxGDx4MJYsWYJt27bh2muvdbvftm3bMHnyZPTr1++iH5vIW37zm9+gX79+eOaZZ/Dggw/ihhtucNu/L6SSN35OCIGZM2fi448/xm9/+1tcffXVeOedd5CWliY1zqysLGRnZ2PRokUYN24camtrceTIERw7dgy33HILAOCDDz5ASkoK4uLikJWVhfPnz2Pt2rW46aabcOzYsU4fQD9x4kQ8+OCDeOmll/DYY4/h6quvBgDXvxc6f/48Jk+ejPz8fCxduhSDBw/G22+/jfT0dFRXV+Ohhx5yi3/99ddRV1eH3/zmN9A0DatXr8add96J77//nt/2tkdQt7Zy5UoBQNxxxx1uy//93/9dABBffPGFKCoqEmazWTz99NNuMV999ZUICAhwWz5p0iQBQGzcuNEttrCwUAAQkZGRorq62rU8MzNTABAJCQnCbre7lqempoqgoCDR1NQkhBCirq5OREREiMWLF7utt6ysTFitVrfljY2NrbbzjTfeEADEgQMHlLa9xcCBA0VaWprr+kcffSQAiI8++si1LC0tTQwcONB1vaKiQgAQK1eubDWe1NRUERsbK5xOp2vZsWPHBACxefPmVvFERtHy3n/77bfdlrfsTy1U8saF+87OnTsFALF69WrXMofDISZMmCC1jyQkJIjbb7+9w5gxY8aIvn37iqqqKteyL774QphMJrFgwQLXss2bNwsAorCw0LWsvf36wjzx9ttvt8oTLSZNmiQmTZrkur5mzRoBQGzdutW1zGaziaSkJBEaGipqa2uFED/l0iuuuEKcPXvWFbtr1y4BQPz973/vcLsvZ/xJyE8sWbLE7frvfvc7AMDu3buxY8cO6LqOuXPnorKy0nWJjo7G0KFD8dFHH7nd12KxYOHChW0+zpw5c2C1Wl3XExMTAQD33HOP2+/fiYmJsNlsrp9G9u7di+rqaqSmprqNwWw2IzEx0W0MISEhrv83NTWhsrISv/jFLwAAx44dU9p2b1mwYAFOnz7tNu5t27YhJCQEd911l9cel+hSUc0bP7d7924EBATggQcecC0zm82ufdOTiIgIfPPNNzhx4kSbt5eWliIvLw/p6eno3bu3a3l8fDxuueUWr+777dm9ezeio6ORmprqWhYYGIgHH3wQ9fX12L9/v1v8vHnz0KtXL9f1CRMmAAC+//77SzPgbog/CfmJoUOHul0fMmQITCYTioqKYDKZIIRoFdPiwq8f+/Xr1+4BZQMGDHC73lK89O/fv83l586dAwBX4rn55pvbXG94eLjr/2fPnsWqVavw5ptvug4QbFFTU9Pqvh1tu7fccsstiImJwbZt2zB16lTouo433ngDM2fORFhYmNcel+hSOXHihFLe+LmTJ08iJiYGoaGhbsuHDx8u9dhPPPEEZs6ciWHDhmH06NGYPn067r33XsTHx7vW3976rr76avzjH/9AQ0MDevbsKfV4XeHkyZMYOnRoq7MRW35CahlziwtzaUvx0pIzqTUWLH7q58d56LoOTdPw3nvvwWw2t4q9MKn8/BuOC7V1/46WCyFcYwB+PI4lOjq6VdzPv52ZO3cuDh48iIcffhhjxoxBaGgodF3H9OnT3U47bs/FNr9SYTab8atf/Qqvvvoq/uu//gv//Oc/cfr0adxzzz1ef2yiS0E1b3SliRMnoqCgALt27cL777+PP/3pT3jxxRexceNGt1Oxu5rT6fTaui/kKWdSayxY/MSJEycwePBg1/X8/Hzouo5BgwbBbDZDCIHBgwdj2LBhPhnfkCFDAAB9+/bFtGnT2o07d+4ccnNzsWrVKqxYscK1vL2vhltua2/bO8NT4bNgwQI8//zz+Pvf/4733nsPkZGRSE5O7tRjEhnFkCFDLjpvDBw4ELm5uaivr3crbI4fPy69jt69e2PhwoVYuHAh6uvrMXHiRGRlZWHRokWuRndtre+7775Dnz59Ovx2pVevXq3OWLLZbCgtLXVbpvLhZ+DAgfjyyy+h67rbtyzfffed63bqHB7D4ifWr1/vdn3t2rUAgJSUFNx5550wm81YtWpVq+pdCNHu6cJdKTk5GeHh4XjmmWdgt9tb3V5RUQHgp08dF45zzZo17a67o23vjJbeM+2dihkfH4/4+Hj86U9/wt/+9jfMnz+ffSzIb3Qmb9x2221wOBzYsGGDa5nT6XTtm55cuO7Q0FBcddVVaG5uBgDExMRgzJgx2LJli9v++fXXX+P999/Hbbfd1uH6hwwZggMHDrgte+WVV1p9w9JS9Micjn3bbbehrKwM27dvdy1zOBxYu3YtQkNDXaeT08VjdvUThYWFuOOOOzB9+nQcOnQIW7duxa9+9SskJCQAAJ566ilkZmaiqKgIs2bNQlhYGAoLC/HOO+/g/vvvx7Jly7w6vvDwcGzYsAH33nsvrrvuOsyfPx+RkZEoLi7Gu+++i5tuugnr1q1DeHg4Jk6ciNWrV8Nut6Nfv354//33UVhYeNHbfrFCQkIwcuRIbN++HcOGDUPv3r0xevRojB492hWzYMEC13PHn4PInwwZMuSi88aMGTNw00034dFHH0VRURFGjhyJHTt2tHkMWltGjhyJyZMnY+zYsejduzeOHDmCv/71r1i6dKkr5rnnnkNKSgqSkpJw3333uU5rtlqtbfZY+blFixbht7/9Le666y7ccsst+OKLL/CPf/wDffr0cYsbM2YMzGYznn32WdTU1MBiseDmm29G3759W63z/vvvx8svv4z09HQcPXoUgwYNwl//+lf885//xJo1a3hsW1fw0dlJ1EVaTkX89ttvxS9/+UsRFhYmevXqJZYuXSrOnz/vFvu3v/1NjB8/XvTs2VP07NlTjBgxQixZskQcP37cFTNp0iQxatSoVo/Tcirec88957a8vVMkW04l/Oyzz1rFJycnC6vVKoKDg8WQIUNEenq6OHLkiCvm1KlTYvbs2SIiIkJYrVYxZ84ccfr06VanIqps+8Wc1iyEEAcPHhRjx44VQUFBbZ4KWVpaKsxmsxg2bFir54zIiGRPa24hkzfa2neqqqrEvffeK8LDw4XVahX33nuv+Pzzz6VOa37qqafEuHHjREREhAgJCREjRowQTz/9tLDZbG5xH3zwgbjppptESEiICA8PFzNmzBDffvutW0xbpzU7nU7xyCOPiD59+ogePXqI5ORkkZ+f3ypPCCHEq6++KuLi4oTZbHbLGRee1iyEEOXl5WLhwoWiT58+IigoSFxzzTWttrW9XCpE+6db0480IXiET3eWlZWFVatWoaKiotWnA/K+yspKxMTEYMWKFVi+fLmvh0NE5Ld4DAtRJ+Tk5MDpdOLee+/19VCIiPwaj2Ehuggffvghvv32Wzz99NOYNWtWp89IIiKijrFgIboITzzxBA4ePIibbrpJ+swHIiK6eDyGhYiIiAyPx7AQERGR4bFgISIiIsPzi2NYdF3H6dOnERYWdknmkSGi1oQQqKurQ2xsbKsJ4IyKuYPIt1Tyhl8ULKdPn241WzAR+UZJSQmuvPJKXw9DCnMHkTHI5A2vFSzr16/Hc889h7KyMiQkJGDt2rUYN25cu/Fvv/02li9fjqKiIgwdOhTPPvusx/kgWrS0PL5y5eMwBQd7jNc8T/jrpmexwqdFhQ9ptVe3nlOnPZbS9qdy76weZ+SPuw6sVztGO6RCfhthkn/ydLPap+GGaPm3evA5+TdIwHn5WE3x+Pag6mbpWD2w7Zlf2xKQX+o56GeclZVScQ7Y8TF2d6oF+aXMG8BPuePksUEID/W8n1c4G6TXDQCVTvnccVZvf5b0C51ztj+x34VqFNYLAOcc8rMwV9nlx1FpU5vdubJZPr7yfA/p2LN18rEAYKvx/DelRUC1fJ6xVMnnsOAqxbxbKT/rdEh5o3SsufSs0jgcZeWeYxTyhlcKlu3btyMjIwMbN25EYmIi1qxZg+TkZBw/frzNORgOHjyI1NRUZGdn49/+7d/w+uuvY9asWTh27JjbvC3tafkq1xQc7JWCxWzxTsFiCpH/I2O2eK9gMQfJ7wwqsQAQECC/jUoFS4BawWIOkn+rBwQqFCwOhYJFV3zuFIoyXeF5DjAFKY1D0yTfe6Il/uJ+WrnUeePnYw0PNSE8zPN+3qRQgKjGN+sKsU7519vmVEvz5x0Kf3Tt8nkpyKb2vgsMkI8PMFmkY81O+QIEAEw2+XhTk/xzZ7bI7yfKeTdQvmAJMMvHmhVzB2Ryh0Le8MoPzS+88AIWL16MhQsXYuTIkdi4cSN69OiBTZs2tRn/xz/+EdOnT8fDDz+Mq6++Gk8++SSuu+46rFu3zhvDIyIDYt4goo50ecFis9lw9OhRTJs27acHMZkwbdo0HDp0qM37HDp0yC0eAJKTk9uNb25uRm1trduFiLqvS5E3AOYOou6sywuWyspKOJ1OREVFuS2PiopCWVlZm/cpKytTis/OzobVanVdeNAcUfd2KfIGwNxB1J11j3MPL5CZmYmamhrXpaSkxNdDIqJugLmDqPvq8oNu+/TpA7PZjPJy96ODy8vLER0d3eZ9oqOjleItFgssFvmDrIjI2C5F3gCYO4i6sy7/hiUoKAhjx45Fbm6ua5mu68jNzUVSUlKb90lKSnKLB4C9e/e2G09E/oV5g4g88cppzRkZGUhLS8P111+PcePGYc2aNWhoaMDChQsBAAsWLEC/fv2QnZ0NAHjooYcwadIkPP/887j99tvx5ptv4siRI3jllVe8MTwiMiDmDSLqiFcKlnnz5qGiogIrVqxAWVkZxowZgz179rgOkCsuLnZrwXvjjTfi9ddfx+OPP47HHnsMQ4cOxc6dO6V7KbQIOK/BJDyfyx1Qr9Ynom6IfK+N4P518iuukm9g1BTrkF8vgB7FCn1H5PsGQQ9UfO4GyJ+3L/HS/RSr0N4FUOu9UzNYfuWawkBCT8n3OwCA5gj5plxhBfLvO9EvUmkcmuSZNJowAfK97lrxVd4g6pC3ZmwwykwQKn2TfDzlhiaEYvtNA6qtrYXVakXciqelGsepFixNUd4pWBoVChboamNWKVhCS+TfAqqFglB4fxulYLGFyw9EU9h7VAsWlUSiUrBoTrVdXvyrQCrOIez4qPkt1NTUIDw8XOkxfKUld5z7vzipxnFnFDvdVig0jqtS6Eh71infBbbaqdbZ9axCp9tKu3xshU2tA3JFk8K6z8sX91W18rEA0KzS6fasQkdtlU63FWr7bI8KhU63ZV7sdPvDac8xwo592CWVN7rlWUJERER0eWHBQkRERIbHgoWIiIgMjwULERERGR4LFiIiIjI8FixERERkeCxYiIiIyPBYsBAREZHhsWAhIiIiw2PBQkRERIbnlbmEfMUW6YApxPOcO7YotTbHVw/x3F64xakaq3RseN966djmLyOkYwG1lvj1/eVbRJvUpjRSoiu029flpygCADiDFaYfUOi3rzJFgB6gNp9Aj3L5cTRFybdf71FwTmkcWoTce9qk24AzSqumS8CpkgwAOBUmuXEqfObVVebe8CKVqXN+vIPC3wuV6UVUpi1R/WpBZRuVnxDf4TcsREREZHgsWIiIiMjwWLAQERGR4bFgISIiIsNjwUJERESGx4KFiIiIDI8FCxERERlelxcs2dnZuOGGGxAWFoa+ffti1qxZOH78eIf3ycnJgaZpbpfg4OCuHhoRGRTzBhF50uUFy/79+7FkyRIcPnwYe/fuhd1ux6233oqGhoYO7xceHo7S0lLX5eTJk109NCIyKOYNIvKkyzvd7tmzx+16Tk4O+vbti6NHj2LixInt3k/TNERHR3f1cIioG2DeICJPvN6av6amBgDQu3fvDuPq6+sxcOBA6LqO6667Ds888wxGjRrVZmxzczOam5td12trawEAWpATWpDT45ii+tbIDv/Hx3PKP02Desm3Pf86b5B0rBaqNp0AIB9vssm3ZlZt4qzUbj9El44VIZ5f558zS0zZ0CIgQH7dtkb5OQIagtV2N5NN/skLPqfQTt0aojQOc5Xke1pXe0064o28AbSfO/yZSvt8ANAV+sB7s92+rpxt5GgqrfYBxXb7CtN6mBRWrPhbiMqUIUovoY/b+Hv1oFtd1/H73/8eN910E0aPHt1u3PDhw7Fp0ybs2rULW7duha7ruPHGG3Hq1Kk247Ozs2G1Wl2X/v37e2sTiOgS81beAJg7iLozrxYsS5Yswddff40333yzw7ikpCQsWLAAY8aMwaRJk7Bjxw5ERkbi5ZdfbjM+MzMTNTU1rktJSYk3hk9EPuCtvAEwdxB1Z177SWjp0qX4n//5Hxw4cABXXnml0n0DAwNx7bXXIj8/v83bLRYLLBZLVwyTiAzEm3kDYO4g6s66/BsWIQSWLl2Kd955Bx9++CEGDx6svA6n04mvvvoKMTExXT08IjIg5g0i8qTLv2FZsmQJXn/9dezatQthYWEoKysDAFitVoSE/Hiw34IFC9CvXz9kZ2cDAJ544gn84he/wFVXXYXq6mo899xzOHnyJBYtWtTVwyMiA2LeICJPurxg2bBhAwBg8uTJbss3b96M9PR0AEBxcTFMpp++3Dl37hwWL16MsrIy9OrVC2PHjsXBgwcxcuTIrh4eERkQ8wYRedLlBYsQnk/r2rdvn9v1F198ES+++GJXD4WIugnmDSLyhHMJERERkeGxYCEiIiLDY8FCREREhuf11vyXUlCIHeYennsS9w+rVlqvSaGV8ydfXSW/YotCK3qzYjtph3wLZV2lU3uA2jhMCi3xQ0ObpGMjQzueFK9VfEi9UryssoZw6diTJX2U1m2LkO+vrQcqTK9gl3/f0aWj0hLfptB7XbV9vtNLrfm91WofUJsyRLU1v6aQe5Va4qtMW6L4l1o3yz8jwqzwvYU/t+YnIiIi6gosWIiIiMjwWLAQERGR4bFgISIiIsNjwUJERESGx4KFiIiIDI8FCxERERkeCxYiIiIyPBYsREREZHgsWIiIiMjw/Ko1v9A16Lrn1sHHK/sqrbfulHz7daXW9YEKLdIltstNkPy6TUFO6djgHjalYUSGybfQHxJeKR07IrRUaRyxgdXSsQ26RTr2E3OcdGxZWJh0LABoepB0rFOhNb+pQX4KBADQAuTShKaz5f+FnApN41ViVdr424Vamrcr9Ix3KnzmdSpOEeAtyt3lTQqt+RXyv0q7fT1AbdBCqTW/wrpNbM1PRERE1CEWLERERGR4XV6wZGVlQdM0t8uIESM6vM/bb7+NESNGIDg4GNdccw12797d1cMiIgNj3iAiT7zyDcuoUaNQWlrqunz88cftxh48eBCpqam477778Pnnn2PWrFmYNWsWvv76a28MjYgMinmDiDrilYIlICAA0dHRrkufPn3ajf3jH/+I6dOn4+GHH8bVV1+NJ598Etdddx3WrVvnjaERkUExbxBRR7xSsJw4cQKxsbGIi4vD3XffjeLi4nZjDx06hGnTprktS05OxqFDh9q9T3NzM2pra90uRNS9eTtvAMwdRN1ZlxcsiYmJyMnJwZ49e7BhwwYUFhZiwoQJqKurazO+rKwMUVFRbsuioqJQVlbW7mNkZ2fDarW6Lv379+/SbSCiS+tS5A2AuYOoO+vygiUlJQVz5sxBfHw8kpOTsXv3blRXV+Ott97qssfIzMxETU2N61JSUtJl6yaiS+9S5A2AuYOoO/N647iIiAgMGzYM+fn5bd4eHR2N8vJyt2Xl5eWIjo5ud50WiwUWi3xzLyLqXryRNwDmDqLuzOt9WOrr61FQUICYmJg2b09KSkJubq7bsr179yIpKcnbQyMig2LeIKILdXnBsmzZMuzfvx9FRUU4ePAgZs+eDbPZjNTUVADAggULkJmZ6Yp/6KGHsGfPHjz//PP47rvvkJWVhSNHjmDp0qVdPTQiMijmDSLypMt/Ejp16hRSU1NRVVWFyMhIjB8/HocPH0ZkZCQAoLi4GCbTT3XSjTfeiNdffx2PP/44HnvsMQwdOhQ7d+7E6NGjlR9biB8vntSdtCqt12SXj9WDFVbskJ+zQ2WOCgCAWT4+QGEuoV49zysNY3hEueeg/y8pvEA69rrg9s8gaUukySEdW+KU/8mguPkK6ViTwpwkAKAy9YrZpvj+UKBX18jFCbV5pn7Ol3nDm3SFF1Flzh+b0nw/avO/6ArxDl1+HCrzH/0Y7515a0wmtTmvNIX9VlfIuyo5XSjOJaQyT5EIUHhdTL5tjt/lBcubb77Z4e379u1rtWzOnDmYM2dOVw+FiLoJ5g0i8oRzCREREZHhsWAhIiIiw2PBQkRERIbHgoWIiIgMjwULERERGR4LFiIiIjI8FixERERkeCxYiIiIyPBYsBAREZHheX225ksp6PNQmC2ee+MHBqqtV1eY3FU0KrRxVhiHHqTWet2peadVe89Atfbr0ZZa6di4oDPysYrv3FBTqHRso6iXjjVr8m2+nU61zwcBClNCaLpCm+9gtR1At8kNRBfy0x90V06ZuT9+xqbwmVCl3b5KG3+7wnoBwK7Sbl+hjb9qq31vteZXpdKaX+UrAIWnWanV/o/x8s+dMCs8z5pvXxN+w0JERESGx4KFiIiIDI8FCxERERkeCxYiIiIyPBYsREREZHgsWIiIiMjwWLAQERGR4XV5wTJo0CBomtbqsmTJkjbjc3JyWsUGB3vupUJE/oW5g4g60uWN4z777DM4nU7X9a+//hq33HIL5syZ0+59wsPDcfz4cdd1zcfNaYjo0mPuIKKOdHnBEhkZ6Xb9P/7jPzBkyBBMmjSp3ftomobo6OiuHgoRdSPMHUTUEa+25rfZbNi6dSsyMjI6/ORTX1+PgQMHQtd1XHfddXjmmWcwatSoduObm5vR3Nzsul5b+2P7d2H+8eKJpVqtvbYwKbSfVmiJrNLy3xGi9slRpau1PVj+bVBvD1IaR6NTPr5R4QlpFA1K44DeJB1a4QyRjq2yybf8d9jUWqQHKXS6N9nl39PaebXpFUxBcq38TUIA8k9zhy517pDl9BziRhfyv7rrCr/QOxVi7Yp93VVa+TsVEo1KG39vUh2FUrxCG3+VGRPUW/OrxKq08fftYa9effSdO3eiuroa6enp7cYMHz4cmzZtwq5du7B161bouo4bb7wRp06davc+2dnZsFqtrkv//v29MHoi8hXmDiK6kFcLlj//+c9ISUlBbGxsuzFJSUlYsGABxowZg0mTJmHHjh2IjIzEyy+/3O59MjMzUVNT47qUlJR4Y/hE5CPMHUR0Ia/9JHTy5El88MEH2LFjh9L9AgMDce211yI/P7/dGIvFAotF4fcUIuo2mDuIqC1e+4Zl8+bN6Nu3L26//Xal+zmdTnz11VeIiYnx0siIyMiYO4ioLV4pWHRdx+bNm5GWloaAAPcvcRYsWIDMzEzX9SeeeALvv/8+vv/+exw7dgz33HMPTp48iUWLFnljaERkYMwdRNQer/wk9MEHH6C4uBi//vWvW91WXFwMk+mnOuncuXNYvHgxysrK0KtXL4wdOxYHDx7EyJEjvTE0IjIw5g4iao9XCpZbb70VQrR9ete+ffvcrr/44ot48cUXvTEMIupmmDuIqD2cS4iIiIgMjwULERERGR4LFiIiIjI8FixERERkeF6dS+hSCy9yIiBQdbYPz5TmEpKbdgUA4LDIr9cUqjYDhqbL16JNZvn5fios8nPnAMC/guUnprMGnJeOtalMxAEgWLNLx37b3E86Nr+uj3Ss3qDw5gBgUpnyR2V6LJPa5xQtWK7Rmia0LptLyKgUpmwCADgVZqJReU+rzPejOoePrjA/kENXGIfKBGeK8YovixJNU1i7yiYqrVdxLjmFXVzl75vqOLoav2EhIiIiw2PBQkRERIbHgoWIiIgMjwULERERGR4LFiIiIjI8FixERERkeCxYiIiIyPBYsBAREZHhsWAhIiIiw2PBQkRERIbnV635g885EBDg8BinByi2uVfooKyybmewfL0Y0Kw2ZpNdoRZVaOPfLHoojeM7PUo6trY5WDr229AYpXGEmOVb85c3hUnHFlX2lo4NqFWbTsBsk3/jqbxHodKKGwACJadu0NVWezlwKnwm1BX6qavEOlX6tENxzAq96FVb8wvFeOn1qsZ7aRwqbfxVh6D0kitNJ8DW/EREREQdUi5YDhw4gBkzZiA2NhaapmHnzp1utwshsGLFCsTExCAkJATTpk3DiRMnPK53/fr1GDRoEIKDg5GYmIhPP/1UdWhEZFDMG0TUWcoFS0NDAxISErB+/fo2b1+9ejVeeuklbNy4EZ988gl69uyJ5ORkNDW1P43r9u3bkZGRgZUrV+LYsWNISEhAcnIyzpw5ozo8IjIg5g0i6izlgiUlJQVPPfUUZs+e3eo2IQTWrFmDxx9/HDNnzkR8fDxee+01nD59utUnqp974YUXsHjxYixcuBAjR47Exo0b0aNHD2zatEl1eERkQMwbRNRZXXoMS2FhIcrKyjBt2jTXMqvVisTERBw6dKjN+9hsNhw9etTtPiaTCdOmTWv3Ps3NzaitrXW7EFH3dKnyBsDcQdSddWnBUlZWBgCIinI/MyQqKsp124UqKyvhdDqV7pOdnQ2r1eq69O/fvwtGT0S+cKnyBsDcQdSddcuzhDIzM1FTU+O6lJSU+HpIRNQNMHcQdV9dWrBER0cDAMrLy92Wl5eXu267UJ8+fWA2m5XuY7FYEB4e7nYhou7pUuUNgLmDqDvr0oJl8ODBiI6ORm5urmtZbW0tPvnkEyQlJbV5n6CgIIwdO9btPrquIzc3t937EJH/YN4gIhnKnW7r6+uRn5/vul5YWIi8vDz07t0bAwYMwO9//3s89dRTGDp0KAYPHozly5cjNjYWs2bNct1n6tSpmD17NpYuXQoAyMjIQFpaGq6//nqMGzcOa9asQUNDAxYuXNj5LSQin2PeIKLOUi5Yjhw5gilTpriuZ2RkAADS0tKQk5ODP/zhD2hoaMD999+P6upqjB8/Hnv27EFw8E9t1wsKClBZWem6Pm/ePFRUVGDFihUoKyvDmDFjsGfPnlYH1HkSWNOMAInu53pwoNJ6Nad8z3FhVmhrfV4+1tyk2tZdJVa+3bLZpvalXHNjiHRsSbVkC3gAp3r2UhqHKVD+NdRtCs91rfwuFFyrOL2CQmt+ldbdIlDtvWQKlZuOQdPNQGXbtxk5b6hQnX1AqYW+whfeTpWW+Eq919Vb6Htrvd5q+6/aal+plb9KsOocAQqUpuroRjQhRLfftNraWlitVky59lEEmC0e4w1TsATJxzpC1P7I2EPl45vD5Xfg5l5qO3tzL/m3lz3CKR2r9fQ8Z9TPGaJgOaP2GvYolX/uQk/LPx8hP9QpjcNUd14qzqE344Oidaipqek2x4a05I5z/xeH8DDP+2OBvV5p/SUO+efhB4d8EV7hkJ/v6pyjp3QsAJyzy88XVueQn/+rxiYfCwB1dvn4epv8h536Js9/I37u/Hn5dTvr5f+2mBXmFgs6p/ZBsccZhdzxg3zu6FFwTmkczuP5HmMcwo592CWVN7rlWUJERER0eWHBQkRERIbHgoWIiIgMjwULERERGR4LFiIiIjI8FixERERkeCxYiIiIyPBYsBAREZHhsWAhIiIiw2PBQkRERIanPJeQPzDXNCnewTtza5gU5nQxNSvOJdQs/9IGnJdfd2CDWo1rqZZ/7uyh8mN2hKi9dXX57tpKc3yYFOZsClTr6q40H5TK3CHOULXW5DDJvea687JMJ4bnVJjPCFCc/0hhXh6V9f4Y7535gVTnElKaqEuXj9WUYuWHAMAwcxp1NX7DQkRERIbHgoWIiIgMjwULERERGR4LFiIiIjI8FixERERkeCxYiIiIyPCUC5YDBw5gxowZiI2NhaZp2Llzp+s2u92ORx55BNdccw169uyJ2NhYLFiwAKdPn+5wnVlZWdA0ze0yYsQI5Y0hImNi3iCizlIuWBoaGpCQkID169e3uq2xsRHHjh3D8uXLcezYMezYsQPHjx/HHXfc4XG9o0aNQmlpqevy8ccfqw6NiAyKeYOIOku501NKSgpSUlLavM1qtWLv3r1uy9atW4dx48ahuLgYAwYMaH8gAQGIjo5WHQ4RdQPMG0TUWV4/hqWmpgaapiEiIqLDuBMnTiA2NhZxcXG4++67UVxc3G5sc3Mzamtr3S5E5D+8kTcA5g6i7syrvbSbmprwyCOPIDU1FeHh4e3GJSYmIicnB8OHD0dpaSlWrVqFCRMm4Ouvv0ZYWFir+OzsbKxatarVclOjDSaZNvqaYmtm1bbIkjSbQz62Sa01v+m8/LoDGhRa89epvWWcwfLrdoQotAS3qL2GziCFFuJqT7U0k/xLAgAIbJR/49l7yj93eoBaa/4gk9xz53B0zY7irbwBtJ876OKotNvXobbPqrTmd+oK41BoiQ8AulOlNb98qEq7fZNTPhYANIV4lWk9IJkLvMVr37DY7XbMnTsXQghs2LChw9iUlBTMmTMH8fHxSE5Oxu7du1FdXY233nqrzfjMzEzU1NS4LiUlJd7YBCK6xLyZNwDmDqLuzCvfsLQknZMnT+LDDz/s8FNSWyIiIjBs2DDk5+e3ebvFYoHFojiBGxEZmrfzBsDcQdSddfk3LC1J58SJE/jggw9wxRVXKK+jvr4eBQUFiImJ6erhEZEBMW8QkSfKBUt9fT3y8vKQl5cHACgsLEReXh6Ki4tht9vxy1/+EkeOHMG2bdvgdDpRVlaGsrIy2Gw21zqmTp2KdevWua4vW7YM+/fvR1FREQ4ePIjZs2fDbDYjNTW181tIRD7HvEFEnaX8k9CRI0cwZcoU1/WMjAwAQFpaGrKysvDf//3fAIAxY8a43e+jjz7C5MmTAQAFBQWorKx03Xbq1CmkpqaiqqoKkZGRGD9+PA4fPozIyEjV4RGRATFvEFFnKRcskydPhhDtH1bc0W0tioqK3K6/+eabqsMgom6EeYOIOotzCREREZHhsWAhIiIiw2PBQkRERIbHgoWIiIgMz6ut+S81zWaHZvJcg4kQxcZRNrvCIBRaF0scaNjC1KwwBgBCpYVyg3zdag5Q61svAuXjVWJ1i9pbVw9SaPuvECsCvNeq2t5T/vlwBsqvN6BJsYW+7HtadcoLuqyotNoHAKHUml8h1qn2OV0otObXHAqxKn9WFKf1UJkGRNMVevM7vTRPjSR+w0JERESGx4KFiIiIDI8FCxERERkeCxYiIiIyPBYsREREZHgsWIiIiMjwWLAQERGR4bFgISIiIsNjwUJERESGx4KFiIiIDM+vWvM3DeiNgIBgj3GWH2rUVhwk3/dcO9+stm5ZTqdSuHZerZW/t8hMleCi0PbfFKj41lV4DUWQ/Lr1IPkx23opTgmhMHWDpUahvbYi2dbdSi2+yS+otNtXbc3v0OVzh0qsrtiaH3b5eJNdfhvVYqVDf4x3yu+LmkNhv9XZmp+IiIioQyxYiIiIyPCUC5YDBw5gxowZiI2NhaZp2Llzp9vt6enp0DTN7TJ9+nSP612/fj0GDRqE4OBgJCYm4tNPP1UdGhEZFPMGEXWWcsHS0NCAhIQErF+/vt2Y6dOno7S01HV54403Olzn9u3bkZGRgZUrV+LYsWNISEhAcnIyzpw5ozo8IjIg5g0i6izlg25TUlKQkpLSYYzFYkF0dLT0Ol944QUsXrwYCxcuBABs3LgR7777LjZt2oRHH31UdYhEZDDMG0TUWV45hmXfvn3o27cvhg8fjgceeABVVVXtxtpsNhw9ehTTpk37aVAmE6ZNm4ZDhw61eZ/m5mbU1ta6XYioe/N23gCYO4i6sy4vWKZPn47XXnsNubm5ePbZZ7F//36kpKTA2c5puZWVlXA6nYiKinJbHhUVhbKysjbvk52dDavV6rr079+/qzeDiC6hS5E3AOYOou6sy/uwzJ8/3/X/a665BvHx8RgyZAj27duHqVOndsljZGZmIiMjw3W9traWiYeoG7sUeQNg7iDqzrx+WnNcXBz69OmD/Pz8Nm/v06cPzGYzysvL3ZaXl5e3+3u2xWJBeHi424WI/Ic38gbA3EHUnXm9YDl16hSqqqoQExPT5u1BQUEYO3YscnNzXct0XUdubi6SkpK8PTwiMiDmDSK6kHLBUl9fj7y8POTl5QEACgsLkZeXh+LiYtTX1+Phhx/G4cOHUVRUhNzcXMycORNXXXUVkpOTXeuYOnUq1q1b57qekZGBV199FVu2bMG//vUvPPDAA2hoaHAd/U9E3RvzBhF1lvIxLEeOHMGUKVNc11t+D05LS8OGDRvw5ZdfYsuWLaiurkZsbCxuvfVWPPnkk7BYfppHpaCgAJWVla7r8+bNQ0VFBVasWIGysjKMGTMGe/bsaXVAnSdBVY0IMHuec8fZq4fSes3nGqVjRXCQdKzWcF5+EE61ORyEQ2HuIbtNPtab88WYFOYaMcvP4QMAWpDC66Iw7xAiQqVDzTaF9QIIqVB4zRWeOkew2nNX39/z/FwA4OhgvhMj5w0Vqp/wTJr8a2iGSqz8fmhWGIM3qc4l5FSZp0hpLiG1cUAhXvPS/EAmlfl+AJgcCrF2+feHpvJ3xQs0IRRmWDOo2tpaWK1W3Dz6YQSYPU8wp/eQ/+MFqBUsMMvvOF4tWJoVihAWLO4UChZdoWCx9w6RHwMAhb9JXi1YbFa5eIe9CUf/+jhqamq6zbEhLbnj3P/FITzM875baK9XWn+JU/79cdreSzq2wiH//FY65McAAFU2+fizNvkPf2eb1T4o1jTLFcoAUHdePvZ8g1r+1+vl80FAjfy+FVQjv9NaqtXybkiVfHyPUvkJe4NK2m810BZHUbHnGGHHPuySyhucS4iIiIgMjwULERERGR4LFiIiIjI8FixERERkeCxYiIiIyPBYsBAREZHhsWAhIiIiw2PBQkRERIbHgoWIiIgMT7k1v5GJQDOERAdUzaHWNVallb/5bIN0rFBpAV+r1mETunwLZWGX7+MsHAo9nwGvdcbVFDoKA4Boku/maAoPkx+HTf75CDoj/94AAD1Ifve0XSHf6VOVPUSuI6fTrNjyvBvy5iaqtPE3Cl2hxbJQbc2v0G7f6VRoze9Qyx2aQ6HdvkLneqVYxbSr0spf6e+hYsf1rsZvWIiIiMjwWLAQERGR4bFgISIiIsNjwUJERESGx4KFiIiIDI8FCxERERkeCxYiIiIyPOWC5cCBA5gxYwZiY2OhaRp27tzpdrumaW1ennvuuXbXmZWV1Sp+xIgRyhtDRMbEvEFEnaVcsDQ0NCAhIQHr169v8/bS0lK3y6ZNm6BpGu66664O1ztq1Ci3+3388ceqQyMig2LeIKLOUu50m5KSgpSUlHZvj46Odru+a9cuTJkyBXFxcR0PJCCg1X2JyD8wbxBRZ3m1NX95eTneffddbNmyxWPsiRMnEBsbi+DgYCQlJSE7OxsDBgxoM7a5uRnNzT+1Wq+trQUANEWGICDQc4tyy1n5Nu0AYGq0yQcHeJ4a4GJoQUFK8eJ8k/y6JaYzcK23We25Eyqt+RWmEwDUng+TwvMn7Hb5FVeclQ7VLBb59QJA7BXSoSqtuGsHqu32gefl1q3Zu2YaBm/lDaD93OEtZninlbk32/irtNvXFdrtq6wXAJy6fLzTqbBulViotebXFFKYt2IBwKSwL5pUWvML70y1IsurB91u2bIFYWFhuPPOOzuMS0xMRE5ODvbs2YMNGzagsLAQEyZMQF1dXZvx2dnZsFqtrkv//v29MXwi8gFv5Q2AuYOoO/NqwbJp0ybcfffdCA7u+FuPlJQUzJkzB/Hx8UhOTsbu3btRXV2Nt956q834zMxM1NTUuC4lJSXeGD4R+YC38gbA3EHUnXntJ6H//d//xfHjx7F9+3bl+0ZERGDYsGHIz89v83aLxQKL6tfrRGR43swbAHMHUXfmtW9Y/vznP2Ps2LFISEhQvm99fT0KCgoQExPjhZERkVExbxBRe5QLlvr6euTl5SEvLw8AUFhYiLy8PBQXF7tiamtr8fbbb2PRokVtrmPq1KlYt26d6/qyZcuwf/9+FBUV4eDBg5g9ezbMZjNSU1NVh0dEBsS8QUSdpfyT0JEjRzBlyhTX9YyMDABAWloacnJyAABvvvkmhBDtJo6CggJUVla6rp86dQqpqamoqqpCZGQkxo8fj8OHDyMyMlJ1eERkQMwbRNRZmhA+Pk+pC9TW1sJqtSLp1lU+P61ZUzlFTIFW16gUL2q8c7qm3qg4Di+d1qwFKp7WHOL5feESFKi0blmqpzU7FU5rdvSUH3P1ELVxyJ7W7LQ34ehbj6Ompgbh4eFKj+ErLbnj3P/FITzM8xfOxY56pfWXOHpIx/7g6CUdW+GQf37P2NVeiwpbmHRsZXNP6dhzzfLPBQCcbQyRjq1rkN+/7fWKuaNO/nN9YJ38KdBBNfKxlnNqf6Z7VMjn0pCy89Kx5lL5Ng4A4Dj1g+cYYcc+7JLKG5xLiIiIiAyPBQsREREZHgsWIiIiMjwWLERERGR4Xp1L6FIz2XWYhOeDXu2hagdVmkLkn6agsvbbgl9Ia1KYo0hlfhsAWi+rdKxeXiEdK5xqk1qYFA40FU75OY00s1qt7axvUIqXZQ6VP+hQ9I5QWrdml3+unQrPs8mhNAxA9ni/bn/4vm+ZvfQEmhTXqxovS2XeIQAQKvHeigXU3tdK45APVZ46yk/3RX7DQkRERIbHgoWIiIgMjwULERERGR4LFiIiIjI8FixERERkeCxYiIiIyPBYsBAREZHhsWAhIiIiw2PBQkRERIbnF51uhfixrZ/D0eyV9Zuc8m0DTU75MWi6QqdblVgA0OXHrAv5detCreOuSaHzY8vrKENT7FapC9X2rnKEwnMHhfcGAAiFpsIOh/xnD6dNrVsx7HKvi9PeBEDtdfS1lrHW1su1Eq1zqLUcbVCIb3TIvy7nnfLv52bFLtk2m3y83Sb//nc0qf25cZ6Xf0/rjfLr1c8rto1tku/A7WxS2A8V0oHTprZPORS6ZDucTdKxQlfLYQ6JvxcO/Bgjkzf8omCpq/uxHf4n//sfPh4JtSK/L6hRy8HeU+ul2G6srq4OVqv81BC+1JI7Bl5X5NuBEF3mZPKGJrrTx6F26LqO06dPIywsDJr20yfv2tpa9O/fHyUlJQgPD/fhCL3D37cP8P9t9KftE0Kgrq4OsbGxMJm6x6/NzB3cvu7KX7ZRJW/4xTcsJpMJV155Zbu3h4eHd+sX1BN/3z7A/7fRX7avu3yz0oK5g9vX3fnDNsrmje7xMYiIiIguayxYiIiIyPD8umCxWCxYuXIlLBaLr4fiFf6+fYD/b6O/b1935e+vC7ev+7sctvFCfnHQLREREfk3v/6GhYiIiPwDCxYiIiIyPBYsREREZHgsWIiIiMjwWLAQERGR4fl1wbJ+/XoMGjQIwcHBSExMxKeffurrIXWJrKwsaJrmdhkxYoSvh9UpBw4cwIwZMxAbGwtN07Bz506324UQWLFiBWJiYhASEoJp06bhxIkTvhnsRfC0fenp6a1e0+nTp/tmsJc5f80bgP/lDuaNyytv+G3Bsn37dmRkZGDlypU4duwYEhISkJycjDNnzvh6aF1i1KhRKC0tdV0+/vhjXw+pUxoaGpCQkID169e3efvq1avx0ksvYePGjfjkk0/Qs2dPJCcno6nJW7Mrdi1P2wcA06dPd3tN33jjjUs4QgL8P28A/pU7mDcus7wh/NS4cePEkiVLXNedTqeIjY0V2dnZPhxV11i5cqVISEjw9TC8BoB45513XNd1XRfR0dHiueeecy2rrq4WFotFvPHGGz4YYedcuH1CCJGWliZmzpzpk/HQT/w5bwjh37mDecP/+eU3LDabDUePHsW0adNcy0wmE6ZNm4ZDhw75cGRd58SJE4iNjUVcXBzuvvtuFBcX+3pIXlNYWIiysjK319NqtSIxMdFvXk8A2LdvH/r27Yvhw4fjgQceQFVVla+HdFm5HPIGcPnkDuYN/+OXBUtlZSWcTieioqLclkdFRaGsrMxHo+o6iYmJyMnJwZ49e7BhwwYUFhZiwoQJqKur8/XQvKLlNfPX1xP48Wvd1157Dbm5uXj22Wexf/9+pKSkwOl0+npolw1/zxvA5ZU7mDf8T4CvB0DqUlJSXP+Pj49HYmIiBg4ciLfeegv33XefD0dGF2v+/Pmu/19zzTWIj4/HkCFDsG/fPkydOtWHIyN/wtzhXy63vOGX37D06dMHZrMZ5eXlbsvLy8sRHR3to1F5T0REBIYNG4b8/HxfD8UrWl6zy+X1BIC4uDj06dPHb19TI7rc8gbg37mDecP/+GXBEhQUhLFjxyI3N9e1TNd15ObmIikpyYcj8476+noUFBQgJibG10PxisGDByM6Otrt9aytrcUnn3zil68nAJw6dQpVVVV++5oa0eWWNwD/zh3MG/7Hb38SysjIQFpaGq6//nqMGzcOa9asQUNDAxYuXOjroXXasmXLMGPGDAwcOBCnT5/GypUrYTabkZqa6uuhXbT6+nq3TwWFhYXIy8tD7969MWDAAPz+97/HU089haFDh2Lw4MFYvnw5YmNjMWvWLN8NWkFH29e7d2+sWrUKd911F6Kjo1FQUIA//OEPuOqqq5CcnOzDUV9+/DlvAP6XO5g3LrO84evTlLxp7dq1YsCAASIoKEiMGzdOHD582NdD6hLz5s0TMTExIigoSPTr10/MmzdP5Ofn+3pYnfLRRx8JAK0uaWlpQogfT1Fcvny5iIqKEhaLRUydOlUcP37ct4NW0NH2NTY2iltvvVVERkaKwMBAMXDgQLF48WJRVlbm62Fflvw1bwjhf7mDeePyyhuaEEJc2hKJiIiISI1fHsNCRERE/oUFCxERERkeCxYiIiIyPBYsREREZHgsWIiIiMjwWLAQERGR4bFgISIiIsNjwUJERESGx4KFiIiIDI8FCxERERkeCxYiIiIyvP8HXODpCG4iMjAAAAAASUVORK5CYII=",
|
|
"text/plain": [
|
|
"<Figure size 640x480 with 2 Axes>"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"plt.subplot(1, 2, 1)\n",
|
|
"plt.title('permeability')\n",
|
|
"plt.imshow(k_train.squeeze(-1)[0])\n",
|
|
"plt.subplot(1, 2, 2)\n",
|
|
"plt.title('field solution')\n",
|
|
"plt.imshow(u_train.squeeze(-1)[0])\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "89a77ff1",
|
|
"metadata": {},
|
|
"source": [
|
|
"We now create the neural operator class. It is a very simple class, inheriting from `AbstractProblem`."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 69,
|
|
"id": "8b27d283",
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"class NeuralOperatorSolver(AbstractProblem):\n",
|
|
" input_variables = ['u_0']\n",
|
|
" output_variables = ['u']\n",
|
|
" conditions = {'data' : Condition(input_points=LabelTensor(k_train, input_variables), \n",
|
|
" output_points=LabelTensor(u_train, input_variables))}\n",
|
|
"\n",
|
|
"# make problem\n",
|
|
"problem = NeuralOperatorSolver()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "1096cc20",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Solving the problem with a FeedForward Neural Network\n",
|
|
"\n",
|
|
"We will first solve the problem using a Feedforward neural network. We will use the `SupervisedSolver` for solving the problem, since we are training using supervised learning."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 78,
|
|
"id": "e34f18b0",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stderr",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"GPU available: False, used: False\n",
|
|
"TPU available: False, using: 0 TPU cores\n",
|
|
"IPU available: False, using: 0 IPUs\n",
|
|
"HPU available: False, using: 0 HPUs\n"
|
|
]
|
|
},
|
|
{
|
|
"name": "stderr",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"\n",
|
|
" | Name | Type | Params\n",
|
|
"----------------------------------------\n",
|
|
"0 | _loss | MSELoss | 0 \n",
|
|
"1 | _neural_net | Network | 481 \n",
|
|
"----------------------------------------\n",
|
|
"481 Trainable params\n",
|
|
"0 Non-trainable params\n",
|
|
"481 Total params\n",
|
|
"0.002 Total estimated model params size (MB)\n"
|
|
]
|
|
},
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Epoch 99: : 1it [00:00, 15.95it/s, v_num=85, mean_loss=0.105]"
|
|
]
|
|
},
|
|
{
|
|
"name": "stderr",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"`Trainer.fit` stopped: `max_epochs=100` reached.\n"
|
|
]
|
|
},
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Epoch 99: : 1it [00:00, 15.53it/s, v_num=85, mean_loss=0.105]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"# make model\n",
|
|
"model=FeedForward(input_dimensions=1, output_dimensions=1)\n",
|
|
"\n",
|
|
"\n",
|
|
"# make solver\n",
|
|
"solver = SupervisedSolver(problem=problem, model=model)\n",
|
|
"\n",
|
|
"# make the trainer and train\n",
|
|
"trainer = Trainer(solver=solver, max_epochs=100)\n",
|
|
"trainer.train()\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "7b2c35be",
|
|
"metadata": {},
|
|
"source": [
|
|
"The final loss is pretty high... We can calculate the error by importing `LpLoss`."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 79,
|
|
"id": "0e2a6aa4",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Final error training 56.06%\n",
|
|
"Final error testing 55.95%\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"from pina.loss import LpLoss\n",
|
|
"\n",
|
|
"# make the metric\n",
|
|
"metric_err = LpLoss(relative=True)\n",
|
|
"\n",
|
|
"\n",
|
|
"err = float(metric_err(u_train.squeeze(-1), solver.models[0](k_train).squeeze(-1)).mean())*100\n",
|
|
"print(f'Final error training {err:.2f}%')\n",
|
|
"\n",
|
|
"err = float(metric_err(u_test.squeeze(-1), solver.models[0](k_test).squeeze(-1)).mean())*100\n",
|
|
"print(f'Final error testing {err:.2f}%')"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "6b5e5aa6",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Solving the problem with a Fuorier Neural Operator (FNO)\n",
|
|
"\n",
|
|
"We will now move to solve the problem using a FNO. Since we are learning operator this approach is better suited, as we shall see."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 70,
|
|
"id": "9af523a5",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stderr",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"GPU available: False, used: False\n",
|
|
"TPU available: False, using: 0 TPU cores\n",
|
|
"IPU available: False, using: 0 IPUs\n",
|
|
"HPU available: False, using: 0 HPUs\n",
|
|
"\n",
|
|
" | Name | Type | Params\n",
|
|
"----------------------------------------\n",
|
|
"0 | _loss | MSELoss | 0 \n",
|
|
"1 | _neural_net | Network | 591 K \n",
|
|
"----------------------------------------\n",
|
|
"591 K Trainable params\n",
|
|
"0 Non-trainable params\n",
|
|
"591 K Total params\n",
|
|
"2.364 Total estimated model params size (MB)\n"
|
|
]
|
|
},
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Epoch 19: : 1it [00:02, 2.65s/it, v_num=84, mean_loss=0.0294]"
|
|
]
|
|
},
|
|
{
|
|
"name": "stderr",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"`Trainer.fit` stopped: `max_epochs=20` reached.\n"
|
|
]
|
|
},
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Epoch 19: : 1it [00:02, 2.67s/it, v_num=84, mean_loss=0.0294]\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"# make model\n",
|
|
"lifting_net = torch.nn.Linear(1, 24)\n",
|
|
"projecting_net = torch.nn.Linear(24, 1)\n",
|
|
"model = FNO(lifting_net=lifting_net,\n",
|
|
" projecting_net=projecting_net,\n",
|
|
" n_modes=16,\n",
|
|
" dimensions=2,\n",
|
|
" inner_size=24,\n",
|
|
" padding=11)\n",
|
|
"\n",
|
|
"\n",
|
|
"# make solver\n",
|
|
"solver = SupervisedSolver(problem=problem, model=model)\n",
|
|
"\n",
|
|
"# make the trainer and train\n",
|
|
"trainer = Trainer(solver=solver, max_epochs=20)\n",
|
|
"trainer.train()\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "84964cb9",
|
|
"metadata": {},
|
|
"source": [
|
|
"We can clearly see that with 1/3 of the total epochs the loss is lower. Let's see in testing.. Notice that the number of parameters is way higher than a `FeedForward` network. We suggest to use GPU or TPU for a speed up in training."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 77,
|
|
"id": "58e2db89",
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"Final error training 26.05%\n",
|
|
"Final error testing 25.58%\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"err = float(metric_err(u_train.squeeze(-1), solver.models[0](k_train).squeeze(-1)).mean())*100\n",
|
|
"print(f'Final error training {err:.2f}%')\n",
|
|
"\n",
|
|
"err = float(metric_err(u_test.squeeze(-1), solver.models[0](k_test).squeeze(-1)).mean())*100\n",
|
|
"print(f'Final error testing {err:.2f}%')"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"id": "26e3a6e4",
|
|
"metadata": {},
|
|
"source": [
|
|
"As we can see the loss is way lower!"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## What's next?\n",
|
|
"\n",
|
|
"We have made a very simple example on how to use the `FNO` for learning neural operator. Currently in **PINA** we implement 1D/2D/3D cases. We suggest to extend the tutorial using more complex problems and train for longer, to see the full potential of neural operators."
|
|
]
|
|
}
|
|
],
|
|
"metadata": {
|
|
"interpreter": {
|
|
"hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49"
|
|
},
|
|
"kernelspec": {
|
|
"display_name": "Python 3.9.0 64-bit",
|
|
"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.16"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 5
|
|
}
|