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()

Directional Derivative

Vector fields act on functions via X: Function $\to$ Function, defined by $X(f) = df(X)$. At a point $p$ this evaluates the derivative $df(X)(p)=Df(x)[X(x)]\in\mathbb{R}$. In Diffusion Geometry, functions are callable objects for vector fields.

"""
Calculate the directional derivative of the funciton f(x,y)=x along the vector fields X=∇_x and Y=∇_y.
"""

L=1
num_points=10

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

dg = DiffusionGeometry.from_point_cloud(data_meshgrid)


# Define f(x,y)=x 
f_pointwise = data_meshgrid[:, 0]
f = Function.from_pointwise_basis(f_pointwise, dg)

# Define a vector field pointing along the x-axis
X_pointwise = np.zeros_like(data_meshgrid)
X_pointwise[:,0] = 10
X = VectorField.from_pointwise_basis(X_pointwise, dg)

# Define a vector field pointing along the y-axis
Y_pointwise = np.zeros_like(data_meshgrid)
Y_pointwise[:,1] = 10
Y = VectorField.from_pointwise_basis(Y_pointwise, dg)

# Calculate the directional derivative simply by calling f with the vector fields
X_f = X(f)
Y_f = Y(f)





fig = make_subplots(
    rows=2, cols=3, 
    subplot_titles=["f(x,y)=x", "X=(x,1)", "Y=(0,1)", "", "X(f)≈1 ", "Y(f)≈0"],
    horizontal_spacing=0.05
)

scale=0.02
line_width=1

fig.add_trace(plot_scatter_2d(data_meshgrid, f.to_ambient()).data[0])

fig.add_trace(plot_quiver_2d(data_meshgrid, X.to_ambient(), scale=scale, line_width=line_width).data[0], col=2, row=1)
fig.add_trace(plot_scatter_2d(data_meshgrid, X_f.to_ambient()).data[0], row=2,col=2)

fig.add_trace(plot_quiver_2d(data_meshgrid, Y.to_ambient(), scale=scale, line_width=line_width).data[0], col=3, row=1)
fig.add_trace(plot_scatter_2d(data_meshgrid, Y_f.to_ambient()).data[0], row=2,col=3)

fig.update_yaxes(scaleanchor="x", scaleratio=1)
clean_fig(fig)

fig.update_layout(
    height=600, 
    width=1000,)
fig.update_layout(margin=dict(l=20, r=20, t=40, b=20))

fig.show()

Lie Bracket

The Lie bracket measures the failure of commutativity of the flows of two vector fields. It is a map $[\cdot\ , \ \cdot]$: VectorField \times VectorField $\to$ VectorField, and it is defined via $$[X,Y](f) = X(Y(f)) - Y(X(f)). $$ Computationally, the Lie bracket can be called with the method lie_bracket on the DiffusionGeometry class.

"""
Example of Lie bracket calculation for X=x ∇_y and Y=y∇_x. The theoretical Lie bracket is then [X,Y]=x\nabla_x - y\nabla_y
"""

def sample_perturbed_grid(range_x, range_y, num_x, num_y, noise_strength=0.3):
    x_lin = np.linspace(range_x[0], range_x[1], num_x)
    y_lin = np.linspace(range_y[0], range_y[1], num_y)
    x_grid, y_grid = np.meshgrid(x_lin, y_lin)
    data = np.column_stack([x_grid.ravel(), y_grid.ravel()])
    data += noise_strength * np.random.randn(*data.shape)
    return data


