A Different Kind of Cube

Among the things that came into my possession last Christmas (thanks, Tom and Kate and Annie and James!) was a collection of wooden puzzles that includes this thing:

The cubes have an elastic string running through them, so that adjacent cubes can be rotated with respect to each other around the axis connecting their centres, but are otherwise fixed.

The aim of the puzzle is to fold this string of cubes into a larger cube. This turns out to be difficult, so I decided to do the traditional thing (c.f. this work from 2012 and everything that arose from it) and write a Haskell program to solve it for me.

Without further ado:

   

import Data.Bits (setBit, testBit)

type Snake = String
snake :: Snake
snake = "bwb/w/b/wb/w/bw/b/w/bw/bw/b/w/b/wb/wb/wb/wb"

This represents the snake that we’re trying to fit into a cube. ‘b’ and ‘w’ represent black and white blocks, and ‘/’ represents a change of direction. The only thing that matters for the solution is the location of the direction changes in this list; the ‘b’ and ‘w’ just make it easier to check that the representation matches the real object.

type Cube = Int
emptyCube = 0

This represents the state of the cube into which we’re trying to fit the snake. For a backtracking search, all we really need to record here is which of the 3x3x3 cells in the cube have already been filled - so that we don’t produce solutions where the snake intersects itself. So we can use the first 27 bits of an Int.

Now we’re ready to define functions for working with the cube: one to mark a cell as filled and one to see whether it’s filled.

bitIndex :: Position -> Int
bitIndex (x, y, z) = x*9 + y*3 + z
fill :: Cube -> Position -> Cube
fill cube pos = setBit cube (bitIndex pos)
isFilledAt :: Cube -> Position -> Bool
isFilledAt cube pos = testBit cube (bitIndex pos)

Our algorithm will work its way through the snake and the cube simultaneously, marking cells in the cube as used as it goes. It will need to know where it is in the cube, and in what direction it’s working.

type Position = (Int, Int, Int)
data Axis = X | Y | Z deriving (Eq, Show)
data Sign = Positive | Negative deriving (Eq, Show)
data Direction = Direction Axis Sign

If the algorithm encounters a change-direction marker in the snake then it can go in any of the four directions that are perpendicular to the current direction.

nextDirections :: Direction -> [Direction]
nextDirections (Direction a _) = [Direction a' d |
                                  a' <- [X, Y, Z],
                                  a' /= a,
                                  d <- [Positive, Negative]]

Otherwise, it will just move to the next cell, in the direction in which it’s already heading.

nextPosition :: Position -> Direction -> Position
nextPosition (x, y, z) (Direction a d) =
  case a
  of X -> (x+delta, y, z)
     Y -> (x, y+delta, z)
     Z -> (x, y, z+delta)
  where delta = if d == Positive then 1 else -1

We’ll need to know when we’ve reached the edge of the space.

outOfBounds :: Position -> Bool
outOfBounds (x, y, z) = or [i < 0 || i > 2 | i <- [x, y, z]]

Now we have all the parts we need to put together a search algorithm. We’ll build up a list of the directions in which the segments of the snake need to lie; that’s enough information to follow and it’s easier to interpret than, say, a list of the cell coordinates.

type PartialSolution = [Direction]
type Solution = PartialSolution

Here’s the core of the algorithm: a function that consumes one character of the snake description, works out what possible next steps there are, and collects the solutions found by each of those steps. At each step, it needs to know:

• which cells of the cube have already been filled;
• what the last-filled position and direction are;
• what part of the snake remains to be placed in the cube;
• what solution has been accumulated so far.

solutions' :: Cube
          -> Position
          -> Direction
          -> Snake
          -> PartialSolution
          -> [Solution]

If there’s nothing of the snake left to place, then we’ve found one solution. This is the base case for the recursive search.

solutions' _ _ dir [] soFar = [dir:soFar]

If the next character in the snake indicates a change of direction, then we have to look for solutions in each of the four possible next directions. We don’t change position or mark any cell as filled, but we’ll need to include the change of direction in any solutions we find.

solutions' cube prevPos dir ('/':snake') soFar =
  concat [solutions' cube prevPos dir' snake' (dir:soFar) |
          dir' <- nextDirections dir]

Otherwise, we must mark our position as filled - or return empty-handed if that’s impossible - and move to the next cell, in the same direction as before.

solutions' cube prevPos dir (c:snake') soFar =
  let pos = nextPosition prevPos dir in
  if outOfBounds pos || cube `isFilledAt` pos then []
  else solutions' (fill cube pos) pos dir snake' soFar

A small wrapper describes the starting conditions for the search.

solutions snake = map reverse $ solutions' emptyCube
                                           (-1, 0, 0)
                                           (Direction X Positive)
                                           snake
                                           []

Finally, we will show the directions in a solution as crude arrows.

instance Show Direction where
  show (Direction X Positive) = ">"
  show (Direction X Negative) = "<"
  show (Direction Y Positive) = "^"
  show (Direction Y Negative) = "v"
  show (Direction Z Positive) = "."
  show (Direction Z Negative) = "o"

