In [1]:
import holoviews as hv
import matplotlib.pyplot as plt
import numpy as np
from swktools.plotting import plot3d
import swktools.TakensEmbed as te
from scipy.integrate import odeint

hv.extension('bokeh', 'matplotlib')


Model system


In [2]:
def run_model(Bxy, Byx, rx, ry, x0, y0, N):
    x = np.zeros(N)
    y = np.zeros(N)
    x[0], y[0] = x0, y0
    for i in range(1,N):
        x[i] = x[i-1]*(rx-rx*x[i-1]-Bxy*y[i-1])
        y[i] = y[i-1]*(ry-ry*y[i-1]-Byx*x[i-1])
    return x,y

In [32]:
%%opts Curve {+axiswise} Curve.time [width=900]
Bxy = 0.02
Byx = 0.2
rx = 3.8
ry = 3.5
x0 = 0.4
y0 = 0.2
N = 3000
x,y = run_model(Bxy, Byx, rx, ry, x0, y0, N)
data = np.array([x,y]).T
fig = hv.Curve(zip(x, y))*hv.Scatter(zip(x, y))+hv.Curve(x, group='time')*hv.Curve(y, group='time')
fig.cols(1)


Out[32]:

Get delay manifolds


In [33]:
tau = 1
ndelay = 3
delay_man = te.get_delayed_manifold(data, tau, ndelay)

In [34]:
%matplotlib nbagg
plot3d(delay_man[0,:], ntraj=10)



In [81]:
%matplotlib inline

Do embedding


In [39]:
rnge = range(10, 3000, 40)
cors = te.do_embedding(delay_man, rnge, timesampling='random')

In [40]:
hv.Curve(cors[0,1], label='x->y')*hv.Curve(cors[1,0], label='y->x')


Out[40]:

Do above for Lorentz attractor


In [18]:
sigma = 10
rho = 28
beta = 8.0/3
theta = 3 * np.pi / 4

def lorenz(xyz, t):
    x, y, z = xyz
    x_dot = sigma * (y - x)
    y_dot = x * rho - x * z - y
    z_dot = x * y - beta* z
    return [x_dot, y_dot, z_dot]

initial = (-10, -7, 35)
t = np.arange(0, 100, 0.006)

solution = odeint(lorenz, initial, t)

x = solution[:, 0]
y = solution[:, 1]
z = solution[:, 2]

xprime = np.cos(theta) * x - np.sin(theta) * y

lorenzian = hv.Overlay([hv.Path(d) for d in zip(np.array_split(xprime, 7), np.array_split(z, 7))])

In [19]:
%matplotlib nbagg
plot3d(solution, ntraj=10)


Get delay manifolds


In [20]:
tau = 10
ndelay = 3
delay_man = te.get_delayed_manifold(solution, tau, ndelay)

In [21]:
%matplotlib nbagg
ntraj = 10
fig = plt.figure()
ax = fig.add_subplot(131, projection='3d')
plot3d(delay_man[0,:], ntraj=ntraj, ax=ax)
ax = fig.add_subplot(132, projection='3d')
plot3d(delay_man[1,:], ntraj=ntraj, ax=ax)
ax = fig.add_subplot(133, projection='3d')
plot3d(delay_man[2,:], ntraj=ntraj, ax=ax)


Do embedding


In [26]:
rnge = range(40, 3000, 40)
cors = te.do_embedding(delay_man, rnge, timesampling='random')

In [27]:
fig = hv.Curve(cors[0,1], label='x->y')*hv.Curve(cors[1,0], label='y->x')
fig+= hv.Curve(cors[0,2], label='x->z')*hv.Curve(cors[2,0], label='z->x')
fig+= hv.Curve(cors[1,2], label='y->z')*hv.Curve(cors[2,1], label='z->y')
fig


Out[27]:

In [ ]: