4 files changed, 195 insertions, 0 deletions
diff --git a/code/sunlab/common/mathlib/__init__.py b/code/sunlab/common/mathlib/__init__.py
new file mode 100644
index 0000000..9b5ed21
--- /dev/null
+++ b/code/sunlab/common/mathlib/__init__.py
@@ -0,0 +1 @@
+from .base import *
diff --git a/code/sunlab/common/mathlib/base.py b/code/sunlab/common/mathlib/base.py
new file mode 100644
index 0000000..38ab14c
--- /dev/null
+++ b/code/sunlab/common/mathlib/base.py
@@ -0,0 +1,57 @@
+import numpy as np
+
+
+def angle(a, b):
+    """# Get Angle Between Row Vectors"""
+    from numpy import arctan2, pi
+
+    theta_a = arctan2(a[:, 1], a[:, 0])
+    theta_b = arctan2(b[:, 1], b[:, 0])
+    d_theta = theta_b - theta_a
+    assert (-pi <= d_theta) and (d_theta <= pi), "Theta difference outside of [-π,π]"
+    return d_theta
+
+
+def normalize(column):
+    """# Normalize Column Vector"""
+    from numpy.linalg import norm
+
+    return column / norm(column, axis=0)
+
+
+def winding(xy_grid, trajectory_start, trajectory_end):
+    """# Get Winding Number on Grid"""
+    from numpy import zeros, cross, clip, arcsin
+
+    trajectories = trajectory_end - trajectory_start
+    winding = zeros((xy_grid.shape[0]))
+    for idx, trajectory in enumerate(trajectories):
+        r = xy_grid - trajectory_start[idx]
+        cross = cross(normalize(trajectory), normalize(r))
+        cross = clip(cross, -1, 1)
+        theta = arcsin(cross)
+        winding += theta
+    return winding
+
+
+def vorticity(xy_grid, trajectory_start, trajectory_end):
+    """# Get Vorticity Number on Grid"""
+    from numpy import zeros, cross
+
+    trajectories = trajectory_end - trajectory_start
+    vorticity = zeros((xy_grid.shape[0]))
+    for idx, trajectory in enumerate(trajectories):
+        r = xy_grid - trajectory_start[idx]
+        vorticity += cross(normalize(trajectory), normalize(r))
+    return vorticity
+
+
+def data_range(data):
+    """# Get the range of values for each row"""
+    from numpy import min, max
+
+    return min(data, axis=0), max(data, axis=0)
+
+
+np.normalize = normalize
+np.range = data_range
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
diff --git a/code/sunlab/common/mathlib/random_walks.py b/code/sunlab/common/mathlib/random_walks.py
new file mode 100644
index 0000000..3aa3bcb
--- /dev/null
+++ b/code/sunlab/common/mathlib/random_walks.py
@@ -0,0 +1,83 @@
+def get_levy_flight(T=50, D=2, t0=0.1, alpha=3, periodic=False):
+    from numpy import vstack
+    from mistree import get_levy_flight as get_flight
+
+    if D == 2:
+        x, y = get_flight(T, mode="2D", periodic=periodic, t_0=t0, alpha=alpha)
+        xy = vstack([x, y]).T
+    elif D == 3:
+        x, y, z = get_flight(T, mode="3D", periodic=periodic, t_0=t0, alpha=alpha)
+        xy = vstack([x, y, z]).T
+    else:
+        raise ValueError(f"Dimension {D} not supported!")
+    return xy
+
+
+def get_levy_flights(N=10, T=50, D=2, t0=0.1, alpha=3, periodic=False):
+    from numpy import moveaxis, array
+
+    trajectories = []
+    for _ in range(N):
+        xy = get_levy_flight(T=T, D=D, t0=t0, alpha=alpha, periodic=periodic)
+        trajectories.append(xy)
+    return moveaxis(array(trajectories), 0, 1)
+
+
+def get_jitter_levy_flights(
+    N=10, T=50, D=2, t0=0.1, alpha=3, periodic=False, noise=5e-2
+):
+    from numpy.random import randn
+
+    trajectories = get_levy_flights(
+        N=N, T=T, D=D, t0=t0, alpha=alpha, periodic=periodic
+    )
+    return trajectories + randn(*trajectories.shape) * noise
+
+
+def get_gaussian_random_walk(T=50, D=2, R=5, step_size=0.5, soft=None):
+    from numpy import array, sin, cos, exp, zeros, pi
+    from numpy.random import randn, uniform, rand
+    from numpy.linalg import norm
+
+    def is_in(x, R=1):
+        from numpy.linalg import norm
+
+        return norm(x) < R
+
+    X = zeros((T, D))
+    for t in range(1, T):
+        while True:
+            if D == 2:
+                angle = uniform(0, pi * 2)
+                step = randn(1) * step_size
+                X[t, :] = X[t - 1, :] + array([cos(angle), sin(angle)]) * step
+            else:
+                X[t, :] = X[t - 1, :] + randn(D) / D * step_size
+            if soft is None:
+                if is_in(X[t, :], R):
+                    break
+            elif rand() < exp(-(norm(X[t, :]) - R) * soft):
+                break
+    return X
+
+
+def get_gaussian_random_walks(N=10, T=50, D=2, R=5, step_size=0.5, soft=None):
+    from numpy import moveaxis, array
+
+    trajectories = []
+    for _ in range(N):
+        xy = get_gaussian_random_walk(T=T, D=D, R=R, step_size=step_size, soft=soft)
+        trajectories.append(xy)
+    return moveaxis(array(trajectories), 0, 1)
+
+
+def get_gaussian_sample(T=50, D=2):
+    from numpy.random import randn
+
+    return randn(T, D)
+
+
+def get_gaussian_samples(N=10, T=50, D=2, R=5, step_size=0.5):
+    from numpy.random import randn
+
+    return randn(T, N, D)