// 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;
typedef int pLI_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;
        }
    };

}



using namespace std;

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

//3 H_3
//4 H_4
#define TYPE 4

#if TYPE == 3 //H3
typedef tuple<LI_t,LI_t,pLI_t,pLI_t> Vertex;
typedef tuple<double,LI_t,LI_t,pLI_t,pLI_t> VVertex;
#endif

#if TYPE == 4 //H4
typedef tuple<LI_t,LI_t,LI_t,pLI_t> Vertex;
typedef tuple<double,LI_t,LI_t,LI_t,pLI_t> VVertex;
#endif

vector<Vertex> vertex_neighbours(Vertex);

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



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

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);
}
//The fractal coordinates are encoded in binary
vector<LI_t> nbs4(LI_t vertex)
{
    short first = vertex &15;
    
    LI_t cleared = vertex;
    LI_t vertex_copy = vertex;
    for(int i = 0;i <4;i++)
    {
        cleared &= ~(1ULL<<i);
    }
    LI_t flipped = cleared;
    int i=0;
    for(i = 0;i <16;i++)
    {
        vertex_copy = vertex_copy >>4;
        flipped &= ~(1ULL << (4*i+4));
        flipped &= ~(1ULL << (4*i+5));
        flipped &= ~(1ULL << (4*i+6));
        flipped &= ~(1ULL << (4*i+7));
        if((vertex_copy &15)!=first)
        {
            break;
        }
    }
    short after = vertex_copy &15;
    short replace_first_1=-1;
    short replace_after_1;
    short replace_first_2;
    short replace_after_2;
    if(first ==0 && after ==15)
    {
        replace_first_1=1;
        replace_after_1=0;
        replace_first_2=15;
        replace_after_2=14;
    }
    else if(first ==15 && after ==14)
    {
        replace_first_1=0;
        replace_after_1=15;
        replace_first_2=14;
        replace_after_2=13;
    }
    else if(first ==14 && after ==13)
    {
        replace_first_1=15;
        replace_after_1=14;
        replace_first_2=13;
        replace_after_2=12;
    }
    else if(first ==13 && after ==12)
    {
        replace_first_1=14;
        replace_after_1=13;
        replace_first_2=12;
        replace_after_2=11;
    }
    else if(first ==12 && after ==11)
    {
        replace_first_1=13;
        replace_after_1=12;
        replace_first_2=11;
        replace_after_2=10;
    }
    else if(first ==11 && after ==10)
    {
        replace_first_1=12;
        replace_after_1=11;
        replace_first_2=10;
        replace_after_2=9;
    }
    else if(first ==10 && after ==9)
    {
        replace_first_1=11;
        replace_after_1=10;
        replace_first_2=9;
        replace_after_2=8;
    }
    else if(first ==9 && after ==8)
    {
        replace_first_1=10;
        replace_after_1=9;
        replace_first_2=8;
        replace_after_2=7;
    }
    else if(first ==8 && after ==7)
    {
        replace_first_1=9;
        replace_after_1=8;
        replace_first_2=7;
        replace_after_2=6;
    }
    else if(first ==7 && after ==6)
    {
        replace_first_1=8;
        replace_after_1=7;
        replace_first_2=6;
        replace_after_2=5;
    }
    else if(first ==6 && after ==5)
    {
        replace_first_1=7;
        replace_after_1=6;
        replace_first_2=5;
        replace_after_2=4;
    }
    else if(first ==5 && after ==4)
    {
        replace_first_1=6;
        replace_after_1=5;
        replace_first_2=4;
        replace_after_2=3;
    }
    else if(first ==4 && after ==3)
    {
        replace_first_1=5;
        replace_after_1=4;
        replace_first_2=3;
        replace_after_2=2;
    }
    else if(first ==3 && after ==2)
    {
        replace_first_1=4;
        replace_after_1=3;
        replace_first_2=2;
        replace_after_2=1;
    }
    else if(first ==2 && after ==1)
    {
        replace_first_1=3;
        replace_after_1=2;
        replace_first_2=1;
        replace_after_2=0;
    }
    else if(first ==1 && after ==0)
    {
        replace_first_1=2;
        replace_after_1=1;
        replace_first_2=0;
        replace_after_2=15;
    }
    else
    {
        vector<LI_t> nbs;
        if(first == 0 )
        {
            nbs={cleared|1,cleared|15};
        }
        else if(first == 1)
        {
            nbs={cleared|2,cleared|0};
        }
        else if(first == 2)
        {
            nbs={cleared|3,cleared|1};
        }
        else if(first == 3)
        {
            nbs={cleared|4,cleared|2};
        }
        else if(first == 4)
        {
            nbs={cleared|5,cleared|3};
        }
        else if(first == 5)
        {
            nbs={cleared|6,cleared|4};
        }
        else if(first == 6)
        {
            nbs={cleared|7,cleared|5};
        }
        else if(first == 7)
        {
            nbs={cleared|8,cleared|6};
        }
        else if(first == 8)
        {
            nbs={cleared|9,cleared|7};
        }
        else if(first == 9)
        {
            nbs={cleared|10,cleared|8};
        }
        else if(first == 10)
        {
            nbs={cleared|11,cleared|9};
        }
        else if(first == 11)
        {
            nbs={cleared|12,cleared|10};
        }
        else if(first == 12)
        {
            nbs={cleared|13,cleared|11};
        }
        else if(first == 13)
        {
            nbs={cleared|14,cleared|12};
        }
        else if(first == 14)
        {
            nbs={cleared|15,cleared|13};
        }
        else if(first == 15)
        {
            nbs={cleared|0,cleared|14};
        }
        return nbs;
    }
    LI_t flipped_2 = flipped;
    for(int j=0;j<=i;j++)
    {
        flipped |=(replace_first_1<<4*j);
    }
    flipped |=(replace_after_1<<4*(i+1));
    for(int j=0;j<=i;j++)
    {
        flipped_2 |=(replace_first_2<<4*j);
    }
    flipped_2 |=(replace_after_2<<4*(i+1));
    vector<LI_t> nbs;
    if(first == 0 )
    {
        nbs={flipped,flipped_2,cleared|1,cleared|15};
    }
    else if(first == 1)
    {
        nbs={flipped,flipped_2,cleared|2,cleared|0};
    }
    else if(first == 2)
    {
        nbs={flipped,flipped_2,cleared|3,cleared|1};
    }
    else if(first == 3)
    {
        nbs={flipped,flipped_2,cleared|4,cleared|2};
    }
    else if(first == 4)
    {
        nbs={flipped,flipped_2,cleared|5,cleared|3};
    }
    else if(first == 5)
    {
        nbs={flipped,flipped_2,cleared|6,cleared|4};
    }
    else if(first == 6)
    {
        nbs={flipped,flipped_2,cleared|7,cleared|5};
    }
    else if(first == 7)
    {
        nbs={flipped,flipped_2,cleared|8,cleared|6};
    }
    else if(first == 8)
    {
        nbs={flipped,flipped_2,cleared|9,cleared|7};
    }
    else if(first == 9)
    {
        nbs={flipped,flipped_2,cleared|10,cleared|8};
    }
    else if(first == 10)
    {
        nbs={flipped,flipped_2,cleared|11,cleared|9};
    }
    else if(first == 11)
    {
        nbs={flipped,flipped_2,cleared|12,cleared|10};
    }
    else if(first == 12)
    {
        nbs={flipped,flipped_2,cleared|13,cleared|11};
    }
    else if(first == 13)
    {
        nbs={flipped,flipped_2,cleared|14,cleared|12};
    }
    else if(first == 14)
    {
        nbs={flipped,flipped_2,cleared|15,cleared|13};
    }
    else if(first == 15)
    {
        nbs={flipped,flipped_2,cleared|0,cleared|14};
    }
    return nbs;
}

