上下文

为了更好地理解x射线衍射,我用python语言编写了代码.对于位置为R_i的点的集合,德拜公式为

enter image description here

其中指数中的i是复数,所有其他i是指数,为了简单起见,现在是b_i = b_j = 1.

Now I tried explicitly calculating this sum for a collection of points of which I have the coordinates enter image description here

import numpy as np
# set up grid
dims = 2
side = 30
points = np.power(side, dims)
coords = np.zeros((dims, points)) 

xc, yc = np.meshgrid(np.arange(side), np.arange(side))
coords[0, :] = xc.reshape((points))
coords[1, :] = yc.reshape((points))


# calculate diffraction
xdist = np.subtract.outer(coords[0], coords[0])
ydist = np.subtract.outer(coords[1], coords[1])
rdist = np.stack((xdist, ydist))
rdist = rdist.reshape(2, rdist.shape[1]*rdist.shape[2])

qs = 200
qspace = np.stack((np.linspace(-2, 8, qs), np.zeros(qs)))
diffrac = np.sum(np.exp(-1j * np.tensordot(qspace.T, rdist, axes=1)), axis=1)

几秒钟后我得到了以下信息

enter image description here

这看起来和预期的一样(周期为2 π,因为点的间距为1).这需要一些时间也是有道理的:对于900个点,必须计算810000个距离.我不使用循环,所以我认为代码在效率方面并不差,但我手动计算这个总和的事实似乎固有地慢.

思想

现在看起来,如果我能用一个离散的快速傅立叶变换来实现这个目标——给定和的形状,事情会大大加快.然而:

  • 对于离散傅里叶变换,我仍然需要对图像进行像素化(据我所知),以便在我的信号点之间包含大量的空白空间.就像我要转换我分享的第一张图片的像素一样.这似乎也不太有效(例如,因为采样).
  • 我想在之后移动点,所以第一个图像是一个网格,因此定期采样的事实并没有特别的帮助.看起来好像非均匀傅立叶变换可以帮助我,但仍然需要我对图像进行"像素化",并将某些值设置为0.

问题

有没有一种方法可以使用FFT(或其他方法)更快地计算和,从np.数组坐标(x,y)的列表开始?(狄拉克德尔塔函数,如果你想...).

特别是相关数学技术/Python函数/Python包的指针将受到赞赏.我对实际应用中使用傅立叶变换并不熟悉,但我在网上找到的大多数material 似乎都无关紧要.所以可能我看错了方向,或者我缺乏理解.所有的帮助是赞赏!

(第一张图片是https://www.ill.eu/fileadmin/user_upload/ILL/6_Careers/1_All_our_vacancies/PhD_recruitment/Student_Seminars/2017/19-2017-05-09_Fischer_Cookies.pdf的截图,因为似乎SO上没有数学符号,或者我没有找到它)

推荐答案

这个答案提供了一个解决方案,使代码更高效,从而充分利用CPU的计算能力,从而使代码更快.


超过90%的时间都花在np.exp上,因为计算复数的经验是非常昂贵的.

一个加快速度的解决方案是使用multiple threads(因为Numpy不使用多线程).除此之外,我们还可以使用faster implementation of 100(通常利用CPU的SIMD单元).这两个都可以用Numexpr轻松完成.

然后,我们可以使用矩阵乘法qspace.T @ rdist来加速np.tensordot运算,因为块执行是低效的.

import numexpr as ne

# Equivalent of the last line of the code:
tmp1 = qspace.T @ rdist
tmp2 = ne.evaluate('exp(-1j * tmp1)')
diffrac = np.sum(tmp2, axis=1)

绩效评价

以下是i5—9600KF CPU(6核)的性能结果:

Initial code:         9.3 s
New proposed code:    1.1 s

因此,新的实现是8.5 times faster.大部分时间仍然花在计算复数的指数上(60%).>

Python相关问答推荐

根据条件将新值添加到下面的行或下面新创建的行中

acme错误-Veritas错误:模块收件箱没有属性linear_util'

使用FASTCGI在IIS上运行Django频道

Python解析整数格式说明符的规则?

在Python中动态计算范围

如何使用Python以编程方式判断和检索Angular网站的动态内容?

Pandas计数符合某些条件的特定列的数量

我想一列Panadas的Rashrame,这是一个URL,我保存为CSV,可以直接点击

如何使用scipy的curve_fit与约束,其中拟合的曲线总是在观测值之下?

删除marplotlib条形图上的底边

SQLAlchemy bindparam在mssql上失败(但在mysql上工作)

无法连接到Keycloat服务器

我的字符串搜索算法的平均时间复杂度和最坏时间复杂度是多少?

可以bcrypts AES—256 GCM加密损坏ZIP文件吗?

如何使regex代码只适用于空的目标单元格

跳过嵌套JSON中的级别并转换为Pandas Rame

如何在Python中使用Iscolc迭代器实现观察者模式?

如何获得3D点的平移和旋转,给定的点已经旋转?

如何将返回引用的函数与pybind11绑定?

PyTorch变压器编码器中的填充掩码问题