new solver - point generation

vis/CMakeLists.txt

@@ -1,3 +1,9 @@
+add_custom_target(resources DEPENDS resources_output)
+add_custom_command(
+    OUTPUT resources_output
+    COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR}/res ${CMAKE_CURRENT_BINARY_DIR}/res
+    COMMENT "Copying Resources")
+
 add_executable(vis
     src/main.cpp
     )
@@ -10,3 +16,5 @@ target_link_libraries(vis
     eigen
     nlohmann_json
     )
+target_include_directories(vis PUBLIC include)
+add_dependencies(vis resources)

vis/include/geo/combo.hpp

@@ -0,0 +1,30 @@
+#pragma once
+
+#include <set>
+#include <algorithm>
+
+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<typename V>
+std::vector<V> combinations(const V &data, const size_t count) {
+    std::vector<V> result;
+
+    std::vector<bool> mask(data.size(), false);
+    std::fill(mask.begin(), mask.begin() + count, true);
+
+    do {
+        result.push_back(select(data, mask, count));
+    } while (std::next_permutation(mask.begin(), mask.end(), std::greater<>()));
+
+    return result;
+}

vis/include/geo/geometry.hpp

@@ -0,0 +1,63 @@
+#pragma once
+
+#include <tc/core.hpp>
+#include <cmath>
+#include <optional>
+#include <numeric>
+#include <iostream>
+
+#include <Eigen/Eigen>
+
+#include "combo.hpp"
+
+template<unsigned N>
+using Prims = Eigen::Matrix<unsigned, N, Eigen::Dynamic>;
+
+template<int N>
+using vec = Eigen::Matrix<float, N, 1>;
+template<int N>
+using mat = Eigen::Matrix<float, N, N>;
+
+using vec1 = vec<1>;
+using vec2 = vec<2>;
+using vec3 = vec<3>;
+using vec4 = vec<4>;
+using vec5 = vec<5>;
+
+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;
+}
+
+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;
+}
+
+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;
+}

vis/include/geo/mirror.hpp

@@ -0,0 +1,128 @@
+#pragma once
+
+#include <tc/core.hpp>
+
+#include <cmath>
+#include <vector>
+#include <algorithm>
+
+#include <geo/geometry.hpp>
+
+template<class V>
+float dot(int n, const V &a, const V &b) {
+    float sum = 0;
+    for (int i = 0; i < n; ++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;
+
+    for (int p = 0; p < group.rank; ++p) {
+        std::vector<float> vp;
+        for (int m = 0; m < p; ++m) {
+            auto &vq = mirrors[m];
+            vp.push_back((cos(M_PI / group.get(p, m)) - dot(m, vp, vq)) / vq[m]);
+        }
+        vp.push_back(std::sqrt(1 - dot(p, vp, vp)));
+
+        for (const auto &v : mirrors) {
+            if (dot(p, vp, vp) > 0) {
+                for (auto &e : vp) {
+                    e *= -1;
+                }
+                break;
+            }
+        }
+
+        mirrors.push_back(vp);
+    }
+
+    std::vector<vec<N>> res;
+    for (const auto &v : mirrors) {
+        vec<N> rv = vec<N>::Zero();
+
+        // ortho proj
+        for (int i = 0; i < std::min(v.size(), (size_t) N); ++i) {
+            rv[i] = v[i];
+        }
+
+        res.push_back(rv);
+    }
+    return res;
+}
+
+template<unsigned N>
+vec<N> stereo(const vec<N + 1> &v) {
+    vec<N> r;
+    for (int i = 0; i < N; ++i) {
+        r[i] = v[i] / (1 - v[N]);
+    }
+    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 vec.dot(target) / target.dot(target) * target;
+}
+
+template<class V>
+V reflect(const V &a, const V &axis) {
+    return a - 2.f * project(a, axis);
+}
+
+template<class V>
+V gram_schmidt_last(std::vector<V> vecs) {
+    for (int i = 0; i < vecs.size(); ++i) {
+        for (int j = 0; j < i; ++j) {
+            vecs[i] -= project(vecs[i], vecs[j]);
+        }
+    }
+
+    return vecs[vecs.size() - 1].normalized();
+}
+
+template<class V, class C>
+V barycentric(const std::vector<V> &basis, const C &coords) {
+    V res = V::Zero();
+
+    int N = std::min((int) basis.size(), (int) coords.rows());
+    for (int i = 0; i < N; ++i) {
+        res += basis[i] * coords[i];
+    }
+    return res;
+}
+
+template<class V>
+std::vector<V> plane_intersections(std::vector<V> normals) {
+    std::vector<V> results(normals.size());
+
+    for (int i = 0; i < normals.size(); ++i) {
+        std::rotate(normals.begin(), normals.begin() + 1, normals.end());
+        results[i] = gram_schmidt_last(normals);
+    }
+
+    return results;
+}
+
+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;
+}

