Золотой метод спирали
Вы сказали, что не можете заставить работать метод золотой спирали, и это позор, потому что это действительно очень хорошо. Я хотел бы дать вам полное понимание этого, чтобы, возможно, вы могли понять, как не допустить, чтобы это было «сгруппировано».
Итак, вот быстрый, неслучайный способ создания решетки, которая приблизительно верна; как обсуждалось выше, ни одна решетка не будет идеальной, но это может быть «достаточно хорошо». По сравнению с другими методами, например на BendWavy.org , но он просто выглядит красиво и красиво, а также гарантирует равномерный интервал в пределе.
Грунтовка: спирали подсолнуха на диске агрегата
Чтобы понять этот алгоритм, я сначала приглашаю вас взглянуть на алгоритм 2D спирали подсолнечника. Это основано на том факте, что самым иррациональным числом является золотое сечение (1 + sqrt(5))/2, и если кто-то выбрасывает точки при подходе «встаньте в центре, поверните золотое сечение целых поворотов, а затем испустите другую точку в этом направлении» естественным образом строит спираль, которая, когда вы получаете все большее и большее количество точек, тем не менее отказывается иметь четко определенные «столбики», на которых эти точки располагаются в ряд. (Примечание 1.)
Алгоритм равномерного разделения на диске:
from numpy import pi, cos, sin, sqrt, arange
import matplotlib.pyplot as pp
num_pts = 100
indices = arange(0, num_pts, dtype=float) + 0.5
r = sqrt(indices/num_pts)
theta = pi * (1 + 5**0.5) * indices
pp.scatter(r*cos(theta), r*sin(theta))
pp.show()
и выдает результаты, которые выглядят следующим образом (n = 100 и n = 1000):
Радиальное расстояние между точками
Ключевая странная вещь - формула r = sqrt(indices / num_pts); как я дошел до этого? (Примечание 2.)
Ну, я использую здесь квадратный корень, потому что я хочу, чтобы они имели равномерное пространство вокруг сферы. Это то же самое, что сказать, что в пределе больших N я хочу маленький регион R ∈ ( r , r + d r ), Θ ∈ ( θ , θ + d θ ), чтобы содержать число точек, пропорциональных его площадь, которая составляет r d r d θ . Теперь, если мы притворимся, что речь идет о случайной переменной здесь, это будет иметь прямую интерпретацию, заключающуюся в том, что общая плотность вероятности для ( R , Θ ) составляет всего к для некоторой константы c . Затем нормализация на диске устройства приводит к принудительному выполнению c = 1 / & pi;.
Теперь позвольте мне представить трюк. Это происходит из теории вероятностей, где она известна как выборка обратного CDF : предположим, что вы хотели сгенерировать случайную величину с плотностью вероятности f ( z ) и у вас есть случайная переменная U ~ Uniform (0, 1), как и random() в большинстве языков программирования. Как ты это делаешь?
- Сначала превратите вашу плотность в кумулятивную функцию распределения F ( z ), которая, помните, монотонно увеличивается от 0 до 1 с производной F (* * г тысяча девяносто-один * +1092 *). * * 1 093
- Затем вычислите обратную функцию CDF F -1 ( z ).
- Вы обнаружите, что Z = F -1 ( U ) распределяется в соответствии с целевой плотностью. (Примечание 3).
Теперь спиральный трюк с золотым сечением расставляет точки равномерно и равномерно для θ , поэтому давайте интегрируем это; для единичного круга у нас осталось F ( r ) = r 2 . Таким образом, обратная функция равна F -1 ( u ) = u 1/2 , и поэтому мы будет генерировать случайные точки на сфере в полярных координатах с r = sqrt(random()); theta = 2 * pi * random().
Теперь яВместо случайной выборки этой обратной функции мы равномерно выбираем ее, и приятная вещь о равномерной выборке состоит в том, что наши результаты о том, как точки распределены в пределе больших N будет вести себя так, как будто мы случайным образом взяли его. Эта комбинация - хитрость. Вместо random() мы используем (arange(0, num_pts, dtype=float) + 0.5)/num_pts, так что, скажем, если мы хотим взять 10 точек, они равны r = 0.05, 0.15, 0.25, ... 0.95. Мы равномерно выбираем r , чтобы получить интервал равной площади, и мы используем приращение подсолнечника, чтобы избежать ужасных "полос" точек на выходе.
Сейчас занимаемся подсолнухом на сфере
Изменения, которые нам нужно внести, чтобы расставить точки по точкам, просто включают в себя переключение полярных координат на сферические. Радиальная координата, конечно, не входит в это, потому что мы находимся на единичной сфере. Чтобы здесь все было более согласованно, хотя я и был физиком, я буду использовать координаты математиков, где 0 & le; φ & le; &число Пи; широта, падающая с полюса и 0 & le; θ & le; 2 & пи; это долгота. Таким образом, отличие от вышеизложенного состоит в том, что мы в основном заменяем переменную r на & phi; .
Наш элемент площади, который был r d r d θ , теперь становится не намного более сложным грехом ( φ ) д φ д θ . Таким образом, наша плотность соединения для равномерного расстояния является грехом ( φ ) / 4π. Интегрируя θ , мы находим f ( φ ) = sin ( φ ) / 2, таким образом F ( φ ) = (1 & минус; cos ( φ )) / 2. Инвертируя это, мы можем видеть, что равномерная случайная величина будет выглядеть как acos (1 - 2 u ), но мы выбираем равномерно, а не случайно, поэтому вместо этого используем φ k = acos (1 & минус 2 ( k + 0,5) / N ). А остальная часть алгоритма просто проецирует это на координаты x, y и z:
from numpy import pi, cos, sin, arccos, arange
import mpl_toolkits.mplot3d
import matplotlib.pyplot as pp
num_pts = 1000
indices = arange(0, num_pts, dtype=float) + 0.5
phi = arccos(1 - 2*indices/num_pts)
theta = pi * (1 + 5**0.5) * indices
x, y, z = cos(theta) * sin(phi), sin(theta) * sin(phi), cos(phi);
pp.figure().add_subplot(111, projection='3d').scatter(x, y, z);
pp.show()
Опять же, при n = 100 и n = 1000 результаты выглядят так:
Примечания
Эти "столбцы" образованы рациональными приближениями к числу, и наилучшие рациональные приближения к числу получаются из выражения его непрерывной дроби, z + 1/(n_1 + 1/(n_2 + 1/(n_3 + ...))), где z - целое число, а n_1, n_2, n_3, ... - либо конечная или бесконечная последовательность натуральных чисел:
def continued_fraction(r):
while r != 0:
n = floor(r)
yield n
r = 1/(r - n)
Поскольку дробная часть 1/(...) всегда находится между нулем и единицей, большое целое число в непрерывной дробной части дает особенно хорошее рациональное приближение: «деленное на что-то между 100 и 101» лучше, чем «деленное на что-то» между 1 и 2 ". Поэтому самым иррациональным числом является число, равное 1 + 1/(1 + 1/(1 + ...)) и не имеющее особо хороших рациональных приближений; можно решить φ = 1 + 1 / φ , умножив на φ , чтобы получить формулу для золотого сечения.
Для людей, которые не очень знакомы с NumPy - все функции «векторизованы», так что sqrt(array) совпадает с тем, что могут написать другие языки map(sqrt, array). Так что это приложение за компонентом sqrt. То же самое относится и к делению на скаляр или сложению со скалярами - они применяются ко всем компонентам параллельно.
рКрыша проста, если вы знаете, что это результат. Если вы спросите, какова вероятность того, что z <<em> Z <<em> z + d z , это то же самое, что спросить, какова вероятность того, что z <<em> F -1 ( U ) <<em> z + d z , примените F ко всем трем выражениям, отметив, что это монотонно возрастающая функция, поэтому F ( z ) <<em> U <<em> F ( z + d z ), разверните правую сторону, чтобы найти F ( z ) + f ( z ) d z , а поскольку U равномерно, эта вероятность равна f ( z ) d z как и было обещано.