// 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"

typedef unsigned long long int LI_t;

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

}

// 0 INVASION PERCOLATION
// 1 LEATH ALGORITHM
#define RUN 0


// 1 H1
// 2 H2
#define TYPE 1


using namespace std;



typedef tuple<LI_t,int> Vertex;
typedef tuple<double,LI_t,int> VVertex;

vector<Vertex> vertex_neighbours(Vertex);

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


void write_to_file(vector<double> vec, string file_name)
{
    std::ofstream outFile(file_name);
    for (const auto &e : vec) outFile << e << "\n";
}

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

LI_t rand_long() {
    thread_local static random_device rd;
    thread_local static mt19937_64 rng(rd());
    
    thread_local static uniform_int_distribution<LI_t> dist;

    return dist(rng);
}
vector<double> rand_uniform(int n) {
    thread_local static random_device rd;
    thread_local static mt19937_64 rng(rd());
    
    thread_local static uniform_real_distribution<double> dist(0.0,1.0);

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

    vector<double> values(n);

    generate(begin(values), end(values), gen);
    return values;
}
//The fractal coordinates are encoded in binary
vector<LI_t> neighbours(LI_t vertex)
{

    short first = vertex & 3;

    LI_t cleared = vertex;
    cleared &= ~1ULL;
    cleared &= ~(1ULL<<1);
    
    
    
    
    
    LI_t vertex_copy = vertex;
    LI_t flipped = vertex;
    
    
    
    
    flipped &= ~1ULL;
    flipped &= ~(1ULL<<1);

    
    int i=0;
    for(i=0;i<31;i++)
    {
        vertex_copy = vertex_copy>>2;
        flipped &= ~(1ULL << (2*i+2));
        flipped &= ~(1ULL << (2*i+3));
        if((vertex_copy &3) != first)
        {
            break;
        }
    }
   if(i==30)
   {
    abort();
   }
         
    short after = vertex_copy & 3;
    
    short replace_first;
    short replace_after;
    
    if(first == 1 && after == 0)
    {
        replace_first = 0;
        replace_after = 1;
    }
    else if(first == 0 && after==1)
    {
        replace_first = 1;
        replace_after = 0;
    }
    else if(first == 1 && after==2)
    {
        replace_first = 2;
        replace_after = 1;
    }
    else if(first == 2 && after==1)
    {
        replace_first = 1;
        replace_after = 2;
    }
    else if(first == 2 && after==3)
    {
        replace_first = 3;
        replace_after = 2;
    }
    else if(first == 3 && after==2)
    {
        replace_first = 2;
        replace_after = 3;
    }
    else
    {
        vector<LI_t> nbs;
        if(first == 0)
        {
            nbs = {cleared|1};
        }
        else if(first == 1)
        {
            nbs = {cleared|0,cleared|2};
        }
        else if(first == 2)
        {
            nbs = {cleared|1,cleared|3};
        }
        else
        {
            nbs = {cleared|2};
        }
        return nbs;
    }

    for (int j = 0;j<=i;j++)
    {
        flipped |=(replace_first<<2*j);
    }
    flipped |=(replace_after<<2*(i+1));
    
    vector<LI_t> nbs;
    if(first == 0)
    {
        nbs = {cleared|1,flipped};
    }
    else if(first == 1)
    {
        nbs = {cleared|0,cleared|2,flipped};
    }
    else if(first == 2)
    {
        nbs = {cleared|1,cleared|3,flipped};
    }
    else
    {
        nbs = {cleared|2,flipped};
    }
    
    
    return nbs;
    
}

