#pragma once #include #include #include #include "spin.hpp" #include "vector.hpp" namespace std { template struct hash> { typedef array argument_type; typedef size_t result_type; result_type operator()(const argument_type& a) const { hash hasher; result_type h = 0; for (result_type i = 0; i < N; ++i) { h = h * 31 + hasher(a[i]); } return h; } }; } // namespace std template class Octree { private: unsigned N; T L; std::unordered_map, std::set*>> data; public: Octree(T L_tmp, unsigned N_tmp) { L = L_tmp; N = N_tmp; } std::array ind(Vector x) const { std::array ind; for (unsigned i = 0; i < D; i++) { ind[i] = std::floor(N * x(i) / L); } return ind; } void insert(Spin* s) { data[ind(s->x)].insert(s); }; void remove(Spin* s) { data[ind(s->x)].erase(s); if (data[ind(s->x)].empty()) { data.erase(ind(s->x)); } }; std::set*> at(const Vector& x) const { auto it = data.find(ind(x)); if (it == data.end()) { return {}; } else { return it->second; } } std::set*> neighbors(const Vector& x) const { std::array i0 = ind(x); std::set*> ns; nearest_neighbors_of(i0, D + 1, ns); return ns; }; void nearest_neighbors_of(std::array i0, unsigned depth, std::set*>& ns) const { if (depth == 0) { auto it = data.find(i0); if (it != data.end()) { ns.insert(it->second.begin(), it->second.end()); } } else { for (signed j : {-1, 0, 1}) { std::array i1 = i0; if (N < 2) { i1[depth - 1] += j; } else { i1[depth - 1] = (N + i1[depth - 1] + j) % N; } nearest_neighbors_of(i1, depth - 1, ns); } } }; };