当生成具有约束的随机数时，while 循环不一致地崩溃

Question

从向量开始，vector0
初始化生成另一个随机向量的while循环，vector1
用点积计算它们之间的夹角
如果vector0和vector1之间的夹角theta太大，继续重新制作vector1直到足够小

看起来像这样：

# initialize the angle
theta = 0
# the first vector:
vector0 = [x0, y0, z0]
# initialize while loop:
while theta <= 0 or theta > np.pi/8:
    # create the second vector using random numbers
    x1 = random.uniform(-maxlen, maxlen)
    y1 = random.uniform(-maxlen, maxlen)
    z1 = random.uniform(-maxlen, maxlen)
    vector1 = [x1, y1, z1]
    # find the angle between the two vectors. The loop will start again if it is too large.
    theta = np.arccos(np.dot(vector0, vector1) / np.linalg.norm(vector0)*np.linalg.norm(vector1)

此过程嵌套在另外两个循环中 - 不是特别大的循环，只有 5 步和 100 步。我认为这是一个足够简单的过程。

这是我的问题：这个 while 循环大约有 70% 的时间崩溃。只是放弃。但有些时候，它工作得很好！

杀死它并重新初始化很容易，但有时我要重复十次才能使代码成功通过运行，这变得难以忍受。

我是不是做了什么蠢事导致了这个？也许有一个错误有时会在我的代码中触发，或者我犯了一个数学错误？也许有更多 memory/CPU-efficient 的方法来实现这个结果？还是我只需要使用更强大的机器？

Answer 1

这是一种生成随机向量而不需要检查它是否在要求的角度内的方法：

import numpy as np
import math

max_phi = np.pi/8

v1 = np.array([1, 1, 1])
phi = np.random.rand()*max_phi
psi = np.random.rand()*2*np.pi

# rotate v1 in the plane created by v1 and [0, 0, 1]
# unless v1 is parallel to [0, 0, 1], then use the plane normal to [1, 0, 0]
if (v1/np.sum(v1**2)**0.5).T @ np.array([0, 0, 1]) == 1:
    axis = np.array([1, 0, 0])
else:
    axis = np.cross(v1, np.array([0, 0, 1]))

def rotation_matrix(axis, theta):
    """
    Return the rotation matrix associated with counterclockwise rotation about
    the given axis by theta radians.
    """
    axis = np.asarray(axis)
    axis = axis / math.sqrt(np.dot(axis, axis))
    a = math.cos(theta / 2.0)
    b, c, d = -axis * math.sin(theta / 2.0)
    aa, bb, cc, dd = a * a, b * b, c * c, d * d
    bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
    return np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
                     [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
                     [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])

# find the rotation matrix about the 'axis' axis
R0 = rotation_matrix(axis, phi)
# find the rotation matrix about the v1
R1 = rotation_matrix(v1  , psi)
# apply random rotations to create a random vector withing an angle of phi 
# radians from v1
v2 = R1@R0@v1

注意随机向量的分布会有所不同。越接近原始向量的向量，生成的概率越高。

当生成具有约束的随机数时，while 循环不一致地崩溃

While loop crashes inconsistently when generating random numbers with constraint

python

numpy

random

while-loop

crash