**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: `"enzyme"` (formerly `"e"`) for the enzymes and `"parasite"` (formerly `"p"`) for the parasites. The values of the keys are the numbers of the 2 variants.
Our parameters are:
`population_size` (formerly `pop_size`): the number of the cells
`cell_size`: the maximal number of molecules (enzymes+parasites) of cells at which cell division takes place
`replication_rate_p` (formerly `a_p`): fitness (replication rate) of the parasites relative to the fitness (replication rate) of the enzymes (for example if `replication_rate_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 `population_size` cells with `cell_size // 2` enzymes 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 `"enzyme" < 2` are discarded as dead. After that if the number of viable cells is bigger than `population_size`, we choose `population_size` of them according to cell fitness (`"enzyme"/("enzyme"+"parasite")`), and they move on to the next generation. On the other hand, if the number of viable cells is `population_size` or less, they all move on to the next generation.
**My request**
I've never studied programming in school. The original program is the result of heavy googling. I'm at a point where I need advice from experienced people. At certain parameter values the program gets quite slow.
So far combinations of these parameter values were examined:
`population_size`: 100; 200; 500; 1000
`cell_size`: 20; 50; 100; 200; 500; 1000
`replication_rate_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 `population_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 `replication_rate_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 `population_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.
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?
In which directions should I improve my programming skills in Python to efficiently implement these kind of models?
Or am I near the full potential of Python 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, `write_out_file()` (formerly `writeoutfile()`) writes the data to file.
Updated code with stylistic changes based on HoboProber's answer. The original code is below the updated one.
```python
# -*- coding: utf-8 -*-
"""
Simulation of an enzyme-parasite system,
explicitly assuming compartmentalization.
"""
from random import choices, random
import csv
import time
import numpy as np
def simulation(population_size, cell_size, replication_rate_p, mutation_rate, gen_max):
"""Function doing the simulation."""
def fitness(population):
"""Function for calculating fitnesses of cells."""
return [i["enzyme"] / (i["enzyme"] + i["parasite"]) for i in population]
def population_stats(population):
"""Function for calculatig statistics of the system."""
gyaklist_e = [i["enzyme"] for i in population]
gyaklist_p = [i["parasite"] for i in population]
fitness_list = fitness(population)
return (
sum(gyaklist_e), sum(gyaklist_p),
sum([1 for i in population if i["enzyme"] > 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(fitness_list), np.var(fitness_list),
np.percentile(fitness_list, 25),
np.percentile(fitness_list, 50),
np.percentile(fitness_list, 75)
)
# Creating the system with the starting state being
# half full cells containing only enzymes.
population = [{"enzyme": cell_size // 2,
"parasite": 0} for _ in range(population_size)]
gen = 0
yield (gen, *population_stats(population), population_size,
cell_size, mutation_rate, replication_rate_p, "aft")
print(f"N = {population_size}, rMax = {cell_size}, "
f"aP = {replication_rate_p}, U = {mutation_rate}")
while population and gen < gen_max:
gen += 1
# Replicator proliferation until cell_size in each cell.
for i in population:
while not i["enzyme"] + i["parasite"] == cell_size:
lucky_replicator = choices(
["enzyme", "parasite"],
[i["enzyme"], replication_rate_p*i["parasite"]]
)
if lucky_replicator[0] == "enzyme" and random() < mutation_rate:
lucky_replicator[0] = "parasite"
i[lucky_replicator[0]] += 1
if gen % 100 == 0:
yield (gen, *population_stats(population), population_size,
cell_size, mutation_rate, replication_rate_p, "bef")
# Each cell divides.
new_population = [{"enzyme": np.random.binomial(i["enzyme"], 0.5),
"parasite": np.random.binomial(i["parasite"], 0.5)}
for i in population]
for i in zip(population, new_population):
i[0]["enzyme"] -= i[1]["enzyme"]
i[0]["parasite"] -= i[1]["parasite"]
population += new_population
# Discarding dead cells.
new_population = [i for i in population if i["enzyme"] > 1]
# Choosing surviving cells according to their fitnesses
# if there are viable cells.
# If there are no viable cells, the simulation is over.
if new_population:
fitness_list = fitness(new_population)
fitness_sum = np.sum(fitness_list)
fitness_list = fitness_list / fitness_sum
population = np.random.choice(
new_population, min(population_size, len(new_population)),
replace=False, p=fitness_list
).tolist()
else:
population = new_population
for i in range(2):
yield (gen+i, *(0, 0)*9, population_size,
cell_size, mutation_rate, replication_rate_p, "aft")
print(f"{gen} generations are done. Cells are extinct.")
if gen % 100 == 0 and population:
yield (gen, *population_stats(population), population_size,
cell_size, mutation_rate, replication_rate_p, "aft")
if gen % 1000 == 0 and population:
print(f"{gen} generations are done.")
def write_out_file(result, n_run):
"""Function writing csv output file."""
local_time = time.strftime("%m_%d_%H_%M_%S_%Y", time.localtime(time.time()))
with open("output_data_" + local_time + ".csv", "w", newline="") as out_file:
out_file.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"
)
out_file = csv.writer(out_file, delimiter=";")
counter = 0
print(counter, "/", n_run)
for i in result:
out_file.writerows(i)
counter += 1
print(counter, "/", n_run)
RESULT = [simulation(100, 20, 1, 0, 10000)]
RESULT.append(simulation(100, 20, 1, 1, 10000))
N_RUN = 2
write_out_file(RESULT, N_RUN)
# Normally I call the functions from another script,
# these last 4 lines are meant to be just an example.
```
The original code:
```python
# -*- 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.
```