# Monte Carlo Simulation of "DNA Mutation" Game Mechanic

The game Space Station 13 (on the Paradise Station codebase) contains a DNA mechanic. Every carbon-based lifeform with DNA has 55 "blocks" that can be altered in a DNA Modifier, giving the mob disabilities or superpowers once they reach a certain threshold.

DNA blocks are 3-digit base-16 numbers from 000 to FFF. The activation thresholds are different for each kind of mutation (disabilities and minor powers have lower thresholds), but if a block is set to a value of hexadecimal DAC or above, all of the mutations on that block are guaranteed to be active.

Here is an example of the DNA modification process:

1. A creature enters the DNA modifier.
2. The Geneticist looks at the creature's first DNA block. It is set to 357.
3. The Geneticist selects the first digit of the block, 3, and irradiates it.
4. The digit 3 is replaced with the digit 9.
5. 957 is below the threshold of DAC. The Geneticist irradiates the first digit again, until it is at least D.
6. If the first digit is above D, the process ends. If it is at D, the geneticist irradiates that block until it is at A or higher. Then, the third block, until it is at C or higher.
7. The process continues to until a value of DAC or above is reached. (If the first irradiation had resulted in E57, that would have ended it early.)

I would like to find out how many attempts, on average, it takes to irradiate a block to DAC. Not every digit can turn into every other digit, and the probability table is not uniform. Fortunately, I have access to the source code of Paradise Station, and thus I can accurately simulate the process.

The game is written in a rather clunky language called DreamMaker - I don't really know how to work with it. The function responsible for irradiating a single is called miniscramble, its source code is below: (modified slightly to be more readable, the defines were not originally there)

#define HIGH_SCRAMBLE prob((rs*10))
#define MED_SCRAMBLE  prob((rs*10)-(rd))
#define LOW_SCRAMBLE  prob((rs*5)+(rd)))
/proc/miniscramble(input,rs,rd)
var/output
output = null
if(input == "C" || input == "D" || input == "E" || input == "F")
output = pick(HIGH_SCRAMBLE;"4",HIGH_SCRAMBLE;"5",HIGH_SCRAMBLE;"6",HIGH_SCRAMBLE;"7",LOW_SCRAMBLE;"0",LOW_SCRAMBLE;"1",MED_SCRAMBLE;"2",MED_SCRAMBLE;"3")
if(input == "8" || input == "9" || input == "A" || input == "B")
output = pick(HIGH_SCRAMBLE;"4",HIGH_SCRAMBLE;"5",HIGH_SCRAMBLE;"A",HIGH_SCRAMBLE;"B",LOW_SCRAMBLE;"C",LOW_SCRAMBLE;"D",LOW_SCRAMBLE;"2",LOW_SCRAMBLE;"3")
if(input == "4" || input == "5" || input == "6" || input == "7")
output = pick(HIGH_SCRAMBLE;"4",HIGH_SCRAMBLE;"5",HIGH_SCRAMBLE;"A",HIGH_SCRAMBLE;"B",LOW_SCRAMBLE;"C",LOW_SCRAMBLE;"D",LOW_SCRAMBLE;"2",LOW_SCRAMBLE;"3")
if(input == "0" || input == "1" || input == "2" || input == "3")
output = pick(HIGH_SCRAMBLE;"8",HIGH_SCRAMBLE;"9",HIGH_SCRAMBLE;"A",HIGH_SCRAMBLE;"B",MED_SCRAMBLE;"C",MED_SCRAMBLE;"D",LOW_SCRAMBLE;"E",LOW_SCRAMBLE;"F")
if(!output) output = "5"
return output


My program simulates trying to get each block from 000-800 (values above 800 do not occur in creatures at the start of a round) to at least DAC, by using the strategy outlined above. It is written in Python 3 and outputs a CSV file with two columns: the block being worked on, and how long it takes to get it to DAC, averaging 4096 attempts. The miniscramble_table function is my direct translation of the DreamMaker code above.

from random import random
from functools import reduce
import itertools

print("This program will calculate how many attempts it takes to get to DAC from every block 000-800, on average.")
block_attempts = input("How many attempts should be made per block? (default 4096, more is more accurate)")

if not block_attempts:
block_attempts = 4096
else:
try:
block_attempts = int(block_attempts)
except ValueError:
print("ERR: Could not parse that number. Using 4096.")
user_bound = input("What should be the upper bound on the dataset? (default 800, lower is faster but less comprehensive)")
if not user_bound:
user_bound = '800'
else:
try:
user_bound = [char for char in user_bound]
assert(len(user_bound) == 3)
except:
print("ERR: Could not parse that bound. Using 800.")

user_target = input("What should be the target value for a block? (default DAC)")
if not user_target:
user_target = ['D','A','C']
else:
try:
user_target = [char for char in user_target]
assert(len(user_target) == 3)
except: # bad, but like, c'mon
print("ERR: Could not parse that bound. Using 800.")

# Generate a probability table. This is ugly because it's based off BYOND code.

picks = (("5"),(1.0)) # default, I guess.
if (letter in ["C",  "D", "E", "F"]):
picks = ("4", "5", "6", "7", "0", "1", "2", "3")
probs = (*[HIGH_SCRAMBLE] * 4, *[LOW_SCRAMBLE] * 2, *[MED_SCRAMBLE] * 2)

if (letter in ["8",  "9", "A", "B", "4", "5", "6", "7"]):
picks = ("4", "5", "A", "B", "C", "D", "2", "3")
probs = (*[HIGH_SCRAMBLE] * 4, *[LOW_SCRAMBLE] * 4)

if (letter in ["0",  "1", "2", "3"]):
picks = ("8", "9", "A", "B", "C", "D", "E", "F")
probs = (*[HIGH_SCRAMBLE] * 4, *[MED_SCRAMBLE] * 2, *[LOW_SCRAMBLE] * 2)

total = sum(probs)
probs = map(lambda n: n/total, probs) # sums to 1

# make the output nicer to work with...
out = []
prev = 0
for pick, prob in zip(picks, probs):
out.append((pick, prob+prev))
prev += prob
return out

r = random()
output = filter(lambda entry: entry[1] >= r, table)
# print(r)
return list(output)[0][0]

# tries to get from initial to at least letters with specified settings
# returns # of attempts to get there.
def scramble_to(initial=('3','5','7'), target=user_target, settings=(10,2), log=False):
current = list(initial) # what are we looking at
# letter-iterable to base10 number
def concat_letters(letters):
return int(reduce(lambda x,y: x+y, letters, ''), 16)

for attempts in enumerate(itertools.repeat(0)):
if log: print(f'Miniscramble #{attempts[0]}:', ''.join(current))
if concat_letters(current) >= concat_letters(target):
if log: print(f'Done with {attempts[0]} miniscrambles!')
return attempts[0] # done, since we're above/at the target!
for i in range(3):
if int(current[i], 16) < int(target[i], 16):
current[i] = miniscramble(current[i], *settings)
break # 1 enumerate per attempt

results = {}
def unactivated(seq):
return int(''.join(seq), 16) < int(user_bound, 16) # blocks never start activated, so default is 800

dataset = filter(unactivated, (seq for seq in itertools.product([_ for _ in '0123456789ABCDEF'], repeat=3)))
for block in dataset: # go through a sample set of blocks, 54 in all
# Give each block lots of attempts for bigger sample size. default=4096
intermediate = []
for _ in range(block_attempts):
intermediate.append(scramble_to(initial=block, target=('D','A','C'), settings=(1,2)))
results[block] = (sum(intermediate)/len(intermediate)) # average it out

# Convert results to CSV

out = []
for k,v in results.items():
out.append(f'{"".join(k)},{v}\n')
filename = input('\nDone. Where should the results be saved? (leave blank to not save) ')
if filename:
with open(filename, 'w') as outfile:
[outfile.write(line) for line in out]
else:
print("Not saving.")


I made this plot of the results using Excel:

### Questions:

• What could I do better?
• Are there libraries that make this kind of simulation easier to write, or faster?
• Generating the output takes a pretty long time right now: are there any obvious optimizations to be made? Would it be beneficial to make this asynchronous, spawn worker threads somehow, compile in Cython?
• How could I generate a similar graph using a library like matplotlib?
• What's your question? What do you want from the review? Apr 7, 2020 at 3:22
• @RootTwo Oh yeah, the question part of asking questions... whoops. Added a section. Thanks! Apr 7, 2020 at 18:32

This problem is an "absorbing Markov chain", and the expected number of steps can be solved analytically.

The Markov chain has a node or state corresponding to each of the DNA blocks. The miniscramble routine, along with the steps of the DNA modification process, can be used to define the transition probablities between states. For example, 0x000 can transistion to 0x100, 0x200, 0x300, ... (only the first digit changes). Similarly 0xD05 can go to 0xD15 ...0xDF5 (only the second digit changes) and so on. Any node >= 0xDAC is an absorbing node.

The code could be cleaner, but it demonstrates the point.

import numpy as np
import matplotlib.pyplot as plt

# Hi, Med, Lo transition weights

transition_probability = []

# for digits 0, 1, 2, 3
# picks    0  1  2  3  4  5  6  7  8  9  A  B  C  D  E  F
weights = [0, 0, 0, 0, 0, 0, 0, 0, H, H, H, H, M, M, L, L]
total = sum(weights)
probabilities = [w/total for w in weights]
transition_probability.extend(probabilities for _ in '0123')

# for digits 4, 5, 6, 7, 8, 9, A, B
# picks    0  1  2  3  4  5  6  7  8  9  A  B  C  D  E  F
weights = [0, 0, L, L, H, H, 0, 0, 0, 0, H, H, L, L, 0, 0]
total = sum(weights)
probabilities = [w/total for w in weights]
transition_probability.extend(probabilities for _ in '456789AB')

# for digits C, D, E, F:
#picks     0  1  2  3  4  5  6  7  8  9  A  B  C  D  E  F
weights = [L, L, M, M, H, H, H, H, 0, 0, 0, 0, 0, 0, 0, 0]
total = sum(weights)
probabilities = [w/total for w in weights]
transition_probability.extend(probabilities for _ in 'CDEF')

