Pythonで書かれたFriends-of-friendsアルゴリズムがFortran 90/95にある必要があります

Question

私は 'friends-of-friends'アルゴリズムのために自分のコードを書こうとしています。このアルゴリズムは、一連の3Dデータポイントに対して作用し、データセット内の「ハロー」の数を返します。各ハローは、距離がプログラムの唯一のパラメータであるリンクの長さbより小さいポイントの集合です。

アルゴリズムの説明： FOFアルゴリズムには、リンク長と呼ばれる単一の空きパラメーターがあります。リンクの長さ以下の距離で区切られた2つのパーティクルは、「フレンド」と呼ばれます。 FOFグループは、セット内の各パーティクルが、友人のネットワークを介してセット内の他のすべてのパーティクルに接続されているパーティクルのセットによって定義されます。

FOFグループカウンタj = 1に設定します。

各粒子について

• 、nは、まだグループに関連付けられていない。

• 割り当てNグループjに、新しいメンバリストを初期化し、mlist、粒子とグループjに対するN個の第1のエントリとして、

• mlist内の各新しいパーティクルpの再帰的、：

、既にグループjに割り当てられていないものをmlistに追加し、以下の連結長さに等しい距離内にあるPの近傍を探します
• グループjのレコードmlist、j = j + 1を設定します。

これはアルゴリズムをコーディングする試みです。これを行う上で私が快適な唯一の言語はPythonです。しかし、私はこのコードをFortranで書いたり、速くしたりする必要があります。誰かが私を助けてくれることを願っています

まずI 3つのハローの存在を模倣しなければならない点の集合を生成：一つはリストで終わる終わり

x=array(x) 
    y=array(y) 
    z=array(z) 
    id=array(id) 

    t0 = time.time()       

    id_grp=[] 
    groups=zeros((len(x),1)).tolist() 
    particles=id 
    b=1 # linking length 
    while len(particles)>0: 
    index = particles[0] 
    # remove the particle from the particles list 
    particles.remove(index) 
    groups[index]=[index] 
    print "#N ", index 
    dx=x-x[index] 
    dy=y-y[index] 
    dz=z-z[index] 
    dr=sqrt(dx**2.+dy**2.+dz**2.) 
    id_to_look = where(dr<b)[0].tolist() 
    id_to_look.remove(index) 
    nlist = id_to_look 
    # remove all the neighbors from the particles list 
    for i in nlist: 
     if (i in particles): 
      particles.remove(i) 
    print "--> neighbors", nlist 
    groups[index]=groups[index]+nlist 
    new_nlist = nlist 
    while len(new_nlist)>0: 
      index_n = new_nlist[0] 
      new_nlist.remove(index_n) 
      print "----> neigh", index_n 
      dx=x-x[index_n] 
      dy=y-y[index_n] 
      dz=z-z[index_n] 
      dr=sqrt(dx**2.+dy**2.+dz**2.) 
      id_to_look = where(dr<b)[0].tolist() 
      id_to_look = list(set(id_to_look) & set(particles)) 
      nlist = id_to_look 
      if (len(nlist)==0): 
      print "No new neighbors found" 
      else: 
      groups[index]=groups[index]+nlist 
      new_nlist=new_nlist+nlist 
      print "------> neigh-neigh", new_nlist 
      for k in nlist: 
       particles.remove(k) 

：私はFOFアルゴリズムを符号化

import random 
from random import * 
import math 
from math import * 
import numpy 
from numpy import * 
import time 

points = 1000 

halos=[0,100.,150.] 

x=[] 
y=[] 
z=[] 
id=[] 
for i in arange(0,points,1): 
    x.append(halos[0]+random()) 
    y.append(halos[0]+random()) 
    z.append(halos[0]+random()) 
    id.append(i) 

for i in arange(points,points*2,1): 
    x.append(halos[1]+random()) 
    y.append(halos[1]+random()) 
    z.append(halos[1]+random()) 
    id.append(i) 

for i in arange(points*2,points*3,1): 
    x.append(halos[2]+random()) 
    y.append(halos[2]+random()) 
    z.append(halos[2]+random()) 
    id.append(i) 

をリストの中のハローの名前groups

