# Population dynamic simulation on biological information maintenance 2

This question is the follow-up to this previous question.

## 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.

I have made the changes advised by HoboProber. Namely correction of style and implementing the model relying on Numpy. So now the system is a 2-dimensional array. Cells are the columns of the array. The values of the first row are the numbers of enzymes and the values of the second row are the numbers of parasites.

## My request

The speed of this newer implementation is much better than that of the previous one. But as I would like to increase population_size and gen_max every bit of performance improvement counts.

So far I examined the system in more detail with population sizes ranging from 100 to 1000 cells and with the maximal number of generations being 10000. The amount of increase in population size depends on performance, a million cells would be a perfectly reasonable assumption concerning the modelled system. The maximal number of generations should be 20-30000.

• Primarily, does the code make use of vectorization and Numpy as effectively as it can?
• Which potential efficiency improvements I missed? For example calculating something multiple times instead of assigning it to a variable.
• Is there a better way performance-wise to write data to file?

Could performance significantly profit from, so would it be worth to learn...

• ...multiprocessing or multithreading? For example when the proliferation of molecules happens, the replication events in different cells are independent of each other. They could happen in parallel maybe.
• ...static typing and compiling? For example with Cython or Numba.

I anticipate this to be borderline opinion based. If it's problematic, I'll remove this section.

## The code

"""
Collect data on an enzyme-parasite system explicitly assuming compartmentalization.

Functions
---------
simulation()
Simulate mentioned system.

write_out_file()
Write data to csv output file.
"""
import csv
import time
import numpy as np

def simulation(population_size, cell_size, replication_rate_p, mutation_rate, gen_max):
"""
Simulate an enzyme-parasite system explicitly assuming compartmentalization.

Parameters
----------
population_size : int
The number of cells.

cell_size : int
The maximal number of replicators of cells at which cell division takes place.

replication_rate_p : float
The fitness (replication rate) of the parasites
relative to the fitness (replication rate) of the enzymes.
Example
-------
\$ replication_rate_p = 2
This means that the parasites' fitness is twice as that of the enzymes.

mutation_rate : float
The probability of mutation during a replication event.

gen_max : int
The maximal number of generations.
A generation corresponds to one outer while cycle.
If the system extincts, the number of generations doesn't reach gen_max.

Yield
-------
generator object
Contains data on the simulated system.
"""

def fitness(population):
"""
Calculate fitnesses of cells.
fitness of a cell = number of enzymes/(number of enzymes + number of parasites)

Parameter
---------
population : ndarray
The system itself.

Return
------
ndarray
The fitness of each cell of the system.
"""
return population[0, :]/population.sum(axis=0)

def population_stats(population):
"""
Calculate statistics of the system.

Parameter
---------
population : ndarray
The system itself.

Return
-------
tuple
Contains statistics of the simulated system.
"""
gyak_sums = population.sum(axis=1)
gyak_means = population.mean(axis=1)
gyak_variances = population.var(axis=1)
gyak_percentiles_25 = np.percentile(population, 25, axis=1)
gyak_medians = np.median(population, axis=1)
gyak_percentiles_75 = np.percentile(population, 75, axis=1)
fitness_list = fitness(population)
return (
gyak_sums[0], gyak_sums[1], (population[0, :] > 1).sum(),
gyak_means[0], gyak_variances[0],
gyak_percentiles_25[0], gyak_medians[0], gyak_percentiles_75[0],
gyak_means[1], gyak_variances[1],
gyak_percentiles_25[1], gyak_medians[1], gyak_percentiles_75[1],
fitness_list.mean(), fitness_list.var(),
np.percentile(fitness_list, 25),
np.median(fitness_list),
np.percentile(fitness_list, 75)
)

# Creating the system with the starting state being
# half full cells containing only enzymes.
population = np.zeros((2, population_size), dtype=int)
population[0, :] = cell_size//2
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.size > 0) & (gen < gen_max):
gen += 1

# Replicator proliferation until cell_size in each cell.
# Calculating probabilites of choosing a parasite to replication.
repl_probs_p.astype(float)[1, :] *= replication_rate_p
repl_probs_p = repl_probs_p[1, :]/repl_probs_p.sum(axis=0)
# Determining if an enzyme or a parasite replicates,
# and if an enzyme replicates, will it mutate to a parasite.
# (Outcome can differ among cells. Parasites don't mutate.)
repl_choices = np.random.random_sample(repl_probs_p.shape)
mut_choices = np.random.random_sample(repl_probs_p.shape)
lucky_replicators = np.zeros(repl_probs_p.shape, dtype=int)
lucky_replicators[
(repl_choices < repl_probs_p) | (mut_choices < mutation_rate)
] = 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 = np.empty_like(population)
new_population[0, :] = np.random.binomial(population[0, :], 0.5)
new_population[1, :] = np.random.binomial(population[1, :], 0.5)
population -= new_population

