Background
Using this simulation I investigate a system in which enzymes proliferate in cells. During the replications of enzymes, parasites can come to be due to mutation. They can drive the system into extinction. I'm interested in where in the parameter space coexistence is possible.
In the program the system is a list, the cells are dictionaries with 2 keys: "e" for the enzymes and "p" for the parasites. The values of the keys are the numbers of the 2 variants.
Our parameters are:
pop_size: the number of the cells
cell_size: the maximal number of molecules (enzymes+parasites) of cells at which cell division takes place
a_p: fitness of the parasites relative to the fitness of the enzymes (for example if a_p = 2, the parasites' fitness is twice as that of the enzymes)
mutation_rate: the probability of mutation during a replication event
gen_max: the maximal number of generations (a generation corresponds to one while cycle; if the system extincts, the program doesn't run until gen_max)
We start with pop_size cells with cell_size // 2 enzimes and 0 parasites. In each cell the molecules proliferate until their number reaches cell_size. Each cell divides, the assortment of the molecules happens according to binomial distributions (\$p=0.5\$). Cells with "e" < 2 are discarded as dead. After that if the number of viable cells is bigger than pop_size, we choose pop_size of them according to cell fitness ("e"/("e"+"p")), and they move on to the next generation. On the other hand, if the number of viable cells is pop_size or less, they all move on to the next generation.
My request
I've never studied programming in school. This program is the result of heavy googling. Now I've reached a point where I need advice from experienced people. At certain parameter values the program gets quite slow.
So, what better solutions exist performance-wise than my solutions for the manipulations of the list's items throughout the program and for writing data to file? And algorithm design-wise?
In which directions should I improve my programming skills in Python to efficiently implement these kind of models?
Or am I near the limit of Python's capabilities in this regard? So should I change to a more appropriate programming language in order to achieve significantly better performance at these kind of tasks? If yes, which languages should I consider? (My guess is C.)
Of course any advice is highly appreciated!
The code
The program consists of two functions. simulation() does the simulation, writeoutfile() writes the data to file.
# -*- coding: utf-8 -*-
from random import choices, random
import csv
import time
import numpy as np
def simulation(pop_size, cell_size, a_p, mutation_rate, gen_max):
def fitness(pop):
return [i["e"] / (i["e"] + i["p"]) for i in pop]
def output(pop, gen, pop_size, cell_size, mutation_rate, a_p, boa_split):
if pop:
gyaklist_e = [i["e"] for i in pop]
gyaklist_p = [i["p"] for i in pop]
fitnesslist = fitness(pop)
return (
gen,
sum(gyaklist_e), sum(gyaklist_p),
sum([1 for i in pop if i["e"] > 1]),
np.mean(gyaklist_e), np.var(gyaklist_e),
np.percentile(gyaklist_e, 25),
np.percentile(gyaklist_e, 50),
np.percentile(gyaklist_e, 75),
np.mean(gyaklist_p), np.var(gyaklist_p),
np.percentile(gyaklist_p, 25),
np.percentile(gyaklist_p, 50),
np.percentile(gyaklist_p, 75),
np.mean(fitnesslist), np.var(fitnesslist),
np.percentile(fitnesslist, 25),
np.percentile(fitnesslist, 50),
np.percentile(fitnesslist, 75),
pop_size, cell_size, mutation_rate, a_p, boa_split
)
return (
gen,
0, 0,
0,
0, 0,
0, 0, 0,
0, 0,
0, 0, 0,
0, 0,
0, 0, 0,
pop_size, cell_size, mutation_rate, a_p, boa_split
)
pop = [{"e": cell_size // 2, "p": 0} for _ in range(pop_size)]
gen = 0
yield output(
pop,
gen, pop_size, cell_size, mutation_rate, a_p, boa_split="aft"
)
print(
"N = {}, rMax = {}, aP = {}, U = {}".format(
pop_size, cell_size, a_p, mutation_rate
)
)
while pop and gen < gen_max:
gen += 1
for i in pop:
while not i["e"] + i["p"] == cell_size:
luckyreplicator = choices(
["e", "p"], [i["e"], a_p*i["p"]]
)
if luckyreplicator[0] == "e" and random() < mutation_rate:
luckyreplicator[0] = "p"
i[luckyreplicator[0]] += 1
if gen % 100 == 0:
yield output(
pop,
gen, pop_size, cell_size, mutation_rate, a_p, boa_split="bef"
)
newpop = [
{"e": np.random.binomial(i["e"], 0.5),
"p": np.random.binomial(i["p"], 0.5)}
for i in pop
]
for i in zip(pop, newpop):
i[0]["e"] -= i[1]["e"]
i[0]["p"] -= i[1]["p"]
pop += newpop
newpop = [i for i in pop if i["e"] > 1]
if newpop:
fitnesslist = fitness(newpop)
fitness_sum = np.sum(fitnesslist)
fitnesslist = fitnesslist / fitness_sum
pop = np.random.choice(
newpop, min(pop_size, len(newpop)),
replace=False, p=fitnesslist
).tolist()
else:
pop = newpop
for i in range(2):
yield output(
pop,
gen+i, pop_size, cell_size, mutation_rate, a_p, boa_split="aft"
)
print("{} generations are done. Cells are extinct.".format(gen))
if gen % 100 == 0 and pop:
yield output(
pop,
gen, pop_size, cell_size, mutation_rate, a_p, boa_split="aft"
)
if gen % 1000 == 0 and pop:
print("{} generations are done.".format(gen))
def writeoutfile(simulationresult, runnumber):
localtime = time.strftime(
"%m_%d_%H_%M_%S_%Y", time.localtime(time.time())
)
with open("output_data_" + localtime + ".csv", "w", newline="") as outfile:
outfile.write(
"gen"+";" +
"eSzamSum"+";"+"pSzamSum"+";" +
"alive"+";" +
"eSzamAtl"+";"+"eSzamVar"+";" +
"eSzamAKv"+";" +
"eSzamMed"+";" +
"eSzamFKv"+";" +
"pSzamAtl"+";" + "pSzamVar" + ";" +
"pSzamAKv"+";" +
"pSzamMed"+";" +
"pSzamFKv"+";" +
"fitAtl"+";"+"fitVar"+";" +
"fitAKv"+";" +
"fitMed"+";" +
"fitFKv"+";" +
"N"+";"+"rMax"+";"+"U"+";"+"aP"+";"+"boaSplit"+"\n"
)
outfile = csv.writer(outfile, delimiter=";")
counter = 0
print(counter, "/", runnumber)
for i in simulationresult:
outfile.writerows(i)
counter += 1
print(counter, "/", runnumber)
RESULT = [simulation(100, 20, 1, 0, 10000)]
RESULT.append(simulation(100, 20, 1, 1, 10000))
N_RUN = 2
writeoutfile(RESULT, N_RUN)
# Normally I call the functions from another script,
# these last 4 lines are meant to be an example.
```