この部分は少し話題ですが、私はあなたにそれを見せてくれるといいと思いました。私は基本的にパーティクルのないグループをすべて削除し、パーティクルの数に従ってソートしていくつかのプロパティを表示しています。

def select(test,list): 
    selected = [] 
    for item in list: 
    if test(item) == True: 
     selected.append(item) 
    return selected 

    groups=select(lambda x: sum(x)>0.,groups) 
    # sorting groups 
    groups.sort(lambda x,y: cmp(len(x),len(y))) 
    groups.reverse() 

    print time.time() - t0, "seconds" 

    mass=x 
    for i in arange(0,len(groups),1): 
    total_mass=sum([mass[j] for j in groups[i]]) 
    x_cm = sum([mass[j]*x[j] for j in groups[i]])/total_mass 
    y_cm = sum([mass[j]*y[j] for j in groups[i]])/total_mass 
    z_cm = sum([mass[j]*z[j] for j in groups[i]])/total_mass 
    dummy_x_cm = [x[j]-x_cm for j in groups[i]] 
    dummy_y_cm = [y[j]-y_cm for j in groups[i]] 
    dummy_z_cm = [z[j]-z_cm for j in groups[i]] 
    dummy_x_cm = array(dummy_x_cm) 
    dummy_y_cm = array(dummy_y_cm) 
    dummy_z_cm = array(dummy_z_cm) 
    dr = max(sqrt(dummy_x_cm**2.+dummy_y_cm**2.+dummy_z_cm**2.)) 
    dummy_x_cm = max(dummy_x_cm) 
    dummy_y_cm = max(dummy_y_cm) 
    dummy_z_cm = max(dummy_z_cm) 
    print i, len(groups[i]), x_cm, y_cm, z_cm,dummy_x_cm,dummy_y_cm,dummy_z_cm 

Answer 1

結果のコードが現在の実装よりも速くなることを願って、Fortranの学習を開始することをお勧めします。最終的にはそうかもしれませんが、他の言語、特に外国語で実装することを考える前にできるだけ速くPythonを実装することをお勧めします。

私はFortranを書いていますが、個人的にはその性能はPython全体にうんざりしますが、これらのことを知っている人は、Python + SciPy + Numpyが慎重に作成した場合、Fortranを計算カーネル多くの科学/工学プログラムお使いのコンピュータのすべてのコアが赤くなってしまうまで、Pythonを最適化していないことを忘れないでください。

So：

第1 - Pythonで動作する実装を取得します。

第2 - できるだけ速く実装してください。

IF（大文字の '大文字'なので大文字）はまだ十分ではなく、コンパイルされた言語に翻訳することのコスト/メリットが有利です。あなたがFortranが広く使われているフィールドにいる場合は、Fortranを覚えておいてください。しかし、それはニッチな言葉のようなもので、C++やその親戚のことを知るためにあなたのキャリアにもっと役立つかもしれません。 （コメントボックスに収めるには長すぎる）

EDIT

なぜ、あなたの質問に私たちを誤解？あなたが快適な唯一の言語はPythonですが、Fortranを知っていると言います。私はあなたがそれに不快でなければならないと思う。そして、あなたのコメントから、あなたが本当に必要とするのは、あなたのPython実装をより速くすることにあるかもしれません。 Sideshow Bobはいくつかのアドバイスを提供しています。それを考慮してから、並列化してください。

Answer 2

より効率的なアルゴリズムへのポインタ。私が間違っていない場合は、ポイントを他のすべてのポイントと比較して、いずれかがリンク長より近いかどうかを確認しています。多数のポイントについて、近くにあるものを見つけるためのより迅速な方法があります - 私の頭の上にある空間インデックスとKDツリーですが、他の方法もあります。

Answer 3

最新のグラフィックスカードをお持ちの場合は、PyOpenCLを使用して、Pythonコードで何百ものプロセッサ（カードに応じて）にパラレル化できます。

あなたはそのアルゴリズムFOFが内部に実装されているかどうかを調べることができ、このvoidfinder F90 code

あなたはSQRT（） の使用を避け、使用することを二乗距離として距離を定義することができ、X * Xの代わりに、X ** 2 ...