Figure 1 — Exploring the Logistic Map. Bifurcation structure of the logistic map as the growth parameter r varies. Cobweb plots illustrate transitions: convergence to fixed points, period doubling around r ≈ 3.2, and onset of chaos near r ≈ 3.8.

import numpy as np
import matplotlib.pyplot as plt

# Logistic map function
def logistic_map(r, x0, steps):
    x = [x0]
    for _ in range(steps - 1):
        x.append(r * x[-1] * (1 - x[-1]))
    return x

# Cobweb plot function
def plot_cobweb(ax, r, x0, steps):
    x = np.linspace(0, 1, 500)
    y = r * x * (1 - x)

    ax.plot(x, y, 'k', label='f(x)')
    ax.plot(x, x, 'g--', label='y = x')

    x_n = x0
    for _ in range(steps):
        x_next = r * x_n * (1 - x_n)
        ax.plot([x_n, x_n], [x_n, x_next], 'b', lw=1)
        ax.plot([x_n, x_next], [x_next, x_next], 'b', lw=1)
        x_n = x_next

    ax.set_title(f'Cobweb Plot (r = {r}, x₀ = {x0})')
    ax.set_xlabel('x[n]')
    ax.set_ylabel('x[n+1]')
    ax.legend()
    ax.grid(True)

# Parameters
r = 3.8
x0 = 0.5
steps = 50

# Generate bifurcation data for one r
x_values = logistic_map(r, x0, steps)

# Create side-by-side plots
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Left: Iterated values over time
ax1.plot(range(steps), x_values, marker='o', linestyle='-', color='black')
ax1.set_title(f'Iterated Logistic Map Values (r = {r})')
ax1.set_xlabel('Iteration (n)')
ax1.set_ylabel('x[n]')
ax1.grid(True)

# Right: Cobweb diagram
plot_cobweb(ax2, r, x0, steps)

plt.tight_layout()
plt.show()

Figure 2 — Delay Embedding with GPR (r = 3.8). Delay embedding reconstructs the state space from a univariate time series, and Gaussian Process Regression (GPR) approximates the one-step map in a chaotic regime.

!pip install -q POT

import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel
from mpl_toolkits.mplot3d import Axes3D  # noqa: F401
import ot  # POT: Python Optimal Transport

# --------------------------------
# Logistic Map Generator
# --------------------------------
def logistic_map(r, x0, n_points):
    x = np.zeros(n_points)
    x[0] = x0
    for i in range(1, n_points):
        x[i] = r * x[i-1] * (1 - x[i-1])
    return x

# --------------------------------
# Delay Embedding Function
# --------------------------------
def delay_embedding(series, tau, Q):
    n_vectors = len(series) - (Q - 1) * tau
    return np.array([series[i:i + Q * tau:tau] for i in range(n_vectors)])

# --------------------------------
# GPR: Train and Iterate
# --------------------------------
def train_and_iterate_gpr(train_x, train_y, x0, n_steps, alpha=1e-4):
    kernel = RBF(length_scale=1.0) + WhiteKernel(noise_level=alpha)
    gpr = GaussianProcessRegressor(kernel=kernel, alpha=alpha)
    gpr.fit(train_x.reshape(-1, 1), train_y)

    generated = [x0]
    for _ in range(n_steps - 1):
        x_next = gpr.predict(np.array([[generated[-1]]]))[0]
        generated.append(x_next)
    return np.array(generated)

# --------------------------------
# Wasserstein Distance (POT)
# Returns W2 (order=2) by default
# --------------------------------
def wasserstein_distance_pot(points1, points2, order=2):
    """
    Compute Wasserstein-p distance between two point clouds using POT.
    points1, points2: arrays of shape (n1, d) and (n2, d)
    order: 1 for W1, 2 for W2 (default)
    """
    X = np.asarray(points1, dtype=float)
    Y = np.asarray(points2, dtype=float)
    assert X.ndim == 2 and Y.ndim == 2, "Point sets must be 2D arrays [n_points, dim]"

    n1, n2 = len(X), len(Y)
    a = np.ones(n1) / n1  # uniform weights
    b = np.ones(n2) / n2

    # Cost matrix
    D = ot.dist(X, Y)  # pairwise Euclidean distances
    if order == 2:
        M = D ** 2
        emd_cost = ot.emd2(a, b, M)   # this is the minimum squared cost
        return np.sqrt(emd_cost)       # return W2
    elif order == 1:
        M = D
        return ot.emd2(a, b, M)        # this is W1
    else:
        raise ValueError("Only order=1 (W1) or order=2 (W2) are supported.")

