Listing matrices up to symmetry

Question

I am interested in the equivalence relation on N x N binary matrices, in which two matrices are equivalent if one can be obtained by rotating/reflecting the other. I would like to obtain a list

binclasses = { {B1, m1},...{Bk, mk} },

where B_i is a representative element of class i, and m_i is the size of that class. (Of course the requirement of binary matrices is somewhat artificial here, so the question could be rephrased for matrices taking values in a finite set). For N =2, I am looking to obtain the list

binclasses = {{{{0, 0}, {0, 0}}, 1}, {{{0, 0}, {0, 1}}, 4}, {{{0, 0}, {1, 1}}, 4}, {{{0, 1}, {1, 0}}, 2}, {{{0, 1}, {1, 1}}, 4}, {{{1, 1}, {1, 1}}, 1}}

My thought process so far has been as follows:

1. Define a list binlist of all N x N binary matrices.
2. Choose an element currmat from binlist, and generate a list currsym of all rotations and reflections of currmat.
3. Define a new list binminsym which is binlist with all elements of currsym removed.
4. Add the pair {currmat, Length[currsym]} to a new list binclasses.
5. Set binlist = binminsym, return to step 1.

I am guessing that this is not a particularly efficient approach, and am curious to know how one could better go about the problem.

EDIT A second update, I have improved on the original code, however, it still runs slowly (original edit left below this for information). In this variation I first group binary matrices by the number of 1s, and the number of adjacent 1s; this means that fewer list elements need be checked against each other. At the moment for the 4 x 4 case, this version takes a little over 4 minutes to identify the 8548 equivalence classes, compared to the 17 minutes of the previous method.

(* sTotal function gives total of all matrix entries *)
sTotal[matrix_] := 
  Sum[matrix[[i, j]], {i, 1, Length[matrix]}, {j, 1, Length[matrix]}];

(* eTotal function gives total number of adjacent non-zero entries *)

halo[matrix_] :=
 Module[{cmat},
  cmat = ConstantArray[0, {Length[matrix] + 1, Length[matrix] + 1}];
  cmat[[1 ;; Length[matrix], 1 ;; Length[matrix] ]] = matrix;
  cmat]

eTotal[matrix_] := 
  Sum[halo[matrix][[i, j]] (halo[matrix][[i, j + 1]] + 
      halo[matrix][[i + 1, j]]), {i, 1, Length[matrix]}, {j, 1, 
    Length[matrix]}];


Timing[
 Clear[mdim, blist, bgrlist, binclasses, binreps, currbg, currmat, 
   currsym]

  mdim = 4; (*matrix dimensions*)
 blist = Tuples[{0, 1}, {mdim, 
    mdim}];  (* All square binary matrices of size mdim *) 

 (* blist grouped into sets of matrices with equal number of 1s, and \
equal adjacent 1s *)
 bgrlist = GatherBy[blist, {sTotal, eTotal}]; 

 (* Initialise list of equivalence classes *)
 binclasses = {};
 binreps = {};
 ]


Timing[
 Do[
  Do[
    currbg = bgrlist[[i, j]];
    While[Length[currbg] >= 1,
     currmat = currbg[[1]];
     currsym = symlist[currmat];
     binclasses = Append[binclasses, {currmat, Length[currsym]}];
     binreps = Append[binreps, currmat];
     If[i != 1 + (mdim^2 /2 ),
      binreps = 
        Append[binreps, ConstantArray[1, {mdim, mdim}] - currmat];
      ];
     currbg = DeleteCases[currbg, x_ /; MemberQ[currsym, x] == True];
     ]
    , {j, 1, Length[bgrlist[[i]]]}];
  , {i, 1, Ceiling[Length[bgrlist]/2]}];
 Length[binreps]
 ]

EDIT Here is my code for the above algorithm. Note, that I have defined the symmetries of the square as functions, even though some of them (eg. reflection in diagonal) are simply transposing the matrix. I did this for clarity (for readers on here), rather than efficiency!

