**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.
1. 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?
2. 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?
3. 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.)
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.
On parameter values
So far combinations of these values were examined:
* `pop_size`: 100; 200; 500; 1000
* `cell_size`: 20; 50; 100; 200; 500; 1000
* `a_p`: 0.75; 1; 1.25; 1.5; 1.75; 2; 3
* `mutation_rate`: 0-1
* `gen_max`: 10000
Primarily I would like to increase `pop_size` and above 1000 cells the program is slower than I would prefer. Of course that's somewhat subjective, but for example a million cells would be a perfectly reasonable assumption and at that order of magnitude I think it's objectively impossibly slow.
The program also gets slower with the increase in `cell_size` and slightly slower with `a_p`, but for the time being I'm happy with the values of the former and the effect of the latter is tolerable.
The effect of the mutation rate on speed is also tolerable.
In addition to `pop_size`, `gen_max` should be increased and has significant effect on run time. I know I don't catch every extinction events with 10000 generations. 20000 would be better, 50000 would be quite enough and 100000 would be like cracking a nut with a sledgehammer.