# --------------------------------
# Parameters
# --------------------------------
r = 3.8
x0 = 0.4
n_points = 1000
train_size = 10
tau = 1
Qs = [2, 3]
alpha = 1e-4

# --------------------------------
# Generate Data
# --------------------------------
original_series = logistic_map(r, x0, n_points)
train_x = original_series[:train_size]
train_y = original_series[1:train_size + 1]
gpr_series = train_and_iterate_gpr(train_x, train_y, train_x[0], n_points, alpha=alpha)

# --------------------------------
# Plot Embeddings and Compute Wasserstein Distance
# --------------------------------
fig = plt.figure(figsize=(12, 8))
colors = ['tab:blue', 'tab:orange']
wasserstein_results = {}

for i, Q in enumerate(Qs):
    emb_orig = delay_embedding(original_series, tau, Q)
    emb_gpr = delay_embedding(gpr_series, tau, Q)

    # Compute W2 distance between point clouds
    w2 = wasserstein_distance_pot(emb_orig, emb_gpr, order=2)
    wasserstein_results[Q] = w2

    # Plotting
    if Q == 2:
        ax1 = fig.add_subplot(2, 2, i * 2 + 1)
        ax2 = fig.add_subplot(2, 2, i * 2 + 2)
        ax1.plot(emb_orig[:, 0], emb_orig[:, 1], '.', alpha=0.5, color=colors[0], label='Original')
        ax2.plot(emb_gpr[:, 0], emb_gpr[:, 1], '.', alpha=0.5, color=colors[1], label='GPR')
        ax1.set_xlabel('x(n)')
        ax1.set_ylabel('x(n+1)')
        ax2.set_xlabel('x(n)')
        ax2.set_ylabel('x(n+1)')
        ax1.set_title(f'Original Logistic Map (Q={Q})')
        ax2.set_title(f'GPR Reconstructed (Q={Q})\nWasserstein-2 Distance: {w2:.4f}')
    elif Q == 3:
        ax1 = fig.add_subplot(2, 2, i * 2 + 1, projection='3d')
        ax2 = fig.add_subplot(2, 2, i * 2 + 2, projection='3d')
        ax1.scatter(emb_orig[:, 0], emb_orig[:, 1], emb_orig[:, 2], s=1, alpha=0.5, color=colors[0])
        ax2.scatter(emb_gpr[:, 0], emb_gpr[:, 1], emb_gpr[:, 2], s=1, alpha=0.5, color=colors[1])
        ax1.set_xlabel('x(n)')
        ax1.set_ylabel('x(n+1)')
        ax1.set_zlabel('x(n+2)')
        ax2.set_xlabel('x(n)')
        ax2.set_ylabel('x(n+1)')
        ax2.set_zlabel('x(n+2)')
        ax1.set_title(f'Original Logistic Map (Q={Q})')
        ax2.set_title(f'GPR Reconstructed (Q={Q})\nWasserstein-2 Distance: {w2:.4f}')

plt.tight_layout()
plt.show()

# --------------------------------
# Print Wasserstein Results
# --------------------------------
print("Wasserstein-2 Distances:")
for Q, dist in wasserstein_results.items():
    print(f"Q={Q}: {dist:.6f}")

   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 0.0/901.7 kB ? eta -:--:--
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸━━━━━━ 757.8/901.7 kB 22.0 MB/s eta 0:00:01
   ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 901.7/901.7 kB 16.6 MB/s eta 0:00:00

