在python中向量化蒙特卡洛模拟

Question

我最近一直在使用python编写一些代码，以使用蒙特卡洛方法模拟二维U（1）规范理论。本质上，我有一个n×n×2的unit复数数组（称为链接）（其大小为1）。我随机选择我的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]

在程序开始时，我调用一个称为“随机化”的例程，以生成K个具有较小虚部的随机unit复数，并将它们存储在一个长度为K的Cnum数组中。在同一例程中，我还将遍历Link数组并进行设置每个元素为一个随机的complex复数。该代码在下面列出。

在迭代例程期间，使用以下例程来获取具有较小虚部的随机复数（通过检索我们先前生成的Cnum数组的随机元素）。

我们需要先定义一些变量，然后再运行任何东西，所以这是我在定义任何例程之前定义的初始变量

K = 20000

n = 50

a = 1.0

Link = np.zeros((n,n,2),dtype = complex)

Cnum = np.zeros((2*K), dtype = complex)

这段代码有效，但是速度很慢。有没有一种方法可以使用多处理或其他方法来加快速度？

Answer 1

您应该使用cython和c数据类型。另一个cython链接。它是为快速计算而构建的。您可以在以下两种情况之一中使用multiprocessing。如果您有一个需要多个进程共享的对象，则需要使用Manager（请参阅多处理链接），Lock和Array在进程之间共享该对象。但是，不能保证这样做会提高速度，因为每个过程都需要锁定链接以保证您的预测，假设预测受到链接中所有元素的影响（如果一个过程同时修改了另一个过程的一个元素）正在对元素进行预测，则该预测将不会基于最新信息。如果您的预测未考虑其他元素的状态，即它只关心一个元素，那么您可以将Link数组划分为多个部分，然后将块划分为一个过程池中的多个过程，并在完成后合并分段回到一个数组。这当然可以节省时间，并且您不必使用任何其他的多处理机制。