me/toddcox-faster
1
0
Fork 0
forked from mirror/toddcox-faster
Code Pull Requests Activity

Un-template complex solver

Browse Source
This commit is contained in:
David Allemang
2021-11-12 22:14:46 -05:00
parent cd8ef050e6
commit a1b0ae31dc
2 changed files with 78 additions and 121 deletions
Download Patch File Download Diff File Expand all files Collapse all files

27
examples/complex.cpp
View File

@@ -5,11 +5,11 @@
#include <tc/complex.hpp>
#include <tc/groups.hpp>
template<unsigned int N, class G>
void test(const G &group) {
template<class G>
void test(const G &group, unsigned int N) {
auto s = std::clock();
auto combos = tc::combinations(group.gens, N - 1);
auto data = tc::merge<N>(tc::hull<N>(group, combos, {}));
auto data = tc::merge(tc::hull(group, combos, {}));
auto e = std::clock();
double diff = (double) (e - s) / CLOCKS_PER_SEC;
@@ -24,17 +24,16 @@ void test(const G &group) {
}
int main() {
test<4>(tc::group::H(4));
test<4>(tc::group::B(4));
test<4>(tc::group::B(5));
test<4>(tc::group::B(6));
test<4>(tc::group::E(6));
test<3>(tc::group::H(3));
test<3>(tc::group::B(4));
test<3>(tc::group::B(5));
test<3>(tc::group::B(6));
test<3>(tc::group::E(6));
test(tc::group::H(4), 4);
test(tc::group::B(4), 4);
test(tc::group::B(5), 4);
test(tc::group::B(6), 4);
test(tc::group::E(6), 4);
test(tc::group::H(3), 3);
test(tc::group::B(4), 3);
test(tc::group::B(5), 3);
test(tc::group::B(6), 3);
test(tc::group::E(6), 3);
return 0;
}

172
include/tc/complex.hpp
View File

@@ -6,6 +6,7 @@
#include <optional>
#include <numeric>
#include <iostream>
#include <utility>
namespace tc {
std::vector<Symbol> combinations(const Symbol &symbol, size_t srank) {
@@ -30,38 +31,41 @@ namespace tc {
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.");
Symbol fan(const Symbol &prim, unsigned root) {
Symbol res(prim.size() + 1);
res << prim, root;
return res;
}
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;
std::vector<Symbol> fan(const std::vector<Symbol>& prims, int root) {
std::vector<Symbol> res;
res.reserve(prims.size());
for (const auto &prim: prims) {
Symbol s(prim.size() + 1);
s << prim, root;
res.push_back(s);
}
return res;
}
~Primitive() = default;
inline void flip() {
if (N > 1) std::swap(inds[0], inds[1]);
void flip(Symbol &prim) {
if (prim.size() > 1) {
std::swap(prim[0], prim[1]);
}
}
void apply(const tc::Cosets &table, unsigned int gen) {
for (auto &ind: inds) {
ind = table.get(ind, gen);
}
flip();
void apply(const tc::Cosets &table, unsigned int gen, Symbol &prim) {
for (auto &ind: prim) {
ind = table.get(ind, gen);
}
};
flip(prim);
}
void apply(const tc::Cosets &table, unsigned int gen, std::vector<Symbol> &prims) {
for (auto &prim: prims) {
apply(table, gen, prim);
}
}
/**
* Produce a list of all generators for the group context. The range [0..group.rank).
@@ -94,35 +98,20 @@ namespace tc {
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) {
void flip(std::vector<Symbol> &prims) {
for (auto &prim: prims) {
prim.flip();
flip(prim);
}
}
/**
* 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,
void recontext(
std::vector<Symbol> &prims,
const tc::Group &context,
const Symbol &g_gens,
const Symbol &sg_gens
@@ -136,30 +125,26 @@ namespace tc {
return table.get(coset, gen);
});
std::vector<Primitive<N>> res(prims);
for (Primitive<N> &prim: res) {
for (auto &ind: prim.inds) {
for (Symbol &prim: prims) {
for (auto &ind: prim) {
ind = map[ind];
}
}
if (get_parity(context, g_gens, sg_gens) == 1)
flip(res);
return res;
flip(prims);
}
/**
* Union several meshes of the same dimension
*/
template<unsigned N>
std::vector<Primitive<N>> merge(const std::vector<std::vector<Primitive<N>>> &meshes) {
std::vector<Symbol> merge(const std::vector<std::vector<Symbol>> &meshes) {
size_t size = 0;
for (const auto &mesh: meshes) {
size += mesh.size();
}
std::vector<Primitive<N>> res;
std::vector<Symbol> res;
res.reserve(size);
for (const auto &mesh: meshes) {
res.insert(res.end(), mesh.begin(), mesh.end());
@@ -168,100 +153,73 @@ namespace tc {
return res;
}
template<unsigned N>
[[nodiscard]]
std::vector<std::vector<Primitive<N>>> each_tile(
std::vector<Primitive<N>> prims,
std::vector<std::vector<Symbol>> each_tile(
std::vector<Symbol> base,
const tc::Group &context,
const Symbol &g_gens,
const Symbol &sg_gens
) {
std::vector<Primitive<N>> base = recontext(prims, context, g_gens, sg_gens);
recontext(base, 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);
auto res = path.walk(base, _gens, [&table](auto from, auto &gen) {
apply(table, gen, from);
return from;
});
return res;
}
template<unsigned N>
[[nodiscard]]
std::vector<Primitive<N>> tile(
std::vector<Primitive<N>> prims,
std::vector<Symbol> tile(
std::vector<Symbol> base,
const tc::Group &context,
const Symbol &g_gens,
const Symbol &sg_gens
) {
auto res = each_tile<N>(prims, context, g_gens, sg_gens);
auto res = each_tile(std::move(base), 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(
std::vector<Symbol> 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");
if (g_gens.size() == 0) {
std::vector<Symbol> res;
Symbol prims(1);
prims.setZero();
res.push_back(prims);
return res;
}
const auto &combos = combinations(g_gens, g_gens.size() - 1);
std::vector<std::vector<Primitive<N>>> meshes;
std::vector<std::vector<Symbol>> meshes;
for (const auto &sg_gens: combos) {
auto base = triangulate<N - 1>(context, sg_gens);
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));
auto fanned = fan(raised, 0);
meshes.push_back(fanned);
}
return merge(meshes);
const std::vector<Symbol> &result = merge(meshes);
return result;
}
/**
* 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>
template<class T>
auto hull(const tc::Group &group, T all_sg_gens, const std::vector<Symbol> &exclude) {
std::vector<std::vector<Primitive<N>>> parts;
std::vector<std::vector<Symbol>> parts;
auto g_gens = group.gens;
for (const Symbol &sg_gens: all_sg_gens) {
bool excluded = false;
@@ -273,7 +231,7 @@ namespace tc {
}
if (excluded) continue;
const auto &base = triangulate<N>(group, sg_gens);
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);