import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import CubicSpline, UnivariateSpline # ============================================================================= # Stochastic 2D Fractal with Iterative Spline Fitting # ============================================================================= # # This program demonstrates an iterative process on a stochastic 2D fractal curve. # At each iteration: # 1. Fit a smoothing spline to the current fractal points. # 2. The spline represents a "generalised" version of the fractal (smoother, # less prone to overfitting the stochastic noise). # 3. Generate a new fractal by adding residual noise (or fresh stochastic detail) # to the spline, mimicking a "fractal iteration" where the spline acts as # an attractor. # # The process is inspired by the idea of a generalisation limit: after enough # iterations, the curve stabilises to a smooth shape that captures the essential # low‑frequency structure while ignoring the high‑frequency random fluctuations. # This is analogous to the matrix fixed point in the previous ridge regression # example – here the attractor is a spline function instead of a weight vector. # # The term "stochastic 2D fractal" is implemented via the random midpoint # displacement algorithm, which produces a curve with fractal (self‑similar) # properties. # ============================================================================= def generate_midpoint_displacement(n_points, roughness=0.7, rng=None): """ Generate a 1D fractal curve using random midpoint displacement. Parameters: ----------- n_points : int Number of points (must be 2^k + 1 for perfect midpoint algorithm, but we pad to nearest power of two and interpolate). roughness : float Controls the fractal dimension (higher = rougher). rng : numpy.random.Generator, optional Random number generator. Returns: -------- x : ndarray x coordinates (linearly spaced from 0 to 1). y : ndarray y coordinates (fractal curve). """ if rng is None: rng = np.random.default_rng() # Find smallest power of two >= n_points-1, then add 1 m = 1 while m < n_points - 1: m <<= 1 full_len = m + 1 # Initialise array with endpoints y_full = np.zeros(full_len) y_full[0] = rng.normal(0, 1) y_full[-1] = rng.normal(0, 1) step = full_len - 1 amplitude = 1.0 while step > 1: half = step // 2 for i in range(half, full_len - 1, step): left = y_full[i - half] right = y_full[i + half] midpoint = 0.5 * (left + right) displacement = amplitude * rng.normal(0, 1) y_full[i] = midpoint + displacement amplitude *= roughness step = half # Interpolate to exact n_points if needed x_full = np.linspace(0, 1, full_len) x = np.linspace(0, 1, n_points) y = np.interp(x, x_full, y_full) return x, y def fit_smoothing_spline(x, y, smoothing_factor=0.1): """ Fit a smoothing spline to the data. Adjust `smoothing_factor` to control the trade‑off between fidelity and smoothness (higher = smoother). """ # UnivariateSpline with smoothing parameter s # s=0 interpolates, larger s gives smoother fit spline = UnivariateSpline(x, y, s=smoothing_factor, ext='raise') return spline def iterative_spline_fitting(x, y0, n_iter=5, smoothing_factor=0.5, noise_scale=0.2): """ Iteratively: - Fit a smoothing spline to the current curve. - Generate a new curve by adding scaled residual noise to the spline. The idea: each iteration refines the "fractal" but the spline fit regularises it, driving the curve towards a generalisation limit. """ history = [y0.copy()] current_y = y0.copy() for i in range(n_iter): # Fit spline to current curve spline = fit_smoothing_spline(x, current_y, smoothing_factor) y_smooth = spline(x) # Compute residuals (high‑frequency fractal details) residuals = current_y - y_smooth # Generate new curve as: smooth part + attenuated residuals + fresh noise # The attenuation (0.5) and fresh noise (noise_scale) mimic a # stochastic iteration that gradually loses high‑frequency components. new_y = y_smooth + 0.5 * residuals + noise_scale * np.random.randn(len(x)) history.append(new_y.copy()) current_y = new_y return history def plot_fractal_iterations(x, history, original_fractal=None): """Visualise the evolution of the curve through iterations.""" n_iter = len(history) - 1 # history[0] is initial, history[1..] are after each fit+noise fig, axes = plt.subplots(2, 2, figsize=(12, 8), sharex=True, sharey=True) axes = axes.flatten() # Plot original fractal if original_fractal is None: original_fractal = history[0] axes[0].plot(x, original_fractal, 'k-', alpha=0.7, label='Initial fractal') axes[0].set_title('Iteration 0: Stochastic Fractal') axes[0].legend() # Show intermediate iterations (1, 2, ...) for idx, i in enumerate([1, 2, 3], start=1): if i <= n_iter: axes[idx].plot(x, history[i], 'b-', alpha=0.7, label=f'Iteration {i} (spline + noise)') axes[idx].set_title(f'Iteration {i}') axes[idx].legend() # Final iteration axes[3].plot(x, history[-1], 'g-', label='Final attractor') axes[3].plot(x, original_fractal, 'k--', alpha=0.4, label='Original fractal') axes[3].set_title(f'Iteration {n_iter} (Generalisation limit)') axes[3].legend() for ax in axes: ax.set_xlabel('x') ax.set_ylabel('y') ax.grid(alpha=0.3) plt.suptitle('Iterative Spline Fitting on a Stochastic 2D Fractal', fontsize=14) plt.tight_layout() plt.show() # ============================================================================= # Main execution # ============================================================================= if __name__ == "__main__": # Parameters n_points = 257 # number of points along the curve (2^8+1 works nicely) roughness = 0.7 # fractal roughness (higher = more detail) n_iterations = 4 # number of spline‑fit iterations smoothing = 0.3 # smoothing factor for spline (larger = smoother) noise_scale = 0.08 # magnitude of fresh stochastic noise each iteration # Create a reproducible random generator rng = np.random.default_rng(42) # Generate initial stochastic fractal curve x, y_fractal = generate_midpoint_displacement(n_points, roughness, rng) # Run iterative spline fitting print("Performing iterative spline fitting on fractal curve...") history = iterative_spline_fitting(x, y_fractal, n_iter=n_iterations, smoothing_factor=smoothing, noise_scale=noise_scale) # Visualise results plot_fractal_iterations(x, history, original_fractal=y_fractal) # Optional: show the fixed‑point (final spline) without added noise # i.e., the pure generalisation limit final_spline = fit_smoothing_spline(x, history[-1], smoothing) y_final_smooth = final_spline(x) plt.figure(figsize=(10, 5)) plt.plot(x, y_fractal, 'k-', alpha=0.5, label='Original stochastic fractal') plt.plot(x, y_final_smooth, 'r-', linewidth=2, label='Generalisation limit (pure spline)') plt.xlabel('x') plt.ylabel('y') plt.title('From Stochastic Fractal to Smooth Generalisation Limit') plt.legend() plt.grid(alpha=0.3) plt.show() # Print diagnostic information print("\n--- Fractal dimension estimate (roughness) ---") print(f"Initial fractal roughness parameter: {roughness}") print("Higher roughness gives a more jagged, space‑filling curve.") print("\n--- Generalisation limit interpretation ---") print("The final spline is the attractor of this iterative process.") print("It represents a 'matrix‑like' fixed point that does not overfit") print("the high‑frequency stochastic details of the fractal.")