#!/usr/bin/env python
# coding: utf-8

# # Smooth manifolds, charts and scalar fields
# 
# This notebook accompanies the lecture
# [Symbolic tensor calculus on manifolds](http://sagemanifolds.obspm.fr/jncf2018/) at JNCF 2018.
# 
# Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Worksheets/JNCF2018/jncf18_scalar.ipynb) to download the notebook file (ipynb format). To run it, you must start SageMath within the Jupyter Notebook server, via the command `sage -n jupyter`

# *NB:* a version of SageMath at least equal to 7.5 (8.2 for the SymPy part) is required to run this notebook.

# First we set up the notebook to use LaTeX display:

# In[1]:


get_ipython().run_line_magic('display', 'latex')


# ## Starting with a manifold
# 
# Manifolds are constructed via the global function `Manifold`:

# In[2]:


M = Manifold(2, 'M')
print(M)


# By default, `Manifold` returns a manifold over $\mathbb{R}$:

# In[3]:


M.base_field()


# Other base fields must be specified with the optional keyword `field`, like
# ```python
# M = Manifold(2, 'M', field='complex')
# ```
# We may check that $M$ is a topological space:

# In[4]:


M in Sets().Topological()


# Actually, $M$ belongs to the following categories:

# In[5]:


M.categories()


# As we can see, by default, `Manifold` constructs a smooth manifold.
# If one would like to stick to the topological level, one should write
# ```python
# M = Manifold(2, 'M', structure='topological')
# ```
# Smooth manifolds are implemented by the Python class `DifferentiableManifold`:

# In[6]:


type(M)


# The type of `M` appears as `DifferentiableManifold-with-category` because it is actually a subclass of `DifferentiableManifold`, which is dynamically generated by SageMath's [category mechanism](http://doc.sagemath.org/html/en/reference/categories/sage/categories/primer.html):

# In[7]:


isinstance(M, 
           sage.manifolds.differentiable.manifold.DifferentiableManifold)


# The class `DifferentiableManifold` inherits from `TopologicalManifold`:

# In[8]:


isinstance(M, sage.manifolds.manifold.TopologicalManifold)


# ## Coordinate charts
# 
# We declare a chart, along with the symbols used to denote the coordinates (here $x=x^0$ and $y=x^1$) by

# In[9]:


U = M.open_subset('U')
XU.<x,y> = U.chart()
XU


# Open subsets are implemented by a (dynamically generated) subclass of
# `DifferentiableManifold`, since they are manifolds in their own:

# In[10]:


isinstance(U, 
           sage.manifolds.differentiable.manifold.DifferentiableManifold)


# Points on $M$ are created by passing their coordinates in a given chart:

# In[11]:


p = U((1,2), chart=XU, name='p')
print(p)


# The syntax `U(...)` used to create $p$ as an element of $U$ reflects the 
# **parent/element pattern** employed in SageMath; indeed $U$ is the parent of $p$:

# In[12]:


p.parent()


# Points are implemented by the class `ManifoldPoint`.
# The principal attribute of this class is a Python dictionary storing the point's coordinates
# in various charts:

# In[13]:


p._coordinates


# Of course, we can recover the point's coordinates by letting the chart act on the point:

# In[14]:


XU(p)


# Let us  introduce a second chart on $M$:

# In[15]:


V = M.open_subset('V')
XV.<xp,yp> = V.chart("xp:x' yp:y'")
XV


# and declare that $M$ is covered by only two charts, i.e. that $M=U\cup V$:

# In[16]:


M.declare_union(U,V)


# In[17]:


M.atlas()


# ### Transition map
# 
# We define the transition map `XU` $\to$ `XV` on $W = U\cap V$ as follows:

# In[18]:


XU_to_XV = XU.transition_map(XV, 
                             (x/(x^2+y^2), y/(x^2+y^2)), 
                             intersection_name='W',
                             restrictions1= x^2+y^2!=0, 
                             restrictions2= xp^2+yp^2!=0)
XU_to_XV.display()


# The value of the argument `restrictions1` means that
# $W = U\setminus \{S\}$, where $S$ is the point of coordinates $(x,y)=(0,0)$,
# while the value of `restrictions2` means that
# $W = V\setminus \{N\}$, where $N$ is the point of coordinates $(x',y')=(0,0)$.
# Since $M=U\cup V$, we have then
# $$
#     U = M \setminus \{N\},\qquad
#     V = M \setminus \{S\},\quad\mbox{and}\quad
#     W = M \setminus \{N, S\} .
# $$
# 
# The transition map `XV` $\to$ `XU`  is obtained by computing the inverse of the one defined above:

# In[19]:


XU_to_XV.inverse().display()


# At this stage, the smooth manifold $M$ is fully specified, being covered by
# one atlas with all transition maps specified. One may have recognized that
# $M$ is nothing but the **2-dimensional sphere**:
# $$
#     M = \mathbb{S}^2 ,
# $$
# with `XU` (resp. `XV`) being
# the chart of \defin{stereographic coordinates}\index{stereographic!coordinates}
# from the North pole $N$ (resp. the South pole $S$).
# 
# Since the transition maps have been defined,
# we can ask for the coordinates $(x',y')$ of the point $p$:

# In[20]:


XV(p)


# This operation has updated the dictionary `_coordinates`:

# In[21]:


p._coordinates


# ## Smooth maps
# 
# As a example of smooth map, let us consider the canonical embedding of $M=\mathbb{S}^2$ 
# into $\mathbb{R}^3$. We need first to introduce $\mathbb{R}^3$ as a 3-dimensional smooth
# manifold, endowed with a single chart:

# In[22]:


R3 = Manifold(3, 'R^3', r'\mathbb{R}^3')
XR3.<X,Y,Z> = R3.chart()
XR3


# The embedding $\Phi: \mathbb{S}^2 \to \mathbb{R}^3$ is then defined in terms of its coordinate expression in the two charts covering $M=\mathbb{S}^2$:

# In[23]:


Phi = M.diff_map(R3, {(XU, XR3): 
                       [2*x/(1+x^2+y^2), 2*y/(1+x^2+y^2),
                        (x^2+y^2-1)/(1+x^2+y^2)],
                      (XV, XR3): 
                       [2*xp/(1+xp^2+yp^2), 2*yp/(1+xp^2+yp^2),
                        (1-xp^2-yp^2)/(1+xp^2+yp^2)]},
                  name='Phi', latex_name=r'\Phi')
Phi.display()


# We may use $\Phi$ for graphical purposes, for instance to display the grids of the stereographic charts `XU` (in red) and `XV` (in green), with the point $p$ atop:

# In[24]:


graph = XU.plot(chart=XR3, mapping=Phi, number_values=25, 
                label_axes=False) + \
        XV.plot(chart=XR3, mapping=Phi, number_values=25, 
                color='green', label_axes=False) + \
        p.plot(chart=XR3, mapping=Phi, label_offset=0.05)
show(graph, viewer='threejs', online=True)


# ## Scalar fields
# 
# Let us define a scalar field, i.e. a smooth map $M\to \mathbb{R}$:

# In[25]:


f = M.scalar_field({XU: 1/(1+x^2+y^2), XV: (xp^2+yp^2)/(1+xp^2+yp^2)},
                   name='f')
f.display()


# Scalar fields are implemented by the class `DiffScalarField`
# and the function chart representations are
# stored in the attribute `_express` of this class, which is a
# Python dictionary whose keys are the various charts defined on $M$:

# In[26]:


f._express


# One may wonder about the compatibility of the two coordinate expressions
# provided in the definition of $f$. Actually, to enforce the compatibility, it is possible to declare the scalar field in a single chart, `XU` say,
# and then to obtain its expression in chart `XV` by analytic continuation
# from the expression in $W=U\cap V$, where both expressions are known, thanks
# to the transition map `XV` $\to$ `XU`:

# In[27]:


f0 = M.scalar_field({XU: 1/(1+x^2+y^2)})
f0.add_expr_by_continuation(XV, U.intersection(V))
f == f0


# The representation of the scalar field in a given chart, i.e. the public access to the directory `_express`, is obtained via the method `coord_function()`:

