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

}

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

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

using namespace std;

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


vector<Vertex> vertex_neighbours(Vertex);
tuple<int,int,int,int> LeathRun(int,double);
vector<bool> rand_bernoulli(float,int);
void print_vertex(Vertex);



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

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<int> vec, string file_name)
{
    std::ofstream outFile(file_name);
    for (const auto &e : vec) outFile << e << "\n";
}


vector<Vertex> vertex_neighbours(Vertex v)
{
    int x = get<0>(v);
    int y = get<1>(v);
    int z = get<2>(v);
    vector <Vertex> ret{make_tuple(x+1,y,z),
              make_tuple(x-1,y,z),
              make_tuple(x,y+1,z+x),
              make_tuple(x,y-1,z-x)};
    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;
}

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



void print_vertex(Vertex v)
{
    std::cout<< "("<<get<0>(v)<<","<<get<1>(v)<<","<<get<2>(v)<<")\n";
}

int dist(int x,int y,int z)
{
    if (x==0 && y==0 && z==0)
    {
        return 0;
    }
    if(z<=0)
    {
        y=-y;
        z=-z;
    }
    if (abs(y)>abs(x))
    {
        int temp = y;
        y=x;
        x=temp;
    }
    if (x<=0)
    {
        x=-x;
        y=-y;
    }

    if(y>=0)
    {
        float sz  = float(sqrt(z));
        if(x<=sz)
        {
            return 2*(ceil(2*sz))-x-y;
        }
        else
        {
            if(x*y>=z)
            {
                return x+y;
            }
            else
            {
                return 2*ceil(float(z)/float(x))+x-y;
            }

        }

    }
    else
    {
        float sz = sqrt(float(z-x*y));
        if (x<=sz)
        {
            return 2*ceil(2*sz)-x+y;
        }
        else
        {
            return 2*ceil(float(z)/float(x))+x-y;
        }

    }

}

tuple<int,int,int,int> LeathRun(int n, double p)
{
    Vertex initial_vertex =  std::make_tuple(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;
    // OR max_dist < n
#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,dist(get<0>(current),get<1>(current),get<2>(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,4);
        int num_open = 0;
        int num = 0;
        for(int i=0; i < 4; 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);
}

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;
    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(0,0,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);

    priority_queue<VVertex, std::deque <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;
    
    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),get<3>(currentE));

        
        
        
        frontier.pop();
        
        //current_marker=1;
        current_marker =  (double)visited/(double)visited_number_running+2.347*pow(visited,-0.4);

        
        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 num_n_v = 0;
        short ss = neighbours.size();
        vector<double> random = rand_uniform(ss);
        
        for (int i=0;i<ss;i++)
        {
            Vertex &nb = neighbours[i];
            if (visited_vertices.find(nb) == visited_vertices.end())
            {

                double &value = random[num_n_v];
                num_n_v += 1;

                if (value<=current_marker)
                {
                frontier.emplace(value,get<0>(nb),get<1>(nb),get<2>(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 = 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);
        }
    }
    

}
#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 = path of file in which to store output
#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