2010-05-13 6 views
1

나는 궤적을 따르는 점들과 각 점과 연관된 좌표들로 구성된 궤도 집합을 가지고 있습니다. 나는 이것을 3 차원 배열 (탄도, 점, 매개 변수)에 저장합니다. 나는이 궤도의 가능한 pairwise 조합 사이에 최대 누적 거리를 갖는 r 궤적 세트를 찾고 싶다. num_traj는 500-1000 주위 수거리 메트릭의 조합 최적화

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 

이 영원히 소요, r은 약 5 ~ 20 일 수있다 : 내가 생각하는 내 첫 번째 시도는이 같은 외모 노력하고 있습니다. 꽤되는 것을

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

외에도 : k는 itertools를 많이 사용하고, 임의이지만, 슈퍼 영리 시도 50

에, 나는 두 개의 중첩 지능형리스트에 모든 것을 넣어 가지고까지 일반적으로 할 수있다 읽을 수없는 (!!!), 그것은 또한 오랜 시간이 걸립니다. 누구든지이 문제를 개선 할 수있는 방법을 제안 할 수 있습니까?

답변

3

필요에 따라 각 궤도 쌍 사이의 거리를 다시 계산하기보다는 모든 궤도 쌍 사이의 거리를 계산하여 시작할 수 있습니다. 사전에 저장할 수 있으며 필요에 따라 검색 할 수 있습니다.

이렇게하면 내부 루프 for (i,j) ...이 상수 시간 조회로 바뀝니다.

+0

이 속도가 빨라졌습니다! 감사! – Jose

2

거리 계산에서 제곱근 계산을 버릴 수 있습니다 ... 최대한의 합계는 최대 제곱 합계를 갖지만 일정한 속도 향상 만 발생합니다.

1

알고리즘이 ~ O(C(N, r) * r^2) (약 C(N, r)) 인 경우 r을 선택하면 영원히 오래 걸릴 것입니다. 더 작은 r (또는 N)에 대해서는 괜찮을 지 모르지만 근사치를 사용하는 것과는 대조적으로 최대 값을 절대적으로 찾아야하는 경우 가지와 다양한 전략으로 묶어보십시오. 이것은 더 작은 r에 대해 작동 할 수 있으며 불필요한 재 계산에 꽤 많은 시간을 절약합니다.

2

여기에 다른 사람들이 언급 한 것 외에도 몇 가지 관심과 제안이 있습니다. (덧붙여서, mathmike가 룩업리스트를 생성하는 제안은 모든 쌍거리를 즉각적으로 배치해야한다는 것입니다. 알고리즘 복잡도에서 O (r^2)를 제거합니다.)

처음 i 및 j는 항상 각 루프에서 동일하기 때문에, 라인

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

로 대체 될 수있다. 여기에서 izip을 사용할 필요가 없습니다. 라인

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

당신은 당신이 그것을 계산하는 방법이있다 확신에 대한

둘째,?

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

하거나

A = numpy.sqrt(sum(A)) 

을 I는 변환 있다고 생각하기 때문이이 벡터 사이의 유클리드 거리의 더 있다면, 행이 일 것이기 때문에 아마 그렇게하지만, 그냥 이상한 나를 강타 numpy의 sum 함수를 사용하려면 numpy 배열을 사용하는 것이 내장 파이썬 합계 함수를 사용하는 것보다 느려지지만 잘못 될 수 있습니다. 또한, 만약 당신이 진정으로 유클리드 거리가 원한다면, 그러면 sqrt를 이렇게 적게 할 것입니다.

셋째, 반복 할 가능성이있는 조합의 수를 알 수 있습니까? num_traj = 1000 및 r = 20 인 최악의 경우, 대략 6.79E42 조합입니다. 그것은 현재의 방법으로는 다루기가 어렵습니다. num_traj = 500 및 r = 5 인 경우에도 1.28E12의 조합이 가능하지만 불가능하지는 않습니다. 이것은 mathmike의 조언을 듣고 내가 언급 한 처음 두 점이 그리 중요하지 않기 때문에 여기에있는 실제 문제입니다.

그럼 어떻게 할 수 있습니까? 음, 좀 더 영리해야합니다. 나에게 아직 명확하지 않은 것은 무엇인지 이것에 대한 훌륭한 방법이 될 것이다. 나는 당신이 알고리즘을 어떤 식 으로든 만들어 낼 필요가 있다고 생각한다. 한 가지 생각은 경험적 방법으로 동적 프로그래밍 방식을 시도하는 것입니다. 각 궤적마다 다른 궤적과의 모든 짝짓기 거리의 합 또는 평균을 구하고이를 체력 측정 도구로 사용할 수 있습니다. 운동량이 가장 낮은 궤적 중 일부는 트리오로 이동하기 전에 떨어질 수 있습니다. 그런 다음 트리오로 같은 일을 할 수 있습니다 : 각 궤적이 포함 된 모든 트리오 (남아있는 가능한 궤적 중에서)에 대한 누적 거리의 합 또는 평균을 찾아 이동하기 전에 어느 것을 떨어 뜨릴 지 결정하기위한 적합성 척도로 사용하십시오 foursomes. 그것은 최적의 솔루션을 보장하지는 않지만 꽤 좋은 것이어야하고, 내가 믿는 솔루션의 시간 복잡성을 크게 낮출 것입니다.

1

"가중치 크리크"문제로 들립니다. 예 : 예 : r = 최대 호환성을 가진 네트워크에서 5 명/최대 C (5,2) 쌍의 합계.
구글 "weighted clique"알고리즘 - "clique percolation"→ 3k 히트.
하지만 나는 그것을 이해하고
제어이기 때문에 (저스틴 껍질의 방법 과 함께 할 것입니다 그들로부터 최고의 N3 트리플을 N2 최고의 쌍을 것입니다 ... 은 ... 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]))