vector<LI_t> neighbours2(LI_t vertex)
{
    short first = vertex & 3;
    LI_t cleared = vertex;
    cleared &= ~1ULL;
    cleared &= ~(1ULL<<1);
    
    vector<LI_t> nbs;
    if(first == 0)
    {
        nbs = {cleared|1};
    }
    else if(first == 1)
    {
        nbs = {cleared|0,cleared|2};
    }
    else if(first == 2)
    {
        nbs = {cleared|1,cleared|3};
    }
    else
    {
        nbs = {cleared|2};
    }
    
    
    
    LI_t vertex_copy = vertex;
    LI_t flipped = vertex;
    
    
    
    
    flipped &= ~1ULL;
    flipped &= ~(1ULL<<1);
    //std::cout<<bitset<8>(flipped)<<"\n";

    
    int i=0;
    for(i=0;i<31;i++)
    {
        vertex_copy = vertex_copy>>2;
        flipped &= ~(1ULL << (2*i+2));
        flipped &= ~(1ULL << (2*i+3));
        if((vertex_copy &3) != first)
        {
            break;
        }
    }
    if(i==30)
    {
        cout<<"note";
        abort();
    }
         
    short after = vertex_copy & 3;
    
    short replace_first;
    short replace_after;
    
    short replace_first_2=-1;
    short replace_after_2=-1;
    
    if(first == 2 && after == 0)
    {
        replace_first = 3;
        replace_after = 1;
    }
    else if(first == 3 && after==1)
    {
        replace_first = 2;
        replace_after = 0;
        
        replace_first_2=0;
        replace_after_2=2;
    }
    else if(first == 0 && after==2)
    {
        replace_first = 3;
        replace_after = 1;
        
        replace_first_2 = 1;
        replace_after_2 = 3;
    }
    else if(first == 1 && after==3)
    {
        replace_first = 0;
        replace_after = 2;
    }
    else
    {
        return nbs;
    }
    
    LI_t flipped_2 = flipped;
    
    for (int j = 0;j<=i;j++)
    {
        flipped |=(replace_first<<2*j);
    }
    flipped |=(replace_after<<2*(i+1));
    
    if(replace_first_2!=-1)
    {
        for (int j = 0;j<=i;j++)
        {
            flipped_2 |=(replace_first_2<<2*j);
        }
        flipped_2 |=(replace_after_2<<2*(i+1));
    }
    
    nbs.push_back(flipped);
    if(replace_first_2!=-1)
    {
        nbs.push_back(flipped_2);
    }
    return nbs;
    
}
vector<Vertex> vertex_neighbours(Vertex v)
{
    LI_t HH = get<0>(v);
    int EE = get<1>(v);
    
    //H1
    #if TYPE == 1
    vector<LI_t> Hneighbours = neighbours(HH);
    #endif
    
    //H2
    #if TYPE == 2
    vector<LI_t> Hneighbours = neighbours2(HH);
    #endif
    short s = Hneighbours.size();
    
    vector <Vertex> ret;
    ret.reserve(2+s);
    ret.emplace_back(HH,EE-1);
    ret.emplace_back(HH,EE+1);
    
    for(int i=0;i<s;i++)
    {
        ret.emplace_back(Hneighbours[i],EE);
    }
    
    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> LeathRun(int n, double p)
{
    Vertex initial_vertex =  std::make_tuple(rand_long(),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;
    
    while(visited < n)
    {


        if (vertex_frontier.empty())
        {
            break;
        }


            Vertex current = vertex_frontier.front();
        vertex_frontier.pop();
                
            if (visited_vertices.find(current) != visited_vertices.end())
        {
                continue;
        }
        
        
        visited_vertices.insert(current);
        
        vector<Vertex> neighbours = vertex_neighbours(current);
        short nn = neighbours.size();

        vector<bool> open = rand_bernoulli(p,nn);
        int num_open = 0;
        int num = 0;
        for(int i=0; i < nn; 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);
}

class Compare
{
public:
    bool operator() (VVertex& v1, VVertex& v2)
    {
        return get<0>(v1)>get<0>(v2);
    }
};

vector<int> invasion_percolation(int n)
{
    int max_t = ceil(4*log(float(n)/100)/log(2))+1;
    //std::cout<<max_t<<"\n";
    vector<int> rec_vals(max_t);
    for (int i=0;i<max_t;i++)
    {
        rec_vals[i] = floor(100*pow(2,float(i)/4));
    }
    int rec_index=0;
    Vertex initial_vertex =  std::make_tuple(rand_long(),0);
    robin_hood::unordered_set <Vertex,hash_tuple::hash<Vertex>> visited_vertices;

    visited_vertices.reserve(n+1);
    visited_vertices.insert(initial_vertex);

    vector<Vertex> frontier_edges = vertex_neighbours(initial_vertex);
    int s = frontier_edges.size();
    vector<double> frontier_values = rand_uniform(s);
    //vector<float> frontier_values{0.13944587, 0.90309675, 0.62820165, 0.47304682};

    priority_queue<VVertex, std::deque <VVertex>, Compare> frontier;
    //std::set<VVertex, Compare> frontier;


    for (int i=0;i<s;i++)
    {
        frontier.push(tuple_cat(make_tuple(frontier_values[i]),frontier_edges[i]));
    }
    
    vector<int> visited_numbers;
    visited_numbers.reserve(max_t);
    int visited_number_running = 4;

    vector<double> random = rand_uniform(6+5*n);

    double current_marker = 1;

    int visited = 0;
    
    while(visited<n)
    {
        int num = frontier.size();
        if (num==0)
        {
            abort();
            break;
        }

        
        /*if(visited>0 && visited%1000000==0)
        {
            auto it=frontier.lower_bound(make_tuple(current_marker,0,0));
            frontier.erase(it,frontier.end());
        }*/
        
        
        
        visited+=1;

        VVertex currentE = frontier.top();
            
        Vertex current = make_tuple(get<1>(currentE),get<2>(currentE));
        frontier.pop();

        //current_marker=1;
        
        //H2 and H1
        current_marker = (double)visited/(double)visited_number_running+2.2221*pow(visited,-0.382);

        

        if (visited_vertices.find(current) != visited_vertices.end())
        {
            if (visited==rec_vals[rec_index])
            {
                rec_index+=1;
                visited_numbers.push_back(visited_number_running);
            }
            continue;
        }
        
        visited_vertices.insert(current);
        
        vector<Vertex> neighbours = vertex_neighbours(current);

        int init = 6+(visited-1)*5;

        int num_n_v = 0;
        short ss = neighbours.size();
        for (int i=0;i<ss;i++)
        {
            Vertex &nb = neighbours[i];
            if (visited_vertices.find(nb) == visited_vertices.end())
            {
                
                double &value = random[init+num_n_v];
                num_n_v += 1;

                if (value<current_marker)
                {
                    frontier.emplace(value,get<0>(nb),get<1>(nb));
                }
                
            }
        }


        visited_number_running += num_n_v;
        if(visited==rec_vals[rec_index])
        {
            rec_index+=1;
            visited_numbers.push_back(visited_number_running);
        }
    }
    return visited_numbers;
    
}

// 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
#if RUN == 1
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 = 200;
    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>> 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);
        }
    }
}
#endif


// RUN INVASION PERCOLATION
// Inputs:
// 1. n = length of each run
// 2. N = total number of samples
// 3. num_runs_per_file = number of runs stored in each file
// 4. filename = output path
#if RUN == 0
 int main(int argc, char** argv)
 {
     int n = stoi(argv[1]);
     int N = stoi(argv[2]);
     int num_runs_per_file = stoi(argv[3]);
     string filename = argv[4];

     
     int num_files = ceil(float(N)/float(num_runs_per_file));
     int runs_last_file = N-(num_files-1) * num_runs_per_file;
     for (int j=0;j<num_files;j++)
     {
         int num_runs_this_file;
         if(j==num_files-1)
         {
             num_runs_this_file = runs_last_file;
         }
         else
         {
             num_runs_this_file = num_runs_per_file;
         }
         
         vector<int> results;
         results.reserve((ceil(4*log(float(n)/100)/log(2))+1)* num_runs_this_file);

         for (int i=0;i<num_runs_this_file;i++)
         {
             vector<int> visited_numbers = invasion_percolation(n);
             copy (visited_numbers.begin(), visited_numbers.end(),back_inserter(results));
         }
         write_to_file(results,filename+"/f"+to_string(j));

     }
         std::cout << "done\n";
 }
#endif