//#ifndef _PolySim
//#define _PolySim

#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/connected_components.hpp>

typedef boost::adjacency_list< boost::vecS, boost::vecS, 
    boost::undirectedS, int> Graph;

//Edge iterator for or graph
typedef boost::graph_traits<Graph>::edge_iterator edge_iterator;

template<class T, int N>
using boostmat = typename boost::multi_array<T,N>;

typedef boost::multi_array<int, 2> array_type;
typedef boost::multi_array_types::index_range range_t;

using namespace std;

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

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

// Function to divide a vector/array by a scalar
template <typename Type>
vector<Type> div_scalar(vector<Type> &V, Type S) {
    vector<Type> Ans;
    Ans.reserve(V.size());
 
    for(typename vector<Type>::iterator it = V.begin(); it != V.end(); ++it) {
        Ans.push_back(*it/S);
    }
    return Ans;
}

// Function to perform a cumulative sum of the vector
template <typename Type>
vector<Type> cumsum(boostmat<Type,1> &V) {
    vector<Type> Ans;
    
    Type running_total = 0;
    for(int i = 0; i != V.shape()[0]; ++i) {
        running_total += V[i];
        Ans.push_back(running_total);
    }

    return Ans;
}

// Function to perform a sum of the vector
template <typename Type>
Type sum(vector<Type> &V) {
    vector<Type> Ans;
    
    Type running_total = 0;
    for(int i = 0; i != V.size(); ++i) {
        running_total += V[i];
    }

    return running_total;
}


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

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

int flip_direction(int direction){
    return((direction+2)%4);
}

void dispNeighbourhoodMap(boostmat<int,2> true_lattice, vector<int> label_vec) {
    boostmat<int,2> label_lattice;
    label_lattice.resize(boost::extents[ true_lattice.shape()[0] ][ true_lattice.shape()[1] ] );
    for ( int i = 0; i < label_lattice.shape()[0]; i++ ) {
        for ( int j = 0; j < label_lattice.shape()[1]; j++ ) {
            if ( true_lattice[i][j] > -1) {
                label_lattice[i][j] = label_vec[ true_lattice[i][j] ];
            } else {
                label_lattice[i][j] = -1;
            }
        }
    }
    PrintBoostMatrix(label_lattice);
}
/********************
 Class: Neighbourhood
 ********************/

class Neighbourhood{
    public:
    array<int,4> neighbours;

    Neighbourhood();
    Neighbourhood(array<int,4> IN_neighbours);
    int unlabelledCount();
    Neighbourhood getFreeNeighbourhood();
    void print();

    bool operator<(const Neighbourhood b) const {
        
        //cout << N << " " << moment_inertia << endl;
        //cout << b.N << " " << b.moment_inertia << endl;
        //cout << endl;
        
        if (this -> neighbours[0] < b.neighbours[0]) {
            return(true);
        } else if (this -> neighbours[0] > b.neighbours[0]) {
            return(false);
        }
        
        if (this -> neighbours[1] < b.neighbours[1]) {
            return(true);
        } else if (this -> neighbours[1] > b.neighbours[1]) {
            return(false);
        }
        
        if (this -> neighbours[2] < b.neighbours[2]) {
            return(true);
        } else if (this -> neighbours[2] > b.neighbours[2]) {
            return(false);
        }
        
        
        if (this -> neighbours[3] < b.neighbours[3]) {
            return(true);
        } else if (this -> neighbours[3] > b.neighbours[3]) {
            return(false);
        }
        
        return(false);
    }

    bool operator>(const Neighbourhood b) const {
        
        //cout << N << " " << moment_inertia << endl;
        //cout << b.N << " " << b.moment_inertia << endl;
        //cout << endl;
        
        if (this -> neighbours[0] > b.neighbours[0]) {
            return(true);
        } else if (this -> neighbours[0] < b.neighbours[0]) {
            return(false);
        }
        
        if (this -> neighbours[1] > b.neighbours[1]) {
            return(true);
        } else if (this -> neighbours[1] < b.neighbours[1]) {
            return(false);
        }
        
        if (this -> neighbours[2] > b.neighbours[2]) {
            return(true);
        } else if (this -> neighbours[2] < b.neighbours[2]) {
            return(false);
        }
        
        
        if (this -> neighbours[3] > b.neighbours[3]) {
            return(true);
        } else if (this -> neighbours[3] < b.neighbours[3]) {
            return(false);
        }
        
        return(false);
    }
    
