Ordinal Partition Network (OPN) approach for Poincaré section estimation of the Thomas attractor based on https://github.com/hinsley/MultimodalMaps/blob/main/maps/thomas.jl

opn_sections_thomas_attractor.py

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from itertools import permutations
from collections import Counter

# Define the Thomas attractor system
def thomas_derivatives(t, u, b):
    x, y, z = u
    dxdt = np.sin(y) - b * x
    dydt = np.sin(z) - b * y
    dzdt = np.sin(x) - b * z
    return [dxdt, dydt, dzdt]


# Parameters
b = 0.208186
u0 = [1.0, 1.0, 0.5]  # Initial conditions
t_span = (0.0, 10_000.0)  # Time span for the simulation
t_eval = np.linspace(t_span[0], t_span[1], 10_000)  # Evaluation times

# Solve the system
sol = solve_ivp(thomas_derivatives, t_span, u0, args=(b,), t_eval=t_eval, method='RK45')

# Extract the solution
x, y, z = sol.y

# Define functions to compute ordinal patterns
def get_ordinal_indices(m):
    perms = list(permutations(range(m)))
    lookup = {perm: i for i, perm in enumerate(perms)}
    return lookup

def ordinal_partition(x, m, tau):
    n = len(x)
    ordinal_indices = []
    for i in range(n - (m - 1) * tau):
        window = x[i:i + m * tau:tau]
        sorted_indices = tuple(np.argsort(window))
        ordinal_indices.append(sorted_indices)
    return ordinal_indices

# Parameters for ordinal pattern analysis
m = 3    # Embedding dimension
tau = 1  # Time delay
w = 1    # Window size

# Compute ordinal patterns
ordinal_lookup = get_ordinal_indices(m)
ordinal_patterns = ordinal_partition(x, m, tau)
ordinal_symbols = [ordinal_lookup[pattern] for pattern in ordinal_patterns]

# Identify the highest entropy ordinal symbol
symbol_counts = Counter(ordinal_symbols)
highest_entropy_symbol = symbol_counts.most_common(1)[0][0]

# Find first-of-consecutive occurrences of the highest entropy symbol
section_indices = []
for i in range(1, len(ordinal_symbols)):
    if ordinal_symbols[i] == highest_entropy_symbol and ordinal_symbols[i - 1] != highest_entropy_symbol:
        section_indices.append((i - 1) * w + 1)

# Extract the points in the section
section_points = [(x[i], y[i], z[i]) for i in section_indices]
section_x, section_y, section_z = zip(*section_points)

# Plot the 3D trajectory and section points
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot(x, y, z, lw=0.5, color='blue', label='Thomas Attractor Trajectory')
ax.scatter(section_x, section_y, section_z, color='red', s=10, label='Section Points')

# Labels and title
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
ax.set_title('Thomas Attractor with Section Points')
ax.legend()
plt.show()