在二维 numpy 数组上广播函数答案

【问题标题】：broadcasting a function on a 2-dimensional numpy array在二维 numpy 数组上广播函数
【发布时间】：2015-05-09 14:16:54
【问题描述】：

我想通过在 numpy 数组而不是 for 循环 上计算一次函数来提高我的代码速度，而不是 thispython 的函数图书馆。如果我有如下功能：

import numpy as np
import galsim
from math import *
M200=1e14
conc=6.9
def func(M200, conc):
    halo_z=0.2
    halo_pos =[1200., 3769.7]
    halo_pos = galsim.PositionD(x=halo_pos_arcsec[0],y=halo_pos_arcsec[1])
    nfw      = galsim.NFWHalo(mass=M200, conc=conc, redshift=halo_z,halo_pos=halo_pos, omega_m = 0.3, omega_lam =0.7)
    for i in range(len(shear_z)):
       shear_pos=galsim.PositionD(x=pos_arcsec[i,0],y=pos_arcsec[i,1])
       model_g1, model_g2  = nfw.getShear(pos=self.shear_pos, z_s=shear_z[i])
    l=np.sum(model_g1-model_g2)/sqrt(np.pi)
    return l

虽然pos_arcsec 是24000x2 的二维数组，而shear_z 也是带有24000 元素的一维数组。 主要问题是我想在M200=np.arange(13., 16., 0.01) 和conc = np.arange(3, 10, 0.01) 的网格上计算这个函数。我不知道如何通过M200 和conc 广播这个函数来估计这个二维数组。运行代码需要很多时间。我正在寻找加速这些计算的最佳方法。

【问题讨论】：

广播这个功能是什么意思？
使用 for 循环会是什么样子？我假设 0 和 1 是变量。
for 循环遍历的是什么？分享for循环代码？
@Divakar 我更新了我的问题！
@Tichodroma 我解释了为什么我需要在新更新中进行广播！

标签： python numpy vectorization

【解决方案1】：

当pos 是一个形状为(n,2) 的数组时，这应该可以工作

import numpy as np

def f(pos, z):
    r=np.sqrt(pos[...,0]**2+pos[...,1]**2)
    return np.log(r)*(z+1)

例子：

z = np.arange(10)
pos = np.arange(20).reshape(10,2)

f(pos,z)
# array([  0.        ,   2.56494936,   5.5703581 ,   8.88530251,
#         12.44183436,  16.1944881 ,  20.11171117,  24.17053133,
#         28.35353608,  32.64709419])

【讨论】：

【解决方案2】：

使用numpy.linalg.norm

如果你有一个数组：

import numpy as np
import numpy.linalg as la

a = np.array([[3, 4], [5, 12], [7, 24]])

然后您可以通过

确定结果向量的大小 (sqrt(a^2 + b^2))

b = np.sqrt(la.norm(a, axis=1)

>>> print b
array([  5.,  15.  25.])

【讨论】：