diff --git a/vis/include/combo_iterator.hpp b/vis/include/combo_iterator.hpp
index 569887a..3a5331f 100644
--- a/vis/include/combo_iterator.hpp
+++ b/vis/include/combo_iterator.hpp
@@ -1,98 +1,30 @@
 #pragma once
 
-#include <vector>
-#include <array>
+#include <set>
 #include <algorithm>
-#include <numeric>
-#include <stdexcept>
 
-size_t choose(size_t n, size_t k) {
-    if (k == 0) return 1;
-    return n * choose(n - 1, k - 1) / k;
+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<class T>
-class ComboIterator {
-private:
-    const std::vector<T> &options;
+template<typename V>
+std::vector<V> combinations(const V &data, const size_t count) {
+    std::vector<V> result;
 
-    std::vector<bool> bits;
-    std::vector<T> curr;
-    int at;
+    std::vector<bool> mask(data.size(), false);
+    std::fill(mask.begin(), mask.begin() + count, true);
 
-    void set_curr() {
-        for (int i = 0, j = 0; i < bits.size(); ++i) {
-            if (bits[i]) curr[j++] = options[i];
-        }
-    }
+    do {
+        result.push_back(select(data, mask, count));
+    } while (std::next_permutation(mask.begin(), mask.end(), std::greater<>()));
 
-public:
-    ComboIterator(const std::vector<T> &options, int k, int at = 0)
-        : options(options), bits(options.size()), curr(k), at(at) {
-        std::fill(bits.begin(), bits.begin() + k, true);
-        set_curr();
-    }
-
-    [[nodiscard]] bool operator==(const ComboIterator<T> &o) const {
-        return at == o.at;
-    }
-
-    [[nodiscard]] bool operator!=(const ComboIterator<T> &o) const {
-        return at != o.at;
-    }
-
-    auto operator*() const {
-        return curr;
-    }
-
-    const auto &operator->() const {
-        return &this;
-    }
-
-    auto operator++(int) {
-        std::prev_permutation(bits.begin(), bits.end());
-        set_curr();
-        ++at;
-        return *this;
-    }
-
-    auto operator++() &{
-        auto res = *this;
-        (*this)++;
-        return res;
-    }
-
-    auto operator--(int) {
-        std::next_permutation(bits.begin(), bits.end());
-        set_curr();
-        --at;
-        return *this;
-    }
-
-    auto operator--() &{
-        auto res = *this;
-        (*this)--;
-        return res;
-    }
-};
-
-template<class T>
-class Combos {
-private:
-    const std::vector<T> options;
-    int k;
-    int size;
-
-public:
-    Combos(const std::vector<T> &options, size_t k)
-        : options(options), k(k), size(choose(options.size(), k)) {
-    }
-
-    ComboIterator<T> begin() const {
-        return ComboIterator<T>(options, k);
-    }
-
-    ComboIterator<T> end() const {
-        return ComboIterator<T>(options, k, size);
-    }
-};
+    return result;
+}
diff --git a/vis/include/geometry.hpp b/vis/include/geometry.hpp
index dabe6c0..972b0a7 100644
--- a/vis/include/geometry.hpp
+++ b/vis/include/geometry.hpp
@@ -5,241 +5,59 @@
 #include <optional>
 #include <numeric>
 #include <iostream>
+
+#include <Eigen/Eigen>
+
 #include "combo_iterator.hpp"
 
-/**
- * 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.");
+using Prims = Eigen::Matrix<unsigned, N, Eigen::Dynamic>;
 
-    std::array<unsigned, N> indices;
+template<int N>
+using vec = Eigen::Matrix<float, N, 1>;
+template<int N>
+using mat = Eigen::Matrix<float, N, N>;
 
-    Primitive() = default;
+using vec1 = vec<1>;
+using vec2 = vec<2>;
+using vec3 = vec<3>;
+using vec4 = vec<4>;
+using vec5 = vec<5>;
 
-    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).
- */
-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
- *
- * 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> 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());
-    }
-
-    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);
-    const auto table = context.sub(g_gens).solve({});
-    const auto cosets = context.sub(sg_gens).solve({});
-
-    tc::Path<size_t> path(cosets, proper_sg_gens);
-
-    std::vector<size_t> map(path.order());
-    path.walk(0, [&table](size_t coset, size_t gen) {
-        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];
-        }
-    }
+using mat1 = mat<1>;
+using mat2 = mat<2>;
+using mat3 = mat<3>;
+using mat4 = mat<4>;
+using mat5 = mat<5>;
 
