在 python 中矢量化蒙特卡罗模拟答案

【问题标题】：Vectorizing a monte carlo simulation in python在 python 中矢量化蒙特卡罗模拟
【发布时间】：2019-01-21 21:07:45
【问题描述】：

我最近一直在编写 Python 中的一些代码，以使用蒙特卡罗方法模拟二维 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]

在程序开始时，我调用一个名为“randomize”的例程来生成 K 个具有小虚部的随机酉复数，并将它们存储在一个名为 Cnum、长度为 K 的数组中。在同一个例程中，我还通过我的链接数组并将每个元素设置为随机单位复数。代码如下。

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

以下是迭代例程的用途示例。我已经编写了一个名为斑块的例程，它计算给定 B0 的平均斑块（链接变量的 1 乘 1 闭环的实部）。迭代例程用于生成独立于先前配置的新字段配置。在我们获得新的现场配置后，我们计算所述配置的模板。然后我们使用 while 循环重复这个过程 j1 次，最后我们得到平均斑块。

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)

此代码有效，但速度非常慢。有没有办法可以使用多处理或其他东西来加快速度？

【问题讨论】：

函数内的代码请适当缩进，方便他人复制
正确缩进。然后向我们展示如何调用此例程的示例，其中包含代码或文件中的完整数据。阅读并关注How to create a Minimal, Complete, and Verifiable example。
对 Bazingaa 和 Rory；我很抱歉。我修复了缩进，并阅读了提供的链接。我将很快编辑帖子的其余部分，以举例说明如何调用例程。 @RoryDaulton，您是否希望我发布在他们自己的机器上运行代码所需的所有代码？
能否请您提供一段可以运行的代码？
是的，显示足够多的代码，以便有人可以复制并粘贴到编辑器中然后运行代码。它不需要是你所有的代码，足以显示你的问题。

标签： python montecarlo

【解决方案1】：

您应该使用cython 和c 数据类型。 Another cython link。它专为快速计算而设计。

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

如果您的预测没有考虑其他元素的状态，即它只关心一个元素，那么您可以将您的 Link 数组分成段并将块分成进程池中的多个进程，以及何时done 将这些段组合回一个数组。这肯定会节省时间，而且您不必使用任何额外的多处理机制。

【讨论】：