From 3fe8cc8d166b67ea81343dc17fb506afba693730 Mon Sep 17 00:00:00 2001
From: David Allemang <david.allemang@kitware.com>
Date: Mon, 6 Feb 2023 12:07:43 -0500
Subject: [PATCH] ENH: Use Eigen array for mesh indices

WIP: Rename "Prims" to "Indices"

WIP: Use Eigen array for mesh indices

WIP: Remvoe Primitive::apply
---
 vis/include/geometry.hpp |   2 +-
 vis/include/solver.hpp   | 195 ++++++++++++++-------------------------
 vis/src/comps.hpp        |   4 +-
 3 files changed, 70 insertions(+), 131 deletions(-)

diff --git a/vis/include/geometry.hpp b/vis/include/geometry.hpp
index 972b0a7..718edac 100644
--- a/vis/include/geometry.hpp
+++ b/vis/include/geometry.hpp
@@ -11,7 +11,7 @@
 #include "combo_iterator.hpp"
 
 template<unsigned N>
-using Prims = Eigen::Matrix<unsigned, N, Eigen::Dynamic>;
+using Indices = Eigen::Array<unsigned, N, Eigen::Dynamic>;
 
 template<int N>
 using vec = Eigen::Matrix<float, N, 1>;
diff --git a/vis/include/solver.hpp b/vis/include/solver.hpp
index 1545d6a..760e7e5 100644
--- a/vis/include/solver.hpp
+++ b/vis/include/solver.hpp
@@ -5,41 +5,13 @@
 #include <optional>
 #include <numeric>
 #include <iostream>
+
+#include <geometry.hpp>
+
 #include "combo_iterator.hpp"
 
-#include "fmt/core.h"
-#include "fmt/ranges.h"
-
 /**
- * An primitive stage N indices.
- * @tparam N
- */
-template<unsigned N>
-struct Primitive {
-    static_assert(N > 0, "Primitives must contain at least one point. Primitive<0> or lower is impossible.");
-
-    std::array<unsigned, N> indices;
-
-    Primitive() = default;
-
-    Primitive(const Primitive<N> &) = default;
-
-    Primitive(const Primitive<N - 1> &sub, unsigned root) {
-        std::copy(sub.indices.begin(), sub.indices.end(), indices.begin());
-        indices[N - 1] = root;
-    }
-
-    ~Primitive() = default;
-
-    void apply(const tc::Cosets &table, int gen) {
-        for (auto &ind: indices) {
-            ind = table.get(ind, gen);
-        }
-    }
-};
-
-/**
- * Produce a list of all generators for the group context. The range [0..group.ngens).
+ * 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());
@@ -49,44 +21,34 @@ std::vector<size_t> generators(const tc::Group &context) {
 
 /**
  * Determine which of g_gens are the correct names for sg_gens within the current context
- *
- * Produces the indexes of sg_gens within g_gens; sorted.
  */
-std::vector<size_t> recontext_gens(
-    const std::vector<size_t> &g_gens,
-    const std::vector<size_t> &sg_gens) {
+std::vector<size_t> recontext_gens(const tc::Group &context, std::vector<size_t> g_gens, std::vector<size_t> sg_gens) {
+
+    std::sort(g_gens.begin(), g_gens.end());
+
+    int inv_gen_map[context.rank()];
+    for (size_t i = 0; i < g_gens.size(); i++) {
+        inv_gen_map[g_gens[i]] = i;
+    }
 
     std::vector<size_t> s_sg_gens;
-    for (const auto &gen: sg_gens) {
-        s_sg_gens.push_back(std::find(g_gens.begin(), g_gens.end(), gen) - g_gens.begin());
+    s_sg_gens.reserve(sg_gens.size());
+    for (const auto gen: sg_gens) {
+        s_sg_gens.push_back(inv_gen_map[gen]);
     }
+    std::sort(s_sg_gens.begin(), s_sg_gens.end());
 
     return s_sg_gens;
 }
 
-/**
- * Apply some context transformation to all primitives of this mesh.
- */
-template<unsigned N>
-std::vector<Primitive<N>> apply(std::vector<Primitive<N>> prims, const tc::Cosets &table, int gen) {
-    for (auto &prim: prims) {
-        prim.apply(table, gen);
-    }
-    return prims;
-}
-
 /**
  * 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>
 [[nodiscard]]
-std::vector<Primitive<N>> recontext(
-    std::vector<Primitive<N>> prims,
-    const tc::Group &context,
-    const std::vector<size_t> &g_gens,
-    const std::vector<size_t> &sg_gens
-) {
-    const auto proper_sg_gens = recontext_gens(g_gens, sg_gens);
+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 auto proper_sg_gens = recontext_gens(context, g_gens, sg_gens);
     const auto table = context.sub(g_gens).solve({});
     const auto cosets = context.sub(sg_gens).solve({});
 
@@ -97,11 +59,10 @@ std::vector<Primitive<N>> recontext(
         return table.get(coset, gen);
     }, map.begin());
 
-    std::vector<Primitive<N>> res(prims);
-    for (Primitive<N> &prim: res) {
-        for (auto &ind: prim.indices) {
-            ind = map[ind];
-        }
+    Indices<N> res(prims);
+    auto data = res.data();
+    for (int i = 0; i < prims.size(); ++i) {
+        data[i] = map[data[i]];
     }
 
     return res;
@@ -111,16 +72,18 @@ std::vector<Primitive<N>> recontext(
  * Union several meshes of the same dimension
  */
 template<unsigned N>
-std::vector<Primitive<N>> merge(const std::vector<std::vector<Primitive<N>>> &meshes) {
-    size_t size = 0;
+Indices<N> merge(const std::vector<Indices<N>> &meshes) {
+    size_t cols = 0;
     for (const auto &mesh: meshes) {
-        size += mesh.size();
+        cols += mesh.cols();
     }
 
-    std::vector<Primitive<N>> res;
-    res.reserve(size);
-    for (const auto &mesh: meshes) {
-        res.insert(res.end(), mesh.begin(), mesh.end());
+    Indices<N> res(N, cols);
+
+    size_t offset = 0;
+    for (const Indices<N> &mesh: meshes) {
+        res.middleCols(offset, mesh.cols()) = mesh;
+        offset += mesh.cols();
     }
 
     return res;
@@ -128,55 +91,39 @@ std::vector<Primitive<N>> merge(const std::vector<std::vector<Primitive<N>>> &me
 
 template<unsigned N>
 [[nodiscard]]
-std::vector<std::vector<Primitive<N>>> each_tile(
-    std::vector<Primitive<N>> prims,
-    const tc::Group &context,
-    const std::vector<size_t> &g_gens,
-    const std::vector<size_t> &sg_gens
-) {
-    std::vector<Primitive<N>> base = recontext(prims, context, g_gens, sg_gens);
-    const auto proper_sg_gens = recontext_gens(g_gens, sg_gens);
+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) {
+    Indices<N> base = recontext<N>(prims, context, g_gens, sg_gens);
+    const auto proper_sg_gens = recontext_gens(context, g_gens, sg_gens);
 
     const auto &table = context.sub(g_gens).solve({});
     const auto &cosets = context.sub(g_gens).solve(proper_sg_gens);
 
     tc::Path<> path(cosets);
 
-    std::vector<std::vector<Primitive<N>>> res(path.order());
-    path.walk(base, [&](auto from, auto to) {
-        return apply(from, table, to);
+    std::vector<Indices<N>> res(path.order());
+
+    path.walk(base, [&](Indices<N> from, auto gen) {
+        for (int i = 0; i < from.size(); ++i) {
+            from(i) = table.get(from(i), gen);
+        }
+        return from;
     }, res.begin());
 
     return res;
 }
 
-template<unsigned N>
-[[nodiscard]]
-std::vector<Primitive<N>> tile(
-    std::vector<Primitive<N>> prims,
-    const tc::Group &context,
-    const std::vector<size_t> &g_gens,
-    const std::vector<size_t> &sg_gens
-) {
-    auto res = each_tile<N>(prims, context, g_gens, sg_gens);
-
-    return merge(res);
-}
-
 /**
  * Produce a mesh of higher dimension by fanning a single point to all primitives in this mesh.
  */
 template<unsigned N>
 [[nodiscard]]
-std::vector<Primitive<N + 1>> fan(std::vector<Primitive<N>> prims, size_t root) {
-    std::vector<Primitive<N + 1>> res(prims.size());
-    std::transform(
-        prims.begin(), prims.end(),
-        res.begin(),
-        [root](const Primitive<N> &prim) {
-            return Primitive<N + 1>(prim, root);
-        }
-    );
+Indices<N + 1> fan(Indices<N> prims, int root) {
+    Indices<N + 1> res(N + 1, prims.cols());
+
+    res.topRows(1) = Indices<1>::Constant(1, prims.cols(), root);
+    res.bottomRows(N) = prims;
+
     return res;
 }
 
@@ -184,50 +131,42 @@ std::vector<Primitive<N + 1>> fan(std::vector<Primitive<N>> prims, size_t root)
  * 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>
-std::vector<Primitive<N>> triangulate(
-    const tc::Group &context,
-    const std::vector<size_t> &g_gens
-) {
-    if (g_gens.size() + 1 != N) {
+Indices<N> triangulate(const tc::Group &context, const std::vector<size_t> &g_gens) {
+    if (g_gens.size() + 1 != N) // todo make static assert
         throw std::logic_error("g_gens size must be one less than N");
-    }
 
     const auto &combos = combinations(g_gens, g_gens.size() - 1);
 
-    std::vector<std::vector<Primitive<N>>> meshes;
+    std::vector<Indices<N>> meshes;
 
     for (const auto &sg_gens: combos) {
         auto base = triangulate<N - 1>(context, sg_gens);
-        auto raised = tile(base, context, g_gens, sg_gens);
-        raised.erase(raised.begin(), raised.begin() + base.size());
-        meshes.push_back(fan(raised, 0));
+        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);
     }
 
-    return merge(meshes);
+    return merge<N>(meshes);
 }
 
-// Single-index primitives should not be further triangulated.
+/**
+ * Single-index primitives should not be further triangulated.
+ */
 template<>
-std::vector<Primitive<1>> triangulate(
-    const tc::Group &,
-    const std::vector<size_t> &g_gens
-) {
-    if (not g_gens.empty()) {
+Indices<1> triangulate<1>(const tc::Group &context, const std::vector<size_t> &g_gens) {
+    if (not g_gens.empty()) // todo make static assert
         throw std::logic_error("g_gens must be empty for a trivial Mesh");
-    }
 
-    std::vector<Primitive<1>> res;
-    res.emplace_back();
-    return res;
+    return Indices<1>::Zero(1, 1);
 }
 
 template<unsigned N, class T>
 auto hull(const tc::Group &group, T all_sg_gens, const std::vector<std::vector<size_t>> &exclude) {
-    std::vector<std::vector<Primitive<N>>> parts;
-
+    std::vector<Indices<N>> parts;
     auto g_gens = generators(group);
-
-    for (std::vector<size_t> sg_gens: all_sg_gens) {
+    for (const std::vector<size_t> &sg_gens: all_sg_gens) {
         bool excluded = false;
         for (const auto &test: exclude) {
             if (sg_gens == test) {
@@ -238,7 +177,7 @@ auto hull(const tc::Group &group, T all_sg_gens, const std::vector<std::vector<s
         if (excluded) continue;
 
         const auto &base = triangulate<N>(group, sg_gens);
-        const auto &tiles = each_tile(base, group, g_gens, sg_gens);
+        const auto &tiles = tile<N>(base, group, g_gens, sg_gens);
         for (const auto &tile: tiles) {
             parts.push_back(tile);
         }
diff --git a/vis/src/comps.hpp b/vis/src/comps.hpp
index 4d3816b..6e3060d 100644
--- a/vis/src/comps.hpp
+++ b/vis/src/comps.hpp
@@ -34,7 +34,7 @@ namespace vis {
         };
 
         cgl::Buffer<vec4> verts;
-        cgl::Buffer<Primitive<4>> ibo;
+        cgl::Buffer<Eigen::Array<unsigned int, 4, 1>> ibo;
         cgl::Buffer<ModelMatrix> ubo;
 
         using Affine4f = Eigen::Transform<float, 4, Eigen::Affine>;
@@ -66,7 +66,7 @@ namespace vis {
 
             auto inds = merge<N>(hull<N>(group.group, group.include, group.exclude));
 
-            vbos.ibo.put(inds.begin(), inds.end());
+            vbos.ibo.put(inds.colwise().begin(), inds.colwise().end());
         }
     }