return transition_probability

# build table of all transitions
# P[i][j] = prob to go from i to j
P = []

for i in range(0xFFF + 1):
d0, rem = divmod(i, 0xFF)
d1, d2 = divmod(rem, 0xF)

row = [0]*4096

if d0 < 0xD:
start = d1*0xF + d2
for c, j in enumerate(range(start, start + 0xF00 + 1, 0x100)):
row[j] = transition[d0][c]

elif d0 == 0xD:
if d1 < 0xA:
start = d0 * 0xFF + d2
for c, j in enumerate(range(start, start + 0xF0 + 1, 0x10)):
row[j] = transition[d1][c]

elif d1 == 0xA:
if d2 < 0xC:
start = d0 * 0xFF + d1 * 0xF
for c, j in enumerate(range(start, start + 0xF + 1, 0x1)):
row[j] = transition[d2][c]

P.append(row)

# convert to numpy array to do to more easily
# select Q and do the matrix math
P = np.array(P)

Q = P[:0xDAB,:0xDAB]

I = np.identity(Q.shape[0])

N = np.linalg.inv(I - Q)

# this is the same a N*1 as shown in the Wikipedia article
avg_steps = np.sum(N, axis=1)

# change indices for avg_steps to view different
# ranges of starting points
plt.plot(avg_steps[:0x801])


The table of probabilities is set as soon as your rad strength and duration are known. You only need to generate it once, for all letters, and then use it as a lookup table.

from collections import defaultdict

picks = defaultdict(lambda: "5")
probs = defaultdict(lambda: 1.)
for letter in ["C",  "D", "E", "F"]:
# picks[letter] = ("4", "5", "6", "7", "0", "1", "2", "3")
# probs[letter] = (*[high] * 4, *[low] * 2, *[medium] * 2)
picks[letter] = ("0", "1", "2", "3", "4", "5", "6", "7")
probs[letter] = (low, low, medium, medium, high, high, high, high)
for letter in ["8",  "9", "A", "B", "4", "5", "6", "7"]:
# picks[letter] = ("4", "5", "A", "B", "C", "D", "2", "3")
# probs[letter] = (*[high] * 4, *[low] * 4)
picks[letter] = ("2", "3", "4", "5", "A", "B", "C", "D")
probs[letter] = (low, low, high, high, high, high, low, low)
for letter in ["0",  "1", "2", "3"]:
picks[letter] = ("8", "9", "A", "B", "C", "D", "E", "F")
# probs[letter] = (*[high] * 4, *[medium] * 2, *[low] * 2)
probs[letter] = (high, high, high, high, medium, medium, low, low)
return picks, probs


Note that I did not normalize the weights. This is because you can use random.choices, which takes a weight argument and normalizes it for you if necessary:

from random import choices

def scramble(letter, picks, probs):
return choices(picks[letter], probs[letter])[0]


And if that seems like not enough for a function anymore, you might be right. You also don't need your concat_letters function. Lists are directly comparable and so are strings (which are larger than the strings of the numbers and sorted lexicographically): ["D", "A", "C"] > ["9", "A", "3"] -> True. Instead I made the scramble function an inner function, so you don't have to pass the picks and probabilities every time:

from itertools import count
from random import choices

def scramble_to(initial, target, picks, probs, log=False):
def scramble(letter):
return choices(picks[letter], probs[letter])[0]

current = initial.copy()
for attempt in count():
if log:
print(f'Miniscramble #{attempt}:', ''.join(current))
if current >= target:
if log:
print(f'Done with {attempt} miniscrambles!')
return attempt
for i in range(3):
if current[i] < target[i]:
current[i] = scramble(current[i])
break   # only scramble one letter per attempt


The ouput of itertools.product is directly iterable, and so are strings. No need to iterate over them in a needless list/generator comprehension, respectively. You can also use statistics.mean instead of doing it yourself.

from statistics import mean
from itertools import product

def inactive(seq):
return seq < ("8", "0", "0") # blocks never start activated, so default is 800

if __name__ == "__main__":
picks, probs = get_probabilities(1, 2)
block_attempts = 4096
target = ["D", "A", "C"]
results = {}
dataset = filter(inactive, product('0123456789ABCDEF', repeat=3))
results = {block: mean(scramble_to(list(block), target, picks, probs)
for _ in range(block_attempts))
for block in dataset}


(Untested for now)

Note that I used a if __name__ == "__main__": guard to ensure this code is not being run when importing from this script and followed Python's official style-guide, PEP8, throughout this answer. PEP8 recommends using 4 spaces as indentation and always using a newline after e.g. an if log:.

As for visualization, here is a quick attempt at replicating your graph:

import matplotlib.pyplot as plt

plt.style.use("dark_background")
plt.figure(figsize=(12, 6))
plt.plot(list(results.values()), 'o-', c="orange")
plt.grid()
plt.title("How many times do you have to press the \"Irradiate Block\" button to get a block to DAC?")
plt.xlabel("Block #")
plt.ylabel("# of attempts")
plt.show()


which results in