Monte Carlo simulation of amoeba population

Question

I wrote a simple Python simulation to answer the "Amoeba" population question posed here:

A population of amoebas starts with 1. After 1 period that amoeba can divide into 1, 2, 3, or 0 (it can die) with equal probability. What is the probability that the entire population dies out eventually?

The script seems to work fine, but then I started playing with it to investigate curious effects.

If the number of trial generations is too large (num_generations_per_trial), I have problems with performance - the population size gets huge, and the simulation either runs slow or I encounter OverflowError on my brute force FOR loops.

I would appreciate feedback on efficiency options, and also on general code improvements. I know the runs are independent and could be run in parallel. But that is still sort of brute force. I am curious more about making a single thread approach faster.

from __future__ import division
import random
import math
import time


def run_trial(max_split, num_generations):
    population = 1
    for generation in xrange(num_generations):
        for amoeba in xrange(population):
            amoeba_split = random.randint(0, max_split)
            population -= 1  # remove current amoeba (she will split or die)
            population += amoeba_split

        if population == 0:
            break
    return population


def main():
    extinct_counter = 0
    trials = 10000

    max_split_per_amoeba = 3
    num_generations_per_trial = 20 # populations can get *massive* as generations increase (memory / overflow errors at 100)

    print '***starting simulation***'
    print 'num trials: %i' % (trials)
    print 'max_split_per_amoeba: %i' % (max_split_per_amoeba)
    print 'num_generations_per_trial: %i' % (num_generations_per_trial)

    for trial in xrange(trials):
        outcome_population = run_trial(max_split_per_amoeba, num_generations_per_trial)
        if outcome_population == 0:
            extinct_counter += 1

        if divmod(trial+1, max(1,int(trials/20)))[1] == 0:
            print 'progress: %i trials complete | %i extinction counter | %.4f extinction probability' %  (trial+1, extinct_counter, extinct_counter/(trial+1))

    print 'extinct outcomes: %i' % (extinct_counter)
    print 'total trials: %i' % (trials)

    extinction_probability = extinct_counter / trials
    print 'extinction probability: %.4f' % (extinction_probability)

    expected_answer = math.sqrt(2)- 1
    print 'expected probability: %.4f' % (expected_answer)
    print 'delta from answer: %.4f' % (extinction_probability - expected_answer)


if __name__ == '__main__':
    start = time.clock()
    main()
    print 'runtime: %.3f s' % (time.clock() - start)
    print 'done'

Answer 1

Performance

Limiting each trial to 20 generations is not really helpful. Let's put it this way: the population is expected to increase exponentially with each generation. Once the population reaches, say, 10 amoebas, it's pretty safe to declare that infection has taken hold and the population will never dwindle to 0.

If you run a trial for 20 generations, you can easily get populations in the thousands. The problem with letting populations grow so large is that you have to decide the fate of each individual amoeba at each time step, and that ends up being a lot of work.

If you change num_generations_per_trial to 10, you'll get identical results, and finish much more quickly. But how do you decide what an appropriately low num_generations_per_trial number is? It would be much more intuitive to impose a cap on the population size than on the number of generations.

Expressiveness

I suggest writing a generator for the population count at each time step. You can easily count elements using the sum() function.

from __future__ import division
from random import randint

def amoeba_generations(max_split=3, population=1):
    yield population
    while population > 0:
        population = sum(randint(0, max_split) for _ in range(population))
        yield population

def trial(generator, success_pred, failure_pred):
    for g in generator:
        if success_pred(g):
            return True
        if failure_pred(g):
            return False
    return False


def main():
    num_trials = 10000
    max_split_per_amoeba = 3
    num_dead_populations = sum(
        trial(
            amoeba_generations(max_split_per_amoeba),
            lambda n: n == 0,
            lambda n: n >= 10     # Assume a population of >= 10 will never die
        )
        for _ in range(num_trials)
    )
    print("""***starting simulation***
num trials: {num_trials}
max_split_per_amoeba: {max_split_per_amoeba}
extinct outcomes: {num_dead_populations}
extinction probability: {extinction_probability:.4f}
expected probability: {expected_probability:.4f}
delta from answer: {delta:.4f}""".format(
        extinction_probability=num_dead_populations/num_trials,
        expected_probability=2**.5 - 1,
        delta=(num_dead_populations/num_trials) - (2**.5 - 1),
        **locals()
    ))

if __name__ == '__main__':
    main()

Answer 2

• with performance problems always profiler is huge help. Sometimes your IDE support it, sometimes you must use it directly

• i replace your __main__ section by code

  if __name__ == "__main__":
      import cProfile # should be at the top of file
      cProfile.run('main()')
  

• when run, i get interesting result

           198906976 function calls in 99.713 seconds
  
     Ordered by: standard name
  
     ncalls  tottime  percall  cumtime  percall filename:lineno(function)
          1    0.000    0.000   99.713   99.713 <string>:1(<module>)
          1    0.011    0.011   99.713   99.713 ameba.py:20(main)
      10000   18.129    0.002   99.697    0.010 ameba.py:7(run_trial)
   66292324   59.878    0.000   63.150    0.000 random.py:175(randrange)
   66292324   18.417    0.000   81.567    0.000 random.py:238(randint)
      10000    0.002    0.000    0.002    0.000 {divmod}
          1    0.000    0.000    0.000    0.000 {math.sqrt}
      10000    0.002    0.000    0.002    0.000 {max}
          1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
   66292324    3.272    0.000    3.272    0.000 {method 'random' of '_random.Random' objects}
  

• this result means that your program 81.5s of 99s total time is inside random.randint(0, max_split) function. This is important information. You should now only optimize this line, nothing more.

• first though: Can i reduce calls of randint()? Unfortunately, without changing brute-force algorithm it's impossible.

• second though: What does to hell with randrange? I don't use it in my code. The answer: randint calls randrange. So, maybe we call it directly. I replace random.randint(0, max_split) by random.randrange(0, max_split + 1). Result:

           132657778 function calls in 69.706 seconds
  
     Ordered by: standard name
  
     ncalls  tottime  percall  cumtime  percall filename:lineno(function)
          1    0.000    0.000   69.706   69.706 <string>:1(<module>)
          1    0.010    0.010   69.706   69.706 ameba.py:20(main)
      10000   16.422    0.002   69.691    0.007 ameba.py:7(run_trial)
   66313887   50.472    0.000   53.269    0.000 random.py:175(randrange)
      10000    0.002    0.000    0.002    0.000 {divmod}
          1    0.000    0.000    0.000    0.000 {math.sqrt}
      10000    0.002    0.000    0.002    0.000 {max}
          1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
   66313887    2.797    0.000    2.797    0.000 {method 'random' of '_random.Random' objects}
  

• it's better, but still this one line takes 53s of 70s total time. We are still need only optimize this line. When i try google this problem, i get similar answers: use numpy. It's not a pure standard library, so it's your choice.

• from my point of view numpy and scipy are "standard libraries" :)