Я пытаюсь написать собственный код для алгоритма «друзей друзей». Этот алгоритм действует на множество трехмерных точек данных и возвращает количество «ореолов» в наборе данных. Каждое гало представляет собой набор точек, расстояние которых меньше длины связывания, b, единственный параметр программы.Алгоритм друзей друзей, написанный на Python, должен быть в Fortran 90/95
Алгоритмическое описание: Алгоритм FOF имеет единственный свободный параметр, называемый длиной связывания. Любые две частицы, которые разделены расстоянием, меньшим или равным длине связывания, называются «друзьями». Затем группа FOF определяется набором частиц, для которых каждая частица внутри множества связана с каждой другой частицей в множестве через сеть друзей.
Комплект группового счетчика FOF j = 1.
Для каждой частицы, п, еще не связана с любой группой:
Присвоить п к группе J, инициализировать новый список членов, mlist, для группы J с частицей п в качестве первой записи,
Рекурсивный, для каждой новой частицы р в mlist:
- Найти соседей р, лежащих на расстоянии, меньшем или равном длине связующей, добавить к mlist, уже не отнесены к группе J,
- Запись для группы j, набор 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)
Тогда я кодированный алгоритм 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
Каковы параметры для решения? Нужно ли минимальное количество ореолов? –
У вас есть проблема с алгоритмом или с его внедрением в код? Сначала получите алгоритм правильно, прежде чем пытаться его кодировать. – eriktous
Я добавил описание алгоритма, вы можете мне помочь? – Brian