HeinrichHartmann.com

# @CX-REMOVE
!pip install Pillow

Requirement already satisfied: Pillow in /opt/conda/lib/python3.11/site-packages (10.1.0)

Convolution Algebras and Discrete Geometry - a Study

We start our exploration of numerical geometry with a function on a plane. We represent this function by an array: F[x,y] where x,y are points on an equidistant grid with step size h.

# Example Function
#@CX-TOGGLE
import numpy as np
import matplotlib.pyplot as plt

def f(x, y):
    return (np.sin(x) + np.cos(y) + 0.3 + 3*x - 0.2*y) / (1 + x**2 + y**2 + 20*(0.2*x*y + 0.03*x + 0.02*y)**2)
    # return np.arctan(1*x - 2*y + 0.1 * (x**4 + y**4) )

# Generate a grid of (x, y) coordinates
N = 100
W = 6
x = np.linspace(-W, W, 2*N)
y = np.linspace(-W, W, 2*N)
X, Y = np.meshgrid(x, y)  # This creates a 2D grid of x, y values
h = x[1]-x[0]

# Compute f(x,y) for each coordinate in the grid
F = f(X, Y)

# Plot the heatmap
def PLOT(F):
    plt.figure(figsize=(5,5))
    plt.imshow(F, extent=[-W, W, -W, W], origin='lower', cmap='RdYlBu_r')
    plt.colorbar(shrink=0.8)
    plt.xlabel('x')
    plt.ylabel('y')

PLOT(F)
plt.title('Heatmap of F');

Directional Derivatives

We compute the "balanced" directional derivatives by convolution on the array F:

DxF[x,y] = (F[x+1] - F[x-1])/2/h = [1, 0, -1]/2/h @ F

Here @ denotes the discrete convolution operation. The formula is obtained as the arithmetic mean of the forward differential (F[x+1]-F[x])/h and the backward differential (F[x]-F[x-1])/h.

# Kernels and Derivatives
#@CX-TOGGLE

import numpy as np
from scipy.signal import convolve2d

class Kernel(np.ndarray):
    def __new__(cls, input_array):
        obj = np.asarray(input_array).view(cls)
        return obj

    def __matmul__(self, other):
        """Overloads @ operator for convolution."""
        if isinstance(other, Kernel):
            return Kernel(convolve2d(self, other, mode='full', boundary='fill'))
        elif isinstance(other, np.ndarray):
            return convolve2d(other, self, mode='same', boundary='symm')
        else:
            raise ValueError("Unsupported type for convolution.")
    
    def __add__(self, other):
        # Get the shape of the larger kernel
        shape = (max(self.shape[0], other.shape[0]), max(self.shape[1], other.shape[1]))
        result = np.zeros(shape)
        # x: [ ----------------------------- ]
        #    [ pad_0 ][    self     ][ pad_1 ]
        padx = (shape[0] - self.shape[0]) // 2
        pady = (shape[1] - self.shape[1]) // 2
        result[ padx : padx + self.shape[0], pady : pady + self.shape[1] ] += self
        padx = (shape[0] - other.shape[0]) // 2
        pady = (shape[1] - other.shape[1]) // 2
        result[ padx : padx + other.shape[0], pady : pady + other.shape[1] ] += other
        return Kernel(result)

from scipy.signal import convolve2d
import numpy as np
import matplotlib.pyplot as plt

# Compute the partial derivatives using convolution
Dx = Kernel(np.array([[1, 0, -1]]))/2/h
Dy = Dx.T
DxF = Dx @ F
DyF = Dy @ F

fig, axes = plt.subplots(1, 3, figsize=(18, 6))

# Plot F
im1 = axes[0].imshow(F, extent=[-W, W, -W, W], origin='lower', cmap='RdYlBu_r')
axes[0].set_title('F')
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
fig.colorbar(im1, ax=axes[0])

# Plot DxF
im2 = axes[1].imshow(DxF, extent=[-W, W, -W, W], origin='lower', cmap='RdYlBu_r')
axes[1].set_title('Dx F')
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
fig.colorbar(im2, ax=axes[1])

# Plot DyF
im3 = axes[2].imshow(DyF, extent=[-W, W, -W, W], origin='lower', cmap='RdYlBu_r')
axes[2].set_title('Dy F')
axes[2].set_xlabel('x')
axes[2].set_ylabel('y')
fig.colorbar(im3, ax=axes[2])

plt.tight_layout()
plt.show()

Discrete Laplace Operator

The Laplace Operator is defined as , can be expressed via convolution:

