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/__init__.py | 1 + code/sunlab/common/mathlib/base.py | 57 ++++++++++++++++++++ code/sunlab/common/mathlib/lyapunov.py | 54 +++++++++++++++++++ code/sunlab/common/mathlib/random_walks.py | 83 ++++++++++++++++++++++++++++++ 4 files changed, 195 insertions(+) create mode 100644 code/sunlab/common/mathlib/__init__.py create mode 100644 code/sunlab/common/mathlib/base.py create mode 100644 code/sunlab/common/mathlib/lyapunov.py create mode 100644 code/sunlab/common/mathlib/random_walks.py (limited to 'code/sunlab/common/mathlib') 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) -- cgit v1.2.1