LoU / perc_Euclidean_4d.cpp
perc_Euclidean_4d.cpp
Raw
// 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_set>
#include <functional>
#include <tuple>
#include <fstream>
#include <sched.h>
#include <stdlib.h>
#include <stdio.h>
#include <cstdlib>
#include "robin_hood.h"

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;
        }
    };

}

// LEATH TERMINATION CONDITION
// 0 TERMINATE AT NUM VERTICES
// 1 TERMINATE AT DISTANCE
#define TERM 1


using namespace std;



typedef tuple<int,int,int,int> Vertex;
vector<Vertex> vertex_neighbours(Vertex);

tuple<int,int,int,int> LeathRun(int,double);
vector<bool> rand_bernoulli(float,int);




void write_to_file(vector<tuple<int,int,int,int>> &vec, string file_name)
{
    std::ofstream f(file_name);
    for(vector<tuple<int,int,int,int>>::const_iterator i = vec.begin(); i != vec.end(); ++i) {
        f << get<0>(*i) << ','<< get<1>(*i) << ','<<get<2>(*i) <<','<<get<3>(*i)<< '\n';
    }
}


vector<Vertex> vertex_neighbours(Vertex v)
{
    int x = get<0>(v);
    int y = get<1>(v);
    int z = get<2>(v);
    int w = get<3>(v);
    vector <Vertex> ret{make_tuple(x+1,y,z,w),
              make_tuple(x-1,y,z,w),
              make_tuple(x,y+1,z,w),
              make_tuple(x,y-1,z,w),
              make_tuple(x,y,z+1,w),
              make_tuple(x,y,z-1,w),
              make_tuple(x,y,z,w+1),
              make_tuple(x,y,z,w-1)};
    return ret;
}


vector<bool> rand_bernoulli(double p,int n) {
    thread_local static random_device rd;
    thread_local static mt19937 rng(rd());
    
    thread_local static bernoulli_distribution dist(p);

    thread_local auto gen = [](){
                   return dist(rng);
           };

    vector<bool> open(n);

    generate(begin(open), end(open), gen);
    return open;
}




tuple<int,int,int,int> LeathRun(int n, double p)
{
    Vertex initial_vertex =  std::make_tuple(0,0,0,0);
    robin_hood::unordered_set <Vertex,hash_tuple::hash<Vertex>> visited_vertices;
    queue<Vertex> vertex_frontier;
    vertex_frontier.push(initial_vertex);
    
    int visited = 1;

    int num_open_edges = 0;
    int num_closed_edges = 0;
    
    int max_dist =  0;
#if TERM == 0
    while(visited<n)
#endif
#if TERM == 1
    while(max_dist<n)
#endif
    {
        if (vertex_frontier.empty())
        {
            break;
        }


            Vertex current = vertex_frontier.front();
        vertex_frontier.pop();
        
        max_dist = max(max_dist,abs(get<0>(current))+abs(get<1>(current))+abs(get<2>(current))+abs(get<3>(current)));
        
            if (visited_vertices.find(current) != visited_vertices.end())
        {
                continue;
        }
        
        
        visited_vertices.insert(current);
        
        vector<Vertex> neighbours = vertex_neighbours(current);
            

        vector<bool> open = rand_bernoulli(p,8);
        int num_open = 0;
        int num = 0;
        for(int i=0; i < 8; i++)
        {
             if(visited_vertices.find(neighbours[i]) == visited_vertices.end())
            {
                if(open[num] == true)
                {
                    num_open += 1;
                    vertex_frontier.push(neighbours[i]);
                }
                num+=1;
            }

        }

        num_open_edges += num_open;
        num_closed_edges += (num-num_open);
        visited+=1;

    }
    return make_tuple(visited,num_open_edges,num_closed_edges,max_dist);
}


// RUN LEATH ALGORITHM
// 1. n = max distance/number of steps (depending on TERM)
// 2. p = percolation probability
// 3. N = number of samples
// 4. filename = output path
int main(int argc, char** argv)
{
    int n = stoi(argv[1]);
    double p = stod(argv[2]);
    int N = stoi(argv[3]);
    string filename = argv[4];
    std::cout << "n="<<n<<", p="<<p<<", N= "<<N<<"\n\n";
    
    int num_per_file = 4000;
    int num_files_per_directory = 50;
    

    int total_per_directory = num_files_per_directory * num_per_file;

    int num_directories = int(ceil(float(N)/float(total_per_directory)));
    
    int num_files_last_directory = int(ceil( float(N-total_per_directory * (num_directories-1))/float(num_per_file)));

    int num_in_last_file = N-total_per_directory * (num_directories-1)-(num_files_last_directory-1)*num_per_file;


    for (int i=0; i<num_directories; i++)
    {
        string directory_name = filename+"/"+to_string(i);
        system(("mkdir "+ directory_name).c_str());
        int num_files_this_dir;
        if(i == num_directories-1)
        {
            num_files_this_dir = num_files_last_directory;
        }
        else
        {
            num_files_this_dir = num_files_per_directory;
        }
        
        for (int j=0; j<num_files_this_dir; j++)
        {
            
            string file_name = directory_name+"/res_"+to_string(j);
            int num_this_file;
            if(i == num_directories-1 && j == num_files_this_dir-1)
            {
                num_this_file = num_in_last_file;
            }
            else
            {
                num_this_file = num_per_file;

            }

            vector<tuple<int,int,int,int>> file_results(num_this_file);

            for (int w=0; w<num_this_file; w++)
            {
                file_results[w] = LeathRun(n,p);
            }
            write_to_file(file_results, file_name);
        }
    }
    

}