LF = Dx @ Dx @ F + Dy @ Dy @ F = (Dx @ Dx + Dy @ Dy) @ F

The convolution kernel L = (Dx @ Dx + Dy @ Dy) looks as follows:

Dxx = Dx @ Dx
Dyy = Dy @ Dy
L0 = (Dxx + Dyy)
print(L0*h*h*4)

[[ 0.  0.  1.  0.  0.]
 [ 0.  0.  0.  0.  0.]
 [ 1.  0. -4.  0.  1.]
 [ 0.  0.  0.  0.  0.]
 [ 0.  0.  1.  0.  0.]]

We can "compress" this into the standard Laplace Kernel:

L = Kernel(np.array([[0,1,0],[1,-4,1],[0,1,0]]))/h/h/4
print(L*h*h*4)

[[ 0.  1.  0.]
 [ 1. -4.  1.]
 [ 0.  1.  0.]]

Note that, convolution with this kernel, compares each F[x,y] to the average over the adjacent values.

# Plot Laplace of F
# @CX_TOGGLE
def PLOT2(F, G, imshow_args={}):
    # Create subplots
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

    # Plot the original matrix F
    ax1.imshow(F, extent=[-W, W, -W, W], origin='lower', cmap='RdYlBu_r', **imshow_args)
    ax1.set_xlabel('x')
    ax1.set_ylabel('y')

    # Plot the result of applying Laplacian to F
    ax2.imshow(G, extent=[-W, W, -W, W], origin='lower', cmap='RdYlBu_r', **imshow_args)
    ax2.set_xlabel('x')
    ax2.set_ylabel('y')
    
    # Adjust spacing between subplots
    plt.tight_layout()

    return (ax1, ax2)

ax1, ax2 = PLOT2(F, L@F, {"vmin":-2, "vmax":3})
ax1.set_title('Original F')
ax2.set_title('Laplacian of F');

Heat Flow

The heat equation is given by $\partial F / \partial t = \Delta F$, heuristically, we changethe value of F along in the direction of the Laplacian at that point.

It's straight forward to simulate the heat flow numerically by subsituting:

F = F + eps * L @ F

in an interative way.

When trying this out, naively, we run into numeric instabilities with checker patterns "blowing up". In order to get a numerically stable version, we need to replace the Laplacian with another version, that takes into account diagonal entries. The heuristic of averaging the neighbours remains the same:

# Stable 2D Laplace Kernel
LS = Kernel(np.array([[0.25,.5,.25],[.5,-3,.5],[.25,.5,.25]]))/h/h/4
print(LS*h*h*4)

[[ 0.25  0.5   0.25]
 [ 0.5  -3.    0.5 ]
 [ 0.25  0.5   0.25]]

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy.signal import convolve2d

G = F.copy()

# Create a figure and axis for the animation
fig, ax = plt.subplots(figsize=(5, 5))
cax = ax.imshow(G, extent=[-W, W, -W, W], origin='lower', cmap='RdYlBu_r', vmin=-2, vmax=2)
plt.colorbar(cax, shrink=0.8)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Heat Equation Simulation')

# Define a function to update the animation
def update(frame):
    global G
    for _ in range(30):
        G = G + 0.001 * convolve2d(G, LS, mode="same", boundary="symm")
    cax.set_array(G)
    return cax,

# Create the animation
animation = FuncAnimation(fig, update, frames=100, repeat=True, blit=True)

# Display the animation in the Jupyter notebook
from IPython.display import HTML
html = animation.to_jshtml(fps=20)
HTML(html)

Once Loop Reflect

G = F.copy()
for _ in range(1000):
    G = G + 0.001 * convolve2d(G, LS, mode="same", boundary="symm")

ax1,ax2 = PLOT2(F, G)
ax1.set_title("F")
ax2.set_title("F dispersed by heat equation")

Text(0.5, 1.0, 'F dispersed by heat equation')

PLOT(G)

Find Solutions to F=0

Checking for F=0 is not numerically stable. Instead we are trying to detect sign changes of F along the x or y direction.

G = F
Gs = np.sign(G)

PLOT(Gs)

from scipy.signal import convolve2d
import numpy as np
import matplotlib.pyplot as plt

# Sobel operator for edge detection
sobel_x = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
sobel_y = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])

Gx = convolve2d(DxF, sobel_x, mode='same', boundary='symm')
Gy = convolve2d(DxF, sobel_y, mode='same', boundary='symm')

# Magnitude of gradient (regions of zero-crossings)
magnitude = np.sqrt(Gx**2 + Gy**2)

