Convergent Cross Mapping

This notebook and package is reproducing the results from Detecting Causality in Complex Ecosystems



In [1]:

    
import sys
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('ticks')
sns.set_context(context='paper',font_scale=1.5)
%matplotlib inline



In [2]:

    
%load_ext autoreload
%autoreload 2



In [3]:

    
#alter the line below to correspond to your file system
nonlinpy_dir = '/Users/nickc/Documents/skccm'
sys.path.append(nonlinpy_dir)



In [4]:

    
import skccm.paper as paper
import skccm.data as data

Work Flow

A coupled logistic map. Here is what we are going to do:

Generate the time series
Calculate the mutual information of the time series
Embed the time series (not going to explore embedding dimension, just lag)
Analyze forecast skill for a range of library lengths

1. Generate the time series



In [5]:

    
rx1 = 3.72 #determines chaotic behavior of the x1 series
rx2 = 3.72 #determines chaotic behavior of the x2 series
b12 = 0.2 #Influence of x1 on x2
b21 = 0.01 #Influence of x2 on x1
ts_length = 1000
x1,x2 = data.coupled_logistic(rx1,rx2,b12,b21,ts_length)



In [6]:

    
fig,ax = plt.subplots(nrows=2,sharex=True,sharey=True)
ax[0].plot(x1[0:100])
ax[1].plot(x2[0:100])
ax[0].set_yticks([.1,.3,.5,.7,.9])
ax[1].set_xlabel('Time')
sns.despine()



In [7]:

    
#fig.savefig('../figures/train_test_split/coupled_logistic.png',bbox_inches='tight')

2. Calculate the mutual information



In [9]:

    
e1 = paper.Embed(x1)
e2 = paper.Embed(x2)



In [10]:

    
mi1 = e1.mutual_information(10)
mi2 = e2.mutual_information(10)



In [11]:

    
fig,ax = plt.subplots(nrows=2,sharex=True,sharey=True)
ax[0].plot(np.arange(1,11),mi1)
ax[1].plot(np.arange(1,11),mi2)
ax[1].set_xlabel('Lag')
sns.despine()



In [12]:

    
#fig.savefig('../figures/mutual_info.png',bbox_inches='tight')

3. Embed the time series



In [13]:

    
lag = 1
embed = 2
X1 = e1.embed_vectors_1d(lag,embed)
X2 = e2.embed_vectors_1d(lag,embed)



In [14]:

    
fig,ax = plt.subplots(ncols=2,sharey=True,sharex=True,figsize=(10,4)) 
ax[0].scatter(X1[:,0],X1[:,1])
ax[1].scatter(X2[:,0],X2[:,1])
ax[0].set_xlabel('X1(t)')
ax[0].set_ylabel('X1(t-1)')
ax[1].set_xlabel('X2(t)')
ax[1].set_ylabel('X2(t-1)')
sns.despine()



In [15]:

    
#fig.savefig('../figures/x_embedded.png',bbox_inches='tight')

4. Forecast skill as a function of library length

Same legend as above.



In [16]:

    
CCM = ccm.CCM()

Split it into a testing set and training set



In [18]:

    
from skccm.utilities import train_test_split



In [19]:

    
x1tr, x1te, x2tr, x2te = train_test_split(X1,X2, percent=.75)



In [22]:

    
len_tr = len(x1tr)
lib_lens = np.linspace(10, len_tr/2, dtype='int')
lib_lens









    Out[22]:





array([ 10,  17,  24,  32,  39,  47,  54,  62,  69,  76,  84,  91,  99,
       106, 114, 121, 129, 136, 143, 151, 158, 166, 173, 181, 188, 195,
       203, 210, 218, 225, 233, 240, 248, 255, 262, 270, 277, 285, 292,
       300, 307, 314, 322, 329, 337, 344, 352, 359, 367, 374])



In [23]:

    
CCM.fit(x1tr,x2tr)
x1p, x2p = CCM.predict(x1te, x2te,lib_lengths=lib_lens)



In [24]:

    
sc1,sc2 = CCM.score()



In [25]:

    
fig,ax = plt.subplots()
ax.plot(lib_lens,sc1,label='X1 xmap X2')
ax.plot(lib_lens,sc2, label='X2 xmap X1')
ax.set_xlabel('Library Length')
ax.set_ylabel('Forecast Skill')
sns.despine()



In [26]:

    
#fig.savefig('../figures/train_test_split/xmap_lib_len.png',bbox_inches='tight')

Reproducing the plot from the paper

The paper explores different values for b12 and b21 in Figure 3B.



In [29]:

    
rx1 = 3.7 #determines chaotic behavior of the x1 series
rx2 = 3.7 #determines chaotic behavior of the x2 series

ts_length = 1000