# In[28]:


fU = f.coord_function(XU)
fU.display()


# In[29]:


fV = f.coord_function(XV)
fV.display()


# Both `fU` and `fV` are instances of the class `ChartFunction`:

# In[30]:


isinstance(fU, sage.manifolds.chart_func.ChartFunction)


# Mathematically, **chart functions** are real-valued functions on the codomain of the considered chart. They map coordinates to elements of the base field (here $\mathbb{R}$):

# In[31]:


fU(1,2)


# In[32]:


fU(*XU(p))


# while scalar fields maps manifold points to $\mathbb{R}$:

# In[33]:


f(p)


# Internally, each chart function stores coordinate expressions with respect to various computational engines:
# 
# - SageMath symbolic engine, based on the [Pynac](http://pynac.org) backend, 
#   with   Maxima used for some simplifications or computation of integrals;
# - [SymPy](http://www.sympy.org}) (Python library for symbolic mathematics);
# - in the future, more symbolic engines (e.g. Giac) or numerical ones will be implemented
# 
# The coordinate expressions are stored in the dictionary `_express`,
# whose keys are strings identifying the computational engines. By default
# only SageMath symbolic expressions, i.e. expressions pertaining
# to the so-called Symbolic Ring (`SR`), are stored:

# In[34]:


fU._express


# The public access to the dictionary `_express` is performed via the
# method `expr()`:

# In[35]:


fU.expr()


# In[36]:


type(fU.expr())


# Actually, `fU.expr()` is a shortcut for `fU.expr('SR')` since
# `SR` is the default symbolic engine. Note that the class
# `Expression` is that devoted to SageMath symbolic expressions.
# The method `expr()` can also be invoked to get the expression in
# another symbolic engine, for instance SymPy:

# In[37]:


fU.expr('sympy')


# In[38]:


type(fU.expr('sympy'))


# This operation has updated the dictionary `_express`:

# In[39]:


fU._express


# The default calculus engine for chart functions of chart `XU` can
# changed thanks to the method `set_calculus_method()`:

# In[40]:


XU.set_calculus_method('sympy')
fU.expr()


# We can have better (i.e. LaTeX) rendering of SymPy objects with
# ```python
# from sympy import init_printing
# init_printing()
# ```
# We don't use it here to make clear the distinction between `SR` objects and SymPy ones.

# Reverting to SageMath's symbolic engine:

# In[41]:


XU.set_calculus_method('SR')
fU.expr()


# Symbolic expressions can be accessed directly from the scalar field,
# `f.expr(XU)` being a shortcut for `f.coord_function(XU).expr()`:

# In[42]:


f.expr(XU)


# In[43]:


f.expr(XV)


# ## Algebra of scalar fields
# 
# The commutative algebra of scalar fields on $M$ is

# In[44]:


CM = M.scalar_field_algebra()
CM


# In[45]:


CM.category()


# As for the manifold classes, the actual Python class implementing
# $C^\infty(M)$ is inherited from `DiffScalarFieldAlgebra` via the
# [category mechanism](http://doc.sagemath.org/html/en/reference/categories/sage/categories/primer.html), hence it bares the name `DiffScalarFieldAlgebra-with-category`:

# In[46]:


type(CM)


# The class `DiffScalarFieldAlgebra-with-category` is dynamically generated as a subclass of `DiffScalarFieldAlgebra` with extra functionalities, like
# for instance the method `is_commutative()`:

# In[47]:


CM.is_commutative()


# Let us have a look at the code of this method; it suffices to use the double question mark:

# In[48]:


get_ipython().run_line_magic('pinfo2', 'CM.is_commutative')


# We see from the `File` field that the code belongs to the category part of SageMath, not to the manifold part, where the class `DiffScalarFieldAlgebra` is defined.
# 
# Regarding the scalar field $f$ introduced above, we have of course

# In[49]:


f in CM


# Actually, in Sage language, `CM`=$C^\infty(M)$ is the parent of `f`:

# In[50]:


f.parent() is CM


# The zero element of the algebra $C^\infty(M)$ is

# In[51]:


