#include <algorithm>
#include <fstream>
#include <random>
#include <set>

// https://github.com/salvipeter/libgeom/
#include <geometry.hh>

using namespace Geometry;

#define USE_TOP

std::vector<std::string> front = {
  "#### #### ####",
  "#    #  # #   ",
  "# ## #### ####",
  "#  # #       #",
  "#### #    ####"
}, side = {
  "### ###",
  "# # # #",
  "  # # #",
  "  # # #",
  "  # ###"
}, top = {
  "### ###",
  "# # # #",
  "### ###",
  "#   # #",
  "       ",
  "#   ###",
  "#   # #",
  "#   # #",
  "### ###",
  "       ",
  "#   ###",
  "#   # #",
  "#   # #",
  "### ###"    
};

PointVector adjacent(const Point3D &p) {
  return {
    { p[0] - 1, p[1], p[2] },
    { p[0] + 1, p[1], p[2] },
    { p[0], p[1] - 1, p[2] },
    { p[0], p[1] + 1, p[2] },
    { p[0], p[1], p[2] - 1 },
    { p[0], p[1], p[2] + 1 }
  };
}

// is `points` connected without `p`?
bool connected(const std::list<Point3D> &points, const Point3D &p) {
  auto cmp = [](const Point3D &a, const Point3D &b) {
    for (size_t i = 0; i < 3; ++i) {
      if (a[i] < b[i])
        return true;
      if (a[i] > b[i])
        return false;
    }
    return false;
  };
  std::set<Point3D, decltype(cmp)> all(points.begin(), points.end(), cmp), remaining(cmp);
  std::list<Point3D> queue;
  if (all.contains(p))
    all.erase(p);
  for (const auto &q : all) {
    if (q[1] == 0)
      queue.push_back(q);
    else
      remaining.insert(q);
  }
  while (!queue.empty() && !remaining.empty()) {
    auto next = queue.front();
    queue.pop_front();
    for (const auto &q : adjacent(next))
      if (remaining.contains(q)) {
        queue.push_back(q);
        remaining.erase(q);
      }
  }
  return remaining.empty();
}

bool removable(const std::list<Point3D> &points, const Point3D &p) {
  // not removable if:
  // - it is the only point on a ray
  // - it is the only connection of an adjacent block to the floor
  // 1. check xy
  bool found = false;
  for (const auto &q : points)
    if (p[0] == q[0] && p[1] == q[1] && p[2] != q[2]) {
      found = true;
      break;
    }
  if (!found)
    return false;
  // 2. check yz
  found = false;
  for (const auto &q : points)
    if (p[0] != q[0] && p[1] == q[1] && p[2] == q[2]) {
      found = true;
      break;
    }
  if (!found)
    return false;
  // 3. check xz
#ifdef USE_TOP
  found = false;
  for (const auto &q : points)
    if (p[0] == q[0] && p[1] != q[1] && p[2] == q[2]) {
      found = true;
      break;
    }
  if (!found)
    return false;
#endif  
  // 4. check connectivity
  return connected(points, p);
};

// What to minimize:
// - # of cubes
// - # of cubes not supported from below
size_t energy(const std::list<Point3D> &points) {
  // return points.size();
  size_t count = 0;
  for (const auto &p : points)
    if (p[1] != 0) {
      bool found = false;
      for (const auto &q : points)
        if (q[0] == p[0] && q[1] == p[1] - 1 && q[2] == p[2]) {
          found = true;
          break;
        }
      if (!found)
        count++;
    }
  return count;
}

std::list<Point3D> findMinimal(const std::list<Point3D> &points, size_t iterations) {
  std::random_device rd;
  std::mt19937 g(rd());

  size_t min = energy(points);
  std::list<Point3D> result = points;

  for (size_t i = 0; i < iterations; ++i) {
    PointVector deck(points.begin(), points.end());
    std::shuffle(deck.begin(), deck.end(), g);
    std::list<Point3D> remaining(deck.begin(), deck.end());

    bool changed;
    do {
      changed = false;
      auto p = remaining.begin();
      while (p != remaining.end())
        if (removable(remaining, *p)) {
          p = remaining.erase(p);
          changed = true;
        } else
          p++;
    } while(changed);

    auto e = energy(remaining);
    if (e < min) {
      min = e;
      result = remaining;
    }
  }

  return result;
}

int main(int argc, char **argv) {
  // Initialize the whole block
  size_t nx = front[0].size(), ny = front.size(), nz = side[0].size();
  std::list<Point3D> points;
  for (size_t i = 0; i < nx; ++i)
    for (size_t j = 0; j < ny; ++j)
      for (size_t k = 0; k < nz; ++k)
        points.emplace_back(i, j, k);

  // Create the maximal solution
  for (size_t i = 0; i < nx; ++i)
    for (size_t j = 0; j < ny; ++j)
      if (front[ny-j-1][i] == ' ') {
        auto p = points.begin();
        while (p != points.end())
          if ((*p)[0] == i && (*p)[1] == j)
            p = points.erase(p);
          else
            p++;
      }
  for (size_t j = 0; j < ny; ++j)
    for (size_t k = 0; k < nz; ++k)
      if (side[ny-j-1][nz-k-1] == ' ') {
        auto p = points.begin();
        while (p != points.end())
          if ((*p)[1] == j && (*p)[2] == k)
            p = points.erase(p);
          else
            p++;
      }
#ifdef USE_TOP
  for (size_t i = 0; i < nx; ++i)
    for (size_t k = 0; k < nz; ++k)
      if (top[nx-i-1][k] == ' ') {
        auto p = points.begin();
        while (p != points.end())
          if ((*p)[0] == i && (*p)[2] == k)
            p = points.erase(p);
          else
            p++;
      }
#endif

  if (!connected(points, Point3D(-1, -1, -1))) {
    std::cerr << "No connected solution!" << std::endl;
    return 1;
  }

  points = findMinimal(points, 100);

  std::cout << "Solution size: " << points.size() << std::endl;
  std::cout << "Solution energy: " << energy(points) << std::endl;
  std::ofstream f("gps70.obj");
  for (const auto &p : points)
    f << "v " << p << std::endl;
}