path: root/lib/src/graph.cpp

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

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

// actual mod for floats
double mod(double a, double b) {
  if (a >= 0) {
    return fmod(a, b);
  } else {
    return fmod(a + b * ceil(-a / b), b);
  }
}

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

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

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

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

  for (unsigned 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));
    signed f = (1 + i / (Nx / 2)) % 2 == 1 ? 1 : -1;
    vertices[i].nd = {i, (i + Nx / 2) % nv, Nx / 2 * (i / (Nx / 2)) + ((i + f) % (Nx / 2)), (nv + i - Nx / 2) % nv};
    vertices[i].polygon = {
      {vertices[i].r.x - 1.0, vertices[i].r.y},
      {vertices[i].r.x, vertices[i].r.y - 1.0},
      {vertices[i].r.x + 1.0, vertices[i].r.y},
      {vertices[i].r.x, vertices[i].r.y + 1.0}
    };

    dual_vertices[i].r.x = (double)((i / (Nx / 2)) % 2 + 2 * (i % (Nx / 2)));
    dual_vertices[i].r.y = (double)(i / (Nx / 2));
    dual_vertices[i].nd = {i, (i + Nx / 2) % nv, Nx * (i / (Nx / 2)) + (i % (Nx / 2)), (nv + i - Nx / 2) % nv};
    dual_vertices[i].polygon = {
      {dual_vertices[i].r.x - 1.0, vertices[i].r.y},
      {dual_vertices[i].r.x, vertices[i].r.y - 1.0},
      {dual_vertices[i].r.x + 1.0, vertices[i].r.y},
      {dual_vertices[i].r.x, vertices[i].r.y + 1.0}
    };
  }

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

      bool crossed_x = x == Nx - 1;
      bool crossed_y = y == Ny - 1;

      edges.push_back({{v1, v2}, {0.5 + (double)x, 0.5 + (double)y}, {crossed_x, crossed_y}});

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

      dual_edges.push_back({{dv1, dv2}, {0.5 + (double)x, 0.5 + (double)y}, {crossed_x, crossed_y}});
    }
  }

  for (vertex& v : vertices) {
    v.ne.resize(v.nd.size());
  }

  for (unsigned i = 0; i < edges.size(); i++) {
    for (unsigned vi : edges[i].v) {
      auto it1 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[0]);
      auto it2 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[1]);

      unsigned d1 = std::distance(vertices[vi].nd.begin(), it1);
      unsigned d2 = std::distance(vertices[vi].nd.begin(), it2);

      unsigned dv1 = d1 < d2 ? d1 : d2;
      unsigned dv2 = d1 < d2 ? d2 : d1;

      if (dv2 - dv1 == 1) {
        vertices[vi].ne[dv1] = i;
      } else {
        vertices[vi].ne[dv2] = i;
      }
    }
  }

}

class eulerException: public std::exception
{
  virtual const char* what() const throw()
  {
    return "The voronoi tessellation generated has the incorrect Euler characteristic for a torus and is malformed.";
  }
} eulerex;

class clippingException: public std::exception
{
  virtual const char* what() const throw()
  {
    return "An interior site has a clipped edge and the tesselation is malformed.";
  }
} clippingex;

class triangleException: public std::exception
{
  virtual const char* what() const throw()
  {
    return "A dual-centered triangle has an impossible geometry and the tesselation is malformed.";
  }
} triex;

unsigned get_triangle_signature(unsigned j1, unsigned j2, unsigned j3) {
  // this yucky function takes three unsignedegers representing the
  // location in the nine periodic copies of each corner of a delauney triangle
  // and returns a signature for that triangle, which is an unsignedeger
  // that uniquely labels the way the triangle crosses boundaries of the
  // copies. This allows us to differentiate delauney triangles with identical
  // vertices but which should be identified with different faces
  unsigned x1 = j1 % 3;
  unsigned y1 = j1 / 3;

  unsigned x2 = j2 % 3;
  unsigned y2 = j2 / 3;

  unsigned x3 = j3 % 3;
  unsigned y3 = j3 / 3;

  if ((j1 == j2) && (j2 == j3)) {
    return 0;
  } else if (((j1 == j2) && (x2 < x3) && (y2 == y3)) || ((j1 == j3) && (x3 < x2) && (y2 == y3)) || ((j2 == j3) && (x3 < x1) && (y1 == y3))) {
    return 1;
  } else if (((j1 == j2) && (x2 > x3) && (y2 == y3)) || ((j1 == j3) && (x3 > x2) && (y2 == y3)) || ((j2 == j3) && (x3 > x1) && (y1 == y3))) {
    return 2;
  } else if (((j1 == j2) && (y2 < y3) && (x2 == x3)) || ((j1 == j3) && (y3 < y2) && (x2 == x3)) || ((j2 == j3) && (y3 < y1) && (x1 == x3))) {
    return 3;
  } else if (((j1 == j2) && (y2 > y3) && (x2 == x3)) || ((j1 == j3) && (y3 > y2) && (x2 == x3)) || ((j2 == j3) && (y3 > y1) && (x1 == x3))) {
    return 4;
  } else if (((j1 == j2) && (x2 < x3) && (y2 < y3)) || ((j1 == j3) && (x3 < x2) && (y3 < y2)) || ((j2 == j3) && (x3 < x1) && (y3 < y1))) {
    return 5;
  } else if (((j1 == j2) && (x2 < x3) && (y2 > y3)) || ((j1 == j3) && (x3 < x2) && (y3 > y2)) || ((j2 == j3) && (x3 < x1) && (y3 > y1))) {
    return 6;
  } else if (((j1 == j2) && (x2 > x3) && (y2 < y3)) || ((j1 == j3) && (x3 > x2) && (y3 < y2)) || ((j2 == j3) && (x3 > x1) && (y3 < y1))) {
    return 7;
  } else if (((j1 == j2) && (x2 > x3) && (y2 > y3)) || ((j1 == j3) && (x3 > x2) && (y3 > y2)) || ((j2 == j3) && (x3 > x1) && (y3 > y1))) {
    return 8;
  } else if (((x1 == x2) && (x2 < x3) && ((y1 < y3) || (y2 < y3))) || ((x1 == x3) && (x3 < x2) && ((y1 < y2) || (y3 < y2))) || ((x2 == x3) && (x2 < x1) && ((y2 < y1) || (y3 < y1)))) {
    return 9;
  } else if (((x1 == x2) && (x2 > x3) && ((y1 < y3) || (y2 < y3))) || ((x1 == x3) && (x3 > x2) && ((y1 < y2) || (y3 < y2))) || ((x2 == x3) && (x2 > x1) && ((y2 < y1) || (y3 < y1)))) {
    return 10;
  } else if (((x1 == x2) && (x2 < x3) && ((y1 > y3) || (y2 > y3))) || ((x1 == x3) && (x3 < x2) && ((y1 > y2) || (y3 > y2))) || ((x2 == x3) && (x2 < x1) && ((y2 > y1) || (y3 > y1)))) {
    return 11;
  } else if (((x1 == x2) && (x2 > x3) && ((y1 > y3) || (y2 > y3))) || ((x1 == x3) && (x3 > x2) && ((y1 > y2) || (y3 > y2))) || ((x2 == x3) && (x2 > x1) && ((y2 > y1) || (y3 > y1)))) {
    return 12;
  } else {
    throw triex;
  }
}

graph::graph(double Lx, double Ly, std::mt19937& rng) {
  // randomly choose N to be floor(Lx * Ly / 2) or ceil(Lx * Ly / 2) with
  // probability proportional to the distance from each
  std::uniform_real_distribution<double> d(0.0, 1.0);
  unsigned N = round(Lx * Ly / 2 + d(rng) - 0.5);

  L = {Lx, Ly};

  this->helper(N, rng);
}

graph::graph(unsigned n, double a, std::mt19937& rng) {
  L = {sqrt(2 * n * a), sqrt(2 * n / a)};

  this->helper(n, rng);
}

void graph::helper(unsigned nv, std::mt19937& rng) {
  std::uniform_real_distribution<double> d(0.0, 1.0);
  vertices.resize(nv);

  // the coordinates of the lattice, from which a delaunay triangulation
  // and corresponding voronoi tessellation will be built. Everyone is in the
  // rectangle (0, 0) (Lx, Ly)
  for (vertex &v : vertices) {
    v = {{L.x * d(rng), L.y * d(rng)}};
  }

  // set up the voronoi objects
  jcv_diagram diagram;
  memset(&diagram, 0, sizeof(jcv_diagram));
  jcv_rect bounds = {{-L.x, -L.y}, {2 * L.x, 2 * L.y}};
  std::vector<jcv_point> points(9 * nv);

  for (unsigned i = 0; i < nv; i++) {
    const vertex& v = vertices[i];
    points[9 * i + 0] = {v.r.x - L.x, v.r.y - L.y};
    points[9 * i + 1] = {v.r.x + 0.0, v.r.y - L.y};
    points[9 * i + 2] = {v.r.x + L.x, v.r.y - L.y};
    points[9 * i + 3] = {v.r.x - L.x, v.r.y + 0.0};
    points[9 * i + 4] = {v.r.x + 0.0, v.r.y + 0.0};
    points[9 * i + 5] = {v.r.x + L.x, v.r.y + 0.0};
    points[9 * i + 6] = {v.r.x - L.x, v.r.y + L.y};
    points[9 * i + 7] = {v.r.x + 0.0, v.r.y + L.y};
    points[9 * i + 8] = {v.r.x + L.x, v.r.y + L.y};
  }

  jcv_diagram_generate(9 * nv, points.data(), &bounds, &diagram);

  const jcv_site* sites = jcv_diagram_get_sites(&diagram);

  std::unordered_map<std::array<unsigned, 4>, unsigned> 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 == 4) {
      bool self_bonded = false;
      unsigned i1 = (unsigned)(site->index / 9);
      unsigned j1 = (unsigned)(site->index % 9);
      const jcv_graphedge* e = site->edges;
      const jcv_graphedge* ep = site->edges;
      while (ep->next) {
        ep = ep->next;
      }
      // for each edge on the site's cell
      while(e) {
        // assess whether the edge needs to be added. only neighboring
        // sites whose index is larger than ours are given edges, so we
        // don't double up later
        const jcv_site* neighbor = e->neighbor;
        if (neighbor == NULL) {
          throw clippingex;
        }
        unsigned i2 = (unsigned)(neighbor->index / 9);
        unsigned j2 = (unsigned)(neighbor->index % 9);
        unsigned x2 = j2 % 3;
        unsigned y2 = j2 / 3;

        vertices[i1].polygon.push_back({e->pos[0].x, e->pos[0].y});

        if (ep->neighbor == NULL) {
          throw clippingex;
        }
        unsigned i3p = (unsigned)(ep->neighbor->index / 9);
        unsigned j3p = (unsigned)(ep->neighbor->index % 9);

        unsigned sig1 = get_triangle_signature(j1, j2, j3p);

        std::array<unsigned, 4> t1 = {i1, i2, i3p, sig1};
        std::sort(t1.begin(), t1.begin() + 3);

        auto it1 = known_vertices.find(t1);

        unsigned vi1;

        if (it1 == known_vertices.end()) {
          vi1 = dual_vertices.size();
          dual_vertices.push_back({{mod(e->pos[0].x, L.x), mod(e->pos[0].y, L.y)},{i1},{},{{site->p.x, site->p.y}}});
          known_vertices[t1] = vi1;
        } else {
          vi1 = it1->second;
          dual_vertices[vi1].nd.push_back(i1);
          dual_vertices[vi1].polygon.push_back({site->p.x, site->p.y});
        }

        vertices[i1].nd.push_back(vi1);
        
        if (i1 < i2 || (i1 == i2 && !self_bonded)) {
          if (i1 == i2) {
            self_bonded = true;
          }
          bool crossed_x = x2 != 1;
          bool crossed_y = y2 != 1;
          edges.push_back({{i1, i2},
              {mod((site->p.x + neighbor->p.x) / 2, L.x),
               mod((site->p.y + neighbor->p.y) / 2, L.y)},
               {crossed_x, crossed_y}
               });

          jcv_graphedge *en;
          if (e->next == NULL) {
            en = site->edges;
          } else {
            en = e->next;
          }

          if (en->neighbor == NULL) {
            throw clippingex;
          }

          unsigned i3n = (unsigned)(en->neighbor->index / 9);
          unsigned j3n = (unsigned)(en->neighbor->index % 9);

          unsigned sig2 = get_triangle_signature(j1, j2, j3n);

          std::array<unsigned, 4> t2 = {i1, i2, i3n, sig2};
          std::sort(t2.begin(), t2.begin() + 3);

          auto it2 = known_vertices.find(t2);

          unsigned vi2;

          if (it2 == known_vertices.end()) {
            vi2 = dual_vertices.size();
            dual_vertices.push_back({{mod(e->pos[1].x, L.x), mod(e->pos[1].y, L.y)},{}});
            known_vertices[t2] = vi2;
          } else {
            vi2 = it2->second;
          }

          bool dcrossed_x = (unsigned)floor(e->pos[0].x / L.x) != (unsigned)floor(e->pos[1].x / L.x);
          bool dcrossed_y = (unsigned)floor(e->pos[0].y / L.y) != (unsigned)floor(e->pos[1].y / L.y);

          dual_edges.push_back({{vi1, vi2},
              {mod((e->pos[0].x + e->pos[1].x) / 2, L.x),
               mod((e->pos[0].y + e->pos[1].y) / 2, L.y)},
               {dcrossed_x, dcrossed_y}
               });
        }

        ep = e;
        e = e->next;
      }
    }
  }

  if (edges.size() != vertices.size() + dual_vertices.size()) {
    throw eulerex;
  }

  for (vertex &v : dual_vertices) {
    if (fabs(v.polygon[0].x - v.polygon[1].x) > L.x / 2) {
      if (v.polygon[0].x < L.x / 2) {
        v.polygon[0].x += L.x;
      } else {
        v.polygon[1].x += L.x;
      }
    }

    if (fabs(v.polygon[0].y - v.polygon[1].y) > L.y / 2) {
      if (v.polygon[0].y < L.y / 2) {
        v.polygon[0].y += L.y;
      } else {
        v.polygon[1].y += L.y;
      }
    }

    if (fabs(v.polygon[2].x - v.polygon[1].x) > L.x / 2) {
      if (v.polygon[2].x < L.x / 2) {
        v.polygon[2].x += L.x;
      } else {
        v.polygon[1].x += L.x;
      }
    }

    if (fabs(v.polygon[2].y - v.polygon[1].y) > L.y / 2) {
      if (v.polygon[2].y < L.y / 2) {
        v.polygon[2].y += L.y;
      } else {
        v.polygon[1].y += L.y;
      }
    }

    if (fabs(v.polygon[2].x - v.polygon[0].x) > L.x / 2) {
      if (v.polygon[2].x < L.x / 2) {
        v.polygon[2].x += L.x;
      } else {
        v.polygon[0].x += L.x;
      }
    }

    if (fabs(v.polygon[2].y - v.polygon[0].y) > L.y / 2) {
      if (v.polygon[2].y < L.y / 2) {
        v.polygon[2].y += L.y;
      } else {
        v.polygon[0].y += L.y;
      }
    }
  }

  for (vertex& v : vertices) {
    v.ne.resize(v.nd.size());
  }

  for (unsigned i = 0; i < edges.size(); i++) {
    for (unsigned vi : edges[i].v) {
      auto it1 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[0]);
      auto it2 = std::find(vertices[vi].nd.begin(), vertices[vi].nd.end(), dual_edges[i].v[1]);

      unsigned d1 = std::distance(vertices[vi].nd.begin(), it1);
      unsigned d2 = std::distance(vertices[vi].nd.begin(), it2);

      unsigned dv1 = d1 < d2 ? d1 : d2;
      unsigned dv2 = d1 < d2 ? d2 : d1;

      if (dv2 - dv1 == 1) {
        vertices[vi].ne[dv1] = i;
      } else {
        vertices[vi].ne[dv2] = i;
      }
    }
  }

  jcv_diagram_free(&diagram);
}

graph::coordinate reverse(const graph::coordinate &x) {
  return {x.y, x.x};
}

graph const graph::rotate() {
  graph g2(*this);

  for (graph::vertex &v : g2.vertices) {
    v.r = reverse(v.r);
    for (graph::coordinate &r : v.polygon) {
      r = reverse(r);
    }
  }

  for (graph::edge &e : g2.edges) {
    e.r = reverse(e.r);
    e.crossings = {e.crossings[1], e.crossings[0]};
  }

  for (graph::vertex &v : g2.dual_vertices) {
    v.r = reverse(v.r);
    for (graph::coordinate &r : v.polygon) {
      r = reverse(r);
    }
  }

  for (graph::edge &e : g2.dual_edges) {
    e.r = reverse(e.r);
    e.crossings = {e.crossings[1], e.crossings[0]};
  }

  g2.L = reverse(g2.L);

  return g2;
}