aboutsummaryrefslogtreecommitdiff
path: root/code/sunlab/common/mathlib/lyapunov.py
blob: 3c747f1c573bc8753c4780b1564d9423c2eb0e5d (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
def trajectory_to_distances(x):
    """X: [N,N_t,N_d]
    ret [N,N_t]"""
    from numpy import zeros
    from numpy.linalg import norm
    from itertools import product, combinations

    x = [x[idx, ...] for idx in range(x.shape[0])]
    pairwise_trajectories = combinations(x, 2)
    _N_COMB = len(list(pairwise_trajectories))
    N_max = x[0].shape[0]
    distances = zeros((_N_COMB, N_max))
    pairwise_trajectories = combinations(x, 2)
    for idx, (a_t, b_t) in enumerate(pairwise_trajectories):
        distances[idx, :] = norm(a_t[:N_max, :] - b_t[:N_max, :], axis=-1)
    return distances


def Lyapunov_d(X):
    """X: [N,N_t]
    λ_n = ln(|dX_n|/|dX_0|)/n; n = [1,2,...]"""
    from numpy import zeros, log, repeat

    Y = zeros((X.shape[0], X.shape[1] - 1))
    Y = log(X[:, 1:] / repeat([X[:, 0]], Y.shape[1], axis=0).T) / (
        repeat([range(Y.shape[1])], Y.shape[0], axis=0) + 1
    )
    return Y


def Lyapunov_t(X):
    """X: [N,N_t,N_d]"""
    return Lyapunov_d(trajectory_to_distances(X))


Lyapunov = Lyapunov_d


def RelativeDistance_d(X):
    """X: [N,N_t]
    λ_n = ln(|dX_n|/|dX_0|)/n; n = [1,2,...]"""
    from numpy import zeros, log, repeat

    Y = zeros((X.shape[0], X.shape[1] - 1))
    Y = log(X[:, 1:] / repeat([X[:, 0]], Y.shape[1], axis=0).T)
    return Y


def RelativeDistance_t(X):
    """X: [N,N_t,N_d]"""
    return RelativeDistance_d(trajectory_to_distances(X))


RelativeDistance = RelativeDistance_d