population = np.concatenate([population[:, population[0, :] > 1],
new_population[:, new_population[0, :] > 1]],
axis=1)

# Choosing survivor cells according to their fitnesses
# if there are more viable cells than population_size.
# Hence population_size or less cells move on to the next generation.
if population.shape[1] > population_size:
fitness_list = fitness(population)
fitness_list = fitness_list/fitness_list.sum()
population = population[:, np.random.choice(population.shape[1],
population_size,
replace=False,
p=fitness_list)]
elif population.size == 0:
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) & (population.size > 0):
yield (gen, *population_stats(population), population_size,
cell_size, mutation_rate, replication_rate_p, "aft")

if (gen % 1000 == 0) & (population.size > 0):
print(f"{gen} generations are done.")

def write_out_file(result, n_run):
"""
Write data to csv output file.

Parameters
----------
result : generator object or list of generator objects
Contains data on the simulated system.

n_run : int
The number of consecutive runs.
"""
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(1000, 100, 1.5, 0.0, 10000)]
RESULT.append(simulation(1000, 100, 1.5, 1.0, 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.


### Profiling sample

Sun Jun  2 14:53:41 2019    profiling

11325844 function calls (11323591 primitive calls) in 116.690 seconds

Ordered by: cumulative time
List reduced from 1032 to 32 due to restriction <32>

ncalls  tottime  percall  cumtime  percall filename:lineno(function)
135/1    0.001    0.000  116.690  116.690 {built-in method builtins.exec}
1    0.000    0.000  116.690  116.690 simulation.py:11(<module>)
1    0.000    0.000  116.551  116.551 simulation.py:172(write_out_file)
2    0.005    0.003  116.549   58.275 {method 'writerows' of '_csv.writer' objects}
206   59.710    0.290  116.543    0.566 simulation.py:17(simulation)
2098070   23.118    0.000   23.118    0.000 {method 'reduce' of 'numpy.ufunc' objects}
1353888    1.480    0.000   22.700    0.000 {method 'sum' of 'numpy.ndarray' objects}
1353888    0.679    0.000   21.219    0.000 _methods.py:34(_sum)
1323254   16.553    0.000   16.553    0.000 {method 'random_sample' of 'mtrand.RandomState' objects}
10005    3.244    0.000    5.589    0.001 {method 'choice' of 'mtrand.RandomState' objects}
732561    0.887    0.000    5.214    0.000 fromnumeric.py:1933(any)
742566    1.602    0.000    4.461    0.000 fromnumeric.py:64(_wrapreduction)
20036    1.879    0.000    2.374    0.000 {method 'binomial' of 'mtrand.RandomState' objects}
662435    1.720    0.000    1.720    0.000 {method 'astype' of 'numpy.ndarray' objects}
662268    1.673    0.000    1.673    0.000 {method 'copy' of 'numpy.ndarray' objects}
31365    0.081    0.000    1.621    0.000 arraysetops.py:121(unique)
31365    0.551    0.000    1.469    0.000 arraysetops.py:268(_unique1d)
661629    1.070    0.000    1.070    0.000 {built-in method numpy.core.multiarray.zeros}
31365    0.811    0.000    0.811    0.000 {method 'argsort' of 'numpy.ndarray' objects}
10207    0.106    0.000    0.531    0.000 simulation.py:51(fitness)
31365    0.046    0.000    0.440    0.000 fromnumeric.py:2092(cumsum)
32981    0.032    0.000    0.406    0.000 fromnumeric.py:49(_wrapfunc)
31365    0.352    0.000    0.352    0.000 {method 'cumsum' of 'numpy.ndarray' objects}
202    0.004    0.000    0.290    0.001 simulation.py:68(population_stats)
1212    0.007    0.000    0.202    0.000 function_base.py:3199(_ureduce)
40072    0.077    0.000    0.198    0.000 {method 'any' of 'numpy.generic' objects}
742592    0.193    0.000    0.193    0.000 {method 'items' of 'dict' objects}
6    0.000    0.000    0.189    0.031 __init__.py:1(<module>)
808    0.004    0.000    0.175    0.000 function_base.py:3398(percentile)
808    0.002    0.000    0.166    0.000 function_base.py:3647(_quantile_unchecked)
808    0.035    0.000    0.155    0.000 function_base.py:3672(_quantile_ureduce_func)
10005    0.018    0.000    0.152    0.000 fromnumeric.py:2478(prod)


Of course any advice is highly appreciated!=)

• If performance is a really big concern you might want to use Pytorch or Tensorflow, not for the machine learning part, but for the usage of CUDA accelerated tensors. – IEatBagels Aug 3 at 22:59