Introduce complex solvers

include/tc/complex.hpp

@@ -0,0 +1,284 @@
+#pragma once
+
+#include <tc/group.hpp>
+#include <tc/solver.hpp>
+#include <cmath>
+#include <optional>
+#include <numeric>
+#include <iostream>
+
+namespace tc {
+    std::vector<Symbol> combinations(const Symbol &symbol, size_t srank) {
+        size_t rank = symbol.size();
+
+        std::vector<bool> mask(rank, false);
+        std::fill_n(mask.begin(), srank, true);
+
+        std::vector<Symbol> combos;
+        combos.reserve(choose(rank, srank));
+
+        Symbol row(srank);
+        do {
+            for (int j = 0, k = 0; j < rank; ++j) {
+                if (mask[j]) {
+                    row(k++) = symbol(j);
+                }
+            }
+            combos.emplace_back(row);
+        } while (std::prev_permutation(mask.begin(), mask.end()));
+
+        return combos;
+    }
+
+/**
+ * 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> inds;
+
+        Primitive() = default;
+
+        Primitive(const Primitive<N> &) = default;
+
+        Primitive(const Primitive<N - 1> &sub, unsigned root) {
+            std::copy(sub.inds.begin(), sub.inds.end(), inds.begin());
+            inds[N - 1] = root;
+        }
+
+        ~Primitive() = default;
+
+        inline void flip() {
+            if (N > 1) std::swap(inds[0], inds[1]);
+        }
+
+        void apply(const tc::Cosets &table, unsigned int gen) {
+            for (auto &ind: inds) {
+                ind = table.get(ind, gen);
+            }
+            flip();
+        }
+    };
+
+/**
+ * Produce a list of all generators for the group context. The range [0..group.rank).
+ */
+    std::vector<unsigned int> generators(const tc::Group &context) {
+        std::vector<unsigned int> g_gens(context.rank());
+        std::iota(g_gens.begin(), g_gens.end(), 0);
+        return g_gens;
+    }
+
+/**
+ * Determine whether the orientation of the group sg_gens is reversed from the group g_gens within group context
+ */
+    int get_parity(
+        const tc::Group &context,
+        const Symbol &g_gens,
+        const Symbol &sg_gens
+    ) {
+        if (g_gens.size() != sg_gens.size() + 1) return 0;
+
+        const auto proper_sg_gens = recontext_gens(context.rank(), g_gens, sg_gens);
+
+        int i = 0;
+        for (; i < sg_gens.size(); ++i) {
+            if (proper_sg_gens[i] != i) {
+                break;
+            }
+        }
+
+        return i & 1;
+    }
+
+
+/**
+ * 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, unsigned int gen) {
+        for (auto &prim: prims) {
+            prim.apply(table, gen);
+        }
+        return prims;
+    }
+
+/**
+ * Reverse the orientation of all primitives in this mesh.
+ */
+    template<unsigned N>
+    void flip(std::vector<Primitive<N>> prims) {
+        for (auto &prim: prims) {
+            prim.flip();
+        }
+    }
+
+/**
+ * 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 Symbol &g_gens,
+        const Symbol &sg_gens
+    ) {
+        const auto proper_sg_gens = recontext_gens(context.rank(), g_gens, sg_gens);
+
+        const auto table = solve(context, g_gens, Symbol(0));
+        const auto path = solve(context, sg_gens, Symbol(0)).path();
+
+        auto map = path.walk(0U, proper_sg_gens, [&table](auto coset, auto gen) {
+            return table.get(coset, gen);
+        });
+
+        std::vector<Primitive<N>> res(prims);
+        for (Primitive<N> &prim: res) {
+            for (auto &ind: prim.inds) {
+                ind = map[ind];
+            }
+        }
+
+        if (get_parity(context, g_gens, sg_gens) == 1)
+            flip(res);
+
+        return res;
+    }
+
+/**
+ * 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;
+        for (const auto &mesh: meshes) {
+            size += mesh.size();
+        }
+
+        std::vector<Primitive<N>> res;
+        res.reserve(size);
+        for (const auto &mesh: meshes) {
+            res.insert(res.end(), mesh.begin(), mesh.end());
+        }
+
+        return res;
+    }
+
+    template<unsigned N>
+    [[nodiscard]]
+    std::vector<std::vector<Primitive<N>>> each_tile(
+        std::vector<Primitive<N>> prims,
+        const tc::Group &context,
+        const Symbol &g_gens,
+        const Symbol &sg_gens
+    ) {
+        std::vector<Primitive<N>> base = recontext(prims, context, g_gens, sg_gens);
+
+        const auto table = solve(context, g_gens, Symbol(0));
+        const auto path = solve(context, g_gens, sg_gens).path();
+
+        auto _gens = generators(context);
+
+        auto res = path.walk(base, _gens, [&table](auto from, auto gen){
+            return apply(from, table, gen);
+        });
+
+        return res;
+    }
+
+    template<unsigned N>
+    [[nodiscard]]
+    std::vector<Primitive<N>> tile(
+        std::vector<Primitive<N>> prims,
+        const tc::Group &context,
+        const Symbol &g_gens,
+        const Symbol &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, int 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);
+            }
+        );
+        return res;
+    }
+
+/**
+ * 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 Symbol &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;
+
+        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));
+        }
+
+        return merge(meshes);
+    }
+
+/**
+ * Single-index primitives should not be further triangulated.
+ */
+    template<>
+    std::vector<Primitive<1>> triangulate(
+        const tc::Group &context,
+        const Symbol &g_gens
+    ) {
+        if (g_gens.size() != 0) // 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;
+    }
+
+    template<unsigned N, class T>
+    auto hull(const tc::Group &group, T all_sg_gens, const std::vector<Symbol> &exclude) {
+        std::vector<std::vector<Primitive<N>>> parts;
+        auto g_gens = group.gens;
+        for (const Symbol &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 = each_tile(base, group, g_gens, sg_gens);
+            for (const auto &tile: tiles) {
+                parts.push_back(tile);
+            }
+        }
+        return parts;
+    }
+}

include/tc/cosets.hpp

@@ -33,6 +33,7 @@ namespace tc {
         }
 
         [[nodiscard]] size_t order() const {
+            if (!_rank) return 0;
             return data.size() / _rank;
         }
 
@@ -60,7 +61,7 @@ namespace tc {
         }
 
         template<class T, class F>
-        std::vector<T> walk(const T &start, const F &op) {
+        std::vector<T> walk(const T &start, const F &op) const {
             std::vector<T> res;
             res.reserve(order());
             res.push_back(start);
@@ -74,7 +75,7 @@ namespace tc {
         }
 
         template<class T, class E, class F>
-        std::vector<T> walk(const T &start, const E &gens, const F &op) {
+        std::vector<T> walk(const T &start, const E &gens, const F &op) const {
             return walk(start, [&](const T &s, const int g) {
                 return op(s, gens[g]);
             });

include/tc/group.hpp

@@ -6,7 +6,7 @@
 
 #include <Eigen/Eigen>
 
-namespace {
+namespace tc {
     template<class T>
     std::string stringify(const T &vec) {
         std::stringstream ss;
@@ -56,16 +56,6 @@ namespace tc {
         return res;
     }
 
-    Symbol inverse(size_t rank, const Symbol &gens) {
-        size_t srank = gens.size();
-        Symbol res(rank);
-        res.fill(0);
-        for (int i = 0; i < srank; ++i) {
-            res(gens(i)) = i;
-        }
-        return res;
-    }
-
     unsigned int factorial(unsigned int n) {
         unsigned int res = 1;
         for (int i = 1; i <= n; ++i) {
@@ -101,6 +91,31 @@ namespace tc {
         return res;
     }
 
+    /**
+     * Determine which of g_gens are the correct names for sg_gens within the current context
+     */
+    Symbol recontext_gens(
+        size_t rank,
+        Symbol g_gens,
+        Symbol sg_gens
+    ) {
+        std::sort(g_gens.begin(),  g_gens.end());
+        std::sort(sg_gens.begin(),  sg_gens.end());
+
+        int inv_gen_map[rank];
+        for (int i = 0; i < g_gens.size(); ++i) {
+            inv_gen_map[g_gens[i]] = i;
+        }
+
+        Symbol s_sg_gens(sg_gens.size());
+        for (int i = 0; i < sg_gens.size(); ++i) {
+            s_sg_gens[i] = inv_gen_map[sg_gens[i]];
+        }
+        std::sort(s_sg_gens.begin(),  s_sg_gens.end());
+
+        return s_sg_gens;
+    }
+
     /**
      * Create a named coxeter matrix from a simplified schlafli symbol
      */

include/tc/solver.hpp

@@ -166,7 +166,7 @@ namespace tc {
     /**
      * Assumes that g is a coxeter group - that is, self-adjoint and the diagonal is 2.
      */
-    tc::Cosets solve(const Group &group, const std::vector<unsigned int> &sub_gens = {}) {
+    tc::Cosets solve(const Group &group, const Symbol &s_gens) {
         size_t rank = group.rank();
 
         tc::Cosets cosets(rank);
@@ -176,7 +176,7 @@ namespace tc {
             return cosets;
         }
 
-        for (unsigned int gen: sub_gens) {
+        for (unsigned int gen: s_gens) {
             if (gen < rank)
                 cosets.put(0, gen, 0);
         }
@@ -223,4 +223,18 @@ namespace tc {
 
         return cosets;
     }
+
+/**
+ * Solve the cosets generated by sg_gens within the subgroup generated by g_gens of the group context
+ */
+    Cosets solve(
+        const Group &context,
+        const Symbol &g_gens,
+        const Symbol &sg_gens
+    ) {
+        const Symbol &proper_sg_gens = recontext_gens(context.rank(), g_gens, sg_gens);
+        const Group &group = subgroup(context, g_gens);
+
+        return solve(group, proper_sg_gens);
+    }
 }