#include "nauty.h"

#include <iostream>

#include <iomanip>

#include <fstream>

#include <string>

#include <map>

#include <algorithm>

#include <chrono>

#define BOOST_UBLAS_NDEBUG 1

#include <boost/dynamic_bitset.hpp> // Use Boost's dynamic_bitsets for sets

#include <boost/numeric/ublas/vector.hpp> // Use Boost's basic linear algebra (ublas) vectors

#include <boost/numeric/ublas/matrix.hpp> // Use Boost's basic linear algebra (ublas) matrices

#include <boost/numeric/ublas/io.hpp> // Use Boost's basic linear algebra (ublas) io routines

#include <boost/numeric/ublas/lu.hpp>

namespace ublas = boost::numeric::ublas;

typedef boost::dynamic_bitset<>::size_type size_type;

size_type npos = boost::dynamic_bitset<>::npos;

std::multimap<int,std::vector<setword>> atoms_cg;

std::multimap<int,std::vector<setword>> LD;

int simplicial_facet = 0; // 2L polytopes with a simplicial facet

int cs = 0; // Centrally-symmetric 2L polytopes

int stab = 0; // Stable set polytopes of a perfect graph

int n_suspensions = 0; // 2L suspensions

int n_polar = 0; // Polar 2L polytopes

std::ofstream myfile;

std::vector<boost::dynamic_bitset<>> slabs_sat_by_point;

std::vector<boost::dynamic_bitset<>> points_sat_by_slab;

bool sqsubseteq(boost::dynamic_bitset<> A, boost::dynamic_bitset<> B) {

}

boost::dynamic_bitset<> inc(boost::dynamic_bitset<> A, int i) {

}

boost::dynamic_bitset<> discreteconvexhull_cl(const boost::dynamic_bitset<> &A, boost::dynamic_bitset<> &B) {

}

boost::dynamic_bitset<> incompatibility_cl(const std::vector<boost::dynamic_bitset<>> &incompatibility_adjM,boost::dynamic_bitset<> &A) {

}

int hash_matrix(const std::vector<boost::dynamic_bitset<>> &S, int D) {

}

std::vector<setword> canonicize(const std::vector<boost::dynamic_bitset<>> &S) {

std::vector<setword> cg_vec;

// Make sure that the slack matrix has at least one row and column

assert (S.size() > 0);

assert (S[0].size() > 0);

// Get the number of rows and columns of current atom

int num_rows = (int) S.size();

int num_cols = (int) S[0].size();

// Initializations for Nauty

DYNALLSTAT(int,lab,lab_sz);

DYNALLSTAT(int,ptn,ptn_sz);

DYNALLSTAT(int,orbits,orbits_sz);

DYNALLSTAT(int,map,map_sz);

DYNALLSTAT(graph,g,g_sz);

DYNALLSTAT(graph,cg,cg_sz);

static DEFAULTOPTIONS_GRAPH(options);

statsblk stats;

// Select option for canonical labelling

options.getcanon = TRUE;

// Select option for custom partition

options.defaultptn = FALSE;

int m, n;

n = num_rows + num_cols;

m = SETWORDSNEEDED(n);

nauty_check(WORDSIZE,m,n,NAUTYVERSIONID);

// Allocate memory for the graph

DYNALLOC1(int,lab,lab_sz,n,"malloc");

DYNALLOC1(int,ptn,ptn_sz,n,"malloc");

DYNALLOC1(int,orbits,orbits_sz,n,"malloc");

DYNALLOC1(int,map,map_sz,n,"malloc");

DYNALLOC2(graph,g,g_sz,n,m,"malloc");

DYNALLOC2(graph,cg,cg_sz,n,m,"malloc");

// Empty the graph

EMPTYGRAPH(g,m,n);

// Build the custom partition for the graph

for (int i = 0; i < num_rows+num_cols; i++) {

lab[i] = i;

if (i != num_rows-1 && i != num_rows+num_cols-1)

ptn[i] = 1;

else

ptn[i] = 0;

}

// Build the edges of the nonincidence graph

// Loop through all the entries of the slack matrix and add an edge when there is a one

for (int i = 0; i < num_rows; i++) {

for (int j = 0; j < num_cols; j++) {

if ((S[i])[j] == 1)

ADDONEEDGE(g,i,j+num_rows,m);

}

}

// Obtain canonical graph from Nauty

densenauty(g,lab,ptn,orbits,&options,&stats,m,n,cg);

// Make std::vector<setword> from canonical graph

for (size_t k = 0; k < m*(size_t)n; ++k)

cg_vec.push_back(cg[k]);

// Clean up

DYNFREE(lab,lab_sz);

DYNFREE(ptn,ptn_sz);

DYNFREE(orbits,orbits_sz);

DYNFREE(map,map_sz);

DYNFREE(g,g_sz);

DYNFREE(cg,cg_sz);

return cg_vec;

}

