big library overhual, fixed all examples, added documentation with sphinx

11 files changed, 139 insertions, 123 deletions

diff --git a/lib/include/wolff.hpp b/lib/include/wolff.hpp
index b78cd4b..7e909c8 100644
--- a/lib/include/wolff.hpp
+++ b/lib/include/wolff.hpp
@@ -3,25 +3,30 @@
 #define WOLFF_H
 
 #include "wolff/cluster.hpp"
+#include "wolff/measurement.hpp"
+
+namespace wolff{
 
 template <class R_t, class X_t>
-void wolff(N_t N, wolff_system<R_t, X_t>& S,
-           std::function <R_t(std::mt19937&, X_t)> r_gen,
-           wolff_measurement<R_t, X_t>& A, std::mt19937& rng) {
+void system<R_t, X_t>::run_wolff(N_t N,
+           std::function <R_t(std::mt19937&, const system<R_t, X_t>&, v_t)> r_gen,
+           measurement<R_t, X_t>& A, std::mt19937& rng) {
 
-  std::uniform_int_distribution<v_t> dist(0, S.nv - 1);
+  std::uniform_int_distribution<v_t> dist(0, nv - 1);
 
   for (N_t n = 0; n < N; n++) {
     v_t i0 = dist(rng);
-    R_t r = r_gen(rng, S.s[i0]);
+    R_t r = r_gen(rng, *this, i0);
 
-    A.pre_cluster(n, N, S, i0, r);
+    A.pre_cluster(n, N, *this, i0, r);
 
-    wolff_cluster_flip<R_t, X_t>(S, i0, r, rng, A);
+    this->flip_cluster(i0, r, rng, A);
 
-    A.post_cluster(n, N, S);
+    A.post_cluster(n, N, *this);
   }
 }
 
+}
+
 #endif
 
diff --git a/lib/include/wolff/cluster.hpp b/lib/include/wolff/cluster.hpp
index 055cdf3..b66f367 100644
--- a/lib/include/wolff/cluster.hpp
+++ b/lib/include/wolff/cluster.hpp
@@ -7,20 +7,19 @@
 #include <vector>
 #include <queue>
 
-#include "types.h"
 #include "system.hpp"
-#include "graph.hpp"
-#include "measurement.hpp"
+
+namespace wolff {
 
 template <class R_t, class X_t>
-void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r,
-                  std::mt19937& rng, wolff_measurement<R_t, X_t>& A) {
+void system<R_t, X_t>::flip_cluster(v_t i0, const R_t& r,
+                  std::mt19937& rng, measurement<R_t, X_t>& A) {
   std::uniform_real_distribution<double> dist(0.0, 1.0);
 
   std::queue<v_t> queue;
   queue.push(i0);
 
-  std::vector<bool> marks(S.G.nv, false);
+  std::vector<bool> marks(G.nv, false);
 
   while (!queue.empty()) {
     v_t i = queue.front();
@@ -33,21 +32,21 @@ void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r,
 #ifndef WOLFF_NO_FIELD
       R_t s0_new;
 
-      bool we_are_ghost = (i == S.nv);
+      bool we_are_ghost = (i == nv);
 
       if (we_are_ghost) {
-        s0_new = r.act(S.s0);
+        s0_new = r.act(s0);
       } else
 #endif
       {
-        si_new = r.act(S.s[i]);
+        si_new = r.act(s[i]);
       }
 
-      for (const v_t &j : S.G.adj[i]) {
+      for (const v_t &j : G.adjacency[i]) {
         double dE, p;
 
 #ifndef WOLFF_NO_FIELD
-        bool neighbor_is_ghost = (j == S.nv);
+        bool neighbor_is_ghost = (j == nv);
 
         if (we_are_ghost || neighbor_is_ghost) {
           X_t s0s_old, s0s_new;
@@ -55,18 +54,18 @@ void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r,
 
           if (neighbor_is_ghost) {
             non_ghost = i;
-            s0s_old = S.s0.act_inverse(S.s[i]);
-            s0s_new = S.s0.act_inverse(si_new);
+            s0s_old = s0.act_inverse(s[i]);
+            s0s_new = s0.act_inverse(si_new);
           } else {
             non_ghost = j;
-            s0s_old = S.s0.act_inverse(S.s[j]);
-            s0s_new = s0_new.act_inverse(S.s[j]);
+            s0s_old = s0.act_inverse(s[j]);
+            s0s_new = s0_new.act_inverse(s[j]);
           }
 
 #ifdef WOLFF_SITE_DEPENDENCE
-          dE = S.B(non_ghost, s0s_old) - S.B(non_ghost, s0s_new);
+          dE = B(non_ghost, s0s_old) - B(non_ghost, s0s_new);
 #else
-          dE = S.B(s0s_old) - S.B(s0s_new);
+          dE = B(s0s_old) - B(s0s_new);
 #endif
 
 #ifdef WOLFF_FINITE_STATES
@@ -75,28 +74,28 @@ void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r,
 #endif
 
           // run measurement hooks for encountering a ghost bond
-          A.ghost_bond_visited(S, non_ghost, s0s_old, s0s_new, dE);
+          A.ghost_bond_visited(*this, non_ghost, s0s_old, s0s_new, dE);
         } else // this is a perfectly normal bond!
 #endif
         {
 #ifdef WOLFF_BOND_DEPENDENCE
-          dE = S.Z(i, S.s[i], j, S.s[j]) - S.Z(i, si_new, j, S.s[j]);
+          dE = Z(i, s[i], j, s[j]) - Z(i, si_new, j, s[j]);
 #else
-          dE = S.Z(S.s[i], S.s[j]) - S.Z(si_new, S.s[j]);
+          dE = Z(s[i], s[j]) - Z(si_new, s[j]);
 #endif
 
 #ifdef WOLFF_FINITE_STATES
-          p = finite_states_Zp[finite_states_enum(S.s[i])]
+          p = finite_states_Zp[finite_states_enum(s[i])]
                               [finite_states_enum(si_new)]
-                              [finite_states_enum(S.s[j])];
+                              [finite_states_enum(s[j])];
 #endif
 
           // run measurement hooks for encountering a plain bond
-          A.plain_bond_visited(S, i, si_new, j, dE);
+          A.plain_bond_visited(*this, i, si_new, j, dE);
         }
 
 #ifndef FINITE_STATES
-        p = 1.0 - exp(-dE / S.T);
+        p = 1.0 - exp(-dE / T);
 #endif
 
         if (dist(rng) < p) {
@@ -106,17 +105,19 @@ void wolff_cluster_flip(wolff_system<R_t, X_t>& S, v_t i0, const R_t& r,
 
 #ifndef WOLFF_NO_FIELD
       if (we_are_ghost) {
-        A.ghost_site_transformed(S, s0_new);
-        S.s0 = s0_new;
+        A.ghost_site_transformed(*this, s0_new);
+        s0 = s0_new;
       } else
 #endif
       {
-        A.plain_site_transformed(S, i, si_new);
-        S.s[i] = si_new;
+        A.plain_site_transformed(*this, i, si_new);
+        s[i] = si_new;
       }
     }
   }
 }
 
+}
+
 #endif
 
diff --git a/lib/include/wolff/finite_states.hpp b/lib/include/wolff/finite_states.hpp
index 4cca69e..7e54eff 100644
--- a/lib/include/wolff/finite_states.hpp
+++ b/lib/include/wolff/finite_states.hpp
@@ -8,13 +8,15 @@
 
 #include "system.hpp"
 
+namespace wolff {
+
 std::array<std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> finite_states_Zp;
 #ifndef WOLFF_NO_FIELD
 std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> finite_states_Bp;
 #endif
 
 template <class R_t, class X_t>
-void finite_states_init(const wolff_system<R_t, X_t>& S) {
+void finite_states_init(const system<R_t, X_t>& S) {
 #ifndef WOLFF_NO_FIELD
   for (q_t i = 0; i < WOLFF_FINITE_STATES_N; i++) {
     for (q_t j = 0; j < WOLFF_FINITE_STATES_N; j++) {
@@ -31,5 +33,7 @@ void finite_states_init(const wolff_system<R_t, X_t>& S) {
   }
 }
 
+}
+
 #endif
 
diff --git a/lib/include/wolff/graph.hpp b/lib/include/wolff/graph.hpp
index fcd73fa..65b3941 100644
--- a/lib/include/wolff/graph.hpp
+++ b/lib/include/wolff/graph.hpp
@@ -5,23 +5,26 @@
 #include <cmath>
 #include <vector>
 
-#include "types.h"
+namespace wolff {
 
-typedef enum lattice_t {
-  SQUARE_LATTICE,
-  DIAGONAL_LATTICE
-} lattice_t;
+#include "types.h"
 
-class graph_t {
+class graph {
   public:
+    D_t D;
+    L_t L;
     v_t ne;
     v_t nv;
-    std::vector<std::vector<v_t>> adj;
+    std::vector<std::vector<v_t>> adjacency;
     std::vector<std::vector<double>> coordinate;
 
-    graph_t(D_t D, L_t L, lattice_t lat = SQUARE_LATTICE);
-    void add_ext();
+    graph();
+    graph(D_t D, L_t L);
+
+    void add_ghost();
 };
 
+}
+
 #endif
 
diff --git a/lib/include/wolff/measurement.hpp b/lib/include/wolff/measurement.hpp
index a6a5f79..604eaf5 100644
--- a/lib/include/wolff/measurement.hpp
+++ b/lib/include/wolff/measurement.hpp
@@ -4,21 +4,25 @@
 
 #include "system.hpp"
 
+namespace wolff {
+
 template <class R_t, class X_t>
-class wolff_measurement {
+class measurement {
   public:
-    virtual void pre_cluster(N_t, N_t, const wolff_system<R_t, X_t>&, v_t, const R_t&) = 0;
+    virtual void pre_cluster(N_t, N_t, const system<R_t, X_t>&, v_t, const R_t&) = 0;
 
-    virtual void plain_bond_visited(const wolff_system<R_t, X_t>&, v_t, const X_t&, v_t, double) = 0;
-    virtual void plain_site_transformed(const wolff_system<R_t, X_t>&, v_t, const X_t&) = 0;
+    virtual void plain_bond_visited(const system<R_t, X_t>&, v_t, const X_t&, v_t, double) = 0;
+    virtual void plain_site_transformed(const system<R_t, X_t>&, v_t, const X_t&) = 0;
 
 #ifndef WOLFF_NO_FIELD
-    virtual void ghost_bond_visited(const wolff_system<R_t, X_t>&, v_t, const X_t&, const X_t&, double) = 0;
-    virtual void ghost_site_transformed(const wolff_system<R_t, X_t>&, const R_t&) = 0;
+    virtual void ghost_bond_visited(const system<R_t, X_t>&, v_t, const X_t&, const X_t&, double) = 0;
+    virtual void ghost_site_transformed(const system<R_t, X_t>&, const R_t&) = 0;
 #endif
 
-    virtual void post_cluster(N_t, N_t, const wolff_system<R_t, X_t>&) = 0;
+    virtual void post_cluster(N_t, N_t, const system<R_t, X_t>&) = 0;
 };
 
+}
+
 #endif
 
diff --git a/lib/include/wolff/models/ising.hpp b/lib/include/wolff/models/ising.hpp
index c6a3a47..824c063 100644
--- a/lib/include/wolff/models/ising.hpp
+++ b/lib/include/wolff/models/ising.hpp
@@ -3,6 +3,9 @@
 #define WOLFF_MODELS_ISING
 
 #include "../types.h"
+#include "../system.hpp"
+
+namespace wolff {
 
 class ising_t {
   public:
@@ -73,9 +76,15 @@ inline ising_t::F_t operator*(double a, const ising_t& s) {
   return s * a;
 }
 
+ising_t gen_ising(std::mt19937&, const system<ising_t, ising_t>&, v_t) {
+  return ising_t(true);
+};
+
 #define WOLFF_FINITE_STATES_N 2
 const ising_t finite_states_possible[2] = {ising_t(0), ising_t(1)};
-q_t finite_states_enum(ising_t state) { return (q_t)state.x; }
+q_t finite_states_enum(const ising_t& state) { return (q_t)state.x; }
+
+}
 
 #endif
 
diff --git a/lib/include/wolff/models/orthogonal.hpp b/lib/include/wolff/models/orthogonal.hpp
index 9dd5ddd..58203b7 100644
--- a/lib/include/wolff/models/orthogonal.hpp
+++ b/lib/include/wolff/models/orthogonal.hpp
@@ -4,7 +4,8 @@
 #include <random>
 #include <cmath>
 
-#include <wolff/types.h>
+#include "../types.h"
+#include "../system.hpp"
 #include "vector.hpp"
 
 template <q_t q, class T>
@@ -113,7 +114,7 @@ class orthogonal_t : public std::array<std::array<T, q>, q> {
 };
 
 template <q_t q>
-orthogonal_t <q, double> generate_rotation_uniform (std::mt19937& r, const vector_t <q, double>& v) {
+orthogonal_t <q, double> generate_rotation_uniform (std::mt19937& r, const system<orthogonal_t<q, double>, vector_t<q, double>>&, v_t) {
   std::normal_distribution<double> dist(0.0,1.0);
   orthogonal_t <q, double> ptr;
   ptr.is_reflection = true;
@@ -135,7 +136,7 @@ orthogonal_t <q, double> generate_rotation_uniform (std::mt19937& r, const vecto
 }
 
 template <q_t q>
-orthogonal_t <q, double> generate_rotation_perturbation (std::mt19937& r, const vector_t <q, double>& v0, double epsilon, unsigned int n) {
+orthogonal_t <q, double> generate_rotation_perturbation (std::mt19937& r, const system<orthogonal_t<q, double>, vector_t<q, double>>& S, v_t i0, double epsilon, unsigned int n) {
   std::normal_distribution<double> dist(0.0,1.0);
   orthogonal_t <q, double> m;
   m.is_reflection = true;
@@ -149,14 +150,14 @@ orthogonal_t <q, double> generate_rotation_perturbation (std::mt19937& r, const
     double cosr = cos(2 * M_PI * rotation / (double)n / 2.0);
     double sinr = sin(2 * M_PI * rotation / (double)n / 2.0);
 
-    v[0] = v0[0] * cosr - v0[1] * sinr;
-    v[1] = v0[1] * cosr + v0[0] * sinr;
+    v[0] = S.s[i0][0] * cosr - S.s[i0][1] * sinr;
+    v[1] = S.s[i0][1] * cosr + S.s[i0][0] * sinr;
 
     for (q_t i = 2; i < q; i++) {
-      v[i] = v0[i];
+      v[i] = S.s[i0][i];
     }
   } else {
-    v = v0;
+    v = S.s[i0];
   }
 
   double m_dot_v = 0;
diff --git a/lib/include/wolff/models/potts.hpp b/lib/include/wolff/models/potts.hpp
index 903f25f..10b440e 100644
--- a/lib/include/wolff/models/potts.hpp
+++ b/lib/include/wolff/models/potts.hpp
@@ -46,5 +46,5 @@ class potts_t {
 
 const potts_t<POTTSQ> states[256] = {{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {9}, {11}, {12}, {13}, {14}, {15}, {16}, {17}, {18}, {19}, {20}, {21}, {22}, {23}, {24}, {25}, {26}, {27}, {28}, {29}, {30}, {31}, {32}, {33}, {34}, {35}, {36}, {37}, {38}, {39}, {40}, {41}, {42}, {43}, {44}, {45}, {46}, {47}, {48}, {49}, {50}, {51}, {52}, {53}, {54}, {55}, {56}, {57}, {58}, {59}, {60}, {61}, {62}, {63}, {64}, {65}, {66}, {67}, {68}, {69}, {70}, {71}, {72}, {73}, {74}, {75}, {76}, {77}, {78}, {79}, {80}, {81}, {82}, {83}, {84}, {85}, {86}, {87}, {88}, {89}, {90}, {91}, {92}, {93}, {94}, {95}, {96}, {97}, {98}, {99}, {100}, {101}, {102}, {103}, {104}, {105}, {106}, {107}, {108}, {109}, {110}, {111}, {112}, {113}, {114}, {115}, {116}, {117}, {118}, {119}, {120}, {121}, {122}, {123}, {124}, {125}, {126}, {127}, {128}, {129}, {130}, {131}, {132}, {133}, {134}, {135}, {136}, {137}, {138}, {139}, {140}, {141}, {142}, {143}, {144}, {145}, {146}, {147}, {148}, {149}, {150}, {151}, {152}, {153}, {154}, {155}, {156}, {157}, {158}, {159}, {160}, {161}, {162}, {163}, {164}, {165}, {166}, {167}, {168}, {169}, {170}, {171}, {172}, {173}, {174}, {175}, {176}, {177}, {178}, {179}, {180}, {181}, {182}, {183}, {184}, {185}, {186}, {187}, {188}, {189}, {190}, {191}, {192}, {193}, {194}, {195}, {196}, {197}, {198}, {199}, {200}, {201}, {202}, {203}, {204}, {205}, {206}, {207}, {208}, {209}, {210}, {211}, {212}, {213}, {214}, {215}, {216}, {217}, {218}, {219}, {220}, {221}, {222}, {223}, {224}, {225}, {226}, {227}, {228}, {229}, {230}, {231}, {232}, {233}, {234}, {235}, {236}, {237}, {238}, {239}, {240}, {241}, {242}, {243}, {244}, {245}, {246}, {247}, {248}, {249}, {250}, {251}, {252}, {253}, {254}, {255}};
 template <q_t q>
-q_t finite_states_enum(potts_t<q> state) { return (q_t)state.x; }
+q_t finite_states_enum(const potts_t<q>& state) { return (q_t)state.x; }
 
diff --git a/lib/include/wolff/system.hpp b/lib/include/wolff/system.hpp
index ac82f44..32ad38e 100644
--- a/lib/include/wolff/system.hpp
+++ b/lib/include/wolff/system.hpp
@@ -4,23 +4,25 @@
 
 #include <functional>
 #include <vector>
+#include <random>
 
-#include "types.h"
 #include "graph.hpp"
 
+namespace wolff {
+
+#include "types.h"
+
+template <class R_t, class X_t>
+class measurement;
+
 template <class R_t, class X_t>
-class wolff_system {
+class system {
   public:
-    D_t D; // dimension
-    L_t L; // linear size
     v_t nv; // number of vertices
     v_t ne; // number of edges
-    graph_t G; // the graph defining the lattice with ghost
+    graph G; // the graph defining the lattice with ghost
     double T; // the temperature
     std::vector<X_t> s; // the state of the ordinary spins
-#ifndef WOLFF_NO_FIELD
-    R_t s0; // the current state of the ghost site
-#endif
 
 #ifdef WOLFF_BOND_DEPENDENCE
     std::function <double(v_t, const X_t&, v_t, const X_t&)> Z; // coupling between sites
@@ -29,6 +31,7 @@ class wolff_system {
 #endif
 
 #ifndef WOLFF_NO_FIELD
+    R_t s0; // the current state of the ghost site
 #ifdef WOLFF_SITE_DEPENDENCE
     std::function <double(v_t, const X_t&)> B; // coupling with the external field
 #else
@@ -36,7 +39,7 @@ class wolff_system {
 #endif
 #endif
 
-    wolff_system(D_t D, L_t L, double T,
+    system(graph g, double T,
 #ifdef WOLFF_BOND_DEPENDENCE
         std::function <double(v_t, const X_t&, v_t, const X_t&)> Z
 #else
@@ -49,7 +52,7 @@ class wolff_system {
         , std::function <double(const X_t&)> B
 #endif
 #endif
-        , lattice_t lat = SQUARE_LATTICE) : D(D), L(L), G(D, L, lat), T(T), Z(Z)
+        ) : G(g), T(T), Z(Z)
 #ifndef WOLFF_NO_FIELD
              , s0(), B(B)
 #endif
@@ -58,13 +61,18 @@ class wolff_system {
       ne = G.ne;
       s.resize(nv);
 #ifndef WOLFF_NO_FIELD
-      G.add_ext();
+      G.add_ghost();
 #endif
 #ifdef WOLFF_FINITE_STATES
       finite_states_init(*this);
 #endif
     }
+
+    void flip_cluster(v_t, const R_t&, std::mt19937&, measurement<R_t, X_t>&);
+    void run_wolff(N_t, std::function <R_t(std::mt19937&, const system<R_t, X_t>&, v_t)> r_gen, measurement<R_t, X_t>& A, std::mt19937& rng);
 };
 
+}
+
 #endif
 
diff --git a/lib/include/wolff/types.h b/lib/include/wolff/types.h
index 68bd09e..1bc8c4d 100644
--- a/lib/include/wolff/types.h
+++ b/lib/include/wolff/types.h
@@ -8,7 +8,7 @@ typedef uint_fast16_t L_t;     // linear size
 typedef uint_fast64_t N_t;     // cycle iterator
 
 #define MAX_v UINT_FAST32_MAX
-#define MAX_Q UINT_FAST8_MAX
+#define MAX_q UINT_FAST8_MAX
 #define MAX_D UINT_FAST8_MAX
 #define MAX_L UINT_FAST16_MAX
 #define MAX_N UINT_FAST64_MAX
diff --git a/lib/src/graph.cpp b/lib/src/graph.cpp
index b42aa78..4f91be4 100644
--- a/lib/src/graph.cpp
+++ b/lib/src/graph.cpp
@@ -1,74 +1,55 @@
 
 #include <wolff/graph.hpp>
 
-graph_t::graph_t(D_t D, L_t L, lattice_t lat) {
-  switch (lat) {
-    case SQUARE_LATTICE: {
-        nv = pow(L, D);
-        ne = D * nv;
+namespace wolff {
 
-        adj.resize(nv);
-        coordinate.resize(nv);
-
-        for (std::vector<v_t> adj_i : adj) {
-          adj_i.reserve(2 * D);
-        }
-
-        for (v_t i = 0; i < nv; i++) {
-          coordinate[i].resize(D);
-          for (D_t j = 0; j < D; j++) {
-            coordinate[i][j] = (i / (v_t)pow(L, D - j - 1)) % L;
-
-            adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1)));
-            adj[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(pow(L, j+1) + i - pow(L, j), pow(L, j + 1)));
-          }
-        }
-        break;
-      }
-    case DIAGONAL_LATTICE: {
-        nv = D * pow(L, D);
-        ne = 2 * nv;
+graph::graph() {
+  D = 0;
+  L = 0;
+  nv = 0;
+  ne = 0;
+}
 
-        adj.resize(nv);
-        coordinate.resize(nv);
+graph::graph(D_t D, L_t L) : D(D), L(L) {
+  nv = pow(L, D);
+  ne = D * nv;
 
-        for (std::vector<v_t> adj_i : adj) {
-          adj_i.reserve(4 * (D - 1));
-        }
+  adjacency.resize(nv);
+  coordinate.resize(nv);
 
-        for (D_t i = 0; i < D; i++) {
-          v_t sb = i * pow(L, D);
+  for (std::vector<v_t> adj_i : adjacency) {
+    adj_i.reserve(2 * D);
+  }
 
-          for (v_t j = 0; j < pow(L, D); j++) {
-            v_t vc = sb + j;
+  for (v_t i = 0; i < nv; i++) {
+    coordinate[i].resize(D);
+    for (D_t j = 0; j < D; j++) {
+      coordinate[i][j] = (i / (v_t)pow(L, D - j - 1)) % L;
 
-            adj[vc].push_back(((i + 1) % D) * pow(L, D) + j);
-            adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * (j / (v_t)pow(L, D - 1)) + (j + 1 - 2 * (i % 2)) % L);
-            adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j - i) % L);
-            adj[vc].push_back(((i + 1) % D) * pow(L, D) + pow(L, D - 1) * ((L + (j/ (v_t)pow(L, D - 1)) - 1 + 2 * (i % 2)) % L) + (j + 1 - i) % L);
-          }
-        }
-        break;
-      }
+      adjacency[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1)));
+      adjacency[i].push_back(pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(pow(L, j+1) + i - pow(L, j), pow(L, j + 1)));
+    }
   }
 }
 
-void graph_t::add_ext() {
-  for (std::vector<v_t>& adj_i : adj) {
+void graph::add_ghost() {
+  for (std::vector<v_t>& adj_i : adjacency) {
     adj_i.push_back(nv);
   }
 
-  adj.resize(nv + 1);
+  adjacency.resize(nv + 1);
   coordinate.resize(nv + 1);
-  adj[nv].reserve(nv);
+  adjacency[nv].reserve(nv);
 
   for (v_t i = 0; i < nv; i++) {
-    adj[nv].push_back(i);
+    adjacency[nv].push_back(i);
   }
 
   coordinate[nv].resize(coordinate[0].size());
 
   ne += nv;
-  nv += 1;
+  nv++;
+}
+
 }