5. Differential Operators

In this chapter, we showcase how diffusion geometry computes differential operators, such as gradients, the Lie bracket, the Hessian, and the Levi-Civita connection.

The Gradient

The gradient is a map $\nabla$: Function $\to$ VectorField. At each point, it points in the direction of fastest increase of a function.

import sys
import os

sys.path.append(os.path.abspath('..'))  # Add parent directory to path
from diffusion_geometry import DiffusionGeometry
from diffusion_geometry.visualisation import *

import numpy as np

from plotly.subplots import make_subplots
from figures.generate_data import gen_3d_data, load_image_point_cloud


# Generate Swiss Roll data
data_swiss_roll, _ = gen_3d_data(kind = "swiss_roll", n = 2000, noise = 0.05)
dg = DiffusionGeometry.from_point_cloud(data_swiss_roll, n_function_basis = 100)

# Pick an eigenfunction as an example
f = dg.function_space.zeros()
f.coeffs[8] = 1

camera = dict(eye=dict(x=-0.4, y=1.85, z=0.5))
plot_scatter_3d(data_swiss_roll, color = f.to_ambient(), camera=camera).show()

# Plot the corresponding gradient vector field
# Computing the gradient is as simple as calling the grad() method on the function
gradient_f = f.grad()

plot_quiver_3d(
    data_swiss_roll, 
    quiver=gradient_f.to_ambient(), 
    scale=4, 
    line_width=1,
    camera=camera
).show()

Divergence

The divergence is defined as the negative adjoint $-\nabla^*$:VectorField$\to$Function. It measures the local expansion and contraction of a space under the flow of a vector field. One can compute it by calling the div() method on vector fields.

# Generate Torus Roll data
data_torus, _ = gen_3d_data(kind = "torus", n = 2000, noise = 0.05)
dg = DiffusionGeometry.from_point_cloud(data_torus, n_function_basis = 100)

# Pick an example vector field
X = dg.vector_field_space.zeros()
X.coeffs[3] = 1

camera = dict(eye=dict(x=2, y=0.4, z=0.5))

# Plot the vector field
plot_quiver_3d(data_torus, quiver = X.to_ambient(), camera=camera, scale=0.1, arrow_scale=0.5, line_width=1.5).show()

# Calculate the divergence
div_X = X.div()

# Plot the divergence as a function on the torus
plot_scatter_3d(data_torus, color = div_X.to_ambient(), camera=camera).show()

Exterior derivative

The exterior derivative $d$ is an operator that takes a k-Form to a (k+1)-Form and generalizes familiar notions like gradients, curls, and divergences in a coordinate-free way.

Locally, and for a smooth function $f$ (a 0-form), the exterior derivative $$df = \sum_i \frac{\partial f}{\partial x^i}\, dx^i,$$ is just the differential of $f$. Moreover, we have that $$df = (\nabla f)^{\flat}. $$ In general, it is defined as the unique operator that is linear, satisfies $$ d(\alpha \wedge \beta) = d\alpha \wedge \beta + (-1)^k \alpha \wedge d\beta,$$ and is nilpotent: $$d^2 = 0.$$

Geometrically, $d$ measures how a differential form fails to be locally exact, and it underlies Stokes’ theorem in all dimensions. The following is a minimal example involving the exterior derivative.

# Example of using the differential d operator on functions and k-forms

random_data = np.random.rand(100, 3)
dg = DiffusionGeometry.from_point_cloud(random_data)

# Define a function
f = dg.function_space.zeros() 
df = f.d()

print(f"The differential of a function is a 1-form, type(df): {type(df)}") 
print(f"df has degree {df.degree}.\n")

# Define a k-form 
k=2
alpha = dg.form_space(k).zeros()

d_alpha = alpha.d()
print(f"The degree of a is {alpha.degree}, and the degree of da is {d_alpha.degree}.")

The differential of a function is a 1-form, type(df): <class 'diffusion_geometry.tensors.forms.form.Form'>
df has degree 1.

The degree of a is 2, and the degree of da is 3.

""" 
This script computes the curl of a rotational vector field rotational vector field representing two counter-rotating vortices.
We compute the curl via the formula curl(X) = (*d X^♭)^♯. The resulting curl vector field then identifies the axis of rotation of the vortices.
"""

from diffusion_geometry.tensors.vector_fields.vector_field import VectorField


# Data = 3D meshgrid points
x = np.linspace(-4, 4, 10)
y = np.linspace(-4, 4, 10)
z = np.linspace(-4, 4, 10)
X, Y, Z = np.meshgrid(x, y, z)
data_meshgrid = np.vstack([X.ravel(), Y.ravel(), Z.ravel()]).T
dg = DiffusionGeometry.from_point_cloud(data_meshgrid)

# ------------------------------------------------------------------
# Define the rotational vector field (two counter-rotating vortices)
# ------------------------------------------------------------------

# Constants for stability
eps = 1e-6 

def get_vortex_field(x, y, z, center_x):
    """Computes a 2D rotational field centered at (center_x, 0)"""
    dx = x - center_x
    dy = y
    r_sq = dx**2 + dy**2 + eps
    
    return np.stack([dy / r_sq, -dx / r_sq, np.zeros_like(z)], axis=1)

Xp = X.ravel()
Yp = Y.ravel()
Zp = Z.ravel()

v1 = get_vortex_field(Xp, Yp, Zp, center_x=1)
v2 = get_vortex_field(Xp, Yp, Zp, center_x=-1)

# The total field is the superposition of both (counter-rotating)
X_pointwise_coeff = v1 - v2  # Note: Use minus if you want them spinning in opposite directions
X = VectorField.from_pointwise_basis(X_pointwise_coeff, dg)

# ------------------------------------------------------------------
# Calculate the curl 
# ------------------------------------------------------------------

# Get corresponding 1-form 
X_flat = X.flat()
d_X_flat = X_flat.d()

# The hodge_star_2_form function takes a 2-form in ambient coordinates and returns the Hodge star of that form as a vector field in ambient coordinates
curl_X = hodge_star_2_form(d_X_flat.to_ambient())

# ------------------------------------------------------------------
# plot the original vector field and its curl
# ------------------------------------------------------------------

camera = dict(eye=dict(x=1.5, y=1.5, z=1.5))
fig = make_subplots(rows=1, cols=2, specs=[[{'type': 'scene'}, {'type': 'scene'}]],
                    subplot_titles=("Original Vector Field", "Curl of the Vector Field"))
fig.add_traces(list(plot_quiver_3d(data_meshgrid, quiver=X.to_ambient(), scale=2, line_width=1, arrow_scale=0.5).data), rows=1, cols=1)
fig.add_traces(list(plot_quiver_3d(data_meshgrid, quiver=curl_X, scale=1, arrow_scale=0.4, line_width=1).data), rows=1, cols=2)
fig.update_layout(height=600, width=1200, scene_camera=camera)

Another way to visualise the vorticity is directly via the 2-form. The coloured signed areas indicate areas of high vorticity. This clearly identifies the two vortices from the above example.

plot_2form_3d(data_meshgrid, d_X_flat.to_ambient(), radius=0.8).show()

Codifferential

The codifferential $\partial^k:$ k-Form $\to$ (k-1)-Form is the dual to the exterior derivative w.r.t. the metric defined on forms. We illustrate the action of the codifferential on 2-forms, where it defines a 1-form whose dual vector field is a rotational flow around the boundaries of the support of the 2-form. The direction of rotation reflects the orientation of the 2-form.

"""
Exterior Calculus on 2D Point Clouds: From Scalar Fields to Rotational Flows

1. Geometry: Defines a 2D meshgrid and initializes the DiffusionGeometry framework.
2. Modulation: Creates a scalar function f = cos(0.5x)sin(0.5y) to serve as a density map.
3. 2-Form (α): Scales the volume form by f, representing localized patches of flux.
4. Codifferential (δα): Maps the 2-form to a 1-form. Geometrically, this extracts 
   the "boundary flow" or circulation induced by the area density α.
5. Result: The output is a divergence-free 1-form that loops around the extrema of f.
"""

# Create a 2D meshgrid 



from diffusion_geometry.tensors.functions.function import Function


L = 15
num_points = 40
x = np.linspace(-L, L, num_points)
y = np.linspace(-L, L, num_points)
X, Y = np.meshgrid(x, y)
data_meshgrid = np.vstack([X.ravel(), Y.ravel()]).T

n_function_basis = 200
n_coefficients = 200

dg = DiffusionGeometry.from_point_cloud(data_meshgrid, n_function_basis=n_function_basis, n_coefficients=n_coefficients)

# create the function cos(x) * sin(y)
freq = 0.5
f_pointwise = np.cos(freq * data_meshgrid[:, 0]) * np.sin(freq * data_meshgrid[:, 1])
f = Function.from_pointwise_basis(f_pointwise, dg)

# Create constant 2-form
alpha = dg.form_space(2).zeros()
alpha.coeffs[0] = 1
alpha = f * alpha  # Scale by the function f

# Calculate the Codifferential
codiff_alpha = alpha.codifferential()

# 1. Create the Figure with 3 columns
fig = make_subplots(
    rows=1, cols=3, 
    subplot_titles=("Scalar Field f(x,y)=cos(0.5*x)*sin(0.5*y)", "Constant 2-form scaled by f", "Codifferential of 2-form (Induced Flow)"),
    horizontal_spacing=0.05
)

# 2. Plot the Function (Scalar Field)
# Using the trace from plot_scatter_2d
f_trace = plot_scatter_2d(data_meshgrid, color=f.to_ambient()).data[0]
fig.add_trace(f_trace, row=1, col=1)

# 3. Plot the 2-Form
# Using the trace from plot_2form_2d
alpha_trace = plot_2form_2d(data_meshgrid, alpha.to_ambient(), radius=1)
for trace in alpha_trace.data:
    fig.add_trace(trace, row=1, col=2)

# 4. Plot the Codifferential (1-Form / Quiver)
# Using the trace from plot_quiver_2d
codiff_trace = plot_quiver_2d(data_meshgrid, quiver=codiff_alpha.to_ambient(), scale=2, line_width=0.5).data[0]
fig.add_trace(codiff_trace, row=1, col=3)

# 5. Update layout for a clean look
fig.update_layout(
    height=500, 
    width=1200, 
    showlegend=False,
    template="plotly_white"
)

# Optional: Ensure all axes are equal for geometric accuracy
fig.update_yaxes(scaleanchor="x", scaleratio=1)

clean_fig(fig)
fig.update_layout(margin=dict(l=20, r=20, t=40, b=20))

fig.show()