Regular polytopes

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 ...

Platonic solids

What this means is that there are only 5 regular or platonic solids in 3 dimensional space.

That’s it!

These are solids where all the sides are equal and all the faces are the same regular polygon.

n-gons

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.

Approximations

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 and Ad

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 p2 computes 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.

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

Final thoughts

The gravity algorithm has complexity O(n²) and can be improved by using a kd-tree for example

Back to posts