notebooks/01_foundations__manifolds.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d62088a0",
   "metadata": {},
   "source": [
    "$\\textbf{Lead Author: Adele Myers}$\n",
    "\n",
    "Inspired by: Guigui, Miolane, Pennec, 2022. Introduction to Riemannian Geometry and Geometric Statistics: from basic theory to implementation with Geomstats."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc0913fc",
   "metadata": {},
   "source": [
    "Note: Before you start to read this notebook, run the following code. This will import packages that will allow later code to run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "ac5b2719",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import warnings\n",
    "\n",
    "sys.path.append(os.path.dirname(os.getcwd()))\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import geomstats.backend as gs\n",
    "\n",
    "gs.random.seed(2020)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff5b438a",
   "metadata": {},
   "source": [
    "# 1. Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f4b7295",
   "metadata": {},
   "source": [
    "Welcome to the `Manifold` class notebook! In this notebook, we will discuss manifolds and the `Manifold` class in Geomstats. Geomstats is a software which analyzes data on manifolds, so naturally the `Manifold` class is quite important.\n",
    "\n",
    "You will learn:\n",
    "\n",
    "1. What is a manifold?\n",
    "2. Why do we care about manifolds?\n",
    "3. What does the `Manifold` class do?\n",
    "4. How is the `Manifold` class structured?\n",
    "    1. What is an open set?\n",
    "    2. What is a closed set?\n",
    "    3. What is a tangent space?\n",
    "        "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48ac2a24",
   "metadata": {},
   "source": [
    "<img src=\"figures/manifold_intro.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4083a61a",
   "metadata": {},
   "source": [
    "# 2. What is a Manifold?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a60222e",
   "metadata": {},
   "source": [
    "$\\textbf{Intuition:}$ When you are first learning, it can be a helpful starting point to think of a manifold as a surface. This surface can have any dimension and any shape as long as it is smooth (in the sense of being continuous and differentiable). (This is not a particularly precise definition, but it can be helpful starting point for building intuition.) For example, a hypersphere is a two dimensional manifold, and we will often use this manifold in examples. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5af7df5",
   "metadata": {},
   "source": [
    "## 2.1 Precise Mathematical Definition"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "486860d6",
   "metadata": {},
   "source": [
    "A manifold is a set of points that satisfy a specific set of constraints. More specifically, a nonempty subset M $\\subseteq \\mathbb{R}^{N}$ is a d-dimensional manifold if and only if ANY of the following conditions hold:\n",
    "\n",
    "(1) (Local parametrization) For every $p \\in M$, there are two open subsets $V \\subseteq \\mathbb{R}^{d}$ and $U \\subseteq \\mathbb{R}^{N}$ with $p \\in U$ and $0 \\in V$. There is also a smooth function $f: V \\to \\mathbb{R}^{N}$ such that $f(0) = p$, where $f$ is a homeomorphism between V and $U \\cap M$, and $f$ is an immersion at 0.\n",
    "\n",
    "(2) (Local implicit function) For every $p \\in M$, there exists an open set $U \\in \\mathbb{R}^{N}$ and a smooth map $f: U \\to \\mathbb{R}^{N-d}$ that is a submersion at p, such that $U \\cap M = f^{-1}$({0}).\n",
    "\n",
    "(3) (Local Graph) For every $x \\in M$, there exists an open neighborhood $U \\subseteq \\mathbb{R}^{N}$ of $x$, a neighborhood $V \\subseteq \\mathbb{R}^{d}$ of 0 and a smooth map $f: V \\to \\mathbb{R}^{N-d}$ such that $U \\cap M = graph(f)^{2}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09f06979",
   "metadata": {},
   "source": [
    "## 2.2 Imprecise \"Layman's Terms\" Definition:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa47546e",
   "metadata": {},
   "source": [
    "A $\\textbf{manifold}$ is a set of points that satisfy a specific set of constraints, and these constraints are given by the mathematical definition above. That means that if a set of points satisfies any one of the conditions, then it is a manifold. However, these conditions can be hard to digest. We can translate the above mathematical conditions into three different ways of describing manifolds:\n",
    "\n",
    "1) (Local Parametrization) a manifold is a topological space that locally resembles Euclidean space near each point. For example, consider a two dimensional grid. We would not be able to deform this grid to have the shape of a sphere under any circumstance, but at each $\\textbf{local}$ point on the manifold, we can approximate the space around the point with a Euclidean grid.\n",
    "\n",
    "2) (Local Implicit Function) a manifold can be understood as the set of points that verify a constraint defined by an implicit equation, given by the function $f$. (see hypersphere example)\n",
    "\n",
    "3) (Local Graph) a manifold can be understood as a d-dimensional surface described by d variables, i.e. by the \"graph\" of a smooth function $f: (x_1,...,x_d) \\to f(x_1,...,x_d)$. This local graph cannot be applied globally (to the whole manifold) because a function $f$ must have only one output per set of inputs, and if we were to try to describe the whole manifold with such a graph, then many manifolds would have more than one output for one set of inputs."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84e49d7c",
   "metadata": {},
   "source": [
    "The figure below provides a visualization of each of these definitions. The first figure is a visualization of \"local parametrization\", the second is a visualization of \"local implicit function\" and the third is a visualization of \"local graph\":"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f1386d3",
   "metadata": {},
   "source": [
    "<img src=\"figures/manifold_definitions1.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf243b6d",
   "metadata": {},
   "source": [
    "<img src=\"figures/manifold_definitions2.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73978259",
   "metadata": {},
   "source": [
    "Each one of these definitions of a manifold is important because each one gives us an additional way to describe a manifold. This is useful because some manifolds are much easier described (or computationally implemented) by one of these definitions than the others. For example, any three of these definitions can be used to describe a hypersphere, but a hypersphere is most easily implemented using definition (2)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1a6966f",
   "metadata": {},
   "source": [
    "## 2.3 Hypersphere example:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e3a40d00",
   "metadata": {},
   "source": [
    "Here, we will consider how to prove that a hypersphere is a manifold using the second condition in the definition of a manifold above (2)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c61b4485",
   "metadata": {},
   "source": [
    "A $\\textbf{hypersphere}$ is any of a set of objects (n-dimensional spheres) resulting from the generalization of a one-dimensional circle and a two-dimensional sphere. The dimension of the manifold is n and is equal to the number of degrees of freedom of a point moving on In n-dimensional space, a hypersphere is the set of all points that are a given distance (called the radius) from a given point (called the center). For example, a 2-dimensional hypersphere in 3-dimensions describes all the points in 3D space that lie on the surface of a sphere. In other words, a hypersphere describes all of the points that lie in three dimensions that have two degrees of freedom $(\\phi, \\theta)$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab26c3ae",
   "metadata": {},
   "source": [
    "$\\textbf{How do we know that a hypersphere is a manifold?}$ Well, we know from the definition of a hypersphere that a hypersphere is the set of all points that are a given distance from the center of your coordinate system. For example in 3-dimensions, the surface $S$ of a sphere with radius 1 can be described by the relation: \n",
    "\n",
    "$|x|^{2} = 1$\n",
    "\n",
    "This simply says that all the points must be a distance of 1 away from the center. (see above figure (b) and take a = 1).\n",
    "\n",
    "We can now define a function \n",
    "\n",
    "$f(x) = |x|^{2} - 1$\n",
    "\n",
    "We know that this function $f$ will always equal zero for all points that lie on the surface $S$ because in order for a point to fall on the surface $S$ the condition $|x|^{2} = 1$ must be true. In other words,\n",
    "\n",
    "$x \\in S \\iff f(x) = 0$\n",
    "\n",
    "which tells us that \n",
    "\n",
    "$x \\in S \\iff x \\in f^{-1}$({0})\n",
    "\n",
    "This last line matches the definition of a manifold: $M = f^{-1}$({0}) where $M$ is the set of points $x$ that satisfy the condition $|x|^{2} = 1$. Therefore, the set of points that satisfy the condition $|x|^{2} = 1$ form a manifold."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6427ab8a",
   "metadata": {},
   "source": [
    "## 2.4 More examples of manifolds:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "deb3926b",
   "metadata": {},
   "source": [
    "A hyperphere is one type of manifold, but there are many other types of manifolds that are commonly used and seen in nature. We will show a few examples here to help you build intuition about what a manifold is and what a manifold can look like."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52f14f02",
   "metadata": {},
   "source": [
    "For example, a torus (the surface of a donut) is a manifold."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2074740",
   "metadata": {},
   "source": [
    "<img src=\"figures/manifold_torus.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c614f981",
   "metadata": {},
   "source": [
    "A klein bottle (the surface of the shape shown below) is also a manifold."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f25a4b02",
   "metadata": {},
   "source": [
    "<img src=\"figures/manifold_klein_bottle.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c07e9b9",
   "metadata": {},
   "source": [
    "Perhaps a more familiar example for those in physics: Minkowski space (or Minkowski spacetime) is a combination of three-dimensional Euclidean space and time into a four-dimensional manifold (where x, y, z, and time are the dimensions of the space). This manifold is theorized to take one of three forms. A 4-dimensional manifold of: (a) flat curvature, (b) positive curvature, or (c) negative curvature. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad68e5ce",
   "metadata": {},
   "source": [
    "<img src=\"figures/manifold_spacetime.jpeg\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1252dfcb",
   "metadata": {},
   "source": [
    "# 3. Why do we care about manifolds?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b70fa3dc",
   "metadata": {},
   "source": [
    "$\\textbf{Manifolds are important because data in nature \"naturally falls on manifolds\"}$, and as we discussed in the introduction section, knowing the manifold that a data set belongs to may give you more predictive power and a better understanding of the data's evolution. \n",
    "\n",
    "$\\textbf{What does it mean for data to \"naturally fall on a manifold\"?}$ In nature, data are often subject to constraints, and these constraints force the data to lie on manifolds. For example, consider the position of cities on the earth."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ec06e21",
   "metadata": {},
   "source": [
    "<img src=\"figures/manifold_cities_on_earth.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab1ca213",
   "metadata": {},
   "source": [
    "The cities are subject to the following constraints: 1) they cannot fly above the surface of the earth because gravity holds them down and 2) they cannot sink down into the earth because the surface of the earth holds them up. Therefore, they are constrained to move (or not move) on the surface of a 3-dimensional sphere. This space that they are confined to exist in should sound familiar-- it is the space of a 2-dimensional hypersphere!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74e0c2fa",
   "metadata": {},
   "source": [
    "Other data falls on manifolds in similar ways: the data is subject to certain constraints, which forces it to fall on a manifold."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d55f5fec",
   "metadata": {},
   "source": [
    "# 4. What does the `Manifold` class do?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cdce00c5",
   "metadata": {},
   "source": [
    "The `Manifold` class describes different types of manifolds. The `Manifold` class and its subclasses implement methods that establish the properties of different types of manifolds, and the manifolds that exist under these subclasses inherit the properties implemented in their respective sublclasses.\n",
    "\n",
    "Note: The words \"class\", \"subclass\", \"methods\" refer to [object oriented programming](https://www.educative.io/blog/object-oriented-programming)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c9888ad8",
   "metadata": {},
   "source": [
    "The manifold class also provides ways of checking whether a vector is $\\textit{tangent}$ to a manifold and ways of projecting a vector to the $\\textit{tangent space}$ of a manifold. In order to better describe what this means, we will now define $\\textit{tangent vectors}$ and $\\textit{tangent spaces}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b80dde22",
   "metadata": {},
   "source": [
    "## 4.1 Tangent Vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "705a6c52",
   "metadata": {},
   "source": [
    "Here we will provide a brief reminder of tangent vectors as a way to lay the groundwork for discussing tangent spaces.\n",
    "\n",
    "A $\\textbf{tangent vector}$ is a vector that is \"tangent\" to a curve or surface at a given point. When a vector is \"tangent\" to a curve, this means that the vector has the same slope as the curve does at that point."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93cf1288",
   "metadata": {},
   "source": [
    "<img src=\"figures/manifold_tangent_vector.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bc23ab0",
   "metadata": {},
   "source": [
    "## 4.2 Tangent Spaces"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17f31482",
   "metadata": {},
   "source": [
    "Now we will introduce the concept of a \"tangent space\". The tangent space at a certain point on a manifold is comprised of all of the possible tangent vectors that exist at that point. For example, if you are considering the tangent space on a curve, then the only possible tangent vectors are tangent vectors that point forward and backward along the line (a). However, if you are considering the tangent space at a point on a surface, then the tangent vectors can point forward, backward, left, right, and everywhere in between, and the set of all tangent vectors forms a plane (b)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5265a8ab",
   "metadata": {},
   "source": [
    "<img src=\"figures/manifold_tangent_space.jpeg\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42a741ef",
   "metadata": {},
   "source": [
    "Thus, the tangent space of a 1-dimensional manifold (curve) is also one dimensional, and the tangent space of a 2-dimensional manifold (a 2-dimensional surface) is also 2-dimensional.\n",
    "\n",
    "Similarly, for every n-dimensional manifold, there exists an n-dimensional tangent space at each point on the manifold, and the tangent space is comprised of all possible tangent vectors on that manifold."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13e555d6",
   "metadata": {},
   "source": [
    "# 5. How is the `Manifold` class structured?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0981e6e3",
   "metadata": {},
   "source": [
    "The hierarchical structure of the classes inheriting from the `Manifold` parent class is as follows (this Figure is a courtesy of Nicolas Guigui):"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "348717a4",
   "metadata": {},
   "source": [
    "<img src=\"figures/manifold_hierarchy.jpeg\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "815af4cb",
   "metadata": {},
   "source": [
    "As discussed in the previous section, one of the primary purposes of the `Manifold` class is to hold information about various types of manifolds. Rules that are universally true for all manifolds are implemented in methods in the parent class `Manifold`. Rules that are true for some types of manifolds are implemented in the subclasses of `Manifold`: `LevelSet`, `OpenSet`, `FiberBundle`, `ProductManifold`, `VectorSpace`, `MatrixLieAlgebra`, and `MatrixLieGroup`. Specific types of manifolds are described in methods within these subclasses.\n",
    "\n",
    "In this notebook, we will focus on describing the subclasses pertenant to the geometry module of geomstats: `LevelSet`, `OpenSet`, `ProductManifold` and `VectorSpace`.\n",
    "\n",
    "In the following subsections, we will discuss the methods and ideas implemented in the parent class and its subclasses."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36dc5f07",
   "metadata": {},
   "source": [
    "## 5.1 The Parent Class: `Manifold`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "578faea8",
   "metadata": {},
   "source": [
    "The Manifold parent class is an abstract base class which provides the minimal skeleton of attributes and methods expected in its subclasses. Note that the methods of the abstract parent class are declared, but they contain no implementation, and they are overridden by the subclasses. The properties that are declared in the `Manifold` class are properties that all types of manifold must posess. For example, the following methods and attributes are implemented in `Manifold`:\n",
    "\n",
    "1. `dim`: $\\textit{attribute}$. the dimension of the manifold. \"How many coordinates are necessary to fully describe the manifold?\"\n",
    "2. `belongs()`: $\\textit{method}$. evaluates whether a given element belongs to that manifold\n",
    "3. `is_tangent()`: $\\textit{method}$. evaluates whether a given vector is a tangent vector at a given point \n",
    "4. `random_point()`: $\\textit{method}$. generates a random point that lies on the manifold"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb8c11e1",
   "metadata": {},
   "source": [
    "While the abstract methods in Manifold do not contain any implementation, the methods of the subclasses of `Manifold`, such as the `Hypersphere` subclass, $\\textit{are}$ implemented and can be run. We will now exemplify this in the following section"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "138bbba4",
   "metadata": {},
   "source": [
    "### 5.1.1 Examples of Using `Manifold`'s Attributes and Methods in the Subclass: `Hypersphere`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c875f005",
   "metadata": {},
   "source": [
    "$\\textbf{Attributes:}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "823b2bce",
   "metadata": {},
   "source": [
    "`dim` : If we build a hypersphere of dimension 2 with the following code, we can check that `sphere.dim` gives back 2. Run the following code to verify this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "23a48d29",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The dimension of the sphere is 2\n"
     ]
    }
   ],
   "source": [
    "from geomstats.geometry.hypersphere import Hypersphere\n",
    "sphere = Hypersphere(dim=2)\n",
    "\n",
    "print(f\"The dimension of the sphere is {sphere.dim}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6325d9ed",
   "metadata": {},
   "source": [
    "$\\textbf{Methods:}$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9bf4f6e",
   "metadata": {},
   "source": [
    "`belongs` : We can re-use the sphere we just built (called \"sphere\") and verify that the point (0, 0, 1) belongs to that sphere (it is the north pole). Run the following code to verify this using Geomstats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ad91b2fa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import geomstats.backend as gs\n",
    "\n",
    "sphere.belongs(gs.array([0, 0, 1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bcfcd2a",
   "metadata": {},
   "source": [
    "`is_tangent` : The vector (1, 1, 0 ) is tangent to the sphere at the north pole, since it does not have a vertical component (last component is equal to 0). Run the following code to verify this using Geomstats."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "f080a042",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sphere.is_tangent(vector=gs.array([1, 1, 0]), base_point=gs.array([0, 0, 1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12764fbf",
   "metadata": {},
   "source": [
    "`random_point` : Now, we will use `random_point` to generate a random point, and then we will use `belongs` to prove that this random point belongs to the sphere."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "d89557c6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from geomstats.geometry.hypersphere import Hypersphere\n",
    "\n",
    "rp=Hypersphere.random_point(sphere)\n",
    "\n",
    "sphere.belongs(rp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc797840",
   "metadata": {},
   "source": [
    "### 5.1.2 The Full `Manifold` Class Code"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fde983a",
   "metadata": {},
   "source": [
    "You can see all of the methods in the `Manifold` parent class by running the following code. Observe the abstract methods denoted with \"@abc.abstractmethod\" that do not contain any implementation, but serve as template for the subclasses to overwrite. The code of `Manifold` can also be found [here](https://github.com/geomstats/geomstats/blob/master/geomstats/geometry/manifold.py)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "5a27e4f9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "class Manifold(abc.ABC):\n",
      "    r\"\"\"Class for manifolds.\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    dim : int\n",
      "        Dimension of the manifold.\n",
      "    shape : tuple of int\n",
      "        Shape of one element of the manifold.\n",
      "        Optional, default : None.\n",
      "    metric : RiemannianMetric\n",
      "        Metric object to use on the manifold.\n",
      "    default_point_type : str, {\\'vector\\', \\'matrix\\'}\n",
      "        Point type.\n",
      "        Optional, default: 'vector'.\n",
      "    default_coords_type : str, {\\'intrinsic\\', \\'extrinsic\\', etc}\n",
      "        Coordinate type.\n",
      "        Optional, default: 'intrinsic'.\n",
      "    \"\"\"\n",
      "\n",
      "    def __init__(\n",
      "        self,\n",
      "        dim,\n",
      "        shape,\n",
      "        metric=None,\n",
      "        default_point_type=None,\n",
      "        default_coords_type=\"intrinsic\",\n",
      "        **kwargs\n",
      "    ):\n",
      "        super(Manifold, self).__init__(**kwargs)\n",
      "        geomstats.errors.check_integer(dim, \"dim\")\n",
      "\n",
      "        if not isinstance(shape, tuple):\n",
      "            raise ValueError(\"Expected a tuple for the shape argument.\")\n",
      "        if default_point_type is None:\n",
      "            default_point_type = POINT_TYPES[len(shape)]\n",
      "\n",
      "        geomstats.errors.check_parameter_accepted_values(\n",
      "            default_point_type, \"default_point_type\", [\"vector\", \"matrix\"]\n",
      "        )\n",
      "\n",
      "        self.dim = dim\n",
      "        self.shape = shape\n",
      "        self.default_point_type = default_point_type\n",
      "        self.default_coords_type = default_coords_type\n",
      "        self.metric = metric\n",
      "\n",
      "    @abc.abstractmethod\n",
      "    def belongs(self, point, atol=gs.atol):\n",
      "        \"\"\"Evaluate if a point belongs to the manifold.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        point : array-like, shape=[..., dim]\n",
      "            Point to evaluate.\n",
      "        atol : float\n",
      "            Absolute tolerance.\n",
      "            Optional, default: backend atol.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        belongs : array-like, shape=[...,]\n",
      "            Boolean evaluating if point belongs to the manifold.\n",
      "        \"\"\"\n",
      "\n",
      "    @abc.abstractmethod\n",
      "    def is_tangent(self, vector, base_point, atol=gs.atol):\n",
      "        \"\"\"Check whether the vector is tangent at base_point.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        vector : array-like, shape=[..., dim]\n",
      "            Vector.\n",
      "        base_point : array-like, shape=[..., dim]\n",
      "            Point on the manifold.\n",
      "        atol : float\n",
      "            Absolute tolerance.\n",
      "            Optional, default: backend atol.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        is_tangent : bool\n",
      "            Boolean denoting if vector is a tangent vector at the base point.\n",
      "        \"\"\"\n",
      "\n",
      "    @abc.abstractmethod\n",
      "    def to_tangent(self, vector, base_point):\n",
      "        \"\"\"Project a vector to a tangent space of the manifold.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        vector : array-like, shape=[..., dim]\n",
      "            Vector.\n",
      "        base_point : array-like, shape=[..., dim]\n",
      "            Point on the manifold.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        tangent_vec : array-like, shape=[..., dim]\n",
      "            Tangent vector at base point.\n",
      "        \"\"\"\n",
      "\n",
      "    @abc.abstractmethod\n",
      "    def random_point(self, n_samples=1, bound=1.0):\n",
      "        \"\"\"Sample random points on the manifold.\n",
      "\n",
      "        If the manifold is compact, a uniform distribution is used.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        n_samples : int\n",
      "            Number of samples.\n",
      "            Optional, default: 1.\n",
      "        bound : float\n",
      "            Bound of the interval in which to sample for non compact manifolds.\n",
      "            Optional, default: 1.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        samples : array-like, shape=[..., {dim, [n, n]}]\n",
      "            Points sampled on the hypersphere.\n",
      "        \"\"\"\n",
      "\n",
      "    def regularize(self, point):\n",
      "        \"\"\"Regularize a point to the canonical representation for the manifold.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        point : array-like, shape=[..., dim]\n",
      "            Point.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        regularized_point : array-like, shape=[..., dim]\n",
      "            Regularized point.\n",
      "        \"\"\"\n",
      "        regularized_point = point\n",
      "        return regularized_point\n",
      "\n",
      "    @property\n",
      "    def metric(self):\n",
      "        \"\"\"Riemannian Metric associated to the Manifold.\"\"\"\n",
      "        return self._metric\n",
      "\n",
      "    @metric.setter\n",
      "    def metric(self, metric):\n",
      "        if metric is not None:\n",
      "            if not isinstance(metric, RiemannianMetric):\n",
      "                raise ValueError(\"The argument must be a RiemannianMetric object\")\n",
      "            if metric.dim != self.dim:\n",
      "                metric.dim = self.dim\n",
      "        self._metric = metric\n"
     ]
    }
   ],
   "source": [
    "import inspect\n",
    "\n",
    "from geomstats.geometry.manifold import Manifold\n",
    "\n",
    "for line in inspect.getsourcelines(Manifold)[0]:\n",
    "    line = line.replace('\\n','')\n",
    "    print(line)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ce847b2",
   "metadata": {},
   "source": [
    "## 5.2 `OpenSet`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bb824a3",
   "metadata": {},
   "source": [
    "Earlier in the notebook, we were able to say that a set of points is a manifold if it satisfied one of three constraints. We also said that every manifold can be described by any three of these definitions, and the choice of defintion is merely a question of which definition is most convenient. `OpenSet` provides a way of describing manifolds with local parametrization, which was labeled (1) on our definition list.\n",
    "\n",
    "One such way to describe a manifold is with the concept of an $\\textbf{Open Set}$: a manifold is the open sets of a d-dimensional vector space, called $\\textbf{ambient space}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6ba6ba0",
   "metadata": {},
   "source": [
    "### 5.2.1 What is an Open Set?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf224618",
   "metadata": {},
   "source": [
    "Intuitively, an open set is a group of numbers that does not include points on the boundary of whatever they are describing. For example, if you were to take the set of all points between a and b but $\\textbf{not}$ include the values a and b, this would be an open set (shown in figure a). If you were to take the set of all points between a and b and include the values a and b, this would be a closed set (shown in figure b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38f8dc93",
   "metadata": {},
   "source": [
    "<img src=\"figures/manifold_openSet_lines.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "822181bc",
   "metadata": {},
   "source": [
    "The above examples showed open and closed sets in one dimension. Similarly, in two dimensions, open sets can be defined as sets which do not contain their boundaries. For example, the inside of the sphere, i.e. the ball without its boundary, is a manifold that is an open set. The image below shows an example of an open set (a) and a closed set (b) in two dimensions."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2991834b",
   "metadata": {},
   "source": [
    "<img src=\"figures/manifold_openSurfaces.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cde0af0",
   "metadata": {},
   "source": [
    "### 5.2.2 What Methods are Implemented in `OpenSet`?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "151d22b8",
   "metadata": {},
   "source": [
    "If we know that a manifold is conveniently described as an open set, then some of the manifold's abstract methods can be rewritten in a specific form. For example, `OpenSet` implements the methods:\n",
    "\n",
    "1. `projection()`: a method to project any d-dimensional vector to the manifold.\n",
    "2. `is_tangent()`: checks whether the input vector is tangent at the input point.\n",
    "3. `to_tangent()`: projects a vector to a tangent space of the manifold.\n",
    "\n",
    "Note that we do not (yet) specify which manifold we are talking about (whether it is a sphere or another surface), we are just saying that we are looking at some manifold that can be described as an open set"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40c86d59",
   "metadata": {},
   "source": [
    "Run the code below to see the contents of the `OpenSet` class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "cd5f0071",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "class OpenSet(Manifold, abc.ABC):\n",
      "    \"\"\"Class for manifolds that are open sets of a vector space.\n",
      "\n",
      "    In this case, tangent vectors are identified with vectors of the ambient\n",
      "    space.\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    dim: int\n",
      "        Dimension of the manifold. It is often the same as the ambient space\n",
      "        dimension but may differ in some cases.\n",
      "    embedding_space: VectorSpace\n",
      "        Ambient space that contains the manifold.\n",
      "    \"\"\"\n",
      "\n",
      "    def __init__(self, dim, embedding_space, **kwargs):\n",
      "        if \"default_point_type\" not in kwargs:\n",
      "            kwargs[\"default_point_type\"] = embedding_space.default_point_type\n",
      "        if \"shape\" not in kwargs:\n",
      "            kwargs[\"shape\"] = embedding_space.shape\n",
      "        super().__init__(dim=dim, **kwargs)\n",
      "        self.embedding_space = embedding_space\n",
      "\n",
      "    def is_tangent(self, vector, base_point=None, atol=gs.atol):\n",
      "        \"\"\"Check whether the vector is tangent at base_point.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        vector : array-like, shape=[..., dim]\n",
      "            Vector.\n",
      "        base_point : array-like, shape=[..., dim]\n",
      "            Point on the manifold.\n",
      "        atol : float\n",
      "            Absolute tolerance.\n",
      "            Optional, default: backend atol.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        is_tangent : bool\n",
      "            Boolean denoting if vector is a tangent vector at the base point.\n",
      "        \"\"\"\n",
      "        return self.embedding_space.belongs(vector, atol)\n",
      "\n",
      "    def to_tangent(self, vector, base_point=None):\n",
      "        \"\"\"Project a vector to a tangent space of the manifold.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        vector : array-like, shape=[..., dim]\n",
      "            Vector.\n",
      "        base_point : array-like, shape=[..., dim]\n",
      "            Point on the manifold.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        tangent_vec : array-like, shape=[..., dim]\n",
      "            Tangent vector at base point.\n",
      "        \"\"\"\n",
      "        return self.embedding_space.projection(vector)\n",
      "\n",
      "    def random_point(self, n_samples=1, bound=1.0):\n",
      "        \"\"\"Sample random points on the manifold.\n",
      "\n",
      "        If the manifold is compact, a uniform distribution is used.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        n_samples : int\n",
      "            Number of samples.\n",
      "            Optional, default: 1.\n",
      "        bound : float\n",
      "            Bound of the interval in which to sample for non compact manifolds.\n",
      "            Optional, default: 1.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        samples : array-like, shape=[..., {dim, [n, n]}]\n",
      "            Points sampled on the hypersphere.\n",
      "        \"\"\"\n",
      "        sample = self.embedding_space.random_point(n_samples, bound)\n",
      "        return self.projection(sample)\n",
      "\n",
      "    @abc.abstractmethod\n",
      "    def projection(self, point):\n",
      "        \"\"\"Project a point in ambient manifold on manifold.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        point : array-like, shape=[..., dim]\n",
      "            Point in ambient manifold.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        projected : array-like, shape=[..., dim]\n",
      "            Projected point.\n",
      "        \"\"\"\n"
     ]
    }
   ],
   "source": [
    "import inspect\n",
    "\n",
    "from geomstats.geometry.base import OpenSet\n",
    "\n",
    "for line in inspect.getsourcelines(OpenSet)[0]:\n",
    "    line = line.replace('\\n','')\n",
    "    print(line)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b35c4d4",
   "metadata": {},
   "source": [
    "## 5.3 `LevelSet`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ba53871",
   "metadata": {},
   "source": [
    "### 5.3.1 What is a Level Set?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2b6b387",
   "metadata": {},
   "source": [
    "Another elementary class of manifolds are $\\textbf{Level Sets}$. A level set is the set of values $x$ for which a function f(x) is equal to a given constant. In other words, a level set is a set of curves for which the fucntion describing a manifold is constant along that curve. \n",
    "\n",
    "In the same way that `OpenSet` is an implementation of the first definition of a manifold (Local Parametrization), `LevelSet` is an implementation of the second definition of a manifold (Local Implicit Function). a level set is a set of points for which the function $f$ takes the exact same value. This value is called the \"level\", and does not need to be a scalar, it could also be a vector.\n",
    "\n",
    "For example, consider a hypersphere in three dimensional space. Each of the concentric spheres is a 2-dimensional manifold, each corresponding to a different level $(r1, r2,r3)$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16bd9d4d",
   "metadata": {},
   "source": [
    "<img src=\"figures/manifold_LevelSet.png\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a93bf48d",
   "metadata": {},
   "source": [
    "You can see [here](https://github.com/geomstats/geomstats/blob/20cccf598e94823ef8d7c0821b2968e6d29660c0/geomstats/geometry/hypersphere.py#L22) that the subclass `Hypersphere` is indeed implemented as a `LevelSet`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8608bea",
   "metadata": {},
   "source": [
    "### 5.3.2 `LevelSet` in Geomstats"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4c2afa0",
   "metadata": {},
   "source": [
    "You can run the code below to see the contents of the `LevelSet` class. The methods of LevelSet are similar to the methods of OpenSet , but their implementation is different. Recall from the (Local Implicit Function) definition, which `LevelSet` is an implementation of, that a manifold of this type must satisfy:\n",
    "\n",
    "\"For every $p \\in M$, there exists an open set $U \\in \\mathbb{R}^{N}$ and a smooth map $f: U \\to \\mathbb{R}^{N-d}$ that is a submersion at p, such that $U \\cap M = f^{-1}$({0}).\"\n",
    "\n",
    "`OpenSet` is an implementation of the first manifold definition (Local Parametrization), and therefore need not follow the (Local Implicit Function) rule above. This is the reason that, `OpenSet` methods and `LevelSet` methods are implemented differently. As an example of these different implementations: observe the implementation of the belongs methods in `LevelSet`. For a general level set, in order to verify if a point belongs to the level set, we should verify that the (Local Implicit Function) definition constraint is met, which is done with the line\n",
    "\n",
    "`constraint = gs.isclose(self.submersion(point), value, atol=atol)`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "48c28386",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "class LevelSet(Manifold, abc.ABC):\n",
      "    \"\"\"Class for manifolds embedded in a vector space by a submersion.\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    dim : int\n",
      "        Dimension of the embedded manifold.\n",
      "    embedding_space : VectorSpace\n",
      "        Embedding space.\n",
      "    default_coords_type : str, {'intrinsic', 'extrinsic', etc}\n",
      "        Coordinate type.\n",
      "        Optional, default: 'intrinsic'.\n",
      "    \"\"\"\n",
      "\n",
      "    def __init__(\n",
      "        self,\n",
      "        dim,\n",
      "        embedding_space,\n",
      "        submersion,\n",
      "        value,\n",
      "        tangent_submersion,\n",
      "        default_coords_type=\"intrinsic\",\n",
      "        **kwargs\n",
      "    ):\n",
      "        if \"shape\" not in kwargs:\n",
      "            kwargs[\"shape\"] = embedding_space.shape\n",
      "        super(LevelSet, self).__init__(\n",
      "            dim=dim,\n",
      "            default_point_type=embedding_space.default_point_type,\n",
      "            default_coords_type=default_coords_type,\n",
      "            **kwargs\n",
      "        )\n",
      "        self.embedding_space = embedding_space\n",
      "        self.embedding_metric = embedding_space.metric\n",
      "        self.submersion = submersion\n",
      "        if isinstance(value, float):\n",
      "            value = gs.array(value)\n",
      "        self.value = value\n",
      "        self.tangent_submersion = tangent_submersion\n",
      "\n",
      "    def belongs(self, point, atol=gs.atol):\n",
      "        \"\"\"Evaluate if a point belongs to the manifold.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        point : array-like, shape=[..., dim]\n",
      "            Point to evaluate.\n",
      "        atol : float\n",
      "            Absolute tolerance.\n",
      "            Optional, default: backend atol.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        belongs : array-like, shape=[...,]\n",
      "            Boolean evaluating if point belongs to the manifold.\n",
      "        \"\"\"\n",
      "        belongs = self.embedding_space.belongs(point, atol)\n",
      "        if not gs.any(belongs):\n",
      "            return belongs\n",
      "        value = self.value\n",
      "        constraint = gs.isclose(self.submersion(point), value, atol=atol)\n",
      "        if value.ndim == 2:\n",
      "            constraint = gs.all(constraint, axis=(-2, -1))\n",
      "        elif value.ndim == 1:\n",
      "            constraint = gs.all(constraint, axis=-1)\n",
      "        return gs.logical_and(belongs, constraint)\n",
      "\n",
      "    def is_tangent(self, vector, base_point, atol=gs.atol):\n",
      "        \"\"\"Check whether the vector is tangent at base_point.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        vector : array-like, shape=[..., dim]\n",
      "            Vector.\n",
      "        base_point : array-like, shape=[..., dim]\n",
      "            Point on the manifold.\n",
      "        atol : float\n",
      "            Absolute tolerance.\n",
      "            Optional, default: backend atol.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        is_tangent : bool\n",
      "            Boolean denoting if vector is a tangent vector at the base point.\n",
      "        \"\"\"\n",
      "        belongs = self.embedding_space.is_tangent(vector, base_point, atol)\n",
      "        tangent_sub_applied = self.tangent_submersion(vector, base_point)\n",
      "        constraint = gs.isclose(tangent_sub_applied, 0.0, atol=atol)\n",
      "        value = self.value\n",
      "        if value.ndim == 2:\n",
      "            constraint = gs.all(constraint, axis=(-2, -1))\n",
      "        elif value.ndim == 1:\n",
      "            constraint = gs.all(constraint, axis=-1)\n",
      "        return gs.logical_and(belongs, constraint)\n",
      "\n",
      "    def intrinsic_to_extrinsic_coords(self, point_intrinsic):\n",
      "        \"\"\"Convert from intrinsic to extrinsic coordinates.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        point_intrinsic : array-like, shape=[..., dim]\n",
      "            Point in the embedded manifold in intrinsic coordinates.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        point_extrinsic : array-like, shape=[..., dim_embedding]\n",
      "            Point in the embedded manifold in extrinsic coordinates.\n",
      "        \"\"\"\n",
      "        raise NotImplementedError(\"intrinsic_to_extrinsic_coords is not implemented.\")\n",
      "\n",
      "    def extrinsic_to_intrinsic_coords(self, point_extrinsic):\n",
      "        \"\"\"Convert from extrinsic to intrinsic coordinates.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        point_extrinsic : array-like, shape=[..., dim_embedding]\n",
      "            Point in the embedded manifold in extrinsic coordinates,\n",
      "            i. e. in the coordinates of the embedding manifold.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        point_intrinsic : array-lie, shape=[..., dim]\n",
      "            Point in the embedded manifold in intrinsic coordinates.\n",
      "        \"\"\"\n",
      "        raise NotImplementedError(\"extrinsic_to_intrinsic_coords is not implemented.\")\n",
      "\n",
      "    @abc.abstractmethod\n",
      "    def projection(self, point):\n",
      "        \"\"\"Project a point in embedding manifold on embedded manifold.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        point : array-like, shape=[..., dim_embedding]\n",
      "            Point in embedding manifold.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        projected : array-like, shape=[..., dim_embedding]\n",
      "            Projected point.\n",
      "        \"\"\"\n",
      "\n",
      "    @abc.abstractmethod\n",
      "    def to_tangent(self, vector, base_point):\n",
      "        \"\"\"Project a vector to a tangent space of the manifold.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        vector : array-like, shape=[..., dim]\n",
      "            Vector.\n",
      "        base_point : array-like, shape=[..., dim]\n",
      "            Point on the manifold.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        tangent_vec : array-like, shape=[..., dim]\n",
      "            Tangent vector at base point.\n",
      "        \"\"\"\n"
     ]
    }
   ],
   "source": [
    "import inspect\n",
    "\n",
    "from geomstats.geometry.base import LevelSet\n",
    "\n",
    "for line in inspect.getsourcelines(LevelSet)[0]:\n",
    "    line = line.replace('\\n','')\n",
    "    print(line)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14edde0c",
   "metadata": {},
   "source": [
    "## 5.4 `VectorSpace`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd5fd5bd",
   "metadata": {},
   "source": [
    "A $\\textit{vector space}$ is a special case of manifold where a vector exists at each point.\n",
    "\n",
    "Many operations defined on manifolds become simpler when the manifold is actually a vector space. Additionally, the ambient spaces used to define manifolds as open or closed sets are often vector spaces. Thus, we provide an implementation `VectorSpace`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "962d4ac2",
   "metadata": {},
   "source": [
    "This class does not provide another way of describing a manifold. This class is an abstract class that makes sure that the ambient/embedding vector space is compatible with the methods that are called in `OpenSet` and `LevelSet`.\n",
    "\n",
    "Actual manifolds are implemented as subclasses of this abstract method and must implement all the abstract methods defined in this class."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "963954ba",
   "metadata": {},
   "source": [
    "You can run the code below to see the contents of the `VectorSpace` class. Observe that, in a `VectorSpace` manifold, the manifold is comprised of vectors -- and more specifically, tangent vectors. Thus, in the implementation, we can see that in the `is_tangent()` method: in order to check that a vector is a tangent to the manifold, we only need to check that the vector belongs to the vector space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "237b9252",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "class VectorSpace(Manifold, abc.ABC):\n",
      "    \"\"\"Abstract class for vector spaces.\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    shape : tuple\n",
      "        Shape of the elements of the vector space. The dimension is the\n",
      "        product of these values by default.\n",
      "    default_point_type : str, {'vector', 'matrix'}\n",
      "        Point type.\n",
      "        Optional, default: 'vector'.\n",
      "    \"\"\"\n",
      "\n",
      "    def __init__(self, shape, **kwargs):\n",
      "        if \"dim\" not in kwargs.keys():\n",
      "            kwargs[\"dim\"] = int(gs.prod(gs.array(shape)))\n",
      "        super(VectorSpace, self).__init__(shape=shape, **kwargs)\n",
      "        self.shape = shape\n",
      "        self._basis = None\n",
      "\n",
      "    def belongs(self, point, atol=gs.atol):\n",
      "        \"\"\"Evaluate if the point belongs to the vector space.\n",
      "\n",
      "        This method checks the shape of the input point.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        point : array-like, shape=[.., {dim, [n, n]}]\n",
      "            Point to test.\n",
      "        atol : float\n",
      "            Unused here.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        belongs : array-like, shape=[...,]\n",
      "            Boolean evaluating if point belongs to the space.\n",
      "        \"\"\"\n",
      "        minimal_ndim = len(self.shape)\n",
      "        belongs = point.shape[-minimal_ndim:] == self.shape\n",
      "        if point.ndim <= minimal_ndim:\n",
      "            return belongs\n",
      "        return gs.tile(gs.array([belongs]), [point.shape[0]])\n",
      "\n",
      "    @staticmethod\n",
      "    def projection(point):\n",
      "        \"\"\"Project a point to the vector space.\n",
      "\n",
      "        This method is for compatibility and returns `point`. `point` should\n",
      "        have the right shape,\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        point: array-like, shape[..., {dim, [n, n]}]\n",
      "            Point.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        point: array-like, shape[..., {dim, [n, n]}]\n",
      "            Point.\n",
      "        \"\"\"\n",
      "        return point\n",
      "\n",
      "    def is_tangent(self, vector, base_point=None, atol=gs.atol):\n",
      "        \"\"\"Check whether the vector is tangent at base_point.\n",
      "\n",
      "        Tangent vectors are identified with points of the vector space so\n",
      "        this checks the shape of the input vector.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        vector : array-like, shape=[..., {dim, [n, n]}]\n",
      "            Vector.\n",
      "        base_point : array-like, shape=[..., {dim, [n, n]}]\n",
      "            Point in the vector space.\n",
      "        atol : float\n",
      "            Absolute tolerance.\n",
      "            Optional, default: backend atol.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        is_tangent : bool\n",
      "            Boolean denoting if vector is a tangent vector at the base point.\n",
      "        \"\"\"\n",
      "        return self.belongs(vector, atol)\n",
      "\n",
      "    def to_tangent(self, vector, base_point=None):\n",
      "        \"\"\"Project a vector to a tangent space of the vector space.\n",
      "\n",
      "        This method is for compatibility and returns vector.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        vector : array-like, shape=[..., {dim, [n, n]}]\n",
      "            Vector.\n",
      "        base_point : array-like, shape=[..., {dim, [n, n]}]\n",
      "            Point in the vector space\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        tangent_vec : array-like, shape=[..., {dim, [n, n]}]\n",
      "            Tangent vector at base point.\n",
      "        \"\"\"\n",
      "        return self.projection(vector)\n",
      "\n",
      "    def random_point(self, n_samples=1, bound=1.0):\n",
      "        \"\"\"Sample in the vector space with a uniform distribution in a box.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        n_samples : int\n",
      "            Number of samples.\n",
      "            Optional, default: 1.\n",
      "        bound : float\n",
      "            Side of hypercube support of the uniform distribution.\n",
      "            Optional, default: 1.0\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        point : array-like, shape=[..., dim]\n",
      "           Sample.\n",
      "        \"\"\"\n",
      "        size = self.shape\n",
      "        if n_samples != 1:\n",
      "            size = (n_samples,) + self.shape\n",
      "        point = bound * (gs.random.rand(*size) - 0.5) * 2\n",
      "        return point\n",
      "\n",
      "    @property\n",
      "    def basis(self):\n",
      "        \"\"\"Basis of the vector space.\"\"\"\n",
      "        if self._basis is None:\n",
      "            self._basis = self._create_basis()\n",
      "\n",
      "        return self._basis\n",
      "\n",
      "    @basis.setter\n",
      "    def basis(self, basis):\n",
      "        if len(basis) < self.dim:\n",
      "            raise ValueError(\n",
      "                \"The basis should have length equal to the \" \"dimension of the space.\"\n",
      "            )\n",
      "        self._basis = basis\n",
      "\n",
      "    @abc.abstractmethod\n",
      "    def _create_basis(self):\n",
      "        \"\"\"Create a canonical basis.\"\"\"\n"
     ]
    }
   ],
   "source": [
    "import inspect\n",
    "\n",
    "from geomstats.geometry.base import VectorSpace\n",
    "\n",
    "for line in inspect.getsourcelines(VectorSpace)[0]:\n",
    "    line = line.replace('\\n','')\n",
    "    print(line)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79923105",
   "metadata": {},
   "source": [
    "## 5.5 `ProductManifold`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "addb9c36",
   "metadata": {},
   "source": [
    "New manifolds can also be created by composing existing manifolds. A product manifold combines two manifolds $M_1$ and $M_2$ to create a new manifold $M_3$, and more generally $N$ manifolds $M1, ..., M_N$ can be combined to create a new manifold $M$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7750cdf",
   "metadata": {},
   "source": [
    "### 5.5.1 Example"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77ac9b35",
   "metadata": {},
   "source": [
    "For example, we can create a product of two spheres as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "fba010bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "from geomstats.geometry.product_manifold import ProductManifold\n",
    "sphere1 = Hypersphere(dim=2)\n",
    "\n",
    "sphere2 = Hypersphere(dim=2)\n",
    "\n",
    "\n",
    "product_of_two_spheres = ProductManifold([sphere1, sphere2])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "645c3d2e",
   "metadata": {},
   "source": [
    "If we generate a random point on this product manifold, we will get 6 coordinates: the first 3 represent a point on the first sphere, and the last 3 represent a point on the second sphere:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "3ffe5cbf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.14573203,  0.15373756,  0.97730596, -0.349178  , -0.48842141,\n",
       "        0.79969947])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "product_of_two_spheres.random_point()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e7557bb",
   "metadata": {},
   "source": [
    "You can run the code below to see the contents of the `ProductManifold` class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "f5009680",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "class ProductManifold(Manifold):\n",
      "    \"\"\"Class for a product of manifolds M_1 x ... x M_n.\n",
      "\n",
      "    In contrast to the classes NFoldManifold, Landmarks, or DiscretizedCurves,\n",
      "    the manifolds M_1, ..., M_n need not be the same, nor of\n",
      "    same dimension, but the list of manifolds needs to be provided.\n",
      "\n",
      "    By default, a point is represented by an array of shape:\n",
      "    [..., dim_1 + ... + dim_n_manifolds]\n",
      "    where n_manifolds is the number of manifolds in the product.\n",
      "    This type of representation is called 'vector'.\n",
      "\n",
      "    Alternatively, a point can be represented by an array of shape:\n",
      "    [..., n_manifolds, dim] if the n_manifolds have same dimension dim.\n",
      "    This type of representation is called `matrix`.\n",
      "\n",
      "    Parameters\n",
      "    ----------\n",
      "    manifolds : list\n",
      "        List of manifolds in the product.\n",
      "    default_point_type : str, {'vector', 'matrix'}\n",
      "        Default representation of points.\n",
      "        Optional, default: 'vector'.\n",
      "    n_jobs : int\n",
      "        Number of jobs for parallel computing.\n",
      "        Optional, default: 1.\n",
      "    \"\"\"\n",
      "\n",
      "    # FIXME (nguigs): This only works for 1d points\n",
      "\n",
      "    def __init__(\n",
      "        self, manifolds, metrics=None, default_point_type=\"vector\", n_jobs=1, **kwargs\n",
      "    ):\n",
      "        geomstats.errors.check_parameter_accepted_values(\n",
      "            default_point_type, \"default_point_type\", [\"vector\", \"matrix\"]\n",
      "        )\n",
      "\n",
      "        self.dims = [manifold.dim for manifold in manifolds]\n",
      "        if metrics is None:\n",
      "            metrics = [manifold.metric for manifold in manifolds]\n",
      "        metric = ProductRiemannianMetric(\n",
      "            metrics, default_point_type=default_point_type, n_jobs=n_jobs\n",
      "        )\n",
      "        dim = sum(self.dims)\n",
      "        shape = (\n",
      "            (dim,) if default_point_type == \"vector\" else (len(manifolds), self.dims[0])\n",
      "        )\n",
      "\n",
      "        super(ProductManifold, self).__init__(\n",
      "            dim=dim,\n",
      "            shape=shape,\n",
      "            metric=metric,\n",
      "            default_point_type=default_point_type,\n",
      "            **kwargs,\n",
      "        )\n",
      "        self.manifolds = manifolds\n",
      "        self.n_jobs = n_jobs\n",
      "\n",
      "    @staticmethod\n",
      "    def _get_method(manifold, method_name, metric_args):\n",
      "        return getattr(manifold, method_name)(**metric_args)\n",
      "\n",
      "    def _iterate_over_manifolds(self, func, args, intrinsic=False):\n",
      "\n",
      "        cum_index = (\n",
      "            gs.cumsum(self.dims)[:-1]\n",
      "            if intrinsic\n",
      "            else gs.cumsum([k + 1 for k in self.dims])\n",
      "        )\n",
      "        arguments = {}\n",
      "        float_args = {}\n",
      "        for key, value in args.items():\n",
      "            if not isinstance(value, float):\n",
      "                arguments[key] = gs.split(value, cum_index, axis=-1)\n",
      "            else:\n",
      "                float_args[key] = value\n",
      "        args_list = [\n",
      "            {key: arguments[key][j] for key in arguments}\n",
      "            for j in range(len(self.manifolds))\n",
      "        ]\n",
      "        pool = joblib.Parallel(n_jobs=self.n_jobs)\n",
      "        out = pool(\n",
      "            joblib.delayed(self._get_method)(\n",
      "                self.manifolds[i], func, {**args_list[i], **float_args}\n",
      "            )\n",
      "            for i in range(len(self.manifolds))\n",
      "        )\n",
      "        return out\n",
      "\n",
      "    def belongs(self, point, atol=gs.atol):\n",
      "        \"\"\"Test if a point belongs to the manifold.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        point : array-like, shape=[..., {dim, [n_manifolds, dim_each]}]\n",
      "            Point.\n",
      "        atol : float,\n",
      "            Tolerance.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        belongs : array-like, shape=[...,]\n",
      "            Boolean evaluating if the point belongs to the manifold.\n",
      "        \"\"\"\n",
      "        point_type = self.default_point_type\n",
      "\n",
      "        if point_type == \"vector\":\n",
      "            intrinsic = self.metric.is_intrinsic(point)\n",
      "            belongs = self._iterate_over_manifolds(\n",
      "                \"belongs\", {\"point\": point, \"atol\": atol}, intrinsic\n",
      "            )\n",
      "            belongs = gs.stack(belongs, axis=-1)\n",
      "\n",
      "        else:\n",
      "            belongs = gs.stack(\n",
      "                [\n",
      "                    space.belongs(point[..., i, :], atol)\n",
      "                    for i, space in enumerate(self.manifolds)\n",
      "                ],\n",
      "                axis=-1,\n",
      "            )\n",
      "\n",
      "        belongs = gs.all(belongs, axis=-1)\n",
      "        return belongs\n",
      "\n",
      "    def regularize(self, point):\n",
      "        \"\"\"Regularize the point into the manifold's canonical representation.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        point : array-like, shape=[..., {dim, [n_manifolds, dim_each]}]\n",
      "            Point to be regularized.\n",
      "        point_type : str, {'vector', 'matrix'}\n",
      "            Representation of point.\n",
      "            Optional, default: None.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        regularized_point : array-like,\n",
      "            shape=[..., {dim, [n_manifolds, dim_each]}]\n",
      "            Point in the manifold's canonical representation.\n",
      "        \"\"\"\n",
      "        point_type = self.default_point_type\n",
      "\n",
      "        if point_type == \"vector\":\n",
      "            intrinsic = self.metric.is_intrinsic(point)\n",
      "            regularized_point = self._iterate_over_manifolds(\n",
      "                \"regularize\", {\"point\": point}, intrinsic\n",
      "            )\n",
      "            regularized_point = gs.concatenate(regularized_point, axis=-1)\n",
      "        elif point_type == \"matrix\":\n",
      "            regularized_point = [\n",
      "                manifold_i.regularize(point[..., i, :])\n",
      "                for i, manifold_i in enumerate(self.manifolds)\n",
      "            ]\n",
      "            regularized_point = gs.stack(regularized_point, axis=1)\n",
      "        return regularized_point\n",
      "\n",
      "    def random_point(self, n_samples=1, bound=1.0):\n",
      "        \"\"\"Sample in the product space from the uniform distribution.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        n_samples : int, optional\n",
      "            Number of samples.\n",
      "        bound : float\n",
      "            Bound of the interval in which to sample for non compact manifolds.\n",
      "            Optional, default: 1.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        samples : array-like, shape=[..., {dim, [n_manifolds, dim_each]}]\n",
      "            Points sampled on the hypersphere.\n",
      "        \"\"\"\n",
      "        point_type = self.default_point_type\n",
      "        geomstats.errors.check_parameter_accepted_values(\n",
      "            point_type, \"point_type\", [\"vector\", \"matrix\"]\n",
      "        )\n",
      "\n",
      "        if point_type == \"vector\":\n",
      "            data = self.manifolds[0].random_point(n_samples, bound)\n",
      "            if len(self.manifolds) > 1:\n",
      "                for space in self.manifolds[1:]:\n",
      "                    samples = space.random_point(n_samples, bound)\n",
      "                    data = gs.concatenate([data, samples], axis=-1)\n",
      "            return data\n",
      "\n",
      "        point = [space.random_point(n_samples, bound) for space in self.manifolds]\n",
      "        samples = gs.stack(point, axis=-2)\n",
      "        return samples\n",
      "\n",
      "    def projection(self, point):\n",
      "        \"\"\"Project a point in product embedding manifold on each manifold.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        point : array-like, shape=[..., {dim, [n_manifolds, dim_each]}]\n",
      "            Point in embedding manifold.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        projected : array-like, shape=[..., {dim, [n_manifolds, dim_each]}]\n",
      "            Projected point.\n",
      "        \"\"\"\n",
      "        point_type = self.default_point_type\n",
      "        geomstats.errors.check_parameter_accepted_values(\n",
      "            point_type, \"point_type\", [\"vector\", \"matrix\"]\n",
      "        )\n",
      "\n",
      "        if point_type == \"vector\":\n",
      "            intrinsic = self.metric.is_intrinsic(point)\n",
      "            projected_point = self._iterate_over_manifolds(\n",
      "                \"projection\", {\"point\": point}, intrinsic\n",
      "            )\n",
      "            projected_point = gs.concatenate(projected_point, axis=-1)\n",
      "        elif point_type == \"matrix\":\n",
      "            projected_point = [\n",
      "                manifold_i.projection(point[..., i, :])\n",
      "                for i, manifold_i in enumerate(self.manifolds)\n",
      "            ]\n",
      "            projected_point = gs.stack(projected_point, axis=-2)\n",
      "        return projected_point\n",
      "\n",
      "    def to_tangent(self, vector, base_point):\n",
      "        \"\"\"Project a vector to a tangent space of the manifold.\n",
      "\n",
      "        The tangent space of the product manifold is the direct sum of\n",
      "        tangent spaces.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        vector : array-like, shape=[..., dim]\n",
      "            Vector.\n",
      "        base_point : array-like, shape=[..., dim]\n",
      "            Point on the manifold.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        tangent_vec : array-like, shape=[..., dim]\n",
      "            Tangent vector at base point.\n",
      "        \"\"\"\n",
      "        point_type = self.default_point_type\n",
      "        geomstats.errors.check_parameter_accepted_values(\n",
      "            point_type, \"point_type\", [\"vector\", \"matrix\"]\n",
      "        )\n",
      "\n",
      "        if point_type == \"vector\":\n",
      "            intrinsic = self.metric.is_intrinsic(base_point)\n",
      "            tangent_vec = self._iterate_over_manifolds(\n",
      "                \"to_tangent\", {\"base_point\": base_point, \"vector\": vector}, intrinsic\n",
      "            )\n",
      "            tangent_vec = gs.concatenate(tangent_vec, axis=-1)\n",
      "        elif point_type == \"matrix\":\n",
      "            tangent_vec = [\n",
      "                manifold_i.to_tangent(vector[..., i, :], base_point[..., i, :])\n",
      "                for i, manifold_i in enumerate(self.manifolds)\n",
      "            ]\n",
      "            tangent_vec = gs.stack(tangent_vec, axis=-2)\n",
      "        return tangent_vec\n",
      "\n",
      "    def is_tangent(self, vector, base_point, atol=gs.atol):\n",
      "        \"\"\"Check whether the vector is tangent at base_point.\n",
      "\n",
      "        The tangent space of the product manifold is the direct sum of\n",
      "        tangent spaces.\n",
      "\n",
      "        Parameters\n",
      "        ----------\n",
      "        vector : array-like, shape=[..., dim]\n",
      "            Vector.\n",
      "        base_point : array-like, shape=[..., dim]\n",
      "            Point on the manifold.\n",
      "        atol : float\n",
      "            Absolute tolerance.\n",
      "            Optional, default: backend atol.\n",
      "\n",
      "        Returns\n",
      "        -------\n",
      "        is_tangent : bool\n",
      "            Boolean denoting if vector is a tangent vector at the base point.\n",
      "        \"\"\"\n",
      "        point_type = self.default_point_type\n",
      "        geomstats.errors.check_parameter_accepted_values(\n",
      "            point_type, \"point_type\", [\"vector\", \"matrix\"]\n",
      "        )\n",
      "\n",
      "        if point_type == \"vector\":\n",
      "            intrinsic = self.metric.is_intrinsic(base_point)\n",
      "            is_tangent = self._iterate_over_manifolds(\n",
      "                \"is_tangent\",\n",
      "                {\"base_point\": base_point, \"vector\": vector, \"atol\": atol},\n",
      "                intrinsic,\n",
      "            )\n",
      "            is_tangent = gs.stack(is_tangent, axis=-1)\n",
      "\n",
      "        else:\n",
      "            is_tangent = gs.stack(\n",
      "                [\n",
      "                    space.is_tangent(\n",
      "                        vector[..., i, :], base_point[..., i, :], atol=atol\n",
      "                    )\n",
      "                    for i, space in enumerate(self.manifolds)\n",
      "                ],\n",
      "                axis=-1,\n",
      "            )\n",
      "\n",
      "        is_tangent = gs.all(is_tangent, axis=-1)\n",
      "        return is_tangent\n"
     ]
    }
   ],
   "source": [
    "import inspect\n",
    "\n",
    "from geomstats.geometry.product_manifold import ProductManifold\n",
    "\n",
    "for line in inspect.getsourcelines(ProductManifold)[0]:\n",
    "    line = line.replace('\\n','')\n",
    "    print(line)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e913791",
   "metadata": {},
   "source": [
    "# Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ff91bfd",
   "metadata": {},
   "source": [
    "Key takeaways from this notebook:\n",
    "\n",
    "1. There are three ways to define a manifold, and each of these three definitions provide a different option for manifold implementation.\n",
    "2. Learning to work with manifolds is useful because many datasets naturally lie on manifolds.\n",
    "3. The `Manifold` class stores information about many types of manifolds. Each of the subclasses implements a different $\\textit{type}$ of manifold (for example, `LevelSet`), and specific manifolds (for example, `Hypersphere`) are implemented within these subclasses."
   ]
  }
 ],
 "metadata": {
  "backends": [
   "numpy",
   "autograd"
  ],
  "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}