【问题标题】:Multiply number of distances in distance matrix prior to histogram binning在直方图分箱之前乘以距离矩阵中的距离数
【发布时间】:2014-04-30 14:10:09
【问题描述】:

我正在使用 scipy.spatial.distance.pdist 来计算与坐标数组的距离,然后使用 numpy.histogram 来对结果进行分类。目前,这会将每个坐标视为一个对象存在,但是我在同一个坐标上有多个对象。一种选择是更改数组,以便每个坐标出现多次,对于该坐标处的每个对象一次,但是这将大大增加数组的大小和 pdist 的计算时间,因为它按 N^2 缩放,并且这是非常昂贵的,而且速度在这个应用程序中很重要。

第二种方法是处理得到的距离矩阵,使得每个距离重复 ninj 次,其中 ni 是坐标 i 处的对象数,nj 是坐标 j 处的对象数。这会将原始的 MxM 距离矩阵转换为 NxN 距离矩阵,其中 M 是数组中坐标的总数,但 N 是对象的总数。但同样,这似乎是不必要的昂贵,因为我真正需要做的就是以某种方式告诉直方图函数将距离 ij 处的事件数乘以 ninj。换句话说,有没有办法告诉 numpy.histogram 在距离 ij 处不只有一个对象,而是有 ni*nj 个对象?

显然欢迎其他想法。

编辑:

这是第一种方法的示例。

import numpy as np
from scipy import spatial
import matplotlib.pyplot as plt

#create array of 5 coordinates in 3D
coords = np.random.random(15).reshape(5,3)
'''array([[ 0.66500534,  0.10145476,  0.92528492],
       [ 0.52677892,  0.07756804,  0.50976737],
       [ 0.50030508,  0.37635556,  0.20828815],
       [ 0.02707651,  0.21878467,  0.55855427],
       [ 0.81564621,  0.82750694,  0.53083443]])'''

#number of objects at each coordinate
objects = np.random.randint(1,10,5)
#array([5, 3, 8, 5, 1])

#create new array with coordinates for each individual object
new_coords = np.zeros((objects.sum(),3))

#there's surely a simpler way to do this
j=0
for coord in range(coords.shape[0]):
    for i in range(objects[coord]):
            new_coords[j] = coords[coord]
            j+=1

'''new_coords
array([[ 0.66500534,  0.10145476,  0.92528492],
       [ 0.66500534,  0.10145476,  0.92528492],
       [ 0.66500534,  0.10145476,  0.92528492],
       [ 0.66500534,  0.10145476,  0.92528492],
       [ 0.66500534,  0.10145476,  0.92528492],
       [ 0.52677892,  0.07756804,  0.50976737],
       [ 0.52677892,  0.07756804,  0.50976737],
       [ 0.52677892,  0.07756804,  0.50976737],
       [ 0.50030508,  0.37635556,  0.20828815],
       [ 0.50030508,  0.37635556,  0.20828815],
       [ 0.50030508,  0.37635556,  0.20828815],
       [ 0.50030508,  0.37635556,  0.20828815],
       [ 0.50030508,  0.37635556,  0.20828815],
       [ 0.50030508,  0.37635556,  0.20828815],
       [ 0.50030508,  0.37635556,  0.20828815],
       [ 0.50030508,  0.37635556,  0.20828815],
       [ 0.02707651,  0.21878467,  0.55855427],
       [ 0.02707651,  0.21878467,  0.55855427],
       [ 0.02707651,  0.21878467,  0.55855427],
       [ 0.02707651,  0.21878467,  0.55855427],
       [ 0.02707651,  0.21878467,  0.55855427],
       [ 0.81564621,  0.82750694,  0.53083443]])''' 

#calculate distance matrix of old and new arrays
distances_old = distance.pdist(coords)
distances_new = distance.pdist(new_coords)

#calculate and plot normalized histograms (typically just use np.histogram without plotting)
plt.hist(distances_old, range=(0,1), alpha=.5, normed=True)
(array([ 0.,  0.,  0.,  0.,  2.,  1.,  2.,  2.,  2.,  1.]), array([ 0. ,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9,  1. ]), <a list of 10 Patch objects>)

plt.hist(distances_new, range=(0,1), alpha=.5, normed=True)
(array([ 2.20779221,  0.        ,  0.        ,  0.        ,  1.68831169,
        0.64935065,  2.07792208,  2.81385281,  0.34632035,  0.21645022]), array([ 0. ,  0.1,  0.2,  0.3,  0.4,  0.5,  0.6,  0.7,  0.8,  0.9,  1. ]), <a list of 10 Patch objects>)

plt.show()

第二种方法是处理距离矩阵而不是坐标矩阵,但我还没有弄清楚那个代码。

这两种方法对我来说似乎效率低下,我认为操纵 np.histogram 的分箱过程更有可能是有效的,因为它只是基本的乘法,但我不确定如何告诉 np.histogram 将每个坐标视为要计数的对象数量可变。

【问题讨论】:

  • 问题是什么?除了用文字来描述你现有的代码,你能用(复制粘贴运行)代码来描述它吗?
  • 你怎么知道每个位置有多少对象?如果你有一个与位置数组形状相同的计数数组,你可以很容易地用它来加权直方图。
  • 查看我的编辑(你们都快)。是的,对象计数数组与坐标数组的形状相同。
  • @MrE,我认为这个问题很清楚,并且编辑显示了一个代码示例。
  • @askewchan 是的,现在好多了

标签: python numpy scipy histogram


【解决方案1】:

这样的事情可能会奏效:

from scipy.spatial import distance

positions = np.random.rand(10, 2)
counts = np.random.randint(1, 5, len(positions))

distances = distance.pdist(positions)
i, j = np.triu_indices(len(positions), 1)

bins = np.linspace(0, 1, 10)
h, b = np.histogram(distances, bins=bins, weights=counts[i]*counts[j])

与重复检查相比,除了0-distances:

repeated = np.repeat(positions, counts, 0)
rdistances_r = distance.pdist(repeated)

hr, br = np.histogram(rdistances, bins=bins)

In [83]: h
Out[83]: array([11, 22, 27, 43, 67, 46, 40,  0, 19,  0])

In [84]: hr
Out[84]: array([36, 22, 27, 43, 67, 46, 40,  0, 19,  0])

【讨论】:

  • 啊,是的,加权。很简单。谢谢!
  • 是否有更快的方法来实现 np.triu_indices 或 np.histogram?我发现它们都是我代码中的瓶颈。我有包含数千个坐标的“位置”数组,事实证明这对于我的需要(十分之一秒)来说太慢了(~10+ 秒)我原以为 pdist 会成为瓶颈,但事实并非如此。对于 10,000 个坐标数组 pdist 只需要约 0.8 秒(慢,但可接受的边界),但索引需要约 8 秒,直方图需要约 5 秒,这太长了。有什么想法吗?
  • 我使用了triu_indices,因为它最容易理解,但是有更快的方法来生成它。看看我最近的一个问题的一些讨论:Pairwise displacement vectors among set of points。我不认为你可以加快histogram,但我认为它已经尽可能快了,但你可以在别处搜索或提出其他问题。
猜你喜欢
  • 2016-12-06
  • 1970-01-01
  • 2011-08-08
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多