def create_figure(fig, row):
    # Create grid data
    num_x = 20
    num_y = 20 
    lim = -4
    data = sample_perturbed_grid(range_x=(-lim,lim), range_y=(-lim,lim), num_x=num_x, num_y=num_y, noise_strength=0.0)

    # Create DiffusionGeometry instance
    # One can play with n_function_basis to increase/decrease precision on the Lie bracket
    dg = DiffusionGeometry.from_point_cloud(data, n_function_basis=100)

    # ----------------------------------------------
    # Create the two vector fields X and  Y
    # ----------------------------------------------
    X_pointwise = np.zeros((dg.n, 2))
    X_pointwise[:, 1] = data[:, 0]
    X = VectorField.from_pointwise_basis(X_pointwise, dg)

    Y_pointwise = np.zeros((dg.n, 2))
    Y_pointwise[:, 0] = data[:, 1]
    Y = VectorField.from_pointwise_basis(Y_pointwise, dg)

    # ----------------------------------------------
    # Compute the Lie bracket
    # ----------------------------------------------
    Lie_XY = dg.lie_bracket(X, Y)

    # ----------------------------------------------
    # Create theoretical vector field
    # ----------------------------------------------
    theoretical_XY_coeffs = np.zeros((dg.n, 2))
    theoretical_XY_coeffs[:, 0] = data[:, 0]
    theoretical_XY_coeffs[:, 1] = -data[:, 1]
    theoretical_XY = VectorField.from_pointwise_basis(theoretical_XY_coeffs, dg)

    scale = 0.12
    line_width = 1

    fig1 = plot_quiver_2d(data, X.to_ambient(), scale=scale, line_width=line_width)
    fig2 = plot_quiver_2d(data, Y.to_ambient(), scale=scale, line_width=line_width)
    fig3 = plot_quiver_2d(data, Lie_XY.to_ambient(), scale=scale, line_width=line_width)
    fig4 = plot_quiver_2d(data, (theoretical_XY-Lie_XY).to_ambient(), scale=scale, line_width=line_width)



    fig.add_traces(fig1.data, rows=row, cols=1)
    fig.add_traces(fig2.data, rows=row, cols=2)
    fig.add_traces(fig3.data, rows=row, cols=3)
    fig.add_traces(fig4.data, rows=row, cols=4)


fig = make_subplots(
    rows=1,
    cols=4, 
    subplot_titles=("X=x∇_y","X=y∇_x","Diffusion Geometry Lie Bracket","Error to Theoretical Value"),
    horizontal_spacing=0.05,
    vertical_spacing=0.01,
    shared_xaxes=True,
    shared_yaxes=True
)

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

Hessian

On a manifold, the Hessian measures the curvature of a function, taking the intrinsic curvature of the manifold into account. Given a function $f$ and two vector fields $X$ and $Y$ providing the direction in which the curvature is measured, we have that $H(f)(X,Y)\in\mathbb{R}$. In Diffusion Geometry, we access it via the CdC: $$ H(f)(\nabla h_1, \nabla h_2) =\frac{1}{2}\Big(\Gamma(h_1, \Gamma(h_2, f))+ \Gamma(h_2, \Gamma(h_1, f))- \Gamma(f, \Gamma(h_1, h_2))\Big).$$

data1 = np.random.uniform(-1, 1, (2000, 2))

# sample data1 from a grid:
x, y = np.meshgrid(np.linspace(-1, 1, 40), np.linspace(-1, 1, 40))
data1 = np.vstack([x.ravel(), y.ravel()]).T
data1 += 0.015 * np.random.randn(*data1.shape)
dg1 = DiffusionGeometry.from_point_cloud(data1)

x, y = data1[:,0], data1[:,1]

# Create function with interesting curvature along x and y
f = dg1.function(np.sin(5*x) + np.cos(5*y))

# ---------------------------------------------------------
# Compute the Hessian in all directions
# ---------------------------------------------------------

# Compute the tensor first
Hf = f.hessian()

# Apply it in all directions
dx = dg1.function(x).grad()
dy = dg1.function(y).grad()
Hxx = Hf(dx, dx)
Hxy = Hf(dx, dy)
Hyy = Hf(dy, dy)


# ---------------------------------------------------------
# Plot the Hessian in all directions
# ---------------------------------------------------------


fig = make_subplots(
    rows=1,
    cols=4, 
    subplot_titles=("f=sin(5x)+cos(5y)","H(f)(∇x,∇x)","H(f)(∇y,∇y)", "H(f)(∇x,∇y)"),
    horizontal_spacing=0.05,
    vertical_spacing=0.01,
    shared_xaxes=True,
    shared_yaxes=True
)

scatter_size = 4
amax = max(np.max(np.abs(Hxx)), np.max(np.abs(Hyy)), np.max(np.abs(Hxy)))

f1 = plot_scatter_2d(data1, f.to_ambient())
fig.add_trace(f1.data[0])

f2 = plot_scatter_2d(data1, Hxx, size=scatter_size, amax=amax)
fig.add_trace(f2.data[0], col=2, row=1)

f3 = plot_scatter_2d(data1, Hyy, size=scatter_size, amax=amax)
fig.add_trace(f3.data[0], col=3, row=1)

f4 = plot_scatter_2d(data1, Hxy, size=scatter_size, amax=amax)
fig.add_trace(f4.data[0], col=4, row=1)



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

Levi-Civita

Connections are the generalisation of directional derivatives from functions to vector fields on a Riemannian manifold. For a given manifold, there is an infinite choice of connections. The Levi-Civita connection is a special one, defined by two conditions:

• Metric Compatibility: $\nabla g=0$ (lengths/angles are preserved under parallel transport).
• Torsion-free: $\nabla_X(​Y)−\nabla_Y(​X)=[X,Y]$ (no "twisting" of the space). In Diffusion Geometry, we get to it via the Hessian $$\nabla X := \sum_{i,j} ( d\phi_i \otimes d x_j + \phi_i H(x_j)).$$ The Levi-Civita connection is precomputed by the DG framework, and can be accesses via the levi_civita method. One important caveat is that, given a VectorField X, we return the operator $\nabla X$ that acts on another vector field $Y$ via $$\nabla X(Y) = \nabla_YX. $$

"""
Example using the Levi-Civita connection.
"""

def cov_der(dg, X, Y):
    return dg.levi_civita(Y)(X)


def sample_perturbed_grid(range_x, range_y, num_x, num_y, noise_strength=0.3):
    x_lin = np.linspace(range_x[0], range_x[1], num_x)
    y_lin = np.linspace(range_y[0], range_y[1], num_y)
    x_grid, y_grid = np.meshgrid(x_lin, y_lin)
    data = np.column_stack([x_grid.ravel(), y_grid.ravel()])
    data += noise_strength * np.random.randn(*data.shape)
    return data



def create_figure(fig, row):
    n_points = 20 
    lim = 4

    data = sample_perturbed_grid(range_x=(-lim, lim), range_y=(-lim, lim), num_x=n_points, num_y=n_points, noise_strength=0.1)
    dg = DiffusionGeometry.from_point_cloud(data, n_function_basis=399)

    X_pointwise = np.zeros((dg.n, 2))
    X_pointwise[:, 0] = 3
    X = VectorField.from_pointwise_basis(X_pointwise,dg) 

    Y_pointwise = np.zeros((dg.n, 2))
    Y_pointwise[:, 1] = data[:, 0] 
    Y = VectorField.from_pointwise_basis( Y_pointwise,dg)
    nabla_X_Y = cov_der(dg, X, Y)

    scale = 0.12
    line_width = 1
    arrow_scale = 0.3


    fig1 = plot_quiver_2d(data, X.to_ambient(), scale=scale, arrow_scale=arrow_scale, line_width=line_width)
    fig2 = plot_quiver_2d(data, Y.to_ambient(), scale=scale, arrow_scale=arrow_scale, line_width=line_width)
    fig3 = plot_quiver_2d(data, nabla_X_Y.to_ambient(), scale=scale, arrow_scale=arrow_scale, line_width=line_width)


    fig.add_traces(list(fig1.data), rows=row, cols=1)
    fig.add_traces(list(fig2.data), rows=row, cols=2)
    fig.add_traces(list(fig3.data), rows=row, cols=3)


fig = make_subplots(
    rows=1,
    cols=3, 
    subplot_titles=("X=∇x","Y=x∇y","∇_X(Y)"),
    horizontal_spacing=0.05,
    vertical_spacing=0.01,
    shared_xaxes=True,
    shared_yaxes=True
)

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

Hodge Decomposition

