From bdecf4034597e7672fa0cedd705847c582f69617 Mon Sep 17 00:00:00 2001 From: David Allemang Date: Fri, 12 Nov 2021 15:06:29 -0500 Subject: [PATCH] Introduce complex solvers --- .gitignore | 2 + Testing/common.hpp | 12 +- examples/CMakeLists.txt | 3 + examples/bench.cpp | 5 +- examples/complex.cpp | 17 +++ examples/path.cpp | 3 +- include/tc/complex.hpp | 284 ++++++++++++++++++++++++++++++++++++++++ include/tc/cosets.hpp | 5 +- include/tc/group.hpp | 37 ++++-- include/tc/solver.hpp | 18 ++- 10 files changed, 363 insertions(+), 23 deletions(-) create mode 100644 examples/complex.cpp create mode 100644 include/tc/complex.hpp diff --git a/.gitignore b/.gitignore index 6447641..da21595 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,5 @@ *.[oa] cmake-build* + +Testing/Temporary diff --git a/Testing/common.hpp b/Testing/common.hpp index c4aa908..f080707 100644 --- a/Testing/common.hpp +++ b/Testing/common.hpp @@ -8,23 +8,23 @@ #include #include -std::vector parse_vec(const std::string &part) { +tc::Symbol parse_vec(const std::string &part) { std::istringstream iss(part); - std::vector res; + std::vector vec; std::string token; while (std::getline(iss, token, ' ')) { - res.push_back(std::stoul(token)); + vec.push_back(std::stoul(token)); } - return res; + return Eigen::Map(vec.data(), vec.size()); } size_t compute( const tc::Group &group, - const std::vector &vgens + const tc::Symbol &gens ) { - auto table = tc::solve(group, vgens); + auto table = tc::solve(group, gens); return table.order(); } diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 63a2aa8..9616e67 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -6,3 +6,6 @@ target_link_libraries(path PRIVATE tc) add_executable(group group.cpp) target_link_libraries(group PRIVATE tc) + +add_executable(complex complex.cpp) +target_link_libraries(complex PRIVATE tc) diff --git a/examples/bench.cpp b/examples/bench.cpp index 08904eb..362e4e2 100644 --- a/examples/bench.cpp +++ b/examples/bench.cpp @@ -7,8 +7,11 @@ template void test(const G &group) { + tc::Symbol gens(1); + gens << 0; + auto s = std::clock(); - auto cosets = tc::solve(group, {0}); + auto cosets = tc::solve(group, gens); auto e = std::clock(); double diff = (double) (e - s) / CLOCKS_PER_SEC; diff --git a/examples/complex.cpp b/examples/complex.cpp new file mode 100644 index 0000000..f9fc881 --- /dev/null +++ b/examples/complex.cpp @@ -0,0 +1,17 @@ +#include +#include +#include + +int main() { + tc::Symbol symbol(3); + symbol << 5, 3, 3; + auto group = tc::schlafli(symbol); + + constexpr int N = 4; + std::vector combos = tc::combinations(group.gens, N - 1); + auto data = tc::merge(tc::hull<4>(group, combos, {})); + + std::cout << data.size() << std::endl; + + return 0; +} diff --git a/examples/path.cpp b/examples/path.cpp index da7c584..6675e2a 100644 --- a/examples/path.cpp +++ b/examples/path.cpp @@ -6,8 +6,9 @@ #include int main() { + tc::Symbol gens(0); auto cube = tc::group::B(3); - auto vars = tc::solve(cube); + auto vars = tc::solve(cube, gens); std::string start; std::vector names = {"a", "b", "c"}; diff --git a/include/tc/complex.hpp b/include/tc/complex.hpp new file mode 100644 index 0000000..cf6028f --- /dev/null +++ b/include/tc/complex.hpp @@ -0,0 +1,284 @@ +#pragma once + +#include +#include +#include +#include +#include +#include + +namespace tc { + std::vector combinations(const Symbol &symbol, size_t srank) { + size_t rank = symbol.size(); + + std::vector mask(rank, false); + std::fill_n(mask.begin(), srank, true); + + std::vector 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 + struct Primitive { + static_assert(N > 0, "Primitives must contain at least one point. Primitive<0> or lower is impossible."); + + std::array inds; + + Primitive() = default; + + Primitive(const Primitive &) = default; + + Primitive(const Primitive &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 generators(const tc::Group &context) { + std::vector 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 + std::vector> apply(std::vector> 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 + void flip(std::vector> 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 + [[nodiscard]] + std::vector> recontext( + std::vector> 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> res(prims); + for (Primitive &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 + std::vector> merge(const std::vector>> &meshes) { + size_t size = 0; + for (const auto &mesh: meshes) { + size += mesh.size(); + } + + std::vector> res; + res.reserve(size); + for (const auto &mesh: meshes) { + res.insert(res.end(), mesh.begin(), mesh.end()); + } + + return res; + } + + template + [[nodiscard]] + std::vector>> each_tile( + std::vector> prims, + const tc::Group &context, + const Symbol &g_gens, + const Symbol &sg_gens + ) { + std::vector> 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 + [[nodiscard]] + std::vector> tile( + std::vector> prims, + const tc::Group &context, + const Symbol &g_gens, + const Symbol &sg_gens + ) { + auto res = each_tile(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 + [[nodiscard]] + std::vector> fan(std::vector> prims, int root) { + std::vector> res(prims.size()); + std::transform(prims.begin(), prims.end(), res.begin(), + [root](const Primitive &prim) { + return Primitive(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 + std::vector> 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>> meshes; + + for (const auto &sg_gens: combos) { + auto base = triangulate(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> 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> res; + res.emplace_back(); + return res; + } + + template + auto hull(const tc::Group &group, T all_sg_gens, const std::vector &exclude) { + std::vector>> 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(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; + } +} diff --git a/include/tc/cosets.hpp b/include/tc/cosets.hpp index fc7002e..5e04529 100644 --- a/include/tc/cosets.hpp +++ b/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 - std::vector walk(const T &start, const F &op) { + std::vector walk(const T &start, const F &op) const { std::vector res; res.reserve(order()); res.push_back(start); @@ -74,7 +75,7 @@ namespace tc { } template - std::vector walk(const T &start, const E &gens, const F &op) { + std::vector 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]); }); diff --git a/include/tc/group.hpp b/include/tc/group.hpp index 7ac586f..08112dc 100644 --- a/include/tc/group.hpp +++ b/include/tc/group.hpp @@ -6,7 +6,7 @@ #include -namespace { +namespace tc { template 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 */ diff --git a/include/tc/solver.hpp b/include/tc/solver.hpp index e0bd4b0..5e664c8 100644 --- a/include/tc/solver.hpp +++ b/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 &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); + } }