#ifndef _PolySim
#define _SequencedPolyomino

#include <iostream>
#include <numeric>
#include <iterator>
#include <vector>
#include <unordered_map>
#include <math.h>
#include <random>
#include <chrono>
#include<unistd.h>

#include "boost/multi_array.hpp"
#include "boost/graph/adjacency_list.hpp"
#include "boost/graph/boyer_myrvold_planar_test.hpp"

#include "sequenced_polyomino_minimizer.hpp"

using namespace std;

typedef boost::adjacency_list<boost::setS, boost::vecS,
    boost::undirectedS> graph;

// cpp file contains definitions for the abstract GillespieSim class, but
// not the specific implementations 

/******************
 GENERIC FUNCTIONS
 ******************/

// Function to print arrays/vectors
template <class Type>
void print_arr_with_idx(vector<Type> V, int idx) {
    for(typename vector<Type>::iterator it = V.begin(); it != V.begin()+idx; ++it) {
        if (it != V.end()-1){
            cout << *it << " ";
	} else {
	    cout << *it;
	}
    }
    cout << flush;
}

template <class Type>
void printMatrix(boostmat<Type,2> lattice) {
    for (int i = 0; i != lattice.shape()[0]; i++){
        for (int j = 0; j != lattice.shape()[1]; j++){
            if ((lattice[i][j] > -1) && (lattice[i][j] < 10)) {
	        cout << " " << lattice[i][j] << " " << flush;
	    } else {
	        cout << lattice[i][j] << " " << flush;
	    }
        }
        cout << endl;
    }
}

void print_2d_cpp_array(array<int,2> V) {
    cout << "(" << V[0] << "," << V[1] << ")";
    cout << flush;
}

bool neighbourhoodSort(SequencedNeighbourhood neighbourhood_1, SequencedNeighbourhood neighbourhood_2){
    if (neighbourhood_1.neighbourhood_type != neighbourhood_2.neighbourhood_type) {
        return(neighbourhood_1.neighbourhood_type < neighbourhood_2.neighbourhood_type);
    } else {
        return(neighbourhood_1.label < neighbourhood_2.label);
    }
}

array<int,2> getLatticeCoordinate(int input){
    array<int,2> coordinate;
    switch (input) {
        case 0:
            coordinate = {-1,0};
            break;
        case 1:
            coordinate = {0,1};
            break;
        case 2:
            coordinate = {1,0};
            break;
        case 3:
            coordinate = {0,-1};
            break;
        default:
            throw runtime_error("getLatticeCoordinate: Unexpected lattice direction index");
   }
   
   return(coordinate);
}

array<int,2> findStartCoordinate(boostmat<int,2> lattice){
    for (int i = 0; i != lattice.shape()[0]; i++){
        for (int j = 0; j != lattice.shape()[1]; j++){
            if (lattice[i][j] == 0) {
                    array<int,2> coord_return = { i, j };
                    return(coord_return);
            }
        }
    }
    throw runtime_error("findStartCoordinate: starting block 0 not found!");
}

/*
boostmat<int,2> GenerateRandomLatticeint(){
    // 1 is forward, 2 is right, 3 is left
    std::random_device rd{}; 
    mt19937_64 mersenne_generator(rd());
}
*/

/****
 k-coloring utility functions
****/

/* A utility function to check if
   the current color assignment
   is safe for vertex v i.e. checks
   whether the edge exists or not
   (i.e, graph[v][i]==1). If exist
   then checks whether the color to
   be filled in the new vertex(c is
   sent in the parameter) is already
   used by its adjacent
   vertices(i-->adj vertices) or
   not (i.e, color[i]==c) */
bool isSafe(int v, boostmat<bool,2> graph,
            vector<int> color, int c, int V)
{
    for(int i = 0; i < V; i++)
        if (graph[v][i] && c == color[i])
            return false;
             
    return true;
}
 
