.. toctree::
   :maxdepth: 2

.. _examples-topology:

################################################################################
Topology Analysis
################################################################################

The :func:`compute_topomesh_property <cellcomplex.property_topomesh.analysis.compute_topomesh_property>` method provides ways to store advanced topological features of a :class:`PropertyTopomesh` as ``wisp`` properties. This avoids unnecessary re-computations of topological relationships that are not diirectly encoded in the structure, and allows powerful global computations, making advantage of the NumPy encoding of the properties.

****************************
Basic topological properties
****************************

Let's consider a rather simple triangular mesh represented by a :class:`PropertyTopomesh`, here generated as a random triangulation inscribed within a circle. The only information that is stored in memory is the incidence graph displayed on the right, that is already somewhat complex.

.. code-block:: python
    :linenos:

    from cellcomplex.property_topomesh.example_topomesh import circle_voronoi_topomesh

    _,topomesh = circle_voronoi_topomesh(2,circle_size=8,return_triangulation=True)

.. plot::

    import numpy as np
    import matplotlib.pyplot as plt

    from cellcomplex.property_topomesh.example_topomesh import circle_voronoi_topomesh

    from cellcomplex.property_topomesh.utils.matplotlib_tools import mpl_draw_topomesh, mpl_draw_incidence_graph

    np.random.seed(108)
    _,topomesh = circle_voronoi_topomesh(2,circle_size=8,return_triangulation=True)

    figure = plt.figure(0)
    figure.clf()

    figure.add_subplot(1,2,1)
    figure.gca().axis('equal')
    mpl_draw_topomesh(topomesh,figure,2,color='g')
    mpl_draw_topomesh(topomesh,figure,1,color='b')
    mpl_draw_topomesh(topomesh,figure,0,color='m')
    figure.gca().axis('off')

    figure.add_subplot(1,2,2)
    figure.gca().axis('equal')
    mpl_draw_incidence_graph(topomesh,figure)
    figure.gca().axis('off')

    figure.set_size_inches(10,5)
    figure.tight_layout()

This kind of combinatorial structure can be very costly to manipulate in Python as traversals might involve several embedded :func:`for` loops. The philosophy behind :mod:`CellComplex <cellcomplex>` is to give the possibility to perform computations as globally as possible, by providing tools to investigate the structure through NumPy arrays. For instance, the :func:`compute_topomesh_property <cellcomplex.property_topomesh.analysis.compute_topomesh_property>` offers the possibility to access the essence of the topological structure (the :func:`borders` and :func:`regions` relationships) as array properties.


.. code-block:: python
    :linenos:

    from cellcomplex.property_topomesh.analysis import compute_topomesh_property

    compute_topomesh_property(topomesh,'borders',1)
    print("Edge borders: "+str(topomesh.wisp_property('borders',1).values()))

.. code-block:: bash

    Edge borders:   [[ 0  2]
                     [ 0  3]
                     [ 0  4]
                     [ 0  5]
                     [ 0  8]
                     [ 0  9]
                     [ 0 10]
                     [ 1  2]
                     [ 1  6]
                     [ 1 12]
                     [ 1 13]
                     [ 1 14]
                     [ 2  3]
                     [ 2  5]
                     [ 2 11]
                     [ 2 12]
                     [ 2 14]
                     [ 3  7]
                     [ 8  3]
                     [ 3 14]
                     [ 8  4]
                     [ 9  4]
                     [10  5]
                     [11  5]
                     [13  6]
                     [ 6 14]
                     [ 8  7]
                     [14  7]
                     [ 9 10]
                     [10 11]
                     [11 12]
                     [12 13]]

The main advantage of accessing data this way is the possibility to combine arrays of properties in order to make fast computations. For instance, to compute the length of all edges requires computing the norm of the difference of positions (``'barycenter'`` property) of the edge borders, which can be written in a **very condensed and efficient** way in NumPy.

.. code-block:: python
    :linenos:

    edge_borders = topomesh.wisp_property('borders',1).values()
    edge_border_positions = topomesh.wisp_property('barycenter',0).values(edge_borders)
    edge_lengths = np.linalg.norm(np.diff(edge_border_positions,axis=1)[:,0],axis=1)
    print("Edge lengths: "+str(edge_lengths))

