Я пытаюсь написать свой собственный код для алгоритма «друзей друзей». Этот алгоритм воздействует на набор трехмерных точек данных и возвращает количество «ореолов» в наборе данных. Каждый гало представляет собой набор точек, расстояние которых меньше длины связывания, 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