PROSAGA码农传奇-雾计算/网格计算-使用大数组使用NumPy在3D中高效的高斯总和

<div class =“post-text”itemprop =“text”>
  <P>
    考虑到你的许多距离计算在指数之后将给出零重量，你可能会减少很多距离。在丢弃大于阈值的距离的同时进行大块距离计算通常会更快
     <code>
      KDTree
    </code>
    ：
  </p>
   <pre>
    <code>
      import numpy as np
from scipy.spatial import cKDTree # so we can get a `coo_matrix` output

def gaussgrid(coords, sigma = 1, n = 64, side = 100, eps = None):
    x_ = np.linspace(-side/2,side/2,n)
    x,y,z = np.meshgrid(x_,x_,x_,indexing='ij')
    xyz = np.column_stack((x.ravel(),y.ravel(),z.ravel()))
    if eps is None:
        eps = np.finfo('float64').eps
    thr = -np.log(eps) * 2 * sigma**2
    data_tree = cKDTree(coords)
    discr = 1000 # you can tweak this to get best results on your system
    values = np.empty(n**3)
    for i in range(n**3//discr + 1):
        slc = slice(i * discr, i * discr + discr)
        grid_tree = cKDTree(xyz[slc])
        dists = grid_tree.sparse_distance_matrix(data_tree, thr, output_type = 'coo_matrix')
        dists.data = 1./np.sqrt(2*np.pi*sigma**2) * np.exp(-dists.data/(2*sigma**2))
        values[slc] = dists.sum(1).squeeze()
    return values.reshape(n,n,n)

</code>
  </pre>
  <P>
    现在，即使你保持
     <code>
      eps = None
    </code>
     当你仍然返回大约10％的距离时，它会快一点，但是当eps = 1e-6左右时，你应该得到一个很大的加速。在我的系统上：
  </p>
   <pre>
    <code>
      %timeit out = sumofgauss(coords, xyz, 1.0)
1 loop, best of 3: 23.7 s per loop

%timeit out = gaussgrid(coords)
1 loop, best of 3: 2.12 s per loop

%timeit out = gaussgrid(coords, eps = 1e-6)
1 loop, best of 3: 382 ms per loop

</code>
  </pre>
</DIV>