Векторизация симуляции Монте-Карло на Python

Недавно я работал над некоторым кодом на Python, чтобы моделировать двумерную калибровочную теорию U (1) с использованием методов Монте-Карло. По сути, у меня есть массив n x n на 2 (назовите его Link) унитарных комплексных чисел (их величина равна единице). Я случайным образом выбираю элемент своего массива ссылок и предлагаю случайное изменение числа на этом сайте. Затем я вычисляю результирующее изменение действия, которое могло бы произойти из-за этого изменения. Затем я принимаю изменение с вероятностью, равной min (1, exp (-dS)), где dS - это изменение действия. Код итератора следующий

def iteration(j1,B0):
    global Link
    Staple = np.zeros((2),dtype=complex)
    for i0 in range(0,j1):
        x1 = np.random.randint(0,n)
        y1 = np.random.randint(0,n)
        u1 = np.random.randint(0,1)
        Linkrxp1 = np.roll(Link,-1, axis = 0)
        Linkrxn1 = np.roll(Link, 1, axis = 0)
        Linkrtp1 = np.roll(Link, -1, axis = 1)
        Linkrtn1 = np.roll(Link, 1, axis = 1)
        Linkrxp1tn1 = np.roll(np.roll(Link, -1, axis = 0),1, axis = 1)
        Linkrxn1tp1 = np.roll(np.roll(Link, 1, axis = 0),-1, axis = 1)


        Staple[0] = Linkrxp1[x1,y1,1]*Linkrtp1[x1,y1,0].conj()*Link[x1,y1,1].conj() + Linkrxp1tn1[x1,y1,1].conj()*Linkrtn1[x1,y1,0].conj()*Linkrtn1[x1,y1,1]
        Staple[1] = Linkrtp1[x1,y1,0]*Linkrxp1[x1,y1,1].conj()*Link[x1,y1,0].conj() + Linkrxn1tp1[x1,y1,0].conj()*Linkrxn1[x1,y1,1].conj()*Linkrxn1[x1,y1,0]


        uni = unitary()
        Linkprop = uni*Link[x1,y1,u1]
        dE3 = (Linkprop - Link[x1,y1,u1])*Staple[u1]
        dE1 = B0*np.real(dE3)
        d1 = np.random.binomial(1, np.minimum(np.exp(dE1),1))


        d = np.random.uniform(low=0,high=1)
        if d1 >= d:
            Link[x1,y1,u1] = Linkprop
        else:
            Link[x1,y1,u1] = Link[x1,y1,u1]

В начале программы я вызываю процедуру под названием "randomize" для генерации K случайных унитарных комплексных чисел, которые имеют небольшие мнимые части, и сохраняю их в массиве под названием Cnum длины K. В той же процедуре я также просматриваю свой массив Link и устанавливаю каждый элемент к случайному унитарному комплексному числу. Код указан ниже.

def randommatrix():
    global Cnum
    global Link
    for i1 in range(0,K):
        C1 = np.random.normal(0,1)
        Cnum[i1] = np.cos(C1) + 1j*np.sin(C1)
        Cnum[i1+K] = np.cos(C1) - 1j*np.sin(C1)
    for i3,i4 in itertools.product(range(0,n),range(0,n)):
        C2 = np.random.uniform(low=0, high = 2*np.pi)
        Link[i3,i4,0] = np.cos(C2) + 1j*np.sin(C2)
        C2 = np.random.uniform(low=0, high = 2*np.pi)
        Link[i3,i4,1] = np.cos(C2) + 1j*np.sin(C2)

Следующая процедура используется во время итерационной процедуры для получения случайного комплексного числа с небольшой мнимой частью (путем извлечения случайного элемента массива Cnum, который мы сгенерировали ранее).

def unitary():
    I1 = np.random.randint((0),(2*K-1))
    mat = Cnum[I1]
    return mat

Вот пример того, для чего будет использоваться процедура итерации. Я написал процедуру, называемую плакеткой, которая вычисляет средний плакат (действительная часть замкнутого цикла переменных связи 1: 1) для заданного B0. Процедура итерации используется для создания новых конфигураций полей, которые не зависят от предыдущих конфигураций. После получения новой конфигурации поля мы рассчитываем табличку для указанной конфигурации. Затем мы повторяем этот процесс j1 раз, используя цикл while, и в конце получаем табличку со средним значением.

def Plq(j1,B0):
    i5 = 0
    Lboot = np.zeros(j1)
    while i5<j1:
        iteration(25000,B0)
        Linkrxp1 = np.roll(Link,-1, axis = 0)
        Linkrtp1 = np.roll(Link, -1, axis = 1)
        c0 = np.real(Link[:,:,0]*Linkrxp1[:,:,1]*Linkrtp1[:,:,0].conj()*Link[:,:,1].conj())
        i5 = i5 + 1

Нам нужно определить некоторые переменные, прежде чем что-либо запускать, поэтому вот начальные переменные, которые я определяю перед определением любых подпрограмм.

K = 20000
n = 50
a = 1.0
Link = np.zeros((n,n,2),dtype = complex)
Cnum = np.zeros((2*K), dtype = complex)

Этот код работает, но он очень медленный. Есть ли способ использовать многопроцессорность или что-то еще, чтобы ускорить это?