{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "de19422d",
   "metadata": {},
   "source": [
    "# Tutorial: Two dimensional Poisson problem using Extra Features Learning\n",
    "\n",
    "This tutorial presents how to solve with Physics-Informed Neural Networks (PINNs) a 2D Poisson problem with Dirichlet boundary conditions. We will train with standard PINN's training, and with extrafeatures. For more insights on extrafeature learning please read [*An extended physics informed neural network for preliminary analysis of parametric optimal control problems*](https://www.sciencedirect.com/science/article/abs/pii/S0898122123002018).\n",
    "\n",
    "First of all, some useful imports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ad0b8dd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "from torch.nn import Softplus\n",
    "\n",
    "from pina.problem import SpatialProblem\n",
    "from pina.operators import laplacian\n",
    "from pina.model import FeedForward\n",
    "from pina.solvers import PINN\n",
    "from pina.trainer import Trainer\n",
    "from pina.plotter import Plotter\n",
    "from pina.geometry import CartesianDomain\n",
    "from pina.equation import Equation, FixedValue\n",
    "from pina import Condition, LabelTensor\n",
    "from pina.callbacks import MetricTracker"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "492a37b4",
   "metadata": {},
   "source": [
    "## The problem definition"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c0b1777",
   "metadata": {},
   "source": [
    "The two-dimensional Poisson problem is mathematically written as:\n",
    "\\begin{equation}\n",
    "\\begin{cases}\n",
    "\\Delta u = \\sin{(\\pi x)} \\sin{(\\pi y)} \\text{ in } D, \\\\\n",
    "u = 0 \\text{ on } \\Gamma_1 \\cup \\Gamma_2 \\cup \\Gamma_3 \\cup \\Gamma_4,\n",
    "\\end{cases}\n",
    "\\end{equation}\n",
    "where $D$ is a square domain $[0,1]^2$, and $\\Gamma_i$, with $i=1,...,4$, are the boundaries of the square.\n",
    "\n",
    "The Poisson problem is written in **PINA** code as a class. The equations are written as *conditions* that should be satisfied in the corresponding domains. The *truth_solution*\n",
    "is the exact solution which will be compared with the predicted one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "82c24040",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Poisson(SpatialProblem):\n",
    "    output_variables = ['u']\n",
    "    spatial_domain = CartesianDomain({'x': [0, 1], 'y': [0, 1]})\n",
    "\n",
    "    def laplace_equation(input_, output_):\n",
    "        force_term = (torch.sin(input_.extract(['x'])*torch.pi) *\n",
    "                      torch.sin(input_.extract(['y'])*torch.pi))\n",
    "        laplacian_u = laplacian(output_, input_, components=['u'], d=['x', 'y'])\n",
    "        return laplacian_u - force_term\n",
    "\n",
    "    # here we write the problem conditions\n",
    "    conditions = {\n",
    "        'gamma1': Condition(location=CartesianDomain({'x': [0, 1], 'y':  1}), equation=FixedValue(0.)),\n",
    "        'gamma2': Condition(location=CartesianDomain({'x': [0, 1], 'y': 0}), equation=FixedValue(0.)),\n",
    "        'gamma3': Condition(location=CartesianDomain({'x':  1, 'y': [0, 1]}), equation=FixedValue(0.)),\n",
    "        'gamma4': Condition(location=CartesianDomain({'x': 0, 'y': [0, 1]}), equation=FixedValue(0.)),\n",
    "        'D': Condition(location=CartesianDomain({'x': [0, 1], 'y': [0, 1]}), equation=Equation(laplace_equation)),\n",
    "    }\n",
    "\n",
    "    def poisson_sol(self, pts):\n",
    "        return -(\n",
    "            torch.sin(pts.extract(['x'])*torch.pi)*\n",
    "            torch.sin(pts.extract(['y'])*torch.pi)\n",
    "        )/(2*torch.pi**2)\n",
    "    \n",
    "    truth_solution = poisson_sol\n",
    "\n",
    "problem = Poisson()\n",
    "\n",
    "# let's discretise the domain\n",
    "problem.discretise_domain(25, 'grid', locations=['D'])\n",
    "problem.discretise_domain(25, 'grid', locations=['gamma1', 'gamma2', 'gamma3', 'gamma4'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7086c64d",
   "metadata": {},
   "source": [
    "## Solving the problem with standard PINNs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72ba4501",
   "metadata": {},
   "source": [
    "After the problem, the feed-forward neural network is defined, through the class `FeedForward`. This neural network takes as input the coordinates (in this case $x$ and $y$) and provides the unkwown field of the Poisson problem. The residual of the equations are evaluated at several sampling points (which the user can manipulate using the method `CartesianDomain_pts`) and the loss minimized by the neural network is the sum of the residuals.\n",
    "\n",
    "In this tutorial, the neural network is composed by two hidden layers of 10 neurons each, and it is trained for 1000 epochs with a learning rate of 0.006 and $l_2$ weight regularization set to $10^{-7}$. These parameters can be modified as desired. We use the `MetricTracker` class to track the metrics during training."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e7d20d6d",
   "metadata": {
    "scrolled": true
   },
   "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": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 999: : 1it [00:00, 158.53it/s, v_num=3, gamma1_loss=5.29e-5, gamma2_loss=4.09e-5, gamma3_loss=4.73e-5, gamma4_loss=4.18e-5, D_loss=0.00134, mean_loss=0.000304]    "
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "`Trainer.fit` stopped: `max_epochs=1000` reached.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 999: : 1it [00:00, 105.33it/s, v_num=3, gamma1_loss=5.29e-5, gamma2_loss=4.09e-5, gamma3_loss=4.73e-5, gamma4_loss=4.18e-5, D_loss=0.00134, mean_loss=0.000304]\n"
     ]
    }
   ],
   "source": [
    "# make model + solver + trainer\n",
    "model = FeedForward(\n",
    "    layers=[10, 10],\n",
    "    func=Softplus,\n",
    "    output_dimensions=len(problem.output_variables),\n",
    "    input_dimensions=len(problem.input_variables)\n",
    ")\n",
    "pinn = PINN(problem, model, optimizer_kwargs={'lr':0.006, 'weight_decay':1e-8})\n",
    "trainer = Trainer(pinn, max_epochs=1000, callbacks=[MetricTracker()], accelerator='cpu', enable_model_summary=False) # we train on CPU and avoid model summary at beginning of training (optional)\n",
    "\n",
    "# train\n",
    "trainer.train()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb83cc7a",
   "metadata": {},
   "source": [
    "Now the `Plotter` class is used to plot the results.\n",
    "The solution predicted by the neural network is plotted on the left, the exact one is represented at the center and on the right the error between the exact and the predicted solutions is showed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1ab83c03",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 78304 characters]",
      "text/plain": [
       "<Figure size 1600x600 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotter = Plotter()\n",
    "plotter.plot(solver=pinn)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20fdf23e",
   "metadata": {},
   "source": [
    "## Solving the problem with extra-features PINNs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1e76351",
   "metadata": {},
   "source": [
    "Now, the same problem is solved in a different way.\n",
    "A new neural network is now defined, with an additional input variable, named extra-feature, which coincides with the forcing term in the Laplace equation. \n",
    "The set of input variables to the neural network is:\n",
    "\n",
    "\\begin{equation}\n",
    "[x, y, k(x, y)], \\text{ with } k(x, y)=\\sin{(\\pi x)}\\sin{(\\pi y)},\n",
    "\\end{equation}\n",
    "\n",
    "where $x$ and $y$ are the spatial coordinates and $k(x, y)$ is the added feature. \n",
    "\n",
    "This feature is initialized in the class `SinSin`, which needs to be inherited by the `torch.nn.Module` class and to have the `forward` method. After declaring such feature, we can just incorporate in the `FeedForward` class thanks to the `extra_features` argument.\n",
    "**NB**: `extra_features` always needs a `list` as input, you you have one feature just encapsulated it in a class, as in the next cell.\n",
    "\n",
    "Finally, we perform the same training as before: the problem is `Poisson`, the network is composed by the same number of neurons and optimizer parameters are equal to previous test, the only change is the new extra feature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ef3ad372",
   "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": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 999: : 1it [00:00, 111.88it/s, v_num=4, gamma1_loss=2.54e-7, gamma2_loss=2.17e-7, gamma3_loss=1.94e-7, gamma4_loss=2.69e-7, D_loss=9.2e-6, mean_loss=2.03e-6]     "
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "`Trainer.fit` stopped: `max_epochs=1000` reached.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 999: : 1it [00:00, 85.62it/s, v_num=4, gamma1_loss=2.54e-7, gamma2_loss=2.17e-7, gamma3_loss=1.94e-7, gamma4_loss=2.69e-7, D_loss=9.2e-6, mean_loss=2.03e-6] \n"
     ]
    }
   ],
   "source": [
    "class SinSin(torch.nn.Module):\n",
    "    \"\"\"Feature: sin(x)*sin(y)\"\"\"\n",
    "    def __init__(self):\n",
    "        super().__init__()\n",
    "\n",
    "    def forward(self, x):\n",
    "        t = (torch.sin(x.extract(['x'])*torch.pi) *\n",
    "             torch.sin(x.extract(['y'])*torch.pi))\n",
    "        return LabelTensor(t, ['sin(x)sin(y)'])\n",
    "\n",
    "\n",
    "# make model + solver + trainer\n",
    "model_feat = FeedForward(\n",
    "    layers=[10, 10],\n",
    "    func=Softplus,\n",
    "    output_dimensions=len(problem.output_variables),\n",
    "    input_dimensions=len(problem.input_variables)+1\n",
    ")\n",
    "pinn_feat = PINN(problem, model_feat, extra_features=[SinSin()], optimizer_kwargs={'lr':0.006, 'weight_decay':1e-8})\n",
    "trainer_feat = Trainer(pinn_feat, max_epochs=1000, callbacks=[MetricTracker()], accelerator='cpu', enable_model_summary=False) # we train on CPU and avoid model summary at beginning of training (optional)\n",
    "\n",
    "# train\n",
    "trainer_feat.train()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9748a13e",
   "metadata": {},
   "source": [
    "The predicted and exact solutions and the error between them are represented below.\n",
    "We can easily note that now our network, having almost the same condition as before, is able to reach additional order of magnitudes in accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2be6b145",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 68072 characters]",
      "text/plain": [
       "<Figure size 1600x600 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotter.plot(solver=pinn_feat)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7bc0577",
   "metadata": {},
   "source": [
    "## Solving the problem with learnable extra-features PINNs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86c1d7b0",
   "metadata": {},
   "source": [
    "We can still do better!\n",
    "\n",
    "Another way to exploit the  extra features is the addition of learnable parameter inside them.\n",
    "In this way, the added parameters are learned during the training phase of the neural network. In this case, we use:\n",
    "\n",
    "\\begin{equation}\n",
    "k(x, \\mathbf{y}) = \\beta \\sin{(\\alpha x)} \\sin{(\\alpha y)},\n",
    "\\end{equation}\n",
    "\n",
    "where $\\alpha$ and $\\beta$ are the abovementioned parameters.\n",
    "Their implementation is quite trivial: by using the class `torch.nn.Parameter` we cam define all the learnable parameters we need, and they are managed by `autograd` module!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ae8716e7",
   "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": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 999: : 1it [00:00, 119.29it/s, v_num=5, gamma1_loss=3.26e-8, gamma2_loss=7.84e-8, gamma3_loss=1.13e-7, gamma4_loss=3.02e-8, D_loss=2.66e-6, mean_loss=5.82e-7]  "
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "`Trainer.fit` stopped: `max_epochs=1000` reached.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 999: : 1it [00:00, 85.94it/s, v_num=5, gamma1_loss=3.26e-8, gamma2_loss=7.84e-8, gamma3_loss=1.13e-7, gamma4_loss=3.02e-8, D_loss=2.66e-6, mean_loss=5.82e-7] \n"
     ]
    }
   ],
   "source": [
    "class SinSinAB(torch.nn.Module):\n",
    "    \"\"\" \"\"\"\n",
    "    def __init__(self):\n",
    "        super().__init__()\n",
    "        self.alpha = torch.nn.Parameter(torch.tensor([1.0]))\n",
    "        self.beta = torch.nn.Parameter(torch.tensor([1.0]))\n",
    "\n",
    "\n",
    "    def forward(self, x):\n",
    "        t =  (\n",
    "            self.beta*torch.sin(self.alpha*x.extract(['x'])*torch.pi)*\n",
    "                      torch.sin(self.alpha*x.extract(['y'])*torch.pi)\n",
    "        )\n",
    "        return LabelTensor(t, ['b*sin(a*x)sin(a*y)'])\n",
    "\n",
    "\n",
    "# make model + solver + trainer\n",
    "model_lean= FeedForward(\n",
    "    layers=[10, 10],\n",
    "    func=Softplus,\n",
    "    output_dimensions=len(problem.output_variables),\n",
    "    input_dimensions=len(problem.input_variables)+1\n",
    ")\n",
    "pinn_lean = PINN(problem, model_lean, extra_features=[SinSinAB()], optimizer_kwargs={'lr':0.006, 'weight_decay':1e-8})\n",
    "trainer_learn = Trainer(pinn_lean, max_epochs=1000, accelerator='cpu', enable_model_summary=False) # we train on CPU and avoid model summary at beginning of training (optional)\n",
    "\n",
    "# train\n",
    "trainer_learn.train()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0319fb3b",
   "metadata": {},
   "source": [
    "Umh, the final loss is not appreciabily better than previous model (with static extra features), despite the usage of learnable parameters. This is mainly due to the over-parametrization of the network: there are many parameter to optimize during the training, and the model in unable to understand automatically that only the parameters of the extra feature (and not the weights/bias of the FFN) should be tuned in order to fit our problem. A longer training can be helpful, but in this case the faster way to reach machine precision for solving the Poisson problem is removing all the hidden layers in the `FeedForward`, keeping only the $\\alpha$ and $\\beta$ parameters of the extra feature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "daa9cf17",
   "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": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 0: : 0it [00:00, ?it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 999: : 1it [00:00, 131.20it/s, v_num=6, gamma1_loss=2.55e-16, gamma2_loss=4.76e-17, gamma3_loss=2.55e-16, gamma4_loss=4.76e-17, D_loss=1.74e-13, mean_loss=3.5e-14] "
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "`Trainer.fit` stopped: `max_epochs=1000` reached.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 999: : 1it [00:00, 98.81it/s, v_num=6, gamma1_loss=2.55e-16, gamma2_loss=4.76e-17, gamma3_loss=2.55e-16, gamma4_loss=4.76e-17, D_loss=1.74e-13, mean_loss=3.5e-14] \n"
     ]
    }
   ],
   "source": [
    "# make model + solver + trainer\n",
    "model_lean= FeedForward(\n",
    "    layers=[],\n",
    "    func=Softplus,\n",
    "    output_dimensions=len(problem.output_variables),\n",
    "    input_dimensions=len(problem.input_variables)+1\n",
    ")\n",
    "pinn_learn = PINN(problem, model_lean, extra_features=[SinSinAB()], optimizer_kwargs={'lr':0.01, 'weight_decay':1e-8})\n",
    "trainer_learn = Trainer(pinn_learn, max_epochs=1000, callbacks=[MetricTracker()], accelerator='cpu', enable_model_summary=False) # we train on CPU and avoid model summary at beginning of training (optional)\n",
    "\n",
    "# train\n",
    "trainer_learn.train()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "150b3e62",
   "metadata": {},
   "source": [
    "In such a way, the model is able to reach a very high accuracy!\n",
    "Of course, this is a toy problem for understanding the usage of extra features: similar precision could be obtained if the extra features are very similar to the true solution. The analyzed Poisson problem shows a forcing term very close to the solution, resulting in a perfect problem to address with such an approach.\n",
    "\n",
    "We conclude here by showing the graphical comparison of the unknown field and the loss trend for all the test cases presented here: the standard PINN, PINN with extra features, and PINN with learnable extra features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "96e51c43",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 83212 characters]",
      "text/plain": [
       "<Figure size 1600x600 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotter.plot(solver=pinn_learn)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c64fcb4",
   "metadata": {},
   "source": [
    "Let us compare the training losses for the various types of training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2855cea1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 38456 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plotter.plot_loss(trainer, logy=True, label='Standard')\n",
    "plotter.plot_loss(trainer_feat, logy=True,label='Static Features')\n",
    "plotter.plot_loss(trainer_learn, logy=True, label='Learnable Features')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a4c8895",
   "metadata": {},
   "source": [
    "## What's next?\n",
    "\n",
    "Nice you have completed the two dimensional Poisson tutorial of **PINA**! There are multiple directions you can go now:\n",
    "\n",
    "1. Train the network for longer or with different layer sizes and assert the finaly accuracy\n",
    "\n",
    "2. Propose new types of extrafeatures and see how they affect the learning\n",
    "\n",
    "3. Exploit extrafeature training in more complex problems\n",
    "\n",
    "4. Many more..."
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "56be7540488f3dc66429ddf54a0fa9de50124d45fcfccfaf04c4c3886d735a3a"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
