Gloss.hs

import Graphics.Gloss

import Hypercubes
import PVector

import qualified Data.Vector as DVector

import GHC.Float

import Data.Fixed

dimensions = 3

main =
    animate 
        (InWindow
            "Your Brain on Hypercubes"
            (1920, 1080)        -- window size
            (10, 10))         -- window position
        white
        (animatedCube dimensions)


conViewPlane n = normPlane (DVector.unzip (DVector.generate n (conViewPlaneGen n)))

conViewPlaneGen :: (Integral n) => n -> n -> (Double, Double)
conViewPlaneGen dimensions i = ((cos ((fromIntegral i :: Double) * pi / (fromIntegral dimensions :: Double))),
                                (sin ((fromIntegral i :: Double) * pi / (fromIntegral dimensions :: Double))))

normPlane (u, v) = (normalize u, normalize v)


-- Up vector for n dimensions.
up :: Int -> PVector -> PVector
up n eye = proj eye ((DVector.replicate (n-1) 0.0) `DVector.snoc` 1.0)

-- Gets a right vector for the given up vector.
right :: Int -> PVector -> PVector
right n up
    | (up DVector.! 0) == 0 = (1.0 `DVector.cons` (DVector.replicate (n-1) 0.0))
    | (up DVector.! 1) == 0 = (0.0 `DVector.cons` (1.0 `DVector.cons` (DVector.replicate (n-1) 0.0)))
    | otherwise = ((-(up DVector.! 1)) `DVector.cons` ((up DVector.! 0) `DVector.cons` (DVector.replicate (n-1) 0.0)))


viewPlane :: Int -> PVector -> (PVector, PVector)
viewPlane n eye =
    (up n eye, right n (up n eye))


animatedCube :: Int -> Float -> Picture
animatedCube _n time = do
    let tDouble = float2Double time
    let rotation = (fromIntegral $ truncate $ (time*200) `mod'` 1000) / 1000 :: Double
    let (planeUp, planeRight) = conViewPlane dimensions
    render dimensions (planeUp, planeRight) (rotCube dimensions rotation (centerCube dimensions (Hypercubes.connected dimensions)))


centerCube :: Int -> [(PVector, PVector)] -> [(PVector, PVector)]
centerCube n cube = Prelude.map (\(p1, p2) -> (DVector.map (0.5-) p1, DVector.map (0.5-) p2)) cube

rotCube n angle cube = Prelude.map (\(p1, p2) -> (PVector.rotate p1 (2*pi*angle), PVector.rotate p2 (2*pi*angle))) cube


render n plane lines =
    Pictures [cubeLine n (projectSeg plane seg) | seg <- lines]

-- Projects the given segment onto the given set of view plane vectors.
projectSeg :: (PVector, PVector) -> (PVector, PVector) -> ((Double, Double), (Double, Double))
projectSeg plane (u, v) = (projectPoint plane u, projectPoint plane v)

-- Projects the given point vector onto the given view plane.
projectPoint :: (PVector, PVector) -> PVector -> (Double, Double)
projectPoint (up, right) p = (p ^. right, p ^. up)

cubeLine :: Int -> ((Double, Double), (Double, Double)) -> Picture
cubeLine n ((x1, y1), (x2, y2)) = Line [(double2Float x1 * 100, double2Float y1 * 100), (double2Float x2 * 100, double2Float y2 * 100)]

Hypercubes.hs

-- Faces in a cube or something... My brain hurts.
module Hypercubes where

import Data.Vector
import PVector
import Data.Bits
import qualified Data.Vector as Vector


-- Getting connections.
connected :: Int -> [(PVector, PVector)]
connected n = [(p1, p2) | (idx, p1) <- Vector.toList (indexed (Vector.fromList (cubePoints n))), p2 <- Prelude.drop idx (cubePoints n), isConnected p1 p2]


isConnected :: PVector -> PVector -> Bool
isConnected p1 p2
    | Vector.length p1 == Vector.length p2 =
        connIterFun  p1 p2 0 0

connIterFun :: PVector -> PVector -> Int -> Double -> Bool
connIterFun  p1 p2 idx diffctr
    | idx < Vector.length p1 =
        connIterFun p1 p2 (idx + 1) (diffctr + (if (p1 ! idx) /= (p2 ! idx)
                                                 then 1
                                                 else 0))
    | otherwise = diffctr == 1


-- Returns a list of points (as Int Vectors) in an n-dimensional unit cube.
cubePoints :: Int -> [PVector]
cubePoints n = [toPointVec n p | p <- [0..2^n-1]]

toPointVec :: Int -> Integer -> PVector
toPointVec len p = Vector.generate len (vecFromIntFun len p)

vecFromIntFun :: Int -> Integer -> Int -> Double
vecFromIntFun len pint idx = fromIntegral (if testBit pint (len - idx - 1) then 1 else 0)

PVector.hs

-- Module for dealing with vectors.
module PVector where

import Data.Vector
import qualified Data.Vector as Vector

import Data.Matrix


class Addable v where
    (^+), (^-) :: v -> v -> v


-- A vector representing a point.
type PVector = Vector Double

-- A matrix to be used with PVectors.
type PMatrix = Matrix Double

instance (Num a) => Addable (Vector a) where
    (^+) = Vector.zipWith (Prelude.+)
    (^-) = Vector.zipWith (Prelude.-)


-- Gets the dot product of the two given vectors.
(^.) :: PVector -> PVector -> Double
(^.) u v = Vector.sum (Vector.zipWith (*) u v)

-- Multiplies a vector by a scalar.
(^*) :: PVector -> Double -> PVector
(^*) v n = Vector.map (* n) v


-- Converts the given Integral vector into a PVector
fromVector :: (Integral a) => Vector a -> PVector
fromVector v =
    Vector.map (fromIntegral) v


-- Normalizes the given vector.
normalize :: PVector -> PVector
normalize v =
    -- Divide the vector by its magnitude.
    v ^* (1/(magnitude v))


-- Gets the magnitude of the given vector.
magnitude :: PVector -> Double
magnitude v =
    -- To normalize the vector, get the square root of its dot product with itself.
    sqrt (v ^. v)

-- Gets the projection of u onto v.
proj :: PVector -> PVector -> PVector
proj v u = v ^* ((u ^. v) / (v ^. v))



-- Rotation

-- Rotation of the point vector p by the given angle around the x axis.
rotate :: PVector -> Double -> PVector
rotate p angle =
    getCol 1 ((rotMatrix (Vector.length p) angle) `multStd` (colVector p))


-- An n-dimensional rotation matrix for rotating by the given angle in radians around the given axis.
rotMatrix :: Int -> Double -> PMatrix
rotMatrix dims angle =
    matrix dims dims (rotMatGen dims angle)

rotMatGen :: Int -> Double -> (Int, Int) -> Double
rotMatGen dims angle (y, x)
-- Put the 2D rotation matrix in the top left.
    | (x, y) == (dims-1, dims-1) =  cos angle
    | (x, y) == (dims-1, dims) =  sin angle
    | (x, y) == (dims, dims-1) = -sin angle
    | (x, y) == (dims, dims) =  cos angle
-- Everywhere else, put the identity matrix. That is, if x == y, 1, otherwise 0
    | x == y    = 1
    | otherwise = 0