toddcox-visualize/include/solver.hpp

#pragma once

#include <tc/core.hpp>
#include <cmath>
#include <optional>
#include <numeric>
#include <iostream>

#include <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) {
    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) {

    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
) {
    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
) {
    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) {
    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
) {
    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) {
    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
) {
    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;
}