vector<LI_t> nbs3(LI_t vertex)
{
    short first = vertex &7;
    
    LI_t cleared = vertex;
    LI_t vertex_copy = vertex;
    for(int i = 0;i <3;i++)
    {
        cleared &= ~(1ULL<<i);
    }
    LI_t flipped = cleared;
    int i=0;
    for(i = 0;i <21;i++)
    {
        vertex_copy = vertex_copy >>3;
        flipped &= ~(1ULL << (3*i+3));
        flipped &= ~(1ULL << (3*i+4));
        flipped &= ~(1ULL << (3*i+5));
        if((vertex_copy &7)!=first)
        {
            break;
        }
    }
    short after = vertex_copy &7;
    short replace_first_1=-1;
    short replace_after_1;
    short replace_first_2;
    short replace_after_2;
    if(first ==0 && after ==7)
    {
        replace_first_1=1;
        replace_after_1=0;
        replace_first_2=7;
        replace_after_2=6;
    }
    else if(first ==7 && after ==6)
    {
        replace_first_1=0;
        replace_after_1=7;
        replace_first_2=6;
        replace_after_2=5;
    }
    else if(first ==6 && after ==5)
    {
        replace_first_1=7;
        replace_after_1=6;
        replace_first_2=5;
        replace_after_2=4;
    }
    else if(first ==5 && after ==4)
    {
        replace_first_1=6;
        replace_after_1=5;
        replace_first_2=4;
        replace_after_2=3;
    }
    else if(first ==4 && after ==3)
    {
        replace_first_1=5;
        replace_after_1=4;
        replace_first_2=3;
        replace_after_2=2;
    }
    else if(first ==3 && after ==2)
    {
        replace_first_1=4;
        replace_after_1=3;
        replace_first_2=2;
        replace_after_2=1;
    }
    else if(first ==2 && after ==1)
    {
        replace_first_1=3;
        replace_after_1=2;
        replace_first_2=1;
        replace_after_2=0;
    }
    else if(first ==1 && after ==0)
    {
        replace_first_1=2;
        replace_after_1=1;
        replace_first_2=0;
        replace_after_2=7;
    }
    else
    {
        vector<LI_t> nbs;
        if(first == 0 )
        {
            nbs={cleared|1,cleared|7};
        }
        else if(first == 1)
        {
            nbs={cleared|2,cleared|0};
        }
        else if(first == 2)
        {
            nbs={cleared|3,cleared|1};
        }
        else if(first == 3)
        {
            nbs={cleared|4,cleared|2};
        }
        else if(first == 4)
        {
            nbs={cleared|5,cleared|3};
        }
        else if(first == 5)
        {
            nbs={cleared|6,cleared|4};
        }
        else if(first == 6)
        {
            nbs={cleared|7,cleared|5};
        }
        else if(first == 7)
        {
            nbs={cleared|0,cleared|6};
        }
        return nbs;
    }
    LI_t flipped_2 = flipped;
    for(int j=0;j<=i;j++)
    {
        flipped |=(replace_first_1<<3*j);
    }
    flipped |=(replace_after_1<<3*(i+1));
    for(int j=0;j<=i;j++)
    {
        flipped_2 |=(replace_first_2<<3*j);
    }
    flipped_2 |=(replace_after_2<<3*(i+1));
    vector<LI_t> nbs;
    if(first == 0 )
    {
        nbs={flipped,flipped_2,cleared|1,cleared|7};
    }
    else if(first == 1)
    {
        nbs={flipped,flipped_2,cleared|2,cleared|0};
    }
    else if(first == 2)
    {
        nbs={flipped,flipped_2,cleared|3,cleared|1};
    }
    else if(first == 3)
    {
        nbs={flipped,flipped_2,cleared|4,cleared|2};
    }
    else if(first == 4)
    {
        nbs={flipped,flipped_2,cleared|5,cleared|3};
    }
    else if(first == 5)
    {
        nbs={flipped,flipped_2,cleared|6,cleared|4};
    }
    else if(first == 6)
    {
        nbs={flipped,flipped_2,cleared|7,cleared|5};
    }
    else if(first == 7)
    {
        nbs={flipped,flipped_2,cleared|0,cleared|6};
    }
    return nbs;
}
vector<LI_t> nbs3o2(LI_t vertex)
{
    short first = vertex &7;
    
    LI_t cleared = vertex;
    LI_t vertex_copy = vertex;
    for(int i = 0;i <3;i++)
    {
        cleared &= ~(1ULL<<i);
    }
    LI_t flipped = cleared;
    int i=0;
    for(i = 0;i <21;i++)
    {
        vertex_copy = vertex_copy >>3;
        flipped &= ~(1ULL << (3*i+3));
        flipped &= ~(1ULL << (3*i+4));
        flipped &= ~(1ULL << (3*i+5));
        if((vertex_copy &7)!=first)
        {
            break;
        }
    }
    short after = vertex_copy &7;
    short replace_first_1;
    short replace_after_1;
    short replace_first_2=-1;
    short replace_after_2;
    if(after==1 && first ==5)
    {
        replace_after_1=0;
        replace_first_1=1;
    }
    else if(after==0 && first ==1)
    {
        replace_after_1=1;
        replace_first_1=5;
    }
    else if(after==7 && first ==3)
    {
        replace_after_1=6;
        replace_first_1=7;
    }
    else if(after==6 && first ==7)
    {
        replace_after_1=7;
        replace_first_1=3;
    }
    else if(after==2 && first ==3)
    {
        replace_after_1=3;
        replace_first_1=7;
    }
    else if(after==3 && first ==7)
    {
        replace_after_1=2;
        replace_first_1=3;
    }
    else if(after==4 && first ==5)
    {
        replace_after_1=5;
        replace_first_1=1;
    }
    else if(after==5 && first ==1)
    {
        replace_after_1=4;
        replace_first_1=5;
    }
    else if(after==0 && first ==5)
    {
        replace_after_1=2;
        replace_first_1=7;
        
        replace_after_2=6;
        replace_first_2=3;
    }
    else if(after==2 && first ==7)
    {
        replace_after_1=0;
        replace_first_1=5;
        
        replace_after_2=4;
        replace_first_2=1;
    }
    else if(after==4 && first ==1)
    {
        replace_after_1=2;
        replace_first_1=7;
        
        replace_after_2=6;
        replace_first_2=3;
    }
    else if(after==6 && first ==3)
    {
        replace_after_1=0;
        replace_first_1=5;
        
        replace_after_2=4;
        replace_first_2=1;
    }
    else
    {
        vector<LI_t> nbs;
        if(first == 0 )
        {
            nbs={cleared|1,cleared|2,cleared|6};
        }
        else if(first == 1)
        {
            nbs={cleared|0};
        }
        else if(first == 2)
        {
            nbs={cleared|3,cleared|4,cleared|0};
        }
        else if(first == 3)
        {
            nbs={cleared|2};
        }
        else if(first == 4)
        {
            nbs={cleared|5,cleared|6,cleared|2};
        }
        else if(first == 5)
        {
            nbs={cleared|4};
        }
        else if(first == 6)
        {
            nbs={cleared|7,cleared|0,cleared|4};
        }
        else if(first == 7)
        {
            nbs={cleared|6};
        }
        return nbs;
    }
    if(replace_first_2!=-1)
    {
        LI_t flipped_2 = flipped;
        for(int j=0;j<=i;j++)
        {
            flipped |=(replace_first_1<<3*j);
        }
        flipped |=(replace_after_1<<3*(i+1));
        for(int j=0;j<=i;j++)
        {
            flipped_2 |=(replace_first_2<<3*j);
        }
        flipped_2 |=(replace_after_2<<3*(i+1));
        vector<LI_t> nbs;
        if(first == 0 )
        {
            nbs={flipped,flipped_2,cleared|1,cleared|2,cleared|6};
        }
        else if(first == 1)
        {
            nbs={flipped,flipped_2,cleared|0};
        }
        else if(first == 2)
        {
            nbs={flipped,flipped_2,cleared|3,cleared|4,cleared|0};
        }
        else if(first == 3)
        {
            nbs={flipped,flipped_2,cleared|2};
        }
        else if(first == 4)
        {
            nbs={flipped,flipped_2,cleared|5,cleared|6,cleared|2};
        }
        else if(first == 5)
        {
            nbs={flipped,flipped_2,cleared|4};
        }
        else if(first == 6)
        {
            nbs={flipped,flipped_2,cleared|7,cleared|0,cleared|4};
        }
        else if(first == 7)
        {
            nbs={flipped,flipped_2,cleared|6};
        }
        return nbs;
    }
    else {
        for(int j=0;j<=i;j++)
        {
            flipped |=(replace_first_1<<3*j);
        }
        flipped |=(replace_after_1<<3*(i+1));
        vector<LI_t> nbs;
        if(first == 0 )
        {
            nbs={flipped,cleared|1,cleared|2,cleared|6};
        }
        else if(first == 1)
        {
            nbs={flipped,cleared|0};
        }
        else if(first == 2)
        {
            nbs={flipped,cleared|3,cleared|4,cleared|0};
        }
        else if(first == 3)
        {
            nbs={flipped,cleared|2};
        }
        else if(first == 4)
        {
            nbs={flipped,cleared|5,cleared|6,cleared|2};
        }
        else if(first == 5)
        {
            nbs={flipped,cleared|4};
        }
        else if(first == 6)
        {
            nbs={flipped,cleared|7,cleared|0,cleared|4};
        }
        else if(first == 7)
        {
            nbs={flipped,cleared|6};
        }
        return nbs;
    }
}







