2017-10-31 3 views
3

同様の質問を作成:
Generating N random points with certain predefined distance between themは、与えられた長さL(PythonとN = 200)よりも、すべてのより遠くのNのランダムな点

choose n most distant points in R

をしかし、彼らはMATLABで、または必要なタスクを気力ないかです。

任意の2点間の距離がデルタより大きい長さのボックス内にN個の点を作成する必要があります。

例: x、y、z軸に長さ10オングストロームのボックスがあるとします。
このボックスの内側に200個のランダムな点があるので、任意の2点の間の最小距離 が3オングストロームより大きくなります。

試み:

#!python 
# -*- coding: utf-8 -*-# 
import numpy as np 
np.random.seed(100) 
np.set_printoptions(2) 

box_length = 10 
d = box_length 
threshold = 6 
num_points = 5 
x1, y1, z1 = np.random.random(num_points)* box_length, np.random.random(num_points)* box_length, np.random.random(num_points)* box_length 
x2, y2, z2 = np.random.random(num_points)* box_length, np.random.random(num_points)* box_length, np.random.random(num_points)* box_length 

# print(len(x1)) 

# just for checking make ponts integers 
for i in range(len(x1)): 
    x1[i] = int(x1[i]) 
    x2[i] = int(x2[i]) 
    y1[i] = int(y1[i]) 
    y2[i] = int(y2[i]) 
    z1[i] = int(z1[i]) 
    z2[i] = int(z2[i]) 


print(x1) 
print(y1) 
print(z1) 
print("\n") 

pt1_lst = [] 
pt2_lst = [] 
for i in range(len(x1)): 
    a, b, c = x1[i], y1[i], z1[i] 
    a2, b2, c2 = x2[i], y2[i], z2[i] 
    dist2  = (a-a2)**2 + (b-b2)**2 + (c-c2)**2 

    print("\n") 
    print(a,b,c) 
    print(a2,b2,c2) 
    print(dist2) 

    if dist2 > threshold**2: 
     pt1 = (a,b,c) 
     pt2 = (a2,b2,c2) 
     pt1_lst.append(pt1) 
     pt2_lst.append(pt2) 



print("points") 
print(pt1_lst) 
print(pt2_lst) 

コード上の問題: このコードはpoints2にポイントからポイントを比較したがポイントとpoints2内自体の中で比較されません。

この問題を解決するためのより良いアルゴリズムがあり、問題を解決するために華麗なアイデアを持っている人 に帽子をかけるかもしれません。

ありがとうございました。

PS: 私はいくつかの研究を行なったし、関連リンクを探してみてください、しかし 問題を解決することができませんでした。 まだいくつかの関連リンクは次のとおりです。

更新 ::私はStefans'コードの下にしようとし

、それはN = 10のために動作しますが、私はN = 200のためにそれを試してみましたが、それは(私を非常に大きな時間を利用しています10分後にコードを停止しました)。

効率的な方法はありますか?

本当にありがとうございます!

All paths of length L from node n using python
Create Random Points Inside Defined Rectangle with Python

答えて

1

のは、私はX、Y、Z軸上の長さが10オングストロームの箱を持っているとしましょう。 このボックスの中に10個のランダムな点がありたいので、任意の2点間の最小距離が3オングストロームより大きくなるようにします。

私は距離が全て十分な大きさになるまで繰り返し、そのボックスに10ランダムな点を生成し、これがうまくいくと思う:

>>> import numpy as np 
>>> from itertools import combinations 

>>> while True: 
     P = np.random.rand(10, 3) * 10 
     if all(np.linalg.norm(p - q) > 3 
       for p, q in combinations(P, 2)): 
      break 

>>> P 
array([[ 9.02322366, 6.13576854, 3.1745708 ], 
     [ 6.48005836, 7.5280536 , 4.66442095], 
     [ 5.78306167, 1.83922896, 9.48337683], 
     [ 0.70507032, 0.20737532, 5.31191608], 
     [ 3.71977864, 6.40278939, 3.81742814], 
     [ 0.03938102, 6.7705456 , 6.28841217], 
     [ 3.27845597, 2.98811665, 4.81792286], 
     [ 7.74422021, 9.30027671, 8.1770998 ], 
     [ 0.28544716, 0.35155801, 9.77847352], 
     [ 4.84536373, 4.21378476, 0.4456017 ]]) 

は、ポイントの良いセットを見つけるために、約50の試みをとります。ここで私は1000回、20回それを試してみました:それは良かった:

>>> sum(all(np.linalg.norm(p - q) > 3 
      for p, q in combinations(np.random.rand(10, 3) * 10, 2)) 
     for _ in range(1000)) 
20 
+0

KDtreeを使用して解決できませんでした。より多くの点については、 'scipy.spatial.distance.pdist'と' scipy.spatial.distance.squareform'を使って計算を高速化することができます。 – Joe

+0

@ジョーええ、より大きい数字の方が良いかもしれません。彼らはすべての距離または最小距離を計算しますが、そうですか?最初も同じことをやっていましたが、小さすぎる距離のペアを見つけると直ちに停止して再試行するほうがはるかに速いことに気付きました。これは、この自作ブルートフォースチェックの利点です。 –

