// Licensed under the MIT License <http://opensource.org/licenses/MIT>. // SPDX-License-Identifier: MIT // Copyright (c) 2021 Noah H. and Tom H. #include <iostream> #include <numeric> #include <map> #include <string> #include <iterator> #include <queue> #include <random> #include <algorithm> #include <vector> #include <unordered_map> #include <functional> #include <tuple> #include <fstream> #include <sched.h> #include <stdlib.h> #include <stdio.h> #include <list> #include <cstdlib> #include <cstdio> #include <memory> #include <stack> #include "robin_hood.h" #define MS 662891 namespace hash_tuple { template <typename TT> struct hash { size_t operator()(TT const& tt) const { return robin_hood::hash<TT>()(tt); } }; } namespace hash_tuple{ namespace { template <class T> inline void hash_combine(std::size_t& seed, T const& v) { seed ^= hash_tuple::hash<T>()(v) + 0x9e3779b9 + (seed<<6) + (seed>>2); } } } namespace hash_tuple{ namespace { // Recursive template code derived from Matthieu M. template <class Tuple, size_t Index = std::tuple_size<Tuple>::value - 1> struct HashValueImpl { static void apply(size_t& seed, Tuple const& tuple) { HashValueImpl<Tuple, Index-1>::apply(seed, tuple); hash_combine(seed, std::get<Index>(tuple)); } }; template <class Tuple> struct HashValueImpl<Tuple,0> { static void apply(size_t& seed, Tuple const& tuple) { hash_combine(seed, std::get<0>(tuple)); } }; } template <typename ... TT> struct hash<std::tuple<TT...>> { size_t operator()(std::tuple<TT...> const& tt) const { size_t seed = 0; HashValueImpl<std::tuple<TT...> >::apply(seed, tt); return seed; } }; } using namespace std; typedef tuple<long long,long long,long long,long long,long long> Vertex; void write_to_file(vector<tuple<int,int,double>> &vec, string file_name) { std::ofstream f(file_name); for(vector<tuple<int,int,double>>::const_iterator i = vec.begin(); i != vec.end(); ++i) { f << get<0>(*i) << ','<< get<1>(*i) << ','<<get<2>(*i) << '\n'; } } class Node { public: Vertex v; Vertex temp_coords; Node *neighbours[6] = {NULL,NULL,NULL,NULL,NULL,NULL}; Node *from = NULL; unsigned long long int counter = 0; int index=0; int is_leaf() { int count = -1; for (int ii=0;ii<8;ii++) { if (neighbours[ii]!=NULL) { if(count!=-1) { return -1; } count=ii; } } return count; } vector<Vertex> coord_neighbours() { long long x1 = get<0>(v); long long x2 = get<1>(v); long long x3 = get<2>(v); long long x4 = get<3>(v); long long x5 = get<4>(v); vector <Vertex> ret{make_tuple(x1,x2+1,x3,x4,x5), make_tuple(x1+x2,x2,x3,x4+1,x5), make_tuple(x1,x2,x3-x4,x4,x5+1), make_tuple(x1,x2,x3+x4,x4,x5-1), make_tuple(x1-x2,x2,x3,x4-1,x5), make_tuple(x1,x2-1,x3,x4,x5)}; return ret; } Vertex coord_neighbour(int i) { long long x1 = get<0>(v); long long x2 = get<1>(v); long long x3 = get<2>(v); long long x4 = get<3>(v); long long x5 = get<4>(v); switch(i) { case 0: return make_tuple(x1,x2+1,x3,x4,x5);break; case 1: return make_tuple(x1+x2,x2,x3,x4+1,x5);break; case 2: return make_tuple(x1,x2,x3-x4,x4,x5+1);break; case 3: return make_tuple(x1,x2,x3+x4,x4,x5-1);break; case 4: return make_tuple(x1-x2,x2,x3,x4-1,x5);break; default: return make_tuple(x1,x2-1,x3,x4,x5);break; } } }; class circular_buffer { private: Node* buf_[MS]; size_t head_ = 0; size_t tail_ = 0; const size_t max_size_ = MS; bool full_ = 0; public: void put(Node* item); Node* get(); void reset(); bool empty() const; bool full() const; size_t capacity(); size_t size() const; }; void circular_buffer::put(Node* item) { buf_[head_] = item; if(full_) { if(++tail_>=max_size_) { tail_ -= max_size_; } } if(++head_>=max_size_) { head_ -= max_size_; } full_ = head_ == tail_; } Node* circular_buffer::get() { if(empty()) { return NULL; } //Read data and advance the tail (we now have a free space) Node* val = buf_[tail_]; full_ = false; if(++tail_>=max_size_) { tail_ -= max_size_; } return val; } void circular_buffer::reset() { head_ = tail_; full_ = false; } bool circular_buffer::empty() const { //if head and tail are equal, we are empty return (!full_ && (head_ == tail_)); } bool circular_buffer::full() const { //If tail is ahead the head by 1, we are full return full_; } size_t circular_buffer::capacity() { return max_size_; } size_t circular_buffer::size() const { size_t size = max_size_; if(!full_) { if(head_ >= tail_) { size = head_ - tail_; } else { size = max_size_ + head_ - tail_; } } return size; } class circular_buffer_short { private: int buf_[MS]; size_t head_ = 0; size_t tail_ = 0; const size_t max_size_ = MS; bool full_ = 0; public: void put(int item); int get(); void reset(); bool empty() const; bool full() const; size_t capacity(); size_t size() const; int operator [](size_t); }; int circular_buffer_short::operator [](size_t i) { size_t get = i +tail_; if(get>=max_size_) { get -= max_size_; } return buf_[get]; } void circular_buffer_short::put(int item) { buf_[head_] = item; if(full_) { if(++tail_>=max_size_) { tail_ -= max_size_; } } if(++head_>=max_size_) { head_ -= max_size_; } full_ = head_ == tail_; } int circular_buffer_short::get() { if(empty()) { return -1; } //Read data and advance the tail (we now have a free space) int val = buf_[tail_]; full_ = false; if(++tail_>=max_size_) { tail_ -= max_size_; } return val; } void circular_buffer_short::reset() { head_ = tail_; full_ = false; } bool circular_buffer_short::empty() const { //if head and tail are equal, we are empty return (!full_ && (head_ == tail_)); } bool circular_buffer_short::full() const { //If tail is ahead the head by 1, we are full return full_; } size_t circular_buffer_short::capacity() { return max_size_; } size_t circular_buffer_short::size() const { size_t size = max_size_; if(!full_) { if(head_ >= tail_) { size = head_ - tail_; } else { size = max_size_ + head_ - tail_; } } return size; } std::function<Vertex (Vertex&)> from_to(Vertex &source, Vertex &dest) { long long x1 = get<0>(source); long long x2 = get<1>(source); long long x3 = get<2>(source); long long x4 = get<3>(source); long long x5 = get<4>(source); long long y1 = get<0>(dest); long long y2 = get<1>(dest); long long y3 = get<2>(dest); long long y4 = get<3>(dest); long long y5 = get<4>(dest); long long a = x2-y2; long long b = x4-y4; return [x1,x2,x3,x4,x5,y1,y2,y3,y4,y5,a,b](Vertex& v) { long long z1 = get<0>(v); long long z2 = get<1>(v); long long z3 = get<2>(v); long long z4 = get<3>(v); long long z5 = get<4>(v); long long c = x5 - z5; return make_tuple(-x1 + y1 + z1 + a*(x4 - z4), z2 - a, -x3 + y3 + z3 - b*c, z4 - b, y5 - c); }; } class Graph { public: vector<Node> nodes; unsigned long long int counter = 0; std::random_device rd; std::mt19937 rng; std::mt19937_64 rng64; circular_buffer *Q1 = new circular_buffer; circular_buffer *Q2 = new circular_buffer; circular_buffer_short *Q3 = new circular_buffer_short; circular_buffer_short *Q4 = new circular_buffer_short; robin_hood::unordered_map <Vertex,int,hash_tuple::hash<Vertex>> coords_to_node; ~Graph() { delete Q1; delete Q2; delete Q3; delete Q4; } Graph(int size): rng(rd()),rng64(rd()) { nodes.resize(size); coords_to_node.reserve(size); int N = size%2==0 ? size/2 : (size+1)/2; int N_ = size%2==0 ? N : N-1; for (int i=0;i<N;i++) { Vertex coord = make_tuple(0,0,0,i,0); nodes[i].v = coord; nodes[i].index=i; coords_to_node[coord] = i; nodes[i].neighbours[1] = &(nodes[i+1]); nodes[i+1].neighbours[4] = &(nodes[i]); } Vertex coord = make_tuple(0,0,0,N,0); nodes[N].v = coord; nodes[N].index = N; coords_to_node[coord]=N; for (int i=1;i<N_;i++) { Vertex coord = make_tuple(0,1,0,i,0); nodes[N+i].v = coord; nodes[N+i].index = N+i; coords_to_node[coord] = N+i; nodes[N+i].neighbours[5] = &(nodes[i]); nodes[i].neighbours[0] = &(nodes[N+i]); } } void increment_counter(int num) { if (counter>= numeric_limits<unsigned long int>::max() - num - 1) { counter = 0; for(std::vector<Node>::iterator node = nodes.begin(); node != nodes.end(); ++node) { node->counter = 0; } } counter = counter + num; } void run_algorithm(int num_samples,int steps_between_samples,int initial_delay,int samples_per_file,string file_name) { long long int so_far = -initial_delay; long long int num_steps = (long long)num_samples*(long long)steps_between_samples+1; int so_far_samples = 0; int siz = nodes.size(); long long mult = 6*siz; std::uniform_int_distribution<long long> uni_nodes(0,(long long)mult*(long long)mult-1); vector<tuple<int,int,double>>ret; ret.resize(samples_per_file); int swap_count = 0; int cut_and_paste_count = 0; int s_and_i_count = 0; int swap_suc_count = 0; int cut_and_paste_suc_count = 0; int k = 0; int z = 0; while(so_far<num_steps) { int div =so_far/steps_between_samples; if(so_far>=0 && div*steps_between_samples==so_far) { if(so_far>0 && k==samples_per_file) { cout<<div<<"/"<<num_samples<<"\n"; write_to_file(ret, file_name+"_"+to_string(z)); k=0; z++; } pair<int,int> p = longest_path(); double b = average_branching(); ret[k] = make_tuple(p.first,p.second,b); so_far_samples+=1; k++; } so_far+=1; long long random_tot = uni_nodes(rng64); int random_s = random_tot/mult; int node_index = random_s%siz; int direction = random_s/siz; Node &node = nodes[node_index]; Node *neighbour = node.neighbours[direction]; if(neighbour != NULL) { cut_and_paste_count+=1; cut_and_paste_suc_count+= cut_and_paste(&node,direction); continue; } Vertex n_coords = node.coord_neighbour(direction); auto n_index_it = coords_to_node.find(n_coords); if(n_index_it != coords_to_node.end()) { int n_index = n_index_it->second; Node *neighbour = &nodes[n_index]; s_and_i_count+=1; search_and_insert(&node,neighbour,direction); continue; } swap_count+=1; random_s = random_tot%mult; int node_index_2 = random_s%siz; int direction_2 = random_s/siz; Node &node_2 = nodes[node_index_2]; Node *node_2_neighbour = node_2.neighbours[direction_2]; if(node_2_neighbour==NULL) { continue; } int ok_dir = node_2.is_leaf(); if(ok_dir!=-1) { if(node_index == node_index_2) { continue; } swap(node, direction, n_coords ,node_2,direction_2); swap_suc_count+=1; continue; } int ok_dir_2 = node_2_neighbour->is_leaf(); if(ok_dir_2 == -1) { continue; } if(node_2_neighbour == &node) { continue; } swap_suc_count+=1; swap(node, direction, n_coords, *node_2_neighbour, ok_dir_2); } } void swap(Node &node1, int direction1,Vertex &n_coords, Node &node2, int direction2) { //cout<<"swap\n"; Node *n2_n = node2.neighbours[direction2]; n2_n->neighbours[5-direction2] = NULL; node1.neighbours[direction1] = &node2; node2.neighbours[direction2] = NULL; node2.neighbours[5-direction1] = &node1; coords_to_node.erase(node2.v); coords_to_node[n_coords]=node2.index; node2.v = n_coords; } bool cut_and_paste(Node *X1,int ii) { //cout<<"cut and paste\n"; Node *Y1 = X1->neighbours[ii]; X1->neighbours[ii] = NULL; Y1->neighbours[5-ii] = NULL; (*Q1).reset(); (*Q2).reset(); (*Q1).put(X1); (*Q2).put(Y1); int dist = 0; increment_counter(2); unsigned long int c1 = counter; unsigned long int c2 = counter-1; X1 -> counter = c1; Y1 -> counter = c2; bool XWon = false; (*Q3).reset(); (*Q4).reset(); while(true) { dist++; Node *X = (*Q1).get(); Node *Y = (*Q2).get(); if(X==NULL) { XWon=true; break; } if(Y==NULL) { break; } (*Q3).put(X->index); (*Q4).put(Y->index); for (int i = 0; i<6; i++) { Node *neighbour = X->neighbours[i]; if(neighbour!=NULL) { if (neighbour->counter == c1) { continue; } neighbour->counter = c1; (*Q1).put(neighbour); } } for (int i = 0; i<6; i++) { Node *neighbour = Y->neighbours[i]; if(neighbour!=NULL) { if (neighbour->counter == c2) { continue; } (*Q2).put(neighbour); neighbour->counter = c2; } } } Node* small_one; circular_buffer_short* small_list = NULL; unsigned long cc; if(XWon) { small_list= &(*Q3); cc=c1; } else { small_list= &(*Q4); cc=c2; } int sss = small_list->size(); std::uniform_int_distribution<int> uni(0,sss*6-1); int r = uni(rng); int r_a = r/6; int dir = r-6*r_a; small_one = &nodes[(*small_list)[r_a]]; std::uniform_int_distribution<int> uni2(0,nodes.size()-1); while(true) { r = uni2(rng); if(nodes[r].counter!=cc) { break; } } Node *large_one =&nodes[r]; Vertex starting_coord = large_one->coord_neighbour(dir); if(coords_to_node.find(starting_coord)!=coords_to_node.end()) { X1->neighbours[ii] = Y1; Y1->neighbours[5-ii] = X1; return false; } increment_counter(1); unsigned long int c = counter; (*Q1).reset(); (*Q1).put(small_one); small_one->counter = c; small_one->temp_coords = starting_coord; std::function<Vertex (Vertex&)> transform = from_to(small_one->v,starting_coord); while(true) { Node *X = (*Q1).get(); if(X==NULL) break; for (int i = 0; i<6; i++) { Node *neighbour = X->neighbours[i]; if(neighbour!=NULL) { if (neighbour->counter == c) { continue; } neighbour->temp_coords = transform(neighbour->v); auto found = coords_to_node.find(neighbour->temp_coords); if(found!=coords_to_node.end() && nodes[found->second].counter!=cc) { X1->neighbours[ii] = Y1; Y1->neighbours[5-ii] = X1; return false; } neighbour->counter = c; (*Q1).put(neighbour); } } } for(int i=0;i<sss;i++) { Node &node = nodes[(*small_list)[i]]; coords_to_node.erase(node.v); node.v = node.temp_coords; } for(int i=0;i<sss;i++) { Node &node = nodes[(*small_list)[i]]; coords_to_node[node.v]=node.index; } large_one->neighbours[dir]=small_one; small_one->neighbours[5-dir]=large_one; return true; } pair<int,int> BFS(int start_index) { Node& start = nodes[start_index]; (*Q1).reset(); (*Q1).put(&start); vector<int> distances; distances.resize(nodes.size()); distances[start_index] = 0; increment_counter(1); unsigned long long c = counter; start.counter = c; while(true) { Node *X = (*Q1).get(); if(X==NULL) { break; } int X_distance = distances[X->index]; for (int i = 0; i<6; i++) { Node *neighbour = X->neighbours[i]; if(neighbour!=NULL) { if (neighbour->counter == c) { continue; } neighbour->counter = c; distances[neighbour->index] =X_distance+1; (*Q1).put(neighbour); } } } int ss = distances.size(); int max_dist = 0; int index = start_index; for (int i=0;i<ss;i++) { if (distances[i]>max_dist) { max_dist = distances[i]; index = i; } } return make_pair(max_dist, index); } pair<int,int> longest_path() { pair<int,int> ret = BFS(0); int start_index = ret.second; pair<int,int> ret_2 = BFS(start_index); int finish_index = ret_2.second; int distance = H_dist(nodes[start_index].v,nodes[finish_index].v); return make_pair(ret_2.first,distance); } double average_branching() { Node& start = nodes[0]; (*Q1).reset(); (*Q1).put(&start); deque<tuple<int,int>> vd_queue; stack<tuple<int,int>> vd_stack; int sss = nodes.size(); vector<int> edge_to_branch_size(6*sss); auto pos_1 = [&sss](tuple<int,int> &pair) {return sss*get<1>(pair)+get<0>(pair) ;}; auto pos_2 = [&sss](int a, int b) {return sss*b+a ;}; increment_counter(1); unsigned long long c =counter; start.counter = c; while(true) { Node *X = (*Q1).get(); if(X==NULL) { break; } for (int i = 0; i<6; i++) { Node *neighbour = X->neighbours[i]; if(neighbour!=NULL) { if (neighbour->counter == c) { continue; } neighbour->counter = c; vd_queue.push_front(make_tuple(X->index,i)); vd_queue.push_back(make_tuple(neighbour->index,5-i)); (*Q1).put(neighbour); } } } for (deque<tuple<int,int>>::reverse_iterator it = vd_queue.rbegin(); it != vd_queue.rend(); ++it ) { vd_stack.push(*it); while(true) { if (vd_stack.empty()) { break; } tuple<int,int> top_vd = vd_stack.top(); int index = pos_1(top_vd); if(edge_to_branch_size[index] > 0) { vd_stack.pop(); continue; } Node* neighbour = nodes[get<0>(top_vd)].neighbours[get<1>(top_vd)]; int n_index = neighbour->index; bool can_calculate = true; int sum = 1; for (int i=0;i<6;i++) { if(i!=5-get<1>(top_vd) && neighbour->neighbours[i]!=NULL) { tuple<int,int> to_h = make_pair(n_index,i); int present = edge_to_branch_size[pos_1(to_h)]; if(present>0) { sum+=present; } else { can_calculate = false; vd_stack.push(std::move(to_h)); } } } if(can_calculate) { edge_to_branch_size[index] = sum; vd_stack.pop(); } } } int sum = 0; for (int i=0;i<nodes.size()-1;i++) { tuple<int,int> top = std::move(vd_queue.front()); int n_index =nodes[get<0>(top)].neighbours[get<1>(top)]->index; sum+= -1+min(edge_to_branch_size[pos_1(top)],edge_to_branch_size[pos_2(n_index,5-get<1>(top))]); vd_queue.pop_front(); } double average = (double)sum/(double)(nodes.size()-1); return average; } int o_dist(Vertex &v) { long long x1 = abs(get<0>(v)); long long x2 = abs(get<1>(v)); long long x3 = abs(get<2>(v)); long long x4 = abs(get<3>(v)); long long x5 = abs(get<4>(v)); return pow(x1,.5)+x2+ pow(x3,.5)+x4+x5; } int H_dist(Vertex v1, Vertex v2) { Vertex origin = make_tuple(0,0,0,0,0); Vertex translated = from_to(v1, origin)(v2); return o_dist(translated); } void check_tree() { int count = 0; for (int i=0;i < nodes.size();i++) { if (coords_to_node[nodes[i].v]!=nodes[i].index) { cout<<"here"; } for (int j=0;j<6;j++) { Node *neighbour = nodes[i].neighbours[j]; if(neighbour!=NULL) { count++; } else { continue; } if(neighbour->neighbours[5-j]!=&nodes[i]) { cout<<"here"; } } } if(count!=2*(nodes.size()-1)) { cout<<"here"; } Node& start = nodes[0]; (*Q1).reset(); (*Q1).put(&start); increment_counter(1); unsigned long long c = counter; start.counter = c; int count_2 = 0; while(true) { Node *X = (*Q1).get(); if(X==NULL) { break; } count_2++; for (int i = 0; i<6; i++) { Node *neighbour = X->neighbours[i]; if(neighbour!=NULL) { if (neighbour->counter == c) { continue; } neighbour->counter = c; (*Q1).put(neighbour); } } } if(count_2!=nodes.size()) { cout<<"here"; } } void search_and_insert(Node *X1, Node *Y1, int ii) { //cout<<"s_ans_i\n"; (*Q1).reset(); (*Q2).reset(); (*Q1).put(X1); (*Q2).put(Y1); increment_counter(2); unsigned long int c1 = counter; unsigned long int c2 = counter-1; bool success = false; int dist = 0; Node *success_X = NULL; Node *success_Y = NULL; int X_d = 0; int Y_d = 0; X1 -> counter = c1; Y1 -> counter = c2; while(success == false) { dist++; Node *X = (*Q1).get(); Node *Y = (*Q2).get(); for (int i = 0; i<6; i++) { Node *neighbour = X->neighbours[i]; if(neighbour!=NULL) { if (neighbour->counter == c1) { continue; } if (neighbour->counter == c2) { success_X = X; success_Y = neighbour; X_d = i; Y_d = 5-i; success = true; break; } (*Q1).put(neighbour); neighbour->from = X; neighbour->counter = c1; } } if(success==true) { break; } for (int i = 0; i<6; i++) { Node *neighbour = Y->neighbours[i]; if(neighbour!=NULL) { if (neighbour->counter == c2) { continue; } if (neighbour->counter == c1) { success_Y = Y; success_X = neighbour; Y_d = i; X_d = 5-i; success = true; break; } (*Q2).put(neighbour); neighbour->from = Y; neighbour->counter = c2; } } } vector <Node *> hist = {}; hist.reserve(2*dist+2); Node *last = success_X; hist.push_back(X1); while(last != X1) { hist.push_back(last); last = last->from; } last = success_Y; while(last != Y1) { hist.push_back(last); last = last->from; } std::uniform_int_distribution<int> uni(0,hist.size()-1); int r = uni(rng); if (r == 0) { success_X->neighbours[X_d] = NULL; success_Y->neighbours[Y_d] = NULL; } else { Node *v1 = hist[r]; Node *v2 = hist[r]->from; for (int i=0;i<6;i++) { if (v1->neighbours[i] == v2) { v1->neighbours[i] = NULL; v2->neighbours[5-i] = NULL; break; } } } X1->neighbours[ii] = Y1; Y1->neighbours[5-ii] = X1; } }; // RUN LATTICE TREE SIMULATION // 1. graph_size = number of vertices in tree // 2. steps_between_samples_to_gs = (number of steps between samples)/(graph size) // 3. num_samples = number of samples // 4. samples_per_file = samples per file // 5. file_name = output path int main(int argc, char** argv) { cout<<"start\n"; int graph_size = stoi(argv[1]); if (graph_size>MS) { fprintf(stderr, "Graph size too large.\n"); return -1; } int steps_between_samples_to_gs = stoi(argv[2]); int steps_between_samples = graph_size*steps_between_samples_to_gs; int num_samples = stoi(argv[3]); int samples_per_file = stoi(argv[4]); string file_name = argv[5]; Graph G(graph_size); G.run_algorithm(num_samples,steps_between_samples,2*steps_between_samples,samples_per_file,file_name); cout<<"done\n"; }