The nine colourings of the Holt graph

autgroup.gi

LoadPackage("digraphs");
# put code describing edges here
edges := ListX([0..8], [0..2], [1, -1], {i, j, s} -> [3*i+j+1, 3*((i + s*2^(j+1)) mod 9) + (j+1) mod 3 + 1]);
edges := Concatenation(edges, List(edges, Reversed));
Gr := DigraphByEdges(edges);
V := DigraphNrVertices(Gr);
c := 3;
G := DirectProduct(AutomorphismGroup(Gr), SymmetricGroup([V+1..V+c]));

ColourShift := function(col, elt)
    local newpos;
    newpos := OnTuples([1..V+c], elt);
    return List(newpos{[1..V]}, i -> newpos[V+1+col[i]]-V-1);
end;

all_cols := [];
file := InputTextFile("cols");
while true do
    col := ReadLine(file);
    if col = fail then
        break;
    fi;
    Add(all_cols, List(NormalizedWhitespace(col), c -> IntChar(c)-48));
od;
CloseStream(file);

ineqs := List(Orbits(G, all_cols, ColourShift), orb -> SortedList(orb)[1]);
Print(ineqs, "\n");
Print(OrbitLengths(G, all_cols, ColourShift) / Factorial(c), "\n");

# [[ 0, 0, 0, 1, 0, 2, 1, 2, 0, 2, 2, 2, 0, 2, 1, 0, 1, 2, 1, 1, 1, 2, 1, 0, 2, 0, 1 ],
#  [ 0, 0, 0, 1, 1, 1, 2, 2, 2, 0, 0, 0, 1, 1, 1, 2, 2, 2, 0, 0, 0, 1, 1, 1, 2, 2, 2 ],
#  [ 0, 0, 0, 1, 1, 2, 1, 2, 2, 0, 0, 0, 1, 2, 1, 2, 1, 2, 0, 0, 0, 2, 1, 1, 2, 2, 1 ],
#  [ 0, 0, 1, 0, 0, 1, 2, 1, 2, 1, 1, 0, 2, 0, 2, 0, 0, 1, 2, 0, 2, 1, 1, 0, 2, 1, 2 ],
#  [ 0, 0, 1, 0, 0, 1, 2, 1, 2, 1, 1, 0, 2, 0, 2, 0, 2, 1, 2, 0, 2, 1, 1, 0, 2, 1, 2 ],
#  [ 0, 0, 1, 0, 0, 1, 2, 2, 1, 2, 2, 0, 1, 0, 1, 0, 0, 2, 1, 0, 1, 2, 2, 0, 2, 2, 1 ],
#  [ 0, 0, 1, 0, 0, 2, 1, 1, 2, 1, 1, 0, 2, 0, 2, 0, 0, 1, 2, 0, 2, 1, 1, 0, 2, 1, 2 ],
#  [ 0, 0, 1, 2, 0, 1, 2, 1, 1, 2, 1, 1, 2, 0, 1, 0, 0, 1, 2, 0, 2, 1, 2, 0, 2, 0, 2 ],
#  [ 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 1, 2 ] ]
# [ 3, 1, 9, 27, 27, 27, 54, 9, 1 ]

sat.py

#!/usr/bin/env python3
from itertools import product
import re
from subprocess import run
lit_re = re.compile(r"-?\d+")
cnfn = "holt.cnf"

# https://community.wolfram.com/groups/-/m/t/2426803

clauses = []

for (i, j, s) in product(range(9), range(3), (1, -1)):
    v1 = 3*i + j
    v2 = 3*((i + s*2**(j+1))%9) + (j+1)%3
    for (b1, b2) in ((1, 1), (1, -1), (-1, 1)):
        clauses.append([b1*(2*v1+1), b2*(2*v1+2), b1*(2*v2+1), b2*(2*v2+2)])
for v in range(27):
    clauses.append([-(2*v+1), -(2*v+2)])

nvars = max({abs(l) for cl in clauses for l in cl})
nclauses = len(clauses)
with open(cnfn, 'w') as f:
    print(f"p cnf {nvars} {nclauses}", file=f)
    for cl in clauses:
        print(" ".join(map(str, cl)), "0", file=f)

nsols = 0
while True:
    proc = run(["./cadical", "-q", cnfn], capture_output=True, encoding="utf-8")
    if proc.returncode != 10:
        break
    nsols += 1
    sol = list(map(int, lit_re.findall(proc.stdout)[:-1]))
    sol_s = "".join('1' if l > 0 else '0' for l in sol)
    print("".join(str(int(sol_s[2*i:2*i+2], 2)) for i in range(27)))
    with open(cnfn, 'r') as f:
        contents = f.read().partition("\n")[2]
    with open(cnfn, 'w') as f:
        print(f"p cnf {nvars} {nclauses+nsols}", file=f)
        f.write(contents)
        print(" ".join(str(-x) for x in sol), "0", file=f)