.. code-block:: bash

    Edge lengths: [1.54168209 1.58307211 1.00382853 2.11521433 1.16696402 1.1110798
                   1.82197209 0.65995687 1.04620625 1.21589874 0.88074694 1.6630285
                   1.66661489 1.73242215 1.84463235 1.50941624 1.87907293 0.61652321
                   1.63064711 1.19783225 0.6710141  0.88910546 1.04256369 0.5089403
                   0.60474705 0.9297893  1.53073373 1.53073373 1.53073373 1.53073373
                   1.53073373 1.53073373]

A more usual way to perform the same computation would have been to iterate over the mesh edges using the :func:`wisps <cellcomplex.property_topomesh.topomesh.Topomesh.wisps>` method and then to use the  :func:`borders <cellcomplex.property_topomesh.topomesh.Topomesh.borders>` method to access each edge's vertices, to look for their value for the ``barycenter`` property, before computing the norm. Running a computational time evaluation on such simple example already reveals the interest to manipulate the topological structure as arrays when possible, and to store low-level information (that may not appear costly to compute at first, but is **systematically used**) as a ready-to-use array. Especially when the computations are to be repeated, as we see that the time gain is more than 3x if we don't need to perform the initial property storage.

.. code-block:: python
    :linenos:

    from time import time

    start_time = time()
    edge_border_positions = []
    for e in topomesh.wisps(1):
        edge_borders = list(topomesh.borders(1,e))
        edge_border_positions.append([topomesh.wisp_property('barycenter',0)[v] for v in edge_borders])
    edge_lengths = np.linalg.norm(np.diff(edge_border_positions,axis=1)[:,0],axis=1)
    print("Loop-based computation (everytime): "+str(np.round(1000.*(time()-start_time),2))+"ms")

    start_time = time()
    compute_topomesh_property(topomesh,'borders',1)
    edge_borders = topomesh.wisp_property('borders',1).values()
    edge_border_positions = topomesh.wisp_property('barycenter',0).values(edge_borders)
    edge_lengths = np.linalg.norm(np.diff(edge_border_positions,axis=1)[:,0],axis=1)
    print("NumPy computation (1st time): "+str(np.round(1000.*(time()-start_time),2))+"ms")

    start_time = time()
    edge_borders = topomesh.wisp_property('borders',1).values()
    edge_border_positions = topomesh.wisp_property('barycenter',0).values(edge_borders)
    edge_lengths = np.linalg.norm(np.diff(edge_border_positions,axis=1)[:,0],axis=1)
    print("NumPy computation (next time(s)): "+str(np.round(1000.*(time()-start_time),2))+"ms")

.. code-block:: bash

    Loop-based computation (everytime): 0.84ms
    NumPy computation (1st time): 0.59ms
    NumPy computation (next time(s)): 0.22ms

******************************
Complex topology as properties
******************************

As we saw, when performing computationally heavy tasks the recommended way to manipulate a :class:`PropertyTopomesh` is through the use of properties that reflect its topology. This is particularly true when *the same computation has to be performed iteratively on a structure with a fixed topology*, for instance whenever geometry is evolving and geometrical values need to be updated regularly. To be more efficient in this situation, the :func:`compute_topomesh_property <cellcomplex.property_topomesh.analysis.compute_topomesh_property>` gives access to other topological relationships that are not simply the boundary relationship between elements. Indirect relationships such as **neighbor** relationships, or incidence between more than one dimension (the **vertices** of a **face** typically) can be stored once as a property and a re-used to perform advanced computations.

.. code-block:: python
    :linenos:

    compute_topomesh_property(topomesh,'vertices',2)
    print("Face vertices: "+str(topomesh.wisp_property('vertices',2).values()))

.. code-block:: bash

    Face vertices:  [[ 0  2  9]
                     [ 0  4  9]
                     [ 5 11 12]
                     [ 0  4  5]
                     [ 6 12 13]
                     [ 0  5  6]
                     [ 5  6 12]
                     [ 2  7  8]
                     [ 2  8  9]
                     [ 1  7 14]
                     [ 1  2  7]
                     [ 0  1  2]
                     [ 0  1  6]
                     [ 1 13 14]
                     [ 1  6 13]
                     [ 4  9 10]
                     [ 3  4 10]
                     [ 3  5 11]
                     [ 3  4  5]]

In the following we will make use of these topological properties to perform a rather complex task: apply an iterative smoothing algorithm to the mesh geometry. The algorithm is a *Laplacian smoothing* algorithm, where the vertices are iteratively shifted towards the barycenter of their neighbors, but with a rather small coefficient. As the process will not affect the topology of the complex but only its geometry, we can compute the ``neighbors`` property on the vertices and store it in a NumPy array.

.. code-block:: python
    :linenos:

    compute_topomesh_property(topomesh,'neighbors',0)
    print(topomesh.wisp_property('neighbors',0))
    vertex_neighbors = topomesh.wisp_property('neighbors',0).values()

.. code-block:: bash

    {0:  [ 2   3   4   5   8   9  10]
     1:  [ 2   6  12  13  14]
     2:  [ 0   1   3   5  11  12  14] 
     3:  [ 0   2   7   8  14] 
     4:  [ 0   8   9]
     5:  [ 0   2  10  11] 
     6:  [ 1  13  14] 
     7:  [ 3   8  14] 
     8:  [ 0   3   4   7]
     9:  [ 0   4  10] 
     10: [ 0   5   9  11] 
     11: [ 2   5  10  12] 
     12: [ 1   2  11  13] 
     13: [ 1   6  12] 
     14: [ 1   2   3   6   7]}

Then to estimate the center of any vertex neighborhood, we will use the ``barycenter`` property on vertices that stores their 3d position, and apply the NumPy :func:`mean <numpy.mean>` function to average positions of the neighbor vertices.

.. code-block:: python
    :linenos:

    vertex_neighbor_positions =  topomesh.wisp_property('barycenter',0).values(vertex_neighbors)
    vertex_neighbor_center = np.array([np.mean(p,axis=0) for p in vertex_neighbor_positions])

The next step is to compute a new position for all the vertices by adding to the current position a vector in the towards the neighborhood center. The vector is multiplied by a small coefficient to ensure the stability of the process. Once computed, the new positions can be affected to the vertices as ``barycenter`` property through the :func:`update_wisp_property <cellcomplex.property_topomesh.PropertyTopomesh.update_wisp_property>` method.

.. code-block:: python
    :linenos:

    coef = 0.2
    vertex_positions = topomesh.wisp_property('barycenter',0).values()
    new_vertex_positions = vertex_positions + 0.33*(vertex_neighbor_center-vertex_positions)

    topomesh.update_wisp_property('barycenter',0,new_vertex_positions)


.. plot::

    import numpy as np
    import matplotlib.pyplot as plt

    from cellcomplex.property_topomesh.example_topomesh import circle_voronoi_topomesh
    from cellcomplex.property_topomesh.analysis import compute_topomesh_property

    from cellcomplex.property_topomesh.utils.matplotlib_tools import mpl_draw_topomesh, mpl_draw_incidence_graph

    from cellcomplex.property_topomesh.utils.delaunay_tools import delaunay_triangulation_topomesh

    np.random.seed(108)
    _,topomesh = circle_voronoi_topomesh(2,circle_size=8,return_triangulation=True)

    figure = plt.figure(0)
    figure.clf()

    figure.add_subplot(1,2,1)
    figure.gca().axis('equal')
    mpl_draw_topomesh(topomesh,figure,2,color='g')
    mpl_draw_topomesh(topomesh,figure,1,color='b')
    mpl_draw_topomesh(topomesh,figure,0,color='m')
    figure.gca().axis('off')

    compute_topomesh_property(topomesh,'neighbors',0)
    vertex_neighbors = topomesh.wisp_property('neighbors',0).values()

    coef = 0.2

    vertex_positions = topomesh.wisp_property('barycenter',0).values()
    vertex_neighbor_positions =  topomesh.wisp_property('barycenter',0).values(vertex_neighbors)
    vertex_neighbor_center = np.array([np.mean(p,axis=0) for p in vertex_neighbor_positions])
    new_vertex_positions = vertex_positions + coef*(vertex_neighbor_center-vertex_positions)

    for p,n_p in zip(vertex_positions,vertex_neighbor_center):
        figure.gca().plot([p[0],n_p[0]],[p[1],n_p[1]],color='gold',zorder=6,alpha=0.33)
    figure.gca().scatter(vertex_neighbor_center[:,0],vertex_neighbor_center[:,1],color='gold',s=10,alpha=0.33)
    figure.gca().scatter(new_vertex_positions[:,0],new_vertex_positions[:,1],color='gold',s=20,zorder=7)

    topomesh.update_wisp_property('barycenter',0,new_vertex_positions)

    figure.add_subplot(1,2,2)
    figure.gca().axis('equal')
    mpl_draw_topomesh(topomesh,figure,2,color='g')
    mpl_draw_topomesh(topomesh,figure,1,color='b')
    mpl_draw_topomesh(topomesh,figure,0,color='m',plot_ids=True)
    figure.gca().axis('off')

    figure.set_size_inches(10,5)
    figure.tight_layout()