/* A recursive utility function
to solve m coloring problem */
bool graphColoringUtil(boostmat<bool,2> graph, int m,
                       vector<int>& color, int v, int V)
{
    /* base case: If all vertices are
       assigned a color then return true */
    if (v == V)
        return true;
 
    /* Consider this vertex v and
       try different colors */
    for(int c = 1; c <= m; c++)
    {
         
        /* Check if assignment of color
           c to v is fine*/
        if (isSafe(v, graph, color, c, V))
        {
            color[v] = c;
 
            /* recur to assign colors to
               rest of the vertices */
            if (graphColoringUtil(
                graph, m, color, v + 1, V) == true)
                return true;
 
            /* If assigning color c doesn't
               lead to a solution then remove it */
            color[v] = 0;
        }
    }
 
    /* If no color can be assigned to
       this vertex then return false */
    return false;
}

/* A utility function to print solution */
void printSolution(vector<int> color)
{
    int V = color.size();
    cout << "Solution Exists:"
         << " Following are the assigned colors"
         << "\n";
    for(int i = 0; i < V; i++)
        cout << " " << color[i] << " ";
         
    cout << "\n";
}

/* This function solves the m Coloring
   problem using Backtracking. It mainly
   uses graphColoringUtil() to solve the
   problem. It returns false if the m
   colors cannot be assigned, otherwise
   return true and prints assignments of
   colors to all vertices. Please note
   that there may be more than one solutions,
   this function prints one of the
   feasible solutions.*/
bool graphColoring(boostmat<bool,2> graph, int m)
{
     
    // Initialize all color values as 0.
    // This initialization is needed
    // correct functioning of isSafe()
    
    int V = graph.shape()[0];
    
    vector<int> color(V,0);
 
    // Call graphColoringUtil() for vertex 0
    if (graphColoringUtil(graph, m, color, 0, V) == false)
    {
        cout << "Solution does not exist";
        return false;
    }
 
    // Print the solution
    // printSolution(color);
    return true;
}


/*****************
  Helper functions to create a boost multi-array
******************/
template <std::size_t... vs> struct ArraySize;
template <std::size_t v, std::size_t... vs> struct ArraySize<v, vs...>
{ static constexpr std::size_t size = v * ArraySize<vs...>::size; };
template <> struct ArraySize<>
{ static constexpr std::size_t size = 1; };

// Creates your multi_array
template <typename T, int... dims>
boost::multi_array<T, sizeof...(dims)>
makeMultiArray(std::initializer_list<T> l)
{
  constexpr std::size_t asize = ArraySize<dims...>::size;
  assert(l.size() == asize); // could be a static assert in C++14

  // Dump data into a vector (because it has the right kind of ctor)
  const std::vector<T> a(l);
  // This can be used in a multi_array_ref ctor.
  boost::const_multi_array_ref<T, sizeof...(dims)> mar(
    &a[0],
    std::array<int, sizeof...(dims)>{dims...});
  // Finally, deep-copy it into the structure we can return.
  return boost::multi_array<T, sizeof...(dims)>(mar);
}

/****************
 Class: SequencedNeighbourhood
 ****************/
 
 /*Stores locations, backbones and neighbourhood details of a specific block*/
 
SequencedNeighbourhood::SequencedNeighbourhood() {
    block_id = -1;
    row = -1;
    col = -1;
    backbone_type = 0;
    neighbourhood_type = 0;
    label = 0;
}
 
void SequencedNeighbourhood::clear_neighbourhood() {
    block_id = -1;
    row = -1;
    col = -1;
    neighbour_locations.clear();
    backbone_type = 0;
    neighbourhood_type = 0;
    label = 0;
}
 
void SequencedNeighbourhood::reset_labels() {
    neighbourhood_type = 0;
    label = 0;
}

void SequencedNeighbourhood::print() {
    cout << block_id << "\t" <<
         "(" << row << "," << col << ")\t" <<
         backbone_type << "\t{";
         
    for(auto it = neighbour_locations.begin(); it != neighbour_locations.end(); ++it) {
        if (it == neighbour_locations.end()-1){
            cout << *it << "}    \t";
	    } else {
	        cout << *it << ",";
	    }
    }
    
    cout << neighbourhood_type << "\t";
    cout << label-1;
}

/****************
 Class: SequencedPolyominoMinimizer
 ****************/

