# Agent based model in Python

My model is based on the Ultimatum Game. I won't go into the intuition behind it but generally speaking it functions as follows:

1. The game consists of a n × n lattice on which an agent is placed at each node.

2. During each time step, each player on each node plays against a random neighbour by playing a particular strategy.

3. Each of their strategies (a value between 1-9) has a propensity attached to it (which is randomly assigned and is just some number). The propensity then in turn determines the probability of playing that strategy. The probability is calculated as the propensity of that strategy over the sum of the propensities of all strategies.

4. If a game results in a positive payoff, then the payoffs from that game get added to the propensities for those strategies.

5. These propensities then determine the probabilities for their strategies in the next time step, and so on.

6. The simulation ends after time step N is reached.

For games with large lattices and large time steps, my code runs really really slowly. I ran cProfiler to check where the bottleneck(s) are, and as I suspected the update_probabilitiesand play_rounds functions seem to be taking up a lot time. I want to be able to run the game with gridsize of about 40x40 for about 100000+ time steps, but right now that is not happening.

So what would be a more efficient way to calculate and update the probabilities/propensities of each player in the grid? I've considered implementing NumPy arrays but I am not sure if it would be worth the hassle here?


import numpy as np
import random
from random import randint
from numpy.random import choice
from numpy.random import multinomial
import cProfile

mew = 0.001
error = 0.05

def create_grid(row, col):
return [[0 for j in range(col)] for i in range(row)]

def create_random_propensities():
propensities = {}
pre_propensities = [random.uniform(0, 1) for i in range(9)]
a = np.sum(pre_propensities)
for i in range(1, 10):
propensities[i] = (pre_propensities[i - 1]/a) * 10 # normalize sum of propensities to 10
return propensities

class Proposer:
def __init__(self):
self.propensities = create_random_propensities()
self.probabilites = []
self.demand = 0 # the amount the proposer demands for themselves

def pick_strat(self, n_trials): # gets strategy, an integer in the interval [1, 9]
results = multinomial(n_trials, self.probabilites)
i, = np.where(results == max(results))
if len(i) > 1:
return choice(i) + 1
else:
return i + 1

def calculate_probability(self, dict_data, index, total_sum): # calculates probability for particular strat, taking propensity
return dict_data[index]/total_sum                           # of that strat as input

def calculate_sum(self, dict_data):
return sum(dict_data.values())

def initialize(self):
init_sum = self.calculate_sum(self.propensities)
for strategy in range(1, 10):
self.probabilites.append(self.calculate_probability(self.propensities, strategy, init_sum))
self.demand = self.pick_strat(1)

def update_strategy(self):
self.demand = self.pick_strat(1)

def update_probablities(self):
for i in range(9):
self.propensities[1 + i] *= 1 - mew
pensity_sum = self.calculate_sum(self.propensities)
for i in range(9):
self.probabilites[i] = self.calculate_probability(self.propensities, 1 + i, pensity_sum)

def update(self):
self.update_probablities()
self.update_strategy()

class Responder: # methods same as proposer class, can skip-over
def __init__(self):
self.propensities = create_random_propensities()
self.probabilites = []
self.max_thresh = 0 # the maximum demand they are willing to accept

def pick_strat(self, n_trials):
results = multinomial(n_trials, self.probabilites)
i, = np.where(results == max(results))
if len(i) > 1:
return choice(i) + 1
else:
return i + 1

def calculate_probability(self, dict_data, index, total_sum):
return dict_data[index]/total_sum

def calculate_sum(self, dict_data):
return sum(dict_data.values())

def initialize(self):
init_sum = self.calculate_sum(self.propensities)
for strategy in range(1, 10):
self.probabilites.append(self.calculate_probability(self.propensities, strategy, init_sum))
self.max_thresh = self.pick_strat(1)

def update_strategy(self):
self.max_thresh = self.pick_strat(1)

def update_probablities(self):
for i in range(9):
self.propensities[1 + i] *= 1 - mew # stops sum of propensites from growing without bound
pensity_sum = self.calculate_sum(self.propensities)
for i in range(9):
self.probabilites[i] = self.calculate_probability(self.propensities, 1 + i, pensity_sum)

def update(self):
self.update_probablities()
self.update_strategy()

class Agent:
def __init__(self):
self.prop_side = Proposer()
self.resp_side = Responder()
self.prop_side.initialize()
self.resp_side.initialize()

def update_all(self):
self.prop_side.update()
self.resp_side.update()

class Grid:
def __init__(self, rowsize, colsize):
self.rowsize = rowsize
self.colsize = colsize

def make_lattice(self):
return [[Agent() for j in range(self.colsize)] for i in range(self.rowsize)]

@staticmethod
def von_neumann_neighbourhood(array, row, col, wrapped=True): # gets up, bottom, left, right neighbours of some node
neighbours = set([])

if row + 1 <= len(array) - 1:

if row - 1 >= 0:

if col + 1 <= len(array) - 1:

if col - 1 >= 0:
#if wrapped is on, conditions for out of bound points
if row - 1 < 0 and wrapped == True:

if col - 1 < 0 and wrapped == True:

if row + 1 > len(array) - 1 and wrapped == True:

if col + 1 > len(array) - 1 and wrapped == True:
return neighbours

def get_error_term(pay, strategy):
index_strat_2, index_strat_8 = 2, 8
if strategy == 1:
return (1 - (error/2)) * pay, error/2 * pay, index_strat_2
if strategy == 9:
return (1 - (error/2)) * pay, error/2 * pay, index_strat_8
else:
return (1 - error) * pay, error/2 * pay, 0

class Games:
def __init__(self, n_rows, n_cols, n_rounds):
self.rounds = n_rounds
self.rows = n_rows
self.cols = n_cols
self.lattice = Grid(self.rows, self.cols).make_lattice()
self.lookup_table = np.full((self.rows, self.cols), False, dtype=bool)  # if player on grid has updated their strat, set to True

def reset_look_tab(self):
self.lookup_table = np.full((self.rows, self.cols), False, dtype=bool)

def run_game(self):
n = 0
while n < self.rounds:
for r in range(self.rows):
for c in range(self.cols):
if n != 0:
self.lattice[r][c].update_all()
self.lookup_table[r][c] = True
self.play_rounds(self.lattice, r, c)
self.reset_look_tab()
n += 1

def play_rounds(self, grid, row, col):
neighbours = Grid.von_neumann_neighbourhood(grid, row, col)
neighbour = random.sample(neighbours, 1).pop()
neighbour_index = [(ix, iy) for ix, row in enumerate(self.lattice) for iy, i in enumerate(row) if i == neighbour]
if self.lookup_table[neighbour_index][neighbour_index] == False: # see if neighbour has already updated their strat
neighbour.update_all()
player = grid[row][col]
coin_toss = randint(0, 1) # which player acts as proposer or responder in game
if coin_toss == 1:
if player.prop_side.demand <= neighbour.resp_side.max_thresh: # postive payoff
payoff, adjacent_payoff, index = get_error_term(player.prop_side.demand, player.prop_side.demand)
if player.prop_side.demand == 1 or player.prop_side.demand == 9: # extreme strategies get bonus payoffs
player.prop_side.propensities[player.prop_side.demand] += payoff
else:
player.prop_side.propensities[player.prop_side.demand] += payoff
else:
return 0 # if demand > max thresh -> both get zero

if coin_toss != 1:
if neighbour.prop_side.demand <= player.resp_side.max_thresh:
payoff, adjacent_payoff, index = get_error_term(10 - neighbour.prop_side.demand, player.resp_side.max_thresh)
if player.resp_side.max_thresh == 1 or player.resp_side.max_thresh == 9:
player.resp_side.propensities[player.resp_side.max_thresh] += payoff
else:
player.resp_side.propensities[player.resp_side.max_thresh] += payoff
else:
return 0

# pr = cProfile.Profile()
# pr.enable()

my_game = Games(10, 10, 2000) # (rowsize, colsize, n_steps)
my_game.run_game()

# pr.disable()
# pr.print_stats(sort='time')