+mat4 orthographic(float left, float right, float bottom, float top, float front, float back) {
+    mat4 res = mat4();
+    res <<
+        2 / (right - left), 0, 0, -(right + left) / (right - left),
+            0, 2 / (top - bottom), 0, -(top + bottom) / (top - bottom),
+            0, 0, 2 / (front - back), -(front + back) / (front - back),
+            0, 0, 0, 1;
     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());
-    }
+mat4 perspective(float fovy, float aspect, float zNear, float zFar) {
+    float tanHalfFovy(std::tan(fovy / 2));
 
+    mat4 res = mat4::Identity();
+    res(0, 0) = 1 / (aspect * tanHalfFovy);
+    res(1, 1) = 1 / (tanHalfFovy);
+    res(2, 2) = -(zFar + zNear) / (zFar - zNear);
+    res(3, 2) = -1;
+    res(2, 3) = -(2 + zFar * zNear) / (zFar - zNear);
     return res;
 }
 
-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);
-
-    const auto table = context.sub(g_gens).solve({});
-
-    const tc::Cosets<> &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);
-    }, res.begin());
-
+mat4 translation(float x, float y, float z) {
+    mat4 res = mat4();
+    res <<
+        1, 0, 0, x,
+        0, 1, 0, y,
+        0, 0, 1, z,
+        0, 0, 0, 1;
     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);
-        }
-    );
-    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 std::vector<size_t> &g_gens
-) {
-    if (g_gens.size() + 1 != N) {
-        throw std::logic_error("g_gens size must be one less than N");
-    }
-
-    const auto &combos = Combos<size_t>(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<> &,
-    const std::vector<size_t> &g_gens
-) {
-    if (not g_gens.empty()) {
-        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<std::vector<size_t>> &exclude) {
-    std::vector<std::vector<Primitive<N>>> parts;
-
-    auto g_gens = generators(group);
-
-    for (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 = each_tile(base, group, g_gens, sg_gens);
-        for (const auto &tile : tiles) {
-            parts.push_back(tile);
-        }
-    }
-    return parts;
-}
diff --git a/vis/include/mirror.hpp b/vis/include/mirror.hpp
index b3af23a..764611e 100644
--- a/vis/include/mirror.hpp
+++ b/vis/include/mirror.hpp
@@ -6,76 +6,7 @@
 #include <vector>
 #include <algorithm>
 
-#include <Eigen/Eigen>
-
-template<unsigned N>
-using vec = std::array<float, N>;
-
-using vec1 = vec<1>;
-using vec2 = vec<2>;
-using vec3 = vec<3>;
-using vec4 = vec<4>;
-using vec5 = vec<5>;
-
-template<class V>
-V operator*(V a, const float &b) {
-    for (auto &e : a) e *= b;
-    return a;
-}
-
-template<class V>
-V operator*(const float &b, V a) {
-    for (auto &e : a) e *= b;
-    return a;
-}
-
-template<class V>
-V operator/(V a, const float &b) {
-    for (auto &e : a) e /= b;
-    return a;
-}
-
-template<class V>
-V operator+(const V &a, V b) {
-    for (int i = 0; i < a.size(); ++i) {
-        a[i] += b[i];
-    }
-    return a;
-}
-
-template<class V>
-V operator-(V a, const V &b) {
-    for (int i = 0; i < a.size(); ++i) {
-        a[i] -= b[i];
-    }
-    return a;
-}
-
-template<class V>
-void operator-=(V &a, const V &b) {
-    for (int i = 0; i < a.size(); ++i) {
-        a[i] -= b[i];
-    }
-}
-
-template<class V>
-void operator+=(V &a, const V &b) {
-    for (int i = 0; i < a.size(); ++i) {
-        a[i] += b[i];
-    }
-}
-
-template<class V>
-float length(const V &a) {
-    float sum = 0;
-    for (const auto &e : a) sum += e * e;
-    return sqrtf(sum);
-}
-
-template<class V>
-V normalized(const V &a) {
-    return a / length(a);
-}
+#include <geometry.hpp>
 
 template<class V>
 float dot(int n, const V &a, const V &b) {
@@ -86,15 +17,6 @@ float dot(int n, const V &a, const V &b) {
     return sum;
 }
 
-template<class V>
-float dot(const V &a, const V &b) {
-    float sum = 0;
-    for (int i = 0; i < a.size(); ++i) {
-        sum += a[i] * b[i];
-    }
-    return sum;
-}
-
 template<unsigned N>
 std::vector<vec<N>> mirror(const tc::Group<> &group) {
     std::vector<std::vector<float>> mirrors;
@@ -121,7 +43,7 @@ std::vector<vec<N>> mirror(const tc::Group<> &group) {
 
     std::vector<vec<N>> res;
     for (const auto &v : mirrors) {
-        vec<N> rv{};
+        vec<N> rv = vec<N>::Zero();
 
         // ortho proj
         for (int i = 0; i < std::min(v.size(), (size_t) N); ++i) {
@@ -142,9 +64,18 @@ vec<N> stereo(const vec<N + 1> &v) {
     return r;
 }
 
+template<unsigned N>
+vec<N> ortho(const vec<N + 1> &v) {
+    vec<N> r;
+    for (int i = 0; i < N; ++i) {
+        r[i] = v[i];
+    }
+    return r;
+}
+
 template<class V>
 V project(const V &vec, const V &target) {
-    return dot(vec, target) / dot(target, target) * target;
+    return vec.dot(target) / target.dot(target) * target;
 }
 
 template<class V>
@@ -160,14 +91,14 @@ V gram_schmidt_last(std::vector<V> vecs) {
         }
     }
 
-    return normalized(vecs[vecs.size() - 1]);
+    return vecs[vecs.size() - 1].normalized();
 }
 
 template<class V, class C>
 V barycentric(const std::vector<V> &basis, const C &coords) {
-    V res{};
+    V res = V::Zero();
 
-    int N = std::min(basis.size(), coords.size());
+    int N = std::min((int) basis.size(), (int) coords.rows());
     for (int i = 0; i < N; ++i) {
         res += basis[i] * coords[i];
     }
@@ -186,29 +117,12 @@ std::vector<V> plane_intersections(std::vector<V> normals) {
     return results;
 }
 
-Eigen::Matrix4f utilRotate(const int u, const int v, const float theta) {
-    Eigen::Matrix4f res;
-    res.setIdentity();
+template<unsigned N>
+mat<N> rot(int u, int v, float theta) {
+    mat<N> res = mat<N>::Identity();
     res(u, u) = std::cos(theta);
     res(u, v) = std::sin(theta);
     res(v, u) = -std::sin(theta);
     res(v, v) = std::cos(theta);
     return res;
 }
-
-Eigen::Matrix4f ortho(
-    float l,
-    float r,
-    float b,
-    float t,
-    float n,
-    float f
-) {
-    Eigen::Matrix4f res;
-    res <<
-        2 / (r - l), 0, 0, -(r + l) / (r - l),
-        0, 2 / (t - b), 0, -(t + b) / (t - b),
-        0, 0, -2 / (f - n), -(f + n) / (f - n),
-        0, 0, 0, 1;
-    return res;
-}
diff --git a/vis/include/solver.hpp b/vis/include/solver.hpp
new file mode 100644
index 0000000..4034892
--- /dev/null
+++ b/vis/include/solver.hpp
@@ -0,0 +1,245 @@
+#pragma once
+
+#include <tc/core.hpp>
+#include <cmath>
+#include <optional>
+#include <numeric>
+#include <iostream>
+#include "combo_iterator.hpp"
+
+/**
+ * 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).
+ */
+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
+ *
+ * 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> 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());
+    }
+
+    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);
+    const auto table = context.sub(g_gens).solve({});
+    const auto cosets = context.sub(sg_gens).solve({});
+
+    tc::Path<size_t> path(cosets, proper_sg_gens);
+
+    std::vector<size_t> map(path.order());
+    path.walk(0, [&table](size_t coset, size_t gen) {
+        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];
+        }
+    }
+
+    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 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);
+
+    const auto table = context.sub(g_gens).solve({});
+
+    const tc::Cosets<> &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);
+    }, 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);
+        }
+    );
+    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 std::vector<size_t> &g_gens
+) {
+    if (g_gens.size() + 1 != N) {
+        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<> &,
+    const std::vector<size_t> &g_gens
+) {
+    if (not g_gens.empty()) {
+        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<std::vector<size_t>> &exclude) {
+    std::vector<std::vector<Primitive<N>>> parts;
+
+    auto g_gens = generators(group);
+
+    for (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 = each_tile(base, group, g_gens, sg_gens);
+        for (const auto &tile : tiles) {
+            parts.push_back(tile);
+        }
+    }
+    return parts;
+}
diff --git a/vis/src/main.cpp b/vis/src/main.cpp
index ae2c606..56230bb 100644
--- a/vis/src/main.cpp
+++ b/vis/src/main.cpp
@@ -7,7 +7,7 @@
 
 #include "util.hpp"
 #include "mirror.hpp"
-#include "geometry.hpp"
+#include "solver.hpp"
 
 #include <cgl/vertexarray.hpp>
 #include <cgl/shaderprogram.hpp>
@@ -50,7 +50,7 @@ Matrices build(GLFWwindow *window, State &state) {
     auto aspect = (float) width / (float) height;
     auto pheight = 1.4f;
     auto pwidth = aspect * pheight;
-    Eigen::Matrix4f proj = ortho(-pwidth, pwidth, -pheight, pheight, -10.0f, 10.0f);
+    Eigen::Matrix4f proj = orthographic(-pwidth, pwidth, -pheight, pheight, -10.0f, 10.0f);
 
     if (!glfwGetKey(window, GLFW_KEY_LEFT_SHIFT)) {
         state.st += state.time_delta / 8;
@@ -60,20 +60,20 @@ Matrices build(GLFWwindow *window, State &state) {
     view.setIdentity();
 
     if (state.dimension < 4) {
-        view *= utilRotate(2, 3, M_PI_2f32 + 0.01f);
+        view *= rot<4>(2, 3, M_PI_2f32 + 0.01f);
     }
 
     if (state.dimension > 1) {
-        view *= utilRotate(0, 1, state.st * .40f);
+        view *= rot<4>(0, 1, state.st * .40f);
     }
     if (state.dimension > 2) {
-        view *= utilRotate(0, 2, state.st * .20f);
-        view *= utilRotate(1, 2, state.st * .50f);
+        view *= rot<4>(0, 2, state.st * .20f);
+        view *= rot<4>(1, 2, state.st * .50f);
     }
     if (state.dimension > 3) {
-        view *= utilRotate(0, 3, state.st * 1.30f);
-        view *= utilRotate(1, 3, state.st * .25f);
-        view *= utilRotate(2, 3, state.st * 1.42f);
+        view *= rot<4>(0, 3, state.st * 1.30f);
+        view *= rot<4>(1, 3, state.st * .25f);
+        view *= rot<4>(2, 3, state.st * 1.42f);
     }
 
     return Matrices(proj, view);
@@ -179,7 +179,7 @@ struct SliceRenderer : public Renderer<N> {
 
     void _draw(const Prop<N> &prop) const override {
         glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, prop.vbo);
-        glProgramUniform3fv(solid, 2, 1, &prop.color.front());
+        glProgramUniform3fv(solid, 2, 1, prop.color.data());
 //        glProgramUniform3f(solid, 2, 1.f, 1.f, 1.f);
         prop.vao.bound([&]() {
             glDrawArrays(GL_POINTS, 0, prop.ibo.count() * N);
@@ -208,7 +208,7 @@ struct DirectRenderer : public Renderer<N> {
 
     void _draw(const Prop<N> &prop) const override {
         glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, prop.vbo);
-        glProgramUniform3fv(solid, 2, 1, &prop.color.front());
+        glProgramUniform3fv(solid, 2, 1, prop.color.data());
         prop.vao.bound([&]() {
             prop.ibo.bound(GL_ELEMENT_ARRAY_BUFFER, [&]() {
                 glDrawElements(GL_LINES, prop.ibo.count() * N, GL_UNSIGNED_INT, nullptr);
@@ -286,8 +286,14 @@ void run(const std::string &config_file, GLFWwindow *window) {
 
         if (group_info["slices"].IsDefined()) {
             for (const auto &slice_info : group_info["slices"]) {
-                auto root = slice_info["root"].as<vec5>();
-                auto color = slice_info["color"].as<vec3>();
+                auto root_arr = slice_info["root"].as<std::array<float, 5>>();
+                auto color_arr = slice_info["color"].as<std::array<float, 3>>();
+
+                vec5 root;
+                std::copy(root_arr.begin(), root_arr.end(), root.begin());
+                vec3 color;
+                std::copy(color_arr.begin(), color_arr.end(), color.begin());
+
                 auto exclude = std::vector<std::vector<size_t>>();
 
                 if (slice_info["exclude"].IsDefined()) {
@@ -300,7 +306,7 @@ void run(const std::string &config_file, GLFWwindow *window) {
                         group, root, color, subgroups, exclude
                     ));
                 } else {
-                    auto combos = Combos<size_t>(gens, 3);
+                    auto combos = combinations(gens, 3);
                     sRen.props.push_back(SliceProp<4>::build(
                         group, root, color, combos, exclude
                     ));
@@ -310,9 +316,16 @@ void run(const std::string &config_file, GLFWwindow *window) {
 
         if (group_info["wires"].IsDefined()) {
             for (const auto &wire_info : group_info["wires"]) {
-                auto root = wire_info["root"].as<vec5>();
-                auto color = wire_info["color"].as<vec3>();
+                auto root_arr = wire_info["root"].as<std::array<float, 5>>();
+                auto color_arr = wire_info["color"].as<std::array<float, 3>>();
+
+                vec5 root;
+                std::copy(root_arr.begin(), root_arr.end(), root.begin());
+                vec3 color;
+                std::copy(color_arr.begin(), color_arr.end(), color.begin());
+
                 auto exclude = std::vector<std::vector<size_t>>();
+
                 auto curve = wire_info["curve"].IsDefined() && wire_info["curve"].as<bool>();
                 auto ortho = wire_info["ortho"].IsDefined() && wire_info["ortho"].as<bool>();
 
@@ -341,7 +354,7 @@ void run(const std::string &config_file, GLFWwindow *window) {
                         ));
                     }
                 } else {
-                    auto combos = Combos<size_t>(gens, 1);
+                    auto combos = combinations(gens, 1);
 
                     if (ortho && curve) {
                         wocRen.props.push_back(WireframeProp::build(