用python实现蒙特卡罗模拟的矢量化

Question

我最近一直在用python编写一些代码，用montecarlo方法模拟二维U（1）规范理论。本质上，我有一个n乘n乘2的酉复数数组（称之为Link）（它们的大小为1）。我随机选择我的链接数组的元素，并建议对该站点的数字进行随机更改。然后，我计算由于这种变化而导致的动作变化。然后我接受变化，概率等于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]

在程序开始时，我调用一个名为“randominize”的例程来生成K个具有较小虚部的随机酉复数，并将它们存储在一个名为Cnum的长度为K的数组中。在同一个例程中，我还遍历了我的Link数组，并将每个元素设置为一个随机酉复数。代码如下所示。在

^{pr2}$

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

^{3}$

下面是一个迭代例程将用于什么的示例。我编写了一个名为plaquette的例程，它计算给定B0的平均plaquette（链接变量的1×1闭环的实际部分）。迭代例程用于生成独立于先前配置的新字段配置。在我们得到一个新的字段配置之后，我们计算所述配置的plaquette。然后我们使用while循环重复这个过程j1次，最后得到mean plaquette。在

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)

这段代码可以工作，但速度非常慢。有没有一种方法可以让我使用多重处理或其他方法来加快速度？在

Answer 1

您应该使用cython和c数据类型。Another cython link。它是为快速计算而建的。在

您可以在两种情况下使用multiprocessing。 如果多个进程需要共享一个对象，则需要使用Manager（参见multiprocessing link）、Lock和Array在进程之间共享该对象。但是，不能保证这会提高速度，因为每个进程都需要锁定链接以保证您的预测，假设预测受到链接中所有元素的影响（如果一个进程同时修改了一个元素，而另一个进程正在对一个元素进行预测，则预测将不会基于最新信息）。在

如果你的预测没有考虑到其他元素的状态，也就是说，它只关心一个元素，那么你可以把你的链接数组分解成段，把块分割成一个进程池中的几个进程，完成后再把这些段合并成一个数组。这当然可以节省时间，而且您不必使用任何额外的多处理机制。在