// Class that performs the self-avoiding walk. Can be made to return energies or validity 
SequencedPolyominoMinimizer::SequencedPolyominoMinimizer(boostmat<int,2> lattice_input){
    canvas_row_num = lattice_input.shape()[0];
    canvas_col_num = lattice_input.shape()[1];
    lattice.resize(boost::extents[canvas_row_num][canvas_col_num]);
    lattice = lattice_input;
    start_coordinate = findStartCoordinate(lattice);
    
    labelled_lattice.resize(boost::extents[lattice_input.shape()[0]][lattice_input.shape()[1]]);
    
    for (int i = 0; i != lattice_input.shape()[0]; i++){
        for (int j = 0; j != lattice_input.shape()[1]; j++){
            if (lattice_input[i][j] > -1) {
                labelled_lattice[i][j] = 1; // The integer 1 is reserved for unlabelled spots, so all labels start from 2 internally! The user always sees them increasing from 2.
            } else {
                labelled_lattice[i][j] = 0;
            }
        }
    }
    
    // Find the number of blocks
    int maximum_sweep = 0;
    for (int i = 0; i != lattice_input.shape()[0]; i++){
        for (int j = 0; j != lattice_input.shape()[1]; j++){
            if (lattice_input[i][j] > maximum_sweep) {
                maximum_sweep = lattice_input[i][j];
            }
        }
    }
    
    block_num = maximum_sweep+1;
    highest_label = 1;
    
    getNeighbourhoodClassifiers();
    // set all to -1
    
    //TODO: ADD FUNCTIONS TO CHECK BOUNDARIES ARE ALL -1
    //TODO: ADD FUNCTIONS TO CHECK CONSISTENCY OF BLOCK NUMBERING
}

SequencedPolyominoMinimizer::SequencedPolyominoMinimizer(boostmat<int,2> lattice_input, array<int,2> start_coord_input){
    canvas_row_num = lattice_input.shape()[0];
    canvas_col_num = lattice_input.shape()[1];
    lattice.resize(boost::extents[canvas_row_num][canvas_col_num]);
    lattice = lattice_input;
    
    start_coordinate = start_coord_input;
    
    labelled_lattice.resize(boost::extents[lattice_input.shape()[0]][lattice_input.shape()[1]]);
    
    for (int i = 0; i != lattice_input.shape()[0]; i++){
        for (int j = 0; j != lattice_input.shape()[1]; j++){
            if (lattice_input[i][j] > -1) {
                labelled_lattice[i][j] = 1; // The integer 1 is reserved for unlabelled spots, so all labels start from 2 internally! The user always sees them increasing from 2.
            } else {
                labelled_lattice[i][j] = 0;
            }
        }
    }
    
    // Find the number of blocks
    int maximum_sweep = 0;
    for (int i = 0; i != lattice_input.shape()[0]; i++){
        for (int j = 0; j != lattice_input.shape()[1]; j++){
            if (lattice_input[i][j] > maximum_sweep) {
                maximum_sweep = lattice_input[i][j];
            }
        }
    }
    
    block_num = maximum_sweep+1;
    highest_label = 1;
    //TODO: ADD FUNCTIONS TO CHECK BOUNDARIES ARE ALL -1
    //TODO: ADD FUNCTIONS TO CHECK CONSISTENCY OF BLOCK NUMBERING
    
    getNeighbourhoodClassifiers();

}

