path: root/doc/h48.md

# The H48 optimal solver

This document contains information on the H48 Rubik's Cube optimal solver.
This solver is occasionally improved and optimized, and this document
is updated accordingly to reflect the current implementation.

I highly encourage the reader to check out Jaap Scherphuis'
[Computer Puzzling page](https://www.jaapsch.net/puzzles/compcube.htm)
before reading this document.

Other backround material that will be referenced throughout the text includes:

* [Cube coordinates](https://sebastiano.tronto.net/speedcubing/coordinates)
  by Sebastiano Tronto
* [Nxopt](https://github.com/rokicki/cube20src/blob/master/nxopt.md)
  by Tomas Rokicki
* [The Mathematics behind Cube Explorer](https://kociemba.org/cube.htm)
  by Herber Kociemba

Initially I intended to give minimum information here, citing the
above resources when needed, but I ended up rewriting even some of the
basic things in this document.

## Optimal solvers

### IDA*

The basic idea behind the solver is an iterative-deepening
depth-first search with pruning tables, also known as IDA*. You
can find a detailed explanation of this idea in
[Jaap's page](https://www.jaapsch.net/puzzles/compcube.htm#tree).

To summarize, the IDA* algorithm finds the shortest solution(s) for a
scrambled puzzle by using multiple depth-first searches at increasing
depth. A breadth-first search is not applicable in this scenario because
of the large amount of memory required.  At each depth, the puzzle's state
is evaluated to obtain an estimated lower bound for how many moves are
required to solve the puzzle; this is done by employing large pruning
tables, as well as other techniques. If the lower bound exceeds the
number of moves required to reach the current target depth, the branch is
pruned. Otherwise, every possible move is tried, except those that would
"cancel" with the preceding moves, obtaining new positions at depth N+1.

### Pruning tables and coordinates

A pruning table associates to a cube position a value that is a
lower bound for the number of moves it takes to solve that position.
To acces this value, the cube must be turned into an index for the
given table. This is done using a
[*coordinate*](https://sebastiano.tronto.net/speedcubing/coordinates),
by which we mean a function from the set of (admissible, solvable) cube
positions to the set of non-negative integers. The maximum value N of a
coordinate must be less than the size of the pruning table; preferably
there are no "gaps", that is the coordinate reaches all values between
0 and N and the pruning table has size N+1.

Coordinates are often derived from the cosets of a *target subgroup*.
As an example, consider the pruning table that has one entry for each
possible corner position, and the stored values denote the length of
an optimal corner solution for that position. In this situation, the
target subgroup is the set of all cube positions with solved corners
(for the Mathematician: this is indeed a subgroup of the group of all
cube positions). This group consists of `12! * 2^10` positions.  The set
of cosets of this subgroup can be identified with the set of corner
configurations, and it has size `8! * 3^7`.  By looking at corners,
from each cube position we can compute a number from `0` to `8! * 3^7 - 1`
and use it as an index for the pruning table.

The pruning table can be filled in many ways. A common way is starting
with the coordinate of the solved cube (often 0) and visiting the set
of all coordinates with a breadth-first search. To do this, one needs
to apply a cube move to a coordinate; this can either be done directly
(for example using transistion tables, as explained in Jaap's page) or
by converting the coordinate back into a valid cube position, applying
the move and then computing the coordinate of the new position.

The largest a pruning table is, the more information it contains,
the faster the solvers that uses it is. It is therefore convenient to
compress the tables as much as possible. This can be achieved by using
4 bits per position to store the lower bound, or as little as 2 bits
per position with some caveats (see Jaap's page or Rokicki's nxopt
description for two possible ways to achieve this).  Tables that use
only 1 bit per position are possible, but the only information they
can provide is wether the given position requires more or less than a
certain number of moves to be solved.

### Symmetries

Another, extremely effective way of reducing the size of a pruning table
is using symmetries. In the example above of the corner table, one can
notice that all corner positions that are one quarter-turn away from being
solved are essentially the same: they can be turned one into the other
by rotating the cube and optionally applying a mirroring transformation.
Indeed, every position belongs to a set of up to 48 positions that are
"essentially the same" in this way.  Keeping a separate entry in the
pruning table for each of these positions is a waste of memory.

There are ways to reduce a coordinate by symmetry, grouping every
transformation-equivalent position into the same index. You can
find a description in my
[cube coordinates page](https://sebastiano.tronto.net/speedcubing/coordinates).
By doing so, the size of the pruning table is reduced without loss of
information.

Some well-known solvers (Cube Explorer, nxopt) do not take advantage
of the full group of 48 symmetries of the cube, but they use only the
16 symmetries that fix the UD axis. However, they make up for this by
doing 3 pruning table lookups, one for each axis of the cube.

Some cube positions are self-symmetric: when applying certain
transformations to them, they remain the same. For example, the
cube obtained by applying U2 to a solved cube is invariant with
respect to the mirroring along the RL-plane. This fact has a couple
of consequences:

* Reducing by a group of N symmetries does not reduce the size of
  the pruning table by a factor of N, but by a slightly smaller
  factor. If the size of the coordinate is large, this fact is
  harmless.
* When combining symmetry-reduced coordinates with other coordinates,
  one has to be extra careful with handling self-symmetric coordinates.
  An informal explanation of this phenomenon can be found in
  the [docs/nasty-symmetries.md document](./nasty-symmetries.md).

## The H48 solver

### The target subgroup: H48 coordinates

The H48 solver uses a target group that is invariant under all 48
symmetries. This group is defined as follows:

* Each corner is placed in the *corner tetrad* it belongs to, and it
  is oriented. Here by corner tetrad we mean one of the two sets of
  corners {UFR, UBL, DFL, DBR} and {UFL, UBR, DFR, DBL}. For a corner
  that is placed in its own tetrad, the orientation does not depend on
  the reference axis chosen, and an oriented corner can be defined
  has having the U or D sticker facing U or D.
* Each edge is placed in the *slice* it belongs to. The three edge slices
  are M = {UF, UB, DB, DF}, S = {UR, UL, DL, DR} and E = {FR, FL, BL, BR}.

Other options are available for corners: for example, we could
impose that the corner permutation is even or that corners are in
*half-turn reduction* state, that is they are in the group generated
by {U2, D2, R2, L2, F2, B2}. We settled on the corner group described
above because it is easy to calculate and it gives enough options
for pruning table sizes, depending on how edge orientation is
considered (see the next section).

We call the coordinate obtained from the target subgroup described
above the **h0** coordinate. This coordinate has `C(8,4) * 3^7 *
C(12,8) * C(8,4) = 5.304.568.500` possible values.  If the corner
part of the coordinate is reduced by symmetry, it consists of only
3393 elements, giving a total of `3393 * C(12,8) * C(8,4) =
117.567.450` values. This is not very large for a pruning table,
as with 2 bits per entry it would occupy less than 30MB. However,
it is possible to take edge orientation (partially) into account
to produce larger tables.

### Edge orientation

Like for corners, the orientation of an edge is well-defined independently
of any axis when said edge is placed in the slice it belongs to. We may
modify the target subgroup defined in the previous section by imposing
that the edges are all oriented. This yields a new coordinate, which we
call **h11**, whose full size with symmetry-reduced corners is `3393 *
C(12,8) * C(8,4) * 2^11 = 240.778.137.600`. A pruning table based on
this coordinate with 2 bits per entry would occupy around 60GB, which
is a little too much for most personal computers.

One can wonder if it is possible to use a coordinate that considers
the orientation of only *some* of the edges, which we may call **h1**
to **h10**. Such coordinates do exist, but they are not invariant under
the full symmetry group: indeed an edge whose orientation we keep track
of could be moved to any of the untracked edges' positions by one of
the symmetries, making the whole coordinate ill-defined.

It is however possible to compute the symmetry-reduced pruning tables
for these coordinates. One way to construct them is by taking the **h11**
pruning table and "forgetting" about some of the edge orientation values,
collapsing 2 (or a power thereof) values to one by taking the minimum.
It is also possible to compute these tables directly, as explained in
the **Pruning table computation** section below.

### Coordinate computation for pruning value estimation

In order to access the pruning value for a given cube position, we first
need to compute its coordinate value. Let's use for simplicity the **h0**
coordinate, but everything will be valid for any other coordinate of
the type described in the previous section.

As described in the symmetric-composite coordinate section of
[my coordinates page](https://sebastiano.tronto.net/speedcubing/coordinates),
we first need to compute the part of the coordinate where the symmetry
is applied, that is the corner separation + orientation (also called
*cocsep*). The value of the coordinate depends on the equivalence
class of the corner configuration, and is memorized in a table called
`cocsepdata`. To avoid lengthy computations, the index for a cube position
in this table is determined by the corner orientation (between `0` and
`3^7-1`) and a binary representation of the corner separation part (an
8-bit number with 4 zero digits and 4 one digits, where zeros denote
corners in one tetrad and 1 denotes corners in the other; the last bit
can be ignored, as it can easily be deduced from the other 7). This is
slightly less space-efficient than computing the actual corner-separation
coordinate as a value between `0` and `C(8,4)`, but it is much faster.

We can now access the value in the `cocsepdata` table corresponding to
our cube position. This table contains the value of the symmetry-reduced
coordinate and the so-called *ttrep* (transformation to representative)
necessary to compute the full **h0** coordinate. For convenience,
this table also stores a preliminary pruning value based on the corner
coordinate.  If this value is large enough, we may skip the computation of
the **h0** coordinate and any further estimate. These three values can be
stored in a single 32 bit integer, so that the table uses less than 12MB.

### Getting the pruning value

Once we have computed the full h0 coordinate, we can access
the correct entry in the full pruning table. We chose to
use 2 bits per entry, as described by Rokicki in the [nxopt
document](https://github.com/rokicki/cube20src/blob/master/nxopt.md).
This means that every pruning table needs to have an associated *base
value*, that determines the offset to be added to each entry (each entry
can only be 0, 1, 2 or 3).  If the base value is `b`, a pruning value
of 1, 2 or 3 can be used directly as a lower bound of b+1, b+2 and b+3
respectively. However, a value of 0 could mean that the actual lower
bound is anything between 0 and b, so we cannot take b as a lower bound.

To be able to still use some sort of pruning value even when we get a
0 read, we use a **fallback table**. Inspired by nxopt, this table is
interleaved with the main table for cache efficiency: every 254 entries
(508 bits), we store in 4 bits the minimum of the pruning values of these
254 entries as a number from 0 to 16, without subtracting the table's
base. Since a cache line is 512 bits long, we never get a cache miss
when looking up the minimum value in the fallback table after a 0 read.
Smaller lines (of 256 or 128 bits) have been tried, but they do not
give any significant improvement over 512 bit lines.

Moreover, as an additional heuristic, in case of a 0 read we also look
up another pruning value in a table that takes into account only the
position of the edges. This table is small (around 1MB), so repeated
accesses to it are not too slow. In practice, this gives a small speed up
of around 5%.  More tables could be used to refine the fallback estimate,
but each additional table leads to longer lookup times, especially if
it is too large to fit in cache.

Previous versions of this implementation (up to commit 6c42463, or
to version 0.2) also included the possibility of a table with 4 bits
per entry, at least for the `h0` case. Such tables did not require
a fallback lookup, but due to their large size they were not less
efficient. Therefore, they have been removed.

### Estimation refinements

After computing the pruning value, there are a number of different tricks
that can be used to improve the estimation.

#### Inverse estimate

A cube position and its inverse will, in general, give a different
pruning value. If the normal estimate is not enough to prune the branch,
the inverse of the position can be computed and a new estimate can be
obtained from there.

#### Reducing the branching factor - tight bounds

*This trick and the next are and adaptation of a similar technique
introduced by Rokicki in nxopt.*

We can take advantage of the fact that the **h0** (and **h11**) coordinate
is invariant under the subgroup `<U2, D2, R2, L2, F2, B2>`. Suppose we
are looking for a solution of length `m`, we are at depth `d` and the
pruning value `p` *for the inverse position* is a strict bound, that is
`d+p = m`.  In this situation, from the inverse point of view the solution
cannot end with a 180° move: if this was the case, since these moves are
contained in the target group for **h0** (or **h11**), it must be that
we were in the target subgroup as well *before the last move*, i.e. we
have found a solution for the **h0** coordinate of length `p-1`. But this
is in contradiction with the fact that the inverse pruning value is `p`.

From all of this, we conclude that trying any 180° move as next move is
useless (because it would be the last move from the inverse position).
We can therefore reduce the branching factor significantly and only try
quarter-turn moves.

#### Reducing the branching factor - switching to inverse

We can expand on the previous trick by using a technique similar to NISS.

Suppose that we have a strict bound for the *normal position*. As above,
we can deduce that the last move cannot be a 180° move, but this does not
tell us anything about the possibilities for the next move. However, if
we *switch to the inverse scramble* we can take advantage of this fact
as described above. For doing this, we need to replace the cube with its
inverse, and keep track of the moves done from now on so that we can
invert them at the end to construct the final solution.

When this technique is used, it is also possible to avoid some table
lookups: if the last move applied to the cube is a 180° move *on the
inverse position*, then the coordinate on the normal position has not
changed. Thus if we keep track of the last computed pruning value,
we can reuse it and avoid an expensive table lookup.

### Other optimizations

The H48 solver uses various other optimizations.

#### Avoiding inverse cube computation

Computing the inverse of a cube position is expensive. We avoid doing
that (for the inverse pruning value estimate) if we bring along both
the normal and the inverse cube during the search, and we use *premoves*
to apply moves to the inverse scramble.

#### Multi-threading

The solution search is parallelized *per-scramble*. Although when solving
multiple positions it would be more efficient to solve them in parallel
by dedicating a single thread to each of them, the current H48 solver
optimizes for solving a single cube at the time.

To do this, we first compute all cube positions that are 4 moves away
from the scramble. This way we prepare up to 43254 *tasks* (depending
on the symmetries of the starting position, which we take advantage of)
which are equally split between the available threads.  Any solution
encountered in this step is of course added to the list of solutions.

#### Filtering out symmetric cases

If the cube position we wish to solve is, say, mirrored along the RL
plane, it is redundant to attempt starting a solution with both R and L',
as they would lead to mirrored versions of the same solution. Therefore,
we filter out some tasks by looking for symmetries before each of the
first 4 moves. Then, unless we are looking for only one solution, we
have to make sure to report back all symmetric variations of a solution,
so we need to keep track of which transformations were available at each
of these 4 points.

This can reduce the number of starting positions drastically: for
example, for the superflip our method produces 1038 starting tasks -
a 41.67x improvement! Of course this optimization has no effect on any
position that isn't 3 moves or fewer away from a symmetric position -
which is the vast majority of them.

#### Heuristically sorting tasks

The tasks described in the previous paragraph (multi-threading) are
initially searched in an arbitrary order. However, after searching
at a sufficient depth, we have gathered some data that allows us to
make some heuristical improvements: the tasks that leads to visiting
more positions (or in other words, where we go over the estimated lower
bounds less often), are more likely to yield the optimal solution. Thus
we sort the tasks based on this, adjusting by a small factor due to the
fact that sequences ending in U, R or F moves have more continuations
than those ending in D, L or B moves - as we don't allow, for example,
both U D and D U, but only the former.

Preliminary benchmarks show a performance improvement of around 40%
when searching a single solution. When searching for multiple optimal
solutions the effects of this optimization will be less pronounced, and
they are obviously inexistent when looking for *all* optimal solutions.

## Pruning table computation

We first describe how the pruning table would be computed if we stored
each pruning value fully, using 4 bits per entry. This algorithm used
to be implemented (up to commit 6c42463), but has since been removed.

### Full tables for h0 and h11 (not implemented)

If we are allowed to store the full pruning value for each entry,
computing the pruning tables is not that difficult.  First, we set the
value of the solved cube's coordinate to 0 and every other entry to 15,
the largest 4-bit integer. Then we iteratively scan through the table and
compute the coordinates at depth n+1 from those at depth n as follows:
for each coordinate at depth n, we compute a valid representative for it,
we apply each of the possible 18 moves to it, and we set their value to
n+1. Once the table is filled or we have reached depth 15, we stop.

Unfortunately what I described above is an oversimplification: one
also must take into account the case where the corner coordinate is
self-symmetric, and handle it accordingly by making sure that every
symmetric variation of the same coordinate has its value set at the
right iteration.

When this algorithm reaches the last few iterations, the coordinates at
depth n are many more than those whose depth is still unknown. Therefore
it is convenient to look at the unset coordinate and check if they have
any neighbor whose depth is i; if so, we can set such a coordinate to i+1.
Thanks to [torchlight](https://github.com/torchlight) for suggesting
this optimization.

Finally, this algorithm can be easily parallelized by dividing the set
of coordinates into separate sections, but one must take into account
that a coordinate and its neighbors are usually not in the same section.

### 2 bit tables

Computing the pruning tables for intermediate coordinates is not as
simple.  The reason is that these coordinates are not invariant under
the 48 transformations of the cube, but we are treating them as such.
Take for example the case of **h10**. Each **h10** coordinate is
represented by either of two **h11** coordinates: the one where the 11th
edge is correctly oriented and the one where this edge is flipped (of
course, the orientation of the 12th edge can be deduced from the parity
of the orientation of the other 11). We can take any cube whose **h11**
coordinate is one of these two to represent our **h10** coordinate, but
the set of neighbors will change depending on which representative we
picked. In other words: if I take one of the two cubes and make a move,
the position reached is not necessarily obtainable by applying a move
to the other cube. This means that the same algorithm that we described
for the "real coordinate" case cannot be applied here.

It is still possible to do a brute-force depth-first seach to fill
these tables.  To make it reasonably fast, we apply a couple of
optimizations. First of all, we restrict ourselves to computing a 2 bits
table, so we can stop the search early.

The second optimization we employ consists in pre-computing all possible
**h11** coordinates at a fixed depth, and storing their depth in a hash
map indexed by their **h11** coordinate value. This way we do not have
to brute-force our way through from depth 0, and it make this method
considerably faster. Unfortunately, since the number of coordinates
at a certain depth increases exponentially with the depth, we are for
now limited at storing the coordinates at depth 8.

This method works for the "real" coordinates **h0** and **h11** as well.
Moreover, in this case one can optimize it further by avoiding to repeat
the search from a coordinate that has already been visited. Further
optimization are possible for **h0** and **h11**, and we may implement
them in the future.