用 Python 编写的朋友之友算法需要在 Fortran 90/95

Bri*_*ian 1 python fortran finder halo

我正在尝试为“朋友的朋友”算法编写自己的代码。该算法作用于一组 3d 数据点,并返回数据集中“晕圈”的数量。每个光环是距离小于链接长度 b 的点的集合,b 是程序的唯一参数。

算法描述: FOF 算法有一个称为链接长度的自由参数。任何两个间隔小于或等于链接长度的粒子都称为“朋友”。然后,FOF 组由一组粒子定义,其中该组中的每个粒子都通过朋友网络连接到该组中的每个其他粒子。

设置 FOF 组计数器 j=1。

  • 对于尚未与任何组关联的每个粒子 n:

  • 将 n 分配给组 j,为组 j 初始化一个新成员列表 mlist,其中粒子 n 作为第一个条目,

  • 递归地,对于 mlist 中的每个新粒子 p:

  • 在小于或等于链接长度的距离内找到 p 的邻居,将那些尚未分配到组 j 的邻居添加到 mlist,
  • 记录组 j 的 mlist,设置 j=j+1。

这是我对算法进行编码的尝试。我擅长的唯一语言是 Python。但是,我需要用 Fortran 编写此代码或使其更快。我真的希望有人能帮助我。

首先,我生成一组点来模拟 3 个光晕的存在:

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)
Run Code Online (Sandbox Code Playgroud)

然后我编码了FOF算法:

  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)
Run Code Online (Sandbox Code Playgroud)

最后得到一个列表中的光环列表 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
Run Code Online (Sandbox Code Playgroud)

Hig*_*ark 5

我认为您开始学习 Fortran 是不明智的,因为希望生成的代码比您当前的实现更快。最终可能是这样,但我认为最好建议您在考虑用另一种语言(尤其是外语)实现之前尽可能快地实现 Python。

我写 Fortran,我个人认为它的性能在整个 Python 中都很糟糕,但是了解这些事情的人提供了令人信服的论点,即如果精心设计,Python+SciPy+Numpy 可以在许多科学/工程的计算内核中与 Fortran 相媲美程式。不要忘记,在您计算机上的所有内核都运行得非常火热之前,您还没有优化您的 Python。

所以:

1st - 在 Python 中获得一个有效的实现。

第二 - 使您的实施尽可能快。

IF(大写字母,因为它是一个很大的“if”)代码仍然不够快,将其翻译成编译语言的成本/收益是有利的 THEN 考虑将其翻译成哪种编译语言。如果你在 Fortran 被广泛使用的领域,那么一定要学习 Fortran,但它是一种小众语言,学习 C++ 或它的一个亲戚可能会让你的职业生涯受益更多。

编辑(太长,无法放入评论框)

为什么在你的问题中误导我们?你说你唯一熟悉的语言是 Python,现在你说你知道 Fortran。我想你一定对此感到不舒服。而且,从您的评论看来,您真正需要的帮助似乎是使您的 Python 实现更快;Sideshow Bob提供了一些建议。考虑到这一点,然后并行化。