void SequencedPolyominoMinimizer::getNeighbourhoodClassifiers(){
    // Need to change this so that the maximum is dependent on the current maximum
    boostmat<int,2> new_labelled_lattice; 
    new_labelled_lattice.resize(boost::extents[labelled_lattice.shape()[0]][labelled_lattice.shape()[1]]);
    
    int block_top_neighbour, block_right_neighbour, block_bot_neighbour, block_left_neighbour;
    
    int i, j; // lattice row and column
    
    int bool_eval; // int flag that is implemented as a sum of boolean expressions
                       // equal to one if the 
    
    // Neighbourhood vector and temp variable
    SequencedNeighbourhood temp_neighbourhood;
    
    // Backbone related variables
    int previous_backbone_dir = 4;
    int next_backbone_dir;
    int next_block;
    int backbone_type;
    int dir_to_penultimate;
    
    // Temporary variables for storing non-backbone neighbour locations of the given block
    int neighbour_1_dir, neighbour_2_dir, neighbour_3_dir;
    array<int,2> n1_coordinates, n2_coordinates, n3_coordinates;

    vector<int> temp_neighbours;
    
    //set current coordinate to starting coordinate and declare variable that tells you where
        // the next block should be
    array<int,2> coordinate = start_coordinate;
    array<int,2> next_coord_shift;
    
    bool previous_empty = false;
    bool empty = false;
    
    int neighbourhood_id; // Used to identify the kind of neighbourhood
    
    // cout << block_num << endl;
    // For loop for identifying backbones only. No classification of blocks yet.
    for (int k = 0; k < block_num; k++) {
        i = coordinate[0]; //row number of the current coordinate
        j = coordinate[1]; // column number of the current coordinate

        if (lattice[i][j] != k) {
            cout << "At block number " << k << ": " << endl;
            throw runtime_error("SequencedPolyominoMinimizer.ResolveLabels: invalid lattice location!");
        }
        
        neighbour_1_dir = 0;
        neighbour_2_dir = 0;
        neighbour_3_dir = 0;
        temp_neighbours.clear(); // clear coordinate and neighbourhood buffers
        temp_neighbours.reserve(3);

        next_block = lattice[i][j] + 1;

        // PART I: IDENTIFY BACKBONE TYPE

        block_top_neighbour = lattice[i-1][j];
        block_right_neighbour = lattice[i][j+1];
        block_bot_neighbour = lattice[i+1][j];
        block_left_neighbour = lattice[i][j-1];

        // a. Check that exactly one lattice neighbour is the next block in the sequence
        // and identify the direction we would need to move on the lattice (lattice direction)
        // to get to this next block

        bool_eval = int(block_top_neighbour == next_block) + int(block_right_neighbour == next_block) 
        + int(block_bot_neighbour == next_block) + int(block_left_neighbour == next_block);
        
        if (lattice[i][j] == block_num - 1) {
            if (bool_eval != 0) {
                throw runtime_error("SequencedPolyominoMinimizer: block numbering error in terminal block!");
            } else {
                next_backbone_dir = 4; 
                    // next backbone dir is used to decide the non-covalent neighbours
                    // The direction to the previous block can be used for this for a terminal block
                next_coord_shift[0] = 0;
                next_coord_shift[1] = 0;
            }
        } else {
        if (bool_eval != 1) {
                throw runtime_error("SequencedPolyominoMinimizer: block numbering error in non-terminal block!");
            } else {
                next_backbone_dir = int(block_top_neighbour == next_block)*0 + int(block_right_neighbour == next_block)*1 
                + int(block_bot_neighbour == next_block)*2 + int(block_left_neighbour == next_block)*3;

                next_coord_shift = getLatticeCoordinate(next_backbone_dir);
                // 4 lattice directions. Indexed as follows:
                // 0 up, 1 right, 2 down, 3 left
            }
        }
        
        // cout << next_backbone_dir << endl;
        // b. Based on the last and next directions taken to move along the sequence, 
        // decide on the type of backbone we need. 

        if (next_backbone_dir < 4 && previous_backbone_dir < 4) {
            backbone_type = (next_backbone_dir + 4 - previous_backbone_dir) % 4 + 1;
            // Backbone type is DISTINCT from the absolute directions on the lattice. 
            // Backbone type is more relative. Starting from the previous block, 
            // connect to the current block with a backbone. Then see what direction
            // to go for the next block assuming the previous block is the 'back'.
            // 0 is ignored because it is not convenient for integer sorting
            // 1 is forbidden here as it would require a backbone that moves up then down,
            // or left then right, (in lattice direction) or vice versa
            // 2 is a left facing backbone
            // 3 is a forward facing backbone
            // 4 is a right facing backbone
            if (backbone_type == 1) {
                throw runtime_error("SequencedPolyominoMinimizer: forbidden backbone detected!");
            }
        } else {
            if (k == 0) {
                backbone_type = 1;
            } else {
                backbone_type = 5;
            } 
            // 1 is therefore reserved for a starting backbone
            // 5 is reserved for the last backbone
        }


        // c. Decide on the lattice directions of NON-BACKBONE 
        //neighbours relative to the current block 
        // To revise, 0 if the neighbour is up
        // 1 if the neighbour is to the right
        // 2 if the neighbour is down
        // 3 if the nieghbour is to the left
        
        if (backbone_type == 1 || backbone_type == 5) {
            if (k == 0) {
                neighbour_1_dir = (next_backbone_dir+1) % 4;
                neighbour_2_dir = (next_backbone_dir+2) % 4;
                neighbour_3_dir = (next_backbone_dir+3) % 4;
            } else {
                dir_to_penultimate = previous_backbone_dir;
                    // Use the direction to the previous backbone to determine directions
                neighbour_1_dir = (dir_to_penultimate+1) % 4;
                neighbour_2_dir = (dir_to_penultimate+2) % 4;
                neighbour_3_dir = (dir_to_penultimate+3) % 4;  
            }
        } else if (backbone_type == 2) {
            neighbour_1_dir = (next_backbone_dir+1) % 4;
            neighbour_2_dir = (next_backbone_dir+2) % 4;
            neighbour_3_dir = 4;
        } else if (backbone_type == 4) {
            neighbour_1_dir = (next_backbone_dir+2) % 4;
            neighbour_2_dir = (next_backbone_dir+3) % 4;
            neighbour_3_dir = 4;
        } else if (backbone_type == 3) {
            neighbour_1_dir = (next_backbone_dir+1) % 4;
            neighbour_2_dir = (next_backbone_dir+3) % 4; 
            neighbour_3_dir = 4;
        } else {
            throw runtime_error("SequencedPolyominoMinimizer: unknown backbone type, cannot find neighbours unfound.");
        }
        
        // Save the previous block in the sequence as a neighbour ONLY IF there are other neighbours
            // This is to prevent connecting with the previous block in the sequence with a non-backbone connection
            // If there are no other neighbours, then it would not matter as the non-backbone connections cannot interact
                // with any other block on the non-backbone connections anyway 
        
        
        // Save all the non-covalent neighbours into a vector
        n1_coordinates = getLatticeCoordinate(neighbour_1_dir);
        temp_neighbours.push_back(lattice[i+n1_coordinates[0]][j+n1_coordinates[1]]);
        n2_coordinates = getLatticeCoordinate(neighbour_2_dir);
        temp_neighbours.push_back(lattice[i+n2_coordinates[0]][j+n2_coordinates[1]]);
        if (neighbour_3_dir < 4) { // if terminal block
            n3_coordinates = getLatticeCoordinate(neighbour_3_dir);
            temp_neighbours.push_back(lattice[i+n3_coordinates[0]][j+n3_coordinates[1]]);
            empty = true;
        } else { // for non-terminal blocks
            empty = ((temp_neighbours[0] == -1) && (temp_neighbours[1] == -1));
            // Save the previous block in the sequence as a neighbour ONLY IF there are other neighbours
            // This is to prevent connecting with the previous block in the sequence with a non-backbone connection
            // This is only relevant for non-terminal blocks because the first terminal block will have no depencies on previous blocks,
                // so it can always be made to match the last assuming the backbone and other neighbours are equivalent
            if ( empty || previous_empty ){
                // If there are no other neighbours, then it would not matter as the non-backbone connections cannot interact
                // with any other block on the non-backbone connections. So we set this uniformly to -1 (no neighbour)
                
                // If the previous block has no neighbours, we set the neighbour to -1 and we will later allow this block to take on the same
                // label as a block with a different previous neighbour but equivalent in everything else. 
                temp_neighbours.push_back(-1);
            } else {
                temp_neighbours.push_back(k-1);
            }
        }
        
        // Fill in the neighbourhood object and push the object into the vector
        temp_neighbourhood.block_id = lattice[i][j];
        temp_neighbourhood.row = i;
        temp_neighbourhood.col = j;
        temp_neighbourhood.neighbour_locations = temp_neighbours;
        temp_neighbourhood.backbone_type = backbone_type;
        temp_neighbourhood.label = 1;
        
        neighbourhood_classifiers.push_back(temp_neighbourhood);
        
        // Update to the next coordinate in the sequence
        coordinate[0] = coordinate[0] + next_coord_shift[0];
        coordinate[1] = coordinate[1] + next_coord_shift[1];
        previous_backbone_dir = (next_backbone_dir+2) % 4;
        
        previous_empty = empty;
    }

}

