In [ ]:
import numpy as np
from swktools.plotting import plot3d
import swktools.TakensEmbed as te
import holoviews as hv
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from sklearn.neighbors import NearestNeighbors
hv.extension('matplotlib','bokeh')
colors= hv.core.options.Cycle.default_cycles['default_colors']
Basis chaotic system from Sugihara et al (2012) - Detecting Causality in Complex Ecosystems
In [ ]:
%%opts Curve {+axiswise}
%%output backend='bokeh'
# parameters
rx = 3.8
ry = 3.5
bxy = .02
byx = .1
x1 = .4
y1 = .2
n = 3000
poisson = False
shift = False
sol = np.zeros((n,2))
sol[0,0] = x1
sol[0,1] = y1
for i in range(1,n):
sol[i,0] = sol[i-1,0]*(rx-rx*sol[i-1,0]-bxy*sol[i-1,1])
sol[i,1] = sol[i-1,1]*(ry-ry*sol[i-1,1]-byx*sol[i-1,0])
if poisson:
measT = 100
sol = np.random.poisson(sol*measT)/measT
if shift:
shift = 500
sol_shifted = np.zeros((n-shift, 2))
sol_shifted[:, 0] = sol[:-shift, 0]
sol_shifted[:, 1] = sol[shift:, 1]
sol = sol_shifted
x = sol[:, 0]
y = sol[:, 1]
lorenzian = hv.Overlay([hv.Path(d) for d in zip(np.array_split(x, 1), np.array_split(y, 1))])
lorenzian(style={'Path': dict(color=hv.Palette('Blues'), linewidth=1)})+hv.Curve(x)*hv.Curve(y)
Show delay manifolds
In [ ]:
%matplotlib nbagg
tau = 1 # how many time steps to go back
ndelay = 3 # how many dimensions to do for the delays
delayed = te.get_delayed_manifold(sol, tau, ndelay)
if ndelay == 3:
fig = plt.figure()
ntraj=20
ax = fig.add_subplot(131, projection='3d')
plot3d(delayed[0,::1,:], ntraj=ntraj, labels=['x(t)','x(t-tau)','x(t-2tau)'], ax=ax)
ax = fig.add_subplot(132, projection='3d')
plot3d(delayed[1,::1,:], ntraj=ntraj, labels=['y(t)','y(t-tau)','y(t-2tau)'], ax=ax)
Do cross-embedding
In [ ]:
%%output backend='bokeh'
reload(te)
cors = te.do_embedding(delayed, range(7, 2500, 40))
hv.Curve(cors[0,1,:], label='y|Mx')*hv.Curve(cors[1,0,:],label='x|My')
With randomized delay coordinates as in Tajima et al (2015)
In [ ]:
%%output backend='bokeh'
cors = te.do_embedding(delayed, range(7, 2500, 40), True)
hv.Curve(cors[0,1,:], label='y|Mx')*hv.Curve(cors[1,0,:],label='x|My')
In [ ]:
%%output backend='bokeh'
%%opts Curve {+axiswise}
sigma = 10
rho = 50
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, 10, 0.006)
N = len(t)
solution = odeint(lorenz, initial, t)
# solution = np.random.poisson((solution+40)*10)
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, 1), np.array_split(z, 1))])
fig = lorenzian(style={'Path': dict(color=hv.Palette('Blues'), linewidth=1)})
fig+= hv.Curve(x, label='x')*hv.Curve(y, label='y')*hv.Curve(z, label='z')
fig
In [ ]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
plot3d(solution[::1,:], ntraj=1, ax=ax)
In [ ]:
tau = 10 # how many time steps to go back
ndelay = 3 # how many dimensions to do for the delays
delayed = te.get_delayed_manifold(solution, tau, ndelay)
if ndelay == 3:
fig = plt.figure()
ntraj=20
ax = fig.add_subplot(131, projection='3d')
plot3d(delayed[0,::1,:], ntraj=ntraj, labels=['x(t)','x(t-tau)','x(t-2tau)'], ax=ax)
ax = fig.add_subplot(132, projection='3d')
plot3d(delayed[1,::1,:], ntraj=ntraj, labels=['y(t)','y(t-tau)','y(t-2tau)'], ax=ax)
ax = fig.add_subplot(133, projection='3d')
plot3d(delayed[2,::1,:], ntraj=ntraj, labels=['z(t)','z(t-tau)','z(t-2tau)'], ax=ax)
Normal cross-embedding
In [ ]:
%%output backend='bokeh'
cors = te.do_embedding(delayed, range(20, 1000, 20))
fig = hv.Curve(cors[0,1,:], label='y|Mx')*hv.Curve(cors[1,0,:], label='x|My')
fig+= hv.Curve(cors[0,2,:], label='z|Mx')*hv.Curve(cors[2,0,:], label='x|Mz')
fig+= hv.Curve(cors[1,2,:], label='z|My')*hv.Curve(cors[2,1,:], label='y|Mz')
fig
With randomized coordinates
In [ ]:
%%output backend='bokeh'
cors = te.do_embedding(delayed, range(20, 1000, 20),True)
fig = hv.Curve(cors[0,1,:], label='y|Mx')*hv.Curve(cors[1,0,:], label='x|My')
fig+= hv.Curve(cors[0,2,:], label='z|Mx')*hv.Curve(cors[2,0,:], label='x|Mz')
fig+= hv.Curve(cors[1,2,:], label='z|My')*hv.Curve(cors[2,1,:], label='y|Mz')
fig
In [ ]: