1

軌跡に沿った点で構成され、各点に関連付けられた座標を持つ一連の軌跡があります。これらを 3D 配列 (trajectory、point、param) に保存します。これらの軌道の可能なペアごとの組み合わせ間の累積距離が最大になる 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 を多用して、すべてを 2 つのネストされたリスト内包表記に入れました。

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

ここでは、他のすべての人が言及したことに加えて、いくつかの興味深い点と提案があります。(ちなみに、すべてのペア距離すべてのルックアップリストを生成するという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関数を使用するよりも遅いと思うからですが、間違っている可能性があります。また、それが本当に必要なユークリッド距離である場合は、この方法でより少ない平方根を実行することになります。

第三に、反復しようとしている可能性のある組み合わせがいくつあるかわかりますか?num_traj=1000およびr=20の最悪の場合、これは私の見積もりでは約6.79E42の組み合わせです。それはあなたの現在の方法ではかなり手に負えません。num_traj=500およびr=5の最良の場合でさえ、それは1.28E12の組み合わせであり、かなりの数ですが、不可能ではありません。これはあなたがここで抱えている本当の問題です。なぜなら、マトミケのアドバイスを取り入れることによって、私が言及した最初の2つのポイントはそれほど重要ではないからです。

それでは何ができますか?さて、あなたはもう少し賢くする必要があります。これに最適な方法が何であるかはまだ私にはわかりません。何らかの方法でアルゴリズムをヒューリスティックにする必要があると思います。私が持っていた1つの考えは、ヒューリスティックを使用して動的計画法のようなアプローチを試すことでした。軌道ごとに、別の軌道とのペアリングごとの距離の合計または平均を求め、これを適合度の尺度として使用できます。適応度の測定値が最も低い軌道の一部は、トリオに移動する前にドロップされる可能性があります。次に、トリオで同じことを行うことができます。各軌道が関与するすべてのトリオ(残りの可能な軌道の中から)の累積距離の合計または平均を見つけ、それをフィットネス測定として使用して、次に進む前にドロップするものを決定します。フォーサムに。それはしません

于 2010-05-13T21:14:03.863 に答える
2

距離計算で平方根計算を捨てることができます...最大和にも最大二乗和がありますが、それは一定のスピードアップしか得られません。

于 2010-05-13T19:59:58.353 に答える
1

これは「重み付けされたクリーク」の問題のように聞こえます。たとえば、最大の互換性 / C(5,2) ペアの重みの最大合計を持つネットワーク内の r=5 人を見つけます。
Google の「加重クリーク」アルゴリズム -「クリーク パーコレーション」 → 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 に答える
1

とにかく永遠にかかる可能性がO( C( N, r ) * r^2 )ありC( N, r )ます.Nはrを選択します. 小さい r (または N) の場合はこれで問題ないかもしれませんが、近似ヒューリスティックを使用するのではなく、絶対に最大値を見つける必要がある場合は、さまざまな戦略で分岐と境界を試す必要があります。これは r が小さい場合に機能する可能性があり、不要な再計算を大幅に節約できます。

于 2010-05-13T20:42:04.283 に答える