#variables for embedding
lag = 1
embed = 2

CCM = ccm.CCM() #intitiate the ccm object

#store values
num_b = 20 #number of bs to test
sc1_store = np.empty((num_b,num_b))
sc2_store = np.empty((num_b,num_b))

#values over which to test
max_b = .4 #maximum b values
b_range = np.linspace(0,max_b,num=num_b)



In [34]:

    
#loop through b values for b12 and b21
for ii,b12 in enumerate(b_range):
    for jj, b21 in enumerate(b_range):
        
        x1,x2 = data.coupled_logistic(rx1,rx2,b12,b21,ts_length)
        
        e1 = ccm.Embed(x1)
        e2 = ccm.Embed(x2)
        
        X1 = e1.embed_vectors_1d(lag,embed)
        X2 = e2.embed_vectors_1d(lag,embed)
        
        x1tr, x1te, x2tr, x2te = train_test_split(X1,X2, percent=.75)
        
        CCM.fit(x1tr,x2tr)
        x1p, x2p = CCM.predict(x1te, x2te,lib_lengths=[500]) #only one prediction
        sc1,sc2 = CCM.score()
        
        
        sc1_store[ii,jj] = sc1[0]
        sc2_store[ii,jj] = sc2[0]
        
sc_diff = sc2_store-sc1_store



In [35]:

    
fig,ax = plt.subplots()
v = np.linspace(-1.6,1.6,21)
cax = ax.contourf(b_range,b_range,sc_diff,v,cmap='magma')
fig.colorbar(cax,ticks=v[::2])
ax.set_xlabel('B12')
ax.set_ylabel('B21')









    Out[35]:





<matplotlib.text.Text at 0x1116bb3c8>



In [26]:

    
#fig.savefig('../figures/train_test_split/xmap_changingB.png',bbox_inches='tight')

Randomly forced



In [36]:

    
rx2 = 3.72 #determines chaotic behavior of the x2 series
b12 = .2 #Influence of x1 on x2

ts_length = 1000
x1,x2 = data.driven_rand_logistic(rx2,b12,ts_length)



In [37]:

    
fig,ax = plt.subplots(nrows=2,sharex=True,sharey=True)
ax[0].plot(x1[0:100])
ax[1].plot(x2[0:100])
ax[0].set_yticks([.1,.3,.5,.7,.9])
ax[1].set_xlabel('Time')
sns.despine()



In [39]:

    
e1 = ccm.Embed(x1)
e2 = ccm.Embed(x2)



In [40]:

    
mi1 = e1.mutual_information(10)
mi2 = e2.mutual_information(10)



In [41]:

    
fig,ax = plt.subplots(nrows=2,sharex=True,sharey=True)
ax[0].plot(np.arange(1,11),mi1)
ax[1].plot(np.arange(1,11),mi2)
ax[1].set_xlabel('Lag')
sns.despine()



In [42]:

    
lag = 1
embed = 3

X1 = e1.embed_vectors_1d(lag,embed)
X2 = e2.embed_vectors_1d(lag,embed)



In [43]:

    
fig,ax = plt.subplots(ncols=2,sharey=True,sharex=True,figsize=(10,4)) 
ax[0].scatter(X1[:,0],X1[:,1])
ax[1].scatter(X2[:,0],X2[:,1])
ax[0].set_xlabel('X1(t)')
ax[0].set_ylabel('X1(t-1)')
ax[1].set_xlabel('X2(t)')
ax[1].set_ylabel('X2(t-1)')
sns.despine()



In [45]:

    
CCM = ccm.CCM()



In [50]:

    
lib_lens = np.linspace(10,ts_length/2,dtype='int')

x1tr, x1te, x2tr, x2te = train_test_split(X1,X2, percent=.75)

CCM.fit(x1tr,x2tr)
x1p, x2p = CCM.predict(x1te, x2te,lib_lengths=lib_lens) #only one prediction
sc1,sc2 = CCM.score()



In [52]:

    
plt.plot(sc1)
plt.plot(sc2)









    Out[52]:





[<matplotlib.lines.Line2D at 0x111c142e8>]

Periodic Forcing



In [53]:

    
rx2 = 3.72 #determines chaotic behavior of the x2 series
b12 = .5 #Influence of x1 on x2

ts_length = 1000
x1,x2 = data.driving_sin(rx2,b12,ts_length)



In [54]:

    
fig,ax = plt.subplots(nrows=2,sharex=True)
ax[0].plot(x1[0:100])
ax[1].plot(x2[0:100])
ax[1].set_xlabel('Time')
sns.despine()



In [55]:

    
em_x1 = ccm.Embed(x1)
em_x2 = ccm.Embed(x2)



In [56]:

    
mi1 = em_x1.mutual_information(10)
mi2 = em_x2.mutual_information(10)



