微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦！

Python中积分循环的有效计算

分类：编程问答

如何解决Python中积分循环的有效计算

我想知道如何加速以下代码，其中我计算了一个涉及数值积分的概率函数，然后我计算了一些置信区间。

我想到的一些可能性是 Numba 或代码的矢量化

编辑： 我做了一些小的修改，因为有一个错误。我正在寻找一些提供重大时间改进的修改（我知道有一些细微的变化会提供一些小的时间改进，例如重复功能，但我并不关心它们） 代码是：

# -*- coding: utf-8 -*-
"""
Created on Tue Jan 26 17:05:46 2021

@author: Ignacio
"""

import numpy as np
from scipy.integrate import simps

def pdf(V,alfa_points):
    alfa=np.linspace(0,2*np.pi,alfa_points)
    return simps(1/np.sqrt(2*np.pi)/np.sqrt(sigma_R2)*np.exp(-(V*np.cos(alfa)-eR)**2/2/sigma_R2)*1/np.sqrt(2*np.pi)/np.sqrt(sigma_I2)*np.exp(-(V*np.sin(alfa)-eI)**2/2/sigma_I2),alfa)
    
def find_nearest(array,value):
    array=np.asarray(array)
    idx = (np.abs(array-value)).argmin()
    return array[idx]

N = 20
n=np.linspace(0,N-1,N)
d=1
sigma_An=0.1
sigma_Pn=0.2

An=np.ones(N)
Pn=np.zeros(N)

Vs=np.linspace(0,30,1000)
inc=np.max(Vs)/len(Vs)
th=np.linspace(0,np.pi/2,250)
R=np.sum(An*np.cos(Pn+2*np.pi*np.sin(th[:,np.newaxis])*n*d),axis=1)
I=np.sum(An*np.sin(Pn+2*np.pi*np.sin(th[:,axis=1)

fmin=np.zeros(len(th))
fmax=np.zeros(len(th))
for tt in range(len(th)):
    eR=np.exp(-sigma_Pn**2/2)*np.sum(An*np.cos(Pn+2*np.pi*np.sin(th[tt])*n*d))
    eI=np.exp(-sigma_Pn**2/2)*np.sum(An*np.sin(Pn+2*np.pi*np.sin(th[tt])*n*d)) 
    
    sigma_R2=1/2*np.sum(An*sigma_An**2)+1/2*(1-np.exp(-sigma_Pn**2))*np.sum(An**2)+1/2*np.sum(np.cos(2*(Pn+2*np.pi*np.sin(th[tt])*n*d))*((An**2+sigma_An**2)*np.exp(-2*sigma_Pn**2)-An**2*np.exp(-sigma_Pn**2)))
    sigma_I2=1/2*np.sum(An*sigma_An**2)+1/2*(1-np.exp(-sigma_Pn**2))*np.sum(An**2)-1/2*np.sum(np.cos(2*(Pn+2*np.pi*np.sin(th[tt])*n*d))*((An**2+sigma_An**2)*np.exp(-2*sigma_Pn**2)-An**2*np.exp(-sigma_Pn**2))) 
    
    PDF=np.zeros(len(Vs))
    for vv in range(len(Vs)):
        PDF[vv]=pdf(Vs[vv],100)
    total=simps(PDF,Vs)
    values=np.cumsum(PDF)*inc/total
    xval_05=find_nearest(values,0.05)
    fmin[tt]=Vs[values==xval_05]
    xval_95=find_nearest(values,0.95)
    fmax[tt]=Vs[values==xval_95]

解决方法

此版本的加速：31 倍

一个简单的分析 (%%prun) 表明大部分时间都花在 simps 上。

您可以控制在 pdf() 中完成的积分：例如，如果您稍微增加 alpha 的分辨率，您可以使用梯形方法代替辛普森，数值差异可以忽略不计。实际上，alpha 的更高采样所获得的更高分辨率足以弥补 simps 和 trapeze 之间的差异（原因见底部图片）。这是迄今为止最高的加速。我们通过自己实现空中飞人方法而不是使用 scipy 来更进一步，因为它非常简单。仅此一项就产生了边际收益，但为更彻底的优化打开了大门（下面，大约 pdf2D。

此外，当它知道 simps(PDF,...) 步是常数时，剩余的 dx 会更快，所以我们可以这样说，而不是传递整个 alpha 数组。

您可以避免执行循环来计算 PDF 并直接在 np.vectorize(pdf) 上使用 Vs，或者更好（如下面的代码所示），执行该计算的二维版本.

还有一些其他的小事（比如直接使用索引 fmin[tt] = Vs[closest(values,0.05)] 而不是查找索引，返回值，然后使用布尔掩码作为 where values == xval_05），或者把所有的函数外的常量（包括 alpha），避免每次都重新计算。

以上为我们提供了 5.2 倍的改进。你的代码中有很多我不明白的地方，例如为什么有 An（一个）和 Pn（零）？

但是，重要的是，另一个约 6 倍的加速来自观察，因为我们正在使用 trapeze 原语实现我们自己的 numpy 方法，我们实际上可以一次性在 2D 中完成它整个PDF。

以下代码的最终加速是 31 倍。我相信，更好地了解您想要做什么的“大局”会带来额外的、也许是实质性的速度提升。

修改后的代码：

import numpy as np
from scipy.integrate import simps


alpha_points = 200  # more points as we'll use trapeze not simps
alpha = np.linspace(0,2*np.pi,alpha_points)
cosalpha = np.cos(alpha)
sinalpha = np.sin(alpha)
d_alpha = np.mean(np.diff(alpha))  # constant dx
coeff = 1 / np.sqrt(2*np.pi)
Vs=np.linspace(0,30,1000)
d_Vs = np.mean(np.diff(Vs))  # constant dx
inc=np.max(Vs)/len(Vs)

def f2D(Vs,eR,sigma_R2,eI,sigma_I2):
    a = coeff / np.sqrt(sigma_R2)
    b = coeff / np.sqrt(sigma_I2)
    y = a * np.exp(-(np.outer(cosalpha,Vs) - eR)**2 / 2 / sigma_R2) * b * np.exp(-(np.outer(sinalpha,Vs) - eI)**2 / 2 / sigma_I2)
    return y

def pdf2D(Vs,sigma_I2):
    y = f2D(Vs,sigma_I2)
    s = y.sum(axis=0) - (y[0] + y[-1]) / 2  # our own impl of trapeze,on 2D y
    return s * d_alpha

def closest(a,val):
    return np.abs(a - val).argmin()

N = 20
n = np.linspace(0,N-1,N)
d = 1
sigma_An = 0.1
sigma_Pn = 0.2

An=np.ones(N)
Pn=np.zeros(N)

th = np.linspace(0,np.pi/2,250)
R = np.sum(An*np.cos(Pn+2*np.pi*np.sin(th[:,np.newaxis])*n*d),axis=1)
I = np.sum(An*np.sin(Pn+2*np.pi*np.sin(th[:,axis=1)

fmin=np.zeros(len(th))
fmax=np.zeros(len(th))
for tt in range(len(th)):
    eR=np.exp(-sigma_Pn**2/2)*np.sum(An*np.cos(Pn+2*np.pi*np.sin(th[tt])*n*d))
    eI=np.exp(-sigma_Pn**2/2)*np.sum(An*np.sin(Pn+2*np.pi*np.sin(th[tt])*n*d)) 
    
    sigma_R2=1/2*np.sum(An*sigma_An**2)+1/2*(1-np.exp(-sigma_Pn**2))*np.sum(An**2)+1/2*np.sum(np.cos(2*(Pn+2*np.pi*np.sin(th[tt])*n*d))*((An**2+sigma_An**2)*np.exp(-2*sigma_Pn**2)-An**2*np.exp(-sigma_Pn**2)))
    sigma_I2=1/2*np.sum(An*sigma_An**2)+1/2*(1-np.exp(-sigma_Pn**2))*np.sum(An**2)-1/2*np.sum(np.cos(2*(Pn+2*np.pi*np.sin(th[tt])*n*d))*((An**2+sigma_An**2)*np.exp(-2*sigma_Pn**2)-An**2*np.exp(-sigma_Pn**2))) 
    
    PDF=pdf2D(Vs,sigma_I2)
    total = simps(PDF,dx=d_Vs)
    values = np.cumsum(PDF) * inc / total
    fmin[tt] = Vs[closest(values,0.05)]
    fmax[tt] = Vs[closest(values,0.95)]

注：大部分的fmin和fmax与原函数相比都是np.allclose()，但也有一些有小错误：经过一番挖掘，发现实现这里更精确，因为函数 f() 可能非常突然，更多的 alpha 点实际上有帮助（并且不仅可以弥补由于使用空中飞人而不是辛普森而造成的微小精度不足）。

例如，在索引 tt=244,vv=400 处：

enter image description here

,

考虑多种方法，提供最大时间改进的方法是 Numba 方法。皮埃尔提出的方法很有意思，不需要安装其他包，是一种资产。

但是，在我计算的示例中，时间改进不如 numba 示例大，特别是当 th 中的点增长到十分之几时（这是我的实际情况）。我在这里发布 Numba 代码，以防有人感兴趣：

import numpy as np
from numba import njit

@njit
def margins(val_min,val_max):
    fmin=np.zeros(len(th))
    fmax=np.zeros(len(th))
    for tt in range(len(th)):    
        eR=np.exp(-sigma_Pn**2/2)*np.sum(An*np.cos(Pn+2*np.pi*np.sin(th[tt])*n*d))
        eI=np.exp(-sigma_Pn**2/2)*np.sum(An*np.sin(Pn+2*np.pi*np.sin(th[tt])*n*d))
    
        sigma_R2=1/2*np.sum(An*sigma_An**2)+1/2*(1-np.exp(-sigma_Pn**2))*np.sum(An**2)+1/2*np.sum(np.cos(2*(Pn+2*np.pi*np.sin(th[tt])*n*d))*((An**2+sigma_An**2)*np.exp(-2*sigma_Pn**2)-An**2*np.exp(-sigma_Pn**2)))
        sigma_I2=1/2*np.sum(An*sigma_An**2)+1/2*(1-np.exp(-sigma_Pn**2))*np.sum(An**2)-1/2*np.sum(np.cos(2*(Pn+2*np.pi*np.sin(th[tt])*n*d))*((An**2+sigma_An**2)*np.exp(-2*sigma_Pn**2)-An**2*np.exp(-sigma_Pn**2)))
    
        Vs=np.linspace(0,1000)
        inc=np.max(Vs)/len(Vs)
        integration_points=200
        PDF=np.zeros(len(Vs))
        for vv in range(len(Vs)):
            PDF[vv]=np.trapz(1/np.sqrt(2*np.pi)/np.sqrt(sigma_R2)*np.exp(-(Vs[vv]*np.cos(np.linspace(0,integration_points))-eR)**2/2/sigma_R2)*1/np.sqrt(2*np.pi)/np.sqrt(sigma_I2)*np.exp(-(Vs[vv]*np.sin(np.linspace(0,integration_points))-eI)**2/2/sigma_I2),np.linspace(0,integration_points))
              
        total=np.trapz(PDF,Vs)
        values=np.cumsum(PDF)*inc/total
        idx = (np.abs(values-val_min)).argmin()
        xval_05=values[idx]            
        fmin[tt]=Vs[np.where(values==xval_05)[0][0]]
        idx = (np.abs(values-val_max)).argmin()   
        xval_95=values[idx]
        fmax[tt]=Vs[np.where(values==xval_95)[0][0]]

    return fmin,fmax


N = 20
n=np.linspace(0,N)
d=1
sigma_An=1/2**6
sigma_Pn=2*np.pi/2**6
An=np.ones(N)
Pn=np.zeros(N)

th=np.linspace(0,250)

R=np.sum(An*np.cos(Pn+2*np.pi*np.sin(th[:,axis=1)
I=np.sum(An*np.sin(Pn+2*np.pi*np.sin(th[:,axis=1)

F=R+1j*I
Fa=np.abs(F)/np.max(np.abs(F))

fmin,fmax = margins(0.05,0.95)

版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。

苹果市值2025年有望达4万亿美元
pythonJavaScriptjavaHTMLPHPreactjsC#AndroidCSSNode.jssqlrpython-3.xMysqLjQueryc++pandasFlutterangularIOSdjangolinuxswifttypescript路由器JSON路由器设置无线路由器h3c华三华三路由器设置华三路由器电脑软件教程arraysdocker软件图文教程Cvue.jslaravelspring-boot