The only thiing left to do is to iterate over this position update loop to perform an interative smoothing of the initial mesh.

.. code-block:: python
    :linenos:

    coef = 0.2

    for iteration in range(10):
        vertex_positions = topomesh.wisp_property('barycenter',0).values()
        vertex_neighbor_positions =  topomesh.wisp_property('barycenter',0).values(vertex_neighbors)
        vertex_neighbor_center = np.array([np.mean(p,axis=0) for p in vertex_neighbor_positions])
        new_vertex_positions = vertex_positions + coef*(vertex_neighbor_center-vertex_positions)

        topomesh.update_wisp_property('barycenter',0,new_vertex_positions)

.. plot::

    import numpy as np
    import matplotlib.pyplot as plt

    from cellcomplex.property_topomesh.example_topomesh import circle_voronoi_topomesh
    from cellcomplex.property_topomesh.analysis import compute_topomesh_property

    from cellcomplex.property_topomesh.utils.matplotlib_tools import mpl_draw_topomesh, mpl_draw_incidence_graph

    from cellcomplex.property_topomesh.utils.delaunay_tools import delaunay_triangulation_topomesh

    np.random.seed(108)
    _,topomesh = circle_voronoi_topomesh(2,circle_size=8,return_triangulation=True)

    figure = plt.figure(0)
    figure.clf()

    figure.add_subplot(1,2,1)
    figure.gca().axis('equal')
    mpl_draw_topomesh(topomesh,figure,2,color='g')
    mpl_draw_topomesh(topomesh,figure,1,color='b')
    mpl_draw_topomesh(topomesh,figure,0,color='m',plot_ids=True)
    figure.gca().axis('off')

    compute_topomesh_property(topomesh,'neighbors',0)
    vertex_neighbors = topomesh.wisp_property('neighbors',0).values()

    coef = 0.2

    for iteration in range(10):
        vertex_positions = topomesh.wisp_property('barycenter',0).values()
        vertex_neighbor_positions =  topomesh.wisp_property('barycenter',0).values(vertex_neighbors)
        vertex_neighbor_center = np.array([np.mean(p,axis=0) for p in vertex_neighbor_positions])
        new_vertex_positions = vertex_positions + coef*(vertex_neighbor_center-vertex_positions)

        topomesh.update_wisp_property('barycenter',0,new_vertex_positions)

    figure.add_subplot(1,2,2)
    figure.gca().axis('equal')
    mpl_draw_topomesh(topomesh,figure,2,color='g')
    mpl_draw_topomesh(topomesh,figure,1,color='b')
    mpl_draw_topomesh(topomesh,figure,0,color='m',plot_ids=True)
    figure.gca().axis('off')

    figure.set_size_inches(10,5)
    figure.tight_layout()

The result of the smoothing is a global regularization of the shapes of the mesh triangles that less elongated than in the initial configuration, at the expense of the overall shape of the mesh. The use of topological properties in this case allows to write the process at a very high level and avoids going into multiple imbricated loops to obtain the same result. Most of the code in :mod:`CellComplex` makes use of this powerful tool for manipulating the structure, notably all the geometrical analysis tools described in :ref:`examples-analysis`.

*********************************
Computable topological properties
*********************************

A list of **topological** properties that can be passed as `property_name` argument to the function :func:`compute_topomesh_property <cellcomplex.property_topomesh.analysis.compute_topomesh_property>`, along with the `wisp` dimension for which they can be computed:

+---------------------+---+---+---+---+----------+
| property_name       | 0 | 1 | 2 | 3 |   type   |
+=====================+===+===+===+===+==========+
|``borders``          |   | x | x | x |   list   |
+---------------------+---+---+---+---+----------+
|``regions``          | x | x | x |   |   list   |
+---------------------+---+---+---+---+----------+
|``vertices``         | x | x | x | x |   list   |
+---------------------+---+---+---+---+----------+
|``triangles``        | x | x | x | x |   list   |
+---------------------+---+---+---+---+----------+
|``cells``            | x | x | x | x |   list   |
+---------------------+---+---+---+---+----------+
|``neighbors``        | x | x | x | x |   list   |
+---------------------+---+---+---+---+----------+