生物信息维护的种群动态模拟

Question

这个问题是前一个问题的后续问题.

背景

通过这个模拟，我研究了一个酶在细胞中增殖的系统。在酶的复制过程中，寄生虫可能是由突变引起的。它们会使系统灭绝。我感兴趣的是参数空间共存在哪里是可能的。

我做了HoboProber建议的更改。即纠正风格，依托Numpy实施模式。所以现在这个系统是一个二维阵列。单元格是数组的列。第一行的数值是酶的数目，第二行的值是寄生虫的数目。

我的请求

这个较新的实现的速度比前一个要好得多。但是，正如我希望增加population_size和gen_max一样，每一点性能改进都很重要。

到目前为止，我更详细地研究了这个系统，种群规模从100个到1000个细胞，最大的世代数为10000。人口数量的增加取决于性能，对于模拟系统来说，一百万个细胞将是一个完全合理的假设。最大世代数应为20~30000代。

• 首先，代码是否尽可能有效地使用向量化和Numpy？
• 我错过了哪些潜在的效率改进？例如，多次计算某物，而不是将其分配给一个变量，或者使(显式和/或隐式)数组不必要地复制很多次。
• 是否有更好的方法--从性能上写数据到文件？

代码

"""
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样本

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()
    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 :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()
     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 :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 :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 :978(_find_and_load)
     4221    0.001    0.000    0.001    0.000 numeric.py:1399()
     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 :56(_path_join)
      179    0.001    0.000    0.034    0.000 :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 :58()
      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 :271(cache_from_source)
      481    0.001    0.000    0.002    0.000 :157(_get_module_lock)
       16    0.001    0.000    0.002    0.000 :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 :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 :1240(_get_spec)
      175    0.001    0.000    0.003    0.000 :504(_init_module_attrs)
    175/2    0.001    0.000    0.164    0.082 :663(_load_unlocked)
   882/71    0.001    0.000    0.146    0.002 :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 :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 :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 :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}

当然，任何建议都非常感谢！

Answer 1

元组返回

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

首先，如果您要麻烦地记录这个函数，那么描述这些值中的每一个都是很重要的。然而，更容易和更易于维护的事情是返回某种类型的对象；选择您的口味--一个简单的老类、一个数据类、一个命名的元组、什么的-您。这些都将允许你返回一个东西，其成员是自我记录，而不是需要魔法知识的位置访问他们。

逻辑，而不是按位排列，运算符

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

在Python中，我唯一见过这样的语法是针对SQLAlchemy的，它使用一些肮脏的技巧来从模糊的布尔式表达式中生成SQL。然而，更有可能的是，您实际上指的是：

while population.size > 0 and gen < gen_max:

因为and是逻辑的，而&是位的。还值得注意的是，您应该像本机一样循环，而不是手动递增gen，请执行

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

类型提示

这是个有理有据的猜测，但是

def write_out_file(result, local_time, n_run):

可以是

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

它看起来(虽然文档中缺少它) local_time实际上是作为一个字符串传入的，但是它不应该是字符串。在这种情况下，字符串应该留给函数本身。

全局代码

这些东西：

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()

有几个问题：

• 这个代码块应该在一个main函数中。
• 一旦发生这种情况，您就可以对这些变量名进行去大写。
• 应该将DEAD_OR_ALIVE放在with块中

使用枚举

这是：

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

应该是

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