#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <limits>
#include <sstream>
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <vector>

namespace {

constexpr double kPi = 3.14159265358979323846;

struct Point {
    double x{};
    double y{};
};

struct Affine {
    double a{1.0}, b{0.0}, c{0.0};
    double d{0.0}, e{1.0}, f{0.0};
};

struct Tile {
    std::string label;
    std::vector<Point> points;
    std::vector<int> neighbors;
    bool alive{false};
    bool nextAlive{false};
    int age{0};
};

struct Config {
    int depth{3};
    int birthFaces{3};
    int surviveMin{2};
    int surviveMax{3};
    int faceWeight{1};
    int steps{120};
    std::string root{"Delta"};
    double density{0.42};
    std::uint64_t seed{1};
    std::string svgPath;
};

struct Bounds {
    double minX{0.0};
    double minY{0.0};
    double maxX{0.0};
    double maxY{0.0};
};

using RuleMap = std::unordered_map<std::string, std::array<std::string, 8>>;

Point pt(double x, double y) { return Point{x, y}; }

Affine mul(const Affine& A, const Affine& B) {
    return {
        A.a * B.a + A.b * B.d,
        A.a * B.b + A.b * B.e,
        A.a * B.c + A.b * B.f + A.c,
        A.d * B.a + A.e * B.d,
        A.d * B.b + A.e * B.e,
        A.d * B.c + A.e * B.f + A.f
    };
}

Point apply(const Affine& M, const Point& p) {
    return {M.a * p.x + M.b * p.y + M.c, M.d * p.x + M.e * p.y + M.f};
}

Affine rotation(double angle) {
    const double c = std::cos(angle);
    const double s = std::sin(angle);
    return {c, -s, 0.0, s, c, 0.0};
}

Affine translation(double tx, double ty) {
    return {1.0, 0.0, tx, 0.0, 1.0, ty};
}

double radians(double degrees) { return degrees * kPi / 180.0; }

Affine transTo(const Point& from, const Point& to) {
    return translation(to.x - from.x, to.y - from.y);
}

std::uint64_t mix64(std::uint64_t x) {
    x += 0x9e3779b97f4a7c15ULL;
    x = (x ^ (x >> 30U)) * 0xbf58476d1ce4e5b9ULL;
    x = (x ^ (x >> 27U)) * 0x94d049bb133111ebULL;
    return x ^ (x >> 31U);
}

double randomUnit(std::uint64_t value) {
    return static_cast<double>(mix64(value) >> 11U) * (1.0 / 9007199254740992.0);
}

std::vector<Point> baseSpectre() {
    return {
        pt(0.0, 0.0),
        pt(1.0, 0.0),
        pt(1.5, -0.8660254037844386),
        pt(2.366025403784439, -0.36602540378443865),
        pt(2.366025403784439, 0.6339745962155614),
        pt(3.366025403784439, 0.6339745962155614),
        pt(3.866025403784439, 1.5),
        pt(3.0, 2.0),
        pt(2.133974596215561, 1.5),
        pt(1.6339745962155614, 2.3660254037844393),
        pt(0.6339745962155614, 2.3660254037844393),
        pt(-0.3660254037844386, 2.3660254037844393),
        pt(-0.866025403784439, 1.5),
        pt(0.0, 1.0)
    };
}

std::array<Point, 4> baseQuad(const std::vector<Point>& spectre) {
    return {spectre[3], spectre[5], spectre[7], spectre[11]};
}

RuleMap superRules() {
    return {
        {"Gamma", {"Pi", "Delta", "null", "Theta", "Sigma", "Xi", "Phi", "Gamma"}},
        {"Delta", {"Xi", "Delta", "Xi", "Phi", "Sigma", "Pi", "Phi", "Gamma"}},
        {"Theta", {"Psi", "Delta", "Pi", "Phi", "Sigma", "Pi", "Phi", "Gamma"}},
        {"Lambda", {"Psi", "Delta", "Xi", "Phi", "Sigma", "Pi", "Phi", "Gamma"}},
        {"Xi", {"Psi", "Delta", "Pi", "Phi", "Sigma", "Psi", "Phi", "Gamma"}},
        {"Pi", {"Psi", "Delta", "Xi", "Phi", "Sigma", "Psi", "Phi", "Gamma"}},
        {"Sigma", {"Xi", "Delta", "Xi", "Phi", "Sigma", "Pi", "Lambda", "Gamma"}},
        {"Phi", {"Psi", "Delta", "Psi", "Phi", "Sigma", "Pi", "Phi", "Gamma"}},
        {"Psi", {"Psi", "Delta", "Psi", "Phi", "Sigma", "Psi", "Phi", "Gamma"}}
    };
}

std::vector<Affine> childTransforms(const std::array<Point, 4>& quad) {
    const std::array<std::array<int, 3>, 7> tRules{{
        {60, 3, 1}, {0, 2, 0}, {60, 3, 1}, {60, 3, 1}, {0, 2, 0}, {60, 3, 1}, {-120, 3, 3}
    }};
    const Affine reflect{-1.0, 0.0, 0.0, 0.0, 1.0, 0.0};

    std::vector<Affine> transforms;
    transforms.push_back({1.0, 0.0, 0.0, 0.0, 1.0, 0.0});

    int totalAngle = 0;
    Affine rot{1.0, 0.0, 0.0, 0.0, 1.0, 0.0};
    auto tquad = quad;
    for (const auto& rule : tRules) {
        totalAngle += rule[0];
        if (rule[0] != 0) {
            rot = rotation(radians(static_cast<double>(totalAngle)));
            for (int i = 0; i < 4; ++i) tquad[i] = apply(rot, quad[i]);
        }
        const auto slide = transTo(tquad[rule[2]], apply(transforms.back(), quad[rule[1]]));
        transforms.push_back(mul(slide, rot));
    }

    for (auto& t : transforms) t = mul(reflect, t);
    return transforms;
}

void emitBaseTile(const std::string& label, const std::vector<Point>& shape, const Affine& xform, std::vector<Tile>& out) {
    Tile tile;
    tile.label = label;
    tile.points.reserve(shape.size());
    for (const auto& p : shape) tile.points.push_back(apply(xform, p));
    out.push_back(std::move(tile));
}

void flattenLabel(const std::string& label,
                  int depth,
                  const Affine& xform,
                  const std::vector<Point>& shape,
                  const std::vector<Affine>& transforms,
                  const RuleMap& rules,
                  std::vector<Tile>& out) {
    if (depth == 0) {
        if (label == "Gamma") {
            emitBaseTile("Gamma1", shape, xform, out);
            const auto gamma2 = mul(xform, mul(translation(shape[8].x, shape[8].y), rotation(kPi / 6.0)));
            emitBaseTile("Gamma2", shape, gamma2, out);
            return;
        }
        emitBaseTile(label, shape, xform, out);
        return;
    }

    const auto it = rules.find(label);
    if (it == rules.end()) throw std::runtime_error("Unknown label: " + label);
    for (std::size_t i = 0; i < it->second.size(); ++i) {
        if (it->second[i] == "null") continue;
        flattenLabel(it->second[i], depth - 1, mul(xform, transforms[i]), shape, transforms, rules, out);
    }
}

std::string edgeKey(const Point& a, const Point& b) {
    const auto ax = static_cast<long long>(std::llround(a.x * 1000000.0));
    const auto ay = static_cast<long long>(std::llround(a.y * 1000000.0));
    const auto bx = static_cast<long long>(std::llround(b.x * 1000000.0));
    const auto by = static_cast<long long>(std::llround(b.y * 1000000.0));
    std::ostringstream oss;
    if (ax < bx || (ax == bx && ay <= by)) {
        oss << ax << ',' << ay << ',' << bx << ',' << by;
    } else {
        oss << bx << ',' << by << ',' << ax << ',' << ay;
    }
    return oss.str();
}

void buildAdjacency(std::vector<Tile>& tiles) {
    struct EdgeRef { int tile; int face; };
    std::unordered_map<std::string, std::vector<EdgeRef>> edgeMap;

    for (int tileIndex = 0; tileIndex < static_cast<int>(tiles.size()); ++tileIndex) {
        auto& tile = tiles[tileIndex];
        tile.neighbors.assign(tile.points.size(), -1);
        for (int face = 0; face < static_cast<int>(tile.points.size()); ++face) {
            const auto& a = tile.points[face];
            const auto& b = tile.points[(face + 1) % tile.points.size()];
            edgeMap[edgeKey(a, b)].push_back({tileIndex, face});
        }
    }

    for (const auto& entry : edgeMap) {
        const auto& refs = entry.second;
        if (refs.size() < 2) continue;
        for (std::size_t i = 0; i < refs.size(); ++i) {
            for (std::size_t j = i + 1; j < refs.size(); ++j) {
                tiles[refs[i].tile].neighbors[refs[i].face] = refs[j].tile;
                tiles[refs[j].tile].neighbors[refs[j].face] = refs[i].tile;
            }
        }
    }
}

int countLiveFaces(const std::vector<Tile>& tiles, const Tile& tile, int faceWeight) {
    std::unordered_map<int, int> counts;
    for (int neighbor : tile.neighbors) {
        if (neighbor < 0 || !tiles[neighbor].alive) continue;
        ++counts[neighbor];
    }
    int total = 0;
    for (const auto& entry : counts) total += std::min(entry.second, faceWeight);
    return total;
}

void stepAutomata(std::vector<Tile>& tiles, const Config& cfg) {
    for (auto& tile : tiles) {
        const int liveFaces = countLiveFaces(tiles, tile, cfg.faceWeight);
        tile.nextAlive = tile.alive
            ? (liveFaces >= cfg.surviveMin && liveFaces <= cfg.surviveMax)
            : (liveFaces == cfg.birthFaces);
    }
    for (auto& tile : tiles) {
        tile.alive = tile.nextAlive;
        tile.age = tile.alive ? (tile.age + 1) : 0;
    }
}

Bounds computeBounds(const std::vector<Tile>& tiles) {
    Bounds bounds;
    bounds.minX = bounds.minY = std::numeric_limits<double>::infinity();
    bounds.maxX = bounds.maxY = -std::numeric_limits<double>::infinity();
    for (const auto& tile : tiles) {
        for (const auto& p : tile.points) {
            bounds.minX = std::min(bounds.minX, p.x);
            bounds.minY = std::min(bounds.minY, p.y);
            bounds.maxX = std::max(bounds.maxX, p.x);
            bounds.maxY = std::max(bounds.maxY, p.y);
        }
    }
    return bounds;
}

void writeSvg(const std::vector<Tile>& tiles, const std::string& path) {
    const auto bounds = computeBounds(tiles);
    const double padding = 24.0;
    const double width = (bounds.maxX - bounds.minX) + padding * 2.0;
    const double height = (bounds.maxY - bounds.minY) + padding * 2.0;

    std::ofstream out(path);
    if (!out) throw std::runtime_error("Unable to write SVG: " + path);

    out << "<svg xmlns=\"http://www.w3.org/2000/svg\" viewBox=\"0 0 " << width << ' ' << height
        << "\" width=\"" << width << "\" height=\"" << height << "\">\n";
    out << "<rect width=\"100%\" height=\"100%\" fill=\"#090a0f\"/>\n";
    for (const auto& tile : tiles) {
        out << "<polygon points=\"";
        for (const auto& p : tile.points) {
            out << std::fixed << std::setprecision(4)
                << (p.x - bounds.minX + padding) << ',' << (p.y - bounds.minY + padding) << ' ';
        }
        out << "\" fill=\"" << (tile.alive ? "#96ff5f" : "#171b26")
            << "\" stroke=\"" << (tile.alive ? "#d9ffb8" : "#394252")
            << "\" stroke-width=\"0.03\"/>\n";
    }
    out << "</svg>\n";
}

std::vector<Tile> buildTiles(const Config& cfg) {
    const auto shape = baseSpectre();
    const auto quad = baseQuad(shape);
    const auto transforms = childTransforms(quad);
    const auto rules = superRules();

    std::vector<Tile> tiles;
    flattenLabel(cfg.root, cfg.depth, {1.0, 0.0, 0.0, 0.0, 1.0, 0.0}, shape, transforms, rules, tiles);
    buildAdjacency(tiles);
    for (std::size_t i = 0; i < tiles.size(); ++i) {
        tiles[i].alive = randomUnit(cfg.seed + static_cast<std::uint64_t>(i) * 7919ULL) < cfg.density;
        tiles[i].age = tiles[i].alive ? 1 : 0;
    }
    return tiles;
}

void printUsage() {
    std::cout
        << "spectre_face_ca [options]\n"
        << "  --depth N          substitution depth (default 3)\n"
        << "  --root NAME        Delta | Gamma | Psi | Sigma (default Delta)\n"
        << "  --birth N          birth face count (default 3)\n"
        << "  --survive-min N    minimum live faces to survive (default 2)\n"
        << "  --survive-max N    maximum live faces to survive (default 3)\n"
        << "  --face-weight N    cap per-neighbor multi-face contribution (default 1)\n"
        << "  --density X        initial density from 0..1 (default 0.42)\n"
        << "  --steps N          automata steps to run (default 120)\n"
        << "  --seed N           deterministic seed (default 1)\n"
        << "  --svg PATH         write current state to an SVG file\n"
        << "  --help             show this help\n";
}

Config parseArgs(int argc, char** argv) {
    Config cfg;
    for (int i = 1; i < argc; ++i) {
        const std::string arg = argv[i];
        auto next = [&]() -> std::string {
            if (i + 1 >= argc) throw std::runtime_error("Missing value after " + arg);
            return argv[++i];
        };

        if (arg == "--depth") cfg.depth = std::stoi(next());
        else if (arg == "--root") cfg.root = next();
        else if (arg == "--birth") cfg.birthFaces = std::stoi(next());
        else if (arg == "--survive-min") cfg.surviveMin = std::stoi(next());
        else if (arg == "--survive-max") cfg.surviveMax = std::stoi(next());
        else if (arg == "--face-weight") cfg.faceWeight = std::stoi(next());
        else if (arg == "--density") cfg.density = std::stod(next());
        else if (arg == "--steps") cfg.steps = std::stoi(next());
        else if (arg == "--seed") cfg.seed = static_cast<std::uint64_t>(std::stoull(next()));
        else if (arg == "--svg") cfg.svgPath = next();
        else if (arg == "--help") { printUsage(); std::exit(0); }
        else throw std::runtime_error("Unknown argument: " + arg);
    }

    cfg.depth = std::max(0, std::min(5, cfg.depth));
    cfg.birthFaces = std::max(0, std::min(14, cfg.birthFaces));
    cfg.surviveMin = std::max(0, std::min(14, cfg.surviveMin));
    cfg.surviveMax = std::max(0, std::min(14, cfg.surviveMax));
    if (cfg.surviveMin > cfg.surviveMax) std::swap(cfg.surviveMin, cfg.surviveMax);
    cfg.faceWeight = std::max(1, std::min(6, cfg.faceWeight));
    cfg.density = std::max(0.0, std::min(1.0, cfg.density));
    cfg.steps = std::max(0, cfg.steps);
    if (cfg.root != "Delta" && cfg.root != "Gamma" && cfg.root != "Psi" && cfg.root != "Sigma" &&
        cfg.root != "Theta" && cfg.root != "Lambda" && cfg.root != "Xi" && cfg.root != "Pi" && cfg.root != "Phi") {
        throw std::runtime_error("Unsupported root tile: " + cfg.root);
    }
    return cfg;
}

}  // namespace

int main(int argc, char** argv) {
    try {
        const Config cfg = parseArgs(argc, argv);
        auto tiles = buildTiles(cfg);

        std::cout << "Spectre Face CA\n";
        std::cout << "tiles=" << tiles.size()
                  << " depth=" << cfg.depth
                  << " root=" << cfg.root
                  << " rule=B" << cfg.birthFaces << "/S" << cfg.surviveMin << '-' << cfg.surviveMax
                  << " weight=" << cfg.faceWeight
                  << " density=" << std::fixed << std::setprecision(2) << cfg.density
                  << " seed=" << cfg.seed << '\n';

        for (int step = 0; step < cfg.steps; ++step) {
            stepAutomata(tiles, cfg);
            int alive = 0;
            int maxAge = 0;
            for (const auto& tile : tiles) {
                if (tile.alive) ++alive;
                maxAge = std::max(maxAge, tile.age);
            }
            std::cout << "step=" << (step + 1) << " alive=" << alive << " max_age=" << maxAge << '\n';
        }

        if (!cfg.svgPath.empty()) {
            writeSvg(tiles, cfg.svgPath);
            std::cout << "svg=" << cfg.svgPath << '\n';
        }
        return 0;
    } catch (const std::exception& ex) {
        std::cerr << "error: " << ex.what() << '\n';
        return 1;
    }
}