Coherent set comparison on Bickley jet¶
First, we set up two methods which allow us to simulate the Bickley jet forward and backward in time.
[1]:
import numpy as np
from deeptime.data import BickleyJet, bickley_jet
state = np.random.RandomState(seed=123) # global random state
def draw_initial_positions(n):
# draws n initial positions in the system's domain
X = np.vstack((state.uniform(0, 20, (n,)),
state.uniform(-3, 3, (n, ))))
return X.T
def forward_transform(X, h=1e-3):
# simulator with positive time-step
simulator = BickleyJet(h=h, n_steps=int(1. / h / 10))
# returns a trajectory evolving X to t=40
return simulator.trajectory(t0=0, x0=X, length=401)
def back_transform(Y, h=1e-3):
# simulator with negative time-step
simulator = BickleyJet(h=-h, n_steps=int(1. / h / 10))
# returns a trajectory evolving Y to t=0
return simulator.trajectory(t0=40, x0=Y, length=401)
Now we can simulate the Bickley jet using particle configurations.
[2]:
Xinit_train = draw_initial_positions(3000)
traj_train = forward_transform(Xinit_train)
print("Simulated training data with {} particles for {} "
"timesteps in {} dimensions.".format(*traj_train.shape))
Xinit_test = draw_initial_positions(15000)
traj_test = forward_transform(Xinit_test)
print("Simulated test data with {} particles for {} "
"timesteps in {} dimensions.".format(*traj_test.shape))
Simulated training data with 3000 particles for 401 timesteps in 2 dimensions.
Simulated test data with 15000 particles for 401 timesteps in 2 dimensions.
[3]:
from IPython.display import HTML
import matplotlib.pyplot as plt
import matplotlib.animation as animation
f, ax = plt.subplots(1, 1, figsize=(12, 4))
ax.set_xlim([0, 20])
ax.set_ylim([-4, 4])
c = traj_train[:, 0, 0] / traj_train[:, 0, 0].max()
handle = ax.scatter(*traj_train[:, 0].T, c=c)
ax.set_aspect('equal')
def step(t):
handle.set_offsets(traj_train[:, t])
return [handle]
ani = animation.FuncAnimation(f, step, interval=80, blit=True, repeat=False,
frames=traj_train.shape[1])
plt.close() # prevent figure from showing in output
HTML(ani.to_html5_video())
[3]: