Алгоритм друзей друзей, написанный на Python, должен быть в Fortran 90/95

Я пытаюсь написать свой собственный код для алгоритма "друзей друзей". Этот алгоритм воздействует на набор трехмерных точек данных и возвращает количество "ореолов" в наборе данных. Каждый гало представляет собой набор точек, расстояние которых меньше длины связывания, b, единственный параметр программы.

Алгоритмическое описание: алгоритм FOF имеет единственный свободный параметр, называемый длиной связывания. Любые две частицы, которые разделены расстоянием, меньшим или равным длине связывания, называются "друзьями". Затем группа FOF определяется набором частиц, для которых каждая частица в наборе связана с каждой другой частицей в наборе через сеть друзей.

Установите счетчик групп FOF j=1.

  • Для каждой частицы n, еще не связанной с какой-либо группой:

  • Присвойте n группе j, инициализируйте новый список членов, mlist, для группы j с частицей n в качестве первой записи,

  • Рекурсивно для каждой новой частицы p в mlist:

  • Найдите соседей p, лежащих на расстоянии, меньшем или равном длине связи, добавьте к mlist тех, кто еще не назначен группе j,
  • Запишите mlist для группы j, установите j=j+1.

Это моя попытка закодировать алгоритм. Единственный язык, на котором мне удобно это Python. Однако мне нужно, чтобы этот код был написан на Фортране или чтобы он был быстрее. Я действительно надеюсь, что кто-нибудь поможет мне.

Сначала я создаю набор точек, которые должны имитировать наличие 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)

Затем я кодировал алгоритм 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)

В конце получается список ореолов в списке 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
Источник

3 ответа

Решение

Я думаю, что вам не следует начинать изучение Фортрана в надежде, что полученный код будет быстрее, чем ваша текущая реализация. Возможно, в конечном итоге, но я думаю, что вам лучше посоветовать сделать реализацию Python как можно быстрее, прежде чем думать о реализации на другом языке, особенно на иностранном.

Я пишу на Fortran, и лично я думаю, что его производительность снижается во всем Python, но люди, которые знают об этих вещах, приводят убедительные аргументы в пользу того, что Python+SciPy+Numpy, если их тщательно продумать, могут сопоставить Fortran с точки зрения скорости в вычислительных ядрах многих научных / инженерных систем программы. Не забывайте, что вы не оптимизировали свой Python до тех пор, пока все ядра на вашем компьютере не станут красными.

Так:

1-й - получить работающую реализацию на Python.

2-й - сделать вашу реализацию максимально быстрой.

ЕСЛИ (заглавные буквы, потому что это большое "если"), код все еще недостаточно быстр, и стоимость / выгода его перевода на скомпилированный язык выгодна, ПОЭТОМУ подумайте, на какой скомпилированный язык его перевести. Если вы находитесь в области, где Фортран широко используется, то изучайте Фортран во что бы то ни стало, но это что-то вроде нишевого языка, и это может принести пользу вашей карьере, чтобы выучить С ++ или одного из его родственников.

РЕДАКТИРОВАТЬ (слишком долго, чтобы поместиться в поле для комментариев)

Зачем вводить нас в заблуждение в вашем вопросе? Вы утверждаете, что единственный язык, с которым вам удобно, это Python, теперь вы говорите, что знаете Fortran. Я полагаю, тебе это неудобно. И из вашего комментария кажется, что, возможно, вам действительно нужна помощь в ускорении реализации Python; Sideshow Bob предложил несколько советов. Примите это во внимание, затем распараллелите.

Источник

Указатель на более эффективный алгоритм. Если я не ошибаюсь, вы сравниваете точку с любой другой точкой, чтобы увидеть, ближе ли она к длине связи. Для большого количества точек есть более быстрые способы поиска ближайших соседей - пространственная индексация и деревья KD на макушке, но, несомненно, есть и другие методы, которые будут работать для вас.

Источник

Если у вас есть современная видеокарта, вы можете парализовать сотни процессоров (в зависимости от вашей карты) в коде Python, используя PyOpenCL.

Вы можете проверить, реализован ли алгоритм FoF в этом коде voidfinder F90.

Вы можете определить расстояние как квадратное расстояние, чтобы избежать использования sqrt() и использовать x*x вместо x**2...

Источник
Другие вопросы по тегам