Analysis of trajectories with regards to system dynamics

I produced sets of trajectories for different values of ecosystem and trade income to see whether they resulted in a stable, trade supported society or predator prey cycles for society and ecosystem.

Therefore, I want to calculate some measure, that tell me whether the trajectories are primarily noise around a stable value or wheter they are oscillating in some sort of chaotic attractor.

First, lets see how the trajectories are saved.


In [2]:
# imports
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import itertools as it
import seaborn as sns

from ipywidgets import IntProgress
from IPython.display import display

# from https://github.com/nikdon/pyEntropy
from pyentrp import entropy as entrp

In [2]:
# load data
path_to_data = '/home/jakob/Project_MayaSim/output_data/X10_trajectories/results/all_trjs.hd5'

with pd.HDFStore(path_to_data) as store:
    data = store.select('d1', where='step>500', columns=['total_population'])
data.index = data.index.droplevel('test')

In [172]:
# generate empty dataframe with same index (minus 'steps')
index = data.index
empty_df = pd.DataFrame(index=index, columns=['e']).sum(level=(0,1,2))
empty_df[['e']] = float('nan')

ann = pd.DataFrame(index=empty_df.index, columns=['ann'], data='').sum(level=(0, 1))

In [173]:
# calculate permutation entropy for trajectories for t>500 (after initial overshoot)
f = IntProgress(min=0, max=len(empty_df)) # instantiate a progress bar
display(f) # display the progress bar
i=0
for ind, row in empty_df.iterrows():
    # increase bar
    f.value += 1
    i += 1
    # check, if already calculated,
    if np.isnan(row['e']):
        # if not, load time series
        dft = data.loc[ind]

        # if timeseries contains data, calculate permutation entropy
        if len(dft.values[:,0])>0:
            e = entrp.permutation_entropy(dft.values[:,0], normalize=True)
        else:
            e = float('nan')
        # and save.
        empty_df.loc[ind, 'e'] = e



In [269]:
# save permutation entropy
path_to_data = '/home/jakob/Project_MayaSim/output_data/X10_trajectories/results/permutation_entropy.hd5'

with pd.HDFStore(path_to_data) as store:
    store.put('d1', 
              empty_df,
             format='table',
             data_columns=True)

In [176]:
tmp = empty_df.xs(key=(6100, 0.0001), level=('r_trade', 'r_es'))
index = pd.MultiIndex.from_product([[6000], [0.0001], list(tmp.index.values)], names=['r_trade', 'r_es', 'run_id'])
empty_df = empty_df.append(pd.DataFrame(index=index, data=tmp.values, columns=['e']))

tmp = empty_df.xs(key=(7100, 0.0001), level=('r_trade', 'r_es'))
index = pd.MultiIndex.from_product([[7000], [0.0001], list(tmp.index.values)], names=['r_trade', 'r_es', 'run_id'])
empty_df = empty_df.append(pd.DataFrame(index=index, data=tmp.values, columns=['e']))

Time series of total population for different values of ecosystem income and a fixed value of trade income:


In [186]:
# plot trajectories for different values of r_es and r_trade
%matplotlib inline
import seaborn as sns
import itertools as it

r_ess = [5e-5, 6e-5, 7e-5, 0.00012, 0.00017]
r_trades=[4500, 6800, 9000]
ann = pd.DataFrame(index=empty_df.index, columns=['ann'], data=' ').sum(level=(0, 1))

fig, ax = plt.subplots(nrows=len(r_trades), ncols=len(r_ess))
fig.set_figwidth(7*len(r_ess))
fig.set_figheight(5*len(r_trades))

with pd.HDFStore('/home/jakob/Project_MayaSim/output_data/X10_trajectories/results/all_trjs.hd5') as store:
    df = store.select('d1', where=f'r_trade = {r_trades} & r_es = {r_ess}', columns=['total_population'])
    df.index = df.index.droplevel('test')