+0

@StefanPochmann、これはN = 10でとても良いですが、N = 200でこのアルゴリズムを実行しようとすると時間がかかりすぎるかもしれません!!! – astro123

1

このような分布の一般的な名前はポアソン球のサンプリングです。これを行うO(n)があります - 確認してくださいhere

+0

偉大な情報をありがとう、あなたが気づいている場合は、この種のデータを生成するための作業コードがある場合、それはよかったです! – astro123

1

これはFORTRANを使用して実行できます。あなたの質問にネームを200に変更してください。

! file: N_random_numbers.f90 
PROGRAM makemodel 
IMPLICIT NONE 
INTEGER :: kwrite(3),i,j,naccepted,failed_attempts 
INTEGER, PARAMETER:: natom1=64,natom2=96,natom3=0,natoms=160 
DOUBLE PRECISION:: lbox,density,molarmass,atoms_per_meter3,box_volume,side_length 
!DOUBLE PRECISION, PARAMETER:: avogadro=6.02214129e23 
DOUBLE PRECISION, DIMENSION(3,3):: lat_vector 
DOUBLE PRECISION, DIMENSION(natoms,3):: aposition 
DOUBLE PRECISION, DIMENSION(1,3):: random 
DOUBLE PRECISION:: local1,local2,local3,x,y,z 
DOUBLE PRECISION:: distance,distance_fraction,distance_thereshold,inner_side 

DATA kwrite/301,302,303/ 

OPEN(kwrite(1),FILE='out_07_Al2O3_ramdom_POSCAR',STATUS='UNKNOWN') 


!########## Files, Initializations ############! 
!##############################################! 


distance_thereshold=1.4       ! (Choose the separation between the atoms thats acceptable) 
lbox= 11.950           ! (Choose what you want it to be) 

! Note:: Box size is given by (latt_con * ncell) 
!#################################BODY OF CALCULATION##############################! 

!******** Finding the box size *********! 
!atoms_per_meter3=5.424958e28       ! From Journal of Non-Crystalline Solids 262(2000) 
box_volume=(lbox)**3 
!box_volume=natoms*1.0e30/atoms_per_meter3   ! In Cubic Angstron 
side_length=box_volume**(1.d0/3.d0)     ! In Angstron 
inner_side=side_length-distance_thereshold 

!write(*,*)side_length 

!####### Writing POSCAR file ########! 
!##########################################! 

write(kwrite(1),*) 'Comment line' 
write(kwrite(1),303)side_length 
write(kwrite(1),303)'1.0000','0.0000','0.0000' 
write(kwrite(1),303)'0.0000','1.0000','0.0000' 
write(kwrite(1),303)'0.0000','0.0000','1.0000' 
write(kwrite(1),303)'AlO' 
write(kwrite(1),303)natoms 
write(kwrite(1),*)'Direct' 

!******** Find and write the atomic positions *********! 
aposition=0.d0 

CALL init_random_seed() 
CALL RANDOM_NUMBER(random) 

aposition(1,:)=random(1,:)*inner_side 
naccepted=1 

105 CONTINUE 
failed_attempts=0 
205 CONTINUE 
CALL RANDOM_NUMBER(random) 
DO i=1,naccepted 
local1=(aposition(i,1)-(random(1,1)*inner_side))**2.d0 
local2=(aposition(i,2)-(random(1,2)*inner_side))**2.d0 
local3=(aposition(i,3)-(random(1,3)*inner_side))**2.d0 
distance=sqrt(local1+local2+local3) 

IF(distance.LT.distance_thereshold)THEN 
failed_attempts=failed_attempts+1 
GO TO 205 
END IF 

END DO 

naccepted=naccepted+1 
aposition(naccepted,:)=random(1,:)*inner_side 
WRITE(*,'(A,I8,A,I10,1x,A)')'found',naccepted,'th position after',failed_attempts,'failed attempts' 

IF (naccepted.LT.natoms)THEN 
GO TO 105 
ELSE 
GO TO 305 
END IF 


305 DO i=1,naccepted 
x=aposition(i,1)+(distance_thereshold/2.0) 
y=aposition(i,2)+(distance_thereshold/2.0) 
z=aposition(i,3)+(distance_thereshold/2.0) 
write(kwrite(1),302) x/lbox,y/lbox,z/lbox 

END DO 

!***************************************************************! 

CLOSE(kwrite(1)) 

301 format(1x,f10.6) 
302 format(3(f16.12,1x)) 
303 format(1x,A,1x,A,1x,A) 

END PROGRAM 

!#################### Subroutine for random numbers generation ###########################! 
!*****************************************************************************************! 

SUBROUTINE init_random_seed() 
INTEGER :: i, n, clock 
INTEGER, DIMENSION(:), ALLOCATABLE :: seed 

CALL RANDOM_SEED(size = n) 
ALLOCATE(seed(n)) 

CALL SYSTEM_CLOCK(COUNT=clock) 

seed = clock + 37 * (/ (i - 1, i = 1, n) /) 
CALL RANDOM_SEED(PUT = seed) 

DEALLOCATE(seed) 
END SUBROUTINE 
関連する問題