{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial: Building custom geometries with PINA `Location` class\n",
    "\n",
    "In this tutorial we will show how to use geometries in PINA. Specifically, the tutorial will include how to create geometries and how to visualize them. The topics covered are:\n",
    "\n",
    "* Creating CartesianDomains and EllipsoidDomains\n",
    "* Getting the Union and Difference of Geometries\n",
    "* Sampling points in the domain (and visualize them)\n",
    "\n",
    "We import the relevant modules first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from pina.geometry import EllipsoidDomain, Difference, CartesianDomain, Union, SimplexDomain\n",
    "from pina.label_tensor import LabelTensor\n",
    "\n",
    "def plot_scatter(ax, pts, title):\n",
    "    ax.title.set_text(title)\n",
    "    ax.scatter(pts.extract('x'), pts.extract('y'), color='blue', alpha=0.5)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Built-in Geometries"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will create one cartesian and two ellipsoids. For the sake of simplicity, we show here the 2-dimensional, but it's trivial the extension to 3D (and higher) cases. The geometries allows also the generation of samples belonging to the boundary. So, we will create one ellipsoid with the border and one without."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "cartesian = CartesianDomain({'x': [0, 2], 'y': [0, 2]})\n",
    "ellipsoid_no_border = EllipsoidDomain({'x': [1, 3], 'y': [1, 3]})\n",
    "ellipsoid_border = EllipsoidDomain({'x': [2, 4], 'y': [2, 4]}, sample_surface=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `{'x': [0, 2], 'y': [0, 2]}` are the bounds of the `CartesianDomain` being created. \n",
    "\n",
    "To visualize these shapes, we need to sample points on them. We will use the `sample` method of the `CartesianDomain` and `EllipsoidDomain` classes. This method takes a `n` argument which is the number of points to sample. It also takes different modes to sample such as random."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "cartesian_samples = cartesian.sample(n=1000, mode='random')\n",
    "ellipsoid_no_border_samples = ellipsoid_no_border.sample(n=1000, mode='random')\n",
    "ellipsoid_border_samples = ellipsoid_border.sample(n=1000, mode='random')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the samples of each of the geometries to see what we are working with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Cartesian Samples: labels(['x', 'y'])\n",
      "LabelTensor([[[0.2300, 1.6698]],\n",
      "             [[1.7785, 0.4063]],\n",
      "             [[1.5143, 1.8979]],\n",
      "             ...,\n",
      "             [[0.0905, 1.4660]],\n",
      "             [[0.8176, 1.7357]],\n",
      "             [[0.0475, 0.0170]]])\n",
      "Ellipsoid No Border Samples: labels(['x', 'y'])\n",
      "LabelTensor([[[1.9341, 2.0182]],\n",
      "             [[1.5503, 1.8426]],\n",
      "             [[2.0392, 1.7597]],\n",
      "             ...,\n",
      "             [[1.8976, 2.2859]],\n",
      "             [[1.8015, 2.0012]],\n",
      "             [[2.2713, 2.2355]]])\n",
      "Ellipsoid Border Samples: labels(['x', 'y'])\n",
      "LabelTensor([[[3.3413, 3.9400]],\n",
      "             [[3.9573, 2.7108]],\n",
      "             [[3.8341, 2.4484]],\n",
      "             ...,\n",
      "             [[2.7251, 2.0385]],\n",
      "             [[3.8654, 2.4990]],\n",
      "             [[3.2292, 3.9734]]])\n"
     ]
    }
   ],
   "source": [
    "print(f\"Cartesian Samples: {cartesian_samples}\")\n",
    "print(f\"Ellipsoid No Border Samples: {ellipsoid_no_border_samples}\")\n",
    "print(f\"Ellipsoid Border Samples: {ellipsoid_border_samples}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice how these are all `LabelTensor` objects. You can read more about these in the [documentation](https://mathlab.github.io/PINA/_rst/label_tensor.html). At a very high level, they are tensors where each element in a tensor has a label that we can access by doing `<tensor_name>.labels`. We can also access the values of the tensor by doing `<tensor_name>.extract(['x'])`.\n",
    "\n",
    "We are now ready to visualize the samples using matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 209468 characters]",
      "text/plain": [
       "<Figure size 1600x400 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 3, figsize=(16, 4))\n",
    "pts_list = [cartesian_samples, ellipsoid_no_border_samples, ellipsoid_border_samples]\n",
    "title_list = ['Cartesian Domain', 'Ellipsoid Domain', 'Ellipsoid Border Domain']\n",
    "for ax, pts, title in zip(axs, pts_list, title_list):\n",
    "    plot_scatter(ax, pts, title)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have now created, sampled, and visualized our first geometries! We can see that the `EllipsoidDomain` with the border has a border around it. We can also see that the `EllipsoidDomain` without the border is just the ellipse. We can also see that the `CartesianDomain` is just a square."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simplex Domain\n",
    "\n",
    "Among the built-in shapes, we quickly show here the usage of `SimplexDomain`, which can be used for polygonal domains!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 212604 characters]",
      "text/plain": [
       "<Figure size 1600x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import torch\n",
    "spatial_domain = SimplexDomain(\n",
    "                    [\n",
    "                        LabelTensor(torch.tensor([[0, 0]]), labels=[\"x\", \"y\"]),\n",
    "                        LabelTensor(torch.tensor([[1, 1]]), labels=[\"x\", \"y\"]),\n",
    "                        LabelTensor(torch.tensor([[0, 2]]), labels=[\"x\", \"y\"]),\n",
    "                    ]\n",
    "                )\n",
    "\n",
    "spatial_domain2 = SimplexDomain(\n",
    "                    [\n",
    "                        LabelTensor(torch.tensor([[ 0., -2.]]), labels=[\"x\", \"y\"]),\n",
    "                        LabelTensor(torch.tensor([[-.5, -.5]]), labels=[\"x\", \"y\"]),\n",
    "                        LabelTensor(torch.tensor([[-2.,  0.]]), labels=[\"x\", \"y\"]),\n",
    "                    ]\n",
    "                )\n",
    "\n",
    "pts = spatial_domain2.sample(100)\n",
    "fig, axs = plt.subplots(1, 2, figsize=(16, 6))\n",
    "for domain, ax in zip([spatial_domain, spatial_domain2], axs):\n",
    "    pts = domain.sample(1000)\n",
    "    plot_scatter(ax, pts, 'Simplex Domain')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Boolean Operations"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create complex shapes we can use the boolean operations, for example to merge two default geometries. We need to simply use the `Union` class: it takes a list of geometries and returns the union of them.\n",
    "\n",
    "Let's create three unions. Firstly, it will be a union of `cartesian` and `ellipsoid_no_border`. Next, it will be a union of `ellipse_no_border` and `ellipse_border`. Lastly, it will be a union of all three geometries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "cart_ellipse_nb_union = Union([cartesian, ellipsoid_no_border])\n",
    "cart_ellipse_b_union = Union([cartesian, ellipsoid_border])\n",
    "three_domain_union = Union([cartesian, ellipsoid_no_border, ellipsoid_border])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can of course sample points over the new geometries, by using the `sample` method as before. We highlihgt that the available sample strategy here is only *random*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "c_e_nb_u_points = cart_ellipse_nb_union.sample(n=2000, mode='random')\n",
    "c_e_b_u_points = cart_ellipse_b_union.sample(n=2000, mode='random')\n",
    "three_domain_union_points = three_domain_union.sample(n=3000, mode='random')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the samples of each of the unions to see what we are working with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 246768 characters]",
      "text/plain": [
       "<Figure size 1600x400 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 3, figsize=(16, 4))\n",
    "pts_list = [c_e_nb_u_points, c_e_b_u_points, three_domain_union_points]\n",
    "title_list = ['Cartesian with Ellipsoid No Border Union', 'Cartesian with Ellipsoid Border Union', 'Three Domain Union']\n",
    "for ax, pts, title in zip(axs, pts_list, title_list):\n",
    "    plot_scatter(ax, pts, title)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we will find the differences of the geometries. We will find the difference of `cartesian` and `ellipsoid_no_border`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 187572 characters]",
      "text/plain": [
       "<Figure size 800x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cart_ellipse_nb_difference = Difference([cartesian, ellipsoid_no_border])\n",
    "c_e_nb_d_points = cart_ellipse_nb_difference.sample(n=2000, mode='random')\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 6))\n",
    "plot_scatter(ax, c_e_nb_d_points, 'Difference')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create Custom Location"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will take a look on how to create our own geometry. The one we will try to make is a heart defined by the function $$(x^2+y^2-1)^3-x^2y^3 \\le 0$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start by importing what we will need to create our own geometry based on this equation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "from pina import Location\n",
    "from pina import LabelTensor\n",
    "import random"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will create the `Heart(Location)` class and initialize it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Heart(Location):\n",
    "    \"\"\"Implementation of the Heart Domain.\"\"\"\n",
    "\n",
    "    def __init__(self, sample_border=False):\n",
    "        super().__init__()\n",
    "        "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the `Location` class we are inherting from requires both a `sample` method and `is_inside` method, we will create them and just add in \"pass\" for the moment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Heart(Location):\n",
    "    \"\"\"Implementation of the Heart Domain.\"\"\"\n",
    "\n",
    "    def __init__(self, sample_border=False):\n",
    "        super().__init__()\n",
    "\n",
    "    def is_inside(self):\n",
    "        pass\n",
    "\n",
    "    def sample(self):\n",
    "        pass"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have the skeleton for our `Heart` class. The `is_inside` method is where most of the work is done so let's fill it out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "class Heart(Location):\n",
    "    \"\"\"Implementation of the Heart Domain.\"\"\"\n",
    "\n",
    "    def __init__(self, sample_border=False):\n",
    "        super().__init__()\n",
    "\n",
    "    def is_inside(self):\n",
    "        pass\n",
    "\n",
    "    def sample(self, n, mode='random', variables='all'):\n",
    "        sampled_points = []\n",
    "\n",
    "        while len(sampled_points) < n:\n",
    "            x = torch.rand(1)*3.-1.5\n",
    "            y = torch.rand(1)*3.-1.5\n",
    "            if ((x**2 + y**2 - 1)**3 - (x**2)*(y**3)) <= 0:\n",
    "                sampled_points.append([x.item(), y.item()])\n",
    "\n",
    "        return LabelTensor(torch.tensor(sampled_points), labels=['x','y'])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create the Heart geometry we simply run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "heart = Heart()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To sample from the Heart geometry we simply run:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 134356 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pts_heart = heart.sample(1500)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "plot_scatter(ax, pts_heart, 'Heart Domain')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What's next?\n",
    "\n",
    "We have made a very simple tutorial on how to build custom geometries and use domain operation to compose base geometries. Now you can play around with different geometries and build your own!  "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