main = mapM_ putStrLn $ map (unwords . map show) $ solutions snake

We get just two solutions:

> ^ < . ^ o > . v < ^ o ^ . v > ^
  > . < ^ . v > ^ o < . v . ^ o > .

which are actually just the unique solution and its mirror image.

Counting to Ten

Earlier I worked out how to fit ten square faces into a single Soma cube. I've now also finished the logic for working out how to assemble and dismantle a particular solution of the cube. It turned out not to be too complicated. And the code for incorporating that into the resulting animation wasn't too tricky either. So, as promised, here's the video:

That is all.

A Ten-Faced Cube

It's been a busy couple of weeks. Among other things, I investigated how many distinct solution faces can be encoded in a single Soma cube. The answer turns out to be ten – the theoretical maximum – which is a surprising result and perhaps even a new one. A little help from Google's math mailing list helped me find that very quickly.

I also:

• stripped out the custom Recipe representation in favour of a more general one built on Data.Vect;
• built most of a representation of 128-bit words;
• finished code for correct assignment of UV coordinates to shape faces, and added those UV coordinates to Haskell export and Python import;
• split the cube face materials into their own Blender library and added scripts to automate their generation.

Excluding comments, the source code current stands at 747 lines of Haskell and 235 lines of Python.

Next, I want to beautify the animation a bit, and then I'll produce another video. In the mean time, here's a big chunk of new code, which computes the tenfold encoding and the corresponding UV coordinates.

> module Stats where

> import Control.Monad
> import Data.Bits
> import Data.Function
> import Data.List
> import Data.Map ((!))
> import qualified Data.Map as Map
> import Data.Maybe
> import Data.Monoid
> import Data.Ord
> import Data.Vect hiding (transpose)
> import Data.Word

> import Cube hiding (main)
> import Mesh hiding (testMeshes)
> import Misc
> import Word128

This counts how many squares cover the surface of a shape in the Soma cube.

A less brute-force way to arrive at the same numbers would be to see that:

• each shape can be built up by starting with a single cube, and sticking successive cubes on by a single face;
• the initial cube has six faces;
• each subsequent cube adds four more faces (it brings six faces into play, but one of its own and one of the existing ones are stuck together and hidden). So, the six four-cube shapes must all have 6 + 3 × 4 = 18 faces, and the one three-cube shape must have 6 + 2 × 4 = 14 faces.

> countFaces = length . allFaces . defaultForm

That makes 122 faces in total. (In any solution, only 6 × 9 = 54 faces can go on the outside, so the other 68 must be inside the cube, in 34 pairs of adjacent faces.)

> totalFaces = sum [countFaces shape | shape <- enumFrom L]

The following defines an offset into a bitfield at which to find the bits for each face. Since there are 122 of them, they pack neatly into 128 bits.

