Decay chain simulated in Python

Question

I am trying to make a computer program which can model the decay of 1,000,000 atoms of Radon-222 (into Polonium-218 into Lead-214 into Bismuth-214 into Lead-210). Each atom has a 0.000126 chance of decaying into Polonium-218 after one unit time. I have a list of 1,000,000 0's m assigned to the variable Rn to represent the Radon-222. Then, to decide if a particular element decays, I choose a random integer between 1 and 1,000,000: if this number is less than 126, the value of Rn[j] switches from 0 to 1 (decay), and we add a 0 to the list Po for polonium, etc. So after 1 unit time, I should expect there to be around 126 Polonium atoms, after another unit time, we should have around 126 extra Polonium atoms, but some of these could have decayed into Lead-214.

We treat time as a discrete entity in this modelling. See here for a list of ideal values after 10 unit times.

Here is a code which works, but is rather slow.

import random

def one(s):
    count=0
    for i in range(len(s)):
        if s[i]==1:
            count+=1
    return count

Rn=[1]*1000000
Po=[0]*1000000
Pb214=[0]*1000000
Bi=[0]*1000000
Pb210=[0]*1000000

print([0,1000000,0,0,0,0])

for i in range(1,100):
    for j in range(len(Rn)):
        rndecay=random.random()
        if rndecay<=0.000126:
            Rn[j]=2        

    for j in range(len(Po)):
        if Po[j]==1:
            podecay=random.random()
            if podecay<=0.203:
                Po[j]=2

        if Rn[j]==2 and Po[j]==0:
            Po[j]=1

    for j in range(len(Pb214)):            
        if Pb214[j]==1:
            decay214=random.random()
            if decay214<=0.0255:
                Pb214[j]=2
        if Po[j]==2 and Pb214[j]==0:
            Pb214[j]=1

    for j in range(len(Bi)):
        if Bi[j]==1:
            bidecay=random.random()
            if bidecay<=0.0346:
                Bi[j]=2
        if Pb214[j]==2 and Bi[j]==0:
            Bi[j]=1

    for j in range(len(Pb210)):
        if Bi[j]==2 and Pb210[j]==0:
            Pb210[j]=1

    print([i,one(Rn),one(Po),one(Pb214),one(Bi),one(Pb210)])

How can I optimize this code?

Answer 1

Thanks for all your help guys! I've been able to come up with this code which uses a simple count and runs quite nicely:

import random

radon=1000000
polonium=0
lead214=0
bismuth=0
lead210=0

for t in range(1,100):
    for i in range(bismuth): 
        pdecay=random.random()
        if pdecay<0.0346:
            bismuth-=1
            lead210+=1

    for i in range(lead214): 
        pdecay=random.random() 
        if pdecay<0.0255:
            lead214-=1
            bismuth+=1

    for i in range(polonium): 
        pdecay=random.random()  
        if pdecay<0.203:
            polonium-=1
            lead214+=1

    for i in range(radon): 
        pdecay=random.random()
        if pdecay<0.000126:
            radon-=1
            polonium+=1

print([t,radon,polonium,lead214,bismuth,lead210])

Answer 2

In this particular case, the binomial law proposed by @o_o is really relevant.

If you want to keep the spirit of your code, here's several way to optimize it :

1. Use as much as possible the built-in function of python.

   one(s) is useless because list.count(value) can do the job.

2. Limit the number of loops.

   You've 5 loops for j in ... but one is enough.

3. Put your data agregation in functions.

   Defining your loops in functions will accelerate execution.

4. Use elif rather than if.

   You can avoid some useless checks with elif.

This code is about 2.5 times faster than yours :

import random

def doit(step,at):
  global Rn,Po,Pb214,Bi,Pb210
  print([0,1000000,0,0,0,0])
  for i in step:
    compute(at,Rn,Po,Pb214,Bi,Pb210)
    print([i,Rn.count(1),Po.count(1),Pb214.count(1),Bi.count(1),Pb210.count(1)])

def compute(at,Rn,Po,Pb214,Bi,Pb210):
  for j in at:
    if random.random()<=0.000126:
      Rn[j]=2

    if Po[j]==1 and random.random()<=0.203:
      Po[j]=2

    elif Rn[j]==2 and Po[j]==0:
      Po[j]=1

    if Pb214[j]==1 and random.random()<=0.0255:
      Pb214[j]=2

    elif Po[j]==2 and Pb214[j]==0:
      Pb214[j]=1

    if Bi[j]==1 and random.random()<=0.0346:
      Bi[j]=2

    elif Pb214[j]==2 and Bi[j]==0:
      Bi[j]=1

    if Bi[j]==2 and Pb210[j]==0:
      Pb210[j]=1

N = 1000000
Rn   =[1]*N
Po   =[0]*N
Pb214=[0]*N
Bi   =[0]*N
Pb210=[0]*N

step = range(1,100)
at = range(N)
doit(step,at)

My attempt to optimize with list comprehensions failed but I think it's possible.

Answer 3

I think you might have better luck approaching the problem differently. Instead of having 1,000,000 element lists, could you just keep track of the number of atoms of each species?

The number of decays in a population of n with a decay probability p follows a binomial distribution with those parameters. Many basic statistics packages will allow you to generate random numbers that follow a distribution (rather than the uniform distribution used by random.random). Thus, you only need to pull one random value for each possible transition and then update the counts with the results. This methodology is exactly equivalent to the one employed in the original post.

Here's an example using counts with transition probabilities that I made up:

from numpy.random import binomial

atoms = {'Rn222': 1000000,
         'Po218': 0,
         'Pb214': 0,
         'Bi214': 0,
         'Pb210': 0}

for _ in range(100):
    Rn222_Po218 = binomial(atoms['Rn222'], 0.000126)
    Po218_Pb214 = binomial(atoms['Po218'], 0.001240)
    Pb214_Bi214 = binomial(atoms['Pb214'], 0.003450)
    Bi214_Pb210 = binomial(atoms['Bi214'], 0.012046)

    atoms['Rn222'] -= Rn222_Po218
    atoms['Po218'] += Rn222_Po218

    atoms['Po218'] -= Po218_Pb214
    atoms['Pb214'] += Po218_Pb214

    atoms['Pb214'] -= Pb214_Bi214
    atoms['Bi214'] += Pb214_Bi214

    atoms['Bi214'] -= Bi214_Pb210
    atoms['Pb210'] += Bi214_Pb210

print atoms

Edited to add example. Edited to add binomial explanation from comment.