for ((i, r_es), (j, r_trade)) in it.product(enumerate(r_ess), enumerate(r_trades)):    
    ann.loc[(r_trade, r_es), 'ann'] = 'X'
    try:
        df.loc[(r_trade, r_es)].unstack('run_id').plot(ax=ax[j,i], legend=False, color='black', alpha=.2)
        ax[j,i].set_title(f'r_es = {r_es}, r_trade = {r_trade}')
    except:
        print('no_data')
fig.savefig('trade_vs_es_trajectories.png')



In [187]:
ann


Out[187]:
ann
r_trade r_es
4000 0.00005
0.00006
0.00007
0.00008
0.00009
0.00010
0.00011
0.00012
0.00013
0.00014
0.00015
0.00016
0.00017
0.00018
4100 0.00005
0.00006
0.00007
0.00008
0.00009
0.00010
0.00011
0.00012
0.00013
0.00014
0.00015
0.00016
0.00017
0.00018
4200 0.00005
0.00006
... ... ...
9200 0.00005
0.00006
0.00007
0.00008
0.00009
0.00010
0.00011
0.00012
0.00013
0.00014
0.00015
0.00016
0.00017
0.00018
9300 0.00005
0.00006
0.00007
0.00008
0.00009
0.00010
0.00011
0.00012
0.00013
0.00014
0.00015
0.00016
0.00017
0.00018
6000 0.00010
7000 0.00010

756 rows × 1 columns

The trajectories show that for different values of ecosystem income and trade icome, different patterns arise. For very low ecosystem income, the population dies out. For intermediate ecosystem income, oscillatory patterns appear, and for a certain (hight) combination of ecosystem income and trade, a set of high population states seems to be stable (except for small noise around them).

Then, I can compare this intuitive analysis to the result of the permutation entropy measure. And I can see that:

  • for trajectories where population goes extinct, entropy is ~0
  • for oscillating trajectories, the permutation entropy is ~1/2
  • for stationary, noisy trajectories, the permutation entropy is ~1
  • and in the range where both dynamic behaviors occur, it shows sort of a bimodal distribution.

In [223]:
# plot heatmaps of permutation entropy for different values of r_es and r_trade
tmp = empty_df
fig, ax = plt.subplots(nrows=1, ncols=3)
fig.set_figwidth(20)
fig.set_figheight(6)

results = tmp.min(level=[0,1])
results.unstack('r_es')
sns.heatmap(results['e'].unstack('r_es'), fmt='s', cmap='viridis', ax=ax[0], annot=ann.unstack('r_es'))
ax[0].set_title('MIN permutation entropy')

# results = tmp.mean(level=[0,1])
# results.unstack('r_es')
# sns.heatmap(results['e'].unstack('r_es'), fmt='s', cmap='viridis', ax=ax[1], annot=ann.unstack('r_es'))
# ax[1].set_title('MEAN permutation entropy')

results = tmp.max(level=[0,1])
results.unstack('r_es')
sns.heatmap(results['e'].unstack('r_es'), fmt='s', cmap='viridis', ax=ax[1], annot=ann.unstack('r_es'))
ax[1].set_title('MAX permutation entropy')


results = tmp.max(level=[0,1]) - tmp.min(level=[0,1])
results.unstack('r_es')
sns.heatmap(results['e'].unstack('r_es'), fmt='s', cmap='viridis', ax=ax[2], annot=ann.unstack('r_es'))
ax[2].set_title('MAX difference in permutation entropy')

def add_at(ax, t, loc=2):
    from mpl_toolkits.axes_grid.anchored_artists import AnchoredText
    fp = dict(size=20)
    _at = AnchoredText(t, loc=loc, prop=fp, )
    ax.add_artist(_at)
    return _at

for x, no in zip(ax, ['A', 'B', 'C']):
    add_at(x, no, loc=4)
    for label in x.get_yaxis().get_ticklabels()[1::2]:
        label.set_visible(False)
        


fig.set_tight_layout(True)
fig.savefig('permutation_entropy.pdf')
fig.savefig('permutation entropy.png')


