<jupyter_start><jupyter_text>Stability of defection, optimisation of strategies and the limits of memory in the PD.----------------------  2. Stability of defection<jupyter_code>import opt_mo

import numpy as np
import sympy as sym
import itertools
import axelrod as axl
import matplotlib.pyplot as plt
from fractions import Fraction
from axelrod.action import Action

C, D = Action.C, Action.D
import tqdm
import matplotlib.transforms as transforms
import matplotlib
import warnings; warnings.simplefilter('ignore')
import operator
import functools
sym.init_printing()
font = {"size": 10, "weight": "bold"}
matplotlib.rc("font", **font)<jupyter_output><empty_output><jupyter_text>Analytical check----------------<jupyter_code>p_1, p_2, p_3, p_4 = sym.symbols("p_1, p_2, p_3, p_4")
q_1, q_2, q_3, q_4 = sym.symbols("q_1, q_2, q_3, q_4")
k_1, k_2, k_3, k_4 = sym.symbols("k_1, k_2, k_3, k_4")

p, q, k = (p_1, p_2, p_3, p_4), (q_1, q_2, q_3, q_4), (k_1, k_2, k_3, k_4)
def get_Q_N_derivative(player, opponent):
    x = np.array(player)
    
    Q = opt_mo.utility.quadratic_term_numerator(opponent)
    c = opt_mo.utility.linear_term_numerator(opponent)
        
    return np.dot(x, Q) + c
def get_Q_N(player, opponent):
    x = np.array(player)

    Q = opt_mo.utility.quadratic_term_numerator(opponent)
    c = opt_mo.utility.linear_term_numerator(opponent)
    a = opt_mo.utility.constant_term_numerator(opponent)

    return np.dot(x, Q.dot(x.T) * 1 / 2) + np.dot(c, x.T) + a
def get_Q_D(player, opponent):
    x = np.array(player)

    Q_bar = opt_mo.utility.quadratic_term_denominator(opponent)
    c_bar = opt_mo.utility.linear_term_denominator(opponent)
    a_bar = opt_mo.utility.constant_term_denominator(opponent)

    return np.dot(x, Q_bar.dot(x.T) * 1 / 2) + np.dot(c_bar, x.T) + a_bar
def get_Q_D_derivative(player, opponent):
    x = np.array(player)

    Q_bar = opt_mo.utility.quadratic_term_denominator(opponent)
    c_bar = opt_mo.utility.linear_term_denominator(opponent)
    a_bar = opt_mo.utility.constant_term_denominator(opponent)

    return np.dot(x, Q_bar) + c_bar<jupyter_output><empty_output><jupyter_text>**Check quadratic derivative of**$$\frac{1}{2}pQp^T + cp +a $$<jupyter_code>expr = get_Q_N(p, q)
diff = [sym.diff(expr, i) for i in p]
derivatives = get_Q_N_derivative(p, q)
for i in range(4):
    assert (diff[i] - derivatives[i]).simplify() == 0<jupyter_output><empty_output><jupyter_text>**Check derivative of utility**<jupyter_code>def get_written_derivative_of_utility(player, opponents):
    sums = 0
    for opponent in opponents:
        numerator = (get_Q_N_derivative(player, opponent) * get_Q_D(player, opponent) -
                      get_Q_D_derivative(player, opponent) * get_Q_N(player, opponent))
        
        denominator = get_Q_D(player, opponent) ** 2
        
        sums += numerator / denominator
    return (sums) * (1 / len(opponents))
for seed in range(100):
    num_players = 5
    np.random.seed(seed)
    opponents = [[np.random.random() for _ in range(4)] for _ in range(num_players)]
    
    np.random.seed(seed + 1000)
    player = [np.random.random() for _ in range(4)]
    
    written_derivative = get_written_derivative_of_utility(player, opponents)
    
    utility = opt_mo.tournament_utility(p, opponents)
    utility_derivative = [sym.diff(utility, i) for i in p]
    utility_derivative = [expr.subs({p_1: player[0], p_2: player[1], p_3: player[2], p_4: player[3]}) for expr
                          in utility_derivative]
    
    differences = written_derivative - utility_derivative
    for difference in differences:
        assert np.isclose(round(difference, 10), 0)
