model.py

import pandas as pd
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import math

import seaborn as sns

def define_strats(h=10):
    return np.linspace(1, h+1, h+1, dtype=int)

def payoff_matrix(h=10, g=1, l=1):
    R = 5
    P = 3
    T = R + g*(R-P)
    S = P - l*(R-P)    

    strats = define_strats(h)

    yv, xv = np.meshgrid(strats, strats)
    fd = np.minimum(xv, yv)

    return R*(fd-1) + S*np.tril(np.ones(h+1),-1) + T*np.triu(np.ones(h+1),1) + \
            P*np.where(np.logical_and(xv == yv, xv < h+1), np.eye(h+1), np.zeros(h+1)) + P*np.maximum(h-fd,0)
            
def softmax(w):    
    maxes = np.amax(w, axis=0)    
    e = np.exp(w - maxes)
    return e / np.sum(e, axis=0)
    
# http://stackoverflow.com/a/34190035/586086
def pick_random(prob_matrix, items):
    s = prob_matrix.cumsum(axis=0)
    r = np.random.rand(prob_matrix.shape[1])
    k = (s < r).sum(axis=0)
    return items[k]
    
def simulate(payoffs, iters, n_resistant=0, nplayers=100, beta=1):
    h = payoffs.shape[0] - 1    
    strats = np.linspace(1, h+1, h+1, dtype=int)
    
    strategies = np.zeros((iters, nplayers))
    opp_strats = np.zeros((iters, nplayers))
        
    start_weights = np.zeros(h+1,)
    start_weights[h] = 1
    
    weights = np.tile(start_weights, (nplayers, 1))
    resistant = np.arange(n_resistant, dtype=int)
    
    player_idx = np.arange(nplayers)
    opp_idx = np.zeros(player_idx.shape, dtype=int)
    opp_idx[0::2] = player_idx[1::2]
    opp_idx[1::2] = player_idx[0::2]        
    
    for j in range(iters):
        # Get prob dist of other strategies
        normed_weights = weights / weights.sum(axis=1, keepdims=True)
        
        # Get strategy payoffs
        current_payoff = np.dot(payoffs, normed_weights.T)

        # Pick strategies
        strat_probs = softmax(current_payoff / beta)        
                
        picked_strats = pick_random(strat_probs, strats)
        picked_strats[resistant] = h+1
        
        strategies[j,:] = picked_strats
            
        # Assign partners to players
        players = np.arange(nplayers)
        np.random.shuffle(players)
        partners = players[opp_idx]
        
        opp_strats[j,:] = picked_strats[partners]
       
        new_weight = np.zeros(weights.shape)
        new_weight[partners, picked_strats-1] += 1
        weights = weights + new_weight
        
        # weights[partners, picked_strats-1] += 1
    
    return strategies, opp_strats