#include "graph.hpp"
#include <cstring>
#include <iostream>

#define JC_VORONOI_IMPLEMENTATION
#define JCV_REAL_TYPE double
#define JCV_ATAN2 atan2
#define JCV_FLT_MAX 1.7976931348623157E+308
#include <jc_voronoi.h>

graph::graph(unsigned int Nx, unsigned int Ny) {
  L = {(double)Nx, (double)Ny};

  unsigned int ne = Nx * Ny;
  unsigned int nv = ne / 2;

  vertices.resize(nv);
  edges.reserve(ne);

  dual_vertices.resize(nv);
  dual_edges.reserve(ne);

  for (unsigned int i = 0; i < nv; i++) {
    vertices[i].r.x = (double)((1 + i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
    vertices[i].r.y = (double)(i / (Nx / 2));

    dual_vertices[i].r.x = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
    dual_vertices[i].r.y = (double)(i / (Nx / 2));
  }

  for (unsigned int x = 0; x < Ny; x++) {
    for (unsigned int y = 0; y < Nx; y++) {
      unsigned int v1 = (Nx * x) / 2 + ((y + x % 2) / 2) % (Nx / 2);
      unsigned int v2 = ((Nx * (x + 1)) / 2 + ((y + (x + 1) % 2) / 2) % (Nx / 2)) % nv;

      edges.push_back({v1, v2});

      unsigned int dv1 = (Nx * x) / 2 + ((y + (x + 1) % 2) / 2) % (Nx / 2);
      unsigned int dv2 = ((Nx * (x + 1)) / 2 + ((y + x % 2) / 2) % (Nx / 2)) % nv;

      dual_edges.push_back({dv1, dv2});
    }
  }
}

graph::graph(unsigned int Nx, unsigned int Ny, std::mt19937& rng, double spread) {
  L = {(double)Nx, (double)Ny};

  unsigned int dnv = Nx * Ny / 2;

  dual_vertices.resize(dnv);
  
  std::normal_distribution<double> d(0.0, spread);

  // the coordinates of the dual lattice, from which a delaunay triangulation
  // and corresponding voronoi tesselation will be built
  for (unsigned int i = 0; i < dnv; i++) {
    double rx = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2))) + d(rng);
    double ry = (double)(i / (Nx / 2)) + d(rng);
    dual_vertices[i] = {{fmod(L.x + rx, L.x), fmod(L.y + ry, L.y)}};
  }

  // to make the resulting tesselation periodic, we tile eight (!) copies of
  // the original dual points for a total of nine. note that the index of each
  // point quotient 9 is equal to the original index (we'll use this to
  // translate back)
  std::vector<jcv_point> points;
  points.reserve(9 * dnv);
  for (const vertex &v : dual_vertices) {
    points.push_back({v.r.x, v.r.y});
    points.push_back({v.r.x + L.x, v.r.y + 0.0});
    points.push_back({v.r.x - L.x, v.r.y + 0.0});
    points.push_back({v.r.x + 0.0, v.r.y + L.y});
    points.push_back({v.r.x + 0.0, v.r.y - L.y});
    points.push_back({v.r.x + L.x, v.r.y + L.y});
    points.push_back({v.r.x - L.x, v.r.y + L.y});
    points.push_back({v.r.x + L.x, v.r.y - L.y});
    points.push_back({v.r.x - L.x, v.r.y - L.y});
  }

  jcv_diagram diagram;
  memset(&diagram, 0, sizeof(jcv_diagram));

  jcv_diagram_generate(9 * dnv, points.data(), NULL, &diagram);

  const jcv_site* sites = jcv_diagram_get_sites(&diagram);

  struct coorcomp {
    bool operator() (const coordinate& lhs, const coordinate& rhs) const
    {return lhs.x<rhs.x;}
  };

  std::map<coordinate, unsigned int, coorcomp> known_vertices;

  for (int i = 0; i < diagram.numsites; i++) {
    const jcv_site* site = &sites[i];

    // we only care about processing the cells of our original, central sites
    if (site->index % 9 == 0) {
      const jcv_graphedge* e = site->edges;
      // for each edge on the site's cell
      while(e) {
        // assess whether the dual edge needs to be added. only neighboring
        // dual sites whose index is larger than ours are given edges, so we
        // don't double up later
        const jcv_site* neighbor = e->neighbor;
        unsigned int index = (unsigned int)(site->index / 9);
        unsigned int real_index = (unsigned int)(neighbor->index / 9);
        
        if (index < real_index) {
          dual_edges.push_back({index, real_index});
          coordinate r1 = {1e-10 * round(1e10 * fmod(L.x + e->pos[0].x, L.x)), 1e-10 * round(1e10 * fmod(L.y + e->pos[0].y, L.y))};
          coordinate r2 = {1e-10 * round(1e10 * fmod(L.x + e->pos[1].x, L.x)), 1e-10 * round(1e10 * fmod(L.y + e->pos[1].y, L.y))};

          std::map<coordinate, unsigned int>::iterator it1 = known_vertices.find(r1);
          std::map<coordinate, unsigned int>::iterator it2 = known_vertices.find(r2);

          unsigned int vi1, vi2;

          if (it1 == known_vertices.end()) {
            vi1 = vertices.size();
            vertices.push_back({r1});
            known_vertices[r1] = vi1;
          } else {
            vi1 = it1->second;
          }
          if (it2 == known_vertices.end()) {
            vi2 = vertices.size();
            vertices.push_back({r2});
            known_vertices[r2] = vi2;
          } else {
            vi2 = it2->second;
          }

          edges.push_back({vi1, vi2});
        }

        e = e->next;
      }
    }
  }

  for (edge &e : edges) { 
    std::cout << e[0] << " " << e[1] << " ";
  }
  std::cout << "\n";
  for (vertex &v : vertices) { 
    std::cout << v.r.x << " " << v.r.y << " ";
  }
  std::cout << "\n";
  for (edge &e : dual_edges) { 
    std::cout << e[0] << " " << e[1] << " ";
  }
  std::cout << "\n";
  for (vertex &v : dual_vertices) { 
    std::cout << v.r.x << " " << v.r.y << " ";
  }
  std::cout << "\n";

  jcv_diagram_free(&diagram);
}