CM.zero().display()


# while its unit element is

# In[52]:


CM.one().display()


# ## Implementation of algebraic operations on scalar fields
# 
# Let us consider some operation in the algebra $C^\infty(M)$:

# In[53]:


h = f + 2*CM.one()
h.display()


# In[54]:


h(p)


# We are going to investigate how the above addition is performed. For the Python interpreter
# `h = f + 2*CM.one()` is equivalent to `h = f.__add__(2*CM.one())`,
# i.e. the `+` operator amounts to calling the method `__add__()` on `f`.
# To have a look at the source code of this method, we use the double question mark:

# In[55]:


get_ipython().run_line_magic('pinfo2', 'f.__add__')


# We see that the method `__add__()` is implemented at the level
# of the class `Element` from which `DiffScalarField` inherits, via
# `CommutativeAlgebraElement`.
# In the present case, `left` = `f` and `right` = `2*CM.one()`
# have the same parent, namely the algebra `CM`, so that the actual
# result is computed in the line
# ```python
# return (<Element>left)._add_(right)
# ```
# This invokes the method `_add_()`
# (note the single underscore on each side of `add`). This operator is
# implemented at the level of `ScalarField`, as it can be checked from the source code:

# In[56]:


get_ipython().run_line_magic('pinfo2', 'f._add_')


# This reflects a [general strategy](http://doc.sagemath.org/html/en/thematic_tutorials/coercion_and_categories.html) in SageMath: the arithmetic Python operators `__add__`, `__sub__`,  etc. are implemented at the
# top-level class `Element`, while the specific element subclasses,
# like `ScalarField` here, implement single-underscore methods `_add_`, `_sub_`, etc., which perform the actual computation
# when both operands have the same parent.
# 
# Looking at the code, we notice that the first step is to search
# for the charts in which both operands of the addition operator have a coordinate
# expression. This is performed by the method `common_charts()`;
# in the current case, we get the two stereographic charts defined on $M$:

# In[57]:


f.common_charts(2*CM.one())


# In general, `common_charts()` returns the charts for which both operands
# have already a known coordinate expression or for which a coordinate
# expression can be computed by a known transition map, as we can see on the source code:

# In[58]:


get_ipython().run_line_magic('pinfo2', 'f.common_charts')


# Once the list of charts in which both operands have a coordinate expression,
# the addition is performed at the chart function level.
# The code for the addition of chart functions defined on the same chart
# is (recall that `fU` is the chart function representing $f$ in chart `XU`):

# In[59]:


get_ipython().run_line_magic('pinfo2', 'fU._add_')


# We notice that the addition is performed on the symbolic expression
# with respect to the symbolic engine currently at work (SageMath/Pynac, SymPy,...), as returned by
# the method `expr()`.
# Let us recall that the user can change the symbolic engine at any time
# by means of the method `set_calculus_method()`, applied either to
# a chart or to an open subset (possibly `M` itself).
# Besides, we notice in the code above that the result of the symbolic addition
# is automatically simplified, by means of the method `_simplify`.
# It invokes a chain of simplification functions, which depends on the
# symbolic engine (See [here](https://github.com/sagemath/sage/blob/develop/src/sage/manifolds/utilities.py)
# for details).

# Let us now discuss the second case in the `__add__` method of
# `Element`, namely the case for which the parents of both operands are
# different. This case is treated via SageMath **coercion model**, which allows one to deal with additions like

# In[60]:


h1 = f + 2
h1.display()


# A priori, `f + 2` is not a well defined operation, since the integer $2$ does not
# belong to the algebra $C^\infty(M)$. However SageMath manages to treat it
# because $2$ can be coerced (i.e. automatically and unambigously converted) via `CM(2)`
# into a element of $C^\infty(M)$, namely the constant scalar field whose value is $2$:

# In[61]:


CM(2).display()


# This happens because there exists a coercion map from the parent of $2$, namely the ring of integers $\mathbb{Z}$
# (denoted `ZZ` in SageMath), to $C^\infty(M)$:

# In[62]:


2.parent()


# In[63]:


CM.has_coerce_map_from(ZZ)