Below we present some animations that illustrate operations on finite patches of Penrose’s Kite and Dart tiles.

These were created using PenroseKiteDart which is a Haskell package available on Hackage making use of the Haskell Diagrams package. For details, see the PenroseKiteDart user guide.

Penrose’s Kite and Dart tiles can produce infinite aperiodic tilings of the plane. There are legal tiling rules to ensure aperiodicity, but these rules do not guarantee that a finite tiling will not get stuck. A legal finite tiling which can be continued to cover the whole plane is called a correct tiling. The rest, which are doomed to get stuck, are called incorrect tilings. (More details can be found in the links at the end of this blog.)

Decomposition Animations

The function decompose is a total operation which is guaranteed to preserve the correctness of a finite tiling represented as a tile graph (or Tgraph). Let us start with a particular Tgraph called sunGraph which is defined in PenroseKiteDart and consists of 5 kites arranged with a common origin vertex. It is drawn using default style in figure 1 on the left. On the right of figure 1 it is drawn with both vertex labels and dotted lines for half-tile join edges.

We can decompose sunGraph three times by selecting index 3 of the infinite list of its decompositions.

    sunD3 :: Tgraph
    sunD3 = decompositions sunGraph !! 3

where we have used

    decompose :: Tgraph -> Tgraph
    
    decompositions :: Tgraph -> [Tgraph]
    decompositions = iterate decompose

The result (sunD3) is drawn in figure 2 (scaled up).

The animation in figure 3 illustrates two further decompositions of sunD3 in two stages.

Figure 4 also illustrates two decompositions, this time starting from forcedKingD.

    forcedKingD :: Tgraph
    forcedKingD = force (decompose kingGraph)

A Composition Animation

An inverse to decomposing (namely composing) has some extra intricacies. In the literature (see for example 1 and 2) versions of the following method are frequently described.

• Firstly, split darts in half.
• Secondly, glue all the short edges of the half-darts where they meet a kite (simultaneously). This will form larger scale complete darts and larger scale half kites.
• Finally join the halves of the larger scale kites.

This works for infinite tilings, but we showed in Graphs,Kites and Darts and Theorems that this method is unsound for finite tilings. There is the trivial problem that a half-dart may not have a complete kite on its short edge. Worse still, the second step can convert a correct finite tiling into an incorrect larger scale tiling. An example of this is given in Graphs, Kites and Darts and Theorems where we also described our own safe method of composing (never producing an incorrect Tgraph when given a correct Tgraph). This composition can leave some boundary half-tiles out of the composition (called remainder half-tiles).

The animation in figure 5 shows such a composition where the remainder half-tiles are indicated with lime green edges.

In general, compose is a partial operation as the resulting half-tiles can break some requirements for Tgraphs (namely, connectedness and no crossing boundaries). However we have shown that it is a total function on forced Tgraphs. (Forcing is discussed next.)

Forcing Animations

The process of forcing a Tgraph adds half-tiles on the boundary where only one legal choice is possible. This continues until either there are no more forced additions possible, or a clash is found showing that the tiling is incorrect. In the latter case it must follow that the initial tiling before forcing was already an incorrect tiling.

The process of forcing is animated in figure 6, starting with a 5 times decomposed kite and in figure 7 with a 5 times decomposed dart.

It is natural to wonder what forcing will do with cut-down (but still correct) Tgraphs. For example, taking just the boundary faces from the final Tgraph shown in the previous animation forms a valid Tgraph (boundaryExample) shown in figure 8.

    boundaryExample :: Tgraph
    boundaryExample = runTry $ tryBoundaryFaceGraph $ force $ decompositions dartGraph !!5

Applying force to boundaryExample just fills in the hole to recreate force (decompositions dartGraph !!5) modulo vertex numbering. To make it more interesting we tried removing further half-tiles from boundaryExample to make a small gap. Forcing this also completes the filling in of the boundary half-tiles to recreate force (decompositions dartGraph !!5). However, we can see that this filling in is constrianed to preserve the required Tgraph property of no crossing boundaries which prevents the tiling closing round a hole.

This is illustrated in the animation shown in figure 9.

As another experiment, we take the boundary faces of a (five times decomposed but not forced) star. When forced this fills in the star and also expands outwards, as illustrated in figure 10.

In the final example, we pick out a shape within a correct Tgraph (ensuring the chosen half-tiles form a valid Tgraph) then animate the force process and then run the animation in both directions (by adding a copy of the frames in reverse order).

The result is shown in figure 11.

Creating Animations

Animations as gif files can be produced by the Haskell Diagrams package using the rasterific back end.

The main module should import both Diagrams.Prelude and Diagrams.Backend.Rasterific.CmdLine. This will expose the type B standing for the imported backend, and diagrams then have type Diagram B.

An animation should have type [(Diagram B, Int)] and consist of a list of frames for the animation, each paired with an integer delay (in one-hundredths of a second).

The animation can then be passed to mainWith.

module Main (main) where
    
import Diagrams.Prelude
import Diagrams.Backend.Rasterific.CmdLine

...

fig::[(Diagram B,Int)]
fig = myExampleAnimation

main :: IO ()
main = mainWith fig

If main is then compiled and run (e.g. with parameters -w 700 -o test.gif) it will produce an output file (test.gif with width 700).

Crossfade tool

The decompose and compose animations were defined using crossfade.

crossfade :: Int -> Diagram B -> Diagram B -> [Diagram B]
crossfade n d1 d2 = map blending ratios 
  where
    blending r = opacity (1-r) d1 <> opacity r d2
    ratios = map ((/ fromIntegral n) . fromIntegral) [0..n]

Thus crossfade n d1 d2 produces n+1 frames, each with d1 overlaid on d2 but with varying opacities (decreasing for d1 and increasing for d2).

Adding the same pause (say 10 hundreths of a second) to every frame can be done by applying map (,10) and this will produce an animation.

Force animation tool

To create force animations it was useful to create a tool to produce frames with stages of forcing.

forceFrames :: Angle Double 
            -> Int
            -> Tgraph 
            -> (Colour Double, Colour Double, Colour Double)
            -> [Diagram B]

This takes as arguments

• an angle argument (to rotate the diagrams in the animation from the default alignment of the Tgraph),
• an Int (for the required number of frames),
• a Tgraph (to be forced),
• a triple of colours for filling darts, kites and grout (edge colour), respectively.

The definition of forceFrames uses stepForce to advance forcing a given number of steps to get the intermediate Tgraphs. The total number of forcing steps will be the number of faces (half-tiles) in the final force g less the number of faces in the initial g. All the Tgraphs are drawn (using colourDKG) but the resulting diagrams must all be aligned properly. The alignment can be achieved by creating a VPatch (vertex patch) from the final Tgraph which is then rotated. All the Tgraphs can then be drawn using sub vertex patches of the final rotated one. (For details see Overlaid examples in the PenroseKiteDart user guide.)

Previous related blogs

• PenroseKiteDart user guide – this explains how to install and use the PenroseKiteDart package.
• Graphs,Kites and Darts and Theorems established some important results relating force, compose, decompose.
• Empires and SuperForce – these new operations were based on observing properties of boundaries of forced Tgraphs.
• Graphs, Kites and Darts introduced Tgraphs. This gave more details of implementation and results of early explorations. (The class Forcible was introduced subsequently).
• Diagrams for Penrose Tiles – the first blog introduced drawing Pieces and Patches (without using Tgraphs) and provided a version of decomposing for Patches (decompPatch).

References

[1] Martin Gardner (1977) MATHEMATICAL GAMES. Scientific American, 236(1), (pages 110 to 121). http://www.jstor.org/stable/24953856

[2] Grünbaum B., Shephard G.C. (1987) Tilings and Patterns. W. H. Freeman and Company, New York. ISBN 0-7167-1193-1 (Hardback) (pages 540 to 542).

(Updated February 2026 for PenroseKiteDart version 1.6.1)

PenroseKiteDart is a Haskell package with tools to experiment with finite tilings of Penrose’s Kites and Darts. It uses the Haskell Diagrams package for drawing tilings. As well as providing drawing tools, this package introduces tile graphs (Tgraphs) for describing finite tilings. (I would like to thank Stephen Huggett for suggesting planar graphs as a way to reperesent the tilings).

This document summarises the design and use of the PenroseKiteDart package.

PenroseKiteDart package is now available on Hackage.

The source files are available on GitHub at https://github.com/chrisreade/PenroseKiteDart.

There is a small art gallery of examples created with PenroseKiteDart here.

Index

1. About Penrose’s Kites and Darts
2. Using the PenroseKiteDart Package (initial set up).
3. Overview of Types and Operations
4. Drawing in more detail
5. Forcing in more detail
6. Advanced Operations
7. Other Reading

1. About Penrose’s Kites and Darts

The Tiles

In figure 1 we show a dart and a kite. All angles are multiples of (a tenth of a full turn). If the shorter edges are of length 1, then the longer edges are of length , where is the golden ratio.

Aperiodic Infinite Tilings

What is interesting about these tiles is:

It is possible to tile the entire plane with kites and darts in an aperiodic way.

Such a tiling is non-periodic and does not contain arbitrarily large periodic regions or patches.

The possibility of aperiodic tilings with kites and darts was discovered by Sir Roger Penrose in 1974. There are other shapes with this property, including a chiral aperiodic monotile discovered in 2023 by Smith, Myers, Kaplan, Goodman-Strauss. (See the Penrose Tiling Wikipedia page for the history of aperiodic tilings)

This package is entirely concerned with Penrose’s kite and dart tilings also known as P2 tilings.

Legal Tilings

In figure 2 we add a temporary green line marking purely to illustrate a rule for making legal tilings. The purpose of the rule is to exclude the possibility of periodic tilings.

If all tiles are marked as shown, then whenever tiles come together at a point, they must all be marked or must all be unmarked at that meeting point. So, for example, each long edge of a kite can be placed legally on only one of the two long edges of a dart. The kite wing vertex (which is marked) has to go next to the dart tip vertex (which is marked) and cannot go next to the dart wing vertex (which is unmarked) for a legal tiling.

Correct Tilings

Unfortunately, having a finite legal tiling is not enough to guarantee you can continue the tiling without getting stuck. Finite legal tilings which can be continued to cover the entire plane are called correct and the others (which are doomed to get stuck) are called incorrect. This means that decomposition and forcing (described later) become important tools for constructing correct finite tilings.

2. Using the PenroseKiteDart Package

You will need the Haskell Diagrams package (See Haskell Diagrams) as well as this package (PenroseKiteDart). When these are installed, you can produce diagrams with a Main.hs module. This should import a chosen backend for diagrams such as the default (SVG) along with Diagrams.Prelude.

    module Main (main) where
    
    import Diagrams.Backend.SVG.CmdLine
    import Diagrams.Prelude

For Penrose’s Kite and Dart tilings, you also need to import the PKD module and (optionally) the TgraphExamples module.

    import PKD
    import TgraphExamples

Then to ouput someExample figure

    fig::Diagram B
    fig = someExample

    main :: IO ()
    main = mainWith fig

Note that the token B is used in the diagrams package to represent the chosen backend for output. So a diagram has type Diagram B. In this case B is bound to SVG by the import of the SVG backend. When the compiled module is executed it will generate an SVG file. (See Haskell Diagrams for more details on producing diagrams and using alternative backends).

3. Overview of Types and Operations

Half-Tiles

In order to implement operations on tilings (decompose in particular), we work with half-tiles. These are illustrated in figure 3 and labelled RD (right dart), LD (left dart), LK (left kite), RK (right kite). The join edges where left and right halves come together are shown with dotted lines, leaving one short edge and one long edge on each half-tile (excluding the join edge). We have shown a red dot at the vertex we regard as the origin of each half-tile (the tip of a half-dart and the base of a half-kite).

The labels are actually data constructors introduced with type operator HalfTile which has an argument type (rep) to allow for more than one representation of the half-tiles.

    data HalfTile rep 
      = LD rep -- Left Dart
      | RD rep -- Right Dart
      | LK rep -- Left Kite
      | RK rep -- Right Kite
      deriving (Show,Eq)

Tgraphs

We introduce tile graphs (Tgraphs) which provide a simple planar graph representation for finite patches of tiles. For Tgraphs we first specialise HalfTile with a triple of vertices (positive integers) to make a TileFace such as RD(1,2,3), where the vertices go clockwise round the half-tile triangle starting with the origin.

    type TileFace  = HalfTile (Vertex,Vertex,Vertex)
    type Vertex    = Int  -- must be positive

The function

    makeTgraph :: [TileFace] -> Tgraph

then constructs a Tgraph from a TileFace list after checking the TileFaces satisfy certain properties (described below). We also have

    faces :: Tgraph -> [TileFace]

to retrieve the TileFace list from a Tgraph.

As an example, the fool (short for fool’s kite and also called an ace in the literature) consists of two kites and a dart (= 4 half-kites and 2 half-darts):

    fool :: Tgraph
    fool = makeTgraph [RD (1,2,3), LD (1,3,4)   -- right and left dart
                      ,LK (5,3,2), RK (5,2,7)   -- left and right kite
                      ,RK (5,4,3), LK (5,6,4)   -- right and left kite
                      ]

To produce a diagram, we simply draw the Tgraph

    foolFigure :: Diagram B
    foolFigure = draw fool

which will produce the diagram on the left in figure 4.

Alternatively,

    foolFigure :: Diagram B
    foolFigure = labelled drawj fool

will produce the diagram on the right in figure 4 (showing vertex labels and dashed join edges).

When any (non-empty) Tgraph is drawn, a default orientation and scale are chosen based on the lowest numbered join edge. This is aligned on the positive x-axis with length 1 (for darts) or length (for kites).

Tgraph Properties

Tgraphs are actually implemented as

    newtype Tgraph = Tgraph [TileFace]
                     deriving (Show)

but the data constructor Tgraph is not exported to avoid accidentally by-passing checks for the required properties. The properties checked by makeTgraph ensure the Tgraph represents a legal tiling as a planar graph with positive vertex numbers, and that the collection of half-tile faces are both connected and have no crossing boundaries (see note below). Finally, there is a check to ensure two or more distinct vertex numbers are not used to represent the same vertex of the graph (a touching vertex check). An error is raised if there is a problem.

Note: If the TileFaces are faces of a planar graph there will also be exterior (untiled) regions, and in graph theory these would also be called faces of the graph. To avoid confusion, we will refer to these only as exterior regions, and unless otherwise stated, face will mean a TileFace. We can then define the boundary of a list of TileFaces as the edges of the exterior regions. There is a crossing boundary if the boundary crosses itself at a vertex. We exclude crossing boundaries from Tgraphs because they prevent us from calculating relative positions of tiles locally and create touching vertex problems.

For convenience, in addition to makeTgraph, we also have

    makeUncheckedTgraph :: [TileFace] -> Tgraph
    checkedTgraph   :: [TileFace] -> Tgraph

The first of these (performing no checks) is useful when you know the required properties hold. The second performs the same checks as makeTgraph except that it omits the touching vertex check. This could be used, for example, when making a Tgraph from a sub-collection of TileFaces of another Tgraph.

Main Tiling Operations

There are three key operations on finite tilings, namely

    decompose :: Tgraph -> Tgraph
    force     :: Tgraph -> Tgraph
    compose   :: Tgraph -> Tgraph

Decompose

Decomposition (also called deflation) works by splitting each half-tile into either 2 or 3 new (smaller scale) half-tiles, to produce a new tiling. The fact that this is possible, is used to establish the existence of infinite aperiodic tilings with kites and darts. Since our Tgraphs have abstracted away from scale, the result of decomposing a Tgraph is just another Tgraph. However if we wish to compare before and after with a drawing, the latter should be scaled by a factor times the scale of the former, to reflect the change in scale.

We can, of course, iterate decompose to produce an infinite list of finer and finer decompositions of a Tgraph

    decompositions :: Tgraph -> [Tgraph]
    decompositions = iterate decompose

Force

Force works by adding any TileFaces on the boundary edges of a Tgraph which are forced. That is, where there is only one legal choice of TileFace addition consistent with the seven possible vertex types. Such additions are continued until either (i) there are no more forced cases, in which case a final (forced) Tgraph is returned, or (ii) the process finds the tiling is stuck, in which case an error is raised indicating an incorrect tiling. [In the latter case, the argument to force must have been an incorrect tiling, because the forced additions cannot produce an incorrect tiling starting from a correct tiling.]

An example is shown in figure 6. When forced, the Tgraph on the left produces the result on the right. The original is highlighted in red in the result to show what has been added.

Compose

Composition (also called inflation) is an opposite to decompose but this has complications for finite tilings, so it is not simply an inverse. (See Graphs,Kites and Darts and Theorems for more discussion of the problems). Figure 7 shows a Tgraph (left) with the result of composing (right) where we have also shown (in pale green) the faces of the original that are not included in the composition – the remainder faces.

Under some circumstances composing can fail to produce a Tgraph because there are crossing boundaries in the resulting TileFaces. However, we have established that

• If g is a forced Tgraph, then compose g is defined and it is also a forced Tgraph.

Try Results

It is convenient to use types of the form Try a for results where we know there can be a failure. For example, compose can fail if the result does not pass the connected and no crossing boundary check, and force can fail if its argument is an incorrect Tgraph. In situations when you would like to continue some computation rather than raise an error when there is a failure, use a try version of a function.

    tryCompose :: Tgraph -> Try Tgraph
    tryForce   :: Tgraph -> Try Tgraph

We define Try as a synonym for Either ShowS (which is a monad) in module Tgraph.Try.

type Try a = Either ShowS a

(Note ShowS is String -> String). Successful results have the form Right r (for some correct result r) and failure results have the form Left (s<>) (where s is a String describing the problem as a failure report).

The function

    runTry:: Try a -> a
    runTry = either error id

will retrieve a correct result but raise an error for failure cases. This means we can always derive an error raising version from a try version of a function by composing with runTry.

    force = runTry . tryForce
    compose = runTry . tryCompose

Elementary Tgraph and TileFace Operations

The module Tgraph.Prelude defines elementary operations on Tgraphs relating vertices, directed edges, and faces. We describe a few of them here.

When we need to refer to particular vertices of a TileFace we use

    originV :: TileFace -> Vertex -- the first vertex - red dot in figure 2
    oppV    :: TileFace -> Vertex -- the vertex at the opposite end of the join edge from the origin
    wingV   :: TileFace -> Vertex -- the vertex not on the join edge

A directed edge is represented as a pair of vertices.

    type Dedge = (Vertex,Vertex)

So (a,b) is regarded as a directed edge from a to b.

When we need to refer to particular edges of a TileFace we use

    joinE  :: TileFace -> Dedge  -- shown dotted in figure 2
    shortE :: TileFace -> Dedge  -- the non-join short edge
    longE  :: TileFace -> Dedge  -- the non-join long edge

which are all directed clockwise round the TileFace. In contrast, joinOfTile is always directed away from the origin vertex, so is not clockwise for right darts or for left kites:

    joinOfTile:: TileFace -> Dedge
    joinOfTile face = (originV face, oppV face)

In the special case that a list of directed edges is symmetrically closed [(b,a) is in the list whenever (a,b) is in the list] we can think of this as an edge list rather than just a directed edge list.

For example,

    internalEdges :: Tgraph -> [Dedge]

produces an edge list, whereas

    boundary :: Tgraph -> [Dedge]

produces single directions. Each directed edge in the resulting boundary will have a TileFace on the left and an exterior region on the right. The function

    dedges :: Tgraph -> [Dedge]

produces all the directed edges obtained by going clockwise round each TileFace so not every edge in the list has an inverse in the list.

Note: There is now a class HasFaces (introduced in version 1.4) which includes instances for both Tgraph and [TileFace] and others. This allows some generalisations. For example

    faces         :: HasFaces a => a -> [TileFace]
    internalEdges :: HasFaces a => a -> [Dedge]
    boundary      :: HasFaces a => a -> [Dedge] 
    dedges        :: HasFaces a => a -> [Dedge] 
    nullFaces     :: HasFaces a => a -> Bool

Patches (Scaled and Positioned Tilings)

Behind the scenes, when a Tgraph is drawn, each TileFace is converted to a Piece. A Piece is another specialisation of HalfTile using a two dimensional vector to indicate the length and direction of the join edge of the half-tile (from the originV to the oppV), thus fixing its scale and orientation. The whole Tgraph then becomes a list of located Pieces called a Patch.

    type Piece = HalfTile (V2 Double)
    type Patch = [Located Piece]

Piece drawing functions derive vectors for other edges of a half-tile piece from its join edge vector. In particular (in the TileLib module) we have

    drawPiece :: Piece -> Diagram B
    darawjPiece :: Piece -> Diagram B
    fillPieceDK :: Colour Double -> Colour Double -> Piece -> Diagram B

where the first draws the non-join edges of a Piece, the second does the same but adds a faint dashed line for the join edge, and the third takes two colours – one for darts and one for kites, which are used to fill the piece as well as using drawPiece.

Patch is an instance of class Transformable so a Patch can be scaled, rotated, and translated.

Vertex Patches

It is useful to have an intermediate form between Tgraphs and Patches, that contains information about both the location of vertices (as 2D points), and the abstract TileFaces. This allows us to introduce labelled drawing functions (to show the vertex labels) which we then extend to Tgraphs. We call the intermediate form a VPatch (short for Vertex Patch).

    type VertexLocMap = IntMap.IntMap (Point V2 Double)
    data VPatch = VPatch {vLocs :: VertexLocMap,  vpFaces::[TileFace]} deriving Show

and

    makeVP :: Tgraph -> VPatch

calculates vertex locations using a default orientation and scale.

VPatch is made an instance of class Transformable so a VPatch can also be scaled and rotated.

