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]:
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:
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')
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.
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')
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]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: