2010-05-13 54 views
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) ] 

除了是相當不可讀(!!!),這也需要很長時間。任何人都可以提出任何改進方法嗎?

回答

3

與其重新計算每對軌跡之間的距離按需,您可以從計算所有軌跡對之間的距離開始。您可以將它們存儲在字典中,並根據需要進行查找。

這樣你的內循環for (i,j) ...將被替換爲恆定時間查找。

+0

這一直是加速的主要來源!謝謝! – Jose 2010-05-14 12:22:11

2

您可以在距離計算中對平方根計算...最大和也將具有最大平方和,但只能產生恆定的加速。

1

無論如何,這可能需要永遠,因爲你的算法需要約O(C(N, r) * r^2),其中C(N, r)是N選擇r。對於較小的r(或N),這可能沒有問題,但是如果您絕對需要找到最大值,而不是使用近似啓發式,則應該嘗試使用不同策略的分支和約束。這可能適用於較小的r,它可以節省您不必要的重新計算。

2

除了其他人提到的,還有幾點興趣點和建議。 (順便說一句,mathmike提出的生成查找表的建議是所有所有對的距離都應該立即設置,從算法複雜度中排除O(r^2)。)

First中,線

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)) 

,因爲我認爲轉換到一個numpy數組使用numpy的sum函數比使用內置的Python sum函數要慢,但我可能是錯的。另外,如果它真的是你想要的歐幾里德距離,那麼你將以這種方式少做sqrt。

第三,你是否意識到有多少潛在的組合可能試圖迭代?對於num_traj = 1000和r = 20的最差情況,按我的估計約爲6.79E42個組合。這與目前的方法相當棘手。即使對於num_traj = 500和r = 5的最佳情況,這是1.28E12的組合,這是很多,但不是不可能的。這是你在這裏遇到的真正問題,因爲通過採用mathmike的建議,我提到的前兩點並不是非常重要。

那你能做什麼?那麼,你需要更聰明一些。目前還不清楚,這將是一個很好的方法。我猜你需要以某種方式使算法成爲啓發式。我想過的一個想法是嘗試一種啓發式的動態編程方法。對於每個軌跡,您可以找到每個軌跡與其他軌跡配對的距離的總和或平均值,並將其用作適應性度量。一些健身措施最低的軌跡可能會在投入三重奏之前被丟棄。然後,您可以對三重奏做同樣的事情:找出每個軌跡所涉及的所有三重奏(剩餘可能的軌跡中)的累計距離的總和或平均值,並將其用作適應性測量,以決定在繼續前進行哪些下降四人組合。它不能保證最佳的解決方案,但它應該是相當不錯的,它會大大降低我相信解決方案的時間複雜度。

1

這聽起來像一個「加權集團」的問題:找到例如 r =網絡中5個人的最大兼容性/ C(5,2)對權重的最大總和。
Google「加權集團」算法 - 「集團滲透」→ 3k次點擊。
而是因爲它是可以理解的,可控的
我將與賈斯汀剝離的方法 去(走N2最好的對,從他們的最好N3三倍...... 調整N2 N3 ......輕鬆地權衡的結果的運行時間/質量。 )

增加了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]))