在numpy阵列上进行集群是一项繁重的任务

Question

我想通过将一个任务划分为不同的线程来加快python/numpy代码的速度（我的笔记本电脑上有四个线程，标准python使用一个线程，大家都知道）。你知道吗

我认为我的程序所做的可以按需要分成多个线程，请告诉我这是否正确：

1. 我们有两张图片A, B（两个大的numpy arrays，形状是（x，y，d））
2. 我们把A乘以一个变量P，然后把结果A'和B相加，得到C = np.add(A*P, B)
3. 我们在C上做了一个测试

这是我想分成四个线程的循环的基本描述。我正在考虑将numpy数组A, B分成四个数组A1, A2, A3, A4, B1, B2, B3, B4,，然后分别执行相同的操作。你知道吗

您能否建议一些起点来学习如何做到这一点（网站、python库等）

谢谢

代码如下：

import numpy as np
import Synth, time


width, height = 1644, 1644

colors_filename = "D65_maxVal16_4096_9900_LYVE_D+50.csv"
maxVal = 15

RGB, XYZ, sRGB = Synth.load_RGB_XYZ_RGB2(colors_filename)
SITE = False
tri, indices, simplices = Synth.delaunay(XYZ, SITE)

wave_start, wave_end, steps = 400, 700, 31
CMF = Synth.CMF_1931_2deg(wave_start, wave_end, steps)
Illuminant = Synth.Illuminant_D65(wave_start, wave_end, steps)
Neutral = Synth.Neutral(CMF, Illuminant)

background = Synth.black_bg(CMF, Illuminant, tri, width, height)
layer = np.zeros_like(background)

scales = 3
stars = 1
radius, power = width*5., 2. 

alive = 0

times = 10

for i in xrange(scales):

    stars *= times
    times *= 1.5
    radius /= 10.

    max_z = Synth.half_power_light_spot_dist(radius, power)
    min_z = max_z / 10.
    z_range = max_z-min_z

    print "scale", i, "max / min", max_z, min_z

    # the loop I described starts here
    for j in xrange(int(stars)):

        x, y = np.random.randint(0, width), np.random.randint(0, height)

        linear_progression = float(j)/float(stars)

        z = z_range * (1- linear_progression)

        L = 60

        C = np.random.randint(0,150)
        H = np.random.uniform(0., 360.)

        cieX, cieY, cieZ = Synth.LCH_to_XYZ(L, C, H, Neutral)
        cieXYZ = cieX/cieY, 1., cieZ/cieY

        power_star = np.random.uniform(2, 10)

        star = Synth.Stars_Light(width, height, x, y, z, power_star)

        print "star", j, "on", int(stars),

        precision = 0.01 

        star_copy = star

        star_copy[star < precision] = 0

        star_XYZ = np.multiply(star_copy, cieXYZ)

        # this is the slowest part, what I called the performed "test" between the two pictures A (here "background"), and B (here "star_XYZ"), see below to see what the function "layer_opt_mask" does.  
        P = Synth.layer_opt_mask(background[star_XYZ != 0], layer[star_XYZ != 0], star_XYZ[star_XYZ != 0], tri)

        print "power", P

        print "x, y, z", x, y, z
        print "L, C, H", L, C, H

        if P > 0.:

            non_linear_progression = Synth.L_to_Y(100 * (1-linear_progression))/100. #before it was not (1-linear_progression), but (linear_progression)

            print "LI, NLI", linear_progression, non_linear_progression

            """
            on pourrait ameliorer tout ce systeme la en effacant les donnees sous la limite apres avoir multiplie, ce serait un poil plus long mais ce serait plus logique
            """

            precision = 0.0001

            star[star < precision] = 0

            star = np.multiply(star, cieXYZ)

            layer[star != 0] = np.add(layer[star != 0], P * star[star != 0] * non_linear_progression * 0.5)

            alive += 1

            LEN = len(layer[star != 0])

            print alive, "ALIVE", "mask size:", LEN


        else:

            print "DEAD"

B, P = Synth.background_opt(background, layer, tri)

print "power Back, Front", B, P

image = np.add(B*background, P*layer)
image = np.reshape(image, (width*height, 3))

fname = "scales_"+str(scales)+"_alive_"+str(alive)+"_"+str(start_time)+"_"+str(version)+".pbm"

Synth.make_pbm_Maxval(tri, image, RGB, "RGB_"+fname, height, width, maxVal)
Synth.make_pbm(tri, image, sRGB, "sRGB_"+fname, height, width)

print fname

最慢的部分是“layer\u opt\u mask”函数，我认为它是集群的：

def layer_opt_mask(background, layer, light, tri):
    """ this function is maximizing the intensity of a new layer "light" on an existing 
background "background" where there's already a layer "layer". The limits of acceptable 
pixels intensities are defined by a set of colorimetric coordinates that form 
a delaunay-tetrahedralized convex hull "tri", the function tri.find_simplex returns -1 
when out of gamut pixels are present.
The objective is to have all the pixels values of the layer "light" inside the gamut 
with a maximal intensity """

    intensity = 0.
    nex = 100.

    #25 iterations might be a lot, it could be less, true
    it = 25

    while it > 0:

        it -= 1
        intensity += nex

        mult = np.multiply(light, intensity)

        test = np.add(mult, layer)
        test = np.add(background, test)

        sh = test.shape[0] / 3

        test = np.reshape(test, (int(sh), 3))

        # This is the slowest part but it is unfortunately absolutely necessary
        tet = np.min(tri.find_simplex(test))


        if tet == -1 and it > 0:
            intensity -= nex
            nex /= 10.

    if intensity > 0.01:
        return intensity
    else:
        return 0