#if TYPE ==3
vector<Vertex> vertex_neighbours(Vertex v)
{
    LI_t nbs3_1 = get<0>(v);
    LI_t nbs3_2 = get<1>(v);
    pLI_t z1 = get<2>(v);
    pLI_t z2 = get<3>(v);
    
    vector<LI_t> nbs3_1_neighbours = nbs3(nbs3_1);
    vector<LI_t> nbs3_2_neighbours = nbs3(nbs3_2);
    //vector<LI_t> quad_neighbours = neighbours_cub(quad);
    
    //std::cout<<tri_1_neighbours.size()<<","<<tri_2_neighbours.size()<<","<<quad_neighbours.size()<<"\n";
    
    vector<Vertex> neighbours;
    //neighbours.reserve(tri_1_neighbours.size()+tri_2_neighbours.size()+quad_neighbours.size());
    neighbours.reserve(4+nbs3_1_neighbours.size()+nbs3_2_neighbours.size());
    

    for (int i=0;i<nbs3_1_neighbours.size();i++)
    {
        neighbours.emplace_back(nbs3_1_neighbours[i],nbs3_2,z1,z2);
    }
    for (int i=0;i<nbs3_2_neighbours.size();i++)
    {
        neighbours.emplace_back(nbs3_1,nbs3_2_neighbours[i],z1,z2);
    }
    neighbours.emplace_back(nbs3_1,nbs3_2,z1-1,z2);
    neighbours.emplace_back(nbs3_1,nbs3_2,z1+1,z2);
    neighbours.emplace_back(nbs3_1,nbs3_2,z1,z2-1);
    neighbours.emplace_back(nbs3_1,nbs3_2,z1,z2+1);
   


    return neighbours;
}
#endif
#if TYPE == 4
vector<Vertex> vertex_neighbours(Vertex v)
{
    LI_t a3o2_1 = get<0>(v);
    LI_t a3o2_2 = get<1>(v);
    LI_t a4_1 = get<2>(v);
    pLI_t z = get<3>(v);
    
    vector<LI_t> a3o2_1_neighbours = nbs3o2(a3o2_1);
    vector<LI_t> a3o2_2_neighbours = nbs3o2(a3o2_2);
    vector<LI_t> a4_1_neighbours = nbs4(a4_1);
    
    
    vector<Vertex> neighbours;
    neighbours.reserve(a3o2_1_neighbours.size()+\
                       a3o2_2_neighbours.size()+\
                       a4_1_neighbours.size()+2);
    
    
    
    for (int i=0;i<a3o2_1_neighbours.size();i++)
    {
        neighbours.emplace_back(a3o2_1_neighbours[i],a3o2_2,a4_1,z);
    }
    for (int i=0;i<a3o2_2_neighbours.size();i++)
    {
        neighbours.emplace_back(a3o2_1,a3o2_2_neighbours[i],a4_1,z);
    }
    for (int i=0;i<a4_1_neighbours.size();i++)
    {
        neighbours.emplace_back(a3o2_1,a3o2_2,a4_1_neighbours[i],z);
    }
    
    neighbours.emplace_back(a3o2_1,a3o2_2,a4_1,z-1);
    neighbours.emplace_back(a3o2_1,a3o2_2,a4_1,z+1);
    
    return neighbours;
}
#endif


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



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


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;
}
class Compare
{
public:
    bool operator() (const VVertex& v1 , const VVertex& v2 )
    {
        return get<0>(v1)>get<0>(v2);
    }
};

