From b85ee9d64a536937912544c7bbd5b98b635b7e8d Mon Sep 17 00:00:00 2001 From: Christian C Date: Mon, 11 Nov 2024 12:29:32 -0800 Subject: Initial commit --- code/sunlab/common/mathlib/lyapunov.py | 54 ++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) create mode 100644 code/sunlab/common/mathlib/lyapunov.py (limited to 'code/sunlab/common/mathlib/lyapunov.py') 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 -- cgit v1.2.1