1

我有一组轨迹，由沿轨迹的点组成，并具有与每个点关联的坐标。我将这些存储在 3d 数组中（轨迹、点、参数）。我想找到在这些轨迹的可能成对组合之间具有最大累积距离的一组 r 轨迹。我认为正在工作的第一次尝试如下所示：

max_dist = 0
for h in itertools.combinations ( xrange(num_traj), r):
    for (m,l) in itertools.combinations (h, 2):
        accum = 0.
        for ( i, j ) in itertools.izip ( range(k), range(k) ):
            A = [ (my_mat[m, i, z] - my_mat[l, j, z])**2 \
                    for z in xrange(k) ]
            A = numpy.array( numpy.sqrt (A) ).sum()
            accum += A
    if max_dist < accum:
        selected_trajectories = h

这需要很长时间，因为 num_traj 可以在 500-1000 左右，而 r 可以在 5-20 左右。k 是任意的，但通常可以达到 50。

为了变得超级聪明，我将所有内容都放入了两个嵌套列表推导中，大量使用了 itertools：

chunk = [[ numpy.sqrt((my_mat[m, i, :] - my_mat[l, j, :])**2).sum() \
        for ((m,l),i,j) in \
        itertools.product ( itertools.combinations(h,2), range(k), range(k)) ]\
        for h in itertools.combinations(range(num_traj), r) ]

除了非常难以阅读（!!!）之外，它还需要很长时间。任何人都可以提出任何改进的方法吗？

4

5 回答 5

3

您可以从计算所有轨迹对之间的距离开始，而不是按需重新计算每对轨迹之间的距离。您可以将它们存储在字典中并根据需要查找它们。

这样，您的内部循环for (i,j) ...将被恒定时间查找所取代。

于 2010-05-13T19:53:54.913 回答
2

您可以在距离计算中放弃平方根计算......最大和也将具有最大平方和，尽管这只会产生恒定的加速。

于 2010-05-13T19:59:58.353 回答
2

除了其他人提到的内容之外，这里还有一些兴趣点和建议。（顺便说一句，mathmike 建议生成所有所有对距离的查找列表是您应该立即实施的建议。它从您的算法复杂性中消除了 O(r^2)。）

一、线条

for ( i, j ) in itertools.izip ( range(k), range(k) ):
    A = [ (my_mat[m, i, z] - my_mat[l, j, z])**2 \
        for z in xrange(k) ]

可以替换为

for i in xrange(k):
    A = [ (my_mat[m, i, z] - my_mat[l, i, z])**2 \
        for z in xrange(k) ]

因为 i 和 j 在每个循环中总是相同的。这里根本不需要使用 izip。

二、关于线路

A = numpy.array( numpy.sqrt (A) ).sum()

您确定这是您要计算的方式吗？可能是这样，但这让我感到很奇怪，因为如果这更像是向量之间的欧几里德距离，那么这条线将是：

A = numpy.sqrt (numpy.array( A ).sum())

要不就

A = numpy.sqrt(sum(A))

因为我认为将 A 转换为 numpy 数组以使用 numpy 的 sum 函数比仅使用内置的 Python sum 函数要慢，但我可能错了。此外，如果它确实是您想要的欧几里德距离，那么您将以这种方式做更少的 sqrt。

第三，您是否意识到您可能会尝试迭代多少种潜在组合？对于 num_traj = 1000 和 r = 20 的最坏情况，据我估计，这大约是 6.79E42 个组合。这对于您当前的方法非常棘手。即使对于 num_traj = 500 和 r = 5 的最佳情况，这也是 1.28E12 组合，这是相当多的，但并非不可能。这是您在这里遇到的真正问题，因为通过接受 mathmike 的建议，我提到的前两点并不是很重要。

那你能做什么？好吧，你需要更聪明一点。我还不清楚什么是一个很好的方法。我猜你需要以某种方式使算法启发式。我的一个想法是尝试一种带有启发式的动态编程方法。对于每个轨迹，您可以找到它与另一个轨迹的每一对的距离的总和或平均值，并将其用作适应度度量。在转向三重奏之前，可以放弃一些具有最低健身措施的轨迹。然后，您可以对三重奏做同样的事情：找到每个轨迹所涉及的所有三重奏（在剩余的可能轨迹中）的累积距离的总和或平均值，并将其用作适应度测量来决定在继续之前放弃哪些到四人一组。它没有

于 2010-05-13T21:14:03.863 回答
1

无论如何，它可能需要永远，因为您的算法大约需要 ~ O( C( N, r ) * r^2 )，其中C( N, r )N 选择 r。对于较小的 r（或 N），这可能没问题，但如果您绝对需要找到最大值，而不是使用近似启发式，您应该尝试使用不同的策略进行分支和绑定。这可能适用于较小的 r，并且可以为您节省很多不必要的重新计算。

于 2010-05-13T20:42:04.283 回答
1

这听起来像是一个“加权集团”问题：在一个网络中找到例如 r=5 人，具有最大兼容性/C(5,2) 对权重的最大总和。
谷歌“加权集团”算法 - “集团渗透”→ 3k 点击。
但我会采用 Justin Peel 的方法，因为它是可以理解和可控的
（取 n2 最好的对，从中选出最好的 n3 三元组……调整 n2 n3 ……以轻松权衡运行时间/结果质量。）

5 月 18 日添加，随后是实施的削减。
@Jose，看看什么 nbest[] 序列适合你会很有趣。

#!/usr/bin/env python
""" cliq.py: grow high-weight 2 3 4 5-cliques, taking nbest at each stage
    weight ab = dist[a,b] -- a symmetric numpy array, diag << 0
    weight abc, abcd ... = sum weight all pairs
    C[2] = [ (dist[j,k], (j,k)) ... ]  nbest[2] pairs
    C[3] = [ (cliqwt(j,k,l), (j,k,l)) ... ]  nbest[3] triples
    ...
    run time ~ N * (N + nbest[2] + nbest[3] ...)

keywords: weighted-clique heuristic python
"""
# cf "graph clustering algorithm"

from __future__ import division
import numpy as np

__version__ = "denis 18may 2010"
me = __file__.split('/') [-1]

def cliqdistances( cliq, dist ):
    return sorted( [dist[j,k] for j in cliq  for k in cliq if j < k], reverse=True )

def maxarray2( a, n ):
    """ -> max n [ (a[j,k], (j,k)) ...]  j <= k, a symmetric """
    jkflat = np.argsort( a, axis=None )[:-2*n:-1]
    jks = [np.unravel_index( jk, a.shape ) for jk in jkflat]
    return [(a[j,k], (j,k)) for j,k in jks if j <= k] [:n]

def _str( iter, fmt="%.2g" ):
    return " ".join( fmt % x  for x in iter )

#...............................................................................

def maxweightcliques( dist, nbest, r, verbose=10 ):

    def cliqwt( cliq, p ):
        return sum( dist[c,p] for c in cliq )  # << 0 if p in c

    def growcliqs( cliqs, nbest ):
        """ [(cliqweight, n-cliq) ...] -> nbest [(cliqweight, n+1 cliq) ...] """
            # heapq the nbest ? here just gen all N * |cliqs|, sort
        all = []
        dups = set()
        for w, c in cliqs:
            for p in xrange(N):
                    # fast gen [sorted c+p ...] with small sorted c ?
                cp = c + [p]
                cp.sort()
                tup = tuple(cp)
                if tup in dups:  continue
                dups.add( tup )
                all.append( (w + cliqwt(c, p), cp ))
        all.sort( reverse=True )
        if verbose:
            print "growcliqs: %s" % _str( w for w,c in all[:verbose] ) ,
            print " best: %s" % _str( cliqdistances( all[0][1], dist )[:10])
        return all[:nbest]

    np.fill_diagonal( dist, -1e10 )  # so cliqwt( c, p in c ) << 0
    C = (r+1) * [(0, None)]  # [(cliqweight, cliq-tuple) ...]
        # C[1] = [(0, (p,)) for p in xrange(N)]
    C[2] = [(w, list(pair)) for w, pair in maxarray2( dist, nbest[2] )]
    for j in range( 3, r+1 ):
        C[j] = growcliqs( C[j-1], nbest[j] )
    return C

#...............................................................................
if __name__ == "__main__":
    import sys

    N = 100
    r = 5  # max clique size
    nbest = 10
    verbose = 0
    seed = 1
    exec "\n".join( sys.argv[1:] )  # N= ...
    np.random.seed(seed)
    nbest = [0, 0, N//2] + (r - 2) * [nbest]  # ?

    print "%s  N=%d  r=%d  nbest=%s"  % (me, N, r, nbest)

        # random graphs w cluster parameters ?
    dist = np.random.exponential( 1, (N,N) )
    dist = (dist + dist.T) / 2
    for j in range( 0, N, r ):
        dist[j:j+r, j:j+r] += 2  # see if we get r in a row
    # dist = np.ones( (N,N) )

    cliqs = maxweightcliques( dist, nbest, r, verbose )[-1]  # [ (wt, cliq) ... ]

    print "Clique weight,  clique,  distances within clique"
    print 50 * "-"
    for w,c in cliqs:
        print "%5.3g  %s  %s" % (
            w, _str( c, fmt="%d" ), _str( cliqdistances( c, dist )[:10]))
于 2010-05-15T18:00:12.607 回答