> shapeBitOffset = Map.fromList $ foldl reserveBits [] (enumFrom L)
>  where
>   reserveBits [] shape = [(shape, 0)]
>   reserveBits os@((shape', o):_) shape = (shape, o + countFaces shape'):os

Given a predicate that selects faces of shapes, the corresponding subset of faces can be represented as 128 bits.

> markingsAsBits :: [(Shape, [Bool])] -> Word128
> markingsAsBits markings = foldl setBit (fromInteger 0) bitsToSet
>  where
>   bitsForShape (shape, marking) = [shapeBitOffset!shape + i |
>                                    (i, marked) <- zip [0..] marking,
>                                    marked]
>   bitsToSet = concatMap bitsForShape markings

> mapFaces :: (a -> b) -> [(Shape, [a])] -> [(Shape, [b])]
> mapFaces f = map (\(shape, faces) -> (shape, map f faces))

> markFaces :: (a -> Bool) -> [(Shape, [a])] -> Word128
> markFaces f = markingsAsBits . mapFaces f

'paintings' takes a solution and lists, for each shape, how the shape's faces would be coloured by that solution.

> paintings :: Solution -> [(Shape, [FacePaint])]
> paintings solution = [(shape, map snd faces) |
>                       (shape, (_, faces)) <- allMeshes solution]

Now there's enough information to work out how many of the shapes' faces are ever on the outside of the cube:

> totalEverExternalFaces = popCount everExternalBits
>  where
>   everExternalBits = foldr1 (.|.) [markFaces (isJust . cubeFace) $
>                                    paintings s | s <- solutions' RotationOnly]

This reveals that 90 of the 122 faces are ever on the outside. Since each face of a solution requires 9 faces of a shape, that means there is a theoretical limit of 90/9 = 10 disjoint solution faces that can simultaneously be assigned to the faces of the pieces.

Obviously, six disjoint solution faces can be simultaneously assigned to the faces of the pieces, just by choosing all the solution faces from the same solution. But that's not very interesting. The interesting problem is to have as many solution faces represented as possible, but where each one comes from a different solution.

In other words, the aim is to find as large a set of solutions as possible, such that there is some choice of one face from each solution having the property that the chosen faces all occupy disjoint sets of piece faces.

This function takes an assignment of cube faces to shape faces, and splits it into a separate list describing, for each cube face, which shape faces it comprises.

> explodePainting :: [(Shape, [FacePaint])] -> [(CubeFace, [(Shape, [Bool])])]
> explodePainting p = [(cf, filterPaint cf p) |
>                      cf <- [CubeFace sign axis |
>                             sign <- [Negative, Positive],
>                             axis <- enumFrom X]]
>  where
>   filterPaint cf shapePaints = [(shape, [match cf fp | fp <- facePaints]) |
>                                 (shape, facePaints) <- shapePaints]
>   match cf (FacePaint cf') = (Just cf == cf')

Applying that splitting function to the face assignments of each solution in turn results in the following function. It produces a list of 6 × 480 = 2880 elements.

> faceInfos :: [((Solution, Int), CubeFace, Word128, [(Shape, [Bool])])]
> faceInfos = [((solution, n), cf, markingsAsBits fs, fs) |
>               (solution, n) <- zip (solutions' RotationOnly) [0..],
>               (cf, fs) <- explodePainting (paintings solution)]
> solutionIndex ((_, n), _, _, _) = n
> bitRep        (_,      _, w, _) = w

Of course, there's some chance that a cube face from one solution will cover exactly the same shape faces as a cube face from another solution (albeit probably with different orientations). Ideally, a set of cube faces would not contain any cube face that is also formed by a different solution.

> distinctFaceInfos = nubBy ((==) `on` bitRep) faceInfos

This reveals that there are only 1176 distinct cube faces. And even fewer will appear on only one solution. That's unsurprising in fact, because while symmetries of individual shapes are discounted, it's often the case that a solution will contain a subunit comprising two or more shapes, where the subunit has a symmetry. In such a solution, the subunit can be removed, transformed by its symmetry, and recombined, to give a new solution; any cube faces not affected by this will then be shared between the two solutions.

> groupedFaces = groupSortBy bitRep faceInfos
> faceFreq = histogram groupSizes
>  where
>   groupSizes = map length groupedFaces

The above function produces this histogram: [(1,619),(2,258),(3,104),(4,72),(5,39),(6,19),(7,22),(8,6),(9,6),(10,7), (11,2),(12,1),(13,4),(14,2),(15,1),(16,2),(18,2),(20,2),(25,1),(28,2), (30,2),(38,1),(46,1),(48,1)]

That is to say, there are 619 cube faces that can only be made by a single solution; there are 258 cube faces that can each only be made by two solutions, and so on ... up to one cube face that appears in 48 different solutions and another that appears in 46.

Here are the solution faces that appear only in a single solution:

> interestingFaceInfos = concat [g | g <- groupedFaces, length g == 1]

These can then be grouped by which solution they appear in:

> interestingSolutionFaces = sortBy (flip $ comparing length) $
>                            groupSortBy solutionIndex interestingFaceInfos

> interestingSolutionFreq = histogram $ map length interestingSolutionFaces

The above function produces this histogram: [(1,226),(2,131),(3,35),(4,4),(5,2)]

So, 226 solutions have one unique face each. But two solutions have five unique faces each!

Although the cube faces between them cover 90 different shape faces, there is some redundant information being encoded, because a shape has several coplanar faces, and when a cube face is assigned to a shape face, it must be assigned to all coplanar faces on the same shape. For example, when a cube face contains one of the L-shaped sides of shape L, it must contain all four of the shape faces comprised by that L-shaped side.

Eliminating this redundant information could be done earlier (at the face painting stage) by comparing faces' geometry, but here will do just as well and keeps things simple, by avoiding the need for a representation of a mapping from cell face to shape face.

By removing this redundancy, the cube-face-to-shape-face mapping can be reduced to just 42 bits per cube face, which is good because it will fit in a native 64-bit word now.

> bitGroups = sortBy (comparing head) $
>             dupBits $ map bitRep interestingFaceInfos

> compressRep :: Word128 -> Word64
> compressRep w = foldl (mapBit w) 0 bitMap
>  where
>   mapBit :: (Bits a, Bits b) => a -> b -> (Int, Int) -> b
>   mapBit w w' (outBit, inBit) = if (testBit w inBit)
>                                 then setBit w' outBit
>                                 else w'
>   bitMap = zip [0..] (map head bitGroups)

> uncompressRep :: Word64 -> Word128
> uncompressRep w' = foldl (unmapBit w') 0 bitUnmap
>  where
>   unmapBit :: (Bits a, Bits b) => a -> b -> (Int, [Int]) -> b
>   unmapBit w' w (inBit, outBits) = if (testBit w' inBit)
>                                    then foldl setBit w outBits
>                                    else w
>   bitUnmap = zip [0..] bitGroups

The following working form now contains, for each interesting solution face, a deduplicated representation of the shape faces that it uses, and also an Integer describing which of the other 618 interesting solution faces cannot appear in conjunction with it.

> compressedFaceInfos = [(n, s, cf, compressRep b, p) |
>                        (n, ((s, _), cf, b, p)) <- zip [0..]
>                                                       interestingFaceInfos]

Each cube face needs 9 cell faces.

> targetComboSize = totalEverExternalFaces `div` 9

The above working provides the basis for a naive depth-first search for solutions.

> search partialSolution n
>        usedFaceBits
>        bestN
>        infos =
>     if n > bestN
>     then let (restSolutions, restBestN) = rest n infos
>          in (partialSolution : restSolutions, max n restBestN)
>     else rest bestN infos
>  where
>   rest bestN' [] = ([], bestN')
>   rest bestN' (info@(_,_,_,faceBits,_):infos) =
>     if (usedFaceBits .&. faceBits == 0)
>     then let (withThisInfo, bestN'W) = search (info:partialSolution) (n+1)
>                                               (usedFaceBits .|. faceBits)
>                                               bestN'
>                                               infos
>              (withoutThisInfo, bestN'') = search partialSolution n
>                                                  usedFaceBits
>                                                  bestN'W
>                                                  infos
>          in (withThisInfo ++ withoutThisInfo, bestN'')
>     else search partialSolution n
>                 usedFaceBits
>                 bestN'
>                 infos

> run' minLength = search [] 0 0 (minLength-1)
> run = run' targetComboSize compressedFaceInfos

> chosenCombo' = head $ fst run

Of course, this simple approach runs straight into the combinatrial explosion of having to choose ten elements from a set of 619; considering each possibility at a time is going to take longer than a lifetime. Luckily, there are other ways to attack the problem and once there's a solution it can be plugged in quickly.

The following was found very quickly (that's to say, in less than a second) with a linear constraint solver. Thanks to Google's math mailing list for that suggestion.

> chosenCombo = map (compressedFaceInfos!!)
>                   [68, 139, 152, 208, 234, 336, 363, 386, 487, 513]

The naive search may not be efficient enough for finding a solution, but it's suitable for checkin a solution. This confirms that the selection made by the linear constraint solver does have the properties required.

> chosenComboIsValid = (not . null . fst) $ run' targetComboSize chosenCombo

Having found a set of ten solution faces that completely cover the shape faces, it's time to condense those faces into one data set that describes, for each shape face, which solution face it belongs in and where in the solution face it appears. That will be enough to asssign a texture to the shape face, and also texture coordinates.

> type ShapeFaceData = Maybe (Int, Recipe)

This elaborate list comprehension fits together all the exploded face information from earlier, into a single structure that collects the facts for each face of each shape.

> implodePainting :: [(Int, Solution, [(Shape, [Bool])])]
>                 -> [(Shape, [ShapeFaceData])]
> implodePainting explodedMarkings = [
>     (shape, superpose paintGroups) |
>     (shape, paintGroups) <- [head $ groupKey grp |
>                              grp <- groupSortBy fst .
>                              concat $ explodedShapeFaceDatas]]
>  where
>   superpose :: [[ShapeFaceData]] -> [ShapeFaceData]
>   superpose = map msum . transpose
>   explodedShapeFaceDatas :: [[(Shape, [ShapeFaceData])]]
>   explodedShapeFaceDatas = [attachMaterials (i, solution, shapeFacePaints) |
>                            (i, solution, shapeFacePaints) <- explodedMarkings]

The painting condensation logic relies on this helper function, whose purpose is to multiply solution face information out onto each face that appears on the solution face.

> attachMaterials :: (Int, Solution, [(Shape, [Bool])]) ->
>                                    [(Shape, [ShapeFaceData])]
> attachMaterials (cubeFaceNumber, solution, shapeFacePaints) = [
>  (shape, [if paint then Just (cubeFaceNumber, recipe) else Nothing |
>           paint <- facePaints]) |
>  (shape, recipe, facePaints) <- merge solution shapeFacePaints]
>   where merge (Solution s) (facePaints) = [
>          (checkEqual "Yikes!" shape shape', recipe, paints) |
>          ((shape, recipe, _), (shape', paints) )<- zip s facePaints]

> checkEqual msg s1 s2 = if s1 == s2 then s1 else error msg

The solution information can now be split into two results.

The first result is an association of solution to cube face, which allows for successive face solution animations to point the same way, even though different chosen solution faces have different orientations.

The second result describes, for each face of each shape, which solution it belongs to (if any) and the position of the shape in that solution (if any), which is necessary for assigning UV coordinates.

> implodeCombo combo = (solutionToCubeFace, implodePainting nsps)
>  where (solutionToCubeFace, nsps) = unzip [((s, cf), (n, s, p)) |
>                                            (n, (_, s, cf, _, p)) <-
>                                            zip [0..] combo]

A couple of small helper functions ease export to a Python-readable form.

> pyVec2 (Vec2 x y) = (x, y)
> pyMaterial Nothing       = "Face" ++ "X"
> pyMaterial (Just (n, _)) = "Face" ++ show n

> painter faceMaterials shape = fromJust . flip lookup faceInfo
>  where 
>   faceInfo = [(face, (pyMaterial mat, uvs mat (faceVerts sign axis face))) |
>               ((sign, axis, face), mat) <-
>               zip (allFaces $ defaultForm shape)
>                   faceMaterials]

Always-internal shape faces don't add up to nice shapes, so they can each display a full texture.

>   uvs Nothing _ = [(0,0),(0,1),(1,1),(1,0)]

For sometimes-external shape faces, the UV coordinates are chosen by transforming the face's vertices into their position in the solution to which the shape face belongs, then seeing whether the vertices appear on the corresponding cube face. At this point there's no information about which cube face it's meant to be, so that is reconstructed from the vertices' normal. (It ought to be possible to keep the cube face information up to this point, then use 'toXaxis' to avoid the normal calculation. But that's for another time.)

>   uvs (Just (_, recipe)) verts = [pyVec2 $ extractUvs vert' &* (1/3) |
>                                   vert' <- verts']
>     where
>       verts' = map (vApplyRecipe recipe) verts  -- into solution space
>       extractUvs = (apply . extractor) verts'
>   apply :: (Vec3 -> Float, Vec3 -> Float) -> Vec3 -> Vec2
>   apply (exU, exV) vec = Vec2 (exU vec) (exV vec)

It's perhaps possible to simplify this mapping, but how is not obvious.

>   extractor [a,b,c,_] =  ex $ pyVert uvNormal
>    where uvNormal = (b &- a) &^ (c &- b)
>          ex (1, 0, 0) = (p _2, p _3)
>          ex (0, 1, 0) = (n _1, p _3)
>          ex (0, 0, 1) = (p _2, n _1)
>          -- For negative normals, just flip the V coordinate
>          ex (a, b, c) = let (eu, ev) = ex (-a, -b, -c) in (p eu, n ev)
>          n = ((3-) .)
>          p = id

> testMeshes = do
>   let (solutionToCubeFace, shapeFaces) = implodeCombo chosenCombo
>   writeMeshes [(shape, mesh (painter faceInfos shape) shape) |
>                (shape, faceInfos) <- shapeFaces]
>   writeSolutions [Solution [(shape, recipe .*. toXaxis sign axis, bitmap) |
>                             (shape, recipe, bitmap) <- s] |
>                   (Solution s, CubeFace sign axis) <- solutionToCubeFace]
>   importIntoBlender

> main = testMeshes

The Moment You Didn't Know You Were Waiting For

Following on from the export side of the Soma Cube work, here's the import side, in Python, to be run as a Blender script. This is more of an exercise in finding one's way around Blender's object model than anything else, so there's not much point in a deep commentary on it. Each shape of the cube is loaded into a separate Mesh object, and each solution of the cube becomes a set of keyframes for those objects.

import math
import bpy
import io

file = io.open

# RGB components are in the range [0,1].
#
# NMesh.materials is a list of Material objects. NMesh.faces[i].mat is the
# index into the mesh's material list of the material to use for the i'th
# face.

# When assigning to a face's vertex index list, a final zero is ignored.
# So, when creating a square face, rotate its vertex list if necessary
# to make sure vertex zero can be specified.
def RawVertices(verts):
  if len(verts) == 4 and verts[-1] == 0:
    return [0] + verts[:-1]
  else:
    return verts

def ReadShapes():
  meshes = eval(file('shapes.txt').read())

  materials = {}
  for (_, (_, faces)) in meshes:
    for i, (_, mat_name) in enumerate(faces):
      if mat_name not in materials:
        n = len(materials)
        material = bpy.data.materials.new(mat_name)
        material.diffuse_color = [(n & (1<<j)) and 1.0 or 0.3 for j in range(3)]
        materials[mat_name] = material

  for i, (name, (verts, faces)) in enumerate(meshes):
    poly = bpy.data.meshes.new(name)

    poly.vertices.add(len(verts))
    for j, v in enumerate(verts):
      poly.vertices[j].co = v

    mesh_materials = []
    for (_, mat_name) in faces:
      if mat_name not in mesh_materials:
        mesh_materials.append(mat_name)
    for mat_name in mesh_materials:
      poly.materials.append(bpy.data.materials[mat_name])

    poly.faces.add(len(faces))
    for j, (verts, mat_name) in enumerate(faces):
      f = poly.faces[j]
      f.vertices_raw = RawVertices(verts)
      f.material_index = mesh_materials.index(mat_name)

    poly.update()
    obj = bpy.data.objects.new(name, poly)
    bpy.context.scene.objects.link(obj)

def ReadSolutions():
  solutions = eval(file('pysolutions.txt').read())
  frames_per_solution = 10

  for i, solution in enumerate(solutions + solutions[:1]):
    for name, (euler, location) in solution:
      props = {'rotation_euler': euler, 'location': location }

      obj = bpy.data.objects[name]

      anim_data = obj.animation_data or obj.animation_data_create()
      action = anim_data.action
      if action is None:
        action = bpy.data.actions.new('%s_Action' % name)
        for prop in 'location', 'rotation_euler':
          for index in range(3):
            action.fcurves.new(prop, index)
        anim_data.action = action

      for curve in action.fcurves:
        value = props[curve.data_path][curve.array_index]
        curve.keyframe_points.insert(i * frames_per_solution, value)

  bpy.context.scene.frame_end = i * frames_per_solution

if __name__ == '__main__':
  ReadShapes()
  ReadSolutions()

... And if you've read this far, here's a little reward - a crude Blender render of the resulting animation. Nothing fancy yet, but it represents progress.

Throwing Shapes

The representation of solutions in the Cube module having been reworked, it's now possible to export data about the shapes and their positions in each solution, in a form that Blender can then load via a short Python script.

Here's the Haskell side of things:

> module Mesh where

> import Data.Map (fromList, keys, lookupIndex)
> import Data.Maybe (catMaybes, fromJust, listToMaybe)

> import System.Process (runCommand)

> import Cube

All the logic for moving a shape around within the cube has equivalents for moving a vertex or face around inside the cube. This will be useful for determining which faces of a shape end up on which faces of the cube in a particular solution.

> vApplyR :: Rotation -> (Int, Int, Int) -> (Int, Int, Int)
> vApplyR (Rotation a n) = times n (r a)
>  where
>   r X (x, y, z) = (x, 3-z, y)
>   r Y (x, y, z) = (z, y, 3-x)
>   r Z (x, y, z) = (3-y, x, z)

> vApplyT :: Translation -> (Int, Int, Int) -> (Int, Int, Int)
> vApplyT (Translation (x', y', z')) (x, y, z) = (x+x', y+y', z+z')

> vApplyRecipe (maj, min, shift) = vApplyR maj . vApplyR min . vApplyT shift

In effect, this function scans all the possible face positions, and determines where the cube is occupied by a form on one side along the X axis but not on the other. These will be the places where the form will need an X-face (i.e. one parallel to the YZ plane). (Locations for the Y- and Z-faces can be found by rotating the form appropriately first.)

> faceCoords :: Form -> [(Int, Int, Int)]
> faceCoords form = [(x, y, z) |
>                    (z, layer)   <- zip [0..] (unform form),
>                    (y, row)     <- zip [0..] layer,
>                    (x, isFace)  <- zip [0..] (markFaces row),
>                    isFace]
>  where markFaces row = zipWith (/=) (' ':row) (row ++ [' '])

A few new data types are needed, to record which side of the solved cube a face belongs to. The 'Show' instances for these have been carefully chosen to make parsing in Python easy.

> data Sign = Negative | Positive deriving (Show)
> data CubeFace = CubeFace Sign Axis
> instance Show CubeFace
>  where show (CubeFace sign axis) = [head $ show axis, head $ show sign]
> newtype FacePaint = FacePaint (Maybe CubeFace)
> instance Show FacePaint
>  where show (FacePaint (Just (CubeFace sign axis))) = [
>         '\"', head $ show axis, head $ show sign, '\"']
>        show (FacePaint Nothing) = "\"xx\""

> type Vert = (Int, Int, Int)

> type Face = (Int, Int, Int)

Now everything necessary is present, to be able to determine the locations of a form's faces for any axis. Note that applyR and vApplyR are used to represent a transformation and its inverse, respectively.

> axisFaces :: Axis -> Form -> [Face]
> axisFaces axis = map (vApplyR $ r axis 3) .
>                  faceCoords .
>                  applyR (r axis 1)
>  where
>   r X _ = Rotation X 0
>   r Y n = Rotation Z n
>   r Z n = Rotation Y n

Given a face on an axis, this function determines the locations of its vertices.

> faceVerts axis (x, y, z) = case axis of
>                              X -> [(x, y+y', z+z') | (y', z') <- square]
>                              Y -> [(x+x', y, z+z') | (x', z') <- square]
>                              Z -> [(x-x', y+y', z) | (x', y') <- square]
>  where square = [(0,0), (0,1), (1,1), (1,0)]

Knowing where a face's vertices lie (after transformation) is enough to determine which face of the cube (if any) it lies on.

> paintFace verts = listToMaybe . catMaybes $
>                   [paint vs axis | (vs, axis) <- [(xs, X), (ys, Y), (zs, Z)]]
>  where
>   (xs, ys, zs) = unzip3 verts
>   paint vs axis = if all (== 0) vs then Just $ CubeFace Negative axis
>                   else if all (==3) vs then Just $ CubeFace Positive axis
>                   else Nothing

Finally, all these things are brought together to give a way to describe how to position and colour all of a shape's faces to suit the space it occupies in a solution.

> axisInfo :: Recipe -> Axis -> Form -> [([Vert], FacePaint)]
> axisInfo recipe axis form = [(faceVerts axis face,
>                               FacePaint $ paintFace $ map (vApplyRecipe recipe) $ faceVerts axis face) |
>                              face <- axisFaces axis form]

> allFaces recipe f = concat [axisInfo recipe axis f | axis <- enumFrom X]

This little helper takes a list of faces (each having a list of vertices) and pulls all the vertices out into a separate (deduplicated) list, replacing them in the original data structure by indices into that vertex list.

> shareVerts :: Ord k => [([k], t)] -> ([k], [([Int], t)])
> shareVerts faces = (keys vertMap,
>                     [([fromJust (lookupIndex v vertMap) | v <- verts], x) |
>                      (verts, x) <- faces])
>  where vertMap = fromList [(v, ()) | v <- concat (map fst faces)]

Because the colouring of a shape depends on the choice of solution, the mesh can't be generated without specifying a recipe for the shape.

> mesh recipe = shareVerts . allFaces recipe . defaultForm
> allMeshes (Solution fs) = [(show shape, mesh recipe shape) | (shape, recipe, _) <- fs]

The file emitted here is one big list of mesh data. Luckily, 'show' for Haskell's built-in types and 'eval' in Python are near enough inverses that, with some careful implementations of 'Show' (as noted above) it's possible to avoid writing a parser for it at all.

> writeMeshes = writeFile "shapes.txt" $ show $ allMeshes (head solutions)

This rather ugly function describes the two rotations and a translation that compose a Recipe, in a form that can easily be parsed in Python.

> pyRecipe recipe@(maj, min, shift) = ((roll, pitch, heading), location)
>  where
>   location = vApplyRecipe recipe (0, 0, 0)
>   heading =
>     case maj of Rotation X _ -> error "X should only appear as a minor axis"
>                 Rotation Y _ -> 0
>                 Rotation Z n -> quarterTurns n
>   pitch =
>     case maj of Rotation X _ -> error "X should only appear as a minor axis"
>                 Rotation Y n -> quarterTurns n
>                 Rotation Z _ -> 0
>   roll =
>     case min of Rotation X n -> quarterTurns n
>                 _            -> error "Only X should appear as a minor axis"
>   quarterTurns n = (fromIntegral n) * pi/2

> pySolution (Solution fs) = [(show shape, pyRecipe recipe) |
>                             (shape, recipe, _) <- fs]

> writeSolutions = writeFile "pysolutions.txt" $
>                  show [pySolution solution | solution <- solutions]

Lastly, a test command glues everything together: writing out the data, then loading it into Blender.

> testMeshes = do
>   writeMeshes
>   writeSolutions
>   _ <- runCommand $ unwords ["blender", "-P", "blender.py"]
>   return ()

Tedious Groundwork

I've been at the cube again. Quite a lot has developed, but it's all been built upon some relatively small but important changes to the way that the original cube solving code models its work. Specifically, each solution now carries with it not only information about which shapes occupy which cells, but also information about how to translate and rotate each piece from its default position into its position within the solution.

This is crucial for tasks such as exporting solution geometry to a 3D modelling package such as Blender ... of which more later. For the moment, here's a very boring rehash of the existing code, for completeness.

> module Cube where

> import Data.Bits
> import Data.Function (on)
> import Data.List (intercalate, nubBy, sortBy, transpose, unfoldr)
> import Data.Monoid
> import Data.Ord (comparing)
> import Data.Word

The shapes remain the same. The last one has been slightly renamed to avoid conflict with the axis of the same name.

> data Shape = L | S | T | R | P | Q | Y'
>   deriving (Enum, Show)

> plan L = "XXX" ++
>          "X  "

> plan S = " XX" ++
>          "XX "

> plan T = "XXX" ++
>          " X "

> plan R = "XX " ++
>          "X  "

> plan P = "X  " ++
>          "X  " ++
>          "   " ++
>          "XX "

> plan Q = " X " ++
>          " X " ++
>          "   " ++
>          "XX "

> plan Y' = "XX " ++
>           "X  " ++
>           "   " ++
>           "X  "

> newtype Form = Form [[[Char]]]
> unform (Form x) = x

> instance Eq Form where
>   (Form x) == (Form y) = x == y

> instance Ord Form where
>   Form x `compare` Form y = x `compare` y

> withForm f = Form . f . unform

The order in which rows of a form are displayed has been reversed, so that Y increases up the screen. With X increasing to the right, this means that each layer of the cube is displayed as if viewed from above (looking towards -Z) in a right-handed coordinate system, so it fits Blender naturally.

> instance Show Form where
>   show (Form f) = unlines [showSlice s | s <- reverse $ transpose f]
>    where
>     showSlice s = intercalate "  " [ showRow r | r <- s ]
>     showRow r = "|" ++ r ++ "|"

> instance Monoid Form where
>   mempty = Form $ (replicate 3 . replicate 3 . replicate 3) ' '
>   mappend (Form x) (Form y) = Form $ (zipWith . zipWith . zipWith)
>                                      superpose x y
>    where
>     superpose ' ' ' ' = ' '
>     superpose  a  ' ' =  a
>     superpose ' '  b  =  b
>     superpose  _   _  = '#'  -- This should never happen!

The following splits a list into chunks, each of which (except maybe the last) has length 3. This would naturally generalise to chunks of length N if needed. This was previously defined within 'defaultForm' but is now more widely useful.

> inThrees = takeWhile (not . null) . unfoldr (Just . splitAt 3)

> flatten = concat . concat . unform

As well as flattening a form into a list, the inverse operation will also be necessary. This takes a flat list and produces a list of triples of triples.

> unflatten = inThrees . inThrees

> asBits :: Form -> Word32
> asBits f = sum [bit n | (n, x) <- zip [0..] (flatten f), x /= ' ']

The opposite of 'asBits', this takes a shape and a bit sequence representing its form, and reconstructs the form. The ability to do this means that computations can take place purely in terms of bit sequences, and the results can still be displayed in a comprehensible manner by expanding them afterwards.

> fromBits shape bs = Form $ unflatten [formChar (testBit bs n) | n <- [0..26]]
>  where
>   formChar False = ' '
>   formChar True  = head (show shape)

> size f = length [x | x <- flatten f, x /= ' ']

> defaultForm :: Shape -> Form
> defaultForm shape = Form $ unflatten fullPlan
>  where
>   fullPlan = [formChar c | c <- plan shape] ++
>              replicate (27 - length (plan shape)) ' '
>   formChar ' ' = ' '
>   formChar  _  = head (show shape)

To be able to produce geometrical descriptions of the solutions, it is useful to record not only which shapes occupy which cell of the cube, but also how the shapes' forms are composed. As mentioned before, each form can be constructed from a translation, followed by a rotation about the X axis and then a rotation about either the Y or Z axes.

These transformations could all be represented as matrices, but they're much easier to decipher if given a higher-level representation. This sequence of three transformations is referred to as a recipe.

> data Axis = X | Y | Z deriving (Enum, Show)
> 
> newtype Translation = Translation (Int, Int, Int) deriving Show
> 
> type QuarterTurns = Int
> data Rotation = Rotation Axis QuarterTurns deriving Show
> 
> type Recipe = (Rotation, Rotation, Translation)

Having given names to these transformations, they can now be listed in terms of those names instead of relying directly on function composition.

> majors = [Rotation Y n | n <- [0..3]] ++
>          [Rotation Z n | n <- [1, 3]]
> 
> minors = [Rotation X n | n <- [0..3]]
> 
> shifts shape = [t | t <- [Translation (x, y, z) |
>                           x <- [0..2],
>                           y <- [0..2],
>                           z <- [0..2]], inBounds t]
>  where
>   inBounds t = size (applyT t $ df) == size df
>   df = defaultForm shape

Translations can now be in either direction; this might prove useful later.

> shiftr def [a, b, _] = [def, a, b]
> shiftl def [_, b, c] = [b, c, def]

> withx op = (map . map) (op ' ')
> withy op = map (op "   ")
> withz op = op (replicate 3 "   ")
> tx  = withx shiftr
> tx' = withx shiftl
> ty  = withy shiftr
> ty' = withy shiftl
> tz  = withz shiftr
> tz' = withz shiftl

> cw = transpose . reverse
> ccw = reverse . transpose

Transformations being defined in terms of Recipes instead of raw functions, a little glue is needed to describe how a recipe is applied to a Form.

Note also that the rotation about X has been changed from clockwise to counter-clockwise. These rotations now match Blender's default settings, which will avoid complication further down the line.

> times 0 f = id
> times n f = f . times (n-1) f
> 
> applyR (Rotation a n) = withForm (times n (r a))
>  where
>   r X = cw
>   r Y = transpose . map ccw . transpose
>   r Z = map cw
> 
> applyT (Translation (x, y, z)) = withForm (times x tx . times y ty . times z tz)
> 
> applyRecipe (maj, min, shift) = applyR maj . applyR min . applyT shift

The special case for the forms of L have been described a little more neatly.

> allForms :: Shape -> [(Recipe, Form)]
> allForms shape = nubBy ((==) `on` snd) $ [(recipe, applyRecipe recipe df) |
>                                           recipe <- recipes shape]
>  where
>   df = defaultForm shape
>   recipes L = [(Rotation Y 0, Rotation X 0, t) |
>                t <- [s | s@(Translation (a, b, c)) <- shifts L, c /= 2]]
>   recipes shape =  [(maj, min, shift) |
>                     maj <- majors,
>                     min <- minors,
>                     shift <- shifts shape]

The list of possible forms for each shape can now carry the recipe for constructing each form, as well as the resulting bit pattern.

> combos :: [(Shape, [(Recipe, Word32)])]
> combos = sortBy (comparing (length . snd))
>                 [(shape, [(r, asBits f) | (r, f) <- allForms shape]) |
>                  shape <- enumFrom L]

> newtype Solution = Solution [(Shape, Recipe, Word32)] deriving (Show)

> solutions :: [Solution]
> solutions = [s | (s, _) <- foldl addForm [(Solution [], 0)] combos] where
>   addForm partialSolutions (shape, formData) = [
>       (Solution ((shape, recipe, bitmap):s), (pb .|. bitmap)) |
>       (Solution s, pb) <- partialSolutions,
>       (recipe, bitmap) <- formData,
>       pb .&. bitmap == 0
>       ]

> expandSolution (Solution s) = foldl1 mappend [fromBits shape bitmap |
>                                               (shape, _, bitmap) <- s]

Finally, this makes it possible to explain what sequence of rotations and transformations was used to place each shape in position in a solution.

> explainSolution (Solution s) = unlines [show (shape, recipe) |
>                                         (shape, recipe, _) <- s]

> main = putStrLn $ intercalate "\n" [show (expandSolution s) | s <- solutions]