我正在做一个与引力透镜有关的项目,我需要计算长度为10^8个复点a = 1+0.48j
和b = 1
的阵列z的汇合超几何函数1F1(a,b,z)。我正在寻找一种有效的方法来计算大阵列大小的函数。scipy implementation速度快,但不接受a和b的复参数。
- mpmath**似乎是计算复杂参数的1F1的最佳方法,但是
mpmath.hyp1f1
不接受数组值。我找到的最佳解决方法是使用np.vectorize
或np.frompyfunc
来允许将NumPy数组作为参数传递。但是,这是非常慢的,需要几天的时间来执行(即使安装了gmpy2)。我想这是因为mpmath函数在大数组大小上总是很慢。
- mpmath**似乎是计算复杂参数的1F1的最佳方法,但是
一个非python的实现也可以,只要我能把结果保存在磁盘上,并把它读入我的python代码中。我见过一些可能工作的实现(例如https://www.math.ucla.edu/~mason/research/pearson_final.pdf),但我不确定。
另一种可能的方法是对函数进行插值(输入数组中的连续点非常接近),但我不确定什么是最好的方法。
谢谢!
1条答案
按热度按时间wsxa1bj11#
我也有和你一样的问题。
我发现mpmath包有一个“隐藏”的函数集(仅限)浮点精度,可以通过写入fp. upfront来访问。hyp1f1中不存在浮点精度,但更一般的hyper中不存在浮点精度。这意味着mpmath包中有一个fp.hyper,它与fp.hyper一起使用([a],[B],z)等价于hyper1f1(a,b,z),但速度快得多。如果使用np.vectorize对其进行矢量化,这将大大加快计算速度。
免责声明:我收到一条错误消息,说在计算这个值时,一些复数值通过去掉虚部转换为真实的,但到目前为止,我得到的结果似乎是合理的,并且与hyper1f1(a,b,z)值兼容。
补充:似乎fp.hyper不喜欢得到numpy数据类型,即使它们是标量,例如在a,B,z是numpy标量的情况下(例如numpy数组的一个元素),它只会返回1,而不会给出与实际输入无关的错误消息。
任何一种方式:在自己的risc使用。