LI_t vector_to_int(vector<short> coords,int mult)
{
    LI_t total = 0;
    int i=0;
    for (auto it = coords.rbegin(); it != coords.rend(); ++it)
    {
        total+=pow(mult,i)* (*it);
        i++;
    }
    return total;
}

vector<short> int_to_vector(LI_t num,short mult)
{
        if(num==0)
        {
            return vector<short>{0};
        }
        std::vector<short> digits (floor(log(num)/log(mult))+1);
        int i=0;
        while(num > 0)
        {
            digits[digits.size()-1-i] = num%mult;
            num /= mult;
            i++;
        }

        return digits;
}



vector<vector<short>> int_to_vector(vector<LI_t> nums,int mult)
{
    vector<vector<short>> all(nums.size());
    for (int i=0;i<nums.size();i++)
    {
        all[i] = int_to_vector(nums[i],mult);
    }
    return all;
}

template<typename T>
ostream& operator<< (ostream& out, const vector<T>& v) {
    out << "{";
    size_t last = v.size() - 1;
    for(size_t i = 0; i < v.size(); ++i) {
        out << v[i];
        if (i != last)
            out << ", ";
    }
    out << "}";
    return out;
}

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;
    #if TYPE==3
    Vertex initial_vertex =  std::make_tuple(rand_long(),rand_long(),0,0);
    #endif
    #if TYPE==4
    Vertex initial_vertex =  std::make_tuple(rand_long(),rand_long(),rand_long(),0);
    #endif
    
    
    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),get<4>(currentE));

        
        
        
        frontier.pop();
        
        //current_marker=1;
        
        //Hmult3
        #if TYPE==3
        current_marker = (double)visited/(double)visited_number_running+1.005*pow(visited,-0.461);
        #endif
        
        //Hmult4
        #if TYPE==4
        current_marker = (double)visited/(double)visited_number_running+0.5578*pow(visited,-0.41583);
        #endif

        
        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),get<3>(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;
    
}


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

tuple<int,int,int> LeathRun(int n, double p)
{
    Vertex initial_vertex =  std::make_tuple(rand_long(),rand_long(),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);
}


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