     bool operator==(const Neighbourhood b) const {
        if ((this -> neighbours[0] == b.neighbours[0]) 
            && (this -> neighbours[1] == b.neighbours[1]) 
            && (this -> neighbours[2] == b.neighbours[2])
            && (this -> neighbours[3] == b.neighbours[3]) ){
            return(true);
        } else {
            return(false);
        }
    }

};

Neighbourhood::Neighbourhood() {
    neighbours = {-1,-1,-1,-1};
}

Neighbourhood::Neighbourhood(array<int,4> IN_neighbours) {
    neighbours = IN_neighbours;
}

void Neighbourhood::print() {
    for ( int i = 0; i < 4; i++ ){
        if ( (neighbours[i] < 10) && (neighbours[i] > -1) ) {
            cout << " " << neighbours[i] << " " << flush;
        } else {
            cout << neighbours[i] << " " << flush;
        }           
    }
}

int Neighbourhood::unlabelledCount(){
    return( int(neighbours[0] == 0) + int(neighbours[1] == 0) +
        int(neighbours[2] == 0) + int(neighbours[3] == 0) );
}

Neighbourhood Neighbourhood::getFreeNeighbourhood(){ 
    //The free Neighbourhood is essentially a neighbourhood descriptor which is independent of rotations
        // It does this by returning the largest orientation under an elementwise sort
    set<Neighbourhood> free_neighbourhoods;

    free_neighbourhoods.insert(Neighbourhood({neighbours[0], neighbours[1], neighbours[2], neighbours[3]}));
	free_neighbourhoods.insert(Neighbourhood({neighbours[3], neighbours[0], neighbours[1], neighbours[2]}));      
	free_neighbourhoods.insert(Neighbourhood({neighbours[2], neighbours[3], neighbours[0], neighbours[1]}));
	free_neighbourhoods.insert(Neighbourhood({neighbours[1], neighbours[2], neighbours[3], neighbours[0]}));  

    return( *free_neighbourhoods.rbegin() );
}

void PrintContactMap(vector<Neighbourhood> contact_map) {
    for( int i = 0; i < contact_map.size(); i++ ){
        contact_map[i].print();
        cout << endl;
    }
}

/********************
 EndClass Neighbourhood
 ********************/

/********************
 Class: Minimizer
 ********************/

class Minimizer{
    public:
    vector<Neighbourhood> contact_map;
    vector<Neighbourhood> label_map;
    vector<int> label_vec;
    boostmat<int,2> true_lattice;
    
    int max_category_one_label; 
    int max_label;

    Minimizer(vector<Neighbourhood> IN_contact_map, boostmat<int,2> IN_true_lattice);
    void UpdateLabelMap();
    void labelCategoryOne();
    void labelCategoryTwo();
    tuple< int, vector<Neighbourhood>, vector<int> > Minimize();

    void dispLabelLattice();
};

Minimizer::Minimizer(vector<Neighbourhood> IN_contact_map, boostmat<int,2> IN_true_lattice) {
    contact_map = IN_contact_map;
    label_vec.resize(contact_map.size());
    label_map.resize(contact_map.size());
    true_lattice.resize(boost::extents[ IN_true_lattice.shape()[0] ][ IN_true_lattice.shape()[1] ]);
    true_lattice = IN_true_lattice;
    //fill(label_vec.begin(),label_vec.end(),0); Not required because integer vectors are always initialized at 0

    UpdateLabelMap();
}

void Minimizer::UpdateLabelMap() {
    for (int i = 0; i < label_map.size() ; i++) {
        for (int j = 0; j < 4; j++) {
            if ( contact_map[i].neighbours[j] > -1  ) {

                label_map[i].neighbours[j] = label_vec[ contact_map[i].neighbours[j] ];
            }
        }
    }
}

void Minimizer::labelCategoryOne() {
    // Find unlabelled locations with only one unlabelled neighbour and determine labels for them 
    int iteration_num = 0;
    bool category_one_solved = false;
    
    while (!category_one_solved) {
        int category_one_count = 0;
        vector<Neighbourhood> FreeNeighbourhoodList;
        for ( int i = 0; i < label_map.size(); i++ ) {
            
            if ( label_vec[i] == 0 && label_map[i].unlabelledCount() == 1 ) {
                category_one_count++;

                Neighbourhood free_neighbourhood = label_map[i].getFreeNeighbourhood();
                auto neighbour_list_iterator = find(FreeNeighbourhoodList.begin(),FreeNeighbourhoodList.end(),free_neighbourhood);
                if ( neighbour_list_iterator == FreeNeighbourhoodList.end() ) {
                    FreeNeighbourhoodList.push_back(free_neighbourhood);
                    label_vec[i] = iteration_num + FreeNeighbourhoodList.size();
                } else {
                    label_vec[i] = iteration_num + int(neighbour_list_iterator-FreeNeighbourhoodList.begin()) + 1;
                }
            }
        }
        UpdateLabelMap();
        // dispLabelLattice();
        // cout << endl;
        if ( category_one_count == 0 ) {
            category_one_solved = true;
        }
        iteration_num = iteration_num + FreeNeighbourhoodList.size();
    }

    max_category_one_label = iteration_num;
}

void Minimizer::labelCategoryTwo() {
    vector<int> unresolved_indices;

    // First step is to find all indexes that are still unlabelled
    for (int i = 0; i < label_vec.size(); i++) {
        if (label_vec[i] == 0) {
            unresolved_indices.push_back(i);
        }
    }
    
    bool solved = false;
    
    while( !solved ) {
        vector<int> last_label_vec = label_vec;
        vector<Neighbourhood> FreeNeighbourhoodList;
        for (int i: unresolved_indices) {
            Neighbourhood free_neighbourhood = label_map[i].getFreeNeighbourhood();
            auto neighbour_list_iterator = find(FreeNeighbourhoodList.begin(),FreeNeighbourhoodList.end(),free_neighbourhood);
            if ( neighbour_list_iterator == FreeNeighbourhoodList.end() ) {
                    FreeNeighbourhoodList.push_back(free_neighbourhood);
                    label_vec[i] = max_category_one_label + FreeNeighbourhoodList.size();
                } else {
                    label_vec[i] = max_category_one_label + int(neighbour_list_iterator-FreeNeighbourhoodList.begin()) + 1;
                }
        }
        
        UpdateLabelMap();

        //dispLabelLattice();
        //cout << endl;

        solved = true;
        for (int i: unresolved_indices) {
            if (label_vec[i] != last_label_vec[i]) {
                solved = false;
                break;
            }
        }
        
        if (solved) {
            max_label = max_category_one_label + FreeNeighbourhoodList.size();
        }
    }
}

tuple< int, vector<Neighbourhood>, vector<int> > Minimizer::Minimize(){
    labelCategoryOne();
    labelCategoryTwo();
    return { max_label, contact_map, label_vec }; 
}

void Minimizer::dispLabelLattice(){
    dispNeighbourhoodMap(true_lattice, label_vec);
}

/********************
 EndClass Minimizer
 ********************/


/********************
 Contact Map Functions
 ********************/

vector<Neighbourhood> latticetoMaximalContactMap(boostmat<int,2> lattice_input, int species_count){
    vector<Neighbourhood> contact_map;
    contact_map.resize(species_count); 
    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 ) {
                array<int,4> temp_neighbourhood;
                temp_neighbourhood[0] = lattice_input[i-1][j];
                temp_neighbourhood[1] = lattice_input[i][j+1];
                temp_neighbourhood[2] = lattice_input[i+1][j];
                temp_neighbourhood[3] = lattice_input[i][j-1];
                contact_map[ lattice_input[i][j] ] = Neighbourhood(temp_neighbourhood);
            }
        }
    }
    return(contact_map);
}

