增长一个numpy数组答案

【问题标题】：Growing a numpy array增长一个numpy数组
【发布时间】：2018-01-24 10:56:48
【问题描述】：

我有一个 3x256 规则的规则集。每个规则都映射到一个 3x3 的值网格，而这些值本身就是规则。

示例规则：

0 -> [[0,0,0],[0,1,0],[0,0,0]]  
1 -> [[1,1,1],[0,0,0],[1,1,1]]

种子：

[[0]]

1 次迭代后：

[[0,0,0],
[0,1,0],
[0,0,0]]

2 次迭代后：

[[0,0,0,0,0,0,0,0,0],
[0,1,0,0,1,0,0,1,0],
[0,0,0,0,0,0,0,0,0],
[0,0,0,1,1,1,0,0,0],
[0,1,0,0,0,0,0,1,0],
[0,0,0,1,1,1,0,0,0],
[0,0,0,0,0,0,0,0,0],
[0,1,0,0,1,0,0,1,0],
[0,0,0,0,0,0,0,0,0]]

现在我有了一个有效的实现，但是，它是我脚本中最慢的函数。我想知道是否有更pythonic和更有效的方法来重写这个函数。

def decode(rules,fractal_iterations,seed):
    final_seed_matrix = np.zeros((3,3**fractal_iterations,3**fractal_iterations))

    for i in range(dimensions):
        seed_matrix = np.array([[seed]])
        for j in range(fractal_iterations):
            size_y = seed_matrix.shape[0]
            size_x = seed_matrix.shape[1]
            new_matrix = np.zeros((size_y*rule_size_sqrt,size_x*rule_size_sqrt))

            for y in range(size_y):
                for x in range(size_x):
                    seed_value = seed_matrix[y,x]
                    new_matrix[y*rule_size_sqrt : y*rule_size_sqrt+rule_size_sqrt, x*rule_size_sqrt : x*rule_size_sqrt+rule_size_sqrt] = rules[int(seed_value),i]

            seed_matrix = new_matrix
        final_seed_matrix[i] = seed_matrix

    return np.moveaxis(final_seed_matrix,0,-1)

【问题讨论】：

如果您不坚持使用纯 Python 而是使用近乎 Pythonic 的代码，也许可以使用 Cython。我不保证你可以使用 Cython 获得比其他 Python 算法更快的代码。我不是 Cython 的专家。 cython.readthedocs.io/en/latest/src/tutorial/numpy.html
一般规则是在追加速度很快的地方使用 python 数组。然后在收集完所有内容后，创建一次数组。
预分配数组很好。嵌套的for 循环是性能杀手，但我不确定是否可以在这里进行大量矢量化。 dimensions 变量和相关外循环的目的是什么？
@kazemakase 维度变量用于 3 个颜色维度。它是全局定义的，但在第 4 行应该是 3。
这是否意味着您独立地迭代每个维度？如果移动到最内层循环，这可能可能被矢量化。

标签： python numpy fractals

【解决方案1】：

这是一个优化版本，它使用高级索引在一个索引步骤中选择和修补所有规则。这将在它替换的像素位置创建一个具有适当规则的 4D 数组。然后将其展平为 2D 是交换中间轴和重塑的问题。它似乎给出了与您相同的结果，但速度明显更快（目前仅针对整数规则进行了测试）：

results equal: True
OP       : 24.883304461836815 ms
optimized: 1.093490980565548 ms

代码：

import numpy as np

dimensions = 3
rule_size_sqrt = 3

def decode(rules,fractal_iterations,seed):
    final_seed_matrix = np.zeros((3,3**fractal_iterations,3**fractal_iterations))

    for i in range(dimensions):
        seed_matrix = np.array([[seed]])
        for j in range(fractal_iterations):
            size_y = seed_matrix.shape[0]
            size_x = seed_matrix.shape[1]
            new_matrix = np.zeros((size_y*rule_size_sqrt,size_x*rule_size_sqrt))

            for y in range(size_y):
                for x in range(size_x):
                    seed_value = seed_matrix[y,x]
                    new_matrix[y*rule_size_sqrt : y*rule_size_sqrt+rule_size_sqrt, x*rule_size_sqrt : x*rule_size_sqrt+rule_size_sqrt] = rules[int(seed_value),i]

            seed_matrix = new_matrix
        final_seed_matrix[i] = seed_matrix

    return np.moveaxis(final_seed_matrix,0,-1)

def decode_fast(rules, fractal_iterations, seed):
    rules_int = rules.astype(int)
    seed = np.array([[seed]], dtype=int)
    res = np.empty((3**fractal_iterations, 3**fractal_iterations, dimensions),
                   dtype=rules.dtype)
    for i in range(dimensions):
        grow = seed
        for j in range(1, fractal_iterations):
            grow = rules_int[grow, i].swapaxes(1, 2).reshape(3**j, -1)
        grow = rules[grow, i].swapaxes(1, 2).reshape(3**fractal_iterations, -1)
        res[..., i] = grow
    return res

rules = np.random.randint(0, 4, (4, dimensions, 3, 3))
seed = 1
fractal_iterations = 5
print('results equal:', np.all(decode(rules, fractal_iterations, seed) == decode_fast(rules, fractal_iterations, seed)))

from timeit import repeat
print('OP       :', min(repeat('decode(rules, fractal_iterations, seed)', globals=globals(), number=50))*20, 'ms')
print('optimized:', min(repeat('decode_fast(rules, fractal_iterations, seed)', globals=globals(), number=50))*20, 'ms')

【讨论】：