One essential use of this intermediate form is to be able to draw a Tgraph with labels, rotated but without the labels themselves being rotated. We can simply convert the Tgraph to a VPatch, and rotate that before drawing with labels.

    labelled draw (rotate someAngle (makeVP g))

We can also align a VPatch using vertex labels.

    alignXaxis :: (Vertex, Vertex) -> VPatch -> VPatch 

So if g is a Tgraph with vertex labels a and b we can align it on the x-axis with a at the origin and b on the positive x-axis (after converting to a VPatch), instead of accepting the default orientation.

    labelled draw (alignXaxis (a,b) (makeVP g))

Another use of VPatches is to share the vertex location map when drawing only subsets of the faces (see Overlaid examples in the next section).

4. Drawing in More Detail

Class Drawable

There is a class Drawable with instances Tgraph, VPatch, Patch. When the token B is in scope standing for a fixed backend then we can assume

    draw   :: Drawable a => a -> Diagram B  -- draws non-join edges
    drawj  :: Drawable a => a -> Diagram B  -- as with draw but also draws dashed join edges
    fillDK :: Drawable a => Colour Double -> Colour Double -> a -> Diagram B -- fills with colours

where fillDK clr1 clr2 will fill darts with colour clr1 and kites with colour clr2 as well as drawing non-join edges.

These are the main drawing tools. However they are actually defined for any suitable backend b so have more general types.

(Update Sept 2024) From version 1.1 onwards of PenroseKiteDart, these are

    draw ::   (Drawable a, OKBackend b) =>
              a -> Diagram b
    drawj ::  (Drawable a, OKBackend) b) =>
              a -> Diagram b
    fillDK :: (Drawable a, OKBackend b) =>
              Colour Double -> Colour Double -> a -> Diagram b

where the class OKBackend is a check to ensure a backend is suitable for drawing 2D tilings with or without labels.

In these notes we will generally use the simpler description of types using B for a fixed chosen backend for the sake of clarity.

The drawing tools are each defined via the class function drawWith using Piece drawing functions.

    class Drawable a where
        drawWith :: (Piece -> Diagram B) -> a -> Diagram B
    
    draw = drawWith drawPiece
    drawj = drawWith drawjPiece
    fillDK clr1 clr2 = drawWith (fillPieceDK clr1 clr2)

To design a new drawing function, you only need to implement a function to draw a Piece, (let us call it newPieceDraw)

    newPieceDraw :: Piece -> Diagram B

This can then be elevated to draw any Drawable (including Tgraphs, VPatches, and Patches) by applying the Drawable class function drawWith:

    newDraw :: Drawable a => a -> Diagram B
    newDraw = drawWith newPieceDraw

Class DrawableLabelled

Class DrawableLabelled is defined with instances Tgraph and VPatch, but Patch is not an instance (because this does not retain vertex label information).

    class DrawableLabelled a where
        labelColourSize :: Colour Double -> Measure Double -> (Patch -> Diagram B) -> a -> Diagram B

So labelColourSize c m modifies a Patch drawing function to add labels (of colour c and size measure m). Measure is defined in Diagrams.Prelude with pre-defined measures tiny, verySmall, small, normal, large, veryLarge, huge. For most of our diagrams of Tgraphs, we use red labels and we also find small is a good default size choice, so we define

    labelSize :: DrawableLabelled a => Measure Double -> (Patch -> Diagram B) -> a -> Diagram B
    labelSize = labelColourSize red

    labelled :: DrawableLabelled a => (Patch -> Diagram B) -> a -> Diagram B
    labelled = labelSize small

and then labelled draw, labelled drawj, labelled (fillDK clr1 clr2) can all be used on both Tgraphs and VPatches as well as (for example) labelSize tiny draw, or labelCoulourSize blue normal drawj.

Further drawing functions

There are a few extra drawing functions built on top of the above ones. The function smart is a modifier to add dashed join edges only when they occur on the boundary of a Tgraph

    smart :: (VPatch -> Diagram B) -> Tgraph -> Diagram B

So smart vpdraw g will draw dashed join edges on the boundary of g before applying the drawing function vpdraw to the VPatch for g. For example the following all draw dashed join edges only on the boundary for a Tgraph g

    smart draw g
    smart (labelled draw) g
    smart (labelSize normal draw) g

When using labels, the function rotating allows a Tgraph to be drawn rotated without rotating the labels.

    rotating :: Angle Double -> (VPatch -> a) -> Tgraph -> a
    rotating angle vpdraw = vpdraw . rotate angle . makeVP

So for example,

    rotating (90@@deg) (labelled draw) g

makes sense for a Tgraph g. Of course if there are no labels we can simply use

    rotate (90@@deg) (draw g)

Similarly aligning allows a Tgraph to be aligned on the X-axis using a pair of vertex numbers before drawing.

    aligning :: (Vertex,Vertex) -> (VPatch -> a) -> Tgraph -> a
    aligning (a,b) vpdraw = vpdraw . alignXaxis (a,b) . makeVP

So, for example, if Tgraph g has vertices a and b, both

    aligning (a,b) draw g
    aligning (a,b) (labelled draw) g

make sense. Note that the following two examples are wrong. Even though they type check, they re-orient g without repositioning the boundary joins.

    smart (labelled draw . rotate angle) g      -- WRONG
    smart (labelled draw . alignXaxis (a,b)) g  -- WRONG

Instead use

    smartRotating angle (labelled draw) g
    smartAligning (a,b) (labelled draw) g

where

    smartRotatinge :: Angle Double -> (VPatch -> Diagram B) -> Tgraph -> Diagram B
    smartAligning  :: (Vertex,Vertex) -> (VPatch -> Diagram B) -> Tgraph -> Diagram B

are defined using

    smartOn :: Tgraph -> (VPatch -> Diagram B) -> VPatch -> Diagram B

Here, smartOn g vpdraw vp uses the given vp for drawing boundary joins and drawing faces of g (with vpdraw) rather than converting g to a new VPatch. This assumes vp has locations for vertices in g.

Overlaid examples (location map sharing)

The function

    drawForce :: Tgraph -> Diagram B

will (smart) draw a Tgraph g in red overlaid (using <>) on the result of force g as in figure 6. Similarly

    drawPCompose  :: Tgraph -> Diagram B

applied to a Tgraph g will draw the result of a partial composition of g as in figure 7. That is a drawing of compose g but overlaid with a drawing of the remainder faces of g shown in pale green.

Both these functions make use of sharing a vertex location map to get correct alignments of overlaid diagrams. In the case of drawForce g, we know that a VPatch for force g will contain all the vertex locations for g since force only adds to a Tgraph (when it succeeds). So when constructing the diagram for g we can use the VPatch created for force g instead of starting afresh. Similarly for drawPCompose g the VPatch for g contains locations for all the vertices of compose g so compose g is drawn using the VPatch for g instead of starting afresh.

The location map sharing is done with

    subFaces :: HasFaces a => 
                a -> VPatch -> VPatch

so that subFaces fcs vp is a VPatch with the same vertex locations as vp, but replacing the faces of vp with fcs. [Of course, this can go wrong if the new faces have vertices not in the domain of the vertex location map so this needs to be used with care. Any errors would only be discovered when a diagram is created.]

For cases where labels are only going to be drawn for certain faces, we need a version of subFaces which also gets rid of vertex locations that are not relevant to the faces. For this situation we have

    restrictTo:: HasFaces a => 
                 a -> VPatch -> VPatch

which filters out un-needed vertex locations from the vertex location map. Unlike subFaces, restrictTo checks for missing vertex locations, so restrictTo fcs vp raises an error if a vertex in fcs is missing from the keys of the vertex location map of vp.

5. Forcing in More Detail

The force rules

The rules used by our force algorithm are local and derived from the fact that there are seven possible vertex types as depicted in figure 8.

Our rules are shown in figure 9 (omitting mirror symmetric versions). In each case the TileFace shown yellow needs to be added in the presence of the other TileFaces shown.

Main Forcing Operations

To make forcing efficient we convert a Tgraph to a BoundaryState to keep track of boundary information of the Tgraph, and then calculate a ForceState which combines the BoundaryState with a record of awaiting boundary edge updates (an update map). Then each face addition is carried out on a ForceState, converting back when all the face additions are complete. It makes sense to apply force (and related functions) to a Tgraph, a BoundaryState, or a ForceState, so we define a class Forcible with instances Tgraph, BoundaryState, and ForceState.

This allows us to define

    force :: Forcible a => a -> a
    tryForce :: Forcible a => a -> Try a

The first will raise an error if a stuck tiling is encountered. The second uses a Try result which produces a Left string for failures and a Right a for successful result a.

There are several other operations related to forcing including

    stepForce :: Forcible a => Int -> a -> a
    tryStepForce  :: Forcible a => Int -> a -> Try a

    addHalfDart, addHalfKite :: Forcible a => Dedge -> a -> a
    tryAddHalfDart, tryAddHalfKite :: Forcible a => Dedge -> a -> Try a

The first two force (up to) a given number of steps (=face additions) and the other four add a half dart/kite on a given boundary edge.

Update Generators

An update generator is used to calculate which boundary edges can have a certain update. There is an update generator for each force rule, but also a combined (all update) generator. The force operations mentioned above all use the default all update generator (defaultAllUGen) but there are more general (with) versions that can be passed an update generator of choice. For example

    forceWith :: Forcible a => UpdateGenerator -> a -> a
    tryForceWith :: Forcible a => UpdateGenerator -> a -> Try a

In fact we defined

    force = forceWith defaultAllUGen
    tryForce = tryForceWith defaultAllUGen

We can also define

    wholeTiles :: Forcible a => a -> a
    wholeTiles = forceWith wholeTileUpdates

where wholeTileUpdates is an update generator that just finds boundary join edges to complete whole tiles.

In addition to defaultAllUGen there is also allUGenerator which does the same thing apart from how failures are reported. The reason for keeping both is that they were constructed differently and so are useful for testing.

In fact UpdateGenerators are functions that take a BoundaryState and a focus (list of boundary directed edges) to produce an update map. Each Update is calculated as either a SafeUpdate (where two of the new face edges are on the existing boundary and no new vertex is needed) or an UnsafeUpdate (where only one edge of the new face is on the boundary and a new vertex needs to be created for a new face).

    type UpdateGenerator = BoundaryState -> [Dedge] -> Try UpdateMap
    type UpdateMap = Map.Map Dedge Update
    data Update = SafeUpdate TileFace 
                | UnsafeUpdate (Vertex -> TileFace)

Completing (executing) an UnsafeUpdate requires a touching vertex check to ensure that the new vertex does not clash with an existing boundary vertex. Using an existing (touching) vertex would create a crossing boundary so such an update has to be blocked.

Forcible Class Operations

The Forcible class operations are higher order and designed to allow for easy additions of further generic operations. They take care of conversions between Tgraphs, BoundaryStates and ForceStates.

    class Forcible a where
      tryFSOpWith :: UpdateGenerator -> (ForceState -> Try ForceState) -> a -> Try a
      tryChangeBoundaryWith :: UpdateGenerator -> (BoundaryState -> Try BoundaryChange) -> a -> Try a
      tryInitFSWith :: UpdateGenerator -> a -> Try ForceState

For example, given an update generator ugen and any f:: ForceState -> Try ForceState , then f can be generalised to work on any Forcible using tryFSOpWith ugen f. This is used to define both tryForceWith and tryStepForceWith.

We also specialize tryFSOpWith to use the default update generator

    tryFSOp :: Forcible a => (ForceState -> Try ForceState) -> a -> Try a
    tryFSOp = tryFSOpWith defaultAllUGen

Similarly given an update generator ugen and any f:: BoundaryState -> Try BoundaryChange , then f can be generalised to work on any Forcible using tryChangeBoundaryWith ugen f. This is used to define tryAddHalfDart and tryAddHalfKite.

We also specialize tryChangeBoundaryWith to use the default update generator

    tryChangeBoundary :: Forcible a => (BoundaryState -> Try BoundaryChange) -> a -> Try a
    tryChangeBoundary = tryChangeBoundaryWith defaultAllUGen

Note that the type BoundaryChange contains a resulting BoundaryState, the single TileFace that has been added, a list of edges removed from the boundary (of the BoundaryState prior to the face addition), and a list of the (3 or 4) boundary edges affected around the change that require checking or re-checking for updates.

The class function tryInitFSWith will use an update generator to create an initial ForceState for any Forcible. If the Forcible is already a ForceState it will do nothing. Otherwise it will calculate updates for the whole boundary. We also have the special case

    tryInitFS :: Forcible a => a -> Try ForceState
    tryInitFS = tryInitFSWith defaultAllUGen

Efficient chains of forcing operations.

Note that (force . force) does the same as force, but we might want to chain other force related steps in a calculation.

For example, consider the following combination which, after decomposing a Tgraph, forces, then adds a half dart on a given boundary edge (d) and then forces again.

    combo :: Dedge -> Tgraph -> Tgraph
    combo d = force . addHalfDart d . force . decompose

Since decompose:: Tgraph -> Tgraph, the instances of force and addHalfDart d will have type Tgraph -> Tgraph so each of these operations, will begin and end with conversions between Tgraph and ForceState. We would do better to avoid these wasted intermediate conversions working only with ForceStates and keeping only those necessary conversions at the beginning and end of the whole sequence.

This can be done using tryFSOp. To see this, let us first re-express the forcing sequence using the Try monad, so

    force . addHalfDart d . force

becomes

    tryForce <=< tryAddHalfDart d <=< tryForce

Note that (<=<) is the Kliesli arrow which replaces composition for Monads (defined in Control.Monad). (We could also have expressed this right to left sequence with a left to right version tryForce >=> tryAddHalfDart d >=> tryForce). The definition of combo becomes

    combo :: Dedge -> Tgraph -> Tgraph
    combo d = runTry . (tryForce <=< tryAddHalfDart d <=< tryForce) . decompose

This has no performance improvement, but now we can pass the sequence to tryFSOp to remove the unnecessary conversions between steps.

    combo :: Dedge -> Tgraph -> Tgraph
    combo d = runTry . tryFSOp (tryForce <=< tryAddHalfDart d <=< tryForce) . decompose

The sequence actually has type Forcible a => a -> Try a but when passed to tryFSOp it specialises to type ForceState -> Try ForseState. This ensures the sequence works on a ForceState and any conversions are confined to the beginning and end of the sequence, avoiding unnecessary intermediate conversions.

A limitation of forcing

To avoid creating touching vertices (or crossing boundaries) a BoundaryState keeps track of locations of boundary vertices. At around 35,000 face additions in a single force operation the calculated positions of boundary vertices can become too inaccurate to prevent touching vertex problems. In such cases it is better to use

    recalibratingForce :: Forcible a => a -> a
    tryRecalibratingForce :: Forcible a => a -> Try a

These work by recalculating all vertex positions at 20,000 step intervals to get more accurate boundary vertex positions. For example, 6 decompositions of the kingGraph has 2,906 faces. Applying force to this should result in 53,574 faces but will go wrong before it reaches that. This can be fixed by calculating either

    recalibratingForce (decompositions kingGraph !!6)

or using an extra force before the decompositions

    force (decompositions (force kingGraph) !!6)

In the latter case, the final force only needs to add 17,864 faces to the 35,710 produced by decompositions (force kingGraph) !!6.

6. Advanced Operations

Guided comparison of Tgraphs

Asking if two Tgraphs are equivalent (the same apart from choice of vertex numbers) is a an np-complete problem. However, we do have an efficient guided way of comparing Tgraphs. In the module Tgraph.Rellabelling we have

    sameGraph :: (Tgraph,Dedge) -> (Tgraph,Dedge) -> Bool

The expression sameGraph (g1,d1) (g2,d2) asks if g2 can be relabelled to match g1 assuming that the directed edge d2 in g2 is identified with d1 in g1. Hence the comparison is guided by the assumption that d2 corresponds to d1.

It is implemented using

    tryRelabelToMatch :: (Tgraph,Dedge) -> (Tgraph,Dedge) -> Try Tgraph

where tryRelabelToMatch (g1,d1) (g2,d2) will either fail with a Left report if a mismatch is found when relabelling g2 to match g1 or will succeed with Right g3 where g3 is a relabelled version of g2. The successful result g3 will match g1 in a maximal tile-connected collection of faces containing the face with edge d1 and have vertices disjoint from those of g1 elsewhere. The comparison tries to grow a suitable relabelling by comparing faces one at a time starting from the face with edge d1 in g1 and the face with edge d2 in g2. (This relies on the fact that Tgraphs are connected with no crossing boundaries, and hence tile-connected.)

The above function is also used to implement

    tryFullUnion:: (Tgraph,Dedge) -> (Tgraph,Dedge) -> Try Tgraph

which tries to find the union of two Tgraphs guided by a directed edge identification. However, there is an extra complexity arising from the fact that Tgraphs might overlap in more than one tile-connected region. After calculating one overlapping region, the full union uses some geometry (calculating vertex locations) to detect further overlaps.

Finally we have

    commonFaces:: (Tgraph,Dedge) -> (Tgraph,Dedge) -> [TileFace]

which will find common regions of overlapping faces of two Tgraphs guided by a directed edge identification. The resulting common faces will be a sub-collection of faces from the first Tgraph. These are returned as a list as they may not be a connected collection of faces and therefore not necessarily a Tgraph.

Empires and SuperForce

In Empires and SuperForce we discussed forced boundary coverings which were used to implement both a superForce operation

    superForce:: Forcible a => a -> a

and operations to calculate empires.

We will not repeat the descriptions here other than to note that

    forcedBoundaryECovering:: Tgraph -> [Forced Tgraph]

finds boundary edge coverings after forcing a Tgraph. That is, forcedBoundaryECovering g will first force g, then (if it succeeds) finds a collection of (forced) extensions to force g such that

• each extension has the whole boundary of force g as internal edges.
• each possible addition to a boundary edge of force g (kite or dart) has been included in the collection.

(possible here means – not leading to a stuck Tgraph when forced.) There is also

    forcedBoundaryVCovering:: Tgraph -> [Forced Tgraph]

which does the same except that the extensions have all boundary vertices internal rather than just the boundary edges. In both cases the result is a list of explicitly forced Tgraphs (discussed next).

Combinations and Explicitly Forced

We introduced a new type Forced (in v 1.3) to enable a forcible to be explictily labelled as being forced. For example

    forceF    :: Forcible a => a -> Forced a 
    tryForceF :: Forcible a => a -> Try (Forced a)
    forgetF   :: Forced a -> a

This allows us to restrict certain functions which expect a forced argument by making this explicit.

    composeF :: Forced Tgraph -> Forced Tgraph

The definition makes use of theorems established in Graphs,Kites and Darts and Theorems that composing a forced Tgraph does not require a check (for connectedness and no crossing boundaries) and the result is also forced. This can then be used to define efficient combinations such as

    compForce:: Tgraph -> Forced Tgraph      -- compose after forcing
    compForce = composeF . forceF

    allCompForce:: Tgraph -> [Forced Tgraph] -- iterated (compose after force) while not emptyTgraph
    maxCompForce:: Tgraph -> Forced Tgraph   -- last item in allCompForce (or emptyTgraph)

Tracked Tgraphs

The type

    data TrackedTgraph = TrackedTgraph
       { tgraph  :: Tgraph
       , tracked :: [[TileFace]] 
       } deriving Show

has proven useful in experimentation as well as in producing artwork with darts and kites. The idea is to keep a record of sub-collections of faces of a Tgraph when doing both force operations and decompositions. A list of the sub-collections forms the tracked list associated with the Tgraph. We make TrackedTgraph an instance of class Forcible by having force operations only affect the Tgraph and not the tracked list. The significant idea is the implementation of

    decomposeTracked :: TrackedTgraph -> TrackedTgraph

Decomposition of a Tgraph involves introducing a new vertex for each long edge and each kite join. These are then used to construct the decomposed faces. For decomposeTracked we do the same for the Tgraph, but when it comes to the tracked collections, we decompose them re-using the same new vertex numbers calculated for the edges in the Tgraph. This keeps a consistent numbering between the Tgraph and tracked faces, so each item in the tracked list remains a sub-collection of faces in the Tgraph.

The function

    drawTrackedTgraph :: [VPatch -> Diagram B] -> TrackedTgraph -> Diagram B

is used to draw a TrackedTgraph. It uses a list of functions to draw VPatches. The first drawing function is applied to a VPatch for any untracked faces. Subsequent functions are applied to VPatches for the tracked list in order. Each diagram is beneath later ones in the list, with the diagram for the untracked faces at the bottom. The VPatches used are all restrictions of a single VPatch for the Tgraph, so will be consistent in vertex locations. When labels are used, there is also a drawTrackedTgraphRotating and drawTrackedTgraphAligning for rotating or aligning the VPatch prior to applying the drawing functions.

Note that the result of calculating empires (see Empires and SuperForce ) is represented as a TrackedTgraph. The result is actually the common faces of a forced boundary covering, but a particular element of the covering (the first one) is chosen as the background Tgraph with the common faces as a tracked sub-collection of faces. Hence we have

    empire1, empire2 :: Tgraph -> TrackedTgraph
    
    drawEmpire :: TrackedTgraph -> Diagram B

Figure 10 was also created using TrackedTgraphs.

7. Other Reading

Previous related blogs are:

• Diagrams for Penrose Tiles – the first blog introduced drawing Pieces and Patches (without using Tgraphs) and provided a version of decomposing for Patches (decompPatch).
• Graphs, Kites and Darts intoduced Tgraphs. This gave more details of implementation and results of early explorations. (The class Forcible was introduced subsequently).
• Empires and SuperForce – these new operations were based on observing properties of boundaries of forced Tgraphs.
• Graphs,Kites and Darts and Theorems established some important results relating force, compose, decompose.

In this blog we discuss some interesting properties we have discovered concerning the , , and operations along with some proofs.

Index

1. Quick Recap (including operations , , on Tgraphs)
2. Composition Problems and a Compose Force Theorem (composition is not a simple inverse to decomposition)
3. Perfect Composition Theorem (establishing relationships between , , )
4. Multiple Compositions (extending the Compose Force theorem for multiple compositions)
5. Proof of the Compose Force Theorem (showing is total on forced Tgraphs)

1. Quick Recap

Haskell diagrams allowed us to render finite patches of tiles easily as discussed in Diagrams for Penrose tiles. Following a suggestion of Stephen Huggett, we found that the description and manipulation of such tilings is greatly enhanced by using planar graphs. In Graphs, Kites and Darts we introduced a specialised planar graph representation for finite tilings of kites and darts which we called Tgraphs (tile graphs). These enabled us to implement operations that use neighbouring tile information and in particular operations , , and .

For ease of reference, we reproduce the half-tiles we are working with here.

Figure 1 shows the right-dart (RD), left-dart (LD), left-kite (LK) and right-kite (RK) half-tiles. Each has a join edge (shown dotted) and a short edge and a long edge. The origin vertex is shown red in each case. The vertex at the opposite end of the join edge from the origin we call the opp vertex, and the remaining vertex we call the wing vertex.

If the short edges have unit length then the long edges have length (the golden ratio) and all angles are multiples of (a tenth turn) with kite halves having two 2s and a 1, and dart halves having a 3 and two 1s. This geometry of the tiles is abstracted away from at the graph representation level but used when checking validity of tile additions and by the drawing functions.

There are rules for how the tiles can be put together to make a legal tiling (see e.g. Diagrams for Penrose tiles). We defined a Tgraph (in Graphs, Kites and Darts) as a list of such half-tiles which are constrained to form a legal tiling but must also be connected with no crossing boundaries (see below).

As a simple example consider kingGraph (2 kites and 3 darts round a king vertex). We represent each half-tile as a TileFace with three vertex numbers, then apply makeTgraph to the list of ten Tilefaces. The function makeTgraph :: [TileFace] -> Tgraph performs the necessary checks to ensure the result is a valid Tgraph.

kingGraph :: Tgraph
kingGraph = makeTgraph 
  [LD (1,2,3),RD (1,11,2),LD (1,4,5),RD (1,3,4),LD (1,10,11)
  ,RD (1,9,10),LK (9,1,7),RK (9,7,8),RK (5,7,1),LK (5,6,7)
  ]

To view the Tgraph we simply form a diagram (in this case 2 diagrams horizontally separated by 1 unit)

  hsep 1 [labelled drawj kingGraph, draw kingGraph]

and the result is shown in figure 2 with labels and dashed join edges (left) and without labels and join edges (right).

The boundary of the Tgraph consists of the edges of half-tiles which are not shared with another half-tile, so they go round untiled/external regions. The no crossing boundary constraint (equivalently, locally tile-connected) means that a boundary vertex has exactly two incident boundary edges and therefore has a single external angle in the tiling. This ensures we can always locally determine the relative angles of tiles at a vertex. We say a collection of half-tiles is a valid Tgraph if it constitutes a legal tiling but also satisfies the connectedness and no crossing boundaries constraints.

Our key operations on Tgraphs are , , and which are illustrated in figure 3.

Figure 3 shows the kingGraph with its decomposition above it (left), the result of forcing the kingGraph (right) and the composition of the forced kingGraph (bottom right).

Decompose

An important property of Penrose dart and kite tilings is that it is possible to divide the half-tile faces of a tiling into smaller half-tile faces, to form a new (smaller scale) tiling.

Figure 4 illustrates the decomposition of a left-dart (top row) and a left-kite (bottom row). With our Tgraph representation we simply introduce new vertices for dart and kite long edges and kite join edges and then form the new faces using these. This does not involve any geometry, because that is taken care of by drawing operations.

Force

Figure 5 illustrates the rules used by our operation (we omit a mirror-reflected version of each rule).

In each case the yellow half-tile is added in the presence of the other half-tiles shown. The yellow half-tile is forced because, by the legal tiling rules and the seven possible vertex types, there is no choice for adding a different half-tile on the edge where the yellow tile is added.

We call a Tgraph correct if it represents a tiling which can be continued infinitely to cover the whole plane without getting stuck, and incorrect otherwise. Forcing involves adding half-tiles by the illustrated rules round the boundary until either no more rules apply (in which case the result is a forced Tgraph) or a stuck tiling is encountered (in which case an incorrect Tgraph error is raised). Hence is a partial function but total on correct Tgraphs.

Compose: This is discussed in the next section.

2. Composition Problems and a Theorem

Compose Choices

For an infinite tiling, composition is a simple inverse to decomposition. However, for a finite tiling with boundary, composition is not so straight forward. Firstly, we may need to leave half-tiles out of a composition because the necessary parts of a composed half-tile are missing. For example, a half-dart with a boundary short edge or a whole kite with both short edges on the boundary must necessarily be excluded from a composition. Secondly, on the boundary, there can sometimes be a problem of choosing whether a half-dart should compose to become a half-dart or a half-kite. This choice in composing only arises when there is a half-dart with its wing on the boundary but insufficient local information to determine whether it should be part of a larger half-dart or a larger half-kite.

In the literature (see for example 1 and 2) there is an often repeated method for composing (also called inflating). This method always make the kite choice when there is a choice. Whilst this is a sound method for an unbounded tiling (where there will be no choice), we show that this is an unsound method for finite tilings as follows.

Clearly composing should preserve correctness. However, figure 6 (left) shows a correct Tgraph which is a forced queen, but the kite-favouring composition of the forced queen produces the incorrect Tgraph shown in figure 6 (centre). Applying our function to this reveals a stuck tiling and reports an incorrect Tgraph.

Our algorithm (discussed in Graphs, Kites and Darts) detects dart wings on the boundary where there is a choice and classifies them as unknowns. Our composition refrains from making a choice by not composing a half dart with an unknown wing vertex. The rightmost Tgraph in figure 6 shows the result of our composition of the forced queen with the half-tile faces left out of the composition (the remainder faces) shown in green. This avoidance of making a choice (when there is a choice) guarantees our composition preserves correctness.

Compose is a Partial Function

A different composition problem can arise when we consider Tgraphs that are not decompositions of Tgraphs. In general, is a partial function on Tgraphs.

Figure 7 shows a Tgraph (left) with its sucessful composition (centre) and the half-tile faces that would result from a second composition (right) which do not form a valid Tgraph because of a crossing boundary (at vertex 6). Thus composition of a Tgraph may fail to produce a Tgraph when the resulting faces are disconnected or have a crossing boundary.

However, we claim that is a total function on forced Tgraphs.

Compose Force Theorem

Theorem: Composition of a forced Tgraph produces a valid Tgraph.

We postpone the proof (outline) for this theorem to section 5. Meanwhile we use the result to establish relationships between , , and in the next section.

3. Perfect Composition Theorem

In Graphs, Kites and Darts we produced a diagram showing relationships between multiple decompositions of a dart and the forced versions of these Tgraphs. We reproduce this here along with a similar diagram for multiple decompositions of a kite.

In figure 8 we show separate (apparently) commuting diagrams for the dart and for the kite. The bottom rows show the decompositions, the middle rows show the result of forcing the decompositions, and the top rows illustrate how the compositions of the forced Tgraphs work by showing both the composed faces (black edges) and the remainder faces (green edges) which are removed in the composition. The diagrams are examples of some commutativity relationships concerning , and which we will prove.

It should be noted that these diagrams break down if we consider only half-tiles as the starting points (bottom right of each diagram). The decomposition of a half-tile does not recompose to its original, but produces an empty composition. So we do not even have in these cases. Forcing the decomposition also results in an empty composition. Clearly there is something special about the depicted cases and it is not merely that they are wholetile complete because the decompositions are not wholetile complete. [Wholetile complete means there are no join edges on the boundary, so every half-tile has its other half.]

Below we have captured the properties that are sufficient for the diagrams to commute as in figure 8. In the proofs we use a partial ordering on Tgraphs (modulo vertex relabelling) which we define next.

Partial ordering of Tgraphs

If and are both valid Tgraphs and consists of a subset of the (half-tile) faces of we have

which gives us a partial order on Tgraphs. Often, though, is only isomorphic to a subset of the faces of , requiring a vertex relabelling to become a subset. In that case we write

which is also a partial ordering and induces an equivalence of Tgraphs defined by

in which case and are isomorphic as Tgraphs.

Both and are monotonic with respect to meaning:

We also have is monotonic, but only when restricted to correct Tgraphs. Also, when restricted to correct Tgraphs, we have is non decreasing because it only adds faces:

and is idempotent (forcing a forced correct Tgraph leaves it the same):

Composing perfectly and perfect compositions

Definition: A Tgraph composes perfectly if all faces of are composable (i.e there are no remainder faces of when composing).

We note that the composed faces must be a valid Tgraph (connected with no crossing boundaries) if all faces are included in the composition because has those properties. Clearly, if composes perfectly then

In general, for arbitrary where the composition is defined, we only have

Definition: A Tgraph is a perfect composition if composes perfectly.

Clearly if is a perfect composition then

(We could use equality here because any new vertex labels introduced by will be removed by ). In general, for arbitrary ,

Lemma 1: is a perfect composition if and only if has the following 2 properties:

1. every half-kite with a boundary join has either a half-dart or a whole kite on the short edge, and
2. every half-dart with a boundary join has a half-kite on the short edge,

(Proof outline:) Firstly note that unknowns in (= ) can only come from boundary joins in . The properties 1 and 2 guarantee that has no unknowns. Since every face of has come from a decomposed face in , there can be no faces in that will not recompose, so will compose perfectly to . Conversely, if is a perfect composition, its decomposition can have no unknowns. This implies boundary joins in must satisfy properties 1 and 2.

(Note: a perfect composition may have unknowns even though its decomposition has none.)

It is easy to see two special cases:

1. If is wholetile complete then is a perfect composition.

   Proof: Wholetile complete implies no boundary joins which implies properties 1 and 2 in lemma 1 which implies is a perfect composition.

2. If is a decomposition then is a perfect composition.

   Proof: If is a decomposition, then every half-dart has a half-kite on the short edge which implies property 2 of lemma 1. Also, any half-kite with a boundary join in must have come from a decomposed half-dart since a decomposed half-kite produces a whole kite with no boundary kite join. So the half-kite must have a half-dart on the short edge which implies property 1 of lemma 1. The two properties imply is a perfect composition.

We note that these two special cases cover all the Tgraphs in the bottom rows of the diagrams in figure 8. So the Tgraphs in each bottom row are perfect compositions, and furthermore, they all compose perfectly except for the rightmost Tgraphs which have empty compositions.

In the following results we make the assumption that a Tgraph is correct, which guarantees that when is applied, it terminates with a correct Tgraph. We also note that preserves correctness as does (provided the composition is defined).

Lemma 2: If is a forced, correct Tgraph then

(Proof outline:) The proof uses a case analysis of boundary and internal vertices of . For internal vertices we just check there is no change at the vertex after using figure 12 (plus an extra case for the forced star). For boundary vertices we check the local context cases shown in figure 9.

This shows two cases for a kite origin, two cases for a kite opp, four cases for a kite wing, and four cases for a dart origin. The only case for a dart wing is one of the two kite opp cases and there are no dart opp cases as these cannot be on the boundary of a forced Tgraph. Actually figure 9 has a repeated case which is both a dart origin and a kite wing, but there are also 3 additional cases when we consider mirror images of those shown. Since there is no local change of the context in each case, and since this is true for all boundary vertices in any forced Tgraph, there can be no non-local change either. (We omit the full details).

Lemma 3: If is a perfect composition and a correct Tgraph, then

(Proof outline:) The proof is by analysis of each possible force rule applicable on a boundary edge of and checking local contexts to establish that (i) the result of applying to the local context must include the added half-tile, and (ii) if the added half tile has a new boundary join, then the result must include both halves of the new half-tile. The two properties of perfect compositions mentioned in lemma 1 are critical for the proof. However, since the result of adding a single half-tile may break the condition of the Tgraph being a pefect composition, we need to arrange that half-tiles are completed first then each subsequent half-tile addition is paired with its wholetile completion. This ensures the perfect composition condition holds at each step for a proof by induction. [A separate proof is needed to show that the ordering of applying force rules makes no difference to a final correct Tgraph (apart from vertex relabelling)].

Lemma 4 If composes perfectly and is a correct Tgraph then

Proof: Assume composes perfectly and is a correct Tgraph. Since is non-decreasing (with respect to on correct Tgraphs)

and since is monotonic

Since composes perfectly, the left hand side is just , so

and since is monotonic (with respect to on correct Tgraphs)

For the opposite direction, we substitute for in lemma 3 to get

Then, since , we have

Apply to both sides (using monotonicity)

For any for which the composition is defined we have so we get

Now apply to both sides and note to get

Combining this with (*) above proves the required equivalence.

Theorem (Perfect Composition): If composes perfectly and is a correct Tgraph then

Proof: Assume composes perfectly and is a correct Tgraph. By lemma 4 we have

Applying to both sides, gives

Now by lemma 2, with , the right hand side is equivalent to

which establishes the result.

Corollaries (of the perfect composition theorem):

1. If is a perfect composition and a correct Tgraph then

   

   Proof: Let (so ) in the theorem.

   [This result generalises lemma 2 because any correct forced Tgraph is necessarily wholetile complete and therefore a perfect composition, and .]

2. If is a perfect composition and a correct Tgraph then

   

   Proof: Apply to both sides of the previous corollary and note that

   

   provided the composition is defined, which it must be for a forced Tgraph by the Compose Force theorem.

3. If is a perfect composition and a correct Tgraph then

   

   Proof: Apply to both sides of the previous corollary noting is monotonic and idempotent for correct Tgraphs

   

   From the fact that is non decreasing and and are monotonic, we also have

   

   Hence combining these two sub-Tgraph results we have

   

   

It is important to point out that if is a correct Tgraph and is a perfect composition then this is not the same as composes perfectly. It could be the case that has more faces than and so could have unknowns. In this case we can only prove that

As an example where this is not an equivalence, choose to be a star. Then its composition is the empty Tgraph (which is still a pefect composition) and so the left hand side is the empty Tgraph, but the right hand side is a sun.

Perfectly composing generators

The perfect composition theorem and lemmas and the three corollaries justify all the commuting implied by the diagrams in figure 8. However, one might ask more general questions like: Under what circumstances do we have (for a correct forced Tgraph )

Definition A generator of a correct forced Tgraph is any Tgraph such that and .

We can now state that

Corollary If a correct forced Tgraph has a generator which composes perfectly, then

Proof: This follows directly from lemma 4 and the perfect composition theorem.

As an example where the required generator does not exist, consider the rightmost Tgraph of the middle row in figure 10. It is generated by the Tgraph directly below it, but it has no generator with a perfect composition. The Tgraph directly above it in the top row is the result of applying which has lost the leftmost dart of the Tgraph.

We could summarise this section by saying that can lose information which cannot be recovered by a subsequent and, similarly, can lose information which cannot be recovered by a subsequent . We have defined perfect compositions which are the Tgraphs that do not lose information when decomposed and Tgraphs which compose perfectly which are those that do not lose information when composed. Forcing does the same thing at each level of composition (that is it commutes with composition) provided information is not lost when composing.

4. Multiple Compositions

We know from the Compose Force theorem that the composition of a Tgraph that is forced is always a valid Tgraph. In this section we use this and the results from the last section to show that composing a forced, correct Tgraph produces a forced Tgraph.

First we note that:

Lemma 5: The composition of a forced, correct Tgraph is wholetile complete.

Proof: Let where is a forced, correct Tgraph. A boundary join in implies there must be a boundary dart wing of the composable faces of . (See for example figure 4 where this would be vertex 2 for the half dart case, and vertex 5 for the half-kite face). This dart wing cannot be an unknown as the half-dart is in the composable faces. However, a known dart wing must be either a large kite centre or a large dart base and therefore internal in the composable faces of (because of the force rules) and therefore not on the boundary in . This is a contradiction showing that can have no boundary joins and is therefore wholetile complete.

Theorem: The composition of a forced, correct Tgraph is a forced Tgraph.

Proof: Let for some forced, correct Tgraph , then is wholetile complete (by lemma 5) and therefore a perfect composition. Let , so composes perfectly (). By the perfect composition theorem we have

We also have

Applying to both sides, noting that is monotonic and the identity on forced Tgraphs, we have

Applying to both sides, noting that is monotonic, we have

By (**) above, the left hand side is equivalent to so we have

but since we also have ( being non-decreasing)

we have established that

which means is a forced Tgraph.

This result means that after forcing once we can repeatedly compose creating valid Tgraphs until we reach the empty Tgraph.

We can also use lemma 5 to establish the converse to a previous corollary:

Corollary If a correct forced Tgraph satisfies:

then has a generator which composes perfectly.

Proof: By lemma 5, is wholetile complete and hence a perfect composition. This means that composes perfectly and it is also a generator for because

5. Proof of the Compose Force theorem

Theorem (Compose Force): Composition of a forced Tgraph produces a valid Tgraph.

Proof: For any forced Tgraph we can construct the composed faces. For the result to be a valid Tgraph we need to show no crossing boundaries and connectedness for the composed faces. These are proved separately by case analysis below.

Proof of no crossing boundaries

Assume is a forced Tgraph and that it has a non-empty set of composed faces (we can ignore cases where the composition is empty as the empty Tgraph is valid). Consider a vertex v in the composed faces of and first take the case that v is on the boundary of . We consider local contexts for a vertex v on a forced Tgraph boundary where the composition is non-empty as shown in figure 11.

In each case v is shown as a red dot, and the composition is shown filled yellow. The cases for v are shown in rows: the first row is for dart origins, the second row is for kite origins, the third row is for kite wings, and the last rows is for kite opps. The dart wing cases are a subset of the kite opp cases, so not repeated, and dart opp vertices are excluded because they cannot be on the boundary of a forced Tgraph. We only show left-hand versions, so there is a mirror symmetric set for right-hand versions.

It is easy to see that there are no crossing boundaries of the composed faces at v in each case. Since any boundary vertex of any forced Tgraph (with a non-empty composition) must match one of these local context cases around the vertex, we can conclude that a boundary vertex of cannot become a crossing boundary in .

Next take the case where v is an internal vertex of .

Figure 12 shows relationships between the forced Tgraphs of the 7 (internal) vertex types (plus a kite at the top right). The red faces are those around the vertex type and the black faces are those produced by forcing (if any). Each forced Tgraph has its composition directly above with empty compositions for the top row. We note that a (forced) star, jack, king, and queen vertex remains an internal vertex in the respective composition so cannot become a crossing boundary vertex. A deuce vertex becomes the centre of a larger kite and is no longer present in the composition (top right). That leaves cases for the sun vertex and ace vertex (=fool vertex). The sun Tgraph (sunGraph) and fool Tgraph (fool) consist of just the red faces at the respective vertex (shown top left and top centre). These both have empty compositions when there is no surrounding context. We thus need to check possible forced local contexts for sunGraph and fool.

The fool case is simple and similar to a duece vertex in that it is never part of a composition. [To see this consider inverting the decomposition arrows shown in figure 4. In both cases we see the half-dart opp vertex (labelled 4 in figure 4) is removed].

For the sunGraph there are only 7 local forced context cases to consider where the sun vertex is on the boundary of the composition.

Six of these are shown in figure 13 (the missing one is just a mirror reflection of the fourth case). Again, the relevant vertex v is shown as a red dot and the composed faces are shown filled yellow, so it is easy to check that there is no crossing boundary of the composed faces at v in each case. Every forced Tgraph containing an internal sun vertex where the vertex is on the boundary of the composition must match one of the 7 cases locally round the vertex.

Thus no vertex from can become a crossing boundary vertex in the composed faces and since the vertices of the composed faces are a subset of those of , we can have no crossing boundary vertex in the composed faces.

Proof of Connectedness

Assume is a forced Tgraph as before. We refer to the half-tile faces of that get included in the composed faces as the composable faces and the rest as the remainder faces. We want to prove that the composable faces are connected as this will imply the composed faces are connected.

As before we can ignore cases where the set of composable faces is empty, and assume this is not the case. We study the nature of the remainder faces of . Firstly, we note:

Lemma (remainder faces)

The remainder faces of are made up entirely of groups of half-tiles which are either:

1. Half-fools (= a half dart and both halves of the kite attached to its short edge) where the other half-fool is entirely composable faces, or
2. Both halves of a kite with both short edges on the () boundary (so they are not part of a half-fool) where (at most) the origin is in common with composable faces, or
3. Whole fools where (at most) the shared kite origin in common with composable faces.

These 3 cases of remainder face groups are shown in figure 14. In each case the possible border in common with composable faces is shown yellow and the red edges are necessarily on the boundary of (the black boundary could be on the boundary of or shared with another reamainder face group). [A mirror symmetric version for the first group is not shown.] Examples can be seen in e.g. figure 13 where the first Tgraph has four examples of case 1, and two of case 2, the second has six examples of case 1 and two of case 2, and the fifth Tgraph has an example of case 3 as well as four of case 1. [We omit the detailed proof of this lemma which reasons about what gets excluded in a composition after forcing. However, all the local context cases are included in figure 15 (left-hand versions), where we only show those contexts where there is a non-empty composition.]

We note from the (remainder faces) lemma that the common boundary of the group of remainder faces with the composable faces (shown yellow in figure 14) is at most a single vertex in cases 2 and 3. In case 1, the common boundary is just a single edge of the composed faces which is made up of 2 adjacent edges of the composable faces that constitute the join of two half-fools.

This means each (remainder face) group shares boundary with exactly one connected component of the composable faces.

Next we establish that if two (remainder face) groups are connected they must share boundary with the same connected component of the composable faces. We need to consider how each (remainder face) group can be connected with a neighbouring such group. It is enough to consider forced contexts of boundary dart long edges (for cases 1 and 3) and boundary kite short edges (for case 2). The cases where the composition is non-empty all appear in figure 15 (left-hand versions) along with boundary kite long edges (middle two rows) which are not relevant for the proof.

We note that, whenever one group of the remainder faces (half-fool, whole-kite, whole-fool) is connected to a neighbouring group of the remainder faces, any common boundary (shared edges and vertices) with the compososable faces is also connected. The combined common boundary forms either 2 adjacent composed face boundary edges (= 4 adjacent edges of the composable faces), or a composed face boundary edge and one of its end vertices, or a single composed face boundary vertex.

It follows that any connected collection of the remainder face groups shares boundary with a unique connected component of the composable faces. Since the collection of composable and remainder faces together is connected ( is connected) the removal of the remainder faces cannot disconnect the composable faces. For this to happen, at least one connected collection of remainder face groups would have to be connected to more than one connected component of composable faces.

This establishes connectedness of any composition of a forced Tgraph, and this completes the proof of the Compose Force theorem.

References

[1] Martin Gardner (1977) MATHEMATICAL GAMES. Scientific American, 236(1), (pages 110 to 121). http://www.jstor.org/stable/24953856

[2] Grünbaum B., Shephard G.C. (1987) Tilings and Patterns. W. H. Freeman and Company, New York. ISBN 0-7167-1193-1 (Hardback) (pages 540 to 542).

Previously in Diagrams for Penrose tiles we implemented tools to draw finite tilings using Haskell diagrams. There we also noted that legal tilings are only correct tilings if they can be continued infinitely and are incorrect otherwise. In Graphs, Kites and Darts we introduced a graph representation for finite tilings (Tgraphs) which enabled us to implement operations that use neighbouring tile information. In particular we implemented a force operation to extend a Tgraph on any boundary edge where there is a unique choice for adding a tile.

In this note we find a limitation of force, show a way to improve on it (superForce), and introduce boundary coverings which are used to implement superForce and calculate empires.

Properties of Tgraphs

A Tgraph is a collection of half-tile faces representing a legal tiling and a half-tile face is either an LD (left dart) , RD (right dart), LK (left kite), or RK (right kite) each with 3 vertices to form a triangle. Faces of the Tgraph which are not half-tile faces are considered external regions and those edges round the external regions are the boundary edges of the Tgraph. The half-tile faces in a Tgraph are required to be connected and locally tile-connected which means that there are exactly two boundary edges at any boundary vertex (no crossing boundaries).

As an example Tgraph we show kingGraph (the three darts and two kites round a king vertex), where

  kingGraph = makeTgraph 
    [LD (1,2,3),RD (1,11,2),LD (1,4,5),RD (1,3,4),LD (1,10,11)
    ,RD (1,9,10),LK (9,1,7),RK (9,7,8),RK (5,7,1),LK (5,6,7)
    ]

This is drawn in figure 1 using

  hsep 1 [labelled drawj kingGraph, draw kingGraph]

which shows vertex labels and dashed join edges (left) and without labels and join edges (right). (hsep 1 provides a horizontal seperator of unit length.)

Properties of forcing

We know there are at most two legal possibilities for adding a half-tile on a boundary edge of a Tgraph. If there are zero legal possibilities for adding a half-tile to some boundary edge, we have a stuck tiling/incorrect Tgraph.

Forcing deals with all cases where there is exactly one possibility for extending on a boundary edge according to the legal tiling rules and consistent with the seven possible vertex types. That means forcing either fails at some stage with a stuck Tgraph (indicating the starting Tgraph was incorrect) or it enlarges the starting Tgraph until every boundary edge has exactly two legal possibilities (consistent with the seven vertex types) for adding a half-tile so a choice would need to be made to grow the Tgraph any further.

Figure 2 shows force kingGraph with kingGraph shown red.

If g is a correct Tgraph, then force g succeeds and the resulting Tgraph will be common to all infinite tilings that extend the finite tiling represented by g. However, we will see that force g is not a greatest lower bound of (infinite) tilings that extend g. Firstly, what is common to all extensions of g may not be a connected collection of tiles. This leads to the concept of empires which we discuss later. Secondly, even if we only consider the connected common region containing g, we will see that we need to go beyond force g to find this, leading to an operation we call superForce.

Our empire and superForce operations are implemented using boundary coverings which we introduce next.

Boundary edge covering

Given a successfully forced Tgraph fg, a boundary edge covering of fg is a list of successfully forced extensions of fg such that

1. no boundary edge of fg remains on the boundary in each extension, and
2. the list takes into account all legal choices for extending on each boundary edge of fg.

[Technically this is a covering of the choices round the boundary, but each extension is also a cover of the boundary edges.] Figure 3 shows a boundary edge covering for a forced kingGraph (force kingGraph is shown red in each extension).

In practice, we do not need to explore both choices for every boundary edge of fg. When one choice is made, it may force choices for other boundary edges, reducing the number of boundary edges we need to consider further.

The main function is boundaryECovering working on a BoundaryState (which is a Tgraph with extra boundary information). It uses covers which works on a list of extensions each paired with the remaining set of the original boundary edges not yet covered. (Initially covers is given a singleton list with the starting boundary state and the full set of boundary edges to be covered.) For each extension in the list, if its uncovered set is empty, that extension is a completed cover. Otherwise covers replaces the extension with further extensions. It picks the (lowest numbered) boundary edge in the uncovered set, tries extending with a half-dart and with a half-kite on that edge, forcing in each case, then pairs each result with its set of remaining uncovered boundary edges before adding the resulting extensions back at the front of the list to be processed again. If one of the choices for a dart/kite leads to an incorrect tiling (a stuck tiling) when forced, that choice is dropped (provided the other choice succeeds). The final list returned consists of all the completed covers.

  boundaryECovering:: BoundaryState -> [BoundaryState]
  boundaryECovering bs = covers [(bs, Set.fromList (boundary bs))]

  covers:: [(BoundaryState, Set.Set Dedge)] -> [BoundaryState]
  covers [] = []
  covers ((bs,es):opens) 
    | Set.null es = bs:covers opens -- bs is complete
    | otherwise   = covers (newcases ++ opens)
       where (de,des) = Set.deleteFindMin es
             newcases = fmap (\b -> (b, commonBdry des b))
                             (atLeastOne $ tryDartAndKite bs de)

Here we have used

  type Try a = Either String a
  tryDartAndKite:: BoundaryState -> Dedge -> [Try BoundaryState]
  atLeastOne    :: [Try a] -> [a]

We frequently use Try as a type for results of partial functions where we need to continue computation if there is a failure. For example we have a version of force (called tryForce) that returns a Try Tgraph so it does not fail by raising an error, but returns a result indicating either an explicit failure situation or a successful result with a final forced Tgraph. The function tryDartAndKite tries adding an appropriate half-dart and half-kite on a given boundary edge, then uses tryForceBoundary (a variant of tryForce which works with boundary states) on each result and returns a list of Try results. The list of Try results is converted with atLeastOne which collects the successful results but will raise an error when there are no successful results.

Boundary vertex covering

You may notice in figure 3 that the top right cover still has boundary vertices of kingGraph on the final boundary. We use a boundary vertex covering rather than a boundary edge covering if we want to exclude these cases. This involves picking a boundary edge that includes such a vertex and continuing the process of growing possible extensions until no boundary vertices of the original remain on the boundary.

Empires

A partial example of an empire was shown in a 1977 article by Martin Gardner 1. The full empire of a finite tiling would consist of the common faces of all the infinite extensions of the tiling. This will include at least the force of the tiling but it is not obviously finite. Here we confine ourselves to the empire in finite local regions.

For example, we can calculate a local empire for a given Tgraph g by finding the common faces of all the extensions in a boundary vertex covering of force g (which we call empire1 g).

This requires an efficient way to compare Tgraphs. We have implemented guided intersection and guided union operations which, when given a common edge starting point for two Tgraphs, proceed to compare the Tgraphs face by face and produce an appropriate relabelling of the second Tgraph to match the first Tgraph only in the overlap where they agree. These operations may also use geometric positioning information to deal with cases where the overlap is not just a single connected region. From these we can return a union as a single Tgraph when it exists, and an intersection as a list of common faces. Since the (guided) intersection of Tgraphs (the common faces) may not be connected, we do not have a resulting Tgraph. However we can arbitrarily pick one of the argument Tgraphs and emphasise which are the common faces in this example Tgraph.

Figure 4 (left) shows empire1 kingGraph where the starting kingGraph is shown in red. The grey-filled faces are the common faces from a boundary vertex covering. We can see that these are not all connected and that the force kingGraph from figure 2 corresponds to the connected set of grey-filled faces around and including the kingGraph in figure 4.

We call this a level 1 empire because we only explored out as far as the first boundary covering. We could instead, find further boundary coverings for each of the extensions in a boundary covering. This grows larger extensions in which to find common faces. On the right of figure 4 is a level 2 empire (empire2 kingGraph) which finds the intersection of the combined boundary edge coverings of each extension in a boundary edge covering of force kingGraph. Obviously this process could be continued further but, in practice, it is too inefficient to go much further.

SuperForce

We might hope that (when not discovering an incorrect tiling), force g produces the maximal connected component containing g of the common faces of all infinite extensions of g. This is true for the kingGraph as noted in figure 4. However, this is not the case in general.

The problem is that forcing will not discover if one of the two legal choices for extending a resulting boundary edge always leads to an incorrect Tgraph. In such a situation, the other choice would be common to all infinite extensions.

We can use a boundary edge covering to reveal such cases, leading us to a superForce operation. For example, figure 5 shows a boundary edge covering for the forced Tgraph shown in red.

This example is particularly interesting because in every case, the leftmost end of the red forced Tgraph has a dart immediately extending it. Why is there no case extending one of the leftmost two red edges with a half-kite? The fact that such cases are missing from the boundary edge covering suggests they are not possible. Indeed we can check this by adding a half-kite to one of the edges and trying to force. This leads to a failure showing that we have an incorrect tiling. Figure 6 illustrates the Tgraph at the point that it is discovered to be stuck (at the bottom left) by forcing.

Our superForce operation starts by forcing a Tgraph. After a successful force, it creates a boundary edge covering for the forced Tgraph and checks to see if there is any boundary edge of the forced Tgraph for which each cover has the same choice. If so, that choice is made to extend the forced Tgraph and the process is repeated by applying superForce to the result. Otherwise, just the result of forcing is returned.

Figure 7 shows a chain of examples (rockets) where superForce has been used. In each case, the starting Tgraph is shown red, the additional faces added by forcing are shown black, and any further extension produced by superForce is shown in blue.

Coda

We still do not know if forcing decides that a Tgraph is correct/incorrect. Can we conclude that if force g succeeds then g (and force g) are correct? We found examples (rockets in figure 7) where force succeeds but one of the 2 legal choices for extending on a boundary edge leads to an incorrect Tgraph. If we find an example g where force g succeeds but both legal choices on a boundary edge lead to incorrect Tgraphs we will have a counter-example. If such a g exists then superForce g will raise an error. [The calculation of a boundary edge covering will call atLeastOne where both branches have led to failure for extending on an edge.]

This means that when superForce succeeds every resulting boundary edge has two legal extensions, neither of which will get stuck when forced.

I would like to thank Stephen Huggett who suggested the idea of using graphs to represent tilings and who is working with me on proof problems relating to the kite and dart tilings.

Reference [1] Martin Gardner (1977) MATHEMATICAL GAMES. Scientific American, 236(1), (pages 110 to 121). http://www.jstor.org/stable/24953856

Non-periodic tilings with Penrose’s kites and darts

(An updated version, since original posting on Jan 6, 2022)

We continue our investigation of the tilings using Haskell with Haskell Diagrams. What is new is the introduction of a planar graph representation. This allows us to define more operations on finite tilings, in particular forcing and composing.

Previously in Diagrams for Penrose Tiles we implemented tools to create and draw finite patches of Penrose kites and darts (such as the samples depicted in figure 1). The code for this and for the new graph representation and tools described here can be found on GitHub https://github.com/chrisreade/PenroseKiteDart.

To describe the tiling operations it is convenient to work with the half-tiles: LD (left dart), RD (right dart), LK (left kite), RK (right kite) using a polymorphic type HalfTile (defined in a module HalfTile)

data HalfTile rep 
 = LD rep | RD rep | LK rep | RK rep   deriving (Show,Eq)

Here rep is a type variable for a representation to be chosen. For drawing purposes, we chose two-dimensional vectors (V2 Double) and called these Pieces.

type Piece = HalfTile (V2 Double)

The vector represents the join edge of the half tile (see figure 2) and thus the scale and orientation are determined (the other tile edges are derived from this when producing a diagram).

Finite tilings or patches are then lists of located pieces.

type Patch = [Located Piece]

Both Piece and Patch are made transformable so rotate, and scale can be applied to both and translate can be applied to a Patch. (Translate has no effect on a Piece unless it is located.)

In Diagrams for Penrose Tiles we also discussed the rules for legal tilings and specifically the problem of incorrect tilings which are legal but get stuck so cannot continue to infinity. In order to create correct tilings we implemented the decompose operation on patches.

The vector representation that we use for drawing is not well suited to exploring properties of a patch such as neighbours of pieces. Knowing about neighbouring tiles is important for being able to reason about composition of patches (inverting a decomposition) and to find which pieces are determined (forced) on the boundary of a patch.

However, the polymorphic type HalfTile allows us to introduce our alternative graph representation alongside Pieces.

Tile Graphs

In the module Tgraph.Prelude, we have the new representation which treats half tiles as triangular faces of a planar graph – a TileFace – by specialising HalfTile with a triple of vertices (clockwise starting with the tile origin). For example

LD (1,3,4)       RK (6,4,3)

type Vertex = Int
type TileFace = HalfTile (Vertex,Vertex,Vertex)

When we need to refer to particular vertices from a TileFace we use originV (the first vertex – red dot in figure 2), oppV (the vertex at the opposite end of the join edge – dashed edge in figure 2), wingV (the remaining vertex not on the join edge).

originV, oppV, wingV :: TileFace -> Vertex

Tgraphs

The Tile Graphs implementation uses a newtype Tgraph which is a list of tile faces.

newtype Tgraph = Tgraph [TileFace]
                 deriving (Show)

faces :: Tgraph -> [TileFace]
faces (Tgraph fcs) = fcs

For example, fool (short for a fool’s kite) is a Tgraph with 6 faces (and 7 vertices), shown in figure 3.

fool = Tgraph [RD (1,2,3),LD (1,3,4),RK (6,2,5)
              ,LK (6,3,2),RK (6,4,3),LK (6,7,4)
              ]

(The fool is also called an ace in the literature)

With this representation we can investigate how composition works with whole patches. Figure 4 shows a twice decomposed sun on the left and a once decomposed sun on the right (both with vertex labels). In addition to decomposing the right Tgraph to form the left Tgraph, we can also compose the left Tgraph to get the right Tgraph.

After implementing composition, we also explore a force operation and an emplace operation to extend tilings.

There are some constraints we impose on Tgraphs.

• No spurious vertices. The vertices of a Tgraph are the vertices that occur in the faces of the Tgraph (and maxV is the largest number occurring).
• Connected. The collection of faces must be a single connected component.
• No crossing boundaries. By this we mean that vertices on the boundary are incident with exactly two boundary edges. The boundary consists of the edges between the Tgraph faces and exterior region(s). This is important for adding faces.
• Tile connected. Roughly, this means that if we collect the faces of a Tgraph by starting from any single face and then add faces which share an edge with those already collected, we get all the Tgraph faces. This is important for drawing purposes.

In fact, if a Tgraph is connected with no crossing boundaries, then it must be tile connected. (We could define tile connected to mean that the dual graph excluding exterior regions is connected.)

Figure 5 shows two excluded graphs which have crossing boundaries at 4 (left graph) and 13 (right graph). The left graph is still tile connected but the right is not tile connected (the two faces at the top right do not have an edge in common with the rest of the faces.)

Although we have allowed for Tgraphs with holes (multiple exterior regions), we note that such holes cannot be created by adding faces one at a time without creating a crossing boundary. They can be created by removing faces from a Tgraph without necessarily creating a crossing boundary.

Important We are using face as an abbreviation for half-tile face of a Tgraph here, and we do not count the exterior of a patch of faces to be a face. The exterior can also be disconnected when we have holes in a patch of faces and the holes are not counted as faces either. In graph theory, the term face would generally include these other regions, but we will call them exterior regions rather than faces.

In addition to the constructor Tgraph we also use

checkedTgraph:: [TileFace] -> Tgraph

which creates a Tgraph from a list of faces, but also performs checks on the required properties of Tgraphs. We can then remove or select faces from a Tgraph and then use checkedTgraph to ensure the resulting Tgraph still satisfies the required properties.

selectFaces, removeFaces  :: [TileFace] -> Tgraph -> Tgraph
selectFaces fcs g = checkedTgraph (faces g `intersect` fcs)
removeFaces fcs g = checkedTgraph (faces g \\ fcs)

Edges and Directed Edges

We do not explicitly record edges as part of a Tgraph, but calculate them as needed. Implicitly we are requiring

• No spurious edges. The edges of a Tgraph are the edges of the faces of the Tgraph.

To represent edges, a pair of vertices (a,b) is regarded as a directed edge from a to b. A list of such pairs will usually be regarded as a directed edge list. In the special case that the list is symmetrically closed [(b,a) is in the list whenever (a,b) is in the list] we will refer to this as an edge list rather than a directed edge list.

The following functions on TileFaces all produce directed edges (going clockwise round a face).

type Dedge = (Vertex,Vertex)

joinE  :: TileFace -> Dedge  -- join edge - dashed in figure 2
shortE :: TileFace -> Dedge  -- the short edge which is not a join edge
longE  :: TileFace -> Dedge  -- the long edge which is not a join edge
faceDedges :: TileFace -> [Dedge]
  -- all three directed edges clockwise from origin

For the whole Tgraph, we often want a list of all the directed edges of all the faces.

graphDedges :: Tgraph -> [Dedge]
graphDedges = concatMap faceDedges . faces

Because our graphs represent tilings they are planar (can be embedded in a plane) so we know that at most two faces can share an edge and they will have opposite directions of the edge. No two faces can have the same directed edge. So from graphDedges g we can easily calculate internal edges (edges shared by 2 faces) and boundary directed edges (directed edges round the external regions).

internalEdges, boundaryDedges :: Tgraph -> [Dedge]

The internal edges of g are those edges which occur in both directions in graphDedges g. The boundary directed edges of g are the missing reverse directions in graphDedges g.

We also refer to all the long edges of a Tgraph (including kite join edges) as phiEdges (both directions of these edges).

phiEdges :: Tgraph -> [Dedge]

This is so named because, when drawn, these long edges are phi times the length of the short edges (phi being the golden ratio which is approximately 1.618).

Drawing Tgraphs (Patches and VPatches)

The module Tgraph.Convert contains functions to convert a Tgraph to our previous vector representation (Patch) defined in TileLib so we can use the existing tools to produce diagrams.

However, it is convenient to have an intermediate stage (a VPatch = Vertex Patch) which contains both faces and calculated vertex locations (a finite map from vertices to locations). This allows vertex labels to be drawn and for faces to be identified and retained/excluded after the location information is calculated.

data VPatch = VPatch { vLocs :: VertexLocMap
                     , vpFaces::[TileFace]
                     } deriving Show

The conversion functions include

makeVP   :: Tgraph -> VPatch

For drawing purposes we introduced a class Drawable which has a means to create a diagram when given a function to draw Pieces.

class Drawable a where
  drawWith :: (Piece -> Diagram B) -> a -> Diagram B

This allows us to make Patch, VPatch and Tgraph instances of Drawable, and we can define special cases for the most frequently used drawing tools.

draw :: Drawable a => a -> Diagram B
draw = drawWith drawPiece

drawj :: Drawable a => a -> Diagram B
drawj = drawWith dashjPiece

We also need to be able to create diagrams with vertex labels, so we use a draw function modifier

class DrawableLabelled a where
  labelSize :: Measure Double -> (VPatch -> Diagram B) -> a -> Diagram B

Both VPatch and Tgraph are made instances (but not Patch as this no longer has vertex information). The type Measure is defined in Diagrams, but we generally use a default measure for labels to define

labelled :: DrawableLabelled a => (VPatch -> Diagram B) -> a -> Diagram B
labelled = labelSize (normalized 0.018)

This allows us to use, for example (where g is a Tgraph or VPatch)

labelled draw g
labelled drawj g

One consequence of using abstract graphs is that there is no unique predefined way to orient or scale or position the VPatch (and Patch) arising from a Tgraph representation. Our implementation selects a particular join edge and aligns it along the x-axis (unit length for a dart, philength for a kite) and tile-connectedness ensures the rest of the VPatch (and Patch) can be calculated from this.

We also have functions to re-orient a VPatch and lists of VPatchs using chosen pairs of vertices. [Simply doing rotations on the final diagrams can cause problems if these include vertex labels. We do not, in general, want to rotate the labels – so we need to orient the VPatch before converting to a diagram]

Decomposing Graphs

We previously implemented decomposition for patches which splits each half-tile into two or three smaller scale half-tiles.

decompPatch :: Patch -> Patch

We now have a Tgraph version of decomposition in the module Tgraph.Decompose:

decompose :: Tgraph -> Tgraph

Graph decomposition is particularly simple. We start by introducing one new vertex for each long edge (the phiEdges) of the Tgraph. We then build the new faces from each old face using the new vertices.

As a running example we take fool (mentioned above) and its decomposition foolD

*Main> foolD = decompose fool

*Main> foolD
Tgraph [LK (1,8,3),RD (2,3,8),RK (1,3,9)
       ,LD (4,9,3),RK (5,13,2),LK (5,10,13)
       ,RD (6,13,10),LK (3,2,13),RK (3,13,11)
       ,LD (6,11,13),RK (3,14,4),LK (3,11,14)
       ,RD (6,14,11),LK (7,4,14),RK (7,14,12)
       ,LD (6,12,14)
       ]

which are best seen together (fool followed by foolD) in figure 6.

Composing Tgraphs, and Unknowns

Composing is meant to be an inverse to decomposing, and one of the main reasons for introducing our graph representation. In the literature, decomposition and composition are defined for infinite tilings and in that context they are unique inverses to each other. For finite patches, however, we will see that composition is not always uniquely determined.

In figure 7 (Two Levels) we have emphasised the larger scale faces on top of the smaller scale faces.

How do we identify the composed tiles? We start by classifying vertices which are at the wing tips of the (smaller) darts as these determine how things compose. In the interior of a graph/patch (e.g in figure 7), a dart wing tip always coincides with a second dart wing tip, and either

1. the 2 dart halves share a long edge. The shared wing tip is then classified as a largeKiteCentre and is at the centre of a larger kite. (See left vertex type in figure 8), or
2. the 2 dart halves touch at their wing tips without sharing an edge. This shared wing tip is classified as a largeDartBase and is the base of a larger dart. (See right vertex type in figure 8)

[We also call these (respectively) a deuce vertex type and a jack vertex type later in figure 10]

Around the boundary of a Tgraph, the dart wing tips may not share with a second dart. Sometimes the wing tip has to be classified as unknown but often it can be decided by looking at neighbouring tiles. In this example of a four times decomposed sun (sunD4), it is possible to classify all the dart wing tips as a largeKiteCentre or a largeDartBase so there are no unknowns.

If there are no unknowns, then we have a function to produce the unique composed Tgraph.

compose:: Tgraph -> Tgraph

Any correct decomposed Tgraph without unknowns will necessarily compose back to its original. This makes compose a left inverse to decompose provided there are no unknowns.

For example, with an (n times) decomposed sun we will have no unknowns, so these will all compose back up to a sun after n applications of compose. For n=4 (sunD4 – the smaller scale shown in figure 7) the dart wing classification returns 70 largeKiteCentres, 45 largeDartBases, and no unknowns.

Similarly with the simpler foolD example, if we classsify the dart wings we get

largeKiteCentres = [14,13]
largeDartBases = [3]
unknowns = []

In foolD (the right hand Tgraph in figure 6), nodes 14 and 13 are new kite centres and node 3 is a new dart base. There are no unknowns so we can use compose safely

*Main> compose foolD
Tgraph [RD (1,2,3),LD (1,3,4),RK (6,2,5)
       ,RK (6,4,3),LK (6,3,2),LK (6,7,4)
       ]

which reproduces the original fool (left hand Tgraph in figure 6).

However, if we now check out unknowns for fool we get

largeKiteCentres = []
largeDartBases = []
unknowns = [4,2]    

So both nodes 2 and 4 are unknowns. It had looked as though fool would simply compose into two half kites back-to-back (sharing their long edge not their join), but the unknowns show there are other possible choices. Each unknown could become a largeKiteCentre or a largeDartBase.

The question is then what to do with unknowns.

Partial Compositions

In fact our compose resolves two problems when dealing with finite patches. One is the unknowns and the other is critical missing faces needed to make up a new face (e.g the absence of any half dart).

It is implemented using an intermediary function for partial composition

partCompose:: Tgraph -> ([TileFace],Tgraph) 

partCompose will compose everything that is uniquely determined, but will leave out faces round the boundary which cannot be determined or cannot be included in a new face. It returns the faces of the argument Tgraph that were not used, along with the composed Tgraph.

Figure 9 shows the result of partCompose applied to two graphs. [These are force kiteD3 and force dartD3 on the left. Force is described later]. In each case, the excluded faces of the starting Tgraph are shown in pale green, overlaid by the composed Tgraph on the right.

Then compose is simply defined to keep the composed faces and ignore the unused faces produced by partCompose.

compose:: Tgraph -> Tgraph
compose = snd . partCompose 

This approach avoids making a decision about unknowns when composing, but it may lose some information by throwing away the uncomposed faces.

For correct Tgraphs g, if decompose g has no unknowns, then compose is a left inverse to decompose. However, if we take g to be two kite halves sharing their long edge (not their join edge), then these decompose to fool which produces an empty Tgraph when recomposed. Thus we do not have g = compose (decompose g) in general. On the other hand we do have g = compose (decompose g) for correct whole-tile Tgraphs g (whole-tile means all half-tiles of g have their matching half-tile on their join edge in g)

Later (figure 21) we show another exception to g = compose (decompose g) with an incorrect tiling.

We make use of

selectFacesVP    :: [TileFace] -> VPatch -> VPatch
removeFacesVP    :: [TileFace] -> VPatch -> VPatch

for creating VPatches from selected tile faces of a Tgraph or VPatch. This allows us to represent and draw a list of faces which need not be connected nor satisfy the no crossing boundaries property provided the Tgraph it was derived from had these properties.

Forcing

When building up a tiling, following the rules, there is often no choice about what tile can be added alongside certain tile edges at the boundary. Such additions are forced by the existing patch of tiles and the rules. For example, if a half tile has its join edge on the boundary, the unique mirror half tile is the only possibility for adding a face to that edge. Similarly, the short edge of a left (respectively, right) dart can only be matched with the short edge of a right (respectively, left) kite. We also make use of the fact that only 7 types of vertex can appear in (the interior of) a patch, so on a boundary vertex we sometimes have enough of the faces to determine the vertex type. These are given the following names in the literature (shown in figure 10): sun, star, jack (=largeDartBase), queen, king, ace (=fool), deuce (=largeKiteCentre).

The function

force :: Tgraph -> Tgraph

will add some faces on the boundary that are forced (i.e new faces where there is exactly one possible choice). For example:

• When a join edge is on the boundary – add the missing half tile to make a whole tile.
• When a half dart has its short edge on the boundary – add the half kite that must be on the short edge.
• When a vertex is both a dart origin and a kite wing (it must be a queen or king vertex) – if there is a boundary short edge of a kite half at the vertex, add another kite half sharing the short edge, (this converts 1 kite to 2 and 3 kites to 4 in combination with the first rule).
• When two half kites share a short edge their common oppV vertex must be a deuce vertex – add any missing half darts needed to complete the vertex.
• …

Figure 11 shows foolDminus (which is foolD with 3 faces removed) on the left and the result of forcing, ie force foolDminus on the right which is the same Tgraph we get from force foolD (modulo vertex renumbering).

foolDminus = 
    removeFaces [RD(6,14,11), LD(6,12,14), RK(5,13,2)] foolD

Figures 12, 13 and 14 illustrate the result of forcing a 5-times decomposed kite, a 5-times decomposed dart, and a 5-times decomposed sun (respectively). The first two figures reproduce diagrams from an article by Roger Penrose illustrating the extent of influence of tiles round a decomposed kite and dart. [Penrose R Tilings and quasi-crystals; a non-local growth problem? in Aperiodicity and Order 2, edited by Jarich M, Academic Press, 1989. (fig 14)].

In figure 15, the bottom row shows successive decompositions of a dart (dashed blue arrows from right to left), so applying compose to each dart will go back (green arrows from left to right). The black vertical arrows are force. The solid blue arrows from right to left are (force . decompose) being applied to the successive forced Tgraphs. The green arrows in the reverse direction are compose again and the intermediate (partCompose) figures are shown in the top row with the remainder faces in pale green.

Figure 16 shows the forced graphs of the seven vertex types (with the starting Tgraphs in red) along with a kite (top right).

These are related to each other as shown in the columns. Each Tgraph composes to the one above (an empty Tgraph for the ones in the top row) and the Tgraph below is its forced decomposition. [The rows have been scaled differently to make the vertex types easier to see.]

Adding Faces to a Tgraph

This is technically tricky because we need to discover what vertices (and implicitly edges) need to be newly created and which ones already exist in the Tgraph. This goes beyond a simple graph operation and requires use of the geometry of the faces. We have chosen not to do a full conversion to vectors to work out all the geometry, but instead we introduce a local representation of relative directions of edges at a vertex allowing a simple equality test.

Edge directions

All directions are integer multiples of 1/10th turn (mod 10) so we use these integers for face internal angles and boundary external angles. The face adding process always adds to the right of a given directed edge (a,b) which must be a boundary directed edge. [Adding to the left of an edge (a,b) would mean that (b,a) will be the boundary direction and so we are really adding to the right of (b,a)]. Face adding looks to see if either of the two other edges already exist in the Tgraph by considering the end points a and b to which the new face is to be added, and checking angles.

This allows an edge in a particular sought direction to be discovered. If it is not found it is assumed not to exist. However, the search will be undermined if there are crossing boundaries. In such a case there will be more than two boundary directed edges at the vertex and there is no unique external angle.

Establishing the no crossing boundaries property ensures these failures cannot occur. We can easily check this property for newly created Tgraphs (with checkedTgraph) and the face adding operations cannot create crossing boundaries.

Touching Vertices and Crossing Boundaries

When a new face to be added on (a,b) has neither of the other two edges already in the Tgraph, the third vertex needs to be created. However it could already exist in the Tgraph – it is not on an edge coming from a or b but from another non-local part of the Tgraph. We call this a touching vertex. If we simply added a new vertex without checking for a clash this would create a non-sensible Tgraph. However, if we do check and find an existing vertex, we still cannot add the face using this because it would create a crossing boundary.

Our version of forcing prevents face additions that would create a touching vertex/crossing boundary by calculating the positions of boundary vertices.

No conflicting edges

There is a final (simple) check when adding a new face, to prevent a long edge (phiEdge) sharing with a short edge. This can arise if we force an incorrect Tgraph (as we will see later).

Implementing Forcing

Our order of forcing prioritises updates (face additions) which do not introduce a new vertex. Such safe updates are easy to recognise and they do not require a touching vertex check. Surprisingly, this pretty much removes the problem of touching vertices altogether.

As an illustration, consider foolDMinus again on the left of figure 11. Adding the left dart onto edge (12,14) is not a safe addition (and would create a crossing boundary at 6). However, adding the right dart RD(6,14,11) is safe and creates the new edge (6,14) which then makes the left dart addition safe. In fact it takes some contrivance to come up with a Tgraph with an update that could fail the check during forcing when safe cases are always done first. Figure 17 shows such a contrived Tgraph formed by removing the faces shown in green from a twice decomposed sun on the left. The forced result is shown on the right. When there are no safe cases, we need to try an unsafe one. The four green faces at the bottom are blocked by the touching vertex check. This leaves any one of 9 half-kites at the centre which would pass the check. But after just one of these is added, the check is not needed again. There is always a safe addition to be done at each step until all the green faces are added.

Boundary information

The implementation of forcing has been made more efficient by calculating some boundary information in advance. This boundary information uses a type BoundaryState

data BoundaryState
  = BoundaryState
    { boundary    :: [Dedge]
    , bvFacesMap  :: Mapping Vertex [TileFace]
    , bvLocMap    :: Mapping Vertex (Point V2 Double)
    , allFaces    :: [TileFace]
    , nextVertex  :: Vertex
    } deriving (Show)

This records the boundary directed edges (boundary) plus a mapping of the boundary vertices to their incident faces (bvFacesMap) plus a mapping of the boundary vertices to their positions (bvLocMap). It also keeps track of all the faces and the vertex number to use when adding a vertex. The boundary information is easily incremented for each face addition without being recalculated from scratch, and a final Tgraph with all the new faces is easily recovered from the boundary information when there are no more updates.

makeBoundaryState  :: Tgraph -> BoundaryState
recoverGraph  :: BoundaryState -> Tgraph

The saving that comes from using boundary information lies in efficient incremental changes to the boundary information and, of course, in avoiding the need to consider internal faces. As a further optimisation we keep track of updates in a mapping from boundary directed edges to updates, and supply a list of affected edges after an update so the update calculator (update generator) need only revise these. The boundary and mapping are combined in a ForceState.

type UpdateMap = Mapping Dedge Update
type UpdateGenerator = BoundaryState -> [Dedge] -> UpdateMap
data ForceState = ForceState 
       { boundaryState:: BoundaryState
       , updateMap:: UpdateMap 
       }

Forcing then involves using a specific update generator (allUGenerator) and initialising the state, then using the recursive forceAll which keeps doing updates until there are no more, before recovering the final Tgraph.

force:: Tgraph -> Tgraph
force = forceWith allUGenerator

forceWith:: UpdateGenerator -> Tgraph -> Tgraph
forceWith uGen = recoverGraph . boundaryState . 
                 forceAll uGen . initForceState uGen

forceAll :: UpdateGenerator -> ForceState -> ForceState
initForceState :: UpdateGenerator -> Tgraph -> ForceState

In addition to force we can easily define

wholeTiles:: Tgraph -> Tgraph
wholeTiles = forceWith wholeTileUpdates 

which just uses the first forcing rule to make sure every half-tile has a matching other half.

We also have a version of force which counts to a specific number of face additions.

stepForce :: Int -> ForceState -> ForceState

This proved essential in uncovering problems of accumulated inaccuracy in calculating boundary positions (now fixed).

Some Other Experiments

Below we describe results of some experiments using the tools introduced above. Specifically: emplacements, sub-Tgraphs, incorrect tilings, and composition choices.

Emplacements

The finite number of rules used in forcing are based on local boundary vertex and edge information only. We thought we may be able to improve on this by considering a composition and forcing at the next level up before decomposing and forcing again. This thus considers slightly broader local information. In fact we can iterate this process to all the higher levels of composition. Some Tgraphs produce an empty Tgraph when composed so we can regard those as maximal compositions. For example compose fool produces an empty Tgraph.

The idea was to take an arbitrary Tgraph and apply (compose . force) repeatedly to find its maximally composed (non-empty) Tgraph, before applying (force . decompose) repeatedly back down to the starting level (so the same number of decompositions as compositions).

We called the function emplace, and called the result the emplacement of the starting Tgraph as it shows a region of influence around the starting Tgraph.

With earlier versions of forcing when we had fewer rules, emplace g often extended force g for a Tgraph g. This allowed the identification of some new rules. However, since adding the new rules we have not found Tgraphs where the result of force had fewer faces than the result of emplace.

[As an important update, we have now found examples where the result of force strictly includes the result of emplace (modulo vertex renumbering).

Sub-Tgraphs

In figure 18 on the left we have a four times decomposed dart dartD4 followed by two sub-Tgraphs brokenDart and badlyBrokenDart which are constructed by removing faces from dartD4 (but retaining the connectedness condition and the no crossing boundaries condition). These all produce the same forced result (depicted middle row left in figure 15).

However, if we do compositions without forcing first we find badlyBrokenDart fails because it produces a graph with crossing boundaries after 3 compositions. So compose on its own is not always safe, where safe means guaranteed to produce a valid Tgraph from a valid correct Tgraph.

In other experiments we tried force on Tgraphs with holes and on incomplete boundaries around a potential hole. For example, we have taken the boundary faces of a forced, 5 times decomposed dart, then removed a few more faces to make a gap (which is still a valid Tgraph). This is shown at the top in figure 19. The result of forcing reconstructs the complete original forced graph. The bottom figure shows an intermediate stage after 2200 face additions. The gap cannot be closed off to make a hole as this would create a crossing boundary, but the channel does get filled and eventually closes the gap without creating a hole.

Incorrect Tilings

When we say a Tgraph g is correct (respectively: incorrect), we mean g represents a correct tiling (respectively: incorrect tiling). A simple example of an incorrect Tgraph is a kite with a dart on each side (referred to as a mistake by Penrose) shown on the left of figure 20.

*Main> mistake
Tgraph [RK (1,2,4),LK (1,3,2),RD (3,1,5)
       ,LD (4,6,1),LD (3,5,7),RD (4,8,6)
       ]

If we try to force (or emplace) this Tgraph it produces an error in construction which is detected by the test for conflicting edge types (a phiEdge sharing with a non-phiEdge).

*Main> force mistake
... *** Exception: doUpdate:(incorrect tiling)
Conflicting new face RK (11,1,6)
with neighbouring faces
[RK (9,1,11),LK (9,5,1),RK (1,2,4),LK (1,3,2),RD (3,1,5),LD (4,6,1),RD (4,8,6)]
in boundary
BoundaryState ...

In figure 20 on the right, we see that after successfully constructing the two whole kites on the top dart short edges, there is an attempt to add an RK on edge (1,6). The process finds an existing edge (1,11) in the correct direction for one of the new edges so tries to add the erroneous RK (11,1,6) which fails a noConflicts test.

So it is certainly true that incorrect Tgraphs may fail on forcing, but forcing cannot create an incorrect Tgraph from a correct Tgraph.

If we apply decompose to mistake it produces another incorrect Tgraph (which is similarly detected if we apply force), but will nevertheless still compose back to mistake if we do not try to force.

Interestingly, though, the incorrectness of a Tgraph is not always preserved by decompose. If we start with mistake1 which is mistake with just two of the half darts (and also incorrect) we still get a similar failure on forcing, but decompose mistake1 is no longer incorrect. If we apply compose to the result or force then compose the mistake is thrown away to leave just a kite (see figure 21). This is an example where compose is not a left inverse to either decompose or (force . decompose).

Composing with Choices

We know that unknowns indicate possible choices (although some choices may lead to incorrect Tgraphs). As an experiment we introduce

makeChoices :: Tgraph -> [Tgraph]

which produces alternatives for the 2 choices of each of unknowns (prior to composing). This uses forceLDB which forces an unknown to be a largeDartBase by adding an appropriate joined half dart at the node, and forceLKC which forces an unknown to be a largeKiteCentre by adding a half dart and a whole kite at the node (making up the 3 pieces for a larger half kite).

Figure 22 illustrates the four choices for composing fool this way. The top row has the four choices of makeChoices fool (with the fool shown embeded in red in each case). The bottom row shows the result of applying compose to each choice.

In this case, all four compositions are correct tilings. The problem is that, in general, some of the choices may lead to incorrect tilings. More specifically, a choice of one unknown can determine what other unknowns have to become with constraints such as

• a and b have to be opposite choices
• a and b have to be the same choice
• a and b cannot both be largeKiteCentres
• a and b cannot both be largeDartBases

This analysis of constraints on unknowns is not trivial. The potential exponential results from choices suggests we should compose and force as much as possible and only consider unknowns of a maximal Tgraph.

For calculating the emplacement of a Tgraph, we first find the forced maximal Tgraph before decomposing. We could also consider using makeChoices at this top step when there are unknowns, i.e a version of emplace which produces these alternative results (emplaceChoices)

The result of emplaceChoices is illustrated for foolD in figure 23. The first force and composition is unique producing the fool level at which point we get 4 alternatives each of which compose further as previously illustrated in figure 22. Each of these are forced, then decomposed and forced, decomposed and forced again back down to the starting level. In figure 23 foolD is overlaid on the 4 alternative results. What they have in common is (as you might expect) emplace foolD which equals force foolD and is the graph shown on the right of figure 11.

Future Work

I am collaborating with Stephen Huggett who suggested the use of graphs for exploring properties of the tilings. We now have some tools to experiment with but we would also like to complete some formalisation and proofs.

It would also be good to establish whether it is true that g is incorrect iff force g fails.

We have other conjectures relating to subgraph ordering of Tgraphs and Galois connections to explore.

Revised version (no longer the full program in this literate Haskell)

Infinite non-periodic tessellations of Roger Penrose’s kite and dart tiles.

As part of a collaboration with Stephen Huggett, working on some mathematical properties of Penrose tilings, I recognised the need for quick renderings of tilings. I thought Haskell diagrams would be helpful here, and that turned out to be an excellent choice. Two dimensional vectors were well-suited to describing tiling operations and these are included as part of the diagrams package.

This literate Haskell uses the Haskell diagrams package to draw tilings with kites and darts. It also implements the main operations of compChoices and decompPatch which are used for constructing tilings (explained below).

Firstly, these 5 lines are needed in Haskell to use the diagrams package:

{-# LANGUAGE NoMonomorphismRestriction #-}
{-# LANGUAGE FlexibleContexts          #-}
{-# LANGUAGE TypeFamilies              #-}

import Diagrams.Prelude
import Diagrams.Backend.SVG.CmdLine

and we will also import a module for half tiles (explained later)

import HalfTile

Legal tilings

These are the kite and dart tiles.

The red line marking here on the right hand copies, is purely to illustrate rules about how tiles can be put together for legal (non-periodic) tilings. Obviously edges can only be put together when they have the same length. If all the tiles are marked with red lines as illustrated on the right, the vertices where tiles meet must all have a red line or none must have a red line at that vertex. This prevents us from forming a simple rombus by placing a kite top at the base of a dart and thus enabling periodic tilings.

All edges are powers of the golden section which we write as phi.

phi::Double
phi = (1.0 + sqrt 5.0) / 2.0

So if the shorter edges are unit length, then the longer edges have length phi. We also have the interesting property of the golden section that and so , and .

All angles in the figures are multiples of tt which is 36 deg or 1/10 turn. We use ttangle to express such angles (e.g 180 degrees is ttangle 5).

ttangle:: Int -> Angle Double
ttangle n = (fromIntegral (n `mod` 10))*^tt
             where tt = 1/10 @@ turn

Pieces

In order to implement compChoices and decompPatch, we need to work with half tiles. We now define these in the separately imported module HalfTile with constructors for Left Dart, Right Dart, Left Kite, Right Kite

data HalfTile rep = LD rep -- defined in HalfTile module
                  | RD rep
                  | LK rep
                  | RK rep

where rep is a type variable allowing for different representations. However, here, we want to use a more specific type which we will call Piece:

type Piece = HalfTile (V2 Double)

where the half tiles have a simple 2D vector representation to provide orientation and scale. The vector represents the join edge of each half tile where halves come together. The origin for a dart is the tip, and the origin for a kite is the acute angle tip (marked in the figure with a red dot).

These are the only 4 pieces we use (oriented along the x axis)

ldart,rdart,lkite,rkite:: Piece
ldart = LD unitX
rdart = RD unitX
lkite = LK (phi*^unitX)
rkite = RK (phi*^unitX)

Perhaps confusingly, we regard left and right of a dart differently from left and right of a kite when viewed from the origin. The diagram shows the left dart before the right dart and the left kite before the right kite. Thus in a complete tile, going clockwise round the origin the right dart comes before the left dart, but the left kite comes before the right kite.

When it comes to drawing pieces, for the simplest case, we just want to show the two tile edges of each piece (and not the join edge). These edges are calculated as a list of 2 new vectors, using the join edge vector v. They are ordered clockwise from the origin of each piece

pieceEdges:: Piece -> [V2 Double]
pieceEdges (LD v) = [v',v ^-^ v'] where v' = phi*^rotate (ttangle 9) v
pieceEdges (RD v) = [v',v ^-^ v'] where v' = phi*^rotate (ttangle 1) v
pieceEdges (RK v) = [v',v ^-^ v'] where v' = rotate (ttangle 9) v
pieceEdges (LK v) = [v',v ^-^ v'] where v' = rotate (ttangle 1) v

Now drawing lines for the 2 outer edges of a piece is simply

drawPiece:: Piece -> Diagram B
drawPiece = strokeLine . fromOffsets . pieceEdges

and drawing all 3 edges round a piece is

drawRoundPiece:: Piece -> Diagram B
drawRoundPiece = strokeLoop . closeLine . fromOffsets . pieceEdges

To fill half tile pieces, we can use fillOnlyPiece which fills without showing edges of a half tile (by using line width none).

fillOnlyPiece:: Colour Double -> Piece -> Diagram B
fillOnlyPiece col piece = drawRoundPiece piece # fc col # lw none

We also use fillPieceDK which fills darts and kites with given colours and also draws edges using drawPiece.

fillPieceDK:: Colour Double -> Colour Double -> Piece -> Diagram B
fillPieceDK dcol kcol piece = drawPiece piece <> fillOnlyPiece col piece where
    col = case piece of (LD _) -> dcol
                        (RD _) -> dcol
                        (LK _) -> kcol
                        (RK _) -> kcol

For an alternative fill operation on whole tiles, it is useful to calculate a list of the 4 tile edges of a completed half-tile piece clockwise from the origin of the tile. (This will allow colour filling a whole tile)

wholeTileEdges:: Piece -> [V2 Double]
wholeTileEdges (LD v) = pieceEdges (RD v) ++ map negated (reverse (pieceEdges (LD v)))
wholeTileEdges (RD v) = wholeTileEdges (LD v)
wholeTileEdges (LK v) = pieceEdges (LK v) ++ map negated (reverse (pieceEdges (RK v)))
wholeTileEdges (RK v) = wholeTileEdges (LK v)

To fill whole tiles with colours, darts with dcol and kites with kcol we can now use leftFillPieceDK. This uses only the left pieces to identify the whole tile and ignores right pieces so that a tile is not filled twice.

leftFillPieceDK:: Colour Double -> Colour Double -> Piece -> Diagram B
leftFillPieceDK dcol kcol c = case c of 
  (LD _) -> (strokeLoop $ glueLine $ fromOffsets $ wholeTileEdges c)  # fc dcol
  (LK _) -> (strokeLoop $ glueLine $ fromOffsets $ wholeTileEdges c)  # fc kcol
  _      -> mempty

By making Pieces transformable we can reuse generic transform operations. These 4 lines of code are required to do this

type instance N (HalfTile a) = N a
type instance V (HalfTile a) = V a
instance Transformable a => Transformable (HalfTile a) where
    transform t ht = fmap (transform t) ht

So we can also scale and rotate a piece by an angle. (Positive rotations are in the anticlockwise direction.)

scale :: Double -> Piece -> Piece
rotate :: Angle Double -> Piece -> Piece

Patches

A patch is a list of located pieces (each with a 2D point)

type Patch = [Located Piece]

To turn a whole patch into a diagram using some function pd for drawing the pieces, we use

drawPatchWith:: (Piece -> Diagram B) -> Patch -> Diagram B 
drawPatchWith pd patch = position $ fmap (viewLoc . mapLoc pd) patch

Here mapLoc applies a function to the piece in a located piece – producing a located diagram in this case, and viewLoc returns the pair of point and diagram from a located diagram. Finally position forms a single diagram from the list of pairs of points and diagrams.

Update: We now use a class for drawable tilings, making Patch an instance

class Drawable a where
 drawWith :: (Piece -> Diagram B) -> a -> Diagram B

instance Drawable Patch where
 drawWith = drawPatchWith

We then introduce special cases:

draw :: Drawable a => a -> Diagram B
draw = drawWith drawPiece

fillDK:: Drawable a => Colour Double -> Colour Double -> a -> Diagram B
fillDK c1 c2 = drawWith (fillPieceDK c1 c2)

Patches are automatically inferred to be transformable now Pieces are transformable, so we can also scale a patch, translate a patch by a vector, and rotate a patch by an angle.

scale :: Double -> Patch -> Patch
rotate :: Angle Double -> Patch -> Patch
translate:: V2 Double -> Patch -> Patch

As an aid to creating patches with 5-fold rotational symmetry, we combine 5 copies of a basic patch (rotated by multiples of ttangle 2 successively).

penta:: Patch -> Patch
penta p = concatMap copy [0..4] 
            where copy n = rotate (ttangle (2*n)) p

This must be used with care to avoid nonsense patches. But two special cases are

sun,star::Patch         
sun =  penta [rkite `at` origin, lkite `at` origin]
star = penta [rdart `at` origin, ldart `at` origin]

This figure shows some example patches, drawn with draw The first is a star and the second is a sun.

The tools so far for creating patches may seem limited (and do not help with ensuring legal tilings), but there is an even bigger problem.

Correct Tilings

Unfortunately, correct tilings – that is, tilings which can be extended to infinity – are not as simple as just legal tilings. It is not enough to have a legal tiling, because an apparent (legal) choice of placing one tile can have non-local consequences, causing a conflict with a choice made far away in a patch of tiles, resulting in a patch which cannot be extended. This suggests that constructing correct patches is far from trivial.

The infinite number of possible infinite tilings do have some remarkable properties. Any finite patch from one of them, will occur in all the others (infinitely many times) and within a relatively small radius of any point in an infinite tiling. (For details of this see links at the end)

This is why we need a different approach to constructing larger patches. There are two significant processes used for creating patches, namely inflate (also called compose) and decompose.

To understand these processes, take a look at the following figure.

Here the small pieces have been drawn in an unusual way. The edges have been drawn with dashed lines, but long edges of kites have been emphasised with a solid line and the join edges of darts marked with a red line. From this you may be able to make out a patch of larger scale kites and darts. This is an inflated patch arising from the smaller scale patch. Conversely, the larger kites and darts decompose to the smaller scale ones.

Decomposition

Since the rule for decomposition is uniquely determined, we can express it as a simple function on patches.

decompPatch :: Patch -> Patch
decompPatch = concatMap decompPiece

where the function decompPiece acts on located pieces and produces a list of the smaller located pieces contained in the piece. For example, a larger right dart will produce both a smaller right dart and a smaller left kite. Decomposing a located piece also takes care of the location, scale and rotation of the new pieces.

decompPiece lp = case viewLoc lp of
  (p, RD vd)-> [ LK vd  `at` p
               , RD vd' `at` (p .+^ v')
               ] where v'  = phi*^rotate (ttangle 1) vd
                       vd' = (2-phi) *^ (negated v') -- (2-phi) = 1/phi^2
  (p, LD vd)-> [ RK vd `at` p
               , LD vd' `at` (p .+^ v')
               ]  where v'  = phi*^rotate (ttangle 9) vd
                        vd' = (2-phi) *^ (negated v')  -- (2-phi) = 1/phi^2
  (p, RK vk)-> [ RD vd' `at` p
               , LK vk' `at` (p .+^ v')
               , RK vk' `at` (p .+^ v')
               ] where v'  = rotate (ttangle 9) vk
                       vd' = (2-phi) *^ v' -- v'/phi^2
                       vk' = ((phi-1) *^ vk) ^-^ v' -- (phi-1) = 1/phi
  (p, LK vk)-> [ LD vd' `at` p
               , RK vk' `at` (p .+^ v')
               , LK vk' `at` (p .+^ v')
               ] where v'  = rotate (ttangle 1) vk
                       vd' = (2-phi) *^ v' -- v'/phi^2
                       vk' = ((phi-1) *^ vk) ^-^ v' -- (phi-1) = 1/phi

This is illustrated in the following figure for the cases of a right dart and a right kite.

The symmetric diagrams for left pieces are easy to work out from these, so they are not illustrated.

With the decompPatch operation we can start with a simple correct patch, and decompose repeatedly to get more and more detailed patches. (Each decomposition scales the tiles down by a factor of but we can rescale at any time.)

This figure illustrates how each piece decomposes with 4 decomposition steps below each one.

thePieces =  [ldart, rdart, lkite, rkite]  
fourDecomps = hsep 1 $ fmap decomps thePieces # lw thin where
        decomps pc = vsep 1 $ fmap draw $ take 5 $ decompositionsP [pc `at` origin] 

We have made use of the fact that we can create an infinite list of finer and finer decompositions of any patch, using:

decompositionsP:: Patch -> [Patch]
decompositionsP = iterate decompPatch

We could get the n-fold decomposition of a patch as just the nth item in a list of decompositions.

For example, here is an infinite list of decomposed versions of sun.

suns = decompositionsP sun

The coloured tiling shown at the beginning is simply 6 decompositions of sun displayed using leftFillPieceDK

leftFilledSun6 :: Diagram B
leftFilledSun6 = drawWith (leftFillPieceDK red blue) (suns !!6) # lw thin

The earlier figure illustrating larger kites and darts emphasised from the smaller ones is also suns!!6 but this time pieces are drawn with experiment.

experimentFig = drawWith experiment (suns!!6) # lw thin

experiment:: Piece -> Diagram B
experiment pc = emph pc <> (drawRoundPiece pc # dashingN [0.002,0.002] 0 # lw ultraThin)
  where emph pc = case pc of
          (LD v) -> (strokeLine . fromOffsets) [v] # lc red   -- emphasise join edge of darts in red
          (RD v) -> (strokeLine . fromOffsets) [v] # lc red 
          (LK v) -> (strokeLine . fromOffsets) [rotate (ttangle 1) v] -- emphasise long edge for kites
          (RK v) -> (strokeLine . fromOffsets) [rotate (ttangle 9) v]

Compose Choices

You might expect composition (also called inflation) to be a kind of inverse to decomposition, but it is a bit more complicated than that. With our current representation of pieces, we can only compose single pieces. This amounts to embedding the piece into a larger piece that matches how the larger piece decomposes. There is thus a choice at each inflation step as to which of several possibilities we select as the larger half-tile. We represent this choice as a list of alternatives. This list should not be confused with a patch. It only makes sense to select one of the alternatives giving a new single piece.

The earlier diagram illustrating how decompositions are calculated also shows the two choices for embedding a right dart into either a right kite or a larger right dart. There will be two symmetric choices for a left dart, and three choices for left and right kites.

Once again we work with located pieces to ensure the resulting larger piece contains the original in its original position in a decomposition.

compChoices :: Located Piece -> [Located Piece]
compChoices lp = case viewLoc lp of
  (p, RD vd)-> [ RD vd' `at` (p .+^ v')
               , RK vk  `at` p
               ] where v'  = (phi+1) *^ vd                  -- vd*phi^2
                       vd' = rotate (ttangle 9) (vd ^-^ v')
                       vk  = rotate (ttangle 1) v'
  (p, LD vd)-> [ LD vd' `at` (p .+^ v')
               , LK vk `at` p
               ] where v'  = (phi+1) *^ vd                  -- vd*phi^2
                       vd' = rotate (ttangle 1) (vd ^-^ v')
                       vk  = rotate (ttangle 9) v'
  (p, RK vk)-> [ LD vk  `at` p
               , LK lvk' `at` (p .+^ lv') 
               , RK rvk' `at` (p .+^ rv')
               ] where lv'  = phi*^rotate (ttangle 9) vk
                       rv'  = phi*^rotate (ttangle 1) vk
                       rvk' = phi*^rotate (ttangle 7) vk
                       lvk' = phi*^rotate (ttangle 3) vk
  (p, LK vk)-> [ RD vk  `at` p
               , RK rvk' `at` (p .+^ rv')
               , LK lvk' `at` (p .+^ lv')
               ] where v0 = rotate (ttangle 1) vk
                       lv'  = phi*^rotate (ttangle 9) vk
                       rv'  = phi*^rotate (ttangle 1) vk
                       rvk' = phi*^rotate (ttangle 7) vk
                       lvk' = phi*^rotate (ttangle 3) vk

As the result is a list of alternatives, we need to select one to do further inflations. We can express all the alternatives after n steps as compNChoices n where

compNChoices :: Int -> Located Piece -> [Located Piece]
compNChoices 0 lp = [lp]
compNChoices n lp = do
    lp' <- inflate lp
    inflations (n-1) lp'

This figure illustrates 5 consecutive choices for inflating a left dart to produce a left kite. On the left, the finishing piece is shown with the starting piece embedded, and on the right the 5-fold decomposition of the result is shown.

fiveCompChoices = hsep 1 $ [ draw [ld] <> draw [lk']
                           , draw (decompositionsP [lk'] !!5)
                           ] where
  ld  = (ldart `at` origin)
  lk  = compChoices ld  !!1
  rk  = compChoices lk  !!1
  rk' = compChoices rk  !!2
  ld' = compChoices rk' !!0
  lk' = compChoices ld' !!1

Finally, at the end of this literate haskell program we choose which figure to draw as output.

fig :: Diagram B
fig = leftFilledSun6

main = mainWith fig

That’s it. But, What about composing whole patches?, I hear you ask. Unfortunately we need to answer questions like what pieces are adjacent to a piece in a patch and whether there is a corresponding other half for a piece. These cannot be done with our simple vector representations. We would need some form of planar graph representation, which is much more involved. That is another story.

Many thanks to Stephen Huggett for his inspirations concerning the tilings. A library version of the above code is available on GitHub

Further reading on Penrose Tilings

As well as the Wikipedia entry Penrose Tilings I recommend two articles in Scientific American from 2005 by David Austin Penrose Tiles Talk Across Miles and Penrose Tilings Tied up in Ribbons.

There is also a very interesting article by Roger Penrose himself: Penrose R Tilings and quasi-crystals; a non-local growth problem? in Aperiodicity and Order 2, edited by Jarich M, Academic Press, 1989.

More information about the diagrams package can be found from the home page Haskell diagrams

Introduction

We implement the multigrid and full multigrid schemes using the Haskell parallel array library Repa. Whilst this file is a literate Haskell program it omits some preliminaries. The full code can be found on GitHub at multiGrid.

Repa (short for ‘regular parallel arrays’) is a Haskell library for efficiently calculating with arrays functionally. It allows parallel computation to be expressed with a single primitive (computeP). This is possible because we only use immutable arrays so calculations are deterministic without any need to worry about locks/mutual exclusion/race conditions/etc.

The multigrid and full multigrid schemes are used for approximating solutions of partial differential equations but they are strategies rather than solvers. They can be used to drastically improve more basic solvers provided convergence and error analysis are considered for the basic solver.

We only consider linear partial differential equations and the Poisson equation in particular where linear partial differential equations can be summarised in the form

Here, is a linear operator (such as the laplace operator or higher order or lower order partial derivatives or any linear combination of these), is given and is the function we want to solve for, and both are defined on a region . The multigrid scheme starts with a fine grid (where is the grid spacing) on which we want to obtain an approximate solution by finite difference methods. The scheme involves the use of a coarser grid (e.g. ) and, recursively, a stack of coarser and coarser grids to apply error correction on approximations. The full multigrid scheme also uses coarser grids to improve initial approximations used by multigrid.

Outline of Coarse Grid Error Correction

Assume we have a basic method to approximate solutions which we call ‘approximately solve’. Then the scheme for coarse grid correction of errors is

1. Take an initial guess on the fine grid () and use it to approximately solve to get a new approximation .

   We want to estimate errors

   

   We do not have to calculate them, but the errors satisfy

   

   The negative of these values are called the residuals () and these we can calculate.

2. Compute residuals

3. Move to the coarse grid and (recursively) solve

   (This is a problem of the same form as we started with, but it requires solving with restrictions of and for the coarse grid.)

   We use zeros as initial guess for errors and this recursion results in an approximation for errors ,

4. Interpolate coarse errors () into the fine grid to get and get a corrected guess

   

5. Approximately solve again (on the fine grid), but now starting with to get a final approximation.

So a basic requirement of the multigrid scheme is to be able to move values to a coarser grid (restriction) and from a coarser grid to a finer grid (interpolation). We will write functions for these below. First we give a quick summary of Repa and operations we will be using, which experienced Repa users may prefer to skip.

Background Repa Summary

Repa array types have an extent (shape) and an indication of representation as well as a content type. For example Array U DIM2 Double is the type of a two dimensional array of Doubles. The U indicates a ‘manifest’ representation which means that the contents are fully calculated rather than delayed. A delayed representation is expressed with D rather than U. The only other representation type we use explicitly is (TR PC5) which is the type of a partitioned, array resulting from a stencil mapping operation. Non-manifest arrays are useful for enabling the compiler to combine (fuse) operations to optimise after inlining code. Some operations require the argument array to be manifest, so to make a delayed array manifest we can use

• computeP which employes parallel evaluation on the array elements. The type of computeP requires that it is used within a monad. This is to prevent us accidentally writing nested calls of computeP.

The extent DIM2 abbreviates Z :. Int :. Int (where Z represents a base shape for the empty zero dimensional array). The integers give the sizes in each dimension. We can have higher dimensions (adding more with :.) but will only be using DIM2 here. Values of an extent type are used to express the sizes of an array in each dimensions and also used as indexes into an array where indexing starts at 0.

We only use a few other Repa operations

• fromListUnboxed takes an extent and a list of values to create a new manifest array.

• traverse takes 3 arguments to produce a new delayed array. The first is the (old) array. The second is a function which produces the extent of the new array when given the extent of the old array. The third is a function to calculate any item in the new array [when supplied with an operation (get) to retrieve values from the old array and an index in the new array for the item to calculate].

• szipWith is analgous to zipWith for lists, taking a binary operation to combine two arrays. (There is also a zipWith for arrays, but the szipWith version is tuned for partitioned arrays)

• smap is analgous to map for lists, mapping an operation over an array. (There is also a map for arrays, but the smap version is tuned for partitioned arrays)

• mapStencil2 takes a boundary handling option (we only use BoundClamp here), a stencil, and a two dimensional array. It maps (convolves) the stencil over the array to produce a new (partititioned) array.

There is a convenient way of writing down stencils which are easy to visualize. We will see examples later, but this format requires pragmas

> {-# LANGUAGE TemplateHaskell     #-}
> {-# LANGUAGE QuasiQuotes         #-}

Interpolation (Prolongation) and Restriction

To interpolate from coarse to fine we want a (linear) mapping that is full rank. E.g identity plus averaging neighbours.

We can restrict by injection or take some weighting of neighbours as well. A common requirement is that this is also linear and a constant times the transpose of the interpolation when these are expressed as matrix operations.

Repa examples of restriction and interpolation in 2D.

We will be using DIM2 grids and work with odd extents so that border points remain on the border when moving between grids. A common method for restriction is to take a weighted average of fine neighbours with each coarse grid point and we can easily express this as a stencil mapping. We calculate one sixteenth of the results from mapping this stencil

> restStencil :: Stencil DIM2 Double
> restStencil =   [stencil2|  1 2 1         
>                             2 4 2 
>                             1 2 1 |]

I.e we map this stencil then divide by 16 and take the coarse array of items where both the row and column are even. Taking the coarse grid items can be done with an array traversal

> {-# INLINE coarsen #-}  
> coarsen :: Array U DIM2 Double -> Array D DIM2 Double
> coarsen !arr 
>  = traverse arr   -- i+1 and j+1 to deal with odd extents for arr correctly
>             (\ (e :. i :. j) -> (e :. (i+1) `div` 2 :. (j+1) `div` 2))
>             (\ get (e :. i :. j) -> get (e :. 2*i :. 2*j))

Here the second argument for traverse – the function to calculate the new extent from the old – adds one before the division to ensures that odd extents are dealt with appropriately.
The inline pragma and bang pattern (!) on the argument are needed for good optimisation.

Coarsening after mapping the above stencil works well with odd extents but is not so good for even extents. For completeness we allow for even extents as well and in this case it is more appropriate to calculate one quarter of the results from mapping this stencil

> sum4Stencil :: Stencil DIM2 Double    
> sum4Stencil = [stencil2|  0 0 0         
>                           0 1 1 
>                           0 1 1 |] 

We treat mixed even and odd extents this way as well so we can express restriction for all cases as

> {-# INLINE restrict #-}       
> restrict :: Array U DIM2 Double -> Array D DIM2 Double
> restrict !arr 
>   | odd n && odd m
>     = coarsen
>       $ smap (/16)  
>       $ mapStencil2 BoundClamp restStencil arr
>   | otherwise
>     = coarsen
>       $ smap (/4)  
>       $ mapStencil2 BoundClamp sum4Stencil arr
>   where _ :. m :. n = extent arr

For interpolation in the case of odd extents we want to distribute coarse to fine values according to the “give-to” stencil

1/4 1/2 1/4
1/2  1  1/2
1/4 1/2 1/4

This means that the fine value at the centre becomes the coarse value at the corresponding position (if you picture the coarse grid spaced out to overlay the fine grid). The fine values around this central value inherit proportions of the coarse value as indicated by the give-to stencil (along with proportions from other neighbouring coarse values).

Odd Extent Interpolation

For the even extent case we simply make four copies of each coarse value using a traverse

> {-# INLINE inject4 #-} 
> inject4 :: Source r a => Array r DIM2 a -> Array D DIM2 a
> inject4 !arr 
>  = traverse arr -- mod 2s deal with odd extents
>             (\ (e :. i :. j) -> (e :. 2*i - (i `mod` 2) :. 2*j - (j `mod` 2)))
>             (\get (e :. i :. j) -> get(e :. i `div` 2 :. j `div` 2))

Even Extent Interpolation

Again, we just use the latter version to cover cases with mixed even and odd extents. There is no primitive for “give-to” stencils but it is easy enough to define what we want with a traverse.

> {-# INLINE interpolate #-} 
> interpolate :: Array U DIM2 Double -> Array D DIM2 Double
> interpolate !arr 
>   | odd m && odd n
>        = traverse arr 
>                   (\ (e :. i :. j) -> (e :. 2*i - (i `mod` 2) :. 2*j - (j `mod` 2)))
>                   (\get (e :. i :. j) -> case () of
>                    _ | even i && even j     -> get(e :. i `div` 2 :. j `div` 2)
>                      | even i               -> (0.5)*(get(e :. i `div` 2 :. j `div` 2)
>                                                       + get(e :. i `div` 2 :. (j `div` 2)+1)) 
>                       -- odd i
>                      | even j               -> (0.5)*(get(e :. i `div` 2 :. j `div` 2)
>                                                       + get(e :. (i `div` 2)+1 :. j `div` 2)) 
>                       -- odd i and j
>                      | otherwise            -> (0.25)*(get(e :. i `div` 2 :. j `div` 2)
>                                                       + get(e :. i `div` 2 :. (j `div` 2)+1)
>                                                       + get(e :. (i `div` 2)+1 :. j `div` 2)
>                                                       + get(e :. (i `div` 2)+1 :. (j `div` 2)+1))
>                   )
>   | otherwise = inject4 arr
>  where _ :. n :. m = extent arr

The uses of mod 2 in the new size calculation are there to ensure that odd extents remain odd.

An alternative interpolation calculation

As an aside, we found that the above interpolation for odd extents can be implemented with a stencil convolution using sum4Stencil (which we defined above) after applying inject4 and then dividing by 4.
So we could define (for odd extents)

interpolate' :: Monad m
                => Array U DIM2 Double -> m (Array (TR PC5) DIM2 Double) 

interpolate' arr = 
   do fineArr <- computeP $ inject4 arr
      return $ smap (/4) 
             $ mapStencil2 BoundClamp sum4Stencil fineArr

We have to make the intermediate fineArr manifest (using computeP) before mapping the stencil over it which is why a monad is used. The final array is not manifest but a structured array resulting from the stencil map. By not making this manifest, we allow the computation to be combined with other computations improving inlining optimisation opportunities.

To illustrate the effect of interpolate', suppose the coarse array content is just

a b c
d e f
g h i

Then the injected array content will be

a a b b c c
a a b b c c
d d e e f f
d d e e f f
g g h h i i
g g h h i i

After mapping the stencil and dividing by 4 we have (assuming top left is (0,0))

at (1,1) (a+b+d+e)/4               (odd row,  odd column)
at (2,1) (d+e+d+e)/4 = (d+e)/2     (even row, odd column)
at (1,2) (b+b+e+e)/4 = (b+e)/2     (odd row,  even column)
at (2,2) (e+e+e+e)/4 = e           (even row, even column)

This even handles the bottom and right boundaries as in the original interpolation.

Slightly more generally, after inject4 and for any w x y z then convolution with the stencil

0 0 0
0 w x
0 y z

will produce the same as interpolating with the “give-to” stencil

  z     (z+y)    y   
(z+x) (z+y+x+w) (y+w)           
  x     (x+w)    w      

Conversely after inject4, the give-to stencil has to have this form for a stencil convolution to produce the same results.

Since the stencil version does require making an intermediate array manifest it is not clear at face value which is more efficient, so we will stick with interpolate.

We note that for even extents, the combination of an interpolation followed by a restriction is an identity (inject4 followed by coarsen). For odd extents, the combination preserves internal values but not values on the boarder. This will not be a problem if boundary conditions of the problem are used to update boarders of the array.

MultiGrid with Coarse Grid Correction Scheme

The problem to solve for is

where is a linear operator. This linear operator could be represented as a matrix or as a matrix multiplication, but this is not the most efficient way of implementing a solver. In the multigrid scheme described above we need knowledge of to implement an approximate solver, and also to calculate residuals. There are ways to implement these latter two operations with some mathematical analysis of the linear operator using finite differences. We will therefore be using an approximator and a residual calculator directly rather than a direct representation of as a parameter.

Let us assume that the approximate solver is a single step improvement which we want to iterate a few times. We can write down the iterator (iterateSolver)

> {-# INLINE iterateSolver #-}
> iterateSolver !opA !steps !arrInit
>  = go steps arrInit
>   where
>     go 0 !arr = return arr
>     go n !arr 
>       = do   arr' <- computeP $ opA arr
>              go (n - 1) arr'

In the definition of multiGrid we make use of a function to create an array of zeros of a given shape (extent)

> {-# INLINE zeroArray #-}
> zeroArray  :: Shape sh => sh -> Array U sh Double
> zeroArray !sh = fromListUnboxed sh $ replicate (size sh) 0.0

and a function to determine when we have the coarsest grid (the bases case for odd extents will be 3 by 3)

> {-# INLINE coarsest #-}
> coarsest :: Array U DIM2 Double -> Bool
> coarsest !arr = m<4 where (Z :. _ :. m) = extent arr

We also use some global parameters to indicate the number of iteration steps (to simplify expressions by passing fewer arguments)

> steps1,steps2 :: Int
> steps1 = 5       -- number of iterations for first step in multiGrid
> steps2 = 5       -- number of iterations for last step in multiGrid

To write down a first version of multiGrid, we temporarily ignore boundary conditions. We will also assume that the grid spacing is the same in each dimension and use just one parameter (h) as the grid spacing.

{- version ignoring boundary conditions -}
multiGrid approxOp residualOp h f uInit 
 = if coarsest uInit
   then 
     do computeP $ approxOp h f uInit
   else
     do v  <- iterateSolver (approxOp h f) steps1 uInit 
        r  <- computeP $ residualOp h f v                 -- calculate fine residuals
        r' <- computeP $ restrict r                       -- move to coarser grid                          
        err <- multiGrid approxOp residualOp (2*h) r' $ zeroArray (extent r')
        vC <- computeP $ szipWith (+) v  
                       $ interpolate err                  -- correct with errors on fine grid
        iterateSolver (approxOp h f) steps2 vC            -- solve again with improved approximation

The parameters are

1. approxOp – an approximate solver to be iterated.
2. residualOp – a residual calculator
3. h – the grid spacing
4. f – a representation for the (negative of) the function on the right hand side of the equation
5. uInit – an array representing a suitable first guess to start from.

Note that both approxOp and residualOp need to be passed h and f as well as an array when they are used. Also, the recursive call of multiGrid requres the two function parameters to work at each coarseness of grid. The grid spacing doubles and the array of residuals has to be converted to the coarser grid for the recursive call.

Next we consider how we deal with the boundary conditions. We will only deal with the Dirichlet type of boundary conditions here.

MultiGrid with Dirichlet boundary conditions

For as above, Dirichlet boundary conditions have the form for some function on the boundary . We will take the standard approach of representing the boundary conditions using a mask array (with 0’s indicating boundary positions and 1’s at non-boundary positions) along with a boundary values array which has 0’s at non-boundary positions. This will enable us to adapt these two arrays for different coarsenesses of grid.

We are assuming we are working in just two dimensions so . We can thus assume two dimensional arrays represent and the boundary mask and boundary values.

We can now define multiGrid with boundary conditions as:

 multiGrid approxOp residualOp h f boundMask boundValues uInit 
  = if coarsest uInit
    then 
      do computeP $ approxOp h f boundMask boundValues uInit
    else
      do v  <- iterateSolver (approxOp h f boundMask boundValues) steps1 uInit 
         r  <- computeP $ residualOp h f boundMask v      -- calculate fine residuals
         boundMask' <- computeP $ coarsen boundMask       -- move to coarser grid
         r' <- computeP $ szipWith (*) boundMask' $ restrict r                            
         let zeros = zeroArray (extent r')
         err <- multiGrid approxOp residualOp (2*h) r' boundMask' zeros zeros
         vC <- computeP $ szipWith (+) v  
                        $ szipWith (*) boundMask 
                        $ interpolate err                 -- correct with errors on fine grid
         iterateSolver (approxOp h f boundMask boundValues) steps2 vC 

In the base case (when the grid size is 3*3) there is only one value at the centre point which needs to be evaluated (assuming the surrounding 8 points are given by the boundary conditions). This will be solved exactly with a single step of the approximate solver. The recursive call uses zeros for both the boundValues (redundantly) as well the initial guess. We note that the residual calculation makes use of the mask but not the boundary values as these should be zero for residuals. The restriction of residuals for the coarse grid also requires applying the mask adapted for the coarse grid. The interpolation of errors uses the (fine version of the) mask but the boundary values will already be set in v which is added to the errors.

Full MultiGrid

Before looking at specific examples, we can now write down the algorithm for the full multigrid scheme which extends multigrid as follows.

Before calling multiGrid we precalculate the initial guess using a coarser grid and interpolate the result. The precalculation on the coarser grid involves the same full multigrid process recursively on coarser and coarser grids.

fullMG approxOp residualOp h f boundMask boundValues
 = if coarsest boundValues
     then do computeP $ approxOp h f boundMask boundValues boundValues
         -- an approximation with  1 interior point will be exact using boundValues as initial
     else do  
         -- recursively precalculate for the coarser level
            f' <- computeP $ restrict f
            boundMask' <- computeP $ coarsen boundMask
            boundValues' <- computeP $ coarsen boundValues
            v' <- fullMG approxOp residualOp (2*h) f' boundMask' boundValues'
         -- move to finer level     
            v <- computeP $ szipWith (+) boundValues  
                           $ szipWith (*) boundMask 
                           $ interpolate v'  
         -- solve for finer level
            multiGrid approxOp residualOp h f boundMask boundValues v

Note that in the base case we need an initial guess for the coarsest initial array. We simply use the boundValues array for this. The interpolation of v requires resetting the boundary values.

Improved inlining using modules

The previous (higher order) versions of multiGrid and fullMG which take approxOp and residualOp as arguments can by substantially improved by importing the operations in a module and specialising the definitions for the imported functions. Thus we define a module which imports a separate module (defining approxOp and residualOp) and which exports

> multiGrid :: Monad m =>
>              Double
>              -> Array U DIM2 Double
>              -> Array U DIM2 Double
>              -> Array U DIM2 Double
>              -> Array U DIM2 Double
>              -> m (Array U DIM2 Double)

> multiGrid !h !f !boundMask !boundValues !uInit 
>  = if coarsest uInit
>    then 
>     do computeP $ approxOp h f boundMask boundValues uInit
>    else
>     do v  <- iterateSolver (approxOp h f boundMask boundValues) steps1 uInit
>        r  <- computeP $ residualOp h f boundMask v
>        boundMask' <- computeP $ coarsen boundMask
>        r'  <- computeP $ szipWith (*) boundMask' $ restrict r 
>        let zeros = zeroArray (extent r')
>        err  <- multiGrid (2*h) r' boundMask' zeros zeros
>        vC   <- computeP $ szipWith (+) v  
>                            $ szipWith (*) boundMask 
>                            $ interpolate err
>        iterateSolver (approxOp h f boundMask boundValues) steps2 vC 

and the main full multigrid operation

> fullMG :: Monad m =>
>           Double
>           -> Array U DIM2 Double
>           -> Array U DIM2 Double
>           -> Array U DIM2 Double
>           -> m (Array U DIM2 Double)

> fullMG !h !f !boundMask !boundValues
>  = if coarsest boundValues
>      then do computeP $ approxOp h f boundMask boundValues boundValues
>      else do 
>             f'  <- computeP $ restrict f
>             boundMask' <- computeP $ coarsen boundMask
>             boundValues' <- computeP $ coarsen boundValues
>             v'  <- fullMG (2*h) f' boundMask' boundValues' 
>             v   <- computeP $ szipWith (*) boundValues  
>                                $ szipWith (*) boundMask 
>                                $ interpolate v'
>             multiGrid h f boundMask boundValues v

These versions perform significantly better when optimised.

Poisson’s Equation

For the two dimensional case, where , Poisson’s equation has the form

for subject to the Dirichlet boundary condition for

So the linear operator here is the Laplace operator

We need to make use of finite difference analysis to implement the approximation step and the residual calculation for Poisson’s equation.

Finite difference analysis with the central difference operator leads to the following (five point difference) formula approximating Poisson’s equation

Rewriting the above as

gives us an algorithm (approxOp) for the approximating iteration step. This is known as the the Jacobi method. We map a suitable stencil to add the four neighbours of each item, then we add corresponding elements of the array f multiplied by , then we divide by 4 and reset the boundary with the mask and values.

> {-# INLINE approxOp #-} 
> approxOp :: Double
>             -> Array U DIM2 Double
>             -> Array U DIM2 Double
>             -> Array U DIM2 Double
>             -> Array U DIM2 Double
>             -> Array (TR PC5) DIM2 Double               

> approxOp !h !f !boundMask !boundValues !arr
>        =   szipWith (+) boundValues
>          $ szipWith (*) boundMask
>          $ smap (/4)
>          $ szipWith (+) (R.map (*hSq) f )
>          $ mapStencil2 (BoundConst 0) 
>             [stencil2|   0 1 0
>                          1 0 1 
>                          0 1 0 |] arr
>          where hSq = h*h

We have chosen to leave the resulting array in a non-manifest form so that inlining can combine this operation with any subsequent operations to optimise the combination.

To calculate residuals from an estimation , we can use the same formula as above to approximate . Rearranging we have

This leads us to implement residualOp with a five point stencil that adds four neighbours and subtracts four times the middle item. After this we divide all items by and then add f and reset the boundary.

> {-# INLINE residualOp #-}
> residualOp :: Double
>               -> Array U DIM2 Double
>               -> Array U DIM2 Double
>               -> Array U DIM2 Double
>               -> Array (TR PC5) DIM2 Double

> residualOp !h !f !boundMask !v
>        =   szipWith (*) boundMask
>          $ szipWith (+) f
>          $ smap (*hFactor)
>          $ mapStencil2 (BoundConst 0) 
>             [stencil2|   0 1  0
>                          1 -4 1 
>                          0 1  0 |] v
>          where hFactor = 1/(h*h)

As noted earlier, the boundary is reset simply with the mask. The boundary values should all be zero for residuals, so we can drop any resetting of boundary values as the mask will have done that already.

The above two functions (approxOP and residualOp) are put in a module PoissonOps to be imported by the module MultiGrid_Poisson.

Setting up to run an example

This example is taken from (Iserles 2009 p.162). On the unit square

The Dirichlet boundary conditions are

and

We define the following to calculate values from the number of grids to be used.

> fineGridSpaces:: Int -> Int
> fineGridSpaces gridStackSize
>   = 2^gridStackSize

> fineGridShape:: Int -> DIM2
> fineGridShape gridStackSize 
>   = (Z :. n :. n) where n = 1+fineGridSpaces gridStackSize

We will set a global value for the sides of the area (the unit square) and also set a stack size of 9 grids so the finest will have by (i.e. 513 by 513) grid points

> distance :: Double
> distance = 1.0        -- length of sides for square area
> gridStackSize :: Int
> gridStackSize = 9     -- number of grids to be used including finest and coarsest

> intervals = fineGridSpaces gridStackSize
> shapeInit = fineGridShape gridStackSize

> hInit :: Double
> hInit = distance / fromIntegral intervals -- initial finest grid spacing

The initial arrays are defined using

> boundMask :: Array U DIM2 Double
> boundMask = 
>  fromListUnboxed shapeInit $ concat 
>   [edgeRow,
>    take ((intervals-1)*(intervals+1)) (cycle mainRow),
>    edgeRow
>   ]
>   where edgeRow = replicate (intervals+1) 0.0
>         mainRow = 0.0: replicate (intervals-1) 1.0  Prelude.++ [0.0]
  

> coordList :: [Double]
> coordList = Prelude.map ((hInit*) . fromIntegral) [0..intervals] 

> boundValues :: Array U DIM2 Double
> boundValues =
>   fromListUnboxed shapeInit $ concat
>      [Prelude.map (\j -> sin (pi*j)) coordList, 
>       concat $ Prelude.map mainRow $ tail $ init coordList,
>       Prelude.map (\j -> exp pi * sin (pi*j) + 0.5*j^2) coordList
>      ] 
>      where mainRow i = replicate intervals 0.0 Prelude.++ [0.5*i^2]

> fInit :: Array U DIM2 Double
> fInit =                        -- negative of RHS of Poisson Equation 
>  fromListUnboxed shapeInit $ concat $ Prelude.map row coordList    
>  where row i = Prelude.map (item i) coordList
>        item i j = -(i^2 + j^2)

> uInit :: Array U DIM2 Double
> uInit = boundValues

We are now in a position to calculate

> test1 = multiGrid hInit fInit boundMask boundValues uInit

and

> test2 = fullMG hInit fInit boundMask boundValues

as well as comparing with iterations of the basic approximation operation on its own

> solverTest :: Monad m => Int -> m(Array U DIM2 Double)
> solverTest n  = iterateSolver (approxOp hInit fInit boundMask boundValues) n uInit

We note that this particular example of Poisson’s equation does have a known exact solution

So we can calculate an array with values for the exact solution and look at the errors. We just look at the maximum error here.

> exact :: Array U DIM2 Double
>  exact =
>  fromListUnboxed shapeInit $
>  concat $ Prelude.map row coordList    
>  where row i = Prelude.map (item i) coordList
>        item i j = exp (pi*i) * sin (pi*j) + 0.5*i^2*j^2

> maxError :: Monad m => m(Array U DIM2 Double) -> m (Double)
> maxError test = 
>   do ans <- test
>      err :: Array U DIM2 Double 
>          <- computeP
>                $ R.map abs
>                $ R.zipWith (-) ans exact
>      foldAllS max 0.0 err

A note on errors

The Jacobi method used for the approximation step is known to converge very slowly, but fullMG significantly improves convergence rates. The literature analysing multigrid and full multigrid shows that the errors are reduced rapidly if the approximation method smooths the high frequency errors quickly (error components with a wavelength less than half the grid size). Unfortunately the Jacobi method does not do this. Gauss-Seidel and SOR methods do smooth high frequency errors but their usual implementation relies on in-place updating of arrays which is problematic for parallel evaluation. However, although multiGrid with the Jacobi approxOp does not reduce errors rapidly, fullMG does (with only gridStackSize=7 and step1=step2=3). This shows that the improved initial guess using fullMG outweighs the slow reduction of errors in multiGrid.

Some results

Clearly this code is designed to be run on multi-cores to use the parallelism. However, even using a single processor, we can see the impressive performance of Repa with optimisations on some simple tests of the Poisson example (running on a 2.13 GHz Intel Core 2 Duo).

So increasing the grid stack size by 1 roughly quadruples the size of the fine grid, increases runtime by a similar factor, and improves the error by a factor of 4 as well.

As a comparison, multiGrid alone (with the Jacobi approxOp) has high errors and increasing grid stack size does not help.

To see the problems with the basic Jacobi method on its own we note that after 400 iteration steps (with grid stack size of 6) the max error is 5.587133822631921. Increasing the grid stack size just makes this worse.

Compiler options

As we had observed in previous work Repa Laplace and SOR the large amount of inline optimisation seems to trigger a ghc compiler bug warning. This is resolved by adding the compiler option -fsimpl-tick-factor=1000 (along with -O2).

Some variations

We also tried two variations on opproxOp.

The first was to add in an extra relaxation step which weights the old value with the new value at each step using a weighting between 0 and 1. This slightly increases runtime, but does not improve errors for fullMG. It’s purpose is to improve smoothing of errors, but as we have seen although this may be significant for multiGrid it is not for fullMG. The improved starting points are much more significant for the latter. Note that successive over relaxing (SOR) with between 1 and 2 is not suitable for use with the Jacobi method in general as this is likely to diverge.

omega :: Double
omega = 2/3

{-# INLINE approxOp #-} 
approxOp :: Double
            -> Array U DIM2 Double
            -> Array U DIM2 Double
            -> Array U DIM2 Double
            -> Array U DIM2 Double
            -> Array (TR PC5) DIM2 Double               

approxOp !h !f !boundMask !boundValues !arr
       = szipWith weightedSum arr         -- extra relaxation step
         $ szipWith (+) boundValues
         $ szipWith (*) boundMask
         $ smap (/4)
         $ szipWith (+) (R.map (*hSq) f )
         $ mapStencil2 (BoundConst 0) 
            [stencil2|   0 1 0
                         1 0 1 
                         0 1 0 |] arr
         where hSq = h*h
               weightedSum old new = (1-omega) * old + omega * new

The second variation is based on the ‘modified 9 point scheme’ method discussed in (Iserles 2009 p.162). This has slower runtimes than the 5 point Jacobi scheme as it involves two stencil mappings (one across f as well as u). It had no noticable improvement on error reduction for fullMG.

-- modified 9 point scheme
approxOp ::  Double
                -> Array U DIM2 Double
                -> Array U DIM2 Double
                -> Array U DIM2 Double
                -> Array U DIM2 Double
                -> Array (TR PC5) DIM2 Double
            
approxOp !h !f !boundMask !boundValues !arr
       =   R.szipWith (+) boundValues
         $ R.szipWith (*) boundMask
         $ R.smap (/20)
         $ R.szipWith (+) 
            (R.map (*hFactor) 
             $ mapStencil2 (BoundConst 0) 
                  [stencil2|   0 1 0
                               1 8 1 
                               0 1 0 |] f )
         $ mapStencil2 (BoundConst 0) 
            [stencil2|   1 4 1
                         4 0 4 
                         1 4 1 |] arr
         where hFactor = h*h/2

This algorithm was obtained from rearranging the formula

Parallel Processing

The important point to note is that this code is ready to run using parallel processing. Results using multi-core processors will be posted in a separate blog.

Acknowledgements and References

Iserles (2009) gives some analysis of the multigrid method and this book also provided the example used (and the modified 9 point scheme). There are numerous other finite methods books covering multigrid and some devoted to it. Online there are some tutorial slides by Briggs part1 and part2 and some by Stüben amongst many others. A paper by Lippmeier, Keller, and Peyton Jones discussing Repa design and optimisation for stencil convolution in Haskell can be found here.

In a previous blog Repa Laplace and SOR, I used Repa to implement a Laplace solver using the Red-Black scheme. The explanation of alternating stencils probably needed a diagram, so here it is.

This diagram illustrates the shapes of the stencils for adding neighbours of red and black cells. It shows that two different stencils are needed (for odd and even rows) and these are swapped over for red and black.

Red-Black Neighbour Stencils

(The diagram was produced with Haskell Diagrams)

This describes the result of some experiments to enhance an existing elegant (haskell parallel array) laplace solver by using successive over-relaxation (SOR) techniques. The starting point was a laplace solver based on stencil convolution and using the Haskell Repa library (Regular Parallel Arrays) as reported in a paper by Lippmeier, Keller and Peyton-Jones. (There is also a 2011 published paper with the first two authors (Lippmeier and Keller 2011))

Background

Laplace’s Equation and Iterative solvers

The Laplace equation

is usually abbreviated to simply , where is a scalar potential (such as temperature) over a smooth surface.

For numerical solutions we will restrict attention to a rectangular surface with some fixed boundary values for u. (The boundary is not necessarily the border of the rectangle, as we can allow for fixed values at internal regions as well). The problem is made discrete for finite difference methods by imposing a grid of points over the surface and here we will assume this is spaced evenly ( in the direction and in the direction). Approximating solutions numerically amounts to approximating a fixedpoint for the grid values satisfying the boundary values and also, for non-boundary points (i,j) satisfying

where

If we also assume then this simplifies to

Iterating to find this fixed poiint from some suitable starting values is called relaxation. The simplest (Jacobi) method involves iterating the calculation of new values from previous values until the iterations are close enough to converging.

This method, whilst very simple, can be extremely slow to converge. One improvement is given by the Gauss-Seidel method which assumes the next set of values in each iteration are calculated row by row in a scan across the columns allowing some new values to be used earlier during the same iteration step (changing to for calculated values in the previous column and in the previouus row).

Unfortunately, this assumption about the order of evaluation of array elements interferes with opportunities to use parallel computation of the array elements (as we wish to do using Repa array primitives).

Successive over-relaxation (SOR)

Successive over-relaxation is a method used to speed up convergence of the iterations by using a weighted sum of the previous iteration with the new one at each relaxation step. Using as the weight we calculate first (substituting for in the above equation), then calculate

Values of less than 1 will slow down the convergence (under-relaxation). A value of between 1 and 2 should speed up convergence. In fact the optimal value for will vary for each iteration, but we will work with constant values only here. It is also important to note that starting with values above 1 will generally cause oscillations (divergence) if just the Jacobi scheme is used, so successive over-relaxation relies on a speed up of the basic iteration such as that provided by the Gauss-Seidel scheme.

The original stencil version using Repa

In the paper a stencil mapping technique is used to get fast performance using the Repa array primitives. The code uses the Repa library version 3

The main relaxation step is expressed as

relaxLaplace arr
 = computeP
 $ szipWith (+) arrBoundValue   -- boundary setting
 $ szipWith (*) arrBoundMask    -- boundary clearing
 $ smap (/ 4)
 $ mapStencil2 (BoundConst 0)
    [stencil2|   0 1 0
                 1 0 1 
                 0 1 0 |] arr

We will be taking this as a starting point, so we review details of this (working from the bottom up).

A stencil is described in a quoted form to express which neighbours are to be added (using zeros and ones) with the element being scanned assumed to be the middle item. Thus this stencil describes adding the north, west, east, and south neighbours for each item. The mapStencil2 is a Repa array primitive which applies the stencil to all elements of the previous iteration array (arr). It is also supplied with a border directive (Boundconst 0) which just says that any stencil item arising from indices outside the array is to be treated as 0. This stencil mapping will result in a non-manifest array, in fact, in a partitioned array split up into border parts and non-border parts to aid later parallel computation. After the stencil map there is a mapping across all the results to divide by 4 using smap (which improves over the basic repa array map for partitioned arrays). These calculations implement the Jacobi scheme.

The subsequent two szipWith operations are used to reset the fixed boundary values. (The primitive szipWith again improves over the basic repa array zipWith by catering explicitly for partitioned arrays.) More specifically, szipWith(*) is used with an array (arrBoundMask) which has zero for boundary positions and 1 for other positions – thus setting the boundary items to 0. After this szipWith(+) is used with an array (arrBoundValues) which has the fixed initial values for the boundary items and zero for all other items. Thus the addition reinstates the initial boundary values.

These array calculations are all delayed so a final computeP is used to force the parallel evaluations of items to produce a manifest array result. This technique allows a single pass with inlining and fusion optimisations of the code to be possible for the calculation of each element. Use of computeP requires the whole calculation to be part of a monad computation. This is a type constraint for the primitive to exclude the possibility of nested calls of computeP overlapping.

A significant advantage of this stencil implementation is that programming directly with array indices is avoided (these are built into the primitives once and for all and are implemented to avoid array bound checks being repeated unnecessarily).

Unfortunately, though, this implementation is using the Jacobi scheme for the iteration steps which is known to have very slow convergence.

We would like to improve on this by using successive over-relaxation, but as pointed out earlier, this will not work with the Jacobi scheme. Furthermore, the Gauss-Seidel improvement will not help because it is based on retrieving values from an array still being updated. This is not compatible with functional array primitives and stencils and prevents simple exploitation of parallel array calculations.

To the rescue – the red-black scheme for calculating iterations.

Red-black Iterations

This scheme is well known and based on the observation that on a red-black checkerboard, for any red square, the north, west, east, and south neighbours are all black and conversely neighbours of black squares are all red. This leads to the simple idea of calculating updates to all the red squares first, then using these updated values to calculate the new black square values.

Implementation using stencils and repa

The set up

We split the original array into a red array and black array with simple traverse operations. We assume the original array (arr) has an even number of columns so the number of red and black cells are the same. As a convention we will take the (0,0) item to be red so we want, for red array (r)

r(i,j) = arr(i, 2*j + (i `mod` 2))

where

arr has i <- 0..n-1, j <- 0..2m-1
r   has i <- 0..n-1, j <- 0..m-1

The traverse operation from the Repa library takes as arguments, the original array, a mapping to express the change in shape, and a lookup function to calculate values in the new array (when given a get operation for the original array and a coordinate (i,j) in the new array)

> projectRed :: Array U DIM2 Double -> Array D DIM2 Double
> projectRed arr =  
>   traverse arr 
>           (\ (e :. i :. j) -> (e :. i :. (j `div` 2)))
>           (\get (e :. i :. j) -> get (e :. i :. 2*j + (i `mod` 2)))

Here (and throughout) we have restricted the type to work with two dimensional arrays of Double, although a more general type is possible. Notice also that the argument array is assumed to be manifest (U) and the result is a delayed array (D).

Similarly for the black array (b) we want

b(i,j) = arr(i, 2*j + ((i+1) `mod` 2))

with the same extent (shape) as for r, hence

> projectBlack :: Array U DIM2 Double -> Array D DIM2 Double
> projectBlack arr =  
>     traverse arr 
>             (\ (e :. i :. j) -> (e :. i :. (j `div` 2)))
>             (\ get (e :. i :. j) -> get(e :. i :. 2*j + ((i+1) `mod` 2)))

We can also use these same functions to set up boundary mask and value arrays separately for red and black from the starting array (arrInit), initial mask (arrBoundMask), and boundary values (arrBoundValue)

do redBoundMask    <- computeP $ projectRed arrBoundMask
   blackBoundMask  <- computeP $ projectBlack arrBoundMask
   redBoundValue   <- computeP $ projectRed arrBoundValue
   blackBoundValue <- computeP $ projectBlack arrBoundValue
   redInit         <- computeP $ projectRed arrInit
   blackInit       <- computeP $ projectBlack arrInit

This is part of a monad computation with each step using a computeP (parallel array computation) to create manifest versions of each of the arrays we need before beginning the iterations. These calculations are independent, so the sequencing is arbitrary.

The finish

At the end of the iterations we will reverse the split into red and black by combining using the traverse2 operation from the Repa library. We need

arr(i,j) = r(i, j `div` 2) when even(i+j)
         = b(i, j `div` 2) otherwise            

where

arr has i <- 0..n-1 , j <- 0..2m-1
r   has i <- 0..n-1 , j <- 0..m-1
b   has i <- 0..n-1 , j <- 0..m-1

The traverse2 operation takes two arrays (here – of the same extent), a mapping to express the new extent (when given the two extents of the argument arrays), and a function to calculate values in the new array (when given the respective get operations for the original arrays (here – get1 and get2) and a coordinate (i,j) in the new array).

> combineRB r b =   
>   traverse2 r b
>            (\ (e :. i :. j) _ -> (e :. i :. 2*j))
>            (\ get1 get2 (e :. i :. j) -> 
>               (if even(i+j) then get1 else get2) (e :. i :. j `div` 2)
>            )
 

Stencils in the iteration step

We use stencils as in the original, but when we consider a stencil on black cells which are to be combined as neighbours of the red cell r(i,j) we need the following shape, where b(i,j) corresponds to the middle item

                    0 1 0         
                    1 1 0
                    0 1 0

BUT this is only when processing EVEN rows from the black array. For odd rows we need a different shape

                    0 1 0 
                    0 1 1 
                    0 1 0

That is, we need to apply one of two different stencils depending on whether the row is even or odd. Similarly in processing the red array to combine red neighbours of b(i,j) we need the same shaped stencils but the first shape above is used on odd rows of red cells and the second shape on even rows. We define and name the stencils as leftSt and rightSt

> leftSt :: Stencil DIM2 Double
> leftSt  =   [stencil2|  0 1 0         
>                         1 1 0 
>                         0 1 0 |]

> rightSt :: Stencil DIM2 Double
> rightSt =   [stencil2|  0 1 0 
>                         0 1 1 
>                         0 1 0 |]

Critically, we will need an efficient way to apply alternate stencils as we map across an array. At first, it may seem that we might have to rebuild a version of the primitive mapStencil2 to accommodate this, but that will get us into some complex array representation handling which is built into that primitive. On reflection, though, the combination of lazy evaluation and smart inlining/fusion optimisation by the compiler should allow us to simply apply both stencils on all elements and then choose the results we actually want. The laziness should ensure that the unwanted stencil applications are not actually calculated, and the administration of choosing should be simplified by compiler optimisations.

To apply our stencils we will use the Repa array primitive mapStencil2 which takes as arguments, a border handling directive, a stencil, and a (manifest) array, to produce a resulting (delayed) array.

 mapStencil2 :: Boundary Double
             -> Stencil DIM2 Double
             -> Array U DIM2 Double
             -> Array D DIM2 Double 

The choice of stencil will depend on the position (actually the evenness of the row). As we cannot refer to the indices when using the stencil mapping operation we are led to using traverse2 after mapping both stencils to select results from the two arrays produced. Our alternate stencil mapping operation will have a similar type to mapStencil2 except that it expects two stencils rather than one.

altMapStencil2 :: Boundary Double
               -> Stencil DIM2 Double
               -> Stencil DIM2 Double
               -> Array U DIM2 Double
               -> Array D DIM2 Double 

altMapStencil2 !bd !s1 !s2 !arr
= traverse2 (mapStencil2 bd s1 arr) (mapStencil2 bd s2 arr)
            (\ e _ -> e)
            (\ get1 get2 (e :. i :. j) -> 
                   (if even i then get1 else get2) (e :. i :. j)
            )

This function needs to be inlined (with a pragma) and the bang annotations on arguments are there to improve optimisation opportunities.

The iteration step

The main relaxLaplace iteration step involves the following (where r and b are the previous red and black arrays)

• Calculating a new version of r using stencils over b

  r1 = smap (/4)
       $ altMapStencil2 (BoundConst 0) leftSt rightSt b

• Over-relaxation of this taking weighted sums (using r and r1)

  r2 = szipWith weightedSum r r1
       where weightedSum old new = (1-omega)*old + omega*new

• Reset the boundaries to get the new version of r (r’)

  r' = szipWith (+) redBoundValue      -- boundary resetting
       $ szipWith (*) redBoundMask r2  -- boundary clearing

• Do similar steps for calculating b’ but starting with stencils over r’ (not r) and switching the left and right stencils

This is combined in the following monad computation (where r and b are the old arrays). The monad is necessary because we want to use computeP to ensure that the final returned arrays are manifest.

 do r' <- computeP
          $ relaxStep r b redBoundValue redBoundMask leftSt rightSt
    b' <- computeP 
          $ relaxStep b r' blackBoundValue blackBoundMask rightSt leftSt
    ...

which uses

relaxStep !arrOld !arrNbs !boundValue !boundMask !stencil1 !stencil2 
   = szipWith (+) boundValue
     $ szipWith (*) boundMask
     $ szipWith weightedSum arrOld
     $ smap (/4)
     $ altMapStencil2 (BoundConst 0) stencil1 stencil2 arrNbs

weightedSum !old !new = (1-omega)*old + omega*new

The first argument for relaxStep is the old (red or black) array we are calculating an update for, and the second argument is the neighbour array to use stencils on. The old array is needed for the over-relaxation step with weightedSum.

It is also worth pointing out that we have the over-relaxation step being done before the two boundary resetting steps, but this could just as easily be done after the boundary resetting as it does not make a difference to the final array produced.

Laplace with Red-Black and Stencils

Finally the main function (solveLaplace) sets up the initial arrays and passes them to the function with the main loop (iterateLaplace)

> solveLaplace::
>   Monad m
>   => Int		   -- ^ Number of iterations to use.
>   -> Double		   -- ^ weight for over relaxing (>0.0 and <2.0)
>   -> Array U DIM2 Double -- ^ Boundary value mask.
>   -> Array U DIM2 Double -- ^ Boundary values.
>   -> Array U DIM2 Double -- ^ Initial state. Should have even number of columns
>   -> m (Array U DIM2 Double)

> solveLaplace !steps !omega !arrBoundMask !arrBoundValue !arrInit
> do redBoundMask    <- computeP $ projectRed arrBoundMask
>    blackBoundMask  <- computeP $ projectBlack arrBoundMask
>    redBoundValue   <- computeP $ projectRed arrBoundValue
>    blackBoundValue <- computeP $ projectBlack arrBoundValue
>    redInit         <- computeP $ projectRed arrInit
>    blackInit       <- computeP $ projectBlack arrInit
>    iterateLaplace steps omega redInit blackInit
>                   redBoundValue blackBoundValue redBoundMask blackBoundMask

where

iterateLaplace !steps !omega !redInit !blackInit
               !redBoundValue !blackBoundValue !redBoundMask !blackBoundMask 
     = go steps redInit blackInit
       where 
         go 0 !r !b = computeP $ combineRB r b -- return final combined array
         go n !r !b 
            = do r' <- computeP
                       $ relaxStep r b redBoundValue redBoundMask leftSt rightSt
                 b' <- computeP 
                       $ relaxStep b r' blackBoundValue blackBoundMask rightSt leftSt
                 go (n - 1) r' b'

         {-# INLINE relaxStep #-}
         relaxStep !arrOld !arrNbs !boundValue !boundMask !stencil1 !stencil2 
            = szipWith (+) boundValue
              $ szipWith (*) boundMask
              $ szipWith weightedSum arrOld
              $ smap (/4)
              $ altMapStencil2 (BoundConst 0) stencil1 stencil2 arrNbs

         {-# INLINE weightedSum #-}
         weightedSum !old !new = (1-omega)*old + omega*new

{-# INLINE iterateLaplace #-}

Performance and a New Version

The number of calculations in an iteration of red-black is comparable to the original stencil implementation (with just the weighting operations addded in). Although there are two arrays to update, they are half the size of the original. We would expect optimised performance to be only fractionally slower than the original. The speedup in progress towards convergence can be dramatic (one iteration per 8 of the original for our test examples). So, in principle, this would be a big improvement. Unfortunately optimisation did not seem to be achieving the same speedups as the original and the code was roughly 12 times slower.

After much experimentation it looked as though the inner traverse2 operation of altMapStencil2 was inhibiting the optimisations (fusions with stencil mapping code and subsequent maps and zips).

A better performance was achieved by separating the stencil mapping from the traverse2 for alternate row selection and delaying the alternate row selection until after all the other operations. The new version drops altMapStencil2 and instead simply uses altRows

> altRows :: forall r1 r2 a . (Source r1 a, Source r2 a)
>         => Array r1 DIM2 a -> Array r2 DIM2 a -> Array D DIM2 a  
>             
> altRows !arr1 !arr2 =  -- assumes argument arrays with the same shape
>     traverse2 arr1 arr2
>               (\ e _ -> e)
>               (\ get1 get2 e@(_ :. i :. _) -> 
>                   if even i then get1 e else get2 e
>               )
>                  
> {-# INLINE altRows #-}

The function relaxStep in the following revised version of iterateLaplace has altRows done last, with mapStencil2 applied first (using a different stencil in each alternative)

> iterateLaplace ::
>   Monad m
>   => Int
>   -> Double
>   -> Array U DIM2 Double	
> 	-> Array U DIM2 Double
> 	-> Array U DIM2 Double	
> 	-> Array U DIM2 Double
> 	-> Array U DIM2 Double
> 	-> Array U DIM2 Double
> 	-> m (Array U DIM2 Double)
>           
> iterateLaplace !steps !omega !redInit !blackInit
>                !redBoundValue !blackBoundValue !redBoundMask !blackBoundMask 
>      = go steps redInit blackInit
>        where 
>          go 0 !r !b = computeP $ combineRB r b -- return final combined array
>          go n !r !b 
>             = do r' <- computeP 
>                        $ relaxStep r b redBoundValue redBoundMask leftSt rightSt
>                  b' <- computeP
>                        $ relaxStep b r' blackBoundValue blackBoundMask rightSt leftSt
>                  go (n - 1) r' b'
>              
>          {-# INLINE relaxStep #-}
>          relaxStep !arrOld !arrNbs !boundValue !boundMask !stencil1 !stencil2 
>             = altRows (f stencil1) (f stencil2)
>               where 
>                     {-# INLINE f #-}
>                     f s = szipWith (+) boundValue
>                           $ szipWith (*) boundMask
>                           $ szipWith weightedSum arrOld
>                           $ smap (/4)
>                           $ mapStencil2 (BoundConst 0) s arrNbs
>
>          {-# INLINE weightedSum #-}
>          weightedSum !old !new = (1-omega)*old + omega*new
>
> {-# INLINE iterateLaplace #-}

This now seems to be only a little slower to run than the original stencil solution (about 4 times slower so about 3 times faster than the previous version of red-black). Thus, with an approximately 8-fold speed up in convergence, this does give an overall improvement.

In order to compile this version, it was necessary to use not just the ghc -O2 flag, but also -fsimpl-tick-factor=1000. The need for this flag was indicated by a ghc compiler bug message.

All the code above with birdfeet ( >) in this (incomplete) literate haskell document is essentially that in the module RedBlackStencilOpt.hs (which has some extra preliminaries and imports and checking of arguments in functions which were elided here for clarity). This code can be found here along with the module RedBlackStencil.hs which contains the first version. These can both be loaded by the wrapper newWrapper.hs which is just an adaptation of the original wrapper to allow for passing the extra parameter.

Acknowledgements and References

There are many numerical methods books covering relevant background mathematics, but online, there are also lecture notes. In particular, on Computational Numerical Analysis of Partial Differential Equations by J.M.McDonough and on Numerical Solution of Laplace Equation by G.E.Urroz. T. Kadin has some tutorial notes with a diagram for the red-black scheme (and an implementation using Fortran 90). There is a useful Repa tutorial online and more examples on parallel array fusion are reported in (Lippmeier et al. 2012).