1import os
  2import time
  3
  4import matplotlib
  5matplotlib.use('Agg')
  6import matplotlib.pyplot as plt
  7import matplotlib.font_manager as fm
  8import matplotlib.colors as mcolors
  9
 10import numpy as np
 11from numba import njit
 12
 13@njit(fastmath=True)
 14def compute_attractor(xs, ys, a, b, c, d):
 15    x, y = 0.1, 0.0
 16    for i in range(xs.shape[0]):
 17        xn = np.sin(a * y) + c * np.cos(a * x)
 18        yn = np.sin(b * x) + d * np.cos(b * y)
 19        x, y = xn, yn
 20        xs[i] = x
 21        ys[i] = y
 22
 23def main():
 24    libertinus_math_path = os.path.expanduser('~/.local/share/fonts/LibertinusMath-Regular.otf')
 25    font_name = 'Libertinus Math'
 26
 27    if os.path.exists(libertinus_math_path):
 28        try:
 29            fm.fontManager.addfont(libertinus_math_path)
 30            print(f"Registered font: {libertinus_math_path}")
 31            matplotlib.rcParams['font.family'] = [font_name, 'serif']
 32            matplotlib.rcParams['mathtext.fontset'] = 'custom'
 33            matplotlib.rcParams['mathtext.rm'] = font_name
 34            matplotlib.rcParams['mathtext.it'] = font_name
 35            matplotlib.rcParams['mathtext.bf'] = font_name
 36            print(f"Matplotlib configured to use '{font_name}' for equations.")
 37        except Exception as e:
 38            print(f"Error registering or configuring {font_name}: {e}")
 39            print("Falling back to 'serif' for equation text.")
 40            matplotlib.rcParams['font.family'] = ['serif']
 41            matplotlib.rcParams['mathtext.fontset'] = 'cm'
 42    else:
 43        print(f"Warning: {font_name} font not found at {libertinus_math_path}. Ensure it's installed or check the path.")
 44        print("Falling back to 'serif' for equation text.")
 45        matplotlib.rcParams['font.family'] = ['serif']
 46        matplotlib.rcParams['mathtext.fontset'] = 'cm'
 47
 48
 49    n_points = 300_000_000
 50    burn_in  =  10_000
 51    a, b, c, d = 1.8, -1.9, 1.0, 1.5
 52
 53    xs = np.empty(n_points, dtype=np.float64)
 54    ys = np.empty(n_points, dtype=np.float64)
 55
 56    t0 = time.perf_counter()
 57    compute_attractor(xs, ys, a, b, c, d)
 58    t1 = time.perf_counter()
 59
 60    xs = xs[burn_in:]
 61    ys = ys[burn_in:]
 62
 63    print(f"Computed {n_points} points in {t1 - t0:.2f} s.")
 64    print(f"x ∈ [{xs.min():.3f}, {xs.max():.3f}], y ∈ [{ys.min():.3f}, {ys.max():.3f}]")
 65
 66    dark_bg = '#111111'
 67    fig = plt.figure(figsize=(12, 8), facecolor=dark_bg)
 68    ax  = fig.add_axes((0, 0, 1, 1), facecolor=dark_bg)
 69    x_min, x_max = -2.5, 2.5
 70    y_min, y_max = -2.5, 2.5
 71    ax.set_xlim(x_min, x_max)
 72    ax.set_ylim(y_min, y_max)
 73    ax.axis('off')
 74
 75    H, xedges, yedges = np.histogram2d(xs, ys, bins=2000, range=[[x_min, x_max], [y_min, y_max]])
 76
 77    ax.imshow(
 78        H.T,
 79        origin='lower',
 80        extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]],
 81        cmap='inferno',
 82        norm=mcolors.LogNorm(vmin=1),
 83        aspect='auto',
 84        interpolation='bilinear'
 85    )
 86
 87    eq_text = (
 88        r"$x_{t+1} = \sin(%.1f y_t) + %.1f \cos(%.1f x_t)$" "\n"
 89        r"$y_{t+1} = \sin(%.1f x_t) + %.1f \cos(%.1f y_t)$"
 90    ) % (a, c, a, b, d, b)
 91
 92    fig.text(
 93        0.01, 0.99, eq_text,
 94        color='#BF40BF',
 95        fontsize=16,
 96        va='top',
 97        ha='left',
 98        alpha=0.8
 99    )
100
101    interval_text = (
102        r"$x \in [%.3f, %.3f]$" "\n"
103        r"$y \in [%.3f, %.3f]$"
104    ) % (xs.min(), xs.max(), ys.min(), ys.max())
105
106    fig.text(
107        0.01, 0.92,
108        interval_text,
109        color='#BF40BF',
110        fontsize=14,
111        va='top',
112        ha='left',
113        alpha=0.8
114    )
115
116    outpath = 'clifford_attractor.png'
117
118    plt.savefig(
119        outpath,
120        dpi=300,
121        bbox_inches='tight',
122        pad_inches=0,
123        facecolor=fig.get_facecolor()
124    )
125    plt.close(fig)
126
127if __name__ == '__main__':
128    main()