vis/include/geo/solver.hpp

@@ -0,0 +1,233 @@
+#pragma once
+
+#include <tc/core.hpp>
+#include <cmath>
+#include <optional>
+#include <numeric>
+#include <iostream>
+
+#include <geo/geometry.hpp>
+
+#include "combo.hpp"
+
+/**
+ * Produce a list of all generators for the group context. The range [0..group.ngens).
+ */
+std::vector<int> generators(const tc::Group &context) {
+    // todo if tc::Group has 'global' generators, then this will be a member of tc::Group.
+    //  std::iota would populate a 'default' list of names, if names are not provided.
+    std::vector<int> g_gens(context.ngens);
+    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
+ */
+std::vector<int> recontext_gens(
+    const tc::Group &context,
+    std::vector<int> g_gens,
+    std::vector<int> sg_gens) {
+    // todo ideally tc::Group will deal in 'global' generators so this stell will be unecessary.
+
+    std::sort(g_gens.begin(), g_gens.end());
+
+    int inv_gen_map[context.ngens];
+    for (size_t i = 0; i < g_gens.size(); i++) {
+        inv_gen_map[g_gens[i]] = i;
+    }
+
+    std::vector<int> s_sg_gens;
+    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;
+}
+
+/**
+ * Solve the cosets generated by sg_gens within the subgroup generated by g_gens of the group context
+ */
+tc::Cosets solve(
+    const tc::Group &context,
+    const std::vector<int> &g_gens,
+    const std::vector<int> &sg_gens
+) {
+    // todo this should also be handled with 'global' generators.
+    const auto proper_sg_gens = recontext_gens(context, g_gens, sg_gens);
+    return context.subgroup(g_gens).solve(proper_sg_gens);
+}
+
+/**
+ * Apply some context transformation to all primitives of this mesh.
+ */
+template<unsigned N>
+void apply(const tc::Cosets &table, int gen, Prims<N> &mat) {
+    auto data = mat.data();
+    for (int i = 0; i < mat.size(); ++i) {
+        data[i] = table.get(data[i], gen);
+    }
+}
+
+/**
+ * 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]]
+Prims<N> recontext(
+    Prims<N> prims,
+    const tc::Group &context,
+    const std::vector<int> &g_gens,
+    const std::vector<int> &sg_gens
+) {
+    // todo this will be simpler with 'global' gens, but it's still not free...
+    const auto proper_sg_gens = recontext_gens(context, g_gens, sg_gens);
+    const auto table = solve(context, g_gens, {});
+    const auto path = solve(context, sg_gens, {}).path;
+
+    auto map = path.template walk<int, int>(0, proper_sg_gens, [table](int coset, int gen) {
+        return table.get(coset, gen);
+    });
+
+    Prims<N> res(prims);
+    auto data = res.data();
+    for (int i = 0; i < prims.size(); ++i) {
+        data[i] = map[data[i]];
+    }
+
+    return res;
+}
+
+/**
+ * Union several meshes of the same dimension
+ */
+template<unsigned N>
+Prims<N> merge(const std::vector<Prims<N>> &meshes) {
+    // todo (?) might be possible with NullaryExpr
+    size_t cols = 0;
+    for (const auto &mesh: meshes) {
+        cols += mesh.cols();
+    }
+
+    Prims<N> res(N, cols);
+
+    size_t offset = 0;
+    for (const Prims<N> &mesh: meshes) {
+        res.middleCols(offset, mesh.cols()) = mesh;
+        offset += mesh.cols();
+    }
+
+    return res;
+}
+
+template<unsigned N>
+[[nodiscard]]
+std::vector<Prims<N>> tile(
+    Prims<N> prims,
+    const tc::Group &context,
+    const std::vector<int> &g_gens,
+    const std::vector<int> &sg_gens
+) {
+    // todo convert to nullaryexpr.
+    //  some stuff will be easier with global generators, but not all.
+    Prims<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 = solve(context, g_gens, {});
+    const auto path = solve(context, g_gens, sg_gens).path;
+
+    std::vector<int> _gens = generators(context);
+
+    std::vector<Prims<N>> res = path.walk<Prims<N>, int>(
+        base, _gens,
+        [&](Prims<N> from, int gen) {
+            apply<N>(table, gen, from);
+            return from;
+        }
+    );
+
+    return res;
+}
+
+/**
+ * Produce a mesh of higher dimension by fanning a single point to all primitives in this mesh.
+ */
+template<unsigned N>
+[[nodiscard]]
+Prims<N + 1> fan(Prims<N> prims, int root) {
+    // todo convert to nullaryexpr.
+    Prims<N + 1> res(N + 1, prims.cols());
+
+    res.topRows(1) = Prims<1>::Constant(1, prims.cols(), root);
+    res.bottomRows(N) = prims;
+
+    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>
+Prims<N> triangulate(
+    const tc::Group &context,
+    const std::vector<int> &g_gens
+) {
+    // todo (?) might be possible with nullaryexpr
+    //  not so sure, though.
+    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<Prims<N>> meshes;
+
+    for (const auto &sg_gens: combos) {
+        auto base = triangulate<N - 1>(context, sg_gens);
+        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<N>(meshes);
+}
+
+/**
+ * Single-index primitives should not be further triangulated.
+ */
+template<>
+Prims<1> triangulate<1>(
+    const tc::Group &context,
+    const std::vector<int> &g_gens
+) {
+    if (not g_gens.empty()) // todo make static assert
+        throw std::logic_error("g_gens must be empty for a trivial Mesh");
+
+    return Prims<1>::Zero(1, 1);
+}
+
+template<unsigned N, class T>
+auto hull(const tc::Group &group, T all_sg_gens, const std::vector<std::vector<int>> &exclude) {
+    std::vector<Prims<N>> parts;
+    auto g_gens = generators(group);
+    for (const std::vector<int> &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 = tile<N>(base, group, g_gens, sg_gens);
+        for (const auto &tile: tiles) {
+            parts.push_back(tile);
+        }
+    }
+    return parts;
+}

vis/include/gl/buffer.hpp

@@ -0,0 +1,53 @@
+#pragma once
+
+#include <glad/glad.h>
+
+#include <utility>
+
+template<typename T_>
+class Buffer {
+public:
+    using Type = T_;
+
+private:
+    GLuint id = 0;
+
+public:
+    Buffer() {
+        glCreateBuffers(1, &id);
+    }
+
+    Buffer(Buffer&& o) noexcept {
+        id = std::exchange(o.id, 0);
+    }
+
+    Buffer(const Buffer&) = delete;  // this is doable, but would be slow.
+
+    operator GLuint() const { // NOLINT(google-explicit-constructor)
+        return id;
+    }
+
+    template<typename RandomIt>
+    GLuint upload(RandomIt first, RandomIt last, GLenum mode = GL_STATIC_DRAW) {
+        size_t count = last - first;
+
+        // todo StaticBuffer that uses BufferStorage
+        glNamedBufferData(id, sizeof(Type) * count, nullptr, mode);
+
+        Type* out = (Type*) glMapNamedBuffer(id, GL_WRITE_ONLY);
+        std::copy(first, last, out);
+        glUnmapNamedBuffer(id);
+
+        return count;
+    }
+
+    template<typename T>
+    GLuint upload(const T& data, GLenum mode = GL_STATIC_DRAW) {
+        return upload(data.begin(), data.end(), mode);
+    }
+
+    ~Buffer() {
+        // delete silently ignores 0.
+        glDeleteBuffers(1, &id);
+    }
+};

vis/include/gl/debug.hpp

@@ -1,5 +1,10 @@
 #pragma once
 
+#include <glad/glad.h>
+#include <fmt/core.h>
+
+#include <string>
+
 #ifndef NDEBUG
 
 void GLAPIENTRY log_gl_debug_callback(

vis/include/gl/shader.hpp

@@ -0,0 +1,79 @@
+#pragma once
+
+#include <glad/glad.h>
+
+#include <string>
+#include <memory>
+#include <vector>
+
+template<GLenum mode>
+class Shader {
+public:
+private:
+    GLuint id = 0;
+
+public:
+    explicit Shader(const std::string &src) {
+        id = glCreateShader(mode);
+
+        const char *str = src.c_str();
+        glShaderSource(id, 1, &str, nullptr);
+        glCompileShader(id);
+
+        // todo throw if compile failed
+    }
+
+    explicit Shader(std::ifstream source)
+        : Shader(std::string(
+        std::istreambuf_iterator<char>(source),
+        std::istreambuf_iterator<char>()
+    )) {}
+
+    Shader(const Shader &) = delete;
+
+    Shader(Shader &&o) noexcept {
+        id = std::exchange(o.id, 0);
+    }
+
+    ~Shader() {
+        glDeleteShader(id);
+    }
+
+    operator GLuint() const { // NOLINT(google-explicit-constructor)
+        return id;
+    }
+};
+
+using VertexShader = Shader<GL_VERTEX_SHADER>;
+using FragmentShader = Shader<GL_FRAGMENT_SHADER>;
+
+class Program {
+private:
+    GLuint id = 0;
+
+public:
+    template<GLenum ...mode>
+    explicit Program(const Shader<mode> &...shader) {
+        id = glCreateProgram();
+
+        (glAttachShader(id, shader), ...);
+
+        glLinkProgram(id);
+
+        // todo throw if link failed
+    }
+
+    Program(const Program &) = delete;
+
+    Program(Program &&o) noexcept {
+        id = std::exchange(o.id, 0);
+    }
+
+    ~Program() {
+        glDeleteProgram(id);
+    }
+
+    operator GLuint() const { // NOLINT(google-explicit-constructor)
+        return id;
+    }
+};

vis/include/gl/types.hpp

@@ -0,0 +1,118 @@
+#pragma once
+
+#include <glad/glad.h>
+
+#include <Eigen/Eigen>
+
+template<typename T_>
+struct Buffer;
+
+template<class ...Fmt_>
+struct VertexArray;
+
+
+template<typename ctype, GLenum type, GLint size, GLuint offset = 0, GLboolean norm = GL_FALSE>
+struct Format {
+    using CType = ctype;
+    static constexpr GLenum Type = type;
+    static constexpr GLint Size = size;
+    static constexpr GLuint Offset = offset;
+    static constexpr GLboolean Norm = norm;
+
+    template<typename ctype_, GLint size_>
+    using As = Format<ctype_, type, size_, offset, norm>;
+
+    static inline constexpr void apply(GLuint vao, GLuint idx) {
+        glEnableVertexArrayAttrib(vao, idx);
+        glVertexArrayAttribFormat(vao, idx, size, type, norm, offset);
+    }
+};
+
+template<typename ctype, GLenum type, GLint size, GLuint offset = 0>
+struct IFormat {
+    using CType = ctype;
+    static constexpr GLenum Type = type;
+    static constexpr GLint Size = size;
+    static constexpr GLuint Offset = offset;
+
+    template<typename ctype_, GLint size_>
+    using As = IFormat<ctype_, type, size_, offset>;
+
+    static inline constexpr void apply(GLuint vao, GLuint idx) {
+        glEnableVertexArrayAttrib(vao, idx);
+        glVertexArrayAttribIFormat(vao, idx, size, type, offset);
+    }
+};
+
+template<typename ctype, GLenum type, GLint size, GLuint offset = 0>
+struct LFormat {
+    using CType = ctype;
+    static constexpr GLenum Type = type;
+    static constexpr GLint Size = size;
+    static constexpr GLuint Offset = offset;
+
+    template<typename ctype_, GLint size_>
+    using As = LFormat<ctype_, type, size_, offset>;
+
+    static inline constexpr void apply(GLuint vao, GLuint idx) {
+        glEnableVertexArrayAttrib(vao, idx);
+        glVertexArrayAttribLFormat(vao, idx, size, type, offset);
+    }
+};
+
+template<typename Fmt_>
+struct AutoFormat {
+    using Fmt = Fmt_;
+};
+
+template<>
+struct AutoFormat<GLdouble> {
+    using Fmt = LFormat<GLdouble, GL_DOUBLE, 1>;
+};
+
+template<>
+struct AutoFormat<GLfloat> {
+    using Fmt = Format<GLfloat, GL_FLOAT, 1>;
+};
+
+template<>
+struct AutoFormat<GLint> {
+    using Fmt = IFormat<GLint, GL_INT, 1>;
+};
+
+template<>
+struct AutoFormat<GLuint> {
+    using Fmt = IFormat<GLuint, GL_UNSIGNED_INT, 1>;
+};
+
+template<typename Scalar, int Dim>
+struct AutoFormat<Eigen::Vector<Scalar, Dim>> {
+    static_assert((Dim >= 1) && (Dim <= 4), "Dim must be in range 1-4");
+
+    using Fmt = typename AutoFormat<Scalar>::Fmt::template As<Eigen::Vector<Scalar, Dim>, Dim>;
+};
+
+template<typename Fmt_>
+struct Binder {
+    using CType = typename AutoFormat<Fmt_>::Fmt::CType;
+
+    const GLuint buf;
+    const GLuint offset;
+    const GLuint stride;
+
+    Binder(const Buffer<CType> &buf) // NOLINT(google-explicit-constructor)
+        : buf((GLuint) buf), offset(0), stride(sizeof(CType)) {}
+
+    template<typename CType_>
+    Binder(const Buffer<CType_> &buf, GLuint offset)
+        : buf((GLuint) buf), offset(offset), stride(sizeof(CType_)) {}
+
+    Binder(const Binder<Fmt_> &o)
+        : buf(o.buf), offset(o.offset), stride(o.stride) {}
+
+    inline void bind(GLuint vao, GLuint idx) const {
+        glVertexArrayVertexBuffer(vao, idx, buf, offset, stride);
+    }
+};
+
+#define ATTR(buf, field) Binder<decltype(decltype(buf)::Type::field)>(buf, offsetof(decltype(buf)::Type, field))

vis/include/gl/vertexarray.hpp

@@ -0,0 +1,67 @@
+#pragma once
+
+#include <glad/glad.h>
+
+#include <utility>
+
+#include "buffer.hpp"
+#include "types.hpp"
+
+template<typename ...Fmt_>
+class VertexArray {
+public:
+    template<size_t idx>
+    using Fmt = std::tuple_element_t<idx, std::tuple<Fmt_...>>;
+
+private:
+    GLuint id = 0;
+
+    template<size_t ...idx>
+    inline void formatall(
+        std::integer_sequence<size_t, idx...>
+    ) {
+        (AutoFormat<Fmt_>::Fmt::apply(id, idx), ...);
+    }
+
+    template<size_t ...idx>
+    inline void bindall(
+        const Binder<Fmt_> &...buf,
+        std::integer_sequence<size_t, idx...>
+    ) {
+        (buf.bind(id, idx), ...);
+    }
+
+public:
+    explicit VertexArray() {
+        glCreateVertexArrays(1, &id);
+        formatall(std::make_index_sequence<sizeof...(Fmt_)>());
+    }
+
+    explicit VertexArray(const Binder<Fmt_> &...buf)
+        : VertexArray() {
+        bind(buf...);
+    }
+
+    VertexArray(VertexArray &&o) noexcept {
+        id = std::exchange(o.id, 0);
+    }
+
+    VertexArray(const VertexArray &) = delete;  // this is doable, but would be slow.
+
+    ~VertexArray() {
+        glDeleteVertexArrays(1, &id);
+    }
+
+    operator GLuint() const { // NOLINT(google-explicit-constructor)
+        return id;
+    }
+
+    void bind(const Binder<Fmt_> &...buf) {
+        bindall(buf..., std::make_index_sequence<sizeof...(Fmt_)>());
+    }
+
+    template<size_t idx>
+    void bind(const Binder<Fmt<idx>> &buf) {
+        buf.bind(id, idx);
+    }
+};

vis/res/shaders/4d.vert.glsl

@@ -0,0 +1,12 @@
+#version 440
+
+layout(location=1) uniform float time;
+layout(location=2) uniform mat4 proj;
+layout(location=3) uniform mat4 rot;
+
+layout(location=0) in vec4 pos;
+
+void main() {
+    vec4 pos = rot * pos;
+    gl_Position = proj * vec4(pos.xyz / (1 - pos.w), 1.0);
+}

vis/res/shaders/main.frag.glsl

@@ -0,0 +1,12 @@
+#version 440
+
+layout(location=0) uniform vec4 ucol;
+
+layout(location=0) out vec4 col;
+
+void main() {
+    float d = 1.0 - gl_FragCoord.z;
+    d = (d - 0.5) / 0.7 + 0.5;
+    col = ucol;
+    col.xyz *= d;
+}

vis/res/shaders/main.vert.glsl

@@ -0,0 +1,15 @@
+#version 440
+
+//layout(location=1) uniform float time;
+layout(location=2) uniform mat4 proj;
+//layout(location=3) uniform mat4 rot;
+
+layout(location=0) in vec4 pos;
+
+void main() {
+//    mat3 rot3 = mat3(rot);
+//    vec3 pos3 = rot3 * pos;
+
+    vec3 pos3 = pos.xyz;
+    gl_Position = proj * vec4(pos3, 1.0);
+}

vis/src/combotest.cpp

@@ -0,0 +1,22 @@
+#include <combo.hpp>
+
+#include <iostream>
+
+std::ostream &operator<<(std::ostream &o, const std::vector<int> &data) {
+    o << "[ ";
+    for (const auto &el: data) {
+        o << el << " ";
+    }
+    o << "]";
+    return o;
+}
+
+int main() {
+    std::vector<int> data{1, 2, 3, 4, 5};
+
+    for (const auto &combo: combinations(data, 3)) {
+        std::cout << combo << std::endl;
+    }
+
+    return EXIT_SUCCESS;
+}

vis/src/geometrytest.cpp

@@ -0,0 +1,74 @@
+#include <tc/groups.hpp>
+#include <geo/geometry.hpp>
+#include <geo/solver.hpp>
+#include <geo/mirror.hpp>
+
+#include <ml/meshlib.hpp>
+#include <ml/meshlib_json.hpp>
+
+#include <fstream>
+#include <iostream>
+
+template<int N>
+using ProjectiveNf = Eigen::Transform<float, N, Eigen::Projective>;
+
+template<int N>
+using VectorNf = Eigen::Vector<float, N>;
+
+template<int N>
+Eigen::Matrix<float, N, Eigen::Dynamic> make_points(
+    const tc::Group &group,
+    const Eigen::Vector<float, N> &root
+) {
+    // todo clean up mirror / plane_intersections / barycentric
+    // ideally barycentric will work in rotors, so that stellations etc. will be possible
+    auto mirrors = mirror<N>(group);
+    auto corners = plane_intersections(mirrors);
+    auto start = barycentric(corners, root);
+
+    auto cosets = group.solve();
+
+    auto verts = cosets.path.walk<VectorNf<N>, VectorNf<N>>(start, mirrors, reflect<vec4>);
+
+    Eigen::Matrix<float, N, Eigen::Dynamic> points(root.size(), verts.size());
+    std::copy(verts.begin(), verts.end(), points.colwise().begin());
+
+    return points;
+}
+
+template<int N>
+Eigen::Matrix<unsigned, N, Eigen::Dynamic> make_cells(
+    const tc::Group &group,
+    const std::vector<std::vector<int>> &exclude = {}
+) {
+    auto gens = generators(group);
+    auto combos = combinations(gens, N - 1);
+
+    // todo clean up merge(hull(...))
+    Eigen::Matrix<unsigned, N, Eigen::Dynamic> cells = merge<N>(hull<N>(group, combos, exclude));
+
+    return cells;
+}
+
+int main() {
+    std::vector<int> symbol = {5, 3, 3};
+    auto group = tc::schlafli(symbol);
+
+    vec4 root{0.80, 0.02, 0.02, 0.02};
+
+    auto points = make_points(group, root);
+    auto cells = make_cells<3>(group, {{0, 1}});
+
+    ml::Mesh mesh(points, cells);
+
+//    auto transform = ProjectiveNf<vec4::RowsAtCompileTime>::Identity();
+//    transform.translation() << 0, 0.5, 0, 0;
+//    points = (transform * points.colwise().homogeneous()).colwise().hnormalized();
+
+    std::cout << points.rows() << " " << points.cols() << std::endl;
+    std::cout << cells.rows() << " " << cells.cols() << std::endl;
+
+    ml::write(mesh, std::ofstream("dodeca.pak", std::ios::out | std::ios::binary));
+
+    return EXIT_SUCCESS;
+}

vis/src/main.cpp

@@ -3,28 +3,193 @@
 #include <imgui.h>
 #include <backends/imgui_impl_glfw.h>
 #include <backends/imgui_impl_opengl3.h>
-
-#include <fmt/core.h>
+#include <iostream>
 #include <fstream>
 
-#include "debug.hpp"
+#include <fmt/core.h>
+#include <fmt/ranges.h>
 
-int run(GLFWwindow *window, ImGuiContext *ctx) {
-    glClearColor(0.1, 0.1, 0.9, 1.0);
+#include <gl/debug.hpp>
+#include <gl/buffer.hpp>
+#include <gl/shader.hpp>
+#include <gl/vertexarray.hpp>
+
+#include <tc/groups.hpp>
+#include <tc/core.hpp>
+
+#include <geo/mirror.hpp>
+
+struct State {
+    Eigen::Vector4f bg{0.169f, 0.169f, 0.169f, 1.00f};
+    Eigen::Vector4f fg{0.71f, 0.53f, 0.94f, 1.00f};
+    Eigen::Vector4f wf{0.95f, 0.95f, 0.95f, 1.00f};
+
+    Eigen::Vector4f R{1.00f, 0.00f, 0.00f, 1.00f};
+    Eigen::Vector4f G{0.00f, 1.00f, 0.00f, 1.00f};
+    Eigen::Vector4f B{0.00f, 0.00f, 1.00f, 1.00f};
+    Eigen::Vector4f Y{1.20f, 1.20f, 0.00f, 1.00f};
+
+    Eigen::Matrix4f rot = Eigen::Matrix4f::Identity();
+
+    bool color_axes = false;
+};
+
+Eigen::Matrix4f rotor(int u, int v, float rad) {
+    Eigen::Matrix4f res = Eigen::Matrix4f::Identity();
+    res(u, u) = res(v, v) = cosf(rad);
+    res(u, v) = res(v, u) = sinf(rad);
+    res(u, v) *= -1;
+    return res;
+}
+
+template<typename T_>
+T_ mix(const T_ &a, const T_ &b, const typename T_::Scalar &x) {
+    return a * (1 - x) + b * x;
+}
+
+void show_overlay(State &state) {
+    static std::string gl_vendor = (const char *) glGetString(GL_VENDOR);
+    static std::string gl_renderer = (const char *) glGetString(GL_RENDERER);
+    static std::string gl_version = (const char *) glGetString(GL_VERSION);
+    static std::string glsl_version = (const char *) glGetString(GL_SHADING_LANGUAGE_VERSION);
+
+    ImGuiWindowFlags window_flags =
+        ImGuiWindowFlags_AlwaysAutoResize |
+        ImGuiWindowFlags_NoSavedSettings |
+        ImGuiWindowFlags_NoFocusOnAppearing |
+        ImGuiWindowFlags_NoNav |
+        ImGuiWindowFlags_NoBringToFrontOnFocus |
+        ImGuiWindowFlags_NoMove;
+
+    ImGuiStyle &style = ImGui::GetStyle();
+    const auto PAD = style.DisplaySafeAreaPadding;
+    auto window_pos = PAD;
+
+    ImGui::SetNextWindowPos(window_pos, ImGuiCond_Always);
+    ImGui::SetNextWindowBgAlpha(0.35f * style.Alpha);
+    ImGui::SetNextWindowCollapsed(true, ImGuiCond_Appearing);
+    ImGui::Begin("Graphics Information", nullptr, window_flags);
+    ImGui::Text("GL Vendor    | %s", gl_vendor.c_str());
+    ImGui::Text("GL Renderer  | %s", gl_renderer.c_str());
+    ImGui::Text("GL Version   | %s", gl_version.c_str());
+    ImGui::Text("GLSL Version | %s", glsl_version.c_str());
+
+    auto v2 = ImGui::GetWindowSize();
+    window_pos.y += v2.y + PAD.y;
+    ImGui::End();
+
+    ImGui::SetNextWindowPos(window_pos, ImGuiCond_Always);
+    ImGui::SetNextWindowBgAlpha(0.35f * style.Alpha);
+    ImGui::Begin("Controls", nullptr, window_flags);
+
+    ImGuiIO &io = ImGui::GetIO();
+    ImGui::Text("FPS          | %.2f", io.Framerate);
+
+    ImGui::Separator();
+
+    ImGui::ColorEdit3("Background", state.bg.data(), ImGuiColorEditFlags_Float);
+    ImGui::ColorEdit3("Foreground", state.fg.data(), ImGuiColorEditFlags_Float);
+    ImGui::ColorEdit3("Wireframe", state.wf.data(), ImGuiColorEditFlags_Float);
+
+    ImGui::Checkbox("Show RGBY axis colors", &state.color_axes);
+
+    if (io.MouseDown[0] && !io.WantCaptureMouse) {
+        Eigen::Matrix4f rot = Eigen::Matrix4f::Identity();
+        Eigen::Vector2f del{io.MouseDelta.x, io.MouseDelta.y};
+        del /= 200.0f;
+
+        if (io.KeyShift) {
+            del /= 5.0f;
+        }
+
+        if (io.KeyCtrl) {
+            Eigen::Matrix4f rx = rotor(0, 3, -del.x());
+            Eigen::Matrix4f ry = rotor(1, 3, del.y());
+            rot = rx * ry;
+        } else {
+            Eigen::Matrix4f rx = rotor(0, 2, -del.x());
+            Eigen::Matrix4f ry = rotor(1, 2, del.y());
+            rot = rx * ry;
+        }
+
+        state.rot = rot * state.rot;
+    }
+
+    ImGui::End();
+}
+
+void set_style() {
+    ImGui::StyleColorsDark();
+
+    ImGuiStyle &style = ImGui::GetStyle();
+    style.WindowRounding = 4;
+    style.FrameRounding = 2;
+    style.DisplaySafeAreaPadding.x = 10;
+    style.DisplaySafeAreaPadding.y = 10;
+}
+
+int run(GLFWwindow *window, ImGuiContext *context) {
+    State state;
+
+//    Buffer<GLuint> ind_buf;
+    Buffer<Eigen::Vector4f> vert_buf;
+
+    VertexArray<Eigen::Vector4f> vao(vert_buf);
+//    glVertexArrayElementBuffer(vao, ind_buf);
+
+    tc::Group group = tc::coxeter("3 4 3");
+
+    auto cosets = solve(group, {}, 1000000);
+    vec4 coords {1, 1, 1, 1};
+
+    auto mirrors = mirror<4>(group);
+
+    auto corners = plane_intersections(mirrors);
+    auto start = barycentric(corners, coords);
+    start.normalize();
+
+    auto points = cosets.path.walk<vec4, vec4>(start, mirrors, reflect<vec4>);
+
+    vert_buf.upload(points);
+
+    VertexShader vs(std::ifstream("res/shaders/main.vert.glsl"));
+    FragmentShader fs(std::ifstream("res/shaders/main.frag.glsl"));
+
+    Program pgm(vs, fs);
+
+    glEnable(GL_DEPTH_TEST);
+    glPointSize(2);
+
+    Eigen::Projective3f proj;
 
     while (!glfwWindowShouldClose(window)) {
         glfwPollEvents();
 
-        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
-        
         ImGui_ImplOpenGL3_NewFrame();
         ImGui_ImplGlfw_NewFrame();
         ImGui::NewFrame();
-        ImGui::Begin("Hello There!", nullptr, ImGuiWindowFlags_None);
-        ImGui::Text("General Kenobi.");
-        ImGui::End();
+        show_overlay(state);
         ImGui::Render();
 
+        int display_w, display_h;
+        glfwGetFramebufferSize(window, &display_w, &display_h);
+        glViewport(0, 0, display_w, display_h);
+        glClearColor(state.bg[0], state.bg[1], state.bg[2], state.bg[3]);
+        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+        auto aspect = (float) display_h / (float) display_w;
+        proj = Eigen::AlignedScaling3f(aspect, 1.0, -0.6);
+
+        glUseProgram(pgm);
+        glBindVertexArray(vao);
+        glUniform4fv(0, 1, state.fg.data());
+        glUniform1f(1, (GLfloat) glfwGetTime());
+        glUniformMatrix4fv(2, 1, false, proj.data());
+        glUniformMatrix4fv(3, 1, false, state.rot.data());
+        glDrawArrays(GL_POINTS, 0, points.size());
+        glBindVertexArray(0);
+        glUseProgram(0);
+
         ImGui_ImplOpenGL3_RenderDrawData(ImGui::GetDrawData());
         glfwSwapBuffers(window);
     }
@@ -34,21 +199,16 @@ int run(GLFWwindow *window, ImGuiContext *ctx) {
 
 int main() {
     if (!glfwInit()) {
-        fmt::print(stderr, "GLFW : Failed Initialization\n");
+        std::cerr << "GLFW:Failed initialization" << std::endl;
         return EXIT_FAILURE;
     }
 
     glfwWindowHint(GLFW_CONTEXT_VERSION_MAJOR, 4);
     glfwWindowHint(GLFW_CONTEXT_VERSION_MINOR, 6);
     glfwWindowHint(GLFW_OPENGL_PROFILE, GLFW_OPENGL_CORE_PROFILE);
-    GLFWwindow *window = glfwCreateWindow(
-        1280, 720,
-        "Cosets Visualization",
-        nullptr, nullptr
-    );
-
+    auto *window = glfwCreateWindow(1280, 720, "Cosets Visualization", nullptr, nullptr);
     if (!window) {
-        fmt::print("GLFW : Failed to create window.\n");
+        std::cerr << "GLFW:Failed to create window" << std::endl;
         return EXIT_FAILURE;
     }
 
@@ -67,17 +227,19 @@ int main() {
 #endif
 
     IMGUI_CHECKVERSION();
-    ImGuiContext *ctx = ImGui::CreateContext();
+    auto *context = ImGui::CreateContext();
     ImGui_ImplGlfw_InitForOpenGL(window, true);
     ImGui_ImplOpenGL3_Init("#version 130");
 
-    int code = EXIT_SUCCESS;
+    set_style();
+
+    int exit_code = EXIT_SUCCESS;
 
     try {
-        code = run(window, ctx);
-    } catch (const std::exception &ex) {
-        fmt::print(stderr, "{}\n", ex.what());
-        code = EXIT_FAILURE;
+        exit_code = run(window, context);
+    } catch (const std::exception &e) {
+        std::cerr << e.what() << std::endl;
+        exit_code = EXIT_FAILURE;
     }
 
     ImGui_ImplOpenGL3_Shutdown();
@@ -87,5 +249,5 @@ int main() {
     glfwDestroyWindow(window);
     glfwTerminate();
 
-    return code;
+    return exit_code;
 }

vis/src/serialtest.cpp

@@ -0,0 +1,45 @@
+#include <Eigen/Eigen>
+
+#include <fstream>
+#include <iostream>
+
+#include <ml/meshlib.hpp>
+#include <ml/meshlib_json.hpp>
+
+auto make_circle(float radius, size_t npoints = 32) {
+    Eigen::Array<float, 1, Eigen::Dynamic> theta(1, npoints);
+    for (int i = 0; i < theta.size(); ++i) {
+        theta(i) = (float) M_PI * 2.0f * (float) i / (float) npoints;
+    }
+
+    Eigen::Array<float, 3, Eigen::Dynamic> points(3, npoints);
+    points.row(0) = theta.sin();
+    points.row(1) = theta.cos();
+    points.row(2).setZero();
+
+    Eigen::Array<unsigned int, 1, Eigen::Dynamic> idx(1, npoints - 1);
+    for (int i = 0; i < idx.size(); ++i) {
+        idx(i) = i;
+    }
+
+    Eigen::Array<unsigned int, 3, Eigen::Dynamic> cells(3, npoints - 1);
+    cells.row(0) = 0;
+    cells.row(1) = idx;
+    cells.row(2) = idx + 1;
+
+    return ml::Mesh(points, cells);
+}
+
+int main() {
+    std::string path = "circle.pak";
+
+    auto omesh = make_circle(1.0f);
+    using MT = decltype(omesh);
+
+    ml::write(omesh, std::ofstream(path, std::ios::out | std::ios::binary));
+
+    auto imesh = ml::read<MT>(std::ifstream(path, std::ios::in | std::ios::binary));
+
+    std::cout << imesh.points << std::endl;
+    std::cout << imesh.cells << std::endl;
+}