vector<Neighbourhood> unforcedContactMap(boostmat<int,2> lattice_input, int species_count){
    //same as contact map, but does not 'overwrite' neighbours 
    vector<Neighbourhood> contact_map;
    contact_map.resize(species_count); 
    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 ) {
                int map_idx = lattice_input[i][j];
                contact_map[map_idx].neighbours[0] = (lattice_input[i-1][j] > -1) ? lattice_input[i-1][j] : contact_map[map_idx].neighbours[0];
                contact_map[map_idx].neighbours[1] = (lattice_input[i][j+1] > -1) ? lattice_input[i][j+1] : contact_map[map_idx].neighbours[1];
                contact_map[map_idx].neighbours[2] = (lattice_input[i+1][j] > -1) ? lattice_input[i+1][j] : contact_map[map_idx].neighbours[2];
                contact_map[map_idx].neighbours[3] = (lattice_input[i][j-1] > -1) ? lattice_input[i][j-1] : contact_map[map_idx].neighbours[3];
            }
        }
    }
    return(contact_map);
}

Graph getAdjGraph(vector<Neighbourhood> contact_map){
    int vertex_num = contact_map.size();
    Graph adjacency_graph(vertex_num);
    for ( int i = 0; i < vertex_num; i++) {
        for (int j = 0; j < 4; j++) {
            int neighbour = contact_map[i].neighbours[j];
            if (  neighbour > -1 && neighbour < i ){
                boost::add_edge (i, neighbour, adjacency_graph);
            }
        }
    }
    return(adjacency_graph);
}

bool isConnected(Graph adjacency_graph){
    vector<int> component(boost::num_vertices(adjacency_graph));
    int num_components = boost::connected_components (adjacency_graph, &component[0]);

    return((num_components == 1));
}

vector<Neighbourhood> removeContactMapEdge( vector<Neighbourhood> contact_map, int vertex_1, int vertex_2) {
    vector<Neighbourhood> revised_contact_map = contact_map; 
    int vertex_2_idx;
    for( int k = 0; k < 4; k++) {
        if (contact_map[vertex_1].neighbours[k] == vertex_2) {
            vertex_2_idx = k;
            break;
        }
    }
    int vertex_1_idx = flip_direction(vertex_2_idx);
    revised_contact_map[vertex_1].neighbours[vertex_2_idx] = -1;
    revised_contact_map[vertex_2].neighbours[vertex_1_idx] = -1;
    return(revised_contact_map);
}

tuple< int, vector<Neighbourhood>, vector<int> > RecDFMinSearch( vector<Neighbourhood> contact_map, boostmat<int,2> lattice, int iter_distance) {
    Minimizer minimizer(contact_map, lattice);
    auto [running_min, running_min_contact, running_min_label] = minimizer.Minimize();
    
    Graph adjacency_graph = getAdjGraph(contact_map);

    //Iterate through all the edges
    std::pair<edge_iterator, edge_iterator> ei = boost::edges(adjacency_graph);
    edge_iterator start = ei.first;
    advance(start, iter_distance);
    for(edge_iterator edge_iter = start; edge_iter != ei.second; ++edge_iter) {
        Graph new_adjacency_graph = adjacency_graph;
        int vertex_src = source((*edge_iter),adjacency_graph);
        int vertex_tgt = target((*edge_iter),adjacency_graph);
        int new_distance = distance(ei.first,edge_iter);
        remove_edge( vertex_src, vertex_tgt, new_adjacency_graph );
        if ( isConnected( new_adjacency_graph )) {
            vector<Neighbourhood> new_contact_map = removeContactMapEdge(contact_map, vertex_src, vertex_tgt);
            
            auto [new_min, new_min_contact, new_min_label] = RecDFMinSearch( new_contact_map, lattice, new_distance );
            if ( new_min < running_min ){ 
                running_min = new_min;
                running_min_contact = new_min_contact;
                running_min_label = new_min_label;
            } 
        } 
    }

    return {running_min, running_min_contact, running_min_label};

    /*for edges in contact_map:
        if (is_connected( contact_map - edges)):
            (min,map) = test_and_minimize( contact_map - edges )
            if ( min < running_min ):
                min = running_min; 
                map = map;
    return( (min, map) )*/        
}

/********************
 End Contact Map Functions
 ********************/

boostmat<int,2> generateLattice(vector< pair<int,int> > NonZeroLoc, 
    int canvas_row_n, int canvas_col_n) {
    
    boostmat<int,2> lattice_space;
    lattice_space.resize(boost::extents[canvas_row_n][canvas_col_n]);
    
    for (int i = 0; i != lattice_space.shape()[0]; i++){
        for (int j = 0; j != lattice_space.shape()[1]; j++){
	        lattice_space[i][j] = -1;
        }
    }

    for (int i = 0; i < NonZeroLoc.size(); i++) {
        lattice_space[ NonZeroLoc[i].first ][ NonZeroLoc[i].second ] = i;
    }
    
    return(lattice_space);
}

tuple< int, vector<Neighbourhood>, vector<int> > polyomino_minimizer_no_neutral( vector<Neighbourhood> maximal_contact_map, boostmat<int,2> lattice ) {
    Minimizer minimizer(maximal_contact_map, lattice);
    auto [min_num, min_contact, min_label] = minimizer.Minimize();

    dispNeighbourhoodMap(lattice, min_label);
    
    return{min_num, min_contact, min_label};

}


tuple< int, vector<Neighbourhood>, vector<int> > polyomino_minimizer( vector< pair< int, int > > generator ) {
    /*int canvas_row_n = 10;
    int canvas_col_n = 10;
    int M = 0;
    boostmat<int,2> lattice;
    lattice.resize(boost::extents[canvas_row_n][canvas_col_n]);
    
    for (int i = 0; i != lattice.shape()[0]; i++){
        for (int j = 0; j != lattice.shape()[1]; j++){
	        lattice[i][j] = -1;
        }
    }
    
    for (int i = 3; i != lattice.shape()[0]-3; i++){
        for (int j = 3; j != lattice.shape()[1]-4; j++){
	        lattice[i][j] = M;
	        M++;
        }
    }*/

    /*vector< pair<int,int> > generator = { {2,1}, {1,2}, {2,2}, {3,2},
                                          {2,3}, {2,4}, {3,4}, {1,5}, {2,5} };
    */
    /*vector< pair<int,int> > generator = { {2,1}, {1,2}, {2,2}, {3,2},
                                          {2,3}, {1,4}, {2,4}, {3,4}, {1,5}, {3,5}, 
                                          {1,6}, {2,6},{3,6}, {2,7}, {1,8}, {2,8}, {2,9}, {2, 10}, 
                                          {3,10}, {4,10}, {5,10}, {4,9}, {4,11}};*/

    /*vector< pair<int,int> > generator = { {2,1}, {1,2}, {2,2}, {3,2},
                                          {2,3}, {1,4}, {2,4}, {3,4}, {1,5}, {3,5}, 
                                          {1,6}, {2,6},{3,6}, {2,7}, {1,8}, {2,8}, 
                                          {3,8}, {2, 9} };
    */

    boostmat<int,2> lattice = generateLattice( generator, 7, 13);
    PrintBoostMatrix(lattice);
    cout << endl;

    vector<Neighbourhood> maximal_contact_map = latticetoMaximalContactMap(lattice,generator.size());
    // PrintContactMap(maximal_contact_map);

    
    auto [min_num, min_contact, min_label] = RecDFMinSearch(maximal_contact_map,lattice,0);

    dispNeighbourhoodMap(lattice, min_label);

    return{min_num, min_contact, min_label};

}

tuple< int, vector<Neighbourhood>, vector<int> > polyomino_minimizer( vector<Neighbourhood> maximal_contact_map, boostmat<int,2> lattice ) {
    
    auto [min_num, min_contact, min_label] = RecDFMinSearch(maximal_contact_map,lattice,0);
    dispNeighbourhoodMap(lattice, min_label);

    return{min_num, min_contact, min_label};
}