/home/jakob/anaconda3/envs/py36/lib/python3.6/site-packages/matplotlib/figure.py:2022: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  warnings.warn("This figure includes Axes that are not compatible "

In the white regions, the runs on the somehow didn't finish whithin 24h. Have to look into this.

These results also show, that the default parameters of r_trade=6000 and r_es=0.0001 would be a good set of parameters to test drought events with an oscillating system.

Bottom Line: The model shows four different regions in the r_es vs. r_trade parameter space with distinct boundaries.

  • convergence to zero total population
  • oscillatory behavior
  • convergence to a set of constant, high values for total population
  • coexistence of the latter two

In [45]:
# plot trajectories for different values of r_es and r_trade
%matplotlib inline
import itertools as it

import matplotlib
import matplotlib.cm as cm

norm = matplotlib.colors.Normalize(vmin=0., vmax=1., clip=True)
mapper = cm.ScalarMappable(norm=norm, cmap=cm.viridis)

txt = ['A', 'B', 'C', 'D']
r_ess = [5e-5, 9e-5, 0.00011, 0.00014]
r_trades=[7400,]

with pd.HDFStore('/home/jakob/Project_MayaSim/output_data/X10_trajectories/results/all_trjs.hd5') as store:
    df = store.select('d1', where=f'r_trade = {r_trades} & r_es = {r_ess}', columns=['total_population'])
    df.index = df.index.droplevel('test')
    
# load permutation entropy
path_to_data = '/home/jakob/Project_MayaSim/output_data/X10_trajectories/results/permutation_entropy.hd5'

with pd.HDFStore(path_to_data) as store:
    tmp = store.select('d1', where='r_trade>5000 & r_es < 0.00016')

ann = pd.DataFrame(index=tmp.index, columns=['ann'], data=' ').sum(level=(0,1))

In [111]:
fig, ax = plt.subplots(nrows=2, ncols=2)
fig.set_figwidth(7)
fig.set_figheight(5)

k = 0

def add_at(ax, t, loc=2):
    from mpl_toolkits.axes_grid.anchored_artists import AnchoredText
    fp = dict(size=24)
    _at = AnchoredText(t, loc=loc, prop=fp, )
    ax.add_artist(_at)
    return _at

for ((j, r_trade), (i, r_es)) in it.product(enumerate(r_trades), enumerate(r_ess)):  
    axij = ax[i%2, i//2]
    add_at(axij, txt[k], loc=1)
    ann.loc[(r_trade, r_es), 'ann'] = f'{txt[k]}'
    k += 1
    rids = tmp.index.levels[2]
    col = tmp.loc[(r_trade, r_es)].unstack('run_id')
    dta = df.loc[(r_trade, r_es)].unstack('run_id')
    dta.columns = dta.columns.droplevel(0)
    for rid in rids:
        color = mapper.to_rgba(col[[rid]].values)
        ln = dta[[rid]].plot(ax=axij, color=color, legend=False, lw=.5)
    if i%2 == 0:
        axij.set_xticklabels([])
        axij.set_xlabel('')
    if i//2 == 0:
        axij.set_ylabel('population')

# add colorbar
cb_ax = fig.add_axes([.9, 0.1, 0.03, 0.85])
mapper.set_array(np.arange(0, 1, 100))
cbar = fig.colorbar(mapper, cax = cb_ax, ax=ax[0])
cb_ax.set_ylabel('permutation entropy')

fig.tight_layout(rect=[0, 0.0, .9, 1.])

fig.savefig('trade_vs_es_trajectories.pdf')


/home/jakob/anaconda3/envs/py36/lib/python3.6/site-packages/ipykernel/__main__.py:20: PerformanceWarning: indexing past lexsort depth may impact performance.

In [110]:
mins = tmp.min(level=(0,1))
    

classification = pd.DataFrame(columns = ['min'], index=mins.index, data=mins.values)
classification['max'] = tmp['e'].max(level=[0,1])
classification['dif'] = classification['max'] - classification['min']

classification['cls'] = float('nan')

classification.loc[classification['min']>.7, 'cls'] = 4
classification.loc[classification['min']<.7, 'cls'] = 2
classification.loc[classification['dif']>.5, 'cls'] = 1
classification.loc[[(.19<x<.5) & (y>.5) & (z>.9) for x, y, z in zip(classification['dif'], classification['min'], classification['max'])], 'cls'] = 3
# classification.loc[[(x>.4) & (y<.77) for x, y in zip(classification['min'], classification['max'])], 'cls'] = 2

fig, ax = plt.subplots()
fig.set_figwidth(6)
fig.set_figheight(5)
ax.set(ylim=(5000, 9200))
cmap = sns.color_palette("RdBu_r", 4)

sns.heatmap(classification['cls'].unstack('r_es'), ax=ax, cmap=cmap, annot=ann.unstack('r_es'), fmt='s', annot_kws={'fontsize':11})

# use every second x and y tick label
for label in ax.get_xaxis().get_ticklabels()[1::2]:
    label.set_visible(False)

for label in ax.get_yaxis().get_ticklabels()[1::2]:
    label.set_visible(False)
        
# invert y axis
ax.invert_yaxis()

# adjust color bar
cax = plt.gcf().axes[-1]
cax.get_yaxis().set_ticks([])
for j, lab in enumerate(['1','2','3', '4']):
    cax.text(.6, (2 * j + 1) / 8.0, lab, ha='center', va='center', rotation=90)

fig.savefig('classified_dynamics.pdf')



In [251]:
# plot heatmaps of permutation entropy for different values of r_es and r_trade
tmp = empty_df
fig, ax = plt.subplots(nrows=3, ncols=1)
fig.set_figwidth(20)
fig.set_figheight(60)

results = tmp.min(level=[0,1])
results.unstack('r_es')
sns.heatmap(results['e'].unstack('r_es'), fmt='g', cmap='viridis', ax=ax[0], annot=results['e'].unstack('r_es'))
ax[0].set_title('MIN permutation entropy')

# results = tmp.mean(level=[0,1])
# results.unstack('r_es')
# sns.heatmap(results['e'].unstack('r_es'), fmt='s', cmap='viridis', ax=ax[1], annot=ann.unstack('r_es'))
# ax[1].set_title('MEAN permutation entropy')

results = tmp.max(level=[0,1])
results.unstack('r_es')
sns.heatmap(results['e'].unstack('r_es'), fmt='g', cmap='viridis', ax=ax[1], annot=results['e'].unstack('r_es'))
ax[1].set_title('MAX permutation entropy')


results = tmp.max(level=[0,1]) - tmp.min(level=[0,1])
results.unstack('r_es')
sns.heatmap(results['e'].unstack('r_es'), fmt='g', cmap='viridis', ax=ax[2], annot=results['e'].unstack('r_es'))
ax[2].set_title('MAX difference in permutation entropy')

def add_at(ax, t, loc=2):
    from mpl_toolkits.axes_grid.anchored_artists import AnchoredText
    fp = dict(size=20)
    _at = AnchoredText(t, loc=loc, prop=fp, )
    ax.add_artist(_at)
    return _at

for x, no in zip(ax, ['A', 'B', 'C']):
    add_at(x, no, loc=4)
    for label in x.get_yaxis().get_ticklabels()[1::2]:
        label.set_visible(False)



In [319]:
ann


Out[319]:
ann
r_trade r_es
5100 0.00005
0.00006
0.00007
0.00008
0.00009
0.00010
0.00011
0.00012
0.00013
5200 0.00005
0.00006
0.00007
0.00008
0.00009
0.00010
0.00011
0.00012
0.00013
5300 0.00005
0.00006
0.00007
0.00008
0.00009
0.00010
0.00011
0.00012
0.00013
5400 0.00005
0.00006
0.00007
... ... ...
9000 0.00013
9100 0.00005
0.00006
0.00007
0.00008
0.00009
0.00010
0.00011
0.00012
0.00013
9200 0.00005
0.00006
0.00007
0.00008
0.00009
0.00010
0.00011
0.00012
0.00013
9300 0.00005
0.00006
0.00007
0.00008
0.00009
0.00010
0.00011
0.00012
0.00013
6000 0.00010
7000 0.00010

387 rows × 1 columns


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]: