Advent of Code 2019 day 12

    Day 12 turned into a learning journey into Haskell's ad-hoc polymophism and typeclasses.

    Part 1

    I started with the position and velocity vectors represented as V3 Integer values from Data.Linear. A Planet was a record that bundled the position and velocity vectors together, and the four satellites were bundled together into Set of such records.

    import Linear (V3(..), (^+^), (^-^))
    import qualified Data.Set as S
    type Vec = V3 Integer
    data Planet = Planet { _pos :: Vec, _vel :: Vec} deriving (Show, Eq, Ord)
    type Planets = S.Set Planet

    The simulation was pretty straightforward. Each simulationStep applies the changes due to gravity, then applies the changes due to velocity. applyGravity turns out to be O(n2) as the gravity change for each planet depends on the gravitational effect of each other planet. That's why there's both the map and fold in there.

    simulate = iterate simulationStep 
    simulationStep planets = planets''
        where   planets' = applyGravity planets
                planets'' = applyVelocity planets'
    applyGravity planets = (applyGravityHere planets) planets
    applyGravityHere planets here = S.foldl' updateGravity here planets
    updateGravity here there = here { _vel = vel'}
        where   vel = _vel here
                vel' = vel ^+^ gravity ((_pos there) ^-^ (_pos here))
    gravity (V3 x y z) = V3 (signum x) (signum y) (signum z)
    applyVelocity = applyVelocityHere
    applyVelocityHere here = here {_pos = (_pos here) ^+^ (_vel here)}

    The full simulation is just an infinite iterate over the simulationSteps. I find the final answer by droping the first 1000 states (the first of which is the zeroth, starting state) than taking the systemEnergy of the next one.

    part1 :: Planets -> Integer
    part1 = systemEnergy . head . drop 1000 . simulate

    The energy calculation for each planet uses absSum to add up a vector's components, then the systemEnergy is another fold over the set of planets.

    absSum (V3 x y z) = (abs x) + (abs y) + (abs z)
    potentalEnergy planet = absSum $ _pos planet
    kineticEnergy planet = absSum $ _vel planet
    totalEnergy planet = (potentalEnergy planet) * (kineticEnergy planet)
    systemEnergy = (S.foldl' (+) 0) . ( totalEnergy)

    Part 2

    Simplifying part 2 relies on two observations.

    1. The motion of each moon in each dimension is independent of its motion in the other dimensions. That allows us to simulate each dimension independently.
    2. The simulation equations are symmetical in time, which means the first repeated state is the initial state.

    These observations mean that I an find the period of each dimension separately, and the period of the whole system is the least common multiple of these three periods.

    But simulating in one dimension is just about the same as simulating in three dimensions. All the vector calculations are the same, apart from the few low-level functions that explicitly refer to the vector's components. But Haskell's type system means that if a function is defined to operate over V3 Integers, it won't operate over V1 Integers.

    One way to do this is to continue to use V3 Integers, but just two of the components to zero. Another (more interesting) way to do it is to use Haskell's polymorphism features to ensure that all the simulation functions operate over both types of vector.

    I did this by defining my own typeclass, NVec, and defining the operations required to support it. I had to use ^+^^ rather than ^+^ to avoid confusion between my operator and the one from the Data.Linear module. nvZero, nvSignum and nvAbsSum are the low-level operations that are used in the simulation.

    Given the typeclass, I can define how V1 Integer and V3 Integer implement that class, and give the definitions of the required functions.

    class (Eq a) => NVec a where
        (^+^^) :: a -> a -> a
        (^-^^) :: a -> a -> a
        nvZero :: a
        nvSignum :: a -> a
        nvAbsSum :: a -> Integer
    instance NVec (V1 Integer) where
        x ^+^^ y  = x ^+^ y
        x ^-^^ y  = x ^-^ y
        nvZero = V1 0
        nvSignum (V1 x) = V1 (signum x)
        nvAbsSum (V1 x) = abs x
    instance NVec (V3 Integer) where
        x ^+^^ y  = x ^+^ y
        x ^-^^ y  = x ^-^ y
        nvZero = V3 0 0 0
        nvSignum (V3 x y z) = V3 (signum x) (signum y) (signum z)
        nvAbsSum (V3 x y z) = (abs x) + (abs y) + (abs z)

    Given this class, I can now define the data types and functions in terms of this class. For example,

    data Planet a = Planet { _pos :: a , _vel :: a} deriving (Show, Eq, Ord)
    type Planets a = V.Vector (Planet a)
    simulate :: (NVec a) => Planets a -> [Planets a]
    simulate = iterate simulationStep 
    applyGravity :: (NVec a) => Planets a -> Planets a
    applyGravity planets = (applyGravityHere planets) planets

    (Another change I made was in the type of the container for the planets. The state only repeats when each planet is in its original state. But a Set is unordered, so it could be that planet a ends up looking like planet b used to, while b now looks like a. The sets representing these state would be identical, even if the planets themselves aren't. I got around that by using a Vector (fixed-size one-dimenional array) of planets; that allows me to distinguish the planets by their position in the Vector.)

    Once the class was defined, I needed functions to split a single vector of three-dimensional planets into a list of three vectors of one-dimensional planets. These functions need not be polymorphic.

    unzipPlanets :: V.Vector (Planet (V3 Integer)) -> [V.Vector (Planet (V1 Integer))]
    unzipPlanets = dimensionSlice . unzipPlanet
    unzipPlanet :: Planet (V3 Integer) -> [Planet (V1 Integer)]
    unzipPlanet planet = map mkPlanet posVecs
        where posVecs = unzipVec $ _pos planet
    unzipVec :: V3 Integer -> [V1 Integer]
    unzipVec (V3 x y z) = [V1 x, V1 y, V1 z]
    dimensionSlice :: (NVec a) => V.Vector [Planet a] -> [Planets a]
    dimensionSlice slicedPlanets = [sliceDim d | d <- [0..2]]
        where sliceDim d = (!!d) slicedPlanets


    The code is available (and on Github).

    Neil Smith

    Read more posts by this author.