1 files changed, 54 insertions, 0 deletions
diff --git a/code/sunlab/common/mathlib/lyapunov.py b/code/sunlab/common/mathlib/lyapunov.py
new file mode 100644
index 0000000..3c747f1
--- /dev/null
+++ b/code/sunlab/common/mathlib/lyapunov.py
@@ -0,0 +1,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