bool SequencedPolyominoMinimizer::getInteractionConnectivity(){
    // Matrix such that M(i,j) is greater than -1 iff block i and block j are interacting
    // They are enumerated from 0 to etc..
    boostmat<int,2> interaction_enum;
    int mat_size = neighbourhood_classifiers.size();
    interaction_enum.resize(boost::extents[mat_size][mat_size]);
    
    // set all to -1
    for (int i = 0; i < mat_size ;i++) {
        for (int j = 0; j < mat_size ;j++) {
            interaction_enum[i][j] = -1;
        }
    }
    
    int counter = 0;
    int neighbour;
    
    for (int i = 0; i < mat_size ;i++) {
        for (int j = 0; j < neighbourhood_classifiers[i].neighbour_locations.size() ;j++) {
            neighbour = neighbourhood_classifiers[i].neighbour_locations[j];
            if( neighbour > -1 && neighbour > i){
                interaction_enum[i][neighbour] = counter;
                interaction_enum[neighbour][i] = counter;
                counter++;
            }
        }
    }
    
    int max_interaction = counter;

    interaction_graph.resize(boost::extents[max_interaction][max_interaction]);
    for (int i = 0; i < max_interaction ;i++) {
        for (int j = 0; j < max_interaction ;j++) {
            interaction_graph[i][j] = false;
        }
    }
    
    
    
    for (int i = 0; i < mat_size ;i++) {
        for (int j = 0; j < mat_size ;j++) {
            if (interaction_enum[i][j] > -1) {
                for (int k = 0; k < max_interaction ;k++){
                    if ( i < mat_size-1){
                        if(interaction_enum[i+1][k] > -1){
                            interaction_graph[interaction_enum[i][j]][interaction_enum[i+1][k]] = true; 
                        }
                    }
                    if ( i > 0){
                        if(interaction_enum[i-1][k] > -1){
                            interaction_graph[interaction_enum[i][j]][interaction_enum[i-1][k]] = true; 
                        }
                    }
                }
                for (int k = 0; k < max_interaction ;k++){
                    if ( j < mat_size-1){
                        if(interaction_enum[k][j+1] > -1){
                            interaction_graph[interaction_enum[k][j+1]][interaction_enum[i][j]] = true; 
                        }
                    } 
                    if (j > 0 ){
                        if(interaction_enum[k][j-1] > -1){
                            interaction_graph[interaction_enum[k][j-1]][interaction_enum[i][j]] = true; 
                        }
                    }
                }
            }
        }
    }
    
    
    
    graph g;
    
    for (int i = 0; i < max_interaction; i++){
        for(int j = 0; j < i; j++){
            if (interaction_graph[i][j] ) {
                boost::add_edge(i,j,g);
            }
        }
    }
    // printMatrix(interaction_enum);
    // printMatrix(interaction_graph);
    bool sol = graphColoring(interaction_graph,4);
    
    std::pair<graph::edge_iterator,
    graph::edge_iterator> es = boost::edges(g);
    
    //cout << endl; 
    
    //std::copy(es.first, es.second,
    //std::ostream_iterator<graph::edge_descriptor>{std::cout, "\n"});
    
    // bool is_planar = boyer_myrvold_planarity_test(g);
    // cout << "Is Planar? " << is_planar << endl;
    return(sol);
}

