A curious topological result in euclidean spaces is the number of possible regular polytopes in n dimensions. The series (starting with dimension 0) goes like this:
1 1 ∞ 5 6 3 3 3 ...
What this means is that there are only 5 regular or platonic solids in 3 dimensional space.
These are solids where all the sides are equal and all the faces are the same regular polygon.
In 2D there’s an infinite number of perfect polygons (ngons) ⬠ ⬡ . Instead of equivalent faces their only constraints are equal line segments and rotational symmetry (it can be rotated by a certain angle < π an arbitrary number of times and it would look exactly the same). Regular polytopes have all their vertices on an n-hypersphere (they are convex), so for 2D a circle, 3D a sphere.
Meaning if m points are on an n-hypersphere and they are all equidistant from their neighbors they would generate a regular n-polytope with m vertices. In 2D it’s easy,
--assuming polar coordinates in radians vertices :: Floating a => Int -> [a] vertices m = take m $ iterate (+angle) 0 where angle = 2*pi / fromIntegral m
>> vertices 7 [0.0,0.8975979010256552,1.7951958020513104,2.6927937030769655,3.5903916041026207,4.487989505128276,5.385587406153931]
This can generate an m sized regular polygon by dividing a circle into m equal pieces. Since m can be any integer this is proof of the infinite number of regular 2-Polytopes.
Points and spheres
However, in 3D such a simple generator cannot be written as there are only 5 possible regular solids. Which means a sphere cannot be divided perfectly into more than 20 sections (the number of faces of the icosahedron)! any division scheme beyond that number is doomed to be imperfect. This also means you cannot put more than 20 points on a sphere without some being closer to each other’s neighbors than others.
What if we needed to slice spheres into 21 sections? or 100000 sections? Let’s say we wanted to divide earth into approximately equal pieces of land and give one to every human, what can we do?
The answer is approximation. We can start with a random distribution of points on a sphere and perform gradient ascent on the aggregate of the distances of all the points. what this will do is try to push all the points away from each other which would give us a good enough solution.
using stack you can follow along with
stack ghci --install-ghc --resolver=lts-9.14 --package ad --package linear --package split
or if you have nix
nix-shell --pure -I http://nixos.org/channels/nixos-17.09/nixexprs.tar.xz -p 'ghc.withPackages (ps: [ ps.ad ps.linear ps.random ])' --run "cabal repl"
or download and run SpheresAndPoints.hs
And get the imports out of the way
import Numeric.AD import Linear import Data.Function import System.Random import Control.Monad import Control.Monad.State import Control.Applicative import Data.List.Split
Linear is a haskell package used for vector manipulation mostly (rotating vectors and such) it’s very powerful and isn’t really constrained to certain dimensions.
1 / qd p1 p2computes gravity between 2 points and it works for any dimension vector without losing type safety
cost :: ( Metric f , Floating b , Foldable t , Eq (f b) , Functor t) => t (f b) -> b cost points = points & fmap (\p1 -> points & fmap (gravity p1) & sum) & sum where gravity p1 p2 = if p1 == p2 then 0 else (1 / qd p1 p2)
Ad allows the use of automatic differentiation. So the function we want to optimise (cost function) can be specified directly, then we apply
gradientDescent over some initial configuration and we get a list of configurations with decreasing cost.
using the cost function above, the points would be pushed away from each other if the gravity between them is to be minimized
polarsToPoints :: (Eq a,Floating a) => [a] -> [V2 a] polarsToPoints = map angle circleCost :: (Eq a,Floating a) => [a] -> a circleCost = cost.polarsToPoints approximatePerfectCircle :: IO [Double] approximatePerfectCircle = do g <- newStdGen let polars = randomRs (-pi,pi) g & take 10 conjugateGradientDescent circleCost polars & take 10 & map circleCost & return
>> approximatePerfectCircle [323.3763294732505,209.29270123788922,131.12276278726247,111.08208354630008,96.25832338806613,91.4574222020872,86.28924395850842,83.7668331716001,82.86677623362112,82.60473141795445]
the cost is going down meaning gradient descent is working, lets compare it to the perfect one
>> vertex 10 & circleCost 82.5
the approximation reached 82.6 after only 10 steps!
instance (Random a) => Random (V2 a) where random = state random & pure & sequence & runState randomR (a,b) = liftA2 (curry randomR) a b & fmap state & sequence & runState approximateImperfectSphere :: IO [Double] approximateImperfectSphere = do g <- newStdGen let polars = randomRs (pure 0,pure 1) g & take 10 & map sphereAdjustment & map (\(V2 θ ϕ) -> [θ,ϕ]) & concat conjugateGradientDescent sphereCost polars & take 10 & map sphereCost & return where sphereAdjustment (V2 u v) = V2 (2*pi*u) (acos (2*v - 1)) sphericalCoordsToCartesian (V2 θ φ) = radius * V3 (sin θ * cos φ) (sin θ * sin φ) (cos θ) where radius = 5 sphereCost config = config & chunksOf 2 & map (\[θ,ϕ] -> V2 θ ϕ) & map sphericalCoordsToCartesian & cost
the Random instance is there to be able to generate 2 dimensional polar coordinates.
sphereAdjustment exists to make random points on a sphere correctly uniformly distributed
>> approximateImperfectSphere [7.527757046105134,3.4083550282036383,2.5320763637944546,2.252188867357618,2.061022321027768,2.025669773833937,2.0233732580013073,2.0226592876636227,2.0220253327423063,2.0209086575769284]
there now all humans can live in imperfect peace
The gravity algorithm has complexity O(n²) and can be improved by using a kd-tree for example