In [57]:

    
fig,ax = plt.subplots(nrows=2,sharex=True,sharey=True)
ax[0].plot(np.arange(1,11),mi1)
ax[1].plot(np.arange(1,11),mi2)
ax[1].set_xlabel('Lag')
sns.despine()



In [58]:

    
lag1 = 4
lag2 = 1
embed = 3

X1 = em_x1.embed_vectors_1d(lag1,embed)
X2 = em_x2.embed_vectors_1d(lag2,embed)



In [59]:

    
fig,ax = plt.subplots(ncols=2,sharey=True,sharex=True,figsize=(10,4)) 
ax[0].scatter(X1[:,0],X1[:,1])
ax[1].scatter(X2[:,0],X2[:,1])
ax[0].set_xlabel('X1(t)')
ax[0].set_ylabel('X1(t-1)')
ax[1].set_xlabel('X2(t)')
ax[1].set_ylabel('X2(t-1)')
sns.despine()



In [60]:

    
CCM = ccm.CCM()



In [63]:

    
print('X1 Shape:', X1.shape)
print('X2 Shape:', X2.shape)









    



X1 Shape: (992, 3)
X2 Shape: (998, 3)

Need to make them the same length before they can be used. This could be corrected in future versions.



In [64]:

    
X2 = X2[0:len(X1)]



In [65]:

    
lib_lens = np.linspace(10,ts_length/2,dtype='int')

x1tr, x1te, x2tr, x2te = train_test_split(X1,X2, percent=.75)

CCM.fit(x1tr,x2tr)
x1p, x2p = CCM.predict(x1te, x2te,lib_lengths=lib_lens) #only one prediction
sc1,sc2 = CCM.score()



In [66]:

    
plt.plot(sc1)
plt.plot(sc2)









    Out[66]:





[<matplotlib.lines.Line2D at 0x1111d7d30>]



In [71]:

    
x1p[-1].shape









    Out[71]:





(248, 3)



In [72]:

    
x1te.shape









    Out[72]:





(248, 3)



In [75]:

    
plt.plot(x1p[-1][:,0])
plt.plot(x1te[:,0])









    Out[75]:





[<matplotlib.lines.Line2D at 0x1112409b0>]



In [80]:

    
plt.plot(x2te[:,0],alpha=.5)
plt.plot(x2p[-1][:,0])









    Out[80]:





[<matplotlib.lines.Line2D at 0x110bbf550>]

Lagged Coupled Logistic



In [81]:

    
rx1 = 3.72 #determines chaotic behavior of the x1 series
rx2 = 3.72 #determines chaotic behavior of the x2 series
b12 = 0.01 #Influence of x1 on x2
b21 = 0.3 #Influence of x2 on x1
ts_length = 8000
x1,x2 = data.lagged_coupled_logistic(rx1,rx2,b12,b21,ts_length)



In [82]:

    
fig,ax = plt.subplots(nrows=2,sharex=True)
ax[0].plot(x1[0:100])
ax[1].plot(x2[0:100])
ax[1].set_xlabel('Time')
sns.despine()



In [84]:

    
em_x1 = ccm.Embed(x1)
em_x2 = ccm.Embed(x2)



In [85]:

    
mi1 = em_x1.mutual_information(10)
mi2 = em_x2.mutual_information(10)



In [86]:

    
fig,ax = plt.subplots(nrows=2,sharex=True,sharey=True)
ax[0].plot(np.arange(1,11),mi1)
ax[1].plot(np.arange(1,11),mi2)
ax[1].set_xlabel('Lag')
sns.despine()



In [87]:

    
lag1 = 2
lag2 = 2
embed = 3

X1 = em_x1.embed_vectors_1d(lag1,embed)
X2 = em_x2.embed_vectors_1d(lag2,embed)



In [88]:

    
fig,ax = plt.subplots(ncols=2,sharey=True,sharex=True,figsize=(10,4)) 
ax[0].scatter(X1[:,0],X1[:,1])
ax[1].scatter(X2[:,0],X2[:,1])
ax[0].set_xlabel('X1(t)')
ax[0].set_ylabel('X1(t-1)')
ax[1].set_xlabel('X2(t)')
ax[1].set_ylabel('X2(t-1)')
sns.despine()



In [89]:

    
CCM = ccm.CCM()



In [90]:

    
lib_lens = np.linspace(10,ts_length/2,dtype='int')

x1tr, x1te, x2tr, x2te = train_test_split(X1,X2, percent=.75)

CCM.fit(x1tr,x2tr)
x1p, x2p = CCM.predict(x1te, x2te,lib_lengths=lib_lens) #only one prediction
sc1,sc2 = CCM.score()