# Apply a threshold to identify regions of potential zeros
threshold = magnitude.max() * 0.02
binary_image = magnitude > threshold

# Plotting
fig, axes = plt.subplots(1, 3, figsize=(18, 6))

# Plot F
axes[0].imshow(F, extent=[-W, W, -W, W], origin='lower', cmap='RdYlBu_r')
axes[0].set_title('F')
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')

# Plot magnitude of gradient
axes[1].imshow(magnitude, extent=[-W, W, -W, W], origin='lower', cmap='gray')
axes[1].set_title('Magnitude of Gradient')
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')

# Plot binary image
axes[2].imshow(binary_image, extent=[-W, W, -W, W], origin='lower', cmap='gray')
axes[2].set_title('Potential Zeros of F')
axes[2].set_xlabel('x')
axes[2].set_ylabel('y')

plt.tight_layout()
plt.show()

def detect_zero_crossings(F):
    # Compute the sign matrix of F
    sign_matrix = np.sign(F)

    # Check for sign changes along rows and columns
    horizontal_changes = np.diff(sign_matrix, axis=1)
    vertical_changes = np.diff(sign_matrix, axis=0)

    # Find where the changes are not zero. This indicates a zero crossing.
    zero_crossing_rows, zero_crossing_cols = np.where(horizontal_changes != 0)
    zero_crossing_cols += 1  # Adjust columns due to the diff operation

    zero_crossing_rows_v, zero_crossing_cols_v = np.where(vertical_changes != 0)
    zero_crossing_rows_v += 1  # Adjust rows due to the diff operation

    # Combine horizontal and vertical zero crossings
    zero_crossing_rows = np.concatenate((zero_crossing_rows, zero_crossing_rows_v))
    zero_crossing_cols = np.concatenate((zero_crossing_cols, zero_crossing_cols_v))

    return zero_crossing_rows, zero_crossing_cols

zero_rows, zero_cols = detect_zero_crossings(DxF)

plt.imshow(DxF, extent=[-W, W, -W, W], origin='lower', cmap='RdYlBu_r')
plt.scatter(x[zero_cols], y[zero_rows], color='black', s=.0001)  # Plot zero crossings
plt.colorbar()
plt.title('Zero Crossings of f(x,y)')
plt.xlabel('x')
plt.ylabel('y')
plt.show()

from scipy.signal import convolve2d

kernel = np.array([-1/2, 0, 1/2])

# Compute DxF by convolving along columns
DxF = convolve2d(F, kernel[:, np.newaxis], mode='same', boundary='symm')
# Compute DyF by convolving along rows
DyF = convolve2d(F, kernel[np.newaxis, :], mode='same', boundary='symm')

import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
from IPython.display import display

def f(x, y):
    return (np.sin(x) + np.cos(y)) / (1 + x**2 + y**2)

def plot_heatmap(C):
    N = 800
    x = np.linspace(-5, 5, N)
    y = np.linspace(-5, 5, N)
    X, Y = np.meshgrid(x, y)
    Z = f(X, Y)

    Z[Z > C] = 1
    Z[Z < C] = -1

    plt.figure(figsize=(8,6))
    plt.imshow(Z, extent=[-5, 5, -5, 5], origin='lower', cmap='RdYlBu_r')
    plt.colorbar()
    plt.title('Heatmap of f(x,y) with C = {:.2f}'.format(C))
    plt.xlabel('x')
    plt.ylabel('y')
    plt.show()

widgets.interactive(plot_heatmap, C=widgets.FloatSlider(value=0.23, min=-2, max=2, step=0.01, description='C:'))

#@CX-CUTOFF

# !pip install -q ipynbname
import ipynbname
nb_name = ipynbname.path().name

%%script bash -s "$nb_name"
#
# Export Notebook as Blogpost
# 
# WARNING: This will operate on the last saved version!
#
./render.sh "./$1" ./nb.html
cat <<EOF > fm.html
---
title: Convolution Algebras and Discrete Geometry - a Study
date: 2023-10-23
tags: post, math
url: /conv23
outputs:
  - RawHTML
---
EOF

cat fm.html nb.html > out.html
cp out.html "/host/root/home/hhartmann/src/HeinrichHartmann.github.io/hugo/content/posts/raw/conv-2023.html"

[NbConvertApp] Converting notebook ./2023-10-21 Geometry.ipynb to HTML
[NbConvertApp] WARNING | Alternative text is missing on 7 image(s).
[NbConvertApp] Writing 6455735 bytes to nb.html