From 92a55f8c746832766e5248ff682c59a0c01a1068 Mon Sep 17 00:00:00 2001
From: David Allemang <david.allemang@kitware.com>
Date: Tue, 7 Feb 2023 16:15:19 -0500
Subject: [PATCH] ENH: Simplify solving logic

This is a combination of multiple WIP commits:

- Semantic group names - G, H
- Simplify fan<>()
- Template int deduction
- Remove _recontext_gens
- Move g_gens logic into caller
- Move subgroup logic into caller
- Move recontext logic into coller
- Move _generators, _recontext_gens to anonymous namespace
- Simplify hull, tile
---
 tc/include/tc/core.hpp         |   7 +-
 tc/src/group.cpp               |  42 ++++++++--
 vis/include/combo_iterator.hpp |  30 --------
 vis/include/geometry.hpp       |   2 -
 vis/include/solver.hpp         | 137 ++++++++++-----------------------
 vis/src/comps.hpp              |  11 +--
 vis/src/main.cpp               |   9 +--
 7 files changed, 91 insertions(+), 147 deletions(-)
 delete mode 100644 vis/include/combo_iterator.hpp

diff --git a/tc/include/tc/core.hpp b/tc/include/tc/core.hpp
index d4852ba..74c11af 100644
--- a/tc/include/tc/core.hpp
+++ b/tc/include/tc/core.hpp
@@ -92,6 +92,7 @@ namespace tc {
     private:
         size_t _rank;
         std::vector<size_t> _mults;
+        std::vector<size_t> _gens;
 
     public:
         Group(Group const &) = default;
@@ -108,7 +109,11 @@ namespace tc {
 
         [[nodiscard]] size_t rank() const;
 
-        [[nodiscard]] Group sub(std::vector<size_t> const &idxs) const;
+        [[nodiscard]] std::vector<size_t> gens() const;
+
+        [[nodiscard]] Group sub(std::vector<size_t> const &gens) const;
+
+        [[nodiscard]] std::vector<Group> subs(size_t rank) const;
 
         [[nodiscard]] Cosets solve(std::vector<size_t> const &idxs = {}, size_t bound = SIZE_MAX) const;
     };
diff --git a/tc/src/group.cpp b/tc/src/group.cpp
index 38ad0cb..831531f 100644
--- a/tc/src/group.cpp
+++ b/tc/src/group.cpp
@@ -1,9 +1,12 @@
 #include <tc/core.hpp>
 
 #include <cassert>
+#include <numeric>
 
 namespace tc {
-    Group::Group(size_t rank) : _rank(rank), _mults(_rank * _rank, 2) {
+    Group::Group(size_t rank) : _rank(rank), _mults(_rank * _rank, 2), _gens(_rank) {
+        std::iota(_gens.begin(), _gens.end(), 0);
+
         for (int idx = 0; idx < rank; ++idx) {
             set(idx, idx, 1);
         }
@@ -28,15 +31,42 @@ namespace tc {
         return _rank;
     }
 
-    [[nodiscard]] Group Group::sub(std::vector<size_t> const &idxs) const {
-        Group res(idxs.size());
+    [[nodiscard]] std::vector<size_t> Group::gens() const {
+        return _gens;
+    }
 
-        for (int i = 0; i < idxs.size(); ++i) {
-            for (int j = i; j < idxs.size(); ++j) {
-                res.set(i, j, get(idxs[i], idxs[j]));
+    [[nodiscard]] Group Group::sub(std::vector<size_t> const &gens) const {
+        Group res(gens.size());
+        res._gens = gens;
+
+        for (int i = 0; i < gens.size(); ++i) {
+            for (int j = i; j < gens.size(); ++j) {
+                res.set(i, j, get(gens[i], gens[j]));
             }
         }
 
         return res;
     }
+
+    [[nodiscard]] std::vector<Group> Group::subs(size_t rank) const {
+        std::vector<bool> mask(_rank, false);
+        std::fill(mask.begin(), mask.begin() + rank, true);
+
+        std::vector<size_t> sub_gens(rank);
+
+        std::vector<Group> res;
+
+        do {
+            for (int i = 0, j = 0; i < _rank; ++i) {
+                if (mask[i]) sub_gens[j++] = i;
+            }
+            res.push_back(sub(sub_gens));
+        } while (std::next_permutation(
+            mask.begin(),
+            mask.end(),
+            std::greater<>()
+        ));
+
+        return res;
+    }
 }
diff --git a/vis/include/combo_iterator.hpp b/vis/include/combo_iterator.hpp
deleted file mode 100644
index 3a5331f..0000000
--- a/vis/include/combo_iterator.hpp
+++ /dev/null
@@ -1,30 +0,0 @@
-#pragma once
-
-#include <set>
-#include <algorithm>
-
-template<typename V, typename M>
-V select(const V &data, const M &mask, size_t count) {
-    V result;
-    result.reserve(count);
-
-    for (int i = 0; i < mask.size(); ++i) {
-        if (mask[i]) result.push_back(data[i]);
-    }
-
-    return result;
-}
-
-template<typename V>
-std::vector<V> combinations(const V &data, const size_t count) {
-    std::vector<V> result;
-
-    std::vector<bool> mask(data.size(), false);
-    std::fill(mask.begin(), mask.begin() + count, true);
-
-    do {
-        result.push_back(select(data, mask, count));
-    } while (std::next_permutation(mask.begin(), mask.end(), std::greater<>()));
-
-    return result;
-}
diff --git a/vis/include/geometry.hpp b/vis/include/geometry.hpp
index 718edac..bd825a2 100644
--- a/vis/include/geometry.hpp
+++ b/vis/include/geometry.hpp
@@ -8,8 +8,6 @@
 
 #include <Eigen/Eigen>
 
-#include "combo_iterator.hpp"
-
 template<unsigned N>
 using Indices = Eigen::Array<unsigned, N, Eigen::Dynamic>;
 
diff --git a/vis/include/solver.hpp b/vis/include/solver.hpp
index c2a2194..67e1e4a 100644
--- a/vis/include/solver.hpp
+++ b/vis/include/solver.hpp
@@ -8,55 +8,20 @@
 
 #include <geometry.hpp>
 
-#include "combo_iterator.hpp"
-
-/**
- * Produce a list of all generators for the group context. The range [0..group.rank()).
- */
-std::vector<size_t> _generators(
-    const tc::Group &context
-) {
-    std::vector<size_t> g_gens(context.rank());
-    std::iota(g_gens.begin(), g_gens.end(), 0);
-    return g_gens;
-}
-
-/**
- * Determine which of g_gens are the correct names for sg_gens within the current context
- *
- * For example if g_gens contains {a, b, c, d} and sg_gens contains {b, d, a} then the result is {1, 3, 0}
- */
-std::vector<size_t> _recontext_gens(
-    std::vector<size_t> g_gens,
-    std::vector<size_t> sg_gens
-) {
-    auto orig = sg_gens;
-
-    for (size_t &gen: sg_gens) {
-        gen = std::distance(
-            g_gens.begin(),
-            std::find(g_gens.begin(), g_gens.end(), gen)
-        );
-    }
-    return sg_gens;
-}
-
 /**
  * Convert the indexes of this mesh to those of a different context, using g_gens to build the parent context and sg_gens to build this context.
  */
-template<unsigned N>
+template<int N>
 [[nodiscard]]
 Indices<N> recontext(
     Indices<N> prims,
-    const tc::Group &context,
-    const std::vector<size_t> &g_gens,
-    const std::vector<size_t> &sg_gens
+    const tc::Group &G,
+    const tc::Group &H
 ) {
-    const auto proper_sg_gens = _recontext_gens(g_gens, sg_gens);
-    const auto table = context.sub(g_gens).solve({});
-    const auto cosets = context.sub(sg_gens).solve({});
+    const auto table = G.solve({});
+    const auto cosets = H.solve({});
 
-    tc::Path<size_t> path(cosets, proper_sg_gens);
+    tc::Path<size_t> path(cosets, H.gens());
 
     std::vector<size_t> map(path.order());
     path.walk(0, [&table](size_t coset, size_t gen) {
@@ -75,7 +40,7 @@ Indices<N> recontext(
 /**
  * Union several meshes of the same dimension
  */
-template<unsigned N>
+template<int N>
 Indices<N> merge(
     const std::vector<Indices<N>> &meshes
 ) {
@@ -95,25 +60,21 @@ Indices<N> merge(
     return res;
 }
 
-template<unsigned N>
+template<int N>
 [[nodiscard]]
 std::vector<Indices<N>> tile(
     Indices<N> prims,
-    const tc::Group &context,
-    const std::vector<size_t> &g_gens,
-    const std::vector<size_t> &sg_gens
+    const tc::Group &G,
+    const tc::Group &H
 ) {
-    Indices<N> base = recontext<N>(prims, context, g_gens, sg_gens);
-    const auto proper_sg_gens = _recontext_gens(g_gens, sg_gens);
-
-    const auto &table = context.sub(g_gens).solve({});
-    const auto &cosets = context.sub(g_gens).solve(proper_sg_gens);
+    const auto &table = G.solve({});
+    const auto &cosets = G.solve(H.gens());
 
     tc::Path<> path(cosets);
 
     std::vector<Indices<N>> res(path.order());
 
-    path.walk(base, [&](Indices<N> from, auto gen) {
+    path.walk(prims, [&](Indices<N> from, auto gen) {
         for (int i = 0; i < from.size(); ++i) {
             from(i) = table.get(from(i), gen);
         }
@@ -126,7 +87,7 @@ std::vector<Indices<N>> tile(
 /**
  * Produce a mesh of higher dimension by fanning a single point to all primitives in this mesh.
  */
-template<unsigned N>
+template<int N>
 [[nodiscard]]
 Indices<N + 1> fan(
     Indices<N> prims,
@@ -134,8 +95,8 @@ Indices<N + 1> fan(
 ) {
     Indices<N + 1> res(N + 1, prims.cols());
 
-    res.topRows(1) = Indices<1>::Constant(1, prims.cols(), root);
-    res.bottomRows(N) = prims;
+    res.topRows(N) = prims;
+    res.bottomRows(1) = root;
 
     return res;
 }
@@ -143,67 +104,53 @@ Indices<N + 1> fan(
 /**
  * Produce a mesh of primitives that fill out the volume of the subgroup generated by generators g_gens within the group context
  */
-template<unsigned N>
-Indices<N> triangulate(
-    const tc::Group &context,
-    const std::vector<size_t> &g_gens
+template<int N>
+Indices<N> cell(
+    const tc::Group &G
 ) {
-    if (g_gens.size() + 1 != N) // todo make static assert
-        throw std::logic_error("g_gens size must be one less than N");
+    assert(G.rank() + 1 == N);  // g_gens size must be one less than N
 
-    const auto &combos = combinations(g_gens, g_gens.size() - 1);
+    std::vector<Indices<N - 1>> facets;
 
-    std::vector<Indices<N>> meshes;
+    for (auto H : G.subs(N - 2)) {
+        auto sub_base = cell<N - 1>(H);
+        auto base = recontext(sub_base, G, H);
+        auto tiles = tile(base, G, H);
 
-    for (const auto &sg_gens: combos) {
-        auto base = triangulate<N - 1>(context, sg_gens);
-        auto parts = tile<N - 1>(base, context, g_gens, sg_gens);
-        parts.erase(parts.begin(), parts.begin() + 1);
-        auto raised = merge<N - 1>(parts);
-        auto fanned = fan<N - 1>(raised, 0);
-        meshes.push_back(fanned);
+        facets.insert(
+            facets.end(),
+            tiles.begin() + 1, // skip "near" facet
+            tiles.end()
+        );
     }
 
-    return merge<N>(meshes);
+    return fan(merge(facets), 0);
 }
 
 /**
  * Single-index primitives should not be further triangulated.
  */
 template<>
-Indices<1> triangulate<1>(
-    const tc::Group &context,
-    const std::vector<size_t> &g_gens
+Indices<1> cell<1>(
+    const tc::Group &G
 ) {
-    if (not g_gens.empty()) // todo make static assert
-        throw std::logic_error("g_gens must be empty for a trivial Mesh");
+    assert(G.rank() == 0); // rank must be 0 for trivial Mesh
 
     return Indices<1>::Zero(1, 1);
 }
 
-template<unsigned N, class T>
+template<int N>
 auto hull(
-    const tc::Group &group,
-    T all_sg_gens,
-    const std::vector<std::vector<size_t>> &exclude
+    const tc::Group &G
 ) {
     std::vector<Indices<N>> parts;
-    auto g_gens = _generators(group);
-    for (const std::vector<size_t> &sg_gens: all_sg_gens) {
-        bool excluded = false;
-        for (const auto &test: exclude) {
-            if (sg_gens == test) {
-                excluded = true;
-                break;
-            }
-        }
-        if (excluded) continue;
 
-        const auto &base = triangulate<N>(group, sg_gens);
-        const auto &tiles = tile<N>(base, group, g_gens, sg_gens);
-        for (const auto &tile: tiles) {
-            parts.push_back(tile);
-        }
+    for (auto H: G.subs(N - 1)) {
+        auto sub_base = cell<N>(H);
+        auto base = recontext(sub_base, G, H);
+        auto tiles = tile(base, G, H);
+        parts.push_back(merge(tiles));
     }
+
     return parts;
 }
diff --git a/vis/src/comps.hpp b/vis/src/comps.hpp
index 559c756..68c4d9b 100644
--- a/vis/src/comps.hpp
+++ b/vis/src/comps.hpp
@@ -22,9 +22,6 @@ namespace vis {
         tc::Group group;
         vec5 root;
         vec3 color;
-
-        std::vector<std::vector<size_t>> exclude {{0, 1, 2}};
-        std::vector<std::vector<size_t>> include = combinations(_generators(group), 3);
     };
 
     struct VBOs {
@@ -62,9 +59,13 @@ namespace vis {
             vbos.verts.put(lower.colwise().begin(), lower.colwise().end());
 
             // todo generate all, then mask using glMultiDraw.
-            const size_t N = 4;
+            const Eigen::Index N = 4;
 
-            auto inds = merge<N>(hull<N>(group.group, group.include, group.exclude));
+            auto tiles = hull<N>(group.group);
+
+            tiles.erase(tiles.begin());  // remove {0, 1, 2} cells
+
+            auto inds = merge<N>(tiles);
 
             vbos.ibo.put(inds.colwise().begin(), inds.colwise().end());
         }
diff --git a/vis/src/main.cpp b/vis/src/main.cpp
index 12edef0..db05ef2 100644
--- a/vis/src/main.cpp
+++ b/vis/src/main.cpp
@@ -140,14 +140,7 @@ int run(GLFWwindow* window, ImGuiContext* ctx) {
         entity,
         tc::schlafli({5, 3, 3, 2}),
         vec5{0.80, 0.09, 0.09, 0.09, 0.09},
-        vec3{0.90, 0.90, 0.90},
-        std::vector<std::vector<size_t>>{
-            {0, 1, 2},
-
-            {0, 3, 4},
-            {1, 3, 4},
-            {2, 3, 4},
-        }
+        vec3{0.90, 0.90, 0.90}
     );
     registry.emplace<vis::VBOs>(entity);