opponents = [q, k]
utility = opt_mo.tournament_utility(p, opponents)
derivative_of_utility = [sym.diff(utility, i) for i in p]
written_derivative = get_written_derivative_of_utility(p, opponents)
for i in tqdm.tqdm(range(4)):
    assert (written_derivative[i] - derivative_of_utility[i]).simplify() == 0<jupyter_output>0%|          | 0/4 [00:00<?, ?it/s][A
 25%|██▌       | 1/4 [00:06<00:18,  6.28s/it][A
 50%|█████     | 2/4 [02:27<01:33, 46.69s/it][A
 75%|███████▌  | 3/4 [02:33<00:34, 34.56s/it][A
100%|██████████| 4/4 [02:34<00:00, 24.63s/it][A<jupyter_text>**Stability of defection** **Check condition for defection stability**<jupyter_code>opponents = [q, k]
utility = opt_mo.tournament_utility(p, opponents)
diff_utility = [sym.diff(utility, i) for i in p]
diff_utility_at_zero = [expr.subs({p_1: 0, p_2: 0, p_3: 0, p_4: 0}).expand() for expr in diff_utility]
def get_derivate_for_p_zeros(opponents):
    sums = 0
    for opponent in opponents:
        lhs = opt_mo.utility.linear_term_numerator(opponent) * opt_mo.utility.constant_term_denominator(opponent)
        rhs = opt_mo.utility.linear_term_denominator(opponent) * opt_mo.utility.constant_term_numerator(opponent)
        
        denominator = opt_mo.utility.constant_term_denominator(opponent) ** 2
        
        sums +=  (lhs - rhs) / denominator

    return (sums) * (1 / len(opponents))
expression = get_derivate_for_p_zeros(opponents)
for i in tqdm.tqdm(range(4)):
    assert (diff_utility_at_zero[i] - expression[i]).simplify() == 0
def stackplot(plot, eco, logscale=True):

    populations = eco.population_sizes
    figure, ax = plt.subplots(figsize=(10, 8))
#     figure = ax.get_figure()
    figure.patch.set_facecolor('#ffffff')
    figure.patch.set_alpha(0.2)

    turns = range(len(populations))
    pops = [
        [populations[iturn][ir] for iturn in turns]
        for ir in plot.result_set.ranking
    ]
    ax.stackplot(turns, *pops)

    ax.yaxis.tick_left()
    ax.yaxis.set_label_position("right")
    ax.yaxis.labelpad = 25.0

    ax.set_ylim([0.0, 1.0])
    #ax.set_xlim([0.0, 10 ** 3])
    ax.set_ylabel("Relative population size", fontweight='bold')
    ax.set_xlabel("Turn", fontweight='bold')

    trans = transforms.blended_transform_factory(ax.transAxes, ax.transData)
    ticks = []
    for i, n in enumerate(plot.result_set.ranked_names):
        x = -0.01
        y = (i + 0.5) * 1 / plot.result_set.num_players
        if n != 'Defector':
            opponent_strings = (n.split(':')[1].replace('[', '').replace(']', '')).split(",")
            opponent = [Fraction(float(op)).limit_denominator() for op in opponent_strings]
            label = '$q^{(%s)}$'% i
#             for p in range(3):
#                 label += r'\frac{' + str(opponent[p].numerator) + '}{' + str(opponent[p].denominator) + '},'
#             label += r'\frac{' + str(opponent[3].numerator) + '}{' + str(opponent[3].denominator) + '})$'
#             print(label)
            n = label
        ax.annotate(
            n,
            xy=(x, y),
            xycoords=trans,
            clip_on=False,
            va="center",
            ha="right",
            fontsize=10,
        )
        ticks.append(y)
    ax.set_yticks(ticks)
    ax.tick_params(direction="out")
    ax.set_yticklabels([])    
    if logscale:
        ax.set_xscale("log")
    for tick in ax.yaxis.get_majorticklabels():  # example for xaxis
        tick.set_fontsize(20) 
        
    plt.tight_layout()
    plt.xlim(0, 10 **2)
    return figure.savefig('../img/population_defection_takes_over.pdf', facecolor=figure.get_facecolor(),
                          edgecolor='none', bbox_inches='tight')
for seed in range(2, 10):
    np.random.seed(seed)
    opponents = [[np.round(np.random.random(), 5) for _ in range(4)] for _ in range(3)]

    derivative = get_derivate_for_p_zeros(opponents)
    
    if all([el < 0 for el in derivative]):
        print('Found at: %s' % seed)
        break
defection_stable_opponents_set = [axl.MemoryOnePlayer(ps, initial=D) for ps in opponents] + [axl.Defector()]
tournament = axl.Tournament(defection_stable_opponents_set)
results = tournament.play(progress_bar=False)
eco = axl.Ecosystem(results)
eco.reproduce(500)
color = '#ffffff'
plot = axl.Plot(results)
p = stackplot(plot, eco)
populations = eco.population_sizes
turns = range(len(populations))
pops = [
    [populations[iturn][ir] for iturn in turns]
    for ir in plot.result_set.ranking
]
figure, ax = plt.subplots(figsize=(10, 8))

ax.stackplot(turns, *pops)
ax.yaxis.tick_left()
ax.yaxis.set_label_position("right")
ax.yaxis.labelpad = 25.0

ax.set_ylim([0.0, 1.0])
# ax.set_ylabel("Relative population size", fontweight='bold', color=color)
ax.set_xlabel("Turn", fontweight='bold', color=color)

trans = transforms.blended_transform_factory(ax.transAxes, ax.transData)

ticks = []
for i, n in enumerate(plot.result_set.ranked_names):
    x = -0.01
    y = (i + 0.5) * 1 / plot.result_set.num_players
    if n != 'Defector':
        opponent_strings = (n.split(':')[1].replace('[', '').replace(']', '')).split(",")
        opponent = [Fraction(float(op)).limit_denominator() for op in opponent_strings]
        label = '$q^{(%s)}$'% i
        n = label
    ax.annotate(
        n,
        xy=(x, y),
        xycoords=trans,
        clip_on=False,
        va="center",
        ha="right",
        fontsize=17,
        color=color
    )
    ticks.append(y)
ax.set_yticks(ticks)
ax.tick_params(direction="out")
ax.set_yticklabels([])    
ax.set_xscale("log")

ax.spines['bottom'].set_color(color)
ax.spines['top'].set_color(color) 
ax.spines['right'].set_color(color)
ax.spines['left'].set_color(color)

ax.tick_params(axis='x', colors=color)
ax.tick_params(axis='y', colors=color)
    
for tick in ax.yaxis.get_majorticklabels():  # example for xaxis
    tick.set_fontsize(20) 

plt.tight_layout()
plt.xlim(0, 10 **2);

figure.savefig('/Users/storm/src/talks/talks/2020-02-26-Max-Planck/static/population_defection_takes_over.png',
            bbox_inches='tight',
            transparent=True, dpi=100)
p.savefig('../img/population_defection_takes_over.pdf', facecolor=fig.get_facecolor(), edgecolor='none', 
          bbox_inches='tight')
for seed in range(5000):
    np.random.seed(4)
    opponents = [[round(np.random.random(), 5) for _ in range(4)] for _ in range(3)]
    
    np.random.seed(seed)
    other_opponent =[[np.random.random() for _ in range(4)]]
    
    derivative = get_derivate_for_p_zeros(opponents + other_opponent)
    if all([el < 0 for el in derivative]):
        print('Found at: %s' % seed)
        break
defection_stable_opponents_set = [axl.MemoryOnePlayer(ps, initial=C) for ps in opponents]
defection_stable_opponents_set += [axl.Defector()]
tournament = axl.Tournament(defection_stable_opponents_set)
results = tournament.play(progress_bar=False)
eco = axl.Ecosystem(results)
eco.reproduce(50000)
populations = eco.population_sizes
turns = range(len(populations))
pops = [
    [populations[iturn][ir] for iturn in turns]
    for ir in plot.result_set.ranking
]
figure, ax = plt.subplots(figsize=(10, 8))

ax.stackplot(turns, *pops)
ax.yaxis.tick_left()
ax.yaxis.set_label_position("right")
ax.yaxis.labelpad = 25.0

ax.set_ylim([0.0, 1.0])
# ax.set_ylabel("Relative population size", fontweight='bold', color=color)
ax.set_xlabel("Turn", fontweight='bold', color=color)

trans = transforms.blended_transform_factory(ax.transAxes, ax.transData)

ticks = []
for i, n in enumerate(plot.result_set.ranked_names):
    x = -0.01
    y = (i + 0.5) * 1 / plot.result_set.num_players
    if n != 'Defector':
        opponent_strings = (n.split(':')[1].replace('[', '').replace(']', '')).split(",")
        opponent = [Fraction(float(op)).limit_denominator() for op in opponent_strings]
        label = '$q^{(%s)}$'% i
        n = label
    ax.annotate(
        n,
        xy=(x, y),
        xycoords=trans,
        clip_on=False,
        va="center",
        ha="right",
        fontsize=17,
        color=color
    )
    ticks.append(y)
ax.set_yticks(ticks)
ax.tick_params(direction="out")
ax.set_yticklabels([])    
ax.set_xscale("log")

ax.spines['bottom'].set_color(color)
ax.spines['top'].set_color(color) 
ax.spines['right'].set_color(color)
ax.spines['left'].set_color(color)

ax.tick_params(axis='x', colors=color)
ax.tick_params(axis='y', colors=color)
    
for tick in ax.yaxis.get_majorticklabels():  # example for xaxis
    tick.set_fontsize(20) 

plt.tight_layout()
plt.xlim(0, 10 **2);

figure.savefig('/Users/storm/src/talks/talks/2020-02-26-Max-Planck/static/population_defection_fails.png',
            bbox_inches='tight',
            transparent=True, dpi=100)
plot = axl.Plot(results)
p = stackplot(plot, eco)

p.savefig('../img/population_defection_fails.pdf', bbox_inches='tight')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.15240234218336407, 'fraction_numerical': 0.023181198363680114, 'mean_word_length': 3.6246290801186944, 'pattern_counts': {'":': 2, '<': 27, '<?xml version=': 0, '>': 20, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '141877', 'license': 'MIT', 'n_tokens_mistral': 4920, 'n_tokens_neox': 4685, 'n_words': 853, 'path': 'nbs/2. Stability of defection.ipynb', 'repo_name': 'Nikoleta-v3/Memory-size-in-the-prisoners-dilemma', 'size': '12467'}