/usr/local/lib/python3.12/dist-packages/sklearn/gaussian_process/kernels.py:442: ConvergenceWarning: The optimal value found for dimension 0 of parameter k2__noise_level is close to the specified lower bound 1e-05. Decreasing the bound and calling fit again may find a better value.
  warnings.warn(

Wasserstein-2 Distances:
Q=2: 0.001143
Q=3: 0.002913

Figure 3 — Sampling strategy across r. At each parameter value r, N = 10 evenly spaced samples are taken for training and evaluation, standardizing comparisons across periodic and chaotic regimes.

import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C
import warnings
from sklearn.exceptions import ConvergenceWarning
warnings.filterwarnings("ignore", category=ConvergenceWarning)
warnings.filterwarnings("ignore", message=".*optimal value.*upper bound.*")


# ------------------------------------------------------------------
# Reproducibility
# ------------------------------------------------------------------
np.random.seed(42)

# ------------------------------------------------------------------
# Logistic map
# ------------------------------------------------------------------
def logistic_map(x, r):
    return r * x * (1 - x)

# ------------------------------------------------------------------
# Figure + simulation params tuned to match the reference image
# ------------------------------------------------------------------
r_values   = np.linspace(2.5, 4.0, 900)   # start at 2.5 (include fixed point / first bifurcation)
iterations = 1000                         # total iterations per r
plot_last  = 100                          # plot only last N iters to show attractor
x0         = 0.2                          # initial condition for both true and GP

# Storage
orbit_true_r, orbit_true_x = [], []
orbit_gp_r,   orbit_gp_x   = [], []

# ------------------------------------------------------------------
# Loop over r: true orbit (red) and GP-learned orbit (black)
# ------------------------------------------------------------------
for r in r_values:
    # ----- True logistic orbit
    x = x0
    for i in range(iterations):
        x = logistic_map(x, r)
        if i >= (iterations - plot_last):
            orbit_true_r.append(r)
            orbit_true_x.append(x)

    # ----- Train GP on noiseless samples of f(x)=r x(1-x)
    # Using a reasonably dense grid improves agreement in chaotic regime
    X_train = np.linspace(0.0, 1.0, 50).reshape(-1, 1)
    y_train = logistic_map(X_train, r).ravel()

    # Smooth kernel (shorter length_scale tracks curvature well)
    kernel = C(1.0, (1e-3, 1e3)) * RBF(length_scale=0.12, length_scale_bounds=(1e-3, 1e3))
    gp = GaussianProcessRegressor(kernel=kernel, alpha=1e-8, normalize_y=False, n_restarts_optimizer=3)
    gp.fit(X_train, y_train)

    # ----- Iterate the learned GP map x_{t+1} = GP(x_t)
    xg = x0
    for i in range(iterations):
        xg = gp.predict(np.array([[xg]], dtype=float))[0]
        # keep in [0,1] (prevents rare numerical drift in chaotic regime)
        if xg < 0.0: xg = 0.0
        elif xg > 1.0: xg = 1.0

        if i >= (iterations - plot_last):
            orbit_gp_r.append(r)
            orbit_gp_x.append(xg)

# ------------------------------------------------------------------
# Plot styling to match the sample
# ------------------------------------------------------------------
plt.figure(figsize=(16, 8), dpi=110)

# True logistic map (red)
plt.plot(orbit_true_r, orbit_true_x, '.', color='crimson', markersize=0.3, label='True logistic map')

# GP learned map (black/gray, slightly transparent)
plt.plot(orbit_gp_r, orbit_gp_x, '.', color='k', alpha=0.7, markersize=0.3, label='GP learned map')

plt.xlim(2.5, 4.0)
plt.ylim(0.0, 1.0)
plt.xlabel(r'$r$', fontsize=14)
plt.ylabel(r'$x$', fontsize=14)
plt.title('Bifurcation diagram of logistic map and GP learned dynamics', fontsize=24, pad=10)
leg = plt.legend(loc='upper left', frameon=True)
leg.get_frame().set_alpha(0.9)
plt.tight_layout()
plt.show()

Figure 4 — GPR on Logistic Map with only 5 data points. Models were trained with extremely sparse samples under increasing noise (σ = 0.0, 0.02, 0.05, 0.1). In (a), α = 0 forces interpolation through training points, leading to overfitting. In (b), a constant α = 0.005 yields smoother fits and better generalization. In real-world settings, true noise variance is usually unknown; testing multiple α values improves robustness.

# Figure 4(a)
import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C

# For reproducibility
np.random.seed(7)

# Logistic map function
def logistic_map(x, r=3.8):
    return r * x * (1 - x)

# Parameters
r = 3.8
N = 5
ㅈ = [0.0, 0.02, 0.05, 0.1]
X_grid = np.linspace(0, 1, 500).reshape(-1, 1)
y_true_grid = logistic_map(X_grid, r=r)

# Precompute random sample points and their noiseless output
X_sample_random = np.sort(np.random.rand(N)).reshape(-1, 1)
y_sample_random_true = logistic_map(X_sample_random, r=r).ravel()

# Precompute evenly spaced samples
X_sample_even = np.linspace(0, 1, N).reshape(-1, 1)
y_sample_even_true = logistic_map(X_sample_even, r=r).ravel()

# Plot
fig, axs = plt.subplots(len(noise_levels), 2, figsize=(12, 10), sharex=True, sharey=True)

for i, noise_std in enumerate(noise_levels):
    for j, sample_type in enumerate(['even', 'random']):
        if sample_type == 'even':
            X_sample = X_sample_even
            y_true_sample = y_sample_even_true
        else:
            X_sample = X_sample_random
            y_true_sample = y_sample_random_true

        # Add the same noise level to the same base samples
        y_sample_noisy = y_true_sample + np.random.normal(0, noise_std, size=y_true_sample.shape)

        # GPR setup
        kernel = C(1.0, (1e-3, 1e3)) * RBF(length_scale=0.2, length_scale_bounds=(1e-2, 1e2))
        gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10, alpha=noise_std**2 + 1e-10) #0.005 + 1e-10
        gp.fit(X_sample, y_sample_noisy)

        # Predict
        y_pred, sigma = gp.predict(X_grid, return_std=True)

        # Plotting
        ax = axs[i, j]
        ax.plot(X_grid, y_true_grid, 'r--', label='$f(x)$ (true)', linewidth=1.2)
        ax.plot(X_grid, y_pred, 'b-', label='GP fit', linewidth=1.5)
        ax.fill_between(X_grid.ravel(), y_pred - 1.96 * sigma, y_pred + 1.96 * sigma,
                        color='b', alpha=0.2, label='95% CI')
        ax.scatter(X_sample, y_sample_noisy, c='k', s=40, zorder=10, label='Noisy samples')

        if j == 0:
            ax.set_ylabel(f'Noise std = {noise_std}', fontsize=12)
        if i == 0:
            ax.set_title('Evenly Spaced Samples' if j == 0 else 'Randomly Selected Samples', fontsize=14)

        ax.set_xlim(0, 1)
        ax.set_ylim(0, 1.3)

        if i == len(noise_levels) - 1 and j == 1:
            ax.legend(loc='lower center', fontsize=8)

plt.tight_layout()
plt.show()

Figure 4 — GPR on Logistic Map with only 5 data points. Models were trained with extremely sparse samples under increasing noise (σ = 0.0, 0.02, 0.05, 0.1). In (a), α = 0 forces interpolation through training points, leading to overfitting. In (b), a constant α = 0.005 yields smoother fits and better generalization. In real-world settings, true noise variance is usually unknown; testing multiple α values improves robustness

# Figure 4(b)

import numpy as np
import matplotlib.pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C

# For reproducibility
np.random.seed(7)

# Logistic map function
def logistic_map(x, r=3.8):
    return r * x * (1 - x)

# Parameters
r = 3.8
N = 5
noise_levels = [0.0, 0.02, 0.05, 0.1]
X_grid = np.linspace(0, 1, 500).reshape(-1, 1)
y_true_grid = logistic_map(X_grid, r=r)

# Precompute random sample points and their noiseless output
X_sample_random = np.sort(np.random.rand(N)).reshape(-1, 1)
y_sample_random_true = logistic_map(X_sample_random, r=r).ravel()

# Precompute evenly spaced samples
X_sample_even = np.linspace(0, 1, N).reshape(-1, 1)
y_sample_even_true = logistic_map(X_sample_even, r=r).ravel()

# Plot
fig, axs = plt.subplots(len(noise_levels), 2, figsize=(12, 10), sharex=True, sharey=True)

for i, noise_std in enumerate(noise_levels):
    for j, sample_type in enumerate(['even', 'random']):
        if sample_type == 'even':
            X_sample = X_sample_even
            y_true_sample = y_sample_even_true
        else:
            X_sample = X_sample_random
            y_true_sample = y_sample_random_true

        # Add the same noise level to the same base samples
        y_sample_noisy = y_true_sample + np.random.normal(0, noise_std, size=y_true_sample.shape)

        # GPR setup
        kernel = C(1.0, (1e-3, 1e3)) * RBF(length_scale=0.2, length_scale_bounds=(1e-2, 1e2))
        gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10, alpha=0.005 + 1e-10) #0.005 + 1e-10
        gp.fit(X_sample, y_sample_noisy)

        # Predict
        y_pred, sigma = gp.predict(X_grid, return_std=True)

        # Plotting
        ax = axs[i, j]
        ax.plot(X_grid, y_true_grid, 'r--', label='$f(x)$ (true)', linewidth=1.2)
        ax.plot(X_grid, y_pred, 'b-', label='GP fit', linewidth=1.5)
        ax.fill_between(X_grid.ravel(), y_pred - 1.96 * sigma, y_pred + 1.96 * sigma,
                        color='b', alpha=0.2, label='95% CI')
        ax.scatter(X_sample, y_sample_noisy, c='k', s=40, zorder=10, label='Noisy samples')

        if j == 0:
            ax.set_ylabel(f'Noise std = {noise_std}', fontsize=12)
        if i == 0:
            ax.set_title('Evenly Spaced Samples' if j == 0 else 'Randomly Selected Samples', fontsize=14)

        ax.set_xlim(0, 1)
        ax.set_ylim(0, 1.3)

        if i == len(noise_levels) - 1 and j == 1:
            ax.legend(loc='lower center', fontsize=8)

plt.tight_layout()
plt.show()