Hi Victor, Thanks for your interest!

We are slowly but surely moving ahead with this project, and FindExactCover should go through WFR pretty soon. The final write-up will be more interesting than these basically public-facing notes I have here, because we're starting to understand why timing statistics look the way they do. Unfortunately, we still haven't succeeded in compiling Knuth's Dancing Links X, which is the main outstanding "To do" item on this mini-project.

That's nice that you figured out how to go faster making the cover matrices. I also improved this code while I was running more tests on the actual solvers (where the bottleneck really is):

PieceRotations[True, True][piece_] := Union[
  ResourceFunction["ArrayRotations"][piece]]

PieceRotations[True, False][piece_] := Union[
  ResourceFunction["ArrayRotations"][piece, False]]

PieceRotations[_, _][piece_] := {piece}

PieceTranslates[piece_, shift_, dims_
  ] := With[{inset = PadRight[PadRight[#,
        dims[[2]], 0] & /@ piece,
     dims[[1]], {ConstantArray[0, dims[[2]] ]}]},
  Catenate[Outer[RotateRight[inset, {#1, #2}] &,
    Range[0, shift[[1]] ], Range[0, shift[[2]]], 1]]]

PartialCoverRows[tfRot_, tfRef_
   ][template_, piece_] := Module[{
   tot = Total[piece, 2],
   zeroes = Position[template, 0],
   rot = PieceRotations[tfRot, tfRef][piece],
   dims = Dimensions[template], shift, candidates},
  shift = Map[dims - Dimensions[#] &, rot];
  candidates = Catenate[MapThread[
     If[MemberQ[Sign[#2], -1],
       Nothing, PieceTranslates[#1, #2, dims]] &,
     {rot, shift}]];
  Select[Outer[#1[[Sequence @@ #2]] &,
    candidates, zeroes, 1],
   Total[#] == tot &]]

MatricesToCoverMatrix[template_, pieces_,
  OptionsPattern[{
    "Rotations" -> True,
    "Reflections" -> True,
    "PlaceOne" -> False
    }]] := If[TrueQ[OptionValue["PlaceOne"]],
  Catenate[MapThread[Function[{piece, ind},
     Join[ind, #] & /@ PartialCoverRows[
        OptionValue["Rotations"], OptionValue["Reflections"]
        ][template, piece]],
    {pieces, IdentityMatrix[Length[pieces]]}]],
  Catenate[PartialCoverRows[
       OptionValue["Rotations"], OptionValue["Reflections"]
       ][template, #] & /@ pieces]] 

Would you want to see this as a WFR? If yes, can you think of a better name for it?

It should be about as fast as what you have, but it's also keyed in with some options and better free space recognition. The design is so that it's easy to make the cover matrix for "Scott's pentomino problem" from Knuth's arxiv article:

pentominoProblem = MatricesToCoverMatrix[
   ArrayPad[ConstantArray[1, {2, 2}], 3, 0],
   upToPentominos[[5]], "PlaceOne" -> True];

Dimensions@pentominoProblem

Out[]= {1568, 72} 

We can find all solutions of this matrix, but our times are still a little slow compared to industry standard (due I think to not having compiled code). Assuming you can make this matrix, what's your fastest time for solving it with Mathematica code?

And now for a joke:

I can't understand what's taking so long. I thought this problem was so easy it would already be solved by now. Inefficiency!

Just to put things into perspective, here's a problem that apparently haunted me for over ten years (admittedly, while I was working on some other things in physics too) : :

Problem Statement. 6. Find all possible tilings ways to fit four copies of the following tile in a 7 by 7 grid, without any holes. Make a graphics function to display the results.

And, again copying from a 10+ year old notebook, the polyomino tile is plotted as follows:

ArrayPlot[{{1, 0, 1, 1}, {1, 1, 1, 0}, {0, 0, 1, 1}}, Mesh -> True]

Graphics[{Green, 
  Polygon[{{0, 1}, {2, 1}, {2, 0}, {4, 0}, {4, 1}, {3, 1}, {3, 2}, {4,
      2}, {4, 3}, {2, 3}, {2, 2}, {1, 2}, {1, 3}, {0, 3}}], Red, 
  Line[{{0, 1}, {2, 1}, {2, 0}, {4, 0}, {4, 1}, {3, 1}, {3, 2}, {4, 
     2}, {4, 3}, {2, 3}, {2, 2}, {1, 2}, {1, 3}, {0, 3}, {0, 1}}]}]

This problem is very similar to the one mentioned above, since it involves placing four tiles. However, it is not an exact cover we're looking for, only a partial. Since partial tilings can have holes, an additional constraint is added that results should be sorted by genus. For definiteness sake, let us state that the tile transforms by rotations and reflections and that a hole is any contiguous empty region not incident on the boundary of the 7x7 grid.

Now rather than getting to the 7x7 case right away, let's just "tool up" and run a series of smallish test case to make sure that the tools are working correctly. The 5x5 case is easy enough to work out by hand, essentially it has only two solutions, one of genus zero and multiplicity eight, and another with genus 1 and multiplicity 4:

Grid[Partition[ArrayPlot[#, ImageSize -> 80,
     Frame -> None,
     ColorRules -> {
       1 -> Blend[{Yellow, Orange}],
       2 -> Lighter[Blend[{Magenta, Red}, .25],
         0.4], 0 -> Darker@Gray}] & /@ Join[
    DeleteDuplicatesBy[
     Union[ResourceFunction["ArrayRotations"][{
        {1, 0, 1, 1, 0},
        {1, 1, 1, 0, 0},
        {2, 2, 1, 1, 0},
        {0, 2, 2, 2, 0},
        {2, 2, 0, 2, 0}
        }]], Sort@{#, # /. {1 -> 2, 2 -> 1}} &],
    DeleteDuplicatesBy[ResourceFunction["ArrayRotations"][{
       {0, 1, 0, 1, 1},
       {0, 1, 1, 1, 0},
       {2, 2, 0, 1, 1},
       {0, 2, 2, 2, 0},
       {2, 2, 0, 2, 0}
       }], Sort@{#, # /. {1 -> 2, 2 -> 1}} &]]
  , 4], Spacings -> {1, 1}]

Now let's try to get the same result, without just typing numbers into a matrix as are seen to fit. First, we transform the matrix polyomino into a mesh, then into a graph:

PolyominoGraph[arrayMesh_] := Construct[With[{edges = ReplaceAll[
       MeshPrimitives[#, 1], 
       Line[x_] :> UndirectedEdge @@ Round[x]]},
    Graph[edges, VertexCoordinates -> (# -> # & /@ Union[
         edges /. UndirectedEdge -> Sequence])]
    ] &, arrayMesh]

form = PolyominoGraph[
  ArrayMesh[{{1, 0, 1, 1}, {1, 1, 1, 0}, {0, 0, 1, 1}}]]

And we can do the same with the grid space where this tile and one other copy will be placed:

emptiness = PolyominoGraph[ArrayMesh[ConstantArray[1, {5,5}]]]

Now for some magic (using 13.1 because we noticed a bug!), sorry you'll have to figure the nuts and bolts of these next two functions on your own:

PolyominoCoverMatrix[forms_, emptiness_] := With[{
   form4Cycles = Map[Sort, FindCycle[#, {4}, All] & /@ Flatten[
       FindIsomorphicSubgraph[emptiness, #, All] & /@ forms], 3],
   emptiness4Cycles = Map[Sort, FindCycle[emptiness, {4}, All], 2]},
  MapThread[Join, {IdentityMatrix[Length@form4Cycles],
    Function[{cycles}, ReplacePart[Array[0 &, Length@emptiness4Cycles],
       # -> 1 & /@ 
        Flatten[Position[emptiness4Cycles, #] & /@ cycles]]
      ] /@ form4Cycles}]]

CoverGraph[coverMatrix_] := ResourceFunction["NestWhileGraph"][
  Function[{state}, 
   Select[state + # & /@ coverMatrix, ! MemberQ[#, 2] &]],
  {ConstantArray[0, Last[Dimensions[coverMatrix]]]}, UnsameQ @@ # &,  2]

Using FindIsomorphicSubgraph (new in V.13 and bug-fixed in V.13.1), and a new WFR from the games post, NestWhileGraph, we can quickly dispatch the test case, and finally solve the difficult problem from 10 years ago. First a cover matrix is constructed, then NestWhileGraph comprehensively sums rows to find non-overlapping placements:

CoverGraph[PolyominoCoverMatrix[{form},emptiness]]

And indeed if we look closely at the 12 vertices with graph distance 2 from the origin, undoubtedly we would again find the 12 relatively obvious solutions listed above. We can ramp up to a 6x6 grid and still find, fairly instantaneously, the following graph:

CoverGraph[PolyominoCoverMatrix[{form},
  PolyominoGraph[ArrayMesh[ConstantArray[1, {6, 6}]]]]]

And counting vertices by numbers of tiles placed, we find 40 configurations with 3 tiles on a 6x6 grid:

Histogram[ Total /@ VertexList[gTest][[All, 1 ;; First[Dimensions[mTest]]]]]

Now when we go to the 7x7 case, the computation takes more noticeable time, and the graph turns out to be very large, too large in fact, to show up nicely as a plot:

AbsoluteTiming[ gPoly = CoverGraph[PolyominoCoverMatrix[{form},
    PolyominoGraph[ArrayMesh[ConstantArray[1, {7, 7}]]]]]]

With a little more analysis, we can immediately solve the counting problem, first sorting the graph vertices by number of tiles, and next by genus:

ord = Plus[Divide[Mean[#], 2], 1/2] & /@ Map[Sort,
     FindCycle[emptiness, {4}, All], 2] /. UndirectedEdge -> Plus;

leveldVertices = Function[{count},
    Select[VertexList[gPoly2],
     Total[#[[1 ;; First[Dimensions[mSolve]]]]] == count &]
    ] /@ Range[0, 4];

leveledMatrices = Map[
   ReplacePart[ConstantArray[0, {7, 7}],
     MapThread[If[#2 == 1, #1 -> 1, #1 -> 0] &,
      {ord, #[[First[Dimensions[mSolve]] + 1 ;; -1]]}]] &,
   leveldVertices, {2}];

Genus[mat_] := Length[Select[WeaklyConnectedComponents[
    NearestNeighborGraph[Position[mat, 0], {All, 1}]],
   ! MemberQ[Flatten[#], 1 | 7] &]]

data = KeySort[CountsBy[#, Genus]] & /@ leveledMatrices

Out[]=<|0 -> 1|>, <|0 -> 160|>, <|0 -> 2628, 1 -> 828, 
 2 -> 232|>, <|0 -> 3384, 1 -> 4640, 2 -> 3072, 3 -> 1008, 
 4 -> 240|>, <|0 -> 96, 1 -> 322, 2 -> 560, 3 -> 632, 4 -> 392, 
 5 -> 142, 6 -> 40, 7 -> 4|>

BarChart[Values /@ data, ChartLayout -> "Stacked"]

It looks like there are 96 placements of 4 tiles on 7x7 grid with genus 0. Now let's plot some of these arrangements to see what they look like, and the four genus 7 cases might also be interesting. First we need to transform the data a bit, and a depict function:

tileVals = MapThread[If[#2 == 1, #1 -> x, #1 -> 0] &, {ord, #}] & /@ 
   mSolve[[All, First[Dimensions[mSolve]] + 1 ;; -1]];

coloredCovers = Total[
     Map[ReplacePart[ConstantArray[0, {7, 7}], #] &, 
      MapIndexed[# /. x -> #2[[1]] &, 
       tileVals[[Position[#, 1][[All, 1]]]]]]] & /@ 
   Select[VertexList[gPoly2], 
     Total[#[[1 ;; First[Dimensions[mSolve]]]] ] == 4 &][[All, 
     1 ;; First[Dimensions[mSolve]]]];

CoverPlot[cover_] := ArrayPlot[cover,
  ColorRules -> 
   Append[MapIndexed[#2[[1]] -> #1 &, {Blend[{Yellow, Orange}],
      Blend[{Green, Cyan}, .5], Lighter[Blue, .5],
      Lighter[Blend[{Magenta, Red}, .25], 0.4]}], 0 -> Darker@Gray],
  Frame -> None, Mesh -> True, MeshStyle -> Black]

A random sample, labeled by genus:

Grid[{Labeled[Show[CoverPlot@#, ImageSize -> 80], 
     Text@Style[Genus[#], Gray]] & /@ 
   RandomChoice[coloredCovers, 6]},
 Spacings -> {1, 1}]

The 96 results for genus 0:

Grid[Partition[Labeled[Show[CoverPlot@#, ImageSize -> 80],
     Text@Style[Genus[#], Gray]] & /@ 
   Select[coloredCovers, Genus[#] == 0 &], 6],
 Spacings -> {1, 1}]

And four more results for genus 7:

Grid[{Labeled[Show[CoverPlot@#, ImageSize -> 80],
     Text@Style[Genus[#], Gray]] & /@ 
   Select[coloredCovers, Genus[#] == 7 &]},
 Spacings -> {1, 1}]

And please notice obvious dihedral symmetry in this set of results!

As a double check, I can run a completely different algorithm based on MultiwayDeletionsGraph for which the input is not a cover matrix, rather is a cover graph derived from the cover matrix:

PairwiseOverlapGraph[coverMatrix_] := 
 With[{len = Length@coverMatrix},
  Graph[Range[len], Flatten[If[
       ! MemberQ[Total[coverMatrix[[#]]], 2],
       UndirectedEdge @@ #, {}] & /@ Subsets[Range[len], {2}]]]]

PairwiseOverlapGraph[PolyominoCoverMatrix[{form},
 PolyominoGraph[ArrayMesh[ConstantArray[1, {7, 7}]]]]]

The augmented MultiwayDeletionsGraph with method "ExactCover" then crawls this graph and can produce an isomorphic solution graph output, but in less than 2 seconds! A factor of 10 faster! However, this is not really scientific reproducibility because in both cases the machine operator is the same person, namely, me. If another scientist wants to verify or dispute the solution above, please add a comment to this thread! Also, if anyone can beat the 2s time, please say how and how fast, thanks.

The last thing I would like to mention is that Knuth's DLX algorithm can probably outperform MultiwayDeletionsGraph by an order of magnitude. What I really hope to do (instead of just copying someone else's implementation) is formulate a graph input similar to what is directly above, that can then be fed into a new method for MultiwayDeletionsGraph. I think this is possible, but I haven't had the time to work out all the details just yet. Anyone?