void SequencedPolyominoMinimizer::printLattice(){
    for (int i = 0; i != lattice.shape()[0]; i++){
        for (int j = 0; j != lattice.shape()[1]; j++){
            if ((lattice[i][j] > -1) && (lattice[i][j] < 10)) {
	        cout << " " << lattice[i][j] << " " << flush;
	    } else {
	        cout << lattice[i][j] << " " << flush;
	    }
        }
        cout << endl;
    }
}

void SequencedPolyominoMinimizer::printLabelledLattice(){
    for (int i = 0; i != labelled_lattice.shape()[0]; i++){
        for (int j = 0; j != labelled_lattice.shape()[1]; j++){
	        if ((labelled_lattice[i][j] > 0) && (labelled_lattice[i][j] < 11)) {
	            cout << " " << labelled_lattice[i][j]-1 << " " << flush;
	        } else if (labelled_lattice[i][j] > 0){
	            cout << labelled_lattice[i][j]-1 << " " << flush;
	        } else {
	            cout << " "<< 0 << " " << flush;
	        } 
        }
        cout << endl;
    }
}



/****************
 Class: AWSimulator
 ****************/

// Class that manages the simulation of avoiding walkers. Each AWSimulator class will 
// simulate all possible sequences and return the equilibrium conformation as well as 
// other outputs such as the 



#endif