(For those who might be wondering, the get_error_term just returns the propensities for strategies that are next to strategies that receive a positive payoff, for example if the strategy 8 works, then 7 and 9's propensities also get adjusted upwards and this is calculated by said function. And the first for loop inside update_probabilities just makes sure that the sum of propensities don't grow without bound).

Proposer and Responder are basically identical. The only difference is a single variable name. So get rid of one of them. This reduces code duplication and makes refactoring easier.

# 1. Profiling

It's always a good idea to analyze thoroughly what it is that makes your code slow. Everything else is premature optimization. In the spirit of "Measure twice, cut once", enter cProfile for (results for Games(10, 10, 2000)):

         19137129 function calls (19134785 primitive calls) in 43.470 seconds

Ordered by: cumulative time

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
146/1    0.002    0.000   43.471   43.471 {built-in method builtins.exec}
1    0.000    0.000   43.470   43.470 cr-222974_pre.py:1(<module>)
1    0.801    0.801   43.192   43.192 cr-222974_pre.py:130(run_game)
299740    0.385    0.000   28.192    0.000 cr-222974_pre.py:69(update_all)
599480    0.597    0.000   27.807    0.000 cr-222974_pre.py:58(update)
200000    3.809    0.000   23.708    0.000 cr-222974_pre.py:142(play_rounds)
599480    0.662    0.000   15.075    0.000 cr-222974_pre.py:48(update_strategy)
599680    5.023    0.000   14.418    0.000 cr-222974_pre.py:28(pick_strat)
599480    8.525    0.000   12.134    0.000 cr-222974_pre.py:51(update_probablities)
599680    5.486    0.000    5.486    0.000 {method 'multinomial' of 'mtrand.RandomState' objects}
200000    1.142    0.000    3.022    0.000 random.py:286(sample)
599686    2.936    0.000    2.936    0.000 {built-in method builtins.max}
174021    2.565    0.000    2.565    0.000 cr-222974_pre.py:110(get_error_term)
200000    2.511    0.000    2.511    0.000 cr-222974_pre.py:145(<listcomp>)
599680    0.508    0.000    2.017    0.000 cr-222974_pre.py:39(calculate_sum)
5397120    1.593    0.000    1.593    0.000 cr-222974_pre.py:36(calculate_probability)
599680    1.355    0.000    1.355    0.000 {built-in method builtins.sum}
403609    0.342    0.000    1.300    0.000 {built-in method builtins.isinstance}
200000    0.974    0.000    1.169    0.000 cr-222974_pre.py:81(von_neumann_neighbourhood)
400000    0.665    0.000    1.055    0.000 random.py:224(_randbelow)
200000    0.168    0.000    1.025    0.000 random.py:218(randint)
400000    0.149    0.000    0.957    0.000 abc.py:137(__instancecheck__)
200000    0.354    0.000    0.857    0.000 random.py:174(randrange)
599680    0.845    0.000    0.845    0.000 {built-in method numpy.where}
...


First of all, that sets our baseline to around $$\43 s\$$ (to be precise: I actually measured several times and all the results where around that time). Also looking at the results IMHO reveals the following "hot spots":

• update_probabilies takes about a quarter of the time (almost $$\12 s\$$) where $$\8 s\$$ are spent in the method itself, which leave about $$\4 s\$$ to be spent in calculate_sum and calculate_probability. That matches the cumtime values for those two methods as measured by the compiler.
• pick_strat accounts for another $$\14 s\$$ of the total time. More than a third of this is spent on numpy.random.multinomial. This is quite a lot of time for picking a (few) random value(s) in such a small range.

# 2. Optimization

Measuring: Done! Not comes the cut...

## Picking a strategy

You will have to verify this, but I think the strategy selection can be simplified heavily. If I understood your code correctly, your basically choosing random values between 1 and 9 following a given probability distribution. This can be done with the following piece of code:

def pick_strat(self, n_trials):
chosen = choices((1, 2, 3, 4, 5, 6, 7, 8, 9), self.probabilites, k=n_trials)
if n_trials > 1:
return chosen
return chosen


This simple change brings the overall execution time down to about $$\28 s\$$ here on my machine.

