Finally fixed selfsim logic

13 files changed, 83 insertions, 63 deletions

diff --git a/TODO.txt b/TODO.txt
index a52945f..b164e23 100644
--- a/TODO.txt
+++ b/TODO.txt
@@ -1,4 +1,7 @@
-Bug in esep table generation
+H48 table generation
+  - fix gendata tool
+    - use to test generation of full h0k4 table
+    - remove ifdef and old code
   - compute all tables for h<11
     x compute visited up to a fixed depth 8
     - compute additional step (if needed) to fill <=base
@@ -10,9 +13,7 @@ Bug in esep table generation
   - Add long-running test for h0k4 (maybe as a tool?)
     - tests for other sizes?
   - optimize
-    - use cached values for invcoord_esep? check if it is faster
-    - Unify and improve performance of finding nasty sim
-      (right now we don't even use selfsim)
+    - use only transform_edges (need compose trans)
     - parallelize with pthread
 
 (OLD:
diff --git a/src/solve_h48.h b/src/solve_h48.h
index ad77156..ad865a7 100644
--- a/src/solve_h48.h
+++ b/src/solve_h48.h
@@ -27,6 +27,24 @@
 
 #define MAX_SOLUTION_LENGTH 20
 
+/*
+TODO: This loop other similar h48 coordinates can be improved by only
+transforming edges, but we need to compose transformations (i.e. conjugate
+_t by _ttrep).
+*/
+#define _foreach_h48sim(_cube, _cocsepdata, _selfsim, _h, _action) \
+        int64_t _cocsep = coord_cocsep(_cube); \
+        uint8_t _ttrep = TTREP(_cocsepdata[_cocsep]); \
+        int64_t _coclass = COCLASS(_cocsepdata[_cocsep]); \
+        cube_t _rep = transform(_cube, _ttrep); \
+        uint64_t _sim = _selfsim[_coclass]; \
+        for (uint8_t _t = 0; _t < 48 && _sim; _t++, _sim >>= 1) { \
+                if (!(_sim & 1)) continue; \
+                _cube = transform(_rep, _t); \
+                _cube = transform(_cube, inverse_trans(_ttrep)); \
+                _action \
+        }
+
 typedef struct {
         uint64_t n;
         uint64_t capacity;
@@ -229,12 +247,12 @@ coord_h48(cube_t c, const uint32_t *cocsepdata, uint8_t h)
 }
 
 _static_inline int64_t
-coord_h48_edges(cube_t c, int64_t coclass, uint8_t t, uint8_t h)
+coord_h48_edges(cube_t c, int64_t coclass, uint8_t ttrep, uint8_t h)
 {
         cube_t d;
         int64_t esep, eo, edges;
 
-        d = transform_edges(c, t);
+        d = transform_edges(c, ttrep);
         esep = coord_esep(d);
         eo = coord_eo(d);
         edges = (esep << 11) + eo;
@@ -394,11 +412,11 @@ gen_h48short(
         const uint64_t *selfsim,
         h48map_t *map
 ) {
-        uint8_t i, m, t;
-        int64_t coord, cc;
-        uint64_t j, oldn, sim;
+        uint8_t i, m;
+        int64_t coord;
+        uint64_t j, oldn;
         kvpair_t kv;
-        cube_t cube, d, e;
+        cube_t cube, d;
 
         cube = solvedcube();
         coord = coord_h48(cube, cocsepdata, 11);
@@ -418,17 +436,10 @@ gen_h48short(
                         cube = invcoord_h48(kv.key, crep, 11);
                         for (m = 0; m < 18; m++) {
                                 d = move(cube, m);
-                                coord = coord_h48(d, cocsepdata, 11);
-                                h48map_insertmin(map, coord, i+1);
-                                cc = coord / H48_ESIZE(11);
-                                sim = selfsim[cc] >> UINT64_C(1);
-                                for (t = 1; t < 48 && sim; t++) {
-                                /* TODO: optimize by using transform_edges,
-                                   corner coordinate is kept */
-                                        e = transform(d, t);
-                                        coord = coord_h48(e, cocsepdata, 11);
+                                _foreach_h48sim(d, cocsepdata, selfsim, 11,
+                                        coord = coord_h48(d, cocsepdata, 11);
                                         h48map_insertmin(map, coord, i+1);
-                                }
+                                )
                         }
                 }
                 LOG("found %" PRIu8 "\n", map->n-oldn);
@@ -512,9 +523,8 @@ gendata_h48h0k4_bfs_fromdone(bfsarg_esep_t *arg)
 {
         uint8_t c, m, x;
         uint32_t cc;
-        int64_t i, j, k, t, cocsep_coord;
-        uint64_t sim;
-        cube_t cube, moved, transd;
+        int64_t i, j, k;
+        cube_t cube, moved;
 
         for (i = 0, cc = 0; i < (int64_t)ESEP_MAX(0); i++) {
                 c = get_esep_pval(arg->buf32, i);
@@ -522,30 +532,16 @@ gendata_h48h0k4_bfs_fromdone(bfsarg_esep_t *arg)
                         continue;
                 cube = invcoord_h48(i, arg->crep, 0);
                 for (m = 0; m < 18; m++) {
-                        /*
-                         * TODO: here we can optimize by computing at first
-                         * only the corner part of the coordinate, and then
-                         * the edge parts for each transformation.
-                         */
                         moved = move(cube, m);
                         j = coord_h48(moved, arg->cocsepdata, 0);
-                        x = get_esep_pval(arg->buf32, j);
-                        if (x <= arg->depth)
+                        if (get_esep_pval(arg->buf32, j) <= arg->depth)
                                 continue;
-                        set_esep_pval(arg->buf32, j, arg->depth);
-                        cc += x != arg->depth;
-                        cocsep_coord = j / H48_ESIZE(0);
-                        sim = arg->selfsim[cocsep_coord] >> UINT64_C(1);
-                        for (t = 1; t < 48 && sim; t++) {
-                                /* TODO: use only selfsim */
-                                transd = transform(moved, t);
-                                k = coord_h48(transd, arg->cocsepdata, 0);
+                        _foreach_h48sim(moved, arg->cocsepdata, arg->selfsim, 0,
+                                k = coord_h48(moved, arg->cocsepdata, 0);
                                 x = get_esep_pval(arg->buf32, k);
-                                if (x <= arg->depth)
-                                        continue;
                                 set_esep_pval(arg->buf32, k, arg->depth);
                                 cc += x != arg->depth;
-                        }
+                        )
                 }
         }
 
@@ -557,8 +553,8 @@ gendata_h48h0k4_bfs_fromnew(bfsarg_esep_t *arg)
 {
         uint8_t c, m, x;
         uint32_t cc;
-        int64_t i, j, t, cocsep_coord, sim;
-        cube_t cube, moved, transd;
+        int64_t i, j;
+        cube_t cube, moved;
 
         for (i = 0, cc = 0; i < (int64_t)ESEP_MAX(0); i++) {
                 c = get_esep_pval(arg->buf32, i);
@@ -569,22 +565,33 @@ gendata_h48h0k4_bfs_fromnew(bfsarg_esep_t *arg)
                         moved = move(cube, m);
                         j = coord_h48(moved, arg->cocsepdata, 0);
                         x = get_esep_pval(arg->buf32, j);
-                        if (x < arg->depth)
-                                goto neighbor_found;
-                }
-                continue;
-neighbor_found:
-                set_esep_pval(arg->buf32, i, arg->depth);
-                cc++;
-                cocsep_coord = i / H48_ESIZE(0);
-                sim = arg->selfsim[cocsep_coord] >> 1;
-                for (t = 1; t < 48 && sim; t++) {
-                        /* TODO: use only selfsim */
-                        transd = transform(cube, t);
-                        j = coord_h48(transd, arg->cocsepdata, 0);
-                        x = get_esep_pval(arg->buf32, j);
-                        set_esep_pval(arg->buf32, j, arg->depth);
-                        cc += x == 0xF;
+                        if (x >= arg->depth)
+                                continue;
+#if 0
+                        cube_t transd;
+                        int64_t t, cocsep_coord, sim;
+
+                        set_esep_pval(arg->buf32, i, arg->depth);
+                        cc++;
+                        cocsep_coord = i / H48_ESIZE(0);
+                        sim = arg->selfsim[cocsep_coord] >> 1;
+                        for (t = 1; t < 48 && sim; t++) {
+                                transd = transform(cube, t);
+                                j = coord_h48(transd, arg->cocsepdata, 0);
+                                x = get_esep_pval(arg->buf32, j);
+                                set_esep_pval(arg->buf32, j, arg->depth);
+                                cc += x == 0xF;
+                        }
+#else
+                        _foreach_h48sim(cube, arg->cocsepdata, arg->selfsim, 0,
+                                j = coord_h48(cube, arg->cocsepdata, 0);
+                                x = get_esep_pval(arg->buf32, j);
+                                set_esep_pval(arg->buf32, j, arg->depth);
+                                cc += x == 0xF;
+                        )
+#endif
+
+                        break;
                 }
         }
 
diff --git a/test/111_gendata_h48_h0/00_h_0.in b/test/111_gendata_h48_h0/00_h_0.in
deleted file mode 100644
index e69de29..0000000
--- a/test/111_gendata_h48_h0/00_h_0.in
+++ /dev/null
diff --git a/test/112_h48map/00_small.in b/test/111_h48map/00_small.in
index ee8e591..ee8e591 100644
--- a/test/112_h48map/00_small.in
+++ b/test/111_h48map/00_small.in
diff --git a/test/112_h48map/00_small.out b/test/111_h48map/00_small.out
index 260e804..260e804 100644
--- a/test/112_h48map/00_small.out
+++ b/test/111_h48map/00_small.out
diff --git a/test/112_h48map/01_large.in b/test/111_h48map/01_large.in
index 1812b40..1812b40 100644
--- a/test/112_h48map/01_large.in
+++ b/test/111_h48map/01_large.in
diff --git a/test/112_h48map/01_large.out b/test/111_h48map/01_large.out
index f5e00d8..f5e00d8 100644
--- a/test/112_h48map/01_large.out
+++ b/test/111_h48map/01_large.out
diff --git a/test/112_h48map/h48map_tests.c b/test/111_h48map/h48map_tests.c
index d70bc3a..d70bc3a 100644
--- a/test/112_h48map/h48map_tests.c
+++ b/test/111_h48map/h48map_tests.c
diff --git a/test/112_gendata_h48/00_h_0.in b/test/112_gendata_h48/00_h_0.in
new file mode 100644
index 0000000..cb6e5ed
--- /dev/null
+++ b/test/112_gendata_h48/00_h_0.in
@@ -0,0 +1,2 @@
+5
+0
diff --git a/test/111_gendata_h48_h0/00_h_0.out b/test/112_gendata_h48/00_h_0.out
index 04e8bf6..04e8bf6 100644
--- a/test/111_gendata_h48_h0/00_h_0.out
+++ b/test/112_gendata_h48/00_h_0.out
diff --git a/test/111_gendata_h48_h0/gendata_h48_tests.c b/test/112_gendata_h48/gendata_h48_tests.c
index e9de3c0..485d9d0 100644
--- a/test/111_gendata_h48_h0/gendata_h48_tests.c
+++ b/test/112_gendata_h48/gendata_h48_tests.c
@@ -1,20 +1,30 @@
 #include "../test.h"
 
-#define MAXDEPTH 5
 #define COCSEPSIZE 1119792
 #define ETABLESIZE ((3393 * 495 * 70) >> 1)
 
 int64_t gendata_h48h0k4(void *, uint8_t);
 
+int64_t gendata_h48_fixture(void *buf, uint8_t maxdepth, uint8_t h) {
+        if (h == 0)
+                return gendata_h48h0k4(buf, maxdepth);
+        fprintf(stderr, "Error: gendata h48 for h>0 not implemented yet\n");
+        exit(1);
+}
+
 void run(void) {
         char str[STRLENMAX];
-        uint8_t i;
+        uint8_t i, maxdepth, h;
         uint32_t *buf, *h48info;
         size_t result;
 
         fgets(str, STRLENMAX, stdin);
+        maxdepth = atoi(str);
+        fgets(str, STRLENMAX, stdin);
+        h = atoi(str);
+
         buf = (uint32_t *)malloc(sizeof(uint32_t) * 60000000);
-        result = gendata_h48h0k4(buf, MAXDEPTH);
+        result = gendata_h48_fixture(buf, maxdepth, h);
         h48info = buf + (ETABLESIZE + COCSEPSIZE) / 4;
 
         printf("%zu\n\n", result);
@@ -26,7 +36,7 @@ void run(void) {
                 printf("%" PRIu32 ": %" PRIu32 "\n", i, buf[COCSEPSIZE/4-10+i]);
 
         printf("\nh48:\n");
-        for (i = 0; i < MAXDEPTH+1; i++)
+        for (i = 0; i < maxdepth+1; i++)
                 printf("%" PRIu32 ": %" PRIu32 "\n", i, h48info[i+1]);
 
         free(buf);
diff --git a/tools/benchmark_gendata_h48/benchmark_gendata_h48.c b/tools/01_gendata_h48/benchmark_gendata_h48.c
index efb6f6f..efb6f6f 100644
--- a/tools/benchmark_gendata_h48/benchmark_gendata_h48.c
+++ b/tools/01_gendata_h48/benchmark_gendata_h48.c
diff --git a/tools/stats_tables_h48/stats_tables_h48.c b/tools/02_stats_tables_h48/stats_tables_h48.c
index 7c6bc3c..7c6bc3c 100644
--- a/tools/stats_tables_h48/stats_tables_h48.c
+++ b/tools/02_stats_tables_h48/stats_tables_h48.c