(*Symmetry functions.*)
vsym[matrix_] := Reverse[matrix, {2}];
hsym[matrix_] := Reverse[matrix];
d1sym[matrix_] := Transpose[matrix];
r1sym[matrix_] := vsym[d1sym[matrix]];
d2sym[matrix_] := r1sym[vsym[matrix]];

(*List of all distinct symmetries of given matrix *)
symlist[matrix_] := 
  DeleteDuplicates[{matrix, vsym[matrix], hsym[matrix], d1sym[matrix], 
    d2sym[matrix], r1sym[matrix], r1sym[r1sym[matrix]], 
    r1sym[r1sym[r1sym[matrix]]]}];

(* All square binary matrices of size mdim *) 
binlist[mdim_] := binlist = Tuples[{0, 1}, {mdim, mdim}];

(* Initialise list of equivalence classes *)
binclasses = {};

binlist[2];

While[ Length[binlist] >= 1,
 (currmat = binlist[[1]];
   currsym = symlist[currmat];
   binclasses = Append[binclasses, {currmat, Length[currsym]}];
   binlist = DeleteCases[binlist, x_ /; MemberQ[currsym, x] == True]);
 ]

Answer 1

SECOND ATTEMPT - 2nd Edition

The following functions are defined as pure functions for speed.

toMatrix = Partition[IntegerDigits[#, 2, 16], 4] &;
toInteger = FromDigits[Flatten[#], 2] &;
allSyms = Module[{s1, s3, s4, s6}, {#, s1 = Reverse[#, {2}], Reverse[#], 
 s3 = Transpose[#], s4 = Reverse[s3, {2}], 
 Reverse[Transpose[s1], {2}], s6 = Reverse[Transpose[s4], {2}], 
 Reverse[Transpose[s6], {2}]}] &;

You can check what they do using for example:

MatrixForm /@ allSyms[toMatrix[253]]

Finally, the following code gives you the answer in 4 seconds 3.8 seconds!

casesToCheck = Range[0, 65535];
Timing[answer = {MatrixForm@toMatrix@First[#], 
 Length[#]} & /@ (Reap[
  NestWhile[
   Complement[#, 
     Sow[Union[toInteger /@ allSyms[toMatrix[#[[1]]]]]]] &, 
   casesToCheck, (n = Length[#]) > 0 &]][[2, 1]]); Length[answer]]
(*{3.796875, 8548}*)

You can check the answers now. For example, the first 3 are answer[[1 ;; 3]]:

The classes with only one element are Cases[answer, {_, 1}, {1}] :

---

Old slow solution

Have you tried the following?

matrices = Partition[IntegerDigits[#, 2, 4], 2] & /@ Range[0, 15];;
vsym[matrix_] := Reverse[matrix, {2}];
hsym[matrix_] := Reverse[matrix];
d1sym[matrix_] := Transpose[matrix];
r1sym[matrix_] := vsym[d1sym[matrix]];
d2sym[matrix_] := r1sym[vsym[matrix]];
Tally[matrices,
   vsym[#1] == #2 || hsym[#1] === #2 || d1sym[#1] == d1sym[#2] || 
   r1sym[#1] == #2 || d2sym[#1] == #2&]

Answer 2

Here is an implementation of the change I suggested to Hector's code:

toMatrix = Partition[IntegerDigits[#, 2, 16], 5] &;
toInteger = FromDigits[Flatten[#], 2] &;
allSyms = 
  Module[{s1, s3, s4, s6}, {#, s1 = Reverse[#, {2}], Reverse[#], 
     s3 = Transpose[#], s4 = Reverse[s3, {2}], 
     Reverse[Transpose[s1], {2}], s6 = Reverse[Transpose[s4], {2}], 
     Reverse[Transpose[s6], {2}]}] &;
casesToCheck = Range[0, 2^25];

Timing[answer = {MatrixForm@toMatrix@First[#], Length[#]} & /@ (Reap[
      With[{syms = Union[toInteger /@ (allSyms@toMatrix[#])]},
         If[# == syms[[1]], Sow[syms]];
         ] & /@ casesToCheck
      ][[2, 1]]); Length[answer]]

(*{6.593587, 8548}*)

For comparison: Four seconds on Hector's machine corresponds to about six seconds on mine. This solution is a little slower than Hector's, but it should scale better to the 5x5 case. I was surprised not to get a speed-up; this is a linear-time algorithm, but apparently the cost of calling Union[toInteger/@(allSyms@toMatrix[#])] for every single matrix is enough to nullify any advantage in the 4x4 case.

I thought it would be interesting to implement a few other (still linear-time) solutions using Hector's matrix functions and compare performance.

Here is the same solution, but with allSyms memoized (i.e. the integers corresponding to the matrices symmetric to the matrix corresponding to a given integer are remembered in the definition of the function sortedSymInts):

clearSortedSymInts[] := (
   Clear[sortedSymInts];
   sortedSymInts[matint_] := 
    With[{foo = Union[toInteger /@ allSyms[toMatrix[matint]]]},
     (sortedSymInts[#] = foo) & /@ foo;
     foo
    ]
  );

clearSortedSymInts[];
Timing[answer2 = {MatrixForm@toMatrix@First[#], 
     Length[#]} & /@ (Reap[
      With[{syms = sortedSymInts[#]},
         If[# == syms[[1]], Sow[syms]];
         ] & /@ casesToCheck
      ][[2, 1]]); Length[answer]
 ]

(*{1.736454, 8548}*)

Now it's much faster, at the expense of requiring more memory.

Here are two reimplementations of my previous solution, using Hector's matrix functions. The first memoizes integers representing matrices symmetric to matrices which have already been seen. The second memoizes the matrices themselves, and only uses toMatrix to generate the initial list.

First:

clearNewSym[] := (
   Clear[newSym];
   newSym[
     matint_] := ((newSym[#] = False) & /@ (sortedSymInts[matint]); 
     True)
   );

clearSortedSymInts[];
clearNewSym[];
Timing[binclasses = {MatrixForm[toMatrix[#]], 
     Length[sortedSymInts[#]]} & /@ Select[casesToCheck, newSym]; 
 Length[binclasses]]

(*{1.534411, 8548}*)

Second:

clearNewerSym[] := (
   Clear[newerSym];
   newerSym[mat_] := ((newerSym[#] = False) & /@ allSyms[mat]; True)
 );
clearSortedSymInts[];
clearNewerSym[];
Timing[binclasses2 = {MatrixForm[#], Length[Union[allSyms[#]]]} & /@ 
   Select[toMatrix /@ casesToCheck, newerSym]; Length[binclasses2]]

(*{1.241104, 8548}*)

Surprisingly (to me), the last is the fastest. Apparently it is faster to memoize the matrices than to call toInteger and then memoize the integers. The hash function used by Mathematica's rule processing code must be very fast for lists (or perhaps slow for integers???).

Answer 3

I am rewriting this answer.

an[u_] := Module[{n, perm, p, tp, ptp}, n = Length[u];
   perm = 
    Flatten[Outer[List, Permutations@Range[n], Permutations@Range[n], 
      1], 1];
   p = Map[u[[#[[1]], #[[2]]]] &, perm];
   tp = Transpose /@ p;
   ptp = Union[p, tp]
   ];

eqc = Union /@ Gather[an /@ test];

Displaying data:

{Map[MatrixForm, #, {2}], Length @@ #} & /@eqc // Grid

Answer 4

The following can replace your second Timing block. The memoizing function newsymrep returns True if no matrix equivalent to the input has previously been passed as input.

symlist[matrix_] := {matrix, vsym[matrix], hsym[matrix], 
   d1sym[matrix], d2sym[matrix], r1sym[matrix], r1sym[r1sym[matrix]], 
   r1sym[r1sym[r1sym[matrix]]]};

Clear[newsymrep];
newsymrep[matrix_] := 
  With[{foo = symlist[matrix]}, (newsymrep[#] = False) & /@ foo;
   True];

Timing[
 binreps = Select[blist, newsymrep];
 binclasses = {MatrixForm[#], Length@Union@symlist[#]} & /@ binreps;
 Length[binreps]]

(*{1.513778, 8548}*)

I removed DeleteCases from symlist; as symlist is used here, it's faster to redundantly memoize.