13
\$\begingroup\$

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 or making (explicit and/or implicit) array copies unnecessarily many times.
  • Is there a better way performance-wise to write data to file?

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 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 = population[0, :]/population.sum(axis=0)
        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=np.int32)
    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}",
          file=DEAD_OR_ALIVE)

    while (population.size > 0) & (gen < gen_max):
        gen += 1

        # Replicator proliferation until cell_size in each cell.
        mask = (population.sum(axis=0) < cell_size).nonzero()
        while mask[0].size > 0:
            # Calculating probabilites of choosing a parasite to replication.
            repl_probs_p = population[:, mask].copy()
            repl_probs_p.view(np.float32)[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=np.int32)
            lucky_replicators[
                (repl_choices < repl_probs_p) | (mut_choices < mutation_rate)
                ] = 1
            population[lucky_replicators, mask] += 1
            mask = (population.sum(axis=0) < cell_size).nonzero()

        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.random.binomial(population, 0.5)
        population -= new_population

        # Discarding dead cells.
        population = np.concatenate((population[:, (population[0, :] > 1).nonzero()[0]],
                                     new_population[:, (new_population[0, :] > 1).nonzero()[0]]),
                                    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 = population[0, :]/population.sum(axis=0)
            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.")
            print("Cells are extinct.", file=DEAD_OR_ALIVE)

        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.")

    print("Simulation ended successfully.\n", file=DEAD_OR_ALIVE)


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

    Parameters
    ----------
    result : list of generator object(s)
        Contains data on the simulated system.

    n_run : int
        The number of consecutive runs.
    """
    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
        for i in result:
            out_file.writerows(i)
            counter += 1
            print(counter, "/", n_run, "\n")


LOCAL_TIME = time.strftime("%m_%d_%H_%M_%S_%Y", time.localtime(time.time()))
DEAD_OR_ALIVE = open("output_data_" + LOCAL_TIME + ".txt", "w")
RESULT = [simulation(1000, 200, 1.5, 0.0, 10000)]
#RESULT.append(simulation(1000, 200, 1.5, 1.0, 10000))
N_RUN = 1
write_out_file(RESULT, LOCAL_TIME, N_RUN)
DEAD_OR_ALIVE.close()
# Normally I call the functions from another script,
# these last 4 lines are meant to be just an example.

line_profiling

Timer unit: 1e-07 s

Total time: 161.05 s
File: simulation.py
Function: simulation at line 16

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    16
    17                                           def simulation(population_size, cell_size, replication_rate_p, mutation_rate, gen_max):
    18                                               """
    19                                               Simulate an enzyme-parasite system explicitly assuming compartmentalization.
    20
    21                                               Parameters
    22                                               ----------
    23                                               population_size : int
    24                                                   The number of cells.
    25
    26                                               cell_size : int
    27                                                   The maximal number of replicators of cells at which cell division takes place.
    28
    29                                               replication_rate_p : float
    30                                                   The fitness (replication rate) of the parasites
    31                                                   relative to the fitness (replication rate) of the enzymes.
    32                                                   Example
    33                                                   -------
    34                                                       $ replication_rate_p = 2
    35                                                   This means that the parasites' fitness is twice as that of the enzymes.
    36
    37                                               mutation_rate : float
    38                                                   The probability of mutation during a replication event.
    39
    40                                               gen_max : int
    41                                                   The maximal number of generations.
    42                                                   A generation corresponds to one outer while cycle.
    43                                                   If the system extincts, the number of generations doesn't reach gen_max.
    44
    45                                               Yield
    46                                               -------
    47                                               generator object
    48                                                   Contains data on the simulated system.
    49                                               """
    50
    51         1         56.0     56.0      0.0      def population_stats(population):
    52                                                   """
    53                                                   Calculate statistics of the system.
    54
    55                                                   Parameter
    56                                                   ---------
    57                                                   population : ndarray
    58                                                       The system itself.
    59
    60                                                   Return
    61                                                   -------
    62                                                   tuple
    63                                                       Contains statistics of the simulated system.
    64                                                   """
    65                                                   gyak_sums = population.sum(axis=1)
    66                                                   gyak_means = population.mean(axis=1)
    67                                                   gyak_variances = population.var(axis=1)
    68                                                   gyak_percentiles_25 = np.percentile(population, 25, axis=1)
    69                                                   gyak_medians = np.median(population, axis=1)
    70                                                   gyak_percentiles_75 = np.percentile(population, 75, axis=1)
    71                                                   fitness_list = population[0, :]/population.sum(axis=0)
    72                                                   return (
    73                                                       gyak_sums[0], gyak_sums[1], (population[0, :] > 1).sum(),
    74                                                       gyak_means[0], gyak_variances[0],
    75                                                       gyak_percentiles_25[0], gyak_medians[0], gyak_percentiles_75[0],
    76                                                       gyak_means[1], gyak_variances[1],
    77                                                       gyak_percentiles_25[1], gyak_medians[1], gyak_percentiles_75[1],
    78                                                       fitness_list.mean(), fitness_list.var(),
    79                                                       np.percentile(fitness_list, 25),
    80                                                       np.median(fitness_list),
    81                                                       np.percentile(fitness_list, 75)
    82                                                       )
    83
    84                                               # Creating the system with the starting state being
    85                                               # half full cells containing only enzymes.
    86         1         68.0     68.0      0.0      population = np.zeros((2, population_size), dtype=np.int32)
    87         1         53.0     53.0      0.0      population[0, :] = cell_size//2
    88         1          9.0      9.0      0.0      gen = 0
    89         1      14828.0  14828.0      0.0      yield (gen, *population_stats(population), population_size,
    90         1         24.0     24.0      0.0             cell_size, mutation_rate, replication_rate_p, "aft")
    91         1         49.0     49.0      0.0      print(f"N = {population_size}, rMax = {cell_size}, "
    92                                                     f"aP = {replication_rate_p}, U = {mutation_rate}",
    93         1        113.0    113.0      0.0            file=DEAD_OR_ALIVE)
    94
    95     10001     140323.0     14.0      0.0      while (population.size > 0) & (gen < gen_max):
    96     10000     123102.0     12.3      0.0          gen += 1
    97
    98                                                   # Replicator proliferation until cell_size in each cell.
    99     10000    3333616.0    333.4      0.2          mask = (population.sum(axis=0) < cell_size).nonzero()
   100   1238245   20308315.0     16.4      1.3          while mask[0].size > 0:
   101                                                       # Calculating probabilites of choosing a parasite to replication.
   102   1228245  239761224.0    195.2     14.9              repl_probs_p = population[:, mask].copy()
   103   1228245   83589799.0     68.1      5.2              repl_probs_p.view(np.float32)[1, :] *= replication_rate_p
   104   1228245  158300271.0    128.9      9.8              repl_probs_p = repl_probs_p[1, :]/repl_probs_p.sum(axis=0)
   105                                                       # Determining if an enzyme or a parasite replicates,
   106                                                       # and if an enzyme replicates, will it mutate to a parasite.
   107                                                       # (Outcome can differ among cells. Parasites don't mutate.)
   108   1228245  132808465.0    108.1      8.2              repl_choices = np.random.random_sample(repl_probs_p.shape)
   109   1228245  117430558.0     95.6      7.3              mut_choices = np.random.random_sample(repl_probs_p.shape)
   110   1228245   35120008.0     28.6      2.2              lucky_replicators = np.zeros(repl_probs_p.shape, dtype=np.int32)
   111                                                       lucky_replicators[
   112                                                           (repl_choices < repl_probs_p) | (mut_choices < mutation_rate)
   113   1228245   76236137.0     62.1      4.7                  ] = 1
   114   1228245  301823109.0    245.7     18.7              population[lucky_replicators, mask] += 1
   115   1228245  357660422.0    291.2     22.2              mask = (population.sum(axis=0) < cell_size).nonzero()
   116
   117     10000     143547.0     14.4      0.0          if gen % 100 == 0:
   118       100    1350075.0  13500.8      0.1              yield (gen, *population_stats(population), population_size,
   119       100       2544.0     25.4      0.0                     cell_size, mutation_rate, replication_rate_p, "bef")
   120
   121                                                   # Each cell divides.
   122     10000   17525435.0   1752.5      1.1          new_population = np.random.binomial(population, 0.5)
   123     10000    1087713.0    108.8      0.1          population -= new_population
   124
   125                                                   # Discarding dead cells.
   126     10000    2526633.0    252.7      0.2          population = np.concatenate((population[:, (population[0, :] > 1).nonzero()[0]],
   127     10000    1979199.0    197.9      0.1                                       new_population[:, (new_population[0, :] > 1).nonzero()[0]]),
   128     10000    1003433.0    100.3      0.1                                      axis=1)
   129
   130                                                   # Choosing survivor cells according to their fitnesses
   131                                                   # if there are more viable cells than population_size.
   132                                                   # Hence population_size or less cells move on to the next generation.
   133     10000     184360.0     18.4      0.0          if population.shape[1] > population_size:
   134     10000    5107803.0    510.8      0.3              fitness_list = population[0, :]/population.sum(axis=0)
   135     10000    1244299.0    124.4      0.1              fitness_list = fitness_list/fitness_list.sum()
   136     10000     213078.0     21.3      0.0              population = population[:, np.random.choice(population.shape[1],
   137     10000     110896.0     11.1      0.0                                                          population_size,
   138     10000     111486.0     11.1      0.0                                                          replace=False,
   139     10000   49497963.0   4949.8      3.1                                                          p=fitness_list)]
   140                                                   elif population.size == 0:
   141                                                       for i in range(2):
   142                                                           yield (gen+i, *(0, 0)*9, population_size,
   143                                                                  cell_size, mutation_rate, replication_rate_p, "aft")
   144                                                       print(f"{gen} generations are done.")
   145                                                       print("Cells are extinct.", file=DEAD_OR_ALIVE)
   146
   147     10000     260742.0     26.1      0.0          if (gen % 100 == 0) & (population.size > 0):
   148       100    1332898.0  13329.0      0.1              yield (gen, *population_stats(population), population_size,
   149       100       2553.0     25.5      0.0                     cell_size, mutation_rate, replication_rate_p, "aft")
   150
   151     10000     147525.0     14.8      0.0          if (gen % 1000 == 0) & (population.size > 0):
   152        10      21265.0   2126.5      0.0              print(f"{gen} generations are done.")
   153
   154         1        226.0    226.0      0.0      print("Simulation ended successfully.\n", file=DEAD_OR_ALIVE)

cProfiling sample

Fri Nov 29 04:53:01 2019    cprofiling

         16375164 function calls (16361694 primitive calls) in 135.937 seconds

   Ordered by: internal time, cumulative time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
      202   72.331    0.358  135.766    0.672 simulation.py:17(simulation)
  2529183   27.246    0.000   27.246    0.000 {method 'reduce' of 'numpy.ufunc' objects}
  2456168   20.346    0.000   20.346    0.000 {method 'random_sample' of 'numpy.random.mtrand.RandomState' objects}
    10000    2.575    0.000    4.456    0.000 {method 'choice' of 'numpy.random.mtrand.RandomState' objects}
  1258084    2.326    0.000    2.326    0.000 {method 'nonzero' of 'numpy.ndarray' objects}
  1228747    2.139    0.000    2.139    0.000 {method 'copy' of 'numpy.ndarray' objects}
  2486771    2.043    0.000   29.905    0.000 {method 'sum' of 'numpy.ndarray' objects}
  1228085    1.420    0.000    1.420    0.000 {built-in method numpy.zeros}
    10000    1.354    0.000    1.683    0.000 {method 'binomial' of 'numpy.random.mtrand.RandomState' objects}
1228088/1228087    0.899    0.000    0.899    0.000 {method 'view' of 'numpy.ndarray' objects}
  2486771    0.783    0.000   27.862    0.000 _methods.py:36(_sum)
    31404    0.585    0.000    0.585    0.000 {method 'argsort' of 'numpy.ndarray' objects}
    31404    0.413    0.000    1.081    0.000 arraysetops.py:297(_unique1d)
    31404    0.262    0.000    0.262    0.000 {method 'cumsum' of 'numpy.ndarray' objects}
134267/124016    0.162    0.000    2.224    0.000 {built-in method numpy.core._multiarray_umath.implement_array_function}
    40804    0.103    0.000    0.334    0.000 fromnumeric.py:73(_wrapreduction)
    31404    0.064    0.000    1.193    0.000 arraysetops.py:151(unique)
    32007    0.039    0.000    0.039    0.000 {method 'flatten' of 'numpy.ndarray' objects}
    31404    0.034    0.000    0.329    0.000 fromnumeric.py:2358(cumsum)
    20000    0.032    0.000    0.092    0.000 {method 'all' of 'numpy.generic' objects}
    31405    0.031    0.000    0.031    0.000 {built-in method numpy.empty}
      804    0.027    0.000    0.111    0.000 function_base.py:3853(_quantile_ureduce_func)
    31404    0.027    0.000    0.382    0.000 <__array_function__ internals>:2(cumsum)
    31404    0.027    0.000    1.256    0.000 <__array_function__ internals>:2(unique)
    68944    0.027    0.000    0.027    0.000 {built-in method numpy.array}
      667    0.025    0.000    0.025    0.000 {built-in method nt.stat}
    33012    0.025    0.000    0.303    0.000 fromnumeric.py:55(_wrapfunc)
    67140    0.025    0.000    0.025    0.000 {built-in method builtins.getattr}
    20000    0.024    0.000    0.029    0.000 getlimits.py:365(__new__)
    40804    0.021    0.000    0.021    0.000 fromnumeric.py:74(<dictcomp>)
    20000    0.021    0.000    0.189    0.000 fromnumeric.py:2277(all)
    24824    0.020    0.000    0.030    0.000 numerictypes.py:293(issubclass_)
    67230    0.020    0.000    0.045    0.000 _asarray.py:88(asanyarray)
    20000    0.019    0.000    0.243    0.000 <__array_function__ internals>:2(all)
    12412    0.019    0.000    0.050    0.000 numerictypes.py:365(issubdtype)
     9045    0.017    0.000    0.025    0.000 numeric.py:1273(normalize_axis_tuple)
      139    0.016    0.000    0.021    0.000 <frozen importlib._bootstrap_external>:914(get_data)
    31404    0.016    0.000    0.021    0.000 arraysetops.py:138(_unpack_tuple)
    10000    0.015    0.000    0.116    0.000 fromnumeric.py:2792(prod)
       19    0.015    0.001    0.017    0.001 {built-in method _imp.create_dynamic}
      317    0.014    0.000    0.014    0.000 {built-in method builtins.compile}
     4221    0.014    0.000    0.043    0.000 numeric.py:1336(moveaxis)
      139    0.014    0.000    0.014    0.000 {built-in method marshal.loads}
    11207    0.012    0.000    0.064    0.000 <__array_function__ internals>:2(concatenate)
    39330    0.011    0.000    0.011    0.000 {built-in method builtins.issubclass}
    10000    0.011    0.000    0.139    0.000 <__array_function__ internals>:2(prod)
    11608    0.011    0.000    0.011    0.000 {built-in method numpy.core._multiarray_umath.count_nonzero}
    11608    0.010    0.000    0.037    0.000 <__array_function__ internals>:2(count_nonzero)
      402    0.010    0.000    0.023    0.000 _methods.py:167(_var)
    10804    0.010    0.000    0.093    0.000 <__array_function__ internals>:2(any)
     1206    0.010    0.000    0.010    0.000 {method 'partition' of 'numpy.ndarray' objects}
    10804    0.009    0.000    0.074    0.000 fromnumeric.py:2189(any)
62590/62386    0.008    0.000    0.008    0.000 {built-in method builtins.len}
    40846    0.007    0.000    0.007    0.000 {method 'items' of 'dict' objects}
    20000    0.007    0.000    0.059    0.000 _methods.py:47(_all)
      804    0.006    0.000    0.017    0.000 _methods.py:134(_mean)
     1608    0.006    0.000    0.006    0.000 {method 'take' of 'numpy.ndarray' objects}
    11608    0.006    0.000    0.017    0.000 numeric.py:409(count_nonzero)
    31404    0.006    0.000    0.006    0.000 fromnumeric.py:2354(_cumsum_dispatcher)
     1206    0.006    0.000    0.145    0.000 function_base.py:3359(_ureduce)
    21762    0.005    0.000    0.005    0.000 {method 'get' of 'dict' objects}
    31404    0.005    0.000    0.005    0.000 arraysetops.py:146(_unique_dispatcher)
      139    0.005    0.000    0.005    0.000 {method 'read' of '_io.FileIO' objects}
  342/339    0.004    0.000    0.006    0.000 {built-in method builtins.__build_class__}
      201    0.004    0.000    0.211    0.001 simulation.py:51(population_stats)
      804    0.004    0.000    0.133    0.000 function_base.py:3569(percentile)
        1    0.004    0.004  135.770  135.770 {method 'writerows' of '_csv.writer' objects}
    20000    0.004    0.000    0.004    0.000 fromnumeric.py:2273(_all_dispatcher)
      804    0.004    0.000    0.009    0.000 function_base.py:3840(_quantile_is_valid)
      402    0.004    0.000    0.025    0.000 function_base.py:3508(_median)
       13    0.003    0.000    0.003    0.000 {built-in method builtins.print}
      642    0.003    0.000    0.003    0.000 {method 'sub' of 're.Pattern' objects}
     9045    0.003    0.000    0.005    0.000 numeric.py:1323(<listcomp>)
     4221    0.003    0.000    0.049    0.000 <__array_function__ internals>:2(moveaxis)
       16    0.003    0.000    0.003    0.000 {built-in method nt.listdir}
      322    0.002    0.000    0.029    0.000 <frozen importlib._bootstrap_external>:1356(find_spec)
    11207    0.002    0.000    0.002    0.000 multiarray.py:145(concatenate)
    10000    0.002    0.000    0.002    0.000 fromnumeric.py:2787(_prod_dispatcher)
     4221    0.002    0.000    0.002    0.000 {method 'transpose' of 'numpy.ndarray' objects}
     4222    0.002    0.000    0.002    0.000 {built-in method builtins.sorted}
     9045    0.002    0.000    0.002    0.000 {built-in method numpy.core._multiarray_umath.normalize_axis_index}
    11608    0.002    0.000    0.002    0.000 numeric.py:405(_count_nonzero_dispatcher)
     1206    0.002    0.000    0.002    0.000 _methods.py:50(_count_reduce_items)
    10804    0.002    0.000    0.002    0.000 fromnumeric.py:2185(_any_dispatcher)
   101/33    0.002    0.000    0.004    0.000 sre_parse.py:469(_parse)
      201    0.002    0.000    0.005    0.000 utils.py:1142(_median_nancheck)
      321    0.002    0.000    0.002    0.000 {method 'findall' of 're.Pattern' objects}
     9499    0.001    0.000    0.001    0.000 {built-in method builtins.isinstance}
    19/14    0.001    0.000    0.011    0.001 {built-in method _imp.exec_dynamic}
    469/1    0.001    0.000  135.938  135.938 {built-in method builtins.exec}
     1608    0.001    0.000    0.009    0.000 fromnumeric.py:97(take)
      614    0.001    0.000    0.002    0.000 _inspect.py:67(getargs)
     1608    0.001    0.000    0.012    0.000 <__array_function__ internals>:2(take)
     3189    0.001    0.000    0.001    0.000 {built-in method builtins.hasattr}
      139    0.001    0.000    0.043    0.000 <frozen importlib._bootstrap_external>:793(get_code)
      804    0.001    0.000    0.119    0.000 function_base.py:3828(_quantile_unchecked)
    182/2    0.001    0.000    0.165    0.083 <frozen importlib._bootstrap>:978(_find_and_load)
     4221    0.001    0.000    0.001    0.000 numeric.py:1399(<listcomp>)
     4226    0.001    0.000    0.001    0.000 {method 'insert' of 'list' objects}
      287    0.001    0.000    0.004    0.000 overrides.py:72(verify_matching_signatures)
      317    0.001    0.000    0.029    0.000 overrides.py:154(decorator)
     1555    0.001    0.000    0.003    0.000 <frozen importlib._bootstrap_external>:56(_path_join)
      179    0.001    0.000    0.034    0.000 <frozen importlib._bootstrap>:882(_find_spec)
      339    0.001    0.000    0.002    0.000 functools.py:37(update_wrapper)
   190/31    0.001    0.000    0.003    0.000 sre_compile.py:71(_compile)
     9045    0.001    0.000    0.001    0.000 {built-in method _operator.index}
       77    0.001    0.000    0.001    0.000 sre_compile.py:276(_optimize_charset)
     1555    0.001    0.000    0.001    0.000 <frozen importlib._bootstrap_external>:58(<listcomp>)
      402    0.001    0.000    0.007    0.000 fromnumeric.py:3153(mean)
      804    0.001    0.000    0.001    0.000 {method 'astype' of 'numpy.ndarray' objects}
      278    0.001    0.000    0.002    0.000 <frozen importlib._bootstrap_external>:271(cache_from_source)
      481    0.001    0.000    0.002    0.000 <frozen importlib._bootstrap>:157(_get_module_lock)
       16    0.001    0.000    0.002    0.000 <frozen importlib._bootstrap_external>:1190(_path_hooks)
      321    0.001    0.000    0.007    0.000 textwrap.py:414(dedent)
        2    0.001    0.000    0.001    0.000 {built-in method _ctypes.LoadLibrary}
      756    0.001    0.000    0.001    0.000 {method 'format' of 'str' objects}
      481    0.001    0.000    0.001    0.000 <frozen importlib._bootstrap>:78(acquire)
      804    0.001    0.000    0.135    0.000 <__array_function__ internals>:2(percentile)
      366    0.001    0.000    0.001    0.000 {built-in method _thread.allocate_lock}
     1608    0.001    0.000    0.001    0.000 {method 'squeeze' of 'numpy.ndarray' objects}
      162    0.001    0.000    0.032    0.000 <frozen importlib._bootstrap_external>:1240(_get_spec)
      175    0.001    0.000    0.003    0.000 <frozen importlib._bootstrap>:504(_init_module_attrs)
    175/2    0.001    0.000    0.164    0.082 <frozen importlib._bootstrap>:663(_load_unlocked)
   882/71    0.001    0.000    0.146    0.002 <frozen importlib._bootstrap>:1009(_handle_fromlist)
      618    0.001    0.000    0.003    0.000 _inspect.py:98(getargspec)
      481    0.001    0.000    0.001    0.000 <frozen importlib._bootstrap>:103(release)
       17    0.001    0.000    0.001    0.000 {built-in method _imp.create_builtin}
      634    0.001    0.000    0.001    0.000 {built-in method __new__ of type object at 0x00007FFFE42159A0}
      455    0.001    0.000    0.010    0.000 re.py:271(_compile)
      278    0.001    0.000    0.001    0.000 <frozen importlib._bootstrap_external>:62(_path_split)
      402    0.001    0.000    0.006    0.000 fromnumeric.py:657(partition)
     4221    0.001    0.000    0.001    0.000 numeric.py:1332(_moveaxis_dispatcher)
    182/2    0.001    0.000    0.165    0.083 <frozen importlib._bootstrap>:948(_find_and_load_unlocked)
       12    0.001    0.000    0.001    0.000 __init__.py:316(namedtuple)
     2064    0.001    0.000    0.001    0.000 {method 'join' of 'str' objects}

Of course any advice is highly appreciated!=)

\$\endgroup\$
1
  • 1
    \$\begingroup\$ 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. \$\endgroup\$
    – IEatBagels
    Aug 3, 2019 at 22:59

1 Answer 1

3
\$\begingroup\$

Tuple returns

    """
    Return
    -------
    tuple
        Contains statistics of the simulated system.
    """
    ...
    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)
        )

First of all - if you're going to bother documenting the function, it would be important to describe every one of these values. However, the easier and significantly more maintainable thing to do is return an object of some kind; pick your flavour - a plain-old class, a data class, a named tuple, what-have-you. These would all allow for you to return one thing whose members are self-documenting, instead of requiring magical knowledge of position to access them.

Logical, not bit-wise, operators

while (population.size > 0) & (gen < gen_max):

The only time I've seen syntax like this in Python is for SQLAlchemy, which does some dirty tricks to produce SQL from vaguely boolean-smelling expressions. However, it's much more likely that you actually mean:

while population.size > 0 and gen < gen_max:

since and is logical and & is bit-wise. It's also worth noting that you should Loop Like a Native, and instead of incrementing gen manually, do

for gen in range(gen_max):
    if population_size <= 0:
        break

Type hints

This is somewhat of an educated guess, but

def write_out_file(result, local_time, n_run):

can be

def write_out_file(result: List[Iterable[int]], local_time: datetime, n_run: int):

It looks (though it is missing from the documentation) that local_time is actually passed in as a string, but it shouldn't be. Stringification should in this case be left to the function itself.

Global code

This stuff:

LOCAL_TIME = time.strftime("%m_%d_%H_%M_%S_%Y", time.localtime(time.time()))
DEAD_OR_ALIVE = open("output_data_" + LOCAL_TIME + ".txt", "w")
RESULT = [simulation(1000, 200, 1.5, 0.0, 10000)]
#RESULT.append(simulation(1000, 200, 1.5, 1.0, 10000))
N_RUN = 1
write_out_file(RESULT, LOCAL_TIME, N_RUN)
DEAD_OR_ALIVE.close()

has a few problems:

  • That code blob should be in a main function
  • Once that happens, you can de-capitalize those variable names.
  • DEAD_OR_ALIVE should be put into a with block

Use enumerate

This:

    counter = 0
    for i in result:
        out_file.writerows(i)
        counter += 1
        print(counter, "/", n_run, "\n")

should be

for counter, i in enumerate(result):
   out_file.writerows(i)
   print(f'{counter}/{n_run}')
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.