In [91]:

    
plt.plot(lib_lens,sc1)
plt.plot(lib_lens,sc2)









    Out[91]:





[<matplotlib.lines.Line2D at 0x111d527b8>]



In [ ]:

Random Linearly Increasing



In [92]:

    
ts_length=1000
x1 = np.random.randn(ts_length,) + np.linspace(0,25,ts_length)
x2 = np.random.randn(ts_length,) + np.linspace(0,10,ts_length)



In [93]:

    
plt.plot(x1)
plt.plot(x2)









    Out[93]:





[<matplotlib.lines.Line2D at 0x1129414a8>]



In [95]:

    
em_x1 = ccm.Embed(x1)
em_x2 = ccm.Embed(x2)



In [96]:

    
mi1 = em_x1.mutual_information(20)
mi2 = em_x2.mutual_information(20)



In [97]:

    
fig,ax = plt.subplots(nrows=2,sharex=True)
ax[0].plot(np.arange(1,21),mi1)
ax[1].plot(np.arange(1,21),mi2)
ax[1].set_xlabel('Lag')
sns.despine()



In [98]:

    
lag1 = 2
lag2 = 2
embed = 8

X1 = em_x1.embed_vectors_1d(lag1,embed)
X2 = em_x2.embed_vectors_1d(lag2,embed)



In [99]:

    
fig,ax = plt.subplots(ncols=2,sharey=True,sharex=True,figsize=(10,4)) 
ax[0].scatter(X1[:,0],X1[:,1])
ax[1].scatter(X2[:,0],X2[:,1])
ax[0].set_xlabel('X1(t)')
ax[0].set_ylabel('X1(t-1)')
ax[1].set_xlabel('X2(t)')
ax[1].set_ylabel('X2(t-1)')
sns.despine()



In [100]:

    
CCM = ccm.CCM()



In [101]:

    
lib_lens = np.linspace(10,ts_length/2,dtype='int')

x1tr, x1te, x2tr, x2te = train_test_split(X1,X2, percent=.75)

CCM.fit(x1tr,x2tr)
x1p, x2p = CCM.predict(x1te, x2te,lib_lengths=lib_lens) #only one prediction
sc1,sc2 = CCM.score()



In [102]:

    
plt.plot(lib_lens,sc1)
plt.plot(lib_lens,sc2)









    Out[102]:





[<matplotlib.lines.Line2D at 0x1128c1e10>]

Cosine and Sine Waves



In [105]:

    
ts_length=500
x1 = np.random.randn(ts_length,) + 10*np.sin(np.linspace(0,10*np.pi,ts_length))
x2 = np.random.randn(ts_length,) + 20*np.cos(np.linspace(0,10*np.pi,ts_length))

#x1[x1<0] = np.random.randn(np.sum(x1<0),)
#x2[x2<0] = np.random.randn(np.sum(x2<0),)



In [106]:

    
plt.plot(x1)
plt.plot(x2,alpha=.5)









    Out[106]:





[<matplotlib.lines.Line2D at 0x11294e6d8>]



In [112]:

    
em_x1 = ccm.Embed(x1)
em_x2 = ccm.Embed(x2)
mi1 = em_x1.mutual_information(100)
mi2 = em_x2.mutual_information(100)



In [114]:

    
fig,ax = plt.subplots(nrows=2,sharex=True,sharey=True)
ax[0].plot(np.arange(1,101),mi1)
ax[1].plot(np.arange(1,101),mi2)
ax[1].set_xlabel('Lag')
sns.despine()



In [115]:

    
lag1 = 20
lag2 = 20
embed = 2

X1 = em_x1.embed_vectors_1d(lag1,embed)
X2 = em_x2.embed_vectors_1d(lag2,embed)



In [116]:

    
fig,ax = plt.subplots(ncols=2,sharey=True,sharex=True,figsize=(10,4)) 
ax[0].scatter(X1[:,0],X1[:,1])
ax[1].scatter(X2[:,0],X2[:,1])
ax[0].set_xlabel('X1(t)')
ax[0].set_ylabel('X1(t-1)')
ax[1].set_xlabel('X2(t)')
ax[1].set_ylabel('X2(t-1)')
sns.despine()



In [117]:

    
CCM = ccm.CCM()



In [118]:

    
lib_lens = np.linspace(10,ts_length/2,dtype='int')

x1tr, x1te, x2tr, x2te = train_test_split(X1,X2, percent=.75)

CCM.fit(x1tr,x2tr)
x1p, x2p = CCM.predict(x1te, x2te,lib_lengths=lib_lens) #only one prediction
sc1,sc2 = CCM.score()



In [119]:

    
plt.plot(lib_lens,sc1)
plt.plot(lib_lens,sc2)









    Out[119]:





[<matplotlib.lines.Line2D at 0x114f23240>]