We can perform the Hodge decomposition of a $k$-form $\alpha \in \Omega^k(M)$ into $$\alpha = d_{k-1}\beta + \gamma + \partial_{k+1}\delta,$$ where $\beta \in \Omega^{k-1}(M)$, $\gamma \in \Omega^k(M)$, and $\delta \in \Omega^{k+1}(M)$, with $d_k \gamma = 0$ and $\partial_k \gamma = 0$.

We call $d_{k-1}\beta$ the exact part of $\alpha$, and $\gamma$ the harmonic part, and $\partial_{k+1}\delta$ the coexact part. We will refer to $\beta$ as the exact potential of $\alpha$, and $\delta$ as the coexact potential. The Hodge decomposition is unique, and the three parts are orthogonal.

The harmonic part $\gamma$ lies in the kernel of the Hodge Laplacian $\Delta_k$. Intuitively, the harmonic and coexact parts of $\alpha$ both capture its rotational behaviour: the coexact part corresponds to vortices whose centre is in the space, and the harmonic part corresponds to global circular flows around holes in the space.

"""
Showcase of how the Hodge decomposition can be used to isolate divergence and rotational parts of a vector field.
"""


def sample_perturbed_grid(range_x, range_y, num_x, num_y, noise_strength=0.3):
    x_lin = np.linspace(range_x[0], range_x[1], num_x)
    y_lin = np.linspace(range_y[0], range_y[1], num_y)
    x_grid, y_grid = np.meshgrid(x_lin, y_lin)
    data = np.column_stack([x_grid.ravel(), y_grid.ravel()])
    data += noise_strength * np.random.randn(*data.shape)
    return data


n_points = 25
lim = 4

data = sample_perturbed_grid(range_x=(-lim, lim), range_y=(-lim, lim), num_x=n_points, num_y=n_points, noise_strength=0.1)
dg = DiffusionGeometry.from_point_cloud(data, n_function_basis=399)

div_vf_pointwise = np.zeros((dg.n, 2))
div_vf_pointwise = data[:]
div_vf = dg.vector_field(div_vf_pointwise)


rot_vf_pointwise = np.zeros((dg.n, 2))
rot_vf_pointwise[:,0] = -data[:,1]
rot_vf_pointwise[:,1] =  data[:,0]
rot_vf = dg.vector_field(rot_vf_pointwise)

exact, coexact, harm = (rot_vf+div_vf).flat().hodge_decomposition()

fig = make_subplots(
    rows=2,
    cols=3, 
    subplot_titles=(
    "Input: Divergence (Source)", 
    "Input: Rotational (Vortex)", 
    "Combined Field (Sum)",
    "Extracted: Exact Part", 
    "Extracted: Coexact Part", 
    "Extracted: Harmonic Part"
),
    horizontal_spacing=0.05,
    vertical_spacing=0.05,
    shared_xaxes=True,
    shared_yaxes=True
)

scale=0.1
line_width = 0.5

fig1 = plot_quiver_2d(data, div_vf.to_ambient(), scale, line_width=line_width).data[0]
fig.add_trace(fig1, row=1, col=1)

fig2 = plot_quiver_2d(data, rot_vf.to_ambient(), scale, line_width=line_width).data[0]
fig.add_trace(fig2, row=1, col=2)

fig3 = plot_quiver_2d(data, (rot_vf+div_vf).to_ambient(), scale, line_width=line_width).data[0]
fig.add_trace(fig3, row=1, col=3)

fig4 = plot_quiver_2d(data, exact.d().sharp().to_ambient(), scale, line_width=line_width).data[0]
fig.add_trace(fig4, row=2, col=1)

fig5 = plot_quiver_2d(data, coexact.codifferential().sharp().to_ambient(), scale, line_width=line_width).data[0]
fig.add_trace(fig5, row=2, col=2)

fig6 = plot_quiver_2d(data, harm.to_ambient(), scale, line_width=line_width).data[0]
fig.add_trace(fig6, row=2, col=3)


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

fig.update_yaxes(scaleanchor="x", scaleratio=1)
fig.update_xaxes(scaleanchor="y", scaleratio=1)

fig.show()