For some reason you picked a dict to store your propensity values, which is a bit unusual in that context, especially considering that the probabilities are stored in a plain list. A dict with consecutive integers as key screams to be an list or something similar. Further looking at the kinds of operation that are be performed on these values supports this assumption even further. For update_probablities this means we go from

def update_probablities(self):
for i in range(9):
self.propensities[1 + i] *= 1 - mew
pensity_sum = self.calculate_sum(self.propensities)
for i in range(9):
self.probabilites[i] = self.calculate_probability(self.propensities, 1 + i, pensity_sum)


to this

mew_complement = 1 - mew

# ... lot of code here

def update_probablities(self):
for i in range(9):
self.propensities[i] *= mew_complement
pensity_sum = sum(self.propensities)
for i in range(9):
self.probabilites[i] = self.propensities[i] / pensity_sum


As you can see, I aggressively removed function calls to improve the performance even further. IMHO the code is still quite readable though. Since I was at it, I also decided to give create_random_propensities and initialize a makeover:

def create_random_propensities(self):
pre_propensities = [random.uniform(0, 1) for i in range(9)]
pensity_sum = sum(pre_propensities) / 10
return [pre / pensity_sum for pre in pre_propensities]

def initialize(self):
init_sum = sum(self.propensities)
self.probabilites = [prop / init_sum for prop in self.propensities]
self.update_strategy()


(Just by the way: I'm not entirely sure what the purpose of * 10 is.) As you can see, I allowed myself to move create_random_propensities into the class as well. Keep in mind that the strategy picking and reward parts later on will have to be slightly adapted to account for the 0-based indexing. With all of these changes in place the runtime is now down to about $$\20 s\$$.

I also decided to implement these changes using NumPy, but at least for Games(10, 10, 2000) using numpy.ndarray instead of list and array operations instead of explicit loops performed considerably worse.

## The game

Now the two mentioned hotspots are handled and allowed the performance to double for the chosen setting. Time to look at the rest of the code.

The first thing that caught my eye was

neighbour = random.sample(neighbours, 1).pop()
neighbour_index = [(ix, iy) for ix, row in enumerate(self.lattice) for iy, i in enumerate(row) if i == neighbour]
if self.lookup_table[neighbour_index][neighbour_index] ==
False: # see if neighbour has already updated their strat
neighbour.update_all()


This seems like a really wasteful approach to check if update_all has already been called on an agent, especially if you take into account that the list comprehension will have to look at all agents even if the correct position has already been found. It would be much more straightforward to have a bool variable in each agent that stores if the update was done or not. This variable could even be set automatically when calling update_all() on an agent.

class Agent:

def __init__(self):
# Participant is the name I chose when I collapsed Proposer and Responder into a single class
self.prop_side = Participant()
self.resp_side = Participant()
self.prop_side.initialize()
self.resp_side.initialize()
self.has_updated_strategy = False

def update_all(self):
self.has_updated_strategy = True
self.prop_side.update()
self.resp_side.update()


With these change the part from above becomes as simple as

if not neighbour.has_updated_strategy:
neighbour.update_all()


The changes to run_game and reset_look_tab should be straightforward. Oh, and while you are at it, the outermost loop in run_game should also be a for loop as well. With all of the changes from above in place we are now at a runtime of about $$\14 s\$$. Not bad, I would say.

# 3. Other considerations

The main part of play_rounds also presents a severe case of code duplication, especially if you follow my advice from above and merge Proposer and Responder into a single class, leaving you with just one name for demand/max_tresh. Then the only difference between the two large code blocks is which agent assumes which role and how the reward is being computed. Sidenote/Bug: In the first block there is payoff, adjacent_payoff, index = get_error_term(player.prop_side.demand, player.prop_side.demand) which is likely a bug.

The code could heavily profit from some documentation. In essence, it's good practice to desribe your methods with a short """documation string""". So instead of

def pick_strat(self, n_trials): # gets strategy, an integer in the interval [1, 9]



do

def pick_strat(self, n_trials):
"""gets strategy, an integer in the interval [1, 9]"""