bool checksimplicialcore(const std::vector<boost::dynamic_bitset<>> &S, int d) {

bool test = true;

// Check that S has at least d rows

assert (S.size() >= d);

for (int i=0; i<=d && test; ++i) {

// Check that the ith row of S has at least d columns

assert (S[i].size() >= d);

for (int j=i; j<=d && test; ++j) {

if (i == j)

test = ((S[i])[j] == 1);

else if (i < j)

test = ((S[i])[j] == 0);

}

}

return test;

}

ublas::matrix<int> extractM(std::vector<boost::dynamic_bitset<>> S, int d) {

ublas::matrix<int> M(d,d);

// Check that S has at least d rows

assert (S.size() >= d);

for (int i=0; i<d; ++i) {

// Check that the ith row of S has at least d columns

assert (S[i].size() >= d);

for (int j=0; j<d; ++j) {

if (i == 0 && j == 0)

M(i,j) = 1;

else if (i == 0 || j == 0)

M(i,j) = 0;

else

M(i,j) = (S[i-1])[j-1];

}

}

return M;

}

bool invertM(const ublas::matrix<int> &M, ublas::matrix<int> &Minv) {

// Create a duplicate of the input matrix

ublas::matrix<size_type> W = M;

assert (W.size1() == W.size2());

size_type d = W.size1();

// Create the permutation matrix

ublas::permutation_matrix<size_type> P(d);

Minv.resize(d,d);

// Assign the identity matrix to the inverse

Minv.assign(ublas::identity_matrix<int>(d));

// LU factorization and substitution

size_type res = lu_factorize(W, P);

if (res != 0)

return false;

lu_substitute(W,P,Minv);

return true;

}

std::vector<boost::dynamic_bitset<>> construct_slack_matrix(const std::vector<ublas::vector<int>> &base_H,const std::vector<ublas::vector<int>> &ground_set_H,const boost::dynamic_bitset<> &A,const boost::dynamic_bitset<> &B,const std::vector<ublas::vector<int>> &slabs, const std::vector<boost::dynamic_bitset<>> &S, std::vector<boost::dynamic_bitset<>> &S_new, int &D) {

std::vector<boost::dynamic_bitset<>> all_ineqs;

for (size_type i = B.find_first(); i != npos; i = B.find_next(i)){

// constuct each row of S with column indexed by e_1, the vertices of the base and then the remaining ones

// in this way the first D+1 columns will be the H-embedding of the canonical affine base of R^D.

boost::dynamic_bitset<> S_row;

int s = ublas::inner_prod(ground_set_H[0],slabs[i]);

bool slack1 = (s != 0);

S_row.push_back(slack1);

for (std::vector<ublas::vector<int>>::const_iterator it=base_H.begin(); it!=base_H.end(); it++) {

int s = ublas::inner_prod(*it,slabs[i]);

bool slack1 = (s != 0);

S_row.push_back(slack1);

}

for (size_type j = A.find_next(0); j != npos; j = A.find_next(j)){

int s = ublas::inner_prod(ground_set_H[j],slabs[i]);

bool slack1 = (s != 0);

S_row.push_back(slack1);

}

if ((~S_row).count() >= D)

all_ineqs.push_back(S_row);

if (S_row.count() >= D)

all_ineqs.push_back(~S_row);

}

// check maximality of rows

std::vector<boost::dynamic_bitset<>>::iterator it1,it2;

for (it1=all_ineqs.begin(); it1!=all_ineqs.end(); ++it1) {

bool is_maximal = true;

for (it2=all_ineqs.begin(); it2!=all_ineqs.end() && is_maximal; ++it2) {

if ((*it2).is_subset_of(*it1) && it2!=it1)

is_maximal = false;

}

if (is_maximal)

S_new.push_back(*it1);

}

// rearranging rows of S

int n_row = 0;

for (std::vector<boost::dynamic_bitset<>>::iterator it = S_new.begin()+n_row+1; it != S_new.end(); it++) {

//for (std::vector<boost::dynamic_bitset<>>::iterator print = S.begin(); print != S.end(); print++)

// std::cout << *print << std::endl;

//std::cout << n_row + 1 << std::endl;

bool accept = true;

for (size_type j = 0; j < S[0].size() && n_row < D && accept; j++)

accept = ((*it).test(j+1) == (S[n_row]).test(j));

if (accept) {

swap(*it, S_new[n_row+1]);

++n_row;

}

}

/*

return S_new;

}

bool istwolevelpolytope(std::vector<boost::dynamic_bitset<>> S, int &D) {

// First test: check that every column contains at least D zeros

// by construction, every row of S contains at least D zeros

bool accept = true;

for (size_type it_cols = 0; it_cols != S[0].size() && accept; it_cols++) {

int n_zeros_col = 0;

for (std::vector<boost::dynamic_bitset<>>::iterator it_rows = S.begin(); it_rows != S.end(); it_rows++) {

n_zeros_col += !(*it_rows).test(it_cols);

}

accept = (n_zeros_col >= D);

}

if (!accept)

return false;

std::vector<boost::dynamic_bitset<>>::iterator it1, it2, it3; // iterators to browse rows of matrix

std::vector<std::vector<int>> neighbors; // Collect the indices of neighbors of each row

bool notadjacent;

// Initialize the neighbors data structure with S.size() empty lists.

std::vector<int> empty_neighborhood;

for (it1=S.begin(); it1!=S.end(); ++it1)

neighbors.push_back(empty_neighborhood);

// Dual graph computation

// std::cout << "Dual graph computation..." << std::endl;

// This loop computes the dual graph of S by going through all unordered pairs

// of rows of S (indexed by it1 and it2), and checking whether the union does

// not contain a third row (indexed by it3)

for (it1=S.begin(); it1!=S.end(); ++it1) {

int i1 = (int)distance(S.begin(),it1);

// Compute the neighbor of first row, with bigger index

//std::cout << "New neighbors of " << distance(S.begin(),it1) << ": ";

for (it2=it1+1; it2!=S.end(); ++it2) {

int i2 = (int)distance(S.begin(),it2);

// Compute the union of the two rows

boost::dynamic_bitset<> u = (*it1) | (*it2);

// Go through all the other rows and check whether union contains some other row.

// If such a row is found, that means that the two rows are not adjacent.

notadjacent = false;

for (it3=S.begin(); it3!=S.end() && !notadjacent; ++it3) {

if (it3 != it1 && it3 != it2 && (*it3).is_subset_of(u))

notadjacent = true;

}

// If an adjacency is detected, add the neighbor list of both rows

if (!notadjacent) {

//std::cout << i2 << ' ';

neighbors[i1].push_back(i2);

neighbors[i2].push_back(i1);

}

}

//std::cout << std::endl;

}

//std::cout << std::endl;

// Computation and verification of submatrices

//std::cout << "Computation of submatrices..." << std::endl;

accept = true;

for (std::vector<std::vector<int>>::iterator it=neighbors.begin(); it!=neighbors.end() && accept; it++)

accept = ((*it).size() >= D);

if (!accept)

return false;

bool found = true;

// Once the dual graph is known, go through all the rows and build the corresponding

// submatrix for each of them. If the input is a slack matrix, this will compute the

// slack matrix of the corresponding facet

for (it1=S.begin(); it1!=S.end() && found; ++it1) {

std::vector<boost::dynamic_bitset<>> SS; // submatrix of S corresponding to the row

int i1 = (int)distance(S.begin(),it1);

std::vector<size_type> vert_idx; // column indices in S of zeroes

size_type num_zeroes = 0; // number of columns

// Compute submatrix for current row

//std::cout << "Submatrix for row #" << i1;

// Count the zeroes in the row and record their positions

for (size_type i = 0; i < (*it1).size(); ++i) {

// If a zero is found, record its position

if (!(*it1).test(i)) {

vert_idx.push_back(i);

++num_zeroes;

}

}

// std::cout << " (" << (*it1) << " - " << "#zeroes = " << num_zeroes << ")" << std::endl;

// Go through all neighbors of the row

for (std::vector<int>::iterator it4=neighbors[i1].begin(); it4!=neighbors[i1].end(); ++it4) {

boost::dynamic_bitset<> new_row(num_zeroes);

// Go through all vertices of the row, and collect the bits from the slack matrix

// to form the slack matrix for the row (which is a submatrix of S)

for (size_type i = 0; i < num_zeroes; ++i) {

size_type j = vert_idx[i];

new_row[i] = (S[*it4])[j];

}

// Add new_row to submatrix

SS.push_back(new_row);

//std::cout << new_row << std::endl;

}

// Turn the submatrix into a nonincidence graph, and canonicize it

std::vector<setword> canonical = canonicize(SS);

found = false;

// Obtain range with atoms_cg corresponding to nonincidence graphs with same hash

std::pair <std::multimap<int,std::vector<setword>>::iterator, std::multimap<int,std::vector<setword>>::iterator> r;

r = atoms_cg.equal_range(hash_matrix(SS,D-1));

// Browse through all nonincidence graphs that have the same hash to see if one of them

// is isomorphic to the current nonincidence graph

for (std::multimap<int,std::vector<setword>>::iterator it5=r.first; it5!=r.second && !found; ++it5)

found = ((*it5).second == canonical);

//if (found) {

// std::cout << "OK: nonincidence graph isomorphic to that of an atom" << std::endl;

//}

//else {

// std::cout << "ERROR: nonincidence graph not isomorphic to that of an atom" << std::endl;

//}

//

//std::cout << "-" << std::endl;

}

return found;

}

bool is_equal(ublas::vector<int> A,ublas::vector<int> B) {

assert(A.size() == B.size());

bool flag = true;

for (size_type it = 0; it != A.size() && flag; it++)

flag = (A[it] == B[it]);

return flag;

}

bool is_susp(std::vector<boost::dynamic_bitset<>> S_new) {

// For all rows i of the slack matrix M

// Partition the columns into F_0 = {j : M(i,j) = 0} and F_1 = {j : M(i,j) = 1}

// For all translation vectors t such that the vertices of F_1 - t are a subset of those of F_0

// Check that F_1 - t is a face by testing whether the intersection of all facets of P that contain it is exactly F_1 - t

bool flag = false;

std::vector<ublas::vector<int>> S_cols; // slack embedding of the polytope having S_new as slack matrix

for (size_type col_idx = 0; col_idx != S_new[0].size(); col_idx++) {

ublas::vector<int> column(S_new.size());

for (int it0 = 0; it0 != S_new.size(); ++it0)

column(it0) = S_new[it0][col_idx];

S_cols.push_back(column);

}

for (std::vector<boost::dynamic_bitset<>>::iterator it0 = S_new.begin(); it0 != S_new.end(); ++it0) {

std::vector<size_type> zeros_idx; // column indices in S of zeroes

std::vector<size_type> ones_idx; // column indices in S of ones

// Count zeroes and ones in the row and record their positions

for (size_type i = 0; i < (*it0).size(); ++i) {

if ((*it0).test(i))

ones_idx.push_back(i);

else

zeros_idx.push_back(i);

}

for (std::vector<size_type>::iterator it1 = ones_idx.begin();it1 != ones_idx.end(); it1++) {

for (std::vector<size_type>::iterator it2 = zeros_idx.begin();it2 != zeros_idx.end(); it2++) {

ublas::vector<int> translation_vect = S_cols[*it2] - S_cols[*it1];

bool is_contained = true;

std::vector<size_type> idx_translated_F1;

for (std::vector<size_type>::iterator it3 = ones_idx.begin();it3 != ones_idx.end() && is_contained; it3++) {

ublas::vector<int> translated_F1_point = S_cols[*it3] + translation_vect;

bool is_found = false;

for (int h = 0; h != zeros_idx.size() && !is_found; h++) {

is_found = (is_equal(translated_F1_point,S_cols[zeros_idx[h]]));

if (is_found)

idx_translated_F1.push_back(zeros_idx[h]);

}

is_contained &= is_found;

}

if (is_contained) {

boost::dynamic_bitset<> char_F1(S_new[0].size());

char_F1.set();

for (int h = 0; h != idx_translated_F1.size(); h++)

char_F1.flip(idx_translated_F1[h]);

boost::dynamic_bitset<> intersect_rows_containing_F1(S_new[0].size());

for (std::vector<boost::dynamic_bitset<>>::iterator it0 = S_new.begin(); it0 != S_new.end(); ++it0) {

if ((*it0).is_subset_of(char_F1))

intersect_rows_containing_F1 |= *it0;

}

flag = (intersect_rows_containing_F1 == char_F1);

}

}

}

}

return flag;

}

int is_polar(std::vector<boost::dynamic_bitset<>> S_new,int &D) {

std::vector<boost::dynamic_bitset<>> S_transpose;

for (size_type col_idx = 0; col_idx != S_new[0].size(); col_idx++) {

boost::dynamic_bitset<> column;

for (size_type row_idx = 0; row_idx != S_new.size(); row_idx++)

column.push_back(S_new[row_idx][col_idx]);

S_transpose.push_back(column);

}

int hash_S = hash_matrix(S_new,D);

int hash_S_T = hash_matrix(S_transpose,D);

std::vector<setword> canonical_S_T = canonicize(S_transpose);

int amount_polar = 0;

if (hash_S == hash_S_T) {

if (canonicize(S_new) == canonical_S_T)

amount_polar = 1; // self-polar

}

else {

bool is_isomorphic = false;

std::pair <std::multimap<int,std::vector<setword>>::iterator, std::multimap<int,std::vector<setword>>::iterator> r;

r = LD.equal_range(hash_S_T);

// Browse through all nonincidence graphs that have the same hash to see if one of them

// is isomorphic to the current nonincidence graph

for (std::multimap<int,std::vector<setword>>::iterator it1=r.first; it1!=r.second && !is_isomorphic; ++it1)

is_isomorphic = ((*it1).second == canonical_S_T);

if (is_isomorphic)

amount_polar = 2;

}

return amount_polar;

}

void to_list(std::vector<boost::dynamic_bitset<>> S_new, int &D, const int &verbose){

int hash_S = hash_matrix(S_new,D);

std::vector<setword> canonical_S = canonicize(S_new);

bool is_isomorphic = false;

// Obtain range with atoms_cg corresponding to nonincidence graphs with same hash

std::pair <std::multimap<int,std::vector<setword>>::iterator, std::multimap<int,std::vector<setword>>::iterator> r;

r = LD.equal_range(hash_S);

// Browse through all nonincidence graphs that have the same hash to see if one of them

// is isomorphic to the current nonincidence graph

for (std::multimap<int,std::vector<setword>>::iterator it1=r.first; it1!=r.second && !is_isomorphic; ++it1)

is_isomorphic = ((*it1).second == canonical_S);

if (!is_isomorphic) {

LD.insert(std::pair<int,std::vector<setword>>(hash_S,canonical_S));

for (std::vector<boost::dynamic_bitset<>>::iterator it2 = S_new.begin(); it2 != S_new.end(); ++it2)

myfile << *it2 << std::endl;

myfile << '-' << std::endl;

if (verbose != 0) {

// check is there exists a simplicial facet

bool has_simplicial = false;

for (std::vector<boost::dynamic_bitset<>>::iterator row = S_new.begin(); row != S_new.end() && !has_simplicial; row++)

has_simplicial = ((~*row).count() == D);

if (has_simplicial)

simplicial_facet++;

// check if there exists a simple vertex

bool STAB = false;

for (size_type col_idx = 0; col_idx != S_new[0].size() && !STAB; col_idx++) {

boost::dynamic_bitset<> column;

for (int it0 = 0; it0 != S_new.size(); ++it0)

column.push_back(S_new[it0][col_idx]);

STAB = ((~column).count() == D);

}

if (STAB)

stab++;

// check if the polytope is centrally symmetric

bool CS = true;

for (std::vector<boost::dynamic_bitset<>>::iterator row = S_new.begin(); row != S_new.end() && CS; row++)

CS = ((*row).count() == S_new[0].size()/2);

if (CS)

cs++;

// tests if the polytope is a suspension

if (is_susp(S_new))

n_suspensions++;

// tests if the polytope has a polar that is a polytope

n_polar+=is_polar(S_new,D);

}

}

}

int main(int argc, const char *argv[]) {

if (argc != 3) {

std::cout << "\nERROR: Wrong input, not enough arguments" << std::endl;

std::cout << "Please insert two arguments: dimension D, verbose_flag." << std::endl;

std::cout << "- D is an integer between 3 and 7." << std::endl;

std::cout << "- verbose_flag is an integer between 0 and 3." << std::endl;

std::cout << "verbose_flag = 0 : minimal output (the benchmark times correspond to this flag)" << std::endl;

std::cout << "verbose_flag = 1 : output includes all closed sets + tests for classes of 2L polytopes" << std::endl;

std::cout << "verbose_flag = 2 : output includes all closed sets + H-embedding of ground set + tests for classes of 2L polytopes" << std::endl;

std::cout << "verbose_flag = 3 : output includes all closed sets, ground sets, slabs_sat_by_point and points_sat_by_slab + tests for classes of 2L polytopes" << std::endl << std::endl;

std::cout << "The input should be of the form ./2L_enum D verbose_flag, e.g.:" << std::endl;

std::cout << "./2L_enum 3 3" << std::endl;

std::cout << "which corresponds to D = 3 and verbose_flag = 3." << std::endl << std::endl;

return 1;

}

int D = atoi(argv[1]);

int verbose = atoi(argv[2]);

// Warnings and bound on D and verbose

if ((D != 1) && (D != 2) &&(D != 3) && (D != 4) && (D != 5) && (D != 6) && (D != 7)) {

std::cout << "\nERROR: Input (dimension D) out of bounds." << std::endl;

std::cout << "Please insert an integer value between 1 and 7 as dimension D." << std::endl;

std::cout << "The input should be of the form: ./2L_enum D verbose_flag." << std::endl << std::endl;

return 1;

}

if ((verbose != 0) && (verbose != 1) && (verbose != 2) && (verbose != 3)) {

std::cout << "\nERROR: Input (verbose_flag) out of bounds." << std::endl;

std::cout << "Please insert an integer value between 0 and 3 as verbose_flag." << std::endl;

std::cout << "verbose_flag = 0 : minimal output (the benchmark times correspond to this flag)" << std::endl;

std::cout << "verbose_flag = 1 : output includes all closed sets + tests for classes of 2L polytopes" << std::endl;

std::cout << "verbose_flag = 2 : output includes all closed sets + H-embedding of ground set + tests for classes of 2L polytopes" << std::endl;

std::cout << "verbose_flag = 3 : output includes all closed sets, ground sets, slabs_sat_by_point and points_sat_by_slab + tests for classes of 2L polytopes" << std::endl << std::endl;

std::cout << "The input should be of the form: ./2L_enum D verbose_flag." << std::endl << std::endl;

return 1;

}

if (D == 7) {

std::cout << "\nWARNING: the computation time for D = 7 is expected to take up to ~61 hours of computation time (on AMD Opteron(TM) 6134 2.3 GHz)." << std::endl;

std::cout << "Please insert 'y' to confirm: ";

char ch = getchar();

if((ch != 'y') && (ch != 'Y')) {

ungetc(ch, stdin);

return 1;

}

}

std::vector<std::vector<boost::dynamic_bitset<>>> atoms; // Vector with all the (D-1)-dim 2L-polytopes

std::vector<boost::dynamic_bitset<>> S; // A slack matrix of a (D-1)-dim polytope

boost::dynamic_bitset<> row; // Some row a slack matrix

std::string line;

// counters for number of 2-level polytopes

int total_2level = 0;

//int index = 0;

// counters for timers

std::chrono::time_point<std::chrono::system_clock> tot_start, tot_end;

std::chrono::time_point<std::chrono::system_clock> start_per_base, end_per_base;

std::chrono::time_point<std::chrono::system_clock> start_next_closure, end_next_closure;

std::chrono::time_point<std::chrono::system_clock> start_slack_matrix, end_slack_matrix;

std::chrono::time_point<std::chrono::system_clock> start_skip_test, end_skip_test;

std::chrono::time_point<std::chrono::system_clock> start_2level_test, end_2level_test;

double tot_next_closure = 0;

double tot_slack_matrix = 0;

double tot_2level_test = 0;

tot_start = std::chrono::system_clock::now();

std::cout << "\nWarning: in dynamic_bitsets, lowest index bits are written to the right!" << std::endl;

std::cout << "Reading all " << (D-1) << "-dimensional 2-level polytopes... ";

if (verbose != 0)

std::cout << std::endl;

// Open the file that contains (D-1)-dim 2L-polytopes

std::ifstream inputfile("./" + std::to_string(D-1) + "d.txt");

if (inputfile.is_open()) {

while (getline(inputfile,line)) {

// If we read some data

if (line.compare(std::string("-")) != 0) {

// Remove spaces

// string::iterator end_pos = std::remove(line.begin(),line.end(),' ');

// line.erase(end_pos, line.end());

// Turn the std::string into a set

row = (boost::dynamic_bitset<>)line;

if (verbose != 0)

std::cout << row << " (" << row.size() << ")" << std::endl;

S.push_back(row);

}

// Otherwise this is the end

else {

// Add the slack matrix to the list

atoms.push_back(S);

// Don't forget to clear S

S.clear();

if (verbose != 0)

std::cout << "-" << std::endl;

}

}

inputfile.close();

std::cout << "OK" << std::endl;

}

else {

std::cout << "Unable to open file." << std::endl;

return 1;

}

std::cout << "Number of polytopes read = " << atoms.size() << std::endl;

std::cout << "Computing canonical forms for all nonincidence graphs... ";

// Loop trough all the atoms, compute canonical form and store it

for (std::vector<std::vector<boost::dynamic_bitset<>>>::iterator it1=atoms.begin(); it1!=atoms.end(); ++it1)

atoms_cg.insert(std::pair<int,std::vector<setword>>(hash_matrix(*it1,D-1),canonicize(*it1)));

std::cout << "OK" << std::endl;

std::cout << "Processing bases..." << std::endl;

myfile.open(std::to_string(D) + "d.txt");

int N_closed_sets = 0;

// Main loop: loop through all the bases

for (std::vector<std::vector<boost::dynamic_bitset<>>>::iterator it1=atoms.begin(); it1!=atoms.end(); ++it1) {

start_per_base = std::chrono::system_clock::now();

std::vector<boost::dynamic_bitset<>> S = *it1;

int n_base = (int)distance(atoms.begin(),it1) +1;

std::cout << "\nBase #" << n_base << ':' << std::endl;

// Check if there is a simplicial core on the "top left" of the matrix

std::cout << "Simplicial core? ";

if (checksimplicialcore(S,D-1))

std::cout << "OK" << std::endl;

// Extract (extended) embedding transformation matrix M_d(0)

ublas::matrix<int> M = extractM(S,D);

if ((verbose==1) || (verbose==2) || (verbose == 3))

std::cout << "M_d(0) = " << M << std::endl;

// Compute inverse of matrix M_d(0)

ublas::matrix<int> Minv;

invertM(M,Minv);

if ((verbose==1) || (verbose==2) || (verbose == 3))

std::cout << "M_d(0)^{-1} = " << Minv << std::endl;

std::cout << "Constructing H-embedding of facets of the base... ";

// computing the facets of the base using the slack matrix S

std::vector<boost::dynamic_bitset<>> facets_base(S.size());

for (std::vector<boost::dynamic_bitset<>>::const_iterator it2 = S.begin(); it2 != S.end(); it2++) {

boost::dynamic_bitset<> E(D);

if ((*it2).test(D-1)) {

for (size_type i = 0; i < D-1; i++)

E[i+1] = !(*it2)[i];

}

else {

for (size_type i = 0; i < D-1; i++)

E[i+1] = (*it2)[i];

}

bool found = false;

for (std::vector<boost::dynamic_bitset<>>::iterator it3 = facets_base.begin(); it3 != facets_base.end() && !found; it3++)

found = ((*it3) == E);

if (!found) {

facets_base.push_back(E);

if (verbose == 3)

std::cout << E << ' ';

}

}

std::cout << "OK" << std::endl;

// Create the set Vert(P_0) (in V-embedding)

std::cout << "Building V-embedding of base... ";

std::vector<ublas::vector<int>> base_V;

// Loop through all vertices

for (int j=0; j<S[0].size(); j++) {

ublas::vector<int> point(D);

// Create a point whose first coordinate is 0, and the others are the D-1 first bits of the jth column of the slack matrix S

point[0] = 0;

for (int i=0; i<D-1; i++)

point[i+1] = (S[i])[j];

// Add point to the V-embedding of the ground set

base_V.push_back(point);

// Print point - for debugging

if (verbose==3)

std::cout << point << ' ';

}

std::cout << "OK" << std::endl;

std::cout << "Building V-embedding of the ground set... ";

std::vector<ublas::vector<int>> ground_set_V;

ublas::vector<int> count(D);

// initialize count to 0

for (int i=0; i<D; i++)

count(i) = 0;

bool carry;

// REDUCED GROUND SET

// Create the set of points in {1} x {-1,0,1}^{D-1} lex greater than e_1

// Its size is 1 + (3^{D-1} - 1)/2

ublas::vector<int> point(D);

point(0) = 1;

for (int i=1; i<D; i++)

point(i) = 0;

ground_set_V.push_back(point);

if (verbose == 3)

std::cout << point << ' ';

for (int i=D-1; i>0 ; i--) {

point(i) = 1;

count.clear();

while (count(D-1) == 0) {

int j;

// Extract a vector in {-1,0,1}^{D-i-1} to fill the vector

for (j=i+1; j<D; j++)

point(j) = count(j-1) - 1;

ground_set_V.push_back(point);

if (verbose == 2)

std::cout << point << ' ';

// Increase counter, by performing mod-3 computation

j = i;

do {

carry = (count(j) == 2);

count(j) = (count(j) + 1) % 3;

j++;

} while (carry && j < D);

}

}

std::cout << "OK" << std::endl;

// Create Vert(P_0) (H-embedding this time), the set of fixed points

std::cout << "Building H-embedding of base... ";

std::vector<ublas::vector<int>> base_H;

for (std::vector<ublas::vector<int>>::iterator it2=base_V.begin(); it2!=base_V.end(); it2++) {

ublas::vector<int> point(D);

point = prod(Minv,*it2);

base_H.push_back(point);

if (verbose == 3)