Metoda złotej spirali
Powiedziałeś, że nie możesz uruchomić metody złotej spirali i szkoda, bo jest naprawdę, naprawdę dobra. Chciałbym dać ci pełne zrozumienie tego, żebyś mógł zrozumieć, jak nie dopuścić do tego, by się „zgniatało”.
Oto szybki, nielosowy sposób tworzenia sieci, która jest w przybliżeniu poprawna; jak omówiono powyżej, żadna krata nie będzie idealna, ale może to wystarczyć. Porównuje się go do innych metod, np. Na BendWavy.org, ale ma po prostu ładny i ładny wygląd oraz gwarancję równych odstępów w limicie.
Podkład: spirale słonecznika na dysku jednostkowym
Aby zrozumieć ten algorytm, najpierw zapraszam do przyjrzenia się algorytmowi spirali słonecznika 2D. Jest to oparte na fakcie, że najbardziej irracjonalną liczbą jest złoty podział (1 + sqrt(5))/2i jeśli ktoś emituje punkty, stosując podejście „stań w środku, obróć złoty stosunek całych obrotów, a następnie emituj kolejny punkt w tym kierunku”, w naturalny sposób konstruuje się spirala, która w miarę jak dochodzi do coraz większej liczby punktów, niemniej jednak odmawia posiadania dobrze zdefiniowanych „słupków”, na których punkty się pokrywają. (Notatka 1.)
Algorytm równych odstępów na dysku to:
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()
i daje wyniki, które wyglądają następująco (n = 100 in = 1000):
Rozmieszczenie punktów promieniowo
Kluczową dziwną rzeczą jest formuła r = sqrt(indices / num_pts); jak doszedłem do tego? (Uwaga 2.)
Cóż, używam tutaj pierwiastka kwadratowego, ponieważ chcę, aby miały równe odstępy wokół dysku. To jest to samo, co powiedzieć, że w granicy dużego N chcę, aby mały obszar R ∈ ( r , r + d r ), Θ ∈ ( θ , θ + d θ ) zawierał liczbę punktów proporcjonalną do jego powierzchni, czyli r d r d θ . Teraz, jeśli udajemy, że mówimy tutaj o zmiennej losowej, ma to prostą interpretację, mówiąc, że łączna gęstość prawdopodobieństwa dla ( R , Θ ) jest po prostu crdla jakiejś stałej c . Normalizacja na dysku jednostkowym wymusiłaby wówczas c = 1 / π.
Pozwólcie, że przedstawię sztuczkę. Pochodzi z teorii prawdopodobieństwa, gdzie jest znane jako próbkowanie odwrotnego CDF : załóżmy, że chcesz wygenerować zmienną losową o gęstości prawdopodobieństwa f ( z ) i masz zmienną losową U ~ Jednolity (0, 1), tak jak wychodzi z random()w większości języków programowania. Jak Ty to robisz?
- Najpierw zamień swoją gęstość w skumulowaną funkcję dystrybucji lub CDF, którą nazwiemy F ( z ). Pamiętaj, że CDF rośnie monotonicznie od 0 do 1 z pochodną f ( z ).
- Następnie oblicz funkcję odwrotną funkcji CDF F -1 ( z ).
- Przekonasz się, że Z = F -1 ( U ) rozkłada się zgodnie z gęstością docelową. (Uwaga 3).
Teraz spirala złotego podziału rozdziela punkty w ładnie równy wzór dla θ, więc zintegrujmy to; dla dysku jednostkowego pozostaje F ( r ) = r 2 . Czyli funkcja odwrotna to F -1 ( u ) = u 1/2 , a zatem wygenerowalibyśmy losowe punkty na dysku we współrzędnych biegunowych z r = sqrt(random()); theta = 2 * pi * random().
Teraz zamiast losowo próbkować tę funkcję odwrotną, próbkujemy ją równomiernie , a fajną rzeczą w próbkowaniu jednorodnym jest to, że nasze wyniki dotyczące rozłożenia punktów na granicy dużego N będą się zachowywać tak, jakbyśmy próbkowali ją losowo. Ta kombinacja jest sztuczką. Zamiast tego random()używamy (arange(0, num_pts, dtype=float) + 0.5)/num_pts, więc powiedzmy, że jeśli chcemy pobrać próbkę 10 punktów, to one są r = 0.05, 0.15, 0.25, ... 0.95. Jednolicie próbujemy r, aby uzyskać równe odstępy, i używamy przyrostu słonecznika, aby uniknąć okropnych „słupków” punktów w wyniku.
Teraz robię słonecznik na kuli
Zmiany, które musimy wprowadzić, aby kropkować kulę punktami, polegają jedynie na zmianie współrzędnych biegunowych na współrzędne sferyczne. Współrzędna promieniowa oczywiście nie wchodzi w to, ponieważ jesteśmy na kuli jednostkowej. Aby zachować trochę spójności w tym miejscu, mimo że byłem wyszkolony jako fizyk, użyję współrzędnych matematyków, gdzie 0 ≤ φ ≤ π to szerokość geograficzna schodząca z bieguna, a 0 ≤ θ ≤ 2π to długość geograficzna. Zatem różnica z powyższego polega na tym, że w zasadzie zastępujemy zmienną r przez φ .
Nasz element obszaru, który był r d r d θ , teraz staje się niewiele bardziej skomplikowanym sin ( φ ) d φ d θ . Więc nasza gęstość spoiny dla równomiernych odstępów wynosi sin ( φ ) / 4π. Całkując θ , otrzymujemy f ( φ ) = sin ( φ ) / 2, więc F ( φ ) = (1 - cos ( φ )) / 2. Odwracając to, możemy zobaczyć, że jednolita zmienna losowa wyglądałaby jak acos (1 - 2 u ), ale próbkujemy jednakowo zamiast losowo, więc zamiast tego używamy φ k = acos (1 - 2 ( k+ 0,5) / N ). Reszta algorytmu po prostu rzutuje to na współrzędne x, y i 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()
Ponownie dla n = 100 in = 1000 wyniki wyglądają następująco:
Dalsze badania
Chciałem się pochwalić blogiem Martina Robertsa. Zauważ, że powyżej utworzyłem przesunięcie moich indeksów, dodając 0,5 do każdego indeksu. To było dla mnie po prostu atrakcyjne wizualnie, ale okazuje się, że wybór offsetu ma duże znaczenie i nie jest stały w przedziale i może oznaczać nawet o 8% lepszą dokładność pakowania, jeśli zostanie wybrany poprawnie. Powinien również istnieć sposób, aby jego sekwencja R 2 pokryła kulę i byłoby interesujące zobaczyć, czy to również stworzyło ładne, równe pokrycie, być może takie, jakie jest, ale być może powinno być, powiedzmy, wzięte tylko z połowy Kwadrat jednostkowy przeciął lub tak po przekątnej i rozciągnął się, aby uzyskać okrąg.
Uwagi
Te „słupki” są tworzone przez racjonalne przybliżenia liczby, a najlepsze racjonalne przybliżenia liczby pochodzą z jej ciągłego wyrażenia ułamkowego, z + 1/(n_1 + 1/(n_2 + 1/(n_3 + ...)))gdzie zjest liczbą całkowitą i n_1, n_2, n_3, ...jest skończoną lub nieskończoną sekwencją dodatnich liczb całkowitych:
def continued_fraction(r):
while r != 0:
n = floor(r)
yield n
r = 1/(r - n)
Ponieważ część ułamkowa 1/(...)zawsze zawiera się w przedziale od zera do jedynki, duża liczba całkowita w ułamku ciągłym pozwala na szczególnie dobre przybliżenie racjonalne: „jedna podzielona przez coś między 100 a 101” jest lepsza niż „podzielona przez coś między 1 a 2”. Dlatego najbardziej irracjonalna liczba jest liczbą, która jest 1 + 1/(1 + 1/(1 + ...))i nie ma szczególnie dobrych racjonalnych przybliżeń; φ = 1 + 1 / φ można rozwiązać mnożąc przez φ, aby otrzymać wzór na złoty współczynnik.
Dla osób, które nie są tak zaznajomione z NumPy - wszystkie funkcje są „wektoryzowane”, więc sqrt(array)jest to to samo, co można napisać w innych językach map(sqrt, array). Więc to jest aplikacja komponent po komponencie sqrt. To samo dotyczy dzielenia przez skalar lub dodawania ze skalarami - te dotyczą wszystkich składników równolegle.
Dowód jest prosty, gdy wiesz, że to jest wynik. Jeśli zapytasz, jakie jest prawdopodobieństwo, że z < Z < z + d z , to jest to to samo, co pytanie, jakie jest prawdopodobieństwo, że z < F -1 ( U ) < z + d z , zastosuj F do wszystkich trzech wyrażeń, zauważając, że jest funkcja monotonicznie rosnąca, stąd F ( z ) < U < F ( z + d z ), rozwiń prawą stronę na zewnątrz, aby znaleźć F ( z ) + f( z ) d z , a ponieważ U jest jednorodne, to prawdopodobieństwo wynosi tylko f ( z ) d z, jak obiecano.
(gdzie rzeczy =