我将二维相空间投影到图上并尝试粗粒度嵌入。我正在尝试执行 Kaplan-Glass 确定性测试。链接：http ://www.medicine.mcgill.ca/physio/glasslab/pub_pdf/coarse-grained_1993.pdf

轨迹通过每个框j一次或多次。我需要使用索引k对通过每个框j的个体进行索引，其中路径包含从R^2 中轨迹进入框j的点到下一个离开框j的点的轨迹段。然后，构造一个单位长度的向量V_(k,j) ，其方向由路径k 进入框j的点到轨迹下一个离开框的点之间的位移给出。轨迹向量V_(k,j) 指向轨迹的割线方向，即切线向量dx/dt的有限差分逼近方向。有没有一种pythonic方法来实现这一点？

到目前为止，这是我的代码：

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from teaspoon.SP.tsa_tools import takens
from sklearn.manifold import SpectralEmbedding
import re


def rossler_system(a, b, c, t, tf, h):
    def derivative(r, t):
        x = r[0]
        y = r[1]
        z = r[2]
        return np.array([- y - z, x + a * y, b + z * (x - c)])

    time = np.array([])
    x = np.array([])
    y = np.array([])
    z = np.array([])
    r = np.array([0.1, 0.1, 0.1])

    while (t <= tf):
        time = np.append(time, t)
        z = np.append(z, r[2])
        y = np.append(y, r[1])
        x = np.append(x, r[0])

        k1 = h * derivative(r, t)
        k2 = h * derivative(r + k1 / 2, t + h / 2)
        k3 = h * derivative(r + k2 / 2, t + h / 2)
        k4 = h * derivative(r + k3, t + h)
        r += (k1 + 2 * k2 + 2 * k3 + k4) / 6

        t = t + h

    return x, y, z

def unit_vector(v):
    return v / np.linalg.norm(v)


if __name__ == "__main__":
    # Rossler System
    a = 0.2
    b = 0.2
    c = 5.7
    t = 0
    tf = 100
    h = 0.01

    x, y, z = rossler_system(a, b, c, t, tf, h)

    # Create Redundant Phase Space 
    taken_matrix = takens(x, n=14, tau=1)  # n=14 is from the correlation dimension

    # Apply Lapalcian EigenMap to get the Reconstructed Phase Space
    taken_matrix = SpectralEmbedding(n_components=3, affinity="rbf").fit_transform(taken_matrix)

    # Apply PCA for a 2-D plot
    taken_matrix = PCA(n_components=2).fit_transform(taken_matrix)

    # COARSE GRAIN PHASE SPACE # 
    x_min, x_max = np.min(taken_matrix[:, 0]), np.max(taken_matrix[:, 0])
    y_min, y_max = np.min(taken_matrix[:, 1]), np.max(taken_matrix[:, 1])

    # - get the number of zeroes after decimal place
    x_min_dec = len(re.search('\d+\.(0*)', str(x_min)).group(1))
    x_max_dec = len(re.search('\d+\.(0*)', str(x_max)).group(1))
    y_min_dec = len(re.search('\d+\.(0*)', str(y_min)).group(1))
    y_max_dec = len(re.search('\d+\.(0*)', str(y_max)).group(1))

    x_dist_dec = len(re.search('\d+\.(0*)', str(x_max - x_min)).group(1))
    y_dist_dec = len(re.search('\d+\.(0*)', str(y_max - y_min)).group(1))

    # - makes equal sized boxes
    if x_dist_dec < y_dist_dec:
        x_min = x_min - (1/(10**(x_min_dec + 1)))
        x_max = x_max + (1/(10**(x_max_dec + 2)))
        x_dist = x_max - x_min
        y_dist = y_max - y_min
        y_min = y_min - (x_dist - y_dist)/2
        y_max = y_max + (x_dist - y_dist)/2
    elif x_dist_dec > y_dist_dec:
        y_min = y_min - (1/(10**(y_min_dec + 1)))
        y_max = y_max + (1/(10**(y_max_dec + 2)))
        y_dist = y_max - y_min
        x_dist = x_max - x_min
        x_min = x_min - (y_dist - x_dist)/2
        x_max = x_max + (y_dist - x_dist)/2

    boxes = 10
    x_vals = np.linspace(x_min, x_max, num=boxes)
    y_vals = np.linspace(y_min, y_max, num=boxes)

    xx, yy = np.meshgrid(x_vals, y_vals)  # box coordinates

    fig = plt.figure(dpi=50) # plot of 2D phase space
    plt.plot(taken_matrix[:, 0], taken_matrix[:, 1], c="red", lw=0.75, alpha=0.8)
    plt.scatter(xx, yy, c="red")
    for i, j in zip(x_vals, y_vals):
        plt.axvline(x=i, c="blue")
        plt.axhline(y=j, c="blue")
    plt.show()

    # Gradient of x and y values
    u = np.gradient(taken_matrix[:, 0])
    v = np.gradient(taken_matrix[:, 1])