{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e80567a6",
   "metadata": {},
   "source": [
    "# Tutorial: Two dimensional Darcy flow using the Fourier Neural Operator\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/tutorial5/tutorial.ipynb)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8762bbe5",
   "metadata": {},
   "source": [
    "In this tutorial we are going to solve the Darcy flow problem in two dimensions, 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 operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5f2744dc",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-09-19T13:35:28.837348Z",
     "start_time": "2024-09-19T13:35:27.611334Z"
    }
   },
   "outputs": [],
   "source": [
    "## routine needed to run the notebook on Google Colab\n",
    "try:\n",
    "  import google.colab\n",
    "  IN_COLAB = True\n",
    "except:\n",
    "  IN_COLAB = False\n",
    "if IN_COLAB:\n",
    "  !pip install \"pina-mathlab\"\n",
    "  !pip install scipy\n",
    "  # get the data\n",
    "  !wget https://github.com/mathLab/PINA/raw/refs/heads/master/tutorials/tutorial5/Data_Darcy.mat\n",
    "\n",
    "import torch\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# !pip install scipy  # install scipy\n",
    "from scipy import io\n",
    "from pina.model import FNO, FeedForward  # let's import some models\n",
    "from pina import Condition, LabelTensor\n",
    "from pina.solver import SupervisedSolver\n",
    "from pina.trainer import Trainer\n",
    "from pina.problem import AbstractProblem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cf5b181",
   "metadata": {},
   "source": [
    "## Data Generation\n",
    "\n",
    "We will focus on solving a specific 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": 2,
   "id": "2ffb8a4c",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-09-19T13:35:28.989631Z",
     "start_time": "2024-09-19T13:35:28.952744Z"
    }
   },
   "outputs": [],
   "source": [
    "# download the dataset\n",
    "data = io.loadmat(\"Data_Darcy.mat\")\n",
    "\n",
    "# extract data (we use only 100 data for train)\n",
    "k_train = LabelTensor(torch.tensor(data['k_train'], dtype=torch.float).unsqueeze(-1), \n",
    "                      labels={3:{'dof': ['u0'], 'name': 'k_train'}})\n",
    "u_train = LabelTensor(torch.tensor(data['u_train'], dtype=torch.float).unsqueeze(-1),\n",
    "                     labels={3:{'dof': ['u'], 'name': 'u_train'}})\n",
    "k_test =  LabelTensor(torch.tensor(data['k_test'], dtype=torch.float).unsqueeze(-1),\n",
    "                      labels={3:{'dof': ['u0'], 'name': 'k_test'}})\n",
    "u_test= LabelTensor(torch.tensor(data['u_test'], dtype=torch.float).unsqueeze(-1),\n",
    "                   labels={3:{'dof': ['u'], 'name': 'u_test'}})\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": 3,
   "id": "c8501b6f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-09-19T13:35:29.108381Z",
     "start_time": "2024-09-19T13:35:29.031076Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAEjCAYAAAARyVqhAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAO95JREFUeJzt3Xt0U2W6P/DvTtqmhbYpSOkFuRW5CNiiKJ0qV0FL9SCgw6Wj0jKCMx4YndWFo3UJFG894lEZgQPqDBQHVHQG4cwRGbEKHAdQAet15NDa0iJtaQu90+ay398f/hoNveR9aUN2w/ezVhZk58nOu5Psp0+SvZ9XE0IIEBERERmYydcDICIiIvKEBQsREREZHgsWIiIiMjwWLERERGR4LFiIiIjI8FiwEBERkeGxYCEiIiLDY8FCREREhseChYiIiAyPBQt1W1lZWdA0DZWVlR5jBw0ahPT0dNf1ffv2QdM07Nu3z7UsPT0dgwYN6vqBEhnMZ599hhtvvBE9e/aEpmnIy8tz7U8XQ3bfKSoqgqZpyMnJuajHuRg5OTnQNA1FRUVdut7Jkydj8uTJXbpO6hgLFqJ2NDY2Iisry62oIeru7HY75syZg7Nnz+LFF1/EX/7yFwwcONDXwzKkb7/9FllZWV1e7NDFCfD1AIguhePHj8Nk6rg+f/XVV6Hruut6Y2MjVq1aBQD8JEV+o6CgACdPnsSrr76KRYsWuZY//vjjePTRR304MuP59ttvsWrVKkyePLnVN0jvv/++bwZ1GWPBQh41NjaiR48evh5Gp1gsFo8xgYGBl2AkRL515swZAEBERITb8oCAAAQE8E+CrKCgIF8P4bLDn4S6uZbfnb/77jvMnTsX4eHhuOKKK/DQQw+hqanJLXbr1q0YO3YsQkJC0Lt3b8yfPx8lJSVuMZMnT8bo0aNx9OhRTJw4ET169MBjjz3m+u35P//zP7F+/XrExcWhR48euPXWW1FSUgIhBJ588klceeWVCAkJwcyZM3H27NlW433vvfcwYcIE9OzZE2FhYbj99tvxzTffuMV8+eWXSE9PR1xcHIKDgxEdHY1f//rXqKqqavM5qKys9LjtFx7D0paf/w5fVFSEyMhIAMCqVaugaRo0TUNWVhY2b94MTdPw+eeft1rHM888A7PZjB9++KHDxyLyhfT0dEyaNAkAMGfOHGia5vr2sL1jWGTyRluqq6uRnp4Oq9WKiIgIpKWlobq6Wmqcdrsdq1atwtChQxEcHIwrrrgC48ePx969e93iPvzwQ1c+iYiIwMyZM/Gvf/3L4/pb9uUL/TxP5OTkYM6cOQCAKVOmuHJAy0/EbR3DcubMGdx3332IiopCcHAwEhISsGXLFreYn+fSV155BUOGDIHFYsENN9yAzz77TOr5uVyxnPYTc+fOxaBBg5CdnY3Dhw/jpZdewrlz5/Daa68BAJ5++mksX74cc+fOxaJFi1BRUYG1a9di4sSJ+Pzzz90+bVVVVSElJQXz58/HPffcg6ioKNdt27Ztg81mw+9+9zucPXsWq1evxty5c3HzzTdj3759eOSRR5Cfn4+1a9di2bJl2LRpk+u+f/nLX5CWlobk5GQ8++yzaGxsxIYNGzB+/Hh8/vnnrmJh7969+P7777Fw4UJER0fjm2++wSuvvIJvvvkGhw8fbpVUPW37xYiMjMSGDRvwwAMPYPbs2bjzzjsBAPHx8Rg8eDCWLFmCbdu24dprr3W737Zt2zB58mT069fvoh+byFt+85vfoF+/fnjmmWfw4IMP4oYbbnDbvy+kkjd+TgiBmTNn4uOPP8Zvf/tbXH311XjnnXeQlpYmNc6srCxkZ2dj0aJFGDduHGpra3HkyBEcO3YMt9xyCwDggw8+QEpKCuLi4pCVlYXz589j7dq1uOmmm3Ds2LFOH0A/ceJEPPjgg3jppZfw2GOP4eqrrwYA178XOn/+PCZPnoz8/HwsXboUgwcPxttvv4309HRUV1fjoYcecot//fXXUVdXh9/85jfQNA2rV6/GnXfeie+//57f9rZHULe2cuVKAUDccccdbsv//d//XQAQX3zxhSgqKhJms1k8/fTTbjFfffWVCAgIcFs+adIkAUBs3LjRLbawsFAAEJGRkaK6utq1PDMzUwAQCQkJwm63u5anpqaKoKAg0dTUJIQQoq6uTkRERIjFixe7rbesrExYrVa35Y2Nja2284033hAAxIEDB5S2vcXAgQNFWlqa6/pHH30kAIiPPvrItSwtLU0MHDjQdb2iokIAECtXrmw1ntTUVBEbGyucTqdr2bFjxwQAsXnz5lbxREbR8t5/++233Za37E8tVPLGhfvOzp07BQCxevVq1zKHwyEmTJggtY8kJCSI22+/vcOYMWPGiL59+4qqqirXsi+++EKYTCaxYMEC17LNmzcLAKKwsNC1rL39+sI88fbbb7fKEy0mTZokJk2a5Lq+Zs0aAUBs3brVtcxms4mkpCQRGhoqamtrhRA/5dIrrrhCnD171hW7a9cuAUD8/e9/73C7L2f8SchPLFmyxO367373OwDA7t27sWPHDui6jrlz56KystJ1iY6OxtChQ/HRRx+53ddisWDhwoVtPs6cOXNgtVpd1xMTEwEA99xzj9vv34mJibDZbK6fRvbu3Yvq6mqkpqa6jcFsNiMxMdFtDCEhIa7/NzU1obKyEr/4xS8AAMeOHVPadm9ZsGABTp8+7Tbubdu2ISQkBHfddZfXHpfoUlHNGz+3e/duBAQE4IEHHnAtM5vNrn3Tk4iICHzzzTc4ceJEm7eXlpYiLy8P6enp6N27t2t5fHw8brnlFq/u++3ZvXs3oqOjkZqa6loWGBiIBx98EPX19di/f79b/Lx589CrVy/X9QkTJgAAvv/++0sz4G6IPwn5iaFDh7pdHzJkCEwmE4qKimAymSCEaBXT4sKvH/v169fuAWUDBgxwu95SvPTv37/N5efOnQMAV+K5+eab21xveHi46/9nz57FqlWr8Oabb7oOEGxRU1PT6r4dbbu33HLLLYiJicG2bdswdepU6LqON954AzNnzkRYWJjXHpfoUjlx4oRS3vi5kydPIiYmBqGhoW7Lhw8fLvXYTzzxBGbOnIlhw4Zh9OjRmD59Ou69917Ex8e71t/e+q6++mr84x//QENDA3r27Cn1eF3h5MmTGDp0aKuzEVt+QmoZc4sLc2lL8dKSM6k1Fix+6ufHeei6Dk3T8N5778FsNreKvTCp/Pwbjgu1df+OlgshXGMAfjyOJTo6ulXcz7+dmTt3Lg4ePIiHH34YY8aMQWhoKHRdx/Tp091OO27PxTa/UmE2m/GrX/0Kr776Kv7rv/4L//znP3H69Gncc889Xn9soktBNW90pYkTJ6KgoAC7du3C+++/jz/96U948cUXsXHjRrdTsbua0+n02rov5ClnUmssWPzEiRMnMHjwYNf1/Px86LqOQYMGwWw2QwiBwYMHY9iwYT4Z35AhQwAAffv2xbRp09qNO3fuHHJzc7Fq1SqsWLHCtby9r4Zbbmtv2zvDU+GzYMECPP/88/j73/+O9957D5GRkUhOTu7UYxIZxZAhQy46bwwcOBC5ubmor693K2yOHz8uvY7evXtj4cKFWLhwIerr6zFx4kRkZWVh0aJFrkZ3ba3vu+++Q58+fTr8dqVXr16tzliy2WwoLS11W6by4WfgwIH48ssvoeu627cs3333net26hwew+In1q9f73Z97dq1AICUlBTceeedMJvNWLVqVavqXQjR7unCXSk5ORnh4eF45plnYLfbW91eUVEB4KdPHReOc82aNe2uu6Nt74yW3jPtnYoZHx+P+Ph4/OlPf8Lf/vY3zJ8/n30syG90Jm/cdtttcDgc2LBhg2uZ0+l07ZueXLju0NBQXHXVVWhubgYAxMTEYMyYMdiyZYvb/vn111/j/fffx2233dbh+ocMGYIDBw64LXvllVdafcPSUvTInI592223oaysDNu3b3ctczgcWLt2LUJDQ12nk9PFY3b1E4WFhbjjjjswffp0HDp0CFu3bsWvfvUrJCQkAACeeuopZGZmoqioCLNmzUJYWBgKCwvxzjvv4P7778eyZcu8Or7w8HBs2LAB9957L6677jrMnz8fkZGRKC4uxrvvvoubbroJ69atQ3h4OCZOnIjVq1fDbrejX79+eP/991FYWHjR236xQkJCMHLkSGzfvh3Dhg1D7969MXr0aIwePdoVs2DBAtdzx5+DyJ8MGTLkovPGjBkzcNNNN+HRRx9FUVERRo4ciR07drR5DFpbRo4cicmTJ2Ps2LHo3bs3jhw5gr/+9a9YunSpK+a5555DSkoKkpKScN9997lOa7ZarW32WPm5RYsW4be//S3uuusu3HLLLfjiiy/wj3/8A3369HGLGzNmDMxmM5599lnU1NTAYrHg5ptvRt++fVut8/7778fLL7+M9PR0HD16FIMGDcJf//pX/POf/8SaNWt4bFtX8NHZSdRFWk5F/Pbbb8Uvf/lLERYWJnr16iWWLl0qzp8/7xb7t7/9TYwfP1707NlT9OzZU4wYMUIsWbJEHD9+3BUzadIkMWrUqFaP03Iq3nPPPee2vL1TJFtOJfzss89axScnJwur1SqCg4PFkCFDRHp6ujhy5Igr5tSpU2L27NkiIiJCWK1WMWfOHHH69OlWpyKqbPvFnNYshBAHDx4UY8eOFUFBQW2eCllaWirMZrMYNmxYq+eMyIhkT2tuIZM32tp3qqqqxL333ivCw8OF1WoV9957r/j888+lTmt+6qmnxLhx40RERIQICQkRI0aMEE8//bSw2WxucR988IG46aabREhIiAgPDxczZswQ3377rVtMW6c1O51O8cgjj4g+ffqIHj16iOTkZJGfn98qTwghxKuvviri4uKE2Wx2yxkXntYshBDl5eVi4cKFok+fPiIoKEhcc801rba1vVwqRPunW9OPNCF4hE93lpWVhVWrVqGioqLVpwPyvsrKSsTExGDFihVYvny5r4dDROS3eAwLUSfk5OTA6XTi3nvv9fVQiIj8Go9hIboIH374Ib799ls8/fTTmDVrVqfPSCIioo6xYCG6CE888QQOHjyIm266SfrMByIiung8hoWIiIgMj8ewEBERkeGxYCEiIiLD84tjWHRdx+nTpxEWFnZJ5pEhotaEEKirq0NsbGyrCeCMirmDyLdU8oZfFCynT59uNVswEflGSUkJrrzySl8PQwpzB5ExyOQNrxUs69evx3PPPYeysjIkJCRg7dq1GDduXLvxb7/9NpYvX46ioiIMHToUzz77rMf5IFq0tDy+cuXjMAUHe4zXPE/466ZnscKnRYUPabVXt55Tpz2W0vancu+sHmfkj7sOrFc7RjukQn4bYZJ/8nSz2qfhhmj5t3rwOfk3SMB5+VhN8fj2oOpm6Vg9sO2ZX9sSkF/qOehnnJWVUnEO2PExdneqBfmlzBvAT7nj5LFBCA/1vJ9XOBuk1w0AlU753HFWb3+W9Audc7Y/sd+FahTWCwDnHPKzMFfZ5cdRaVOb3bmyWT6+8nwP6dizdfKxAGCr8fw3pUVAtXyesVTJ57DgKsW8Wyk/63RIeaN0rLn0rNI4HGXlnmMU8oZXCpbt27cjIyMDGzduRGJiItasWYPk5GQcP368zTkYDh48iNTUVGRnZ+Pf/u3f8Prrr2PWrFk4duyY27wt7Wn5KtcUHOyVgsVs8U7BYgqR/yNjtnivYDEHye8MKrEAEBAgv41KBUuAWsFiDpJ/qwcEKhQsDoWCRVd87hSKMl3heQ4wBSmNQ9Mk33uiJf7iflq51Hnj52MNDzUhPMzzft6kUICoxjfrCrFO+dfb5lRL8+cdCn907fJ5Kcim9r4LDJCPDzBZpGPNTvkCBABMNvl4U5P8c2e2yO8nynk3UL5gCTDLx5oVcwdkcodC3vDKD80vvPACFi9ejIULF2LkyJHYuHEjevTogU2bNrUZ/8c//hHTp0/Hww8/jKuvvhpPPvkkrrvuOqxbt84bwyMiA2LeIKKOdHnBYrPZcPToUUybNu2nBzGZMG3aNBw6dKjN+xw6dMgtHgCSk5PbjW9ubkZtba3bhYi6r0uRNwDmDqLurMsLlsrKSjidTkRFRbktj4qKQllZWZv3KSsrU4rPzs6G1Wp1XXjQHFH3dinyBsDcQdSddY9zDy+QmZmJmpoa16WkpMTXQyKiboC5g6j76vKDbvv06QOz2Yzycvejg8vLyxEdHd3mfaKjo5XiLRYLLBb5g6yIyNguRd4AmDuIurMu/4YlKCgIY8eORW5urmuZruvIzc1FUlJSm/dJSkpyiweAvXv3thtPRP6FeYOIPPHKac0ZGRlIS0vD9ddfj3HjxmHNmjVoaGjAwoULAQALFixAv379kJ2dDQB46KGHMGnSJDz//PO4/fbb8eabb+LIkSN45ZVXvDE8IjIg5g0i6ohXCpZ58+ahoqICK1asQFlZGcaMGYM9e/a4DpArLi52a8F744034vXXX8fjjz+Oxx57DEOHDsXOnTuleym0CDivwSQ8n8sdUK/WJ6JuiHyvjeD+dfIrrpJvYNQU65BfL4AexQp9R+T7BkEPVHzuBsifty/x0v0Uq9DeBVDrvVMzWH7lmsJAQk/J9zsAgOYI+aZcYQXy7zvRL1JpHJrkmTSaMAHyve5a8VXeIOqQt2ZsMMpMECp9k3w85YYmhGL7TQOqra2F1WpF3IqnpRrHqRYsTVHeKVgaFQoW6GpjVilYQkvk3wKqhYJQeH8bpWCxhcsPRFPYe1QLFpVEolKwaE61XV78q0AqziHs+Kj5LdTU1CA8PFzpMXylJXec+784qcZxZxQ73VYoNI6rUuhIe9Yp3wW22qnW2fWsQqfbSrt8bIVNrQNyRZPCus/LF/dVtfKxANCs0un2rEJHbZVOtxVq+2yPCoVOt2Ve7HT7w2nPMcKOfdgllTe65VlCREREdHlhwUJERESGx4KFiIiIDI8FCxERERkeCxYiIiIyPBYsREREZHgsWIiIiMjwWLAQERGR4bFgISIiIsNjwUJERESG55W5hHzFFumAKcTznDu2KLU2x1cP8dxeuMWpGqt0bHjfeunY5i8jpGMBtZb49f3lW0Sb1KY0UqIrtNvX5acoAgA4gxWmH1Dot68yRYAeoDafQI9y+XE0Rcm3X+9RcE5pHFqE3HvapNuAM0qrpkvAqZIMADgVJrlxKnzm1VXm3vAilalzfryDwt8LlelFVKYtUf1qQWUblZ8Q3+E3LERERGR4LFiIiIjI8FiwEBERkeGxYCEiIiLDY8FCREREhseChYiIiAyPBQsREREZXpcXLNnZ2bjhhhsQFhaGvn37YtasWTh+/HiH98nJyYGmaW6X4ODgrh4aERkU8wYRedLlBcv+/fuxZMkSHD58GHv37oXdbsett96KhoaGDu8XHh6O0tJS1+XkyZNdPTQiMijmDSLypMs73e7Zs8ftek5ODvr27YujR49i4sSJ7d5P0zRER0d39XCIqBtg3iAiT7zemr+mpgYA0Lt37w7j6uvrMXDgQOi6juuuuw7PPPMMRo0a1WZsc3MzmpubXddra2sBAFqQE1qQ0+OYovrWyA7/x8dzyj9Ng3rJtz3/Om+QdKwWqjadACAfb7LJt2ZWbeKs1G4/RJeOFSGeX+efM0tM2dAiIEB+3bZG+TkCGoLVdjeTTf7JCz6n0E7dGqI0DnOV5HtaV3tNOuKNvAG0nzv8mUr7fADQFfrAe7Pdvq6cbeRoKq32AcV2+wrTepgUVqz4W4jKlCFKL6GP2/h79aBbXdfx+9//HjfddBNGjx7dbtzw4cOxadMm7Nq1C1u3boWu67jxxhtx6tSpNuOzs7NhtVpdl/79+3trE4joEvNW3gCYO4i6M68WLEuWLMHXX3+NN998s8O4pKQkLFiwAGPGjMGkSZOwY8cOREZG4uWXX24zPjMzEzU1Na5LSUmJN4ZPRD7grbwBMHcQdWde+0lo6dKl+J//+R8cOHAAV155pdJ9AwMDce211yI/P7/N2y0WCywWS1cMk4gMxJt5A2DuIOrOuvwbFiEEli5dinfeeQcffvghBg8erLwOp9OJr776CjExMV09PCIyIOYNIvKky79hWbJkCV5//XXs2rULYWFhKCsrAwBYrVaEhPx4sN+CBQvQr18/ZGdnAwCeeOIJ/OIXv8BVV12F6upqPPfcczh58iQWLVrU1cMjIgNi3iAiT7q8YNmwYQMAYPLkyW7LN2/ejPT0dABAcXExTKafvtw5d+4cFi9ejLKyMvTq1Qtjx47FwYMHMXLkyK4eHhEZEPMGEXnS5QWLEJ5P69q3b5/b9RdffBEvvvhiVw+FiLoJ5g0i8oRzCREREZHhsWAhIiIiw2PBQkRERIbn9db8l1JQiB3mHp57EvcPq1Zar0mhlfMnX10lv2KLQit6s2I7aYd8C2VdpVN7gNo4TAot8UNDm6RjI0M7nhSvVXxIvVK8rLKGcOnYkyV9lNZti5Dvr60HKkyvYJd/39Glo9IS36bQe121fb7TS635vdVqH1CbMkS1Nb+mkHuVWuKrTFui+JdaN8s/I8Ks8L2FP7fmJyIiIuoKLFiIiIjI8FiwEBERkeGxYCEiIiLDY8FCREREhseChYiIiAyPBQsREREZHgsWIiIiMjwWLERERGR4LFiIiIjI8PyqNb/QNei659bBxyv7Kq237pR8+3Wl1vWBCi3SJbbLTZD8uk1BTunY4B42pWFEhsm30B8SXikdOyK0VGkcsYHV0rENukU69hNznHRsWViYdCwAaHqQdKxToTW/qUF+CgQA0ALk0oSms+X/hZwKTeNVYlXa+NuFWpq3K/SMdyp85nUqThHgLcrd5U0KrfkV8r9Ku309QG3QQqk1v8K6TWzNT0RERNQhFixERERkeF1esGRlZUHTNLfLiBEjOrzP22+/jREjRiA4OBjXXHMNdu/e3dXDIiIDY94gIk+88g3LqFGjUFpa6rp8/PHH7cYePHgQqampuO+++/D5559j1qxZmDVrFr7++mtvDI2IDIp5g4g64pWCJSAgANHR0a5Lnz592o394x//iOnTp+Phhx/G1VdfjSeffBLXXXcd1q1b542hEZFBMW8QUUe8UrCcOHECsbGxiIuLw913343i4uJ2Yw8dOoRp06a5LUtOTsahQ4favU9zczNqa2vdLkTUvXk7bwDMHUTdWZcXLImJicjJycGePXuwYcMGFBYWYsKECairq2szvqysDFFRUW7LoqKiUFZW1u5jZGdnw2q1ui79+/fv0m0gokvrUuQNgLmDqDvr8oIlJSUFc+bMQXx8PJKTk7F7925UV1fjrbfe6rLHyMzMRE1NjetSUlLSZesmokvvUuQNgLmDqDvzeuO4iIgIDBs2DPn5+W3eHh0djfLycrdl5eXliI6ObnedFosFFot8cy8i6l68kTcA5g6i7szrfVjq6+tRUFCAmJiYNm9PSkpCbm6u27K9e/ciKSnJ20MjIoNi3iCiC3V5wbJs2TLs378fRUVFOHjwIGbPng2z2YzU1FQAwIIFC5CZmemKf+ihh7Bnzx48//zz+O6775CVlYUjR45g6dKlXT00IjIo5g0i8qTLfxI6deoUUlNTUVVVhcjISIwfPx6HDx9GZGQkAKC4uBgm00910o033ojXX38djz/+OB577DEMHToUO3fuxOjRo5UfW4gfL57UnbQqrddkl4/VgxVW7JCfs0NljgoAgFk+PkBhLqFePc8rDWN4RLnnoP8vKbxAOva64PbPIGlLpMkhHVvilP/JoLj5CulYk8KcJACgMvWK2ab4/lCgV9fIxQm1eaZ+zpd5w5t0hRdRZc4fm9J8P2rzv+gK8Q5dfhwq8x/9GO+deWtMJrU5rzSF/VZXyLsqOV0oziWkMk+RCFB4XUy+bY7f5QXLm2++2eHt+/bta7Vszpw5mDNnTlcPhYi6CeYNIvKEcwkRERGR4bFgISIiIsNjwUJERESGx4KFiIiIDI8FCxERERkeCxYiIiIyPBYsREREZHgsWIiIiMjwWLAQERGR4Xl9tuZLKejzUJgtnnvjBwaqrVdXmNxVNCq0cVYYhx6k1nrdqXmnVXvPQLX269GWWunYuKAz8rGK79xQU6h0bKOol441a/Jtvp1Otc8HAQpTQmi6QpvvYLUdQLfJDUQX8tMfdFdOmbk/fsam8JlQpd2+Sht/u8J6AcCu0m5foY2/aqt9b7XmV6XSml/lKwCFp1mp1f6P8fLPnTArPM+ab18TfsNCREREhseChYiIiAyPBQsREREZHgsWIiIiMjwWLERERGR4LFiIiIjI8FiwEBERkeF1ecEyaNAgaJrW6rJkyZI243NyclrFBgd77qVCRP6FuYOIOtLljeM+++wzOJ1O1/Wvv/4at9xyC+bMmdPufcLDw3H8+HHXdc3HzWmI6NJj7iCijnR5wRIZGel2/T/+4z8wZMgQTJo0qd37aJqG6Ojorh4KEXUjzB1E1BGvtua32WzYunUrMjIyOvzkU19fj4EDB0LXdVx33XV45plnMGrUqHbjm5ub0dzc7LpeW/tj+3dh/vHiiaVarb22MCm0n1ZoiazS8t8RovbJUaWrtT1Y/m1Qbw9SGkejUz6+UeEJaRQNSuOA3iQdWuEMkY6tssm3/HfY1FqkByl0ujfZ5d/T2nm16RVMQXKt/E1CAPJPc4cude6Q5fQc4kYX8r+66wq/0DsVYu2Kfd1VWvk7FRKNSht/b1IdhVK8Qht/lRkT1Fvzq8SqtPH37WGvXn30nTt3orq6Gunp6e3GDB8+HJs2bcKuXbuwdetW6LqOG2+8EadOnWr3PtnZ2bBara5L//79vTB6IvIV5g4iupBXC5Y///nPSElJQWxsbLsxSUlJWLBgAcaMGYNJkyZhx44diIyMxMsvv9zufTIzM1FTU+O6lJSUeGP4ROQjzB1EdCGv/SR08uRJfPDBB9ixY4fS/QIDA3HttdciPz+/3RiLxQKLReH3FCLqNpg7iKgtXvuGZfPmzejbty9uv/12pfs5nU589dVXiImJ8dLIiMjImDuIqC1eKVh0XcfmzZuRlpaGgAD3L3EWLFiAzMxM1/UnnngC77//Pr7//nscO3YM99xzD06ePIlFixZ5Y2hEZGDMHUTUHq/8JPTBBx+guLgYv/71r1vdVlxcDJPppzrp3LlzWLx4McrKytCrVy+MHTsWBw8exMiRI70xNCIyMOYOImqPVwqWW2+9FUK0fXrXvn373K6/+OKLePHFF70xDCLqZpg7iKg9nEuIiIiIDI8FCxERERkeCxYiIiIyPBYsREREZHhenUvoUgsvciIgUHW2D8+U5hKSm3YFAOCwyK/XFKo2A4amy9eiTWb5+X4qLPJz5wDAv4LlJ6azBpyXjrWpTMQBIFizS8d+29xPOja/ro90rN6g8OYAYFKZ8kdleiyT2ucULViu0ZomtC6bS8ioFKZsAgA4FWaiUXlPq8z3ozqHj64wP5BDVxiHygRnivGKL4sSTVNYu8omKq1XcS45hV1c5e+b6ji6Gr9hISIiIsNjwUJERESGx4KFiIiIDI8FCxERERkeCxYiIiIyPBYsREREZHgsWIiIiMjwWLAQERGR4bFgISIiIsNjwUJERESG51et+YPPORAQ4PAYpwcotrlX6KCssm5nsHy9GNCsNmaTXaEWVWjj3yx6KI3jOz1KOra2OVg69tvQGKVxhJjlW/OXN4VJxxZV9paODahVm07AbJN/46m8R6HSihsAAiWnbtDVVns5cCp8JtQV+qmrxDpV+rRDccwKvehVW/MLxXjp9arGe2kcKm38VYeg9JIrTSfA1vxEREREHVIuWA4cOIAZM2YgNjYWmqZh586dbrcLIbBixQrExMQgJCQE06ZNw4kTJzyud/369Rg0aBCCg4ORmJiITz/9VHVoRGRQzBtE1FnKBUtDQwMSEhKwfv36Nm9fvXo1XnrpJWzcuBGffPIJevbsieTkZDQ1tT+N6/bt25GRkYGVK1fi2LFjSEhIQHJyMs6cOaM6PCIyIOYNIuos5YIlJSUFTz31FGbPnt3qNiEE1qxZg8cffxwzZ85EfHw8XnvtNZw+fbrVJ6qfe+GFF7B48WIsXLgQI0eOxMaNG9GjRw9s2rRJdXhEZEDMG0TUWV16DEthYSHKysowbdo01zKr1YrExEQcOnSozfvYbDYcPXrU7T4mkwnTpk1r9z7Nzc2ora11uxBR93Sp8gbA3EHUnXVpwVJWVgYAiIpyPzMkKirKdduFKisr4XQ6le6TnZ0Nq9XquvTv378LRk9EvnCp8gbA3EHUnXXLs4QyMzNRU1PjupSUlPh6SETUDTB3EHVfXVqwREdHAwDKy8vdlpeXl7tuu1CfPn1gNpuV7mOxWBAeHu52IaLu6VLlDYC5g6g769KCZfDgwYiOjkZubq5rWW1tLT755BMkJSW1eZ+goCCMHTvW7T66riM3N7fd+xCR/2DeICIZyp1u6+vrkZ+f77peWFiIvLw89O7dGwMGDMDvf/97PPXUUxg6dCgGDx6M5cuXIzY2FrNmzXLdZ+rUqZg9ezaWLl0KAMjIyEBaWhquv/56jBs3DmvWrEFDQwMWLlzY+S0kIp9j3iCizlIuWI4cOYIpU6a4rmdkZAAA0tLSkJOTgz/84Q9oaGjA/fffj+rqaowfPx579uxBcPBPbdcLCgpQWVnpuj5v3jxUVFRgxYoVKCsrw5gxY7Bnz55WB9R5EljTjACJ7ud6cKDSejWnfM9xYVZoa31ePtbcpNrWXSVWvt2y2ab2pVxzY4h0bEm1ZAt4AKd69lIahylQ/jXUbQrPda38LhRcqzi9gkJrfpXW3SJQ7b1kCpWbjkHTzUBl27cZOW+oUJ19QKmFvsIX3k6VlvhKvdfVW+h7a73eavuv2mpfqZW/SrDqHAEKlKbq6EY0IUS337Ta2lpYrVZMufZRBJgtHuMNU7AEycc6QtT+yNhD5eObw+V34OZeajt7cy/5t5c9wikdq/X0PGfUzxmiYDmj9hr2KJV/7kJPyz8fIT/UKY3DVHdeKs6hN+ODonWoqanpNseGtOSOc/8Xh/Awz/tjgb1eaf0lDvnn4QeHfBFe4ZCf7+qco6d0LACcs8vPF1bnkJ//q8YmHwsAdXb5+Hqb/Ied+ibPfyN+7vx5+XU76+X/tpgV5hYLOqf2QbHHGYXc8YN87uhRcE5pHM7j+R5jHMKOfdgllTe65VlCREREdHlhwUJERESGx4KFiIiIDI8FCxERERkeCxYiIiIyPBYsREREZHgsWIiIiMjwWLAQERGR4bFgISIiIsNjwUJERESGpzyXkD8w1zQp3sE7c2uYFOZ0MTUrziXULP/SBpyXX3dgg1qNa6mWf+7sofJjdoSovXV1+e7aSnN8mBTmbApU6+quNB+UytwhzlC11uQwyb3muvOyTCeG51SYzwhQnP9IYV4elfX+GO+d+YFU5xJSmqhLl4/VlGLlhwDAMHMadTV+w0JERESGx4KFiIiIDI8FCxERERkeCxYiIiIyPBYsREREZHgsWIiIiMjwlAuWAwcOYMaMGYiNjYWmadi5c6frNrvdjkceeQTXXHMNevbsidjYWCxYsACnT5/ucJ1ZWVnQNM3tMmLECOWNISJjYt4gos5SLlgaGhqQkJCA9evXt7qtsbERx44dw/Lly3Hs2DHs2LEDx48fxx133OFxvaNGjUJpaanr8vHHH6sOjYgMinmDiDpLudNTSkoKUlJS2rzNarVi7969bsvWrVuHcePGobi4GAMGDGh/IAEBiI6OVh0OEXUDzBtE1FleP4alpqYGmqYhIiKiw7gTJ04gNjYWcXFxuPvuu1FcXNxubHNzM2pra90uROQ/vJE3AOYOou7Mq720m5qa8MgjjyA1NRXh4eHtxiUmJiInJwfDhw9HaWkpVq1ahQkTJuDrr79GWFhYq/js7GysWrWq1XJTow0mmTb6mmJrZtW2yJI0m0M+tkmtNb/pvPy6AxoUWvPXqb1lnMHy63aEKLQEt6i9hs4ghRbiak+1NJP8SwIACGyUf+PZe8o/d3qAWmv+IJPcc+dwdM2O4q28AbSfO+jiqLTb16G2z6q05nfqCuNQaIkPALpTpTW/fKhKu32TUz4WADSFeJVpPSCZC7zFa9+w2O12zJ07F0IIbNiwocPYlJQUzJkzB/Hx8UhOTsbu3btRXV2Nt956q834zMxM1NTUuC4lJSXe2AQiusS8mTcA5g6i7swr37C0JJ2TJ0/iww8/7PBTUlsiIiIwbNgw5Ofnt3m7xWKBxaI4gRsRGZq38wbA3EHUnXX5NywtSefEiRP44IMPcMUVVyivo76+HgUFBYiJienq4RGRATFvEJEnygVLfX098vLykJeXBwAoLCxEXl4eiouLYbfb8ctf/hJHjhzBtm3b4HQ6UVZWhrKyMthsNtc6pk6dinXr1rmuL1u2DPv370dRUREOHjyI2bNnw2w2IzU1tfNbSEQ+x7xBRJ2l/JPQkSNHMGXKFNf1jIwMAEBaWhqysrLw3//93wCAMWPGuN3vo48+wuTJkwEABQUFqKysdN126tQppKamoqqqCpGRkRg/fjwOHz6MyMhI1eERkQExbxBRZykXLJMnT4YQ7R9W3NFtLYqKityuv/nmm6rDIKJuhHmDiDqLcwkRERGR4bFgISIiIsNjwUJERESGx4KFiIiIDM+rrfkvNc1mh2byXIOJEMXGUTa7wiAUWhdLHGjYwtSsMAYAQqWFcoN83WoOUOtbLwLl41VidYvaW1cPUmj7rxArArzXqtreU/75cAbKrzegSbGFvux7WnXKC7qsqLTaBwCh1JpfIdap9jldKLTm1xwKsSp/VhSn9VCZBkTTFXrzO700T40kfsNCREREhseChYiIiAyPBQsREREZHgsWIiIiMjwWLERERGR4LFiIiIjI8FiwEBERkeGxYCEiIiLDY8FCREREhseChYiIiAzPr1rzNw3ojYCAYI9xlh9q1FYcJN/3XDvfrLZuWU6nUrh2Xq2Vv7fITJXgotD23xSo+NZVeA1FkPy69SD5Mdt6KU4JoTB1g6VGob22ItnW3UotvskvqLTbV23N79Dlc4dKrK7Ymh92+XiTXX4b1WKlQ3+Md8rvi5pDYb/V2ZqfiIiIqEMsWIiIiMjwlAuWAwcOYMaMGYiNjYWmadi5c6fb7enp6dA0ze0yffp0j+tdv349Bg0ahODgYCQmJuLTTz9VHRoRGRTzBhF1lnLB0tDQgISEBKxfv77dmOnTp6O0tNR1eeONNzpc5/bt25GRkYGVK1fi2LFjSEhIQHJyMs6cOaM6PCIyIOYNIuos5YNuU1JSkJKS0mGMxWJBdHS09DpfeOEFLF68GAsXLgQAbNy4Ee+++y42bdqERx99VHWIRGQwzBtE1FleOYZl37596Nu3L4YPH44HHngAVVVV7cbabDYcPXoU06ZN+2lQJhOmTZuGQ4cOtXmf5uZm1NbWul2IqHvzdt4AmDuIurMuL1imT5+O1157Dbm5uXj22Wexf/9+pKSkwNnOabmVlZVwOp2IiopyWx4VFYWysrI275OdnQ2r1eq69O/fv6s3g4guoUuRNwDmDqLurMv7sMyfP9/1/2uuuQbx8fEYMmQI9u3bh6lTp3bJY2RmZiIjI8N1vba2lomHqBu7FHkDYO4g6s68flpzXFwc+vTpg/z8/DZv79OnD8xmM8rLy92Wl5eXt/t7tsViQXh4uNuFiPyHN/IGwNxB1J15vWA5deoUqqqqEBMT0+btQUFBGDt2LHJzc13LdF1Hbm4ukpKSvD08IjIg5g0iupBywVJfX4+8vDzk5eUBAAoLC5GXl4fi4mLU19fj4YcfxuHDh1FUVITc3FzMnDkTV111FZKTk13rmDp1KtatW+e6npGRgVdffRVbtmzBv/71LzzwwANoaGhwHf1PRN0b8wYRdZbyMSxHjhzBlClTXNdbfg9OS0vDhg0b8OWXX2LLli2orq5GbGwsbr31Vjz55JOwWH6aR6WgoACVlZWu6/PmzUNFRQVWrFiBsrIyjBkzBnv27Gl1QJ0nQVWNCDB7nnPH2auH0nrN5xqlY0VwkHSs1nBefhBOtTkchENh7iG7TT7Wm/PFmBTmGjHLz+EDAFqQwuuiMO8QIkKlQ802hfUCCKlQeM0VnjpHsNpzV9/f8/xcAODoYL4TI+cNFaqf8Eya/Gtohkqs/H5oVhiDN6nOJeRUmadIaS4htXFAIV7z0vxAJpX5fgCYHAqxdvn3h6byd8ULNCEUZlgzqNraWlitVtw8+mEEmD1PMKf3kP/jBagVLDDL7zheLViaFYoQFizuFAoWXaFgsfcOkR8DAIW/SV4tWGxWuXiHvQlH//o4ampqus2xIS2549z/xSE8zPO+W2ivV1p/iVP+/XHa3ks6tsIh//xWOuTHAABVNvn4szb5D39nm9U+KNY0yxXKAFB3Xj72fINa/tfr5fNBQI38vhVUI7/TWqrV8m5IlXx8j1L5CXuDStpvNdAWR1Gx5xhhxz7sksobnEuIiIiIDI8FCxERERkeCxYiIiIyPBYsREREZHgsWIiIiMjwWLAQERGR4bFgISIiIsNjwUJERESGx4KFiIiIDE+5Nb+RiUAzhEQHVM2h1jVWpZW/+WyDdKxQaQFfq9ZhE7p8C2Vhl+/jLBwKPZ8Br3XG1RQ6CgOAaJLv5mgKD5Mfh03++Qg6I//eAAA9SH73tF0h3+lTlT1EriOn06zY8rwb8uYmqrTxNwpdocWyUG3Nr9Bu3+lUaM3vUMsdmkOh3b5C53qlWMW0q9LKX+nvoWLH9a7Gb1iIiIjI8FiwEBERkeGxYCEiIiLDY8FCREREhseChYiIiAyPBQsREREZHgsWIiIiMjzlguXAgQOYMWMGYmNjoWkadu7c6Xa7pmltXp577rl215mVldUqfsSIEcobQ0TGxLxBRJ2lXLA0NDQgISEB69evb/P20tJSt8umTZugaRruuuuuDtc7atQot/t9/PHHqkMjIoNi3iCizlLudJuSkoKUlJR2b4+Ojna7vmvXLkyZMgVxcXEdDyQgoNV9icg/MG8QUWd5tTV/eXk53n33XWzZssVj7IkTJxAbG4vg4GAkJSUhOzsbAwYMaDO2ubkZzc0/tVqvra0FADRFhiAg0HOLcstZ+TbtAGBqtMkHB3ieGuBiaEFBSvHifJP8uiWmM3Ctt1ntuRMqrfkVphMA1J4Pk8LzJ+x2+RVXnJUO1SwW+fUCQOwV0qEqrbhrB6rt9oHn5dat2btmGgZv5Q2g/dzhLWZ4p5W5N9v4q7Tb1xXa7ausFwCcuny806mwbpVYqLXm1xRSmLdiAcCksC+aVFrzC+9MtSLLqwfdbtmyBWFhYbjzzjs7jEtMTEROTg727NmDDRs2oLCwEBMmTEBdXV2b8dnZ2bBara5L//79vTF8IvIBb+UNgLmDqDvzasGyadMm3H333QgO7vhbj5SUFMyZMwfx8fFITk7G7t27UV1djbfeeqvN+MzMTNTU1LguJSUl3hg+EfmAt/IGwNxB1J157Seh//3f/8Xx48exfft25ftGRERg2LBhyM/Pb/N2i8UCi+rX60RkeN7MGwBzB1F35rVvWP785z9j7NixSEhIUL5vfX09CgoKEBMT44WREZFRMW8QUXuUC5b6+nrk5eUhLy8PAFBYWIi8vDwUFxe7Ympra/H2229j0aJFba5j6tSpWLdunev6smXLsH//fhQVFeHgwYOYPXs2zGYzUlNTVYdHRAbEvEFEnaX8k9CRI0cwZcoU1/WMjAwAQFpaGnJycgAAb775JoQQ7SaOgoICVFZWuq6fOnUKqampqKqqQmRkJMaPH4/Dhw8jMjJSdXhEZEDMG0TUWZoQPj5PqQvU1tbCarUi6dZVPj+tWVM5RUyBVteoFC9qvHO6pt6oOA4vndasBSqe1hzi+X3hEhSotG5Zqqc1OxVOa3b0lB9z9RC1ccie1uy0N+HoW4+jpqYG4eHhSo/hKy2549z/xSE8zPMXzsWOeqX1lzh6SMf+4OglHVvhkH9+z9jVXosKW5h0bGVzT+nYc83yzwUAnG0MkY6ta5Dfv+31irmjTv5zfWCd/CnQQTXysZZzan+me1TI59KQsvPSseZS+TYOAOA49YPnGGHHPuySyhucS4iIiIgMjwULERERGR4LFiIiIjI8FixERERkeF6dS+hSM9l1mITng17toWoHVZpC5J+moLL224JfSGtSmKNIZX4bAFovq3SsXl4hHSucapNamBQONBVO+TmNNLNare2sb1CKl2UOlT/oUPSOUFq3Zpd/rp0Kz7PJoTQMQPZ4v25/+L5vmb30BJoU16saL0tl3iEAECrx3ooF1N7XSuOQD1WeOspP90V+w0JERESGx4KFiIiIDI8FCxERERkeCxYiIiIyPBYsREREZHgsWIiIiMjwWLAQERGR4bFgISIiIsNjwUJERESG5xedboX4sa2fw9HslfWbnPJtA01O+TFoukKnW5VYANDlx6wL+XXrQq3jrkmh82PL6yhDU+xWqQvV9q5yhMJzB4X3BgAIhabCDof8Zw+nTa1bMexyr4vT3gRA7XX0tZax1tbLtRKtc6i1HG1QiG90yL8u553y7+dmxS7ZNpt8vN0m//53NKn9uXGel39P643y69XPK7aNbZLvwO1sUtgPFdKB06a2TzkUumQ7nE3SsUJXy2EOib8XDvwYI5M3/KJgqav7sR3+J//7Hz4eCbUivy+oUcvB3lPrpdhurK6uDlar/NQQvtSSOwZeV+TbgRBd5mTyhia608ehdui6jtOnTyMsLAya9tMn79raWvTv3x8lJSUIDw/34Qi9w9+3D/D/bfSn7RNCoK6uDrGxsTCZusevzcwd3L7uyl+2USVv+MU3LCaTCVdeeWW7t4eHh3frF9QTf98+wP+30V+2r7t8s9KCuYPb1935wzbK5o3u8TGIiIiILmssWIiIiMjw/LpgsVgsWLlyJSwWi6+H4hX+vn2A/2+jv29fd+Xvrwu3r/u7HLbxQn5x0C0RERH5N7/+hoWIiIj8AwsWIiIiMjwWLERERGR4LFiIiIjI8FiwEBERkeH5dcGyfv16DBo0CMHBwUhMTMSnn37q6yF1iaysLGia5nYZMWKEr4fVKQcOHMCMGTMQGxsLTdOwc+dOt9uFEFixYgViYmIQEhKCadOm4cSJE74Z7EXwtH3p6emtXtPp06f7ZrCXOX/NG4D/5Q7mjcsrb/htwbJ9+3ZkZGRg5cqVOHbsGBISEpCcnIwzZ874emhdYtSoUSgtLXVdPv74Y18PqVMaGhqQkJCA9evXt3n76tWr8dJLL2Hjxo345JNP0LNnTyQnJ6OpyVuzK3YtT9sHANOnT3d7Td94441LOEIC/D9vAP6VO5g3LrO8IfzUuHHjxJIlS1zXnU6niI2NFdnZ2T4cVddYuXKlSEhI8PUwvAaAeOedd1zXdV0X0dHR4rnnnnMtq66uFhaLRbzxxhs+GGHnXLh9QgiRlpYmZs6c6ZPx0E/8OW8I4d+5g3nD//nlNyw2mw1Hjx7FtGnTXMtMJhOmTZuGQ4cO+XBkXefEiROIjY1FXFwc7r77bhQXF/t6SF5TWFiIsrIyt9fTarUiMTHRb15PANi3bx/69u2L4cOH44EHHkBVVZWvh3RZuRzyBnD55A7mDf/jlwVLZWUlnE4noqKi3JZHRUWhrKzMR6PqOomJicjJycGePXuwYcMGFBYWYsKECairq/P10Lyi5TXz19cT+PFr3ddeew25ubl49tlnsX//fqSkpMDpdPp6aJcNf88bwOWVO5g3/E+ArwdA6lJSUlz/j4+PR2JiIgYOHIi33noL9913nw9HRhdr/vz5rv9fc801iI+Px5AhQ7Bv3z5MnTrVhyMjf8Lc4V8ut7zhl9+w9OnTB2azGeXl5W7Ly8vLER0d7aNReU9ERASGDRuG/Px8Xw/FK1pes8vl9QSAuLg49OnTx29fUyO63PIG4N+5g3nD//hlwRIUFISxY8ciNzfXtUzXdeTm5iIpKcmHI/OO+vp6FBQUICYmxtdD8YrBgwcjOjra7fWsra3FJ5984pevJwCcOnUKVVVVfvuaGtHlljcA/84dzBv+x29/EsrIyEBaWhquv/56jBs3DmvWrEFDQwMWLlzo66F12rJlyzBjxgwMHDgQp0+fxsqVK2E2m5GamurroV20+vp6t08FhYWFyMvLQ+/evTFgwAD8/ve/x1NPPYWhQ4di8ODBWL58OWJjYzFr1izfDVpBR9vXu3dvrFq1CnfddReio6NRUFCAP/zhD7jqqquQnJzsw1Fffvw5bwD+lzuYNy6zvOHr05S8ae3atWLAgAEiKChIjBs3Thw+fNjXQ+oS8+bNEzExMSIoKEj069dPzJs3T+Tn5/t6WJ3y0UcfCQCtLmlpaUKIH09RXL58uYiKihIWi0VMnTpVHD9+3LeDVtDR9jU2Nopbb71VREZGisDAQDFw4ECxePFiUVZW5uthX5b8NW8I4X+5g3nj8sobmhBCXNoSiYiIiEiNXx7DQkRERP6FBQsREREZHgsWIiIiMjwWLERERGR4LFiIiIjI8FiwEBERkeGxYCEiIiLDY8FCREREhseChYiIiAyPBQsREREZHgsWIiIiMrz/B1zg6QhuIjIwAAAAAElFTkSuQmCC",
      "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).tensor[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": 4,
   "id": "8b27d283",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-09-19T13:35:29.136572Z",
     "start_time": "2024-09-19T13:35:29.134124Z"
    }
   },
   "outputs": [],
   "source": [
    "class NeuralOperatorSolver(AbstractProblem):\n",
    "    input_variables = k_train.full_labels[3]['dof']\n",
    "    output_variables = u_train.full_labels[3]['dof']\n",
    "    conditions = {'data' : Condition(input=k_train, \n",
    "                                     target=u_train)}\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": 5,
   "id": "e34f18b0",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-09-19T13:35:31.245429Z",
     "start_time": "2024-09-19T13:35:29.154937Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "GPU available: False, used: False\n",
      "TPU available: False, using: 0 TPU cores\n",
      "HPU available: False, using: 0 HPUs\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 9: 100%|██████████| 70/70 [00:01<00:00, 69.54it/s, v_num=14, data_loss_step=0.109, train_loss_step=0.109, val_loss_step=0.109, data_loss_epoch=0.105, val_loss_epoch=0.104, train_loss_epoch=0.105]  "
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "`Trainer.fit` stopped: `max_epochs=10` reached.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 9: 100%|██████████| 70/70 [00:01<00:00, 69.13it/s, v_num=14, data_loss_step=0.109, train_loss_step=0.109, val_loss_step=0.109, data_loss_epoch=0.105, val_loss_epoch=0.104, train_loss_epoch=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=10, accelerator='cpu', enable_model_summary=False, batch_size=10) \n",
    "# We train on CPU and avoid model summary at the beginning of training (optional)\n",
    "trainer.train()"
   ]
  },
  {
   "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": 6,
   "id": "0e2a6aa4",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-09-19T13:35:31.295336Z",
     "start_time": "2024-09-19T13:35:31.256308Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final error training 56.26%\n",
      "Final error testing 56.15%\n"
     ]
    }
   ],
   "source": [
    "from pina.loss import LpLoss\n",
    "\n",
    "# make the metric\n",
    "metric_err = LpLoss(relative=True)\n",
    "\n",
    "model = solver.model\n",
    "err = float(metric_err(u_train.squeeze(-1), model(k_train).squeeze(-1)).mean())*100\n",
    "print(f'Final error training {err:.2f}%')\n",
    "\n",
    "err = float(metric_err(u_test.squeeze(-1), model(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 Fourier 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": 7,
   "id": "9af523a5",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-09-19T13:35:44.717807Z",
     "start_time": "2024-09-19T13:35:31.306689Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "GPU available: False, used: False\n",
      "TPU available: False, using: 0 TPU cores\n",
      "HPU available: False, using: 0 HPUs\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 9: 100%|██████████| 70/70 [00:02<00:00, 26.49it/s, v_num=15, data_loss_step=0.00535, train_loss_step=0.00358, val_loss_step=0.00535, data_loss_epoch=0.00372, val_loss_epoch=0.00392, train_loss_epoch=0.00372]"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "`Trainer.fit` stopped: `max_epochs=10` reached.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 9: 100%|██████████| 70/70 [00:02<00:00, 26.33it/s, v_num=15, data_loss_step=0.00535, train_loss_step=0.00358, val_loss_step=0.00535, data_loss_epoch=0.00372, val_loss_epoch=0.00392, train_loss_epoch=0.00372]\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=8,\n",
    "            dimensions=2,\n",
    "            inner_size=24,\n",
    "            padding=8)\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=10, accelerator='cpu', enable_model_summary=False, batch_size=10) # we train on CPU and avoid model summary at beginning of training (optional)\n",
    "trainer.train()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84964cb9",
   "metadata": {},
   "source": [
    "We can clearly see that the final 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, when many data samples are used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "58e2db89",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2024-09-19T13:35:45.259819Z",
     "start_time": "2024-09-19T13:35:44.729042Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Final error training 9.37%\n",
      "Final error testing 9.25%\n"
     ]
    }
   ],
   "source": [
    "model = solver.model\n",
    "\n",
    "err = float(metric_err(u_train.squeeze(-1), model(k_train).squeeze(-1)).mean())*100\n",
    "print(f'Final error training {err:.2f}%')\n",
    "\n",
    "err = float(metric_err(u_test.squeeze(-1), model(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",
   "id": "ba1dfa4b",
   "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": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}