Refactor: Add SIMD policy option to all templates

src/include/lattice.h

@@ -13,7 +13,7 @@
 
 namespace pivot {
 
-template <int Dim> struct box;
+template <int Dim, bool Simd> struct box;
 
 /**
  * @brief Represents a closed interval [left_, right_].
@@ -39,7 +39,7 @@ struct interval {
 /**
  * @brief Represents a Dim-dimensional point.
  */
-template <int Dim> class point : boost::multipliable<point<Dim>, int> {
+template <int Dim, bool Simd = false> class point : boost::multipliable<point<Dim, Simd>, int> {
 
 public:
   point() = default;
@@ -77,10 +77,10 @@ template <int Dim> class point : boost::multipliable<point<Dim>, int> {
   std::array<int, Dim> coords_{};
 };
 
-template <typename S, typename T, T Dim> point(std::array<S, Dim>) -> point<Dim>;
+template <typename S, typename T, T Dim> point(std::array<S, Dim>) -> point<Dim, false>;
 
 /** @brief Represents a Dim-dimensional box. */
-template <int Dim> struct box : boost::additive<box<Dim>, point<Dim>> {
+template <int Dim, bool Simd = false> struct box : boost::additive<box<Dim, Simd>, point<Dim, Simd>> {
   std::array<interval, Dim> intervals_;
 
   box() = delete;
@@ -94,7 +94,7 @@ template <int Dim> struct box : boost::additive<box<Dim>, point<Dim>> {
    * in the sense that the input sequences of points (p0, ...) is effectively translated by e0 - p0
    * prior to the box's construction. In other words, the box will always have a vertex at e0.
    */
-  box(std::span<const point<Dim>> points);
+  box(std::span<const point<Dim, Simd>> points);
 
   bool operator==(const box &b) const;
 
@@ -105,9 +105,9 @@ template <int Dim> struct box : boost::additive<box<Dim>, point<Dim>> {
   bool empty() const;
 
   /** @brief Action of a point (understood as a vector) on a box */
-  box<Dim> &operator+=(const point<Dim> &b);
+  box<Dim, Simd> &operator+=(const point<Dim, Simd> &b);
 
-  box<Dim> &operator-=(const point<Dim> &b);
+  box<Dim, Simd> &operator-=(const point<Dim, Simd> &b);
 
   /**
    * @brief Returns the "union" of two boxes.
@@ -129,7 +129,7 @@ template <int Dim> struct box : boost::additive<box<Dim>, point<Dim>> {
   std::string to_string() const;
 };
 
-template <typename S, typename T, S Dim, T N> box(std::array<point<Dim>, N>) -> box<Dim>;
+template <typename S, typename T, S Dim, T N, bool Simd> box(std::array<point<Dim, Simd>, N>) -> box<Dim, Simd>;
 
 struct point_hash {
   int num_steps_;
@@ -138,7 +138,7 @@ struct point_hash {
 
   // This hashing method, which exploits the known range of values that can be taken by the
   // sequence of points in a walk, ppears to result in better performance than other methods tested.
-  template <int Dim> std::size_t operator()(const point<Dim> &p) const;
+  template <int Dim, bool Simd = false> std::size_t operator()(const point<Dim, Simd> &p) const;
 };
 
 /**
@@ -156,7 +156,7 @@ struct point_hash {
  * entries of S by S(i). Note that P e(i) = e(P(i)) and S e(i) = S(i) e(i), where e(i)
  * is the i-th standard unit vector.
  */
-template <int Dim> class transform {
+template <int Dim, bool Simd = false> class transform {
 
 public:
   /** @brief Constructs the identity transformation. */
@@ -174,7 +174,7 @@ template <int Dim> class transform {
    *
    * Note that the resulting pivot transformation is not uniquely defined by the inputs.
    */
-  transform(const point<Dim> &p, const point<Dim> &q);
+  transform(const point<Dim, Simd> &p, const point<Dim, Simd> &q);
 
   /** @brief Produce a uniformly random transfom.*/
   template <typename Gen> static transform rand(Gen &gen) {
@@ -203,7 +203,7 @@ template <int Dim> class transform {
    *
    * which means the P(i)-th component of the result is S(P(i)) p[i].
    */
-  point<Dim> operator*(const point<Dim> &p) const;
+  point<Dim, Simd> operator*(const point<Dim, Simd> &p) const;
 
   /**
    * @brief Composes two transforms.
@@ -238,7 +238,7 @@ template <int Dim> class transform {
    * which means the P(i)-th interval from which SP(B) is constructed has bounds S(P(i)) a[i] and
    * S(P(i)) b[i].
    */
-  box<Dim> operator*(const box<Dim> &b) const;
+  box<Dim, Simd> operator*(const box<Dim, Simd> &b) const;
 
   /** @brief Returns true if the transform is the identity. */
   bool is_identity() const;

src/include/utils.h

@@ -4,10 +4,10 @@
 
 namespace pivot {
 
-template <int Dim> std::vector<point<Dim>> from_csv(const std::string &path);
+template <int Dim, bool Simd> std::vector<point<Dim, Simd>> from_csv(const std::string &path);
 
-template <int Dim> void to_csv(const std::string &path, const std::vector<point<Dim>> &points);
+template <int Dim, bool Simd> void to_csv(const std::string &path, const std::vector<point<Dim, Simd>> &points);
 
-template <int Dim> std::vector<point<Dim>> line(int num_steps);
+template <int Dim, bool Simd> std::vector<point<Dim, Simd>> line(int num_steps);
 
 } // namespace pivot

src/include/walk.h

@@ -11,10 +11,10 @@
 
 namespace pivot {
 
-template <int Dim> class walk : public walk_base<Dim> {
+template <int Dim, bool Simd = false> class walk : public walk_base<Dim, Simd> {
 
 public:
-  walk(const std::vector<point<Dim>> &steps, std::optional<unsigned int> seed = std::nullopt);
+  walk(const std::vector<point<Dim, Simd>> &steps, std::optional<unsigned int> seed = std::nullopt);
   walk(int num_steps, std::optional<unsigned int> seed = std::nullopt);
   walk(const std::string &path, std::optional<unsigned int> seed = std::nullopt);
 
@@ -24,15 +24,15 @@ template <int Dim> class walk : public walk_base<Dim> {
 
   ~walk() = default;
 
-  point<Dim> operator[](int i) const { return steps_[i]; }
+  point<Dim, Simd> operator[](int i) const { return steps_[i]; }
 
   int num_steps() const { return steps_.size(); }
 
-  point<Dim> endpoint() const override { return steps_.back(); }
+  point<Dim, Simd> endpoint() const override { return steps_.back(); }
 
-  std::optional<std::vector<point<Dim>>> try_pivot(int step, const transform<Dim> &trans) const;
+  std::optional<std::vector<point<Dim, Simd>>> try_pivot(int step, const transform<Dim, Simd> &trans) const;
 
-  std::pair<int, std::optional<std::vector<point<Dim>>>> try_rand_pivot() const;
+  std::pair<int, std::optional<std::vector<point<Dim, Simd>>>> try_rand_pivot() const;
 
   bool rand_pivot(bool fast = false) override;
 
@@ -43,15 +43,15 @@ template <int Dim> class walk : public walk_base<Dim> {
   void export_csv(const std::string &path) const override;
 
 protected:
-  std::vector<point<Dim>> steps_;
-  boost::unordered_flat_map<point<Dim>, int, point_hash> occupied_;
+  std::vector<point<Dim, Simd>> steps_;
+  boost::unordered_flat_map<point<Dim, Simd>, int, point_hash> occupied_;
 
   mutable std::mt19937 rng_;
   mutable std::uniform_int_distribution<int> dist_;
 
-  void do_pivot(int step, std::vector<point<Dim>> &new_points);
+  void do_pivot(int step, std::vector<point<Dim, Simd>> &new_points);
 
-  point<Dim> pivot_point(int step, int i, const transform<Dim> &trans) const;
+  point<Dim, Simd> pivot_point(int step, int i, const transform<Dim, Simd> &trans) const;
 };
 
 } // namespace pivot

src/include/walk_base.h

@@ -4,7 +4,7 @@
 
 namespace pivot {
 
-template <int Dim> class walk_base {
+template <int Dim, bool Simd = false> class walk_base {
 
 public:
   virtual ~walk_base() = default;
@@ -15,7 +15,7 @@ template <int Dim> class walk_base {
 
   virtual void export_csv(const std::string &path) const = 0;
 
-  virtual point<Dim> endpoint() const = 0;
+  virtual point<Dim, Simd> endpoint() const = 0;
 };
 
 } // namespace pivot

src/include/walk_node.h

@@ -10,18 +10,19 @@ namespace pivot {
 
 /* FORWARD REFERENCES */
 
-template <int Dim> class walk_node;
+template <int Dim, bool Simd> class walk_node;
 
-template <int Dim>
-bool intersect(const walk_node<Dim> *l_walk, const walk_node<Dim> *r_walk, const point<Dim> &l_anchor,
-               const point<Dim> &r_anchor, const transform<Dim> &l_symm, const transform<Dim> &r_symm);
+template <int Dim, bool Simd = false>
+bool intersect(const walk_node<Dim, Simd> *l_walk, const walk_node<Dim, Simd> *r_walk, const point<Dim, Simd> &l_anchor,
+               const point<Dim, Simd> &r_anchor, const transform<Dim, Simd> &l_symm,
+               const transform<Dim, Simd> &r_symm);
 
-template <int Dim> class walk_tree;
+template <int Dim, bool Simd> class walk_tree;
 
 /* WALK NODE */
 
 /** @brief Represents a node in a walk tree. */
-template <int Dim> class walk_node {
+template <int Dim, bool Simd = false> class walk_node {
 
 public:
   /* CONSTRUCTORS, DESTRUCTOR */
@@ -33,7 +34,7 @@ template <int Dim> class walk_node {
    *
    * @return The root of the walk tree.
    */
-  static walk_node *pivot_rep(const std::vector<point<Dim>> &steps, walk_node *buf = nullptr);
+  static walk_node *pivot_rep(const std::vector<point<Dim, Simd>> &steps, walk_node *buf = nullptr);
 
   /** @brief Returns the root of a walk tree for the balanced representation of a walk given by a sequence of points.
    *
@@ -42,7 +43,7 @@ template <int Dim> class walk_node {
    *
    * @return The root of the walk tree.
    */
-  static walk_node *balanced_rep(const std::vector<point<Dim>> &steps, walk_node *buf = nullptr);
+  static walk_node *balanced_rep(const std::vector<point<Dim, Simd>> &steps, walk_node *buf = nullptr);
 
   /** @brief Copies the given node but none of the nodes it links to. */
   walk_node(const walk_node &w) = default;
@@ -57,11 +58,11 @@ template <int Dim> class walk_node {
 
   int id() const { return id_; }
 
-  const box<Dim> &bbox() const { return bbox_; }
+  const box<Dim, Simd> &bbox() const { return bbox_; }
 
-  const point<Dim> &endpoint() const { return end_; }
+  const point<Dim, Simd> &endpoint() const { return end_; }
 
-  const transform<Dim> &symm() const { return symm_; }
+  const transform<Dim, Simd> &symm() const { return symm_; }
 
   walk_node *left() const { return left_; }
 
@@ -125,7 +126,7 @@ template <int Dim> class walk_node {
    *
    * @return Whether the transform creates an intersection.
    */
-  bool shuffle_intersect(const transform<Dim> &t, std::optional<bool> is_left_child);
+  bool shuffle_intersect(const transform<Dim, Simd> &t, std::optional<bool> is_left_child);
 
   /**
    * @brief Checks if the current walk has an intersection via a top-down algorithm.
@@ -136,7 +137,7 @@ template <int Dim> class walk_node {
 
   /* OTHER FUNCTIONS */
 
-  std::vector<point<Dim>> steps() const;
+  std::vector<point<Dim, Simd>> steps() const;
 
   void todot(const std::string &path) const;
 
@@ -146,15 +147,16 @@ template <int Dim> class walk_node {
   walk_node *parent_{};
   walk_node *left_{};
   walk_node *right_{};
-  transform<Dim> symm_;
-  box<Dim> bbox_;
-  point<Dim> end_;
+  transform<Dim, Simd> symm_;
+  box<Dim, Simd> bbox_;
+  point<Dim, Simd> end_;
 
-  friend class walk_tree<Dim>;
+  friend class walk_tree<Dim, Simd>;
 
   /* CONVENIENCE METHODS */
 
-  walk_node(int id, int num_sites, const transform<Dim> &symm, const box<Dim> &bbox, const point<Dim> &end);
+  walk_node(int id, int num_sites, const transform<Dim, Simd> &symm, const box<Dim, Simd> &bbox,
+            const point<Dim, Simd> &end);
 
   void set_left(walk_node *left) {
     left_ = left;
@@ -185,15 +187,15 @@ template <int Dim> class walk_node {
   Agnode_t *todot(Agraph_t *g, const cgraph_t &cgraph) const;
 
   // recursive helper
-  static walk_node *balanced_rep(std::span<const point<Dim>> steps, int start, const transform<Dim> &glob_symm,
-                                 walk_node *buf);
+  static walk_node *balanced_rep(std::span<const point<Dim, Simd>> steps, int start,
+                                 const transform<Dim, Simd> &glob_symm, walk_node *buf);
 
-  bool shuffle_intersect(const transform<Dim> &t, std::optional<bool> was_left_child,
+  bool shuffle_intersect(const transform<Dim, Simd> &t, std::optional<bool> was_left_child,
                          std::optional<bool> is_left_child);
 
-  template <int D>
-  friend bool intersect(const walk_node<D> *l_walk, const walk_node<D> *r_walk, const point<D> &l_anchor,
-                        const point<D> &r_anchor, const transform<D> &l_symm, const transform<D> &r_symm);
+  template <int D, bool S>
+  friend bool intersect(const walk_node<D, S> *l_walk, const walk_node<D, S> *r_walk, const point<D, S> &l_anchor,
+                        const point<D, S> &r_anchor, const transform<D, S> &l_symm, const transform<D, S> &r_symm);
 };
 
 } // namespace pivot

src/include/walk_tree.h

@@ -10,10 +10,10 @@
 
 namespace pivot {
 
-template <int Dim> class walk_node;
+template <int Dim, bool Simd> class walk_node;
 
 /** @brief Represents an entire saw-tree (as per Clisby's 2010 paper). */
-template <int Dim> class walk_tree : public walk_base<Dim> {
+template <int Dim, bool Simd = false> class walk_tree : public walk_base<Dim, Simd> {
 
 public:
   /* CONSTRUCTORS, DESTRUCTOR */
@@ -56,17 +56,17 @@ template <int Dim> class walk_tree : public walk_base<Dim> {
    *
    * @warning It is not recommended to set balanced=false.
    */
-  walk_tree(const std::vector<point<Dim>> &steps, std::optional<unsigned int> seed = std::nullopt,
+  walk_tree(const std::vector<point<Dim, Simd>> &steps, std::optional<unsigned int> seed = std::nullopt,
             bool balanced = true);
 
   /**@brief Deallocates the entire tree and every node it contains. */
   ~walk_tree();
 
   /* GETTERS, SETTERS, SIMPLE UTILITIES */
 
-  walk_node<Dim> *root() const;
+  walk_node<Dim, Simd> *root() const;
 
-  point<Dim> endpoint() const override;
+  point<Dim, Simd> endpoint() const override;
 
   bool is_leaf() const;
 
@@ -84,7 +84,7 @@ template <int Dim> class walk_tree : public walk_base<Dim> {
    *
    * @return reference to the node
    */
-  walk_node<Dim> &find_node(int n);
+  walk_node<Dim, Simd> &find_node(int n);
 
   /* HIGH-LEVEL FUNCTIONS (see Clisby (2010), Section 2.7) */
 
@@ -99,7 +99,7 @@ template <int Dim> class walk_tree : public walk_base<Dim> {
    *
    * @return Whether the pivot was successful.
    */
-  bool try_pivot(int n, const transform<Dim> &r);
+  bool try_pivot(int n, const transform<Dim, Simd> &r);
 
   /**
    * @brief Attempt to pivot the walk about the given lattice site using Clisby's Attempt_pivot_fast function.
@@ -112,7 +112,7 @@ template <int Dim> class walk_tree : public walk_base<Dim> {
    *
    * @return Whether the pivot was successful.
    */
-  bool try_pivot_fast(int n, const transform<Dim> &r);
+  bool try_pivot_fast(int n, const transform<Dim, Simd> &r);
 
   /**
    * @brief Attempts to pivot the walk about a randomly chosen lattice site with a random transform.
@@ -130,7 +130,7 @@ template <int Dim> class walk_tree : public walk_base<Dim> {
    *
    * @return Sequence of lattice sites.
    */
-  std::vector<point<Dim>> steps() const;
+  std::vector<point<Dim, Simd>> steps() const;
 
   /**
    * @brief Check whether the walk is self-avoiding using a naive algorithm.
@@ -149,10 +149,10 @@ template <int Dim> class walk_tree : public walk_base<Dim> {
   void todot(const std::string &path) const;
 
 private:
-  std::unique_ptr<walk_node<Dim>> root_;
+  std::unique_ptr<walk_node<Dim, Simd>> root_;
   std::mt19937 rng_;
   std::uniform_int_distribution<int> dist_; // distribution for choosing a random lattice site
-  walk_node<Dim> *buf_;                     // buffer into which nodes are allocated (used for fast node lookup by id)
+  walk_node<Dim, Simd> *buf_;               // buffer into which nodes are allocated (used for fast node lookup by id)
 };
 
 } // namespace pivot