//============================================================================
// Name        : OSS_start.cpp
// Author      : Efrain H. Gonzalez
// Version     : 07/
// Compile	   : g++ -x c++ -std=c++11 -o OSSCode Downloads/OSS_start.cpp
// Copyright   : Your copyright notice
// Description : Hello World in C++, Ansi-style
//============================================================================

#include <iostream>
using namespace std;

#include <algorithm>
#include <boost/lexical_cast.hpp>
#include <boost/math/special_functions.hpp>
#include <boost/multiprecision/cpp_int.hpp>
#include <boost/random.hpp>
#include <cstdint>
#include <cstdlib>
#include <ctime>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <numeric>
#include <string>
#include <utility>
#include <vector>

using namespace std;
//Function for calculating the binomial coefficient
boost::multiprecision::cpp_int BinomialCoefficient(unsigned int n, unsigned int k) {
	 if (k == 0) { return 1; }
	else { return (n * BinomialCoefficient(n - 1, k - 1)) / k; }
}
//for sorting a vector of integer values
bool SortVecInt(int primero1, int segundo2) {
	return primero1 < segundo2;
}
//for sorting a vector of pairs based on the values in the vector of integers
bool SortVecsOnVecs(const std::pair < double ,vector <int> >& mefirst, const std::pair < double , vector <int> >& mefirst2) {
	bool testingvecs = false;

	size_t vecsize = mefirst.second.size();
	for(size_t frslst = 0; frslst < vecsize; ++frslst) {
		if(mefirst.second[frslst] < mefirst2.second[frslst]){
			testingvecs = true;
			break;
		}
		else if(mefirst.second[frslst] > mefirst2.second[frslst]){
			break;
		}
	}
	//std::cout << testingvecs << std::endl;
	//std::cout << std::endl;
	return testingvecs;
}
//Function for scoring an order given data
void OrderScore(vector < vector <int> > DATA, int totaltuples, int totalvariables,vector < int > nextorder, int variables, int maximumparents, vector < int > mincatval, vector <int> maxcatval,vector <int> numocats, vector <boost::multiprecision::cpp_int> FAMILY, vector< vector <double> >& VariableScores, double &havemeScore, vector < vector <string> >& VariableParentLabels) {
	int COLS=variables;
	//Your data here
	vector < vector <int> > DAT=DATA; // 2d array as a vector of vectors
	int totvars = totalvariables;
	int tottuples = totaltuples;

	//Ask User for the order they would like to test Order starts at 0 e.g. 0,1,2 is valid order for three variables
	vector <int> orden=nextorder;

	//Set the maximum for the amount of parents for any given variable
	unsigned int maxparents=maximumparents;

	//Lets convert to counts for every variable combination which would be 2^n in the case of binary variables starting with the minimum in each category:
	//categories in i
	vector <int> catsi=numocats;

	//Total Families Ui,alpha for a particular variable in the order
	vector <boost::multiprecision::cpp_int> families=FAMILY;


	//How many parent combinations for each step? As well as there counts
	vector< vector <int> > ParentCombos;
	vector< vector <int> > fullNijkvector;
	//Vector for first switching positions based on global index
	vector< vector <string> > switchvector;
	for (size_t i = 0; i < orden.size(); ++i) {
		//i represents the orden of the variable
		if (i == 0) {
			vector <int> tmp, Nijkovercombos1;
			tmp.push_back(1);
			ParentCombos.push_back(tmp);
			//counting the amount of times that a value of the first variable in the orden occurs
			//this starts with the maximum value for that variable
			for (int hello = maxcatval[orden[0]];hello >= mincatval[orden[0]];--hello) {
				//hello cycles through the categories of the first variable in the orden
				int Nijk1 = 0;
				//green cycles through tuples
				for (int green = 0; green < tottuples; ++green) {
					if (DAT[green][orden[0]] == hello) {
						Nijk1 += 1;
					}
				}
				Nijkovercombos1.push_back(Nijk1);
			}
			fullNijkvector.push_back(Nijkovercombos1);
		}
		else {
			vector <string> inswchvector;
			vector <int> tmp, Nijkovercombos1;
			tmp.push_back(1);
			int numparnts = i;
			//counting the amount of times that a value of the last variable in the current orden size occurs
			//this starts with the maximum value for that variable
			for (int hello = maxcatval[orden[numparnts]];hello >= mincatval[orden[numparnts]];--hello) {
				//hello cycles through the categories of the first variable in the orden
				int Nijk1 = 0;
				//green cycles through tuples
				for (int green = 0; green < tottuples; ++green) {
					if (DAT[green][orden[numparnts]] == hello) {
						Nijk1 += 1;
					}
				}
				Nijkovercombos1.push_back(Nijk1);
			}
			fullNijkvector.push_back(Nijkovercombos1);

			//j representing the number of parents
			for (int it = 1; it <= numparnts; ++it) {
				//(333)Creating a vector that uses the right combination
				boost::multiprecision::cpp_int Nloopy = 0;

				//Accounting for the limit of parent quantity
				if (it > maxparents) {
					break;
				}
				else {
					vector <int> NewMat(numparnts, 0);
					for (int p = 0; p < it; ++p) {
						NewMat[p] = 1;
					}

					//Nloopy represents the result of numparnts choose i e.g. numparnts choose 1 equals numparnts
					Nloopy = BinomialCoefficient(numparnts,it);
					for (int iNloopy = 0; iNloopy < Nloopy; ++iNloopy) {
						int combsparents = 1;

						//(444)This sets up the process for changing
						//PosOne tells me the position of the last one in the vector
						//We want to change when the position is the last position available in the vector
						int SumOnes = 0, PosOne = 0, SumOnes2 = 0, PosOne2, NxtOne = 0, FrstOne = 0;
						int SumOnes3 = 0, SumOnes4 = 0, SumY = 0;
						for (PosOne = (numparnts - 1); PosOne >= 0; --PosOne) {
							if (NewMat[PosOne] == 1) {
								break;
							}
						}
						for (int y = (numparnts - 1); y >= (numparnts - it); --y) {
							//SumOnes tells you the amount of ones in the last i columns
							//These are the last columns being considered
							SumOnes += NewMat[y];
						}
						for (PosOne2 = (numparnts - 1); PosOne2 >= 0; --PosOne2) {
							//SumOnes2 tells you the amount of ones before you reach the next zero
							//PosOne2 keeps track of the position of the coming zero
							SumOnes2 += NewMat[PosOne2];
							if ((SumOnes2 > 0) && (NewMat[PosOne2] == 0)) {
								break;
							}
						}
						for (FrstOne = 0; FrstOne < numparnts; ++FrstOne) {
							//FrstOne tells you the position of the first number 1 starting from the left hand side
							if (NewMat[FrstOne] == 1) {
								break;
							}
						}
						for (int x = (numparnts - 1); x >= (numparnts - it + 1); --x) {
							//SumOnes4 helps keep track of the sum of all ones located in the last i - 1 positions
							SumOnes4 += NewMat[x];
						}


						//Adding in the code that will allow counts parent combinations for this particular variable
						vector <int> iswvect;
						for (int q = 0; q < i; ++q) {
							if (NewMat[q] == 1) {
								combsparents *= catsi[orden[q]];
								iswvect.push_back(orden[q]);
							}
						}
						string stringOswvect;
						if(iswvect.size() > 1) {
							stringstream swvectstring;
							std::vector<int>::const_iterator findme;
							for(findme = iswvect.begin(); findme != iswvect.end(); ++findme) {
								if(findme != iswvect.begin()) {
									swvectstring << " ";
								}
								swvectstring << *findme;
							}
							stringOswvect = swvectstring.str();
						}
						else {
							stringOswvect = std::to_string(iswvect[0]);
						}
						inswchvector.push_back(stringOswvect);
						tmp.push_back(combsparents);

						vector <int> hvect;
						//hvect tells us which variables are being considered always the last variable is being considered
						//e.g if ABC is our orden and we are on i equals 1 then we are looking at relationships between A and B only
						//continued: A is the only one that is either a parent or isn't a parent so hvect will be < 0 1 >
						//for A C hvect will be < 0 2 >
						for (int h = 0; h < i; ++h) {
							if (NewMat[h] == 1) {
								hvect.push_back(h);
							}
						}
						hvect.push_back(numparnts);
						size_t shvect = hvect.size();
						//Prints out hvect
						/*for (int u = 0; u < shvect; ++u) {
							std::cout << hvect[u] << " ";
						}
						std::cout << endl;*/
						//Counting the amount of values in the data that have that particular parent combination
						vector <int> Nijkovercombos;
						for (int last = mincatval[orden[numparnts]]; last <= maxcatval[orden[numparnts]]; ++last) {
							//(333)Creating a vector that uses the right combination
							//std::cout << "This is for " << i << " place in the orden with value of variable equal to" << last << endl;
							vector <int> Test(shvect, last), maxtest;
							for (int p = 0; p < (shvect - 1); ++p) {
								Test[p] = maxcatval[orden[hvect[p]]];
							}
							maxtest = Test;
							for (int i2Nloopy = 0; i2Nloopy < combsparents; ++i2Nloopy) {

								//(444)This sets up the process for changing
								//NMpos tells me the position of the last non minimum value in the vector
								//We want to change when the position is the last position available in the vector
								int NMpos = 0, minpos = 0;
								for (NMpos = (shvect - 2); NMpos >= 0; --NMpos) {
									if (Test[NMpos] != mincatval[orden[hvect[NMpos]]]) {
										break;
									}
								}
								for (minpos = (shvect - 2); minpos >= 0; --minpos) {
									//minpos tells you the position of the last minimum value
									if (Test[minpos] == mincatval[orden[hvect[minpos]]]) {
										break;
									}
								}
								//Prints out Test
								/*for (int u = 0; u < shvect; ++u) {
									std::cout << Test[u] << " ";
								}
								std::cout << endl;
								std::cout << endl;
								std::cout << endl;*/
								//Count how many occurrences of the value are present in the data
								int Nijk = 0;
								for (int num2size = 0; num2size < tottuples; ++num2size) {
									int countcorrect = 0;
									for (size_t g = 0; g < Test.size(); ++g) {
										//num2size cycles through tuples
										//orden[hvect[g]] represents the variable in the orden that we are considering as a parent
										if (DAT[num2size][orden[hvect[g]]] == Test[g]) {
											countcorrect += 1;
										}
									}
									if (countcorrect == Test.size()) {
										Nijk += 1;
									}
								}
								//Nijkovercombos displays data as follows
								//it starts with the smallest value for the last variable in hvect
								//and the largest values in the first n-1 variables in hvect
								//max,max-1,max-2,max-3 e.g. 2, 1, 0, 2, 1, 0
								//count,count,count,count e.g. 13, 2, 2, 3, 4, 10
								Nijkovercombos.push_back(Nijk);
								//(666)Now that the values have been calculated find out what the next combination of variables should be
								if ((NMpos == -1) && (minpos == (shvect - 2))) {
									//break when the 1st non minimum does not exist and the first minimum is found in the last position e.g. 0000
									break;
								}
								if (minpos < NMpos) {
									Test[NMpos] = Test[NMpos] - 1;
								}
								else if (NMpos < minpos) {
									Test[NMpos] = Test[NMpos] - 1;
									for (int filler = NMpos + 1; filler < (shvect - 1); ++filler) {
										Test[filler] = maxtest[filler];
									}
								}
							}
						}
						fullNijkvector.push_back(Nijkovercombos);
						//(666)Now that the unique values have been calculated find out what the next combination of variables should be
						if ((PosOne == (numparnts - 1)) && (SumOnes == it)) {
							break;
						}
						else if ((PosOne == (numparnts - 1)) && (SumOnes != it)) {
							for (NxtOne = (numparnts - 1); NxtOne >= 0; --NxtOne) {
								//NxtOne tells you the position of the next closest number 1 that we would
								//like to change the position of (we will call it the important number one)
								//SumOnes3 helps keep track of the sum of all ones between now and the next important number one
								SumOnes3 += NewMat[NxtOne];
								if (SumOnes3 == (SumOnes2 + 1)) {
									break;
								}
							}
							if (SumOnes4 == (it - 1)) {
								//If all except one of the 1's are found in the last it - 1 columns
								for (int x = 0; x < numparnts; ++x) {
									if (((x <= (NxtOne + SumOnes3)) && (x > NxtOne)) || (x == (FrstOne + 1))) {
										//If
										NewMat[x] = 1;
									}
									else {
										NewMat[x] = 0;
									}
								}
							}
							else {
								for (int x = 0; x < numparnts; ++x) {
									if (((x <= (NxtOne + SumOnes3)) && (x > NxtOne)) || (x == FrstOne)) {
										//If the position is that of the first 1 or it falls between the changed number one and the total
										//amount of ones that are on that side of the zero 10111
										NewMat[x] = 1;
									}
									else if ((x != FrstOne) && (x != NxtOne) && (NewMat[x] == 1) && (x < PosOne2)) {
										//If it is not the position of the first 1 and it is not the position of the 1 whose position we are interested in changing
										//and the previous value at this position was 1 and the postion is below the value of the first zero spotted from the right
										NewMat[x] = 1;
									}
									else {
										NewMat[x] = 0;
									}
								}
							}
						}
						else if ((PosOne != (numparnts - 1)) && (SumOnes != it)) {
							for (NxtOne = (numparnts - 1); NxtOne >= 0; --NxtOne) {
								//NxtOne tells you the position of the next closest number 1 that we would
								//like to change the position of (we will call it the important number one)
								//SumOnes3 helps keep track of the sum of all ones between now and the next important number one
								SumOnes3 += NewMat[NxtOne];
								if (SumOnes3 == 1) {
									break;
								}
							}
							if (it != 1) {
								for (int x = 0; x < numparnts; ++x) {
									if (x == (NxtOne + 1)) {
										NewMat[x] = 1;
									}
									else if (x == NxtOne) {
										NewMat[x] = 0;
									}
									else if ((NewMat[x] == 1) && (x != NxtOne)) {
										NewMat[x] = 1;
									}
									else {
										NewMat[x] = 0;
									}
								}
							}
							else {
								for (int x = 0; x < numparnts; ++x) {
									if ((x == (NxtOne + 1))) {
										NewMat[x] = 1;
									}
									else {
										NewMat[x] = 0;
									}
								}
							}
						}
					}

				}
			}
			ParentCombos.push_back(tmp);
			switchvector.push_back(inswchvector);
		}
	}

	//Obtaining the actual score from this information
	//varinorder cycles through families (the amount of parent families that should be considered for the variable with a particular orden starting
	//the first variable in the orden)
	//keeping track of the position within the fullNijkvector associated with the varinorder and the qi_Uialpha
	int posinfull = 0;
	//finlogscore is the final score in natural log format
	double finlogscore = 0.0;
	vector< vector <double> > myvecvarparset;
	for (size_t varinorder = 0; varinorder < families.size(); ++varinorder) {
		//sumovUialpha is the the sum over all parent sets for a particular variable
		double sumovUialpha = 0.0;
		//vector of all values of seclastgamma
		vector <double> vec2ndlastgamma;
		double maxseclastgamma;
		//Uialpha cycles through all the parent sets for a particular family
		for (int Uialpha = 0; Uialpha < families[varinorder]; ++Uialpha) {
			// nijkprime represents the value of 1/(ri * qi)
			double nijkprime, nijprime;
			double rij = catsi[orden[varinorder]], PCs = ParentCombos[varinorder][Uialpha];

			nijprime = 1.0 / (PCs);
			nijkprime = 1.0 / (rij * PCs);
			//seclastgamma is the sum over all combinations for the parents in a set sum because it is logarithmic
			double seclastgamma = 0.0;
			//qi_Uialpha cycles through the combinations for the parents in a set
			for (int qi_Uialpha = 0; qi_Uialpha < ParentCombos[varinorder][Uialpha];++qi_Uialpha) {
				double lastgamma = 0.0;
				double nij = 0.0;
				//countijk cycles through the categories of the variable with a particular orden
				//catsi is in the orden that data is input and so one must use the orden[varinorder] to first obtain the variable that we are referring to
				//and then find the categories for it
				for (int countijk = 0; countijk < catsi[orden[varinorder]]; ++countijk) {
					double topy;
					//rightcol lets you find the right column/position of the value that you need for a particular category within the
					int rightcol = qi_Uialpha + (countijk * ParentCombos[varinorder][Uialpha]);
					nij += fullNijkvector[posinfull][rightcol];
					topy = (nijkprime + fullNijkvector[posinfull][rightcol]);

					//Using boost lgamma function for the product over categories and parent combinations
					lastgamma += boost::math::lgamma(topy) - boost::math::lgamma(nijkprime);

				}
				double boty = nij + nijprime;
				seclastgamma += lastgamma + boost::math::lgamma(nijprime) - boost::math::lgamma(boty);

			}
			vec2ndlastgamma.push_back(seclastgamma);
			//Calculate sumovUialpha based on the logsumexp concept
			if (Uialpha + 1 == families[varinorder]) {

				for (size_t que = 0; que < vec2ndlastgamma.size(); ++que) {
					//change the value of maxseclastgamma if new value is larger than the previous value
					if (que == 0) {
						maxseclastgamma = vec2ndlastgamma[0];
					}
					else {
						if (maxseclastgamma < vec2ndlastgamma[que]) {
							maxseclastgamma = vec2ndlastgamma[que];
						}
					}
				}
				for (size_t what = 0; what < vec2ndlastgamma.size(); ++what) {
					sumovUialpha += exp(vec2ndlastgamma[what] - maxseclastgamma);
				}
				//add info on parent set scores for each variable to this vector of vectors
				myvecvarparset.push_back(vec2ndlastgamma);

			}

			posinfull += 1;

		}
		finlogscore += log(sumovUialpha) + maxseclastgamma;
	}

	VariableScores = myvecvarparset;
	havemeScore = finlogscore;
	VariableParentLabels = switchvector;
}

//OSS Function for scoring using parent specifications
void OSS(vector < vector <int> > DATA, vector < vector < vector <int> > > ParentInfo,int totaltuples, int totalvariables,vector < int > nextorder, int variables, int maximumparents, vector < int > mincatval, vector <int> maxcatval,vector <int> numocats, vector< vector <double> >& VariableScores, double &havemeScore, vector < vector <string> >& VariableParentLabels) {
	int COLS=variables;
	//Your data here
	vector < vector <int> > DAT=DATA; // 2d array as a vector of vectors
	int totvars = totalvariables;
	int tottuples = totaltuples;

	//Ask User for the order they would like to test Order starts at 0 e.g. 0,1,2 is valid order for three variables
	vector <int> orden=nextorder;

	//Set the maximum for the amount of parents for any given variable
	unsigned int maxparents=maximumparents;

	//Lets convert to counts for every variable combination which would be 2^n in the case of binary variables starting with the minimum in each category:
	//categories in i
	vector <int> catsi=numocats;


	//How many parent combinations for each step? As well as there counts
	vector< vector <int> > ParentCombos;
	vector< vector <int> > fullNijkvector;

	//Vector for first switching positions based on global index
	vector< vector <string> > switchvector;

	//must use for loop to go through order
	for(int ordervariable = 0; ordervariable < orden.size(); ++ordervariable) {

		//use users input to determine the variables that should be used
		int globalvariable = orden[ordervariable];
		//These are the possible parents for the global variable based on its position in the current order which is orden
		vector <int> possibleparents;
		for (int pp1 = 0; pp1 < orden.size(); ++pp1) {
			if(orden[pp1]== globalvariable) {
				break;
			}
			else {
				possibleparents.push_back(orden[pp1]);
			}
		}
		// vector of combinations for this variable that were used to obtain counts
		vector <int> tmps;
		for(int checking = 0; checking < ParentInfo[globalvariable].size(); ++checking) {
			//The above focuses the search to a particular variables parent combinations
			int countpp = 0;
			for (int checking2 = 0; checking2 < ParentInfo[globalvariable][checking].size(); ++checking2){
				//The above focuses the search to a particular variable in the parent combination
				std::vector <int>::iterator myitpp;
				myitpp = std::find(possibleparents.begin(),possibleparents.end(),ParentInfo[globalvariable][checking][checking2]);
				if((myitpp != possibleparents.end()) || ((ParentInfo[globalvariable][checking][checking2] == globalvariable) && (ParentInfo[globalvariable][checking].size() == 1))) {
					//Check if the value is a possible parent and check if the value is the current variable that we are looking at
					// We are only interested in the current variable when its by itself in the file otherwise it does not make sense for it to
					// be its own parent as well as having another parent.
					// Ex: if the current variable is 4 (globalvariable = 4) then it would not make sense for the PARENTFILE to have a line that reads 1,4 which would mean that both 1 and 4 are
					// a parent for the variable 4.
					// But 4 by itself within the PARENTFILE would mean that 4 is its own parent and thus has no parents.
					countpp += 1;
				}
				else if (myitpp == possibleparents.end()){
					//erase the count because the variable cannot be found in the possible parents vector therefore we should not use this combination
					countpp = 0;
					break;
				}
			}
			int nucombsparents = 1;
			//We reorder so that the variables in the current ParentInfo[globalvariable][checking] vector is in the order of the global variables.
			vector <int> reconfigpar = ParentInfo[globalvariable][checking];
			std::sort(reconfigpar.begin(),reconfigpar.end());
			//added the current variable to the end of reconfigpar
			reconfigpar.push_back(globalvariable);
			vector <int> Nijkovercombosnu;



			//Check to make sure that all of the variables in the parent combination can be found in the order as possible parents of this variable
			// Also check to see if the only variable in the ParentInfo[globalvariable][checking] vector is the current variable
			if((countpp == ParentInfo[globalvariable][checking].size()) && (ParentInfo[globalvariable][checking][0] != globalvariable) ){
				//How many configurations are there for this parent combination
				for (int q1 = 0; q1 < countpp; ++q1) {
					nucombsparents *= catsi[reconfigpar[q1]];
				}
				tmps.push_back(nucombsparents);
				//now we must perform the counts for all configurations of these variables.
				for (int last1 = mincatval[globalvariable]; last1 <= maxcatval[globalvariable]; ++last1){
					//Start with the lowest possible category for the current variable and loop through all of the categories as you are ascending
					vector <int> Test1(countpp+1,last1), maxtest1, mintest1;
					for (int p1 = 0; p1 < countpp; ++p1) {
						//setting the first countpp values to the maximum possible categorical values for those parents
						Test1[p1] = maxcatval[reconfigpar[p1]];
					}
					maxtest1 = Test1;

					//Looping through all possible configurations
					for (int n2loop = 0; n2loop < nucombsparents; ++n2loop) {
						//(444)This sets up the process for changing
						//NMpos tells me the position of the last non minimum value in the vector
						//We want to change when the position is the last position available in the vector
						int NMposs = 0, minposs = 0;
						//countpp - 1 is so that I start from the right Ex: if reconfigpar is 0 1 then start from 1
						for (NMposs = (countpp - 1); NMposs >= 0; --NMposs) {
							if (Test1[NMposs] != mincatval[reconfigpar[NMposs]]) {
								break;
							}
						}

						for (minposs = (countpp - 1); minposs >= 0; --minposs) {
							//minpos tells you the position of the last minimum value
							if (Test1[minposs] == mincatval[reconfigpar[minposs]]) {
								break;
							}
						}

						//Count how many occurrences of the value are present in the data
						int Nijks = 0;
						for (int num2sizes = 0; num2sizes < tottuples; ++num2sizes) {
							int countcorrects = 0;
							for (size_t gs = 0; gs < Test1.size(); ++gs) {
								//num2size cycles through tuples
								//reconfigpar[gs] represents the variable in the order that we are considering as a parent
								if (DAT[num2sizes][reconfigpar[gs]] == Test1[gs]) {
									countcorrects += 1;
								}
							}
							if (countcorrects == Test1.size()) {
								Nijks += 1;
							}
						}
						//Nijkovercombos displays data as follows
						//it starts with the smallest value for the last variable in hvect
						//and the largest values in the first n-1 variables in hvect
						//max,max-1,max-2,max-3 e.g. 2, 1, 0, 2, 1, 0
						//count,count,count,count e.g. 13, 2, 2, 3, 4, 10
						Nijkovercombosnu.push_back(Nijks);
						//(666)Now that the values have been calculated find out what the next combination of variables should be
						if ((NMposs == -1) && (minposs == (countpp - 1))) {
							//break when the 1st non minimum does not exist and the first minimum is found in the last position e.g. 0000
							break;
						}
						if (minposs < NMposs) {
							Test1[NMposs] = Test1[NMposs] - 1;
						}
						else if (NMposs < minposs) {
							Test1[NMposs] = Test1[NMposs] - 1;
							for (int fillers = NMposs + 1; fillers < countpp ; ++fillers) {
								Test1[fillers] = maxtest1[fillers];
							}
						}
					}


				}
				fullNijkvector.push_back(Nijkovercombosnu);
			}
			else if ((countpp == ParentInfo[globalvariable][checking].size()) && (ParentInfo[globalvariable][checking][0] == globalvariable) ){
				//In this case the only parent was the current variable meaning that only counts for the current variable should be calculated
				tmps.push_back(1);
				//counting the amount of times that each category for the current variable occurs
				// starting at the maximum value for that variable
				for (int hellos = maxcatval[globalvariable]; hellos >= mincatval[globalvariable]; --hellos) {
					//hellos cycles through the categories of the variable
					int Nijk1s = 0;
					//greens cycles through tuples/rows
					for (int greens = 0; greens < tottuples; ++greens) {
						if (DAT[greens][globalvariable] == hellos) {
							Nijk1s += 1;
						}
					}
					Nijkovercombosnu.push_back(Nijk1s);
				}
				fullNijkvector.push_back(Nijkovercombosnu);
			}

		}
		//get counts for the combinations
		ParentCombos.push_back(tmps);
	}
	//Obtaining the actual score from this information
	int posinfulls = 0;
	//finlogscore is the final score in natural log format
	double finlogscores = 0.0;
	vector< vector <double> > myvecvarparsets;
	//varinorders cycles through all variables in the order
	for (size_t varinorders = 0; varinorders < orden.size(); ++varinorders){
		//sumovUialpha is the sum over all parent sets for a particular variable
		double sumovUialphas = 0.0;
		//vector of all values of seclastgamma
		vector <double> vec2ndlastgammas;
		double maxseclastgammas;
		//Uialpha cycles through all the parent sets for a particular family
		for (int Uialphas = 0; Uialphas < ParentCombos[varinorders].size(); ++Uialphas) {
			// nijkprime represents the value of 1/(ri * qi)
			double nijkprimes, nijprimes;
			double rijs = catsi[orden[varinorders]], PCs = ParentCombos[varinorders][Uialphas];

			nijprimes = 1.0 / (PCs);
			nijkprimes = 1.0 / (rijs * PCs);
			//seclastgamma is the sum over all combinations for the parents in a set sum because it is logarithmic
			double seclastgammas = 0.0;
			//qi_Uialpha cycles through the configurations for the parents in a set
			for (int qi_Uialphas = 0; qi_Uialphas < ParentCombos[varinorders][Uialphas];++qi_Uialphas) {
				double lastgammas = 0.0;
				double nijs = 0.0;
				//countijk cycles through the categories of the variable with a particular orden
				//catsi is in the orden that data is input and so one must use the orden[varinorder] to first obtain the variable that we are referring to
				//and then find the categories for it
				for (int countijks = 0; countijks < catsi[orden[varinorders]]; ++countijks) {
					double topys;
					//rightcol lets you find the right column/position of the value that you need for a particular category within the
					int rightcols = qi_Uialphas + (countijks * ParentCombos[varinorders][Uialphas]);
					nijs += fullNijkvector[posinfulls][rightcols];
					topys = (nijkprimes + fullNijkvector[posinfulls][rightcols]);

					//Using boost lgamma function for the product over categories and parent combinations
					lastgammas += boost::math::lgamma(topys) - boost::math::lgamma(nijkprimes);
				}
				double boty = nijs + nijprimes;
				seclastgammas += lastgammas + boost::math::lgamma(nijprimes) - boost::math::lgamma(boty);
			}
			vec2ndlastgammas.push_back(seclastgammas);
			//Calculate sumovUialpha based on the logsumexp concept
			if (Uialphas + 1 == ParentCombos[varinorders].size()) {

				for (size_t ques = 0; ques < vec2ndlastgammas.size(); ++ques) {
					//change the value of maxseclastgamma if new value is larger than the previous value
					if (ques == 0) {
						maxseclastgammas = vec2ndlastgammas[0];
					}
					else {
						if (maxseclastgammas < vec2ndlastgammas[ques]) {
							maxseclastgammas = vec2ndlastgammas[ques];
						}
					}
				}
				for (size_t whats = 0; whats < vec2ndlastgammas.size(); ++whats) {
					sumovUialphas += exp(vec2ndlastgammas[whats] - maxseclastgammas);
				}
				//add info on parent set scores for each variable to this vector of vectors
				myvecvarparsets.push_back(vec2ndlastgammas);

			}

			posinfulls += 1;
		}
		finlogscores += log(sumovUialphas) + maxseclastgammas;
	}

	VariableScores = myvecvarparsets;
	havemeScore = finlogscores;
	VariableParentLabels = switchvector;
}

//Function for swapping two positions in the order
vector<int> swapping(vector<int> order) {
	//Randomly swaping two numbers in the order vector
	int size = order.size();
	vector<int> diff_order;
	if (size > 1) {
		int rand1;
		int rand2;
		do {
			rand1 = rand() % (size - 1);

			rand2 = rand() % (size - 1);

			if (rand1 != rand2) {
				std::swap(order[rand1], order[rand2]);
				diff_order = order;
			}

		} while (rand1 == rand2);

	}
	else {
		std::cout << "Error: size of order is one" << '\n';
	}
	return diff_order;
}
//Function for cutting the deck
vector <int> cutdeck(vector<int> order) {
	//Randomly select where to cut
	int size2 = order.size();
	vector <int> diff_order2;
	if (size2 > 1) {
		int randnum;
		do {
			randnum = rand() % (size2 - 1);
			if (randnum != (size2 - 1)) {
				std::rotate(order.begin(),order.begin()+randnum,order.end());
				diff_order2 = order;
			}
		} while ((randnum == (size2 - 1))||(randnum == 0));

	}
	else {
		std::cout << "Error: size of order is one" << '\n';
	}
	return diff_order2;
}
//New function for calculating scores when only two positions have been swapped
void Swappin_Score(vector < vector <int> > DATA, int totaltuples, vector < int > prevOrder, vector < int > currOrder, int maximumparents, vector < int > mincatval, vector < int > maxcatval, vector < int > numocats, vector < boost::multiprecision::cpp_int> FAMILY, vector < vector <double> >& VariableScores1, vector < vector <string> >& YeOlswvector) {

	//Your data here
	vector < vector <int> > DAT=DATA; // 2d array as a vector of vectors

	int tottuples = totaltuples;

	//Ask for the previous orden that was scored
	vector <int> PrevOrder=prevOrder;

	//Ask for the new orden that we want to score
	vector <int> CurrOrder = currOrder;

	//Set the maximum for the amount of parents for any given variable
	unsigned int maxparents=maximumparents;

	//Lets convert to counts for every variable combination which would be 2^n in the case of binary variables starting with the minimum in each category:
	//categories in i
	vector <int> catsi=numocats;

	//Total Families Ui,alpha for a particular variable in the orden
	vector <boost::multiprecision::cpp_int> families2=FAMILY;

	vector < vector <double> > OGvecparset = VariableScores1;

	vector < vector <string> > Oldswitchvector = YeOlswvector;
	//First find which two variables were swapped and what the position of these two variables was
	//Which two variables?
	vector <int> varchange;
	//What was their position?
	vector <int> PositionChange;
	//Search
	for(size_t i = 0; i < PrevOrder.size();++i){
		if(PrevOrder[i] != CurrOrder[i]){
			//The first position will represent the old value and the second will be the new value
			varchange.push_back(PrevOrder[i]);
			PositionChange.push_back(i);
		}
	}

	//Do we expect there to be any variables between these two swapped positions?
	if( ((PositionChange[0] + 1) != PositionChange[1]) && ((PositionChange[0] - 1) != PositionChange[1]) ) {
		//Which are the variables inbetween the swapped positions?
		//Which variables?
		vector <int> betweeners;
		//What was their positon?
		vector <int> PositionBetw;
		for(int j = (PositionChange[0] + 1); j < PositionChange[1];++j){
			betweeners.push_back(PrevOrder[j]);
			PositionBetw.push_back(j);
		}

		vector< vector <int> > NuNijkvector1;
		vector< vector <int> > NuRentCombos1;
		//Tracevec where the changes are occuring in the original vecparset by using traceme
		vector < vector <int> > nTracevec1;
		//Vector for first switching positions based on global index
		vector< vector <string> > switchvector1;
		for(size_t lines1 = 0; lines1 < PositionBetw.size(); ++lines1){
			//for(int genCounts = 0;genCounts < (2^(PositionBetw[lines]-1)); ++genCounts){
			int numparnts1 = PositionBetw[lines1];
			//std::cout << "This is the number of parents: " << numparnts <<std::endl;
			vector <int> tmp1;
			vector <int> nTracin1;
			vector <string> inswchvector1;
				for(int it1 = 1; it1 <= numparnts1; ++it1) {
					boost::multiprecision::cpp_int Nloopy1 = 0;//, NcolFactorial1 = 1, iFactorial1 = 1, NiFactorial1 = 1;

					if(it1 > maxparents) {
						break;
					}
					else {
						vector <int> NewRents1(numparnts1,0);
						for(int p1 = 0; p1 < it1; ++p1){
							NewRents1[p1] = 1;
						}
						//Nloopy represents the result of numparnts choose i e.g. numparnts choose 1 equals numparnts
						//Nloopy1 = NcolFactorial1 / (iFactorial1 * NiFactorial1);
						Nloopy1 = BinomialCoefficient(numparnts1, it1);
						//std::cout << "Nloopy equals: " << Nloopy << std::endl;

						//Position from first combination

						for (int iNloopy1 = 0; iNloopy1 < Nloopy1; ++iNloopy1) {
							int combsparents1 = 1;
							//(444)This sets up the process for changing
							//PosOne tells me the position of the last one in the vector
							//We want to change when the position is the last position available in the vector
							int SumOnes1 = 0, PosOne1 = 0, SumOnes21 = 0, PosOne21, NxtOne1 = 0, FrstOne1 = 0;
							int SumOnes31 = 0, SumOnes41 = 0;
							for (PosOne1 = (numparnts1 - 1); PosOne1 >= 0; --PosOne1) {
								if (NewRents1[PosOne1] == 1) {
									break;
								}
							}
							for (int y1 = (numparnts1 - 1); y1 >= (numparnts1 - it1); --y1) {
								//SumOnes tells you the amount of ones in the last i columns
								//These are the last columns being considered
								SumOnes1 += NewRents1[y1];
							}

							for (PosOne21 = (numparnts1 - 1); PosOne21 >= 0; --PosOne21) {
								//SumOnes2 tells you the amount of ones before you reach the next zero
								//PosOne2 keeps track of the position of the coming zero
								SumOnes21 += NewRents1[PosOne21];
								if ((SumOnes21 > 0) & (NewRents1[PosOne21] == 0)) {
									break;
								}
							}

							for (FrstOne1 = 0; FrstOne1 < numparnts1; ++FrstOne1) {
								//FrstOne tells you the position of the first number 1 starting from the left hand side
								if (NewRents1[FrstOne1] == 1) {
									break;
								}
							}

							for (int x1 = (numparnts1 - 1); x1 >= (numparnts1 - it1 + 1); --x1) {
								//SumOnes4 helps keep track of the sum of all ones located in the last i - 1 positions
								SumOnes41 += NewRents1[x1];
							}
							//Prints out NewRents
							/*std::cout << "This is the NewRents vector: " << std::endl;
							for (int u = 0; u < numparnts; ++u) {
								std::cout << NewRents[u] << " ";
							}
							std::cout << endl;*/

							if(NewRents1[PositionChange[0]] == 1) {
								//Adding in the code that will allow counts parent combinations for this particular variable
								vector <int> iswvect1;
								for (int q1 = 0; q1 < numparnts1; ++q1) {
									if (NewRents1[q1] == 1) {
										combsparents1 *= catsi[CurrOrder[q1]];
										iswvect1.push_back(CurrOrder[q1]);
									}
								}
								string stringOswvect1;
								if(iswvect1.size() > 1) {
									stringstream swvectstring1;
									std::vector<int>::const_iterator findme1;
									for(findme1 = iswvect1.begin(); findme1 != iswvect1.end(); ++findme1) {
										if(findme1 != iswvect1.begin()) {
												swvectstring1 << " ";
										}
										swvectstring1 << *findme1;
									}
									stringOswvect1 = swvectstring1.str();
								}
								else {
									stringOswvect1 = std::to_string(iswvect1[0]);
								}
								inswchvector1.push_back(stringOswvect1);
								//string version of variable that was swapped to lower order
								for(size_t mesearch = 0; mesearch < Oldswitchvector[numparnts1 - 1].size(); ++mesearch) {
									//convert from string to integer vector (new)
									vector <int> Swap_vec;
									int svecin;
									std::stringstream svein(Oldswitchvector[numparnts1 - 1][mesearch]);
									while(svein >> svecin){
										Swap_vec.push_back(svecin);
									}
									//compare two vectors so that we can add to new tracevec the right position
									// We are only interested in checking vectors of the same size
									if( Swap_vec.size() == iswvect1.size() ) {
										//first find out if the vector has the element that was changed in this case the old element
										std::vector<int>::iterator pos1Swap;
										pos1Swap = std::find(Swap_vec.begin(), Swap_vec.end(), varchange[0]);
										//Does the new vector have the new element?
										std::vector<int>::iterator posNuvar;
										posNuvar = std::find(iswvect1.begin(),iswvect1.end(),varchange[1]);
										if((pos1Swap != Swap_vec.end()) && (posNuvar != iswvect1.end())) {
											//Now check if the only difference between the two vectors is the old element only
											int countnosame = 0; //the amount that were not equal if 0 then all except the new variable are equal
											for(size_t nsearch = 0; nsearch < iswvect1.size(); ++nsearch) {
												if (iswvect1[nsearch] != varchange[1]){
													std::vector<int>::iterator posSearch;
													posSearch = std::find(Swap_vec.begin(),Swap_vec.end(),iswvect1[nsearch]);
													if(posSearch == Swap_vec.end()) {
														countnosame += 1;
														break;
													}
												}
											}
											if(countnosame == 0){
												//It is! Lets record the position and move on
												nTracin1.push_back(mesearch + 1);
												Oldswitchvector[numparnts1 - 1][mesearch] = stringOswvect1;
												break;
											}
										}
									}
									else if(Swap_vec.size() < iswvect1.size()){
										continue;
									}
									else {
										break;
									}

								}
								tmp1.push_back(combsparents1);
								vector <int> hvect1;
								//hvect tells us which variables are being considered always the last variable is being considered
								//e.g if ABC is our orden and we are on i equals 1 then we are looking at relationships between A and B only
								//continued: A is the only one that is either a parent or isn't a parent so hvect will be < 0 1 >
								//for A C hvect will be < 0 2 >
								for (int h1 = 0; h1 < numparnts1; ++h1) {
									if (NewRents1[h1] == 1) {
										hvect1.push_back(h1);
									}
								}
								hvect1.push_back(numparnts1);
								size_t shvect1 = hvect1.size();
								//Prints out hvect
								/*for (int u = 0; u < shvect; ++u) {
									std::cout << hvect[u] << " ";
								}
								std::cout << endl;*/
								//Counting the amount of values in the data that have that particular parent combination
								vector <int> Nijkovercombos1;
								for (int last1 = mincatval[CurrOrder[numparnts1]]; last1 <= maxcatval[CurrOrder[numparnts1]]; ++last1) {
									//(333)Creating a vector that uses the right combination

									vector <int> Test1(shvect1, last1), maxtest1;
									for (int p1 = 0; p1 < (shvect1 - 1); ++p1) {
										Test1[p1] = maxcatval[CurrOrder[hvect1[p1]]];
									}
									maxtest1 = Test1;
									for (int i2Nloopy1 = 0; i2Nloopy1 < combsparents1; ++i2Nloopy1) {

										//(444)This sets up the process for changing
										//NMpos tells me the position of the last non minimum value in the vector
										//We want to change when the position is the last position available in the vector
										int NMpos1 = 0, minpos1 = 0;
										for (NMpos1 = (shvect1 - 2); NMpos1 >= 0; --NMpos1) {
											if (Test1[NMpos1] != mincatval[CurrOrder[hvect1[NMpos1]]]) {
												break;
											}
										}
										for (minpos1 = (shvect1 - 2); minpos1 >= 0; --minpos1) {
											//minpos tells you the position of the last minimum value
											if (Test1[minpos1] == mincatval[CurrOrder[hvect1[minpos1]]]) {
												break;
											}
										}
										//Prints out Test
										/*for (int u = 0; u < shvect; ++u) {
											std::cout << Test[u] << " ";
										}
										std::cout << endl;
										std::cout << endl;
										std::cout << endl;*/
										//Count how many occurrences of the value are present in the data
										int Nijk1 = 0;
										for (int num2size1 = 0; num2size1 < tottuples; ++num2size1) {
											int countcorrect1 = 0;
											for (size_t g1 = 0; g1 < Test1.size(); ++g1) {
												//num2size cycles through tuples
												//CurrOrder[hvect[g]] represents the variable in the orden that we are considering as a parent
												if (DAT[num2size1][CurrOrder[hvect1[g1]]] == Test1[g1]) {
													countcorrect1 += 1;
												}
											}
											if (countcorrect1 == Test1.size()) {
												Nijk1 += 1;
											}
										}
										//Nijkovercombos displays data as follows
										//it starts with the smallest value for the last variable in hvect
										//and the largest values in the first n-1 variables in hvect
										//max,max-1,max-2,max-3 e.g. 2, 1, 0, 2, 1, 0
										//count,count,count,count e.g. 13, 2, 2, 3, 4, 10
										Nijkovercombos1.push_back(Nijk1);
										//(666)Now that the values have been calculated find out what the next combination of variables should be
										if ((NMpos1 == -1) & (minpos1 == (shvect1 - 2))) {
											//break when the 1st non minimum does not exist and the first minimum is found in the last position e.g. 0000
											break;
										}
										if (minpos1 < NMpos1) {
											Test1[NMpos1] = Test1[NMpos1] - 1;
										}
										else if (NMpos1 < minpos1) {
											Test1[NMpos1] = Test1[NMpos1] - 1;
											for (int filler1 = NMpos1 + 1; filler1 < (shvect1 - 1); ++filler1) {
												Test1[filler1] = maxtest1[filler1];
											}
										}
									}
								}
								NuNijkvector1.push_back(Nijkovercombos1);
							}
							//Add so that we can keep track of position that was changed in The Previous vector of vectors called nuvecvarparset
							//(666)Now that the unique values have been calculated find out what the next combination of variables should be
							if ((PosOne1 == (numparnts1 - 1)) & (SumOnes1 == it1)) {
								break;
							}
							else if ((PosOne1 == (numparnts1 - 1)) & (SumOnes1 != it1)) {
								for (NxtOne1 = (numparnts1 - 1); NxtOne1 >= 0; --NxtOne1) {
									//NxtOne tells you the position of the next closest number 1 that we would
									//like to change the position of (we will call it the important number one)
									//SumOnes3 helps keep track of the sum of all ones between now and the next important number one
									SumOnes31 += NewRents1[NxtOne1];
									if (SumOnes31 == (SumOnes21 + 1)) {
										break;
									}
								}
								if (SumOnes41 == (it1 - 1)) {
									//If all except one of the 1's are found in the last it - 1 columns
									for (int x1 = 0; x1 < numparnts1; ++x1) {
										if (((x1 <= (NxtOne1 + SumOnes31)) & (x1 > NxtOne1)) | (x1 == (FrstOne1 + 1))) {
											//If
											NewRents1[x1] = 1;
										}
										else {
											NewRents1[x1] = 0;
										}
									}
								}
								else {
									for (int x1 = 0; x1 < numparnts1; ++x1) {
										if (((x1 <= (NxtOne1 + SumOnes31)) & (x1 > NxtOne1)) | (x1 == FrstOne1)) {
											//If the position is that of the first 1 or it falls between the changed number one and the total
											//amount of ones that are on that side of the zero 10111
											NewRents1[x1] = 1;
										}
										else if ((x1 != FrstOne1) & (x1 != NxtOne1) & (NewRents1[x1] == 1) & (x1 < PosOne21)) {
											//If it is not the position of the first 1 and it is not the position of the 1 whose position we are interested
											// in changing
										//and the previous value at this position was 1 and the postion is below the value of the first zero spotted from the
										// right
											NewRents1[x1] = 1;
										}
										else {
											NewRents1[x1] = 0;
										}
									}
								}
							}

							else if ((PosOne1 != (numparnts1 - 1)) & (SumOnes1 != it1)) {
								for (NxtOne1 = (numparnts1 - 1); NxtOne1 >= 0; --NxtOne1) {
									//NxtOne tells you the position of the next closest number 1 that we would
									//like to change the position of (we will call it the important number one)
									//SumOnes3 helps keep track of the sum of all ones between now and the next important number one
									SumOnes31 += NewRents1[NxtOne1];
									if (SumOnes31 == 1) {
										break;
									}
								}
								if (it1 != 1) {
									for (int x1 = 0; x1 < numparnts1; ++x1) {
										if (x1 == (NxtOne1 + 1)) {
											NewRents1[x1] = 1;
										}
										else if (x1 == NxtOne1) {
											NewRents1[x1] = 0;
										}
										else if ((NewRents1[x1] == 1) & (x1 != NxtOne1)) {
											NewRents1[x1] = 1;
										}
										else {
											NewRents1[x1] = 0;
										}
									}
								}
								else {
									for (int x1 = 0; x1 < numparnts1; ++x1) {
										if ((x1 == (NxtOne1 + 1))) {
											NewRents1[x1] = 1;
										}
										else {
											NewRents1[x1] = 0;
										}
									}
								}
							}
						}

					}

				}
				nTracevec1.push_back(nTracin1);
				switchvector1.push_back(inswchvector1);
				NuRentCombos1.push_back(tmp1);
		}
		//Obtaining the actual score from this information
		//varinorder cycles through families2 (the amount of parent families2 that should be considered for the variable with a particular orden starting
		//the first variable in the orden)
		//keeping track of the position within the NuNijkvector associated with the varinorder and the qi_Uialpha
		int posinfull1 = 0;

		vector< vector <double> > nuvecvarparset1;
		int icu1 = 0;
		for (int varinorder1 = PositionBetw[0]; varinorder1 < PositionChange[1]; ++varinorder1) {
			//vector of all values of seclastgamma
			vector <double> vec2ndlastgamma21;
			//Uialpha cycles through all the parent sets for a particular family
			for (size_t Uialpha1 = 0; Uialpha1 < NuRentCombos1[icu1].size(); ++Uialpha1) {
				// nijkprime represents the value of 1/(ri * qi)
				double nijkprime1, nijprime1;
				double rij1 = catsi[CurrOrder[varinorder1]], PCs1 = NuRentCombos1[icu1][Uialpha1];

				nijprime1 = 1.0 / (PCs1);
				nijkprime1 = 1.0 / (rij1 * PCs1);
				//seclastgamma is the sum over all combinations for the parents in a set sum because it is logarithmic
				double seclastgamma1 = 0.0;
				//qi_Uialpha cycles through the combinations for the parents in a set
				for (int qi_Uialpha1 = 0; qi_Uialpha1 < NuRentCombos1[icu1][Uialpha1];++qi_Uialpha1) {
					double lastgamma1 = 0.0;
					double nij1 = 0.0;
					//countijk cycles through the categories of the variable with a particular orden
					//catsi is in the orden that data is input and so one must use the orden[varinorder] to first obtain the variable that we are referring to
					//and then find the categories for it
					for (int countijk1 = 0; countijk1 < catsi[CurrOrder[varinorder1]]; ++countijk1) {
						double topy1;
						//rightcol lets you find the right column/position of the value that you need for a particular category within the
						int rightcol1 = qi_Uialpha1 + (countijk1 * NuRentCombos1[icu1][Uialpha1]);
						nij1 += NuNijkvector1[posinfull1][rightcol1];
						topy1 = (nijkprime1 + NuNijkvector1[posinfull1][rightcol1]);

						//Using boost lgamma function for the product over categories and parent combinations
						lastgamma1 += boost::math::lgamma(topy1) - boost::math::lgamma(nijkprime1);

					}
					double boty1 = nij1 + nijprime1;
					seclastgamma1 += lastgamma1 + boost::math::lgamma(nijprime1) - boost::math::lgamma(boty1);

				}
				vec2ndlastgamma21.push_back(seclastgamma1);
				//Calculate sumovUialpha based on the logsumexp concept
				if ((Uialpha1 + 1) == NuRentCombos1[icu1].size()) {
					//add infor on parent set scores for each variable to this vector of vectors
					nuvecvarparset1.push_back(vec2ndlastgamma21);

				}
				posinfull1 += 1;

			}
			icu1 += 1;
		}
		//Now make the changes in the original vector of Variable parent set scores
		for(size_t rep11 = 0; rep11 < nTracevec1.size(); ++rep11) {
			//rep1 will refer to the rows in Tracevec
			//each row starts at a different variable
			//the first row represents the first variable that was affected by the change so PositionBetw[0]
			for(size_t rep21 = 0; rep21 < nTracevec1[rep11].size(); ++rep21) {
				//rep2 will refer to the columns of Tracevec
				int makechangecol1 = nTracevec1[rep11][rep21];
				int makechangerow1 = PositionBetw[rep11];
				OGvecparset[makechangerow1][makechangecol1]=nuvecvarparset1[rep11][rep21];
			}
		}
	}
	//Second Part goes here:>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	//Store changes for the first and second variable
	vector <vector <int> > NuNijkvector2;
	vector <vector <int> > Tracevec2;
	vector < vector <int> > nTracevec2;
	vector <vector <int> > NuRentCombos2;
	vector< vector <string> > switchvector2;
	for(int lines2 = 0; lines2 < 1; ++lines2){
		//previous number of parents before swap for a variable
		int prevnumparnts = PositionChange[lines2];

		if(lines2 == 0) {
			//First variable that will now be placed in back

			//Number of parents for variable after swap
			int newnumparnts = PositionChange[1];
			//adding in the null parent case to the

			//Replace the values within the null cases first
			double VarUNO = OGvecparset[prevnumparnts][0];
			double VarDOS = OGvecparset[newnumparnts][0];
			OGvecparset[prevnumparnts][0] = VarDOS;
			OGvecparset[newnumparnts][0] = VarUNO;
			int posinOG2 = 1;
			vector <int> tmp2;
			vector <int> nTracin2;
			vector <string> inswchvector2;
			for(int it2 = 1; it2 <= newnumparnts; ++it2) {
				boost::multiprecision::cpp_int Nloopy2 = 0;

				if(it2 > maxparents) {
					break;
				}
				else {
					vector <int> NewRents2(newnumparnts,0);
					for(int p2 = 0; p2 < it2; ++p2){
						NewRents2[p2] = 1;
					}
					//Nloopy represents the result of newnumparnts choose i e.g. newnumparnts choose 1 equals numparnts
					Nloopy2 = BinomialCoefficient(newnumparnts,it2);

					//Position from first combination

					for (int iNloopy2 = 0; iNloopy2 < Nloopy2; ++iNloopy2) {
						int combsparents2 = 1;
						//(444)This sets up the process for changing
						//PosOne tells me the position of the last one in the vector
						//We want to change when the position is the last position available in the vector
						int SumOnes2 = 0, PosOne2 = 0, SumOnes22 = 0, PosOne22, NxtOne2 = 0, FrstOne2 = 0;
						int SumOnes32 = 0, SumOnes42 = 0, SumY2 = 0;
						for (PosOne2 = (newnumparnts - 1); PosOne2 >= 0; --PosOne2) {
							if (NewRents2[PosOne2] == 1) {
								break;
							}
						}
						for (int y2 = (newnumparnts - 1); y2 >= (newnumparnts - it2); --y2) {
							//SumOnes tells you the amount of ones in the last i columns
							//These are the last columns being considered

							SumOnes2 += NewRents2[y2];
						}
						for (PosOne22 = (newnumparnts - 1); PosOne22 >= 0; --PosOne22) {
							//SumOnes2 tells you the amount of ones before you reach the next zero
							//PosOne2 keeps track of the position of the coming zero
							SumOnes22 += NewRents2[PosOne22];
							if ((SumOnes22 > 0) & (NewRents2[PosOne22] == 0)) {
								break;
							}
						}

						for (FrstOne2 = 0; FrstOne2 < newnumparnts; ++FrstOne2) {
							//FrstOne tells you the position of the first number 1 starting from the left hand side
							if (NewRents2[FrstOne2] == 1) {
								break;
							}
						}

						for (int x2 = (newnumparnts - 1); x2 >= (newnumparnts - it2 + 1); --x2) {
							//SumOnes4 helps keep track of the sum of all ones located in the last i - 1 positions
							SumOnes42 += NewRents2[x2];
						}

						if(((PosOne2 + 1) > prevnumparnts)||(it2 > prevnumparnts)) {
							//Adding in the code that will allow counts parent combinations for this particular variable
							vector <int> iswvect2;
							for (int q2 = 0; q2 < newnumparnts; ++q2) {
								if (NewRents2[q2] == 1) {
									combsparents2 *= catsi[CurrOrder[q2]];
									iswvect2.push_back(CurrOrder[q2]);
								}
							}
							string stringOswvect2;
							if(iswvect2.size() > 1) {
								stringstream swvectstring2;
								std::vector<int>::const_iterator findme2;
								for(findme2 = iswvect2.begin(); findme2 != iswvect2.end(); ++findme2) {
									if(findme2 != iswvect2.begin()) {
										swvectstring2 << " ";
									}
									swvectstring2 << *findme2;
								}
								stringOswvect2 = swvectstring2.str();
							}
							else {
								stringOswvect2 = std::to_string(iswvect2[0]);
							}
							inswchvector2.push_back(stringOswvect2);
							for(size_t mesearch2 = 0; mesearch2 < Oldswitchvector[newnumparnts - 1].size(); ++mesearch2) {
								//convert from string to integer vector (new)
								vector <int> Swap_vec2;
								int svecin2;
								std::stringstream svein2(Oldswitchvector[newnumparnts - 1][mesearch2]);
								while(svein2 >> svecin2){
									Swap_vec2.push_back(svecin2);
								}
								//compare two vectors so that we can add to new tracevec the right position
								// We are only interested in checking vectors of the same size
								if( Swap_vec2.size() == iswvect2.size() ) {
									//Do they have the same elements?
									int countnosame2 = 0; //the amount that were not equal, if 0 then all are equal
									for(size_t nsearch2 = 0; nsearch2 < iswvect2.size(); ++nsearch2) {
										std::vector<int>::iterator posSearch2;
										posSearch2 = std::find(Swap_vec2.begin(),Swap_vec2.end(),iswvect2[nsearch2]);
										if(posSearch2 == Swap_vec2.end()) {
											countnosame2 += 1;
											break;
										}
									}
									if(countnosame2 == 0){
										//There is no difference! Lets record the position and move on
										nTracin2.push_back(mesearch2 + 1);
										Oldswitchvector[newnumparnts - 1][mesearch2] = stringOswvect2;
										break;
									}
									else {
										//first find out if the old vector has the element that was changed in this case the old element
										std::vector<int>::iterator posOlvar;
										posOlvar = std::find(Swap_vec2.begin(),Swap_vec2.end(),varchange[0]);
										//Does the new vector have the new element
										std::vector<int>::iterator posNuvar;
										posNuvar = std::find(iswvect2.begin(),iswvect2.end(),varchange[1]);
										if((posOlvar != Swap_vec2.end()) && (posNuvar != iswvect2.end())) {
											//Check if the only difference between the two is the old element
											int countnosame3 = 0; //the amount that were not equal, if 0 then all except new variable are equal
											for(size_t nsearch3 = 0; nsearch3 < iswvect2.size(); ++nsearch3) {
												if (iswvect2[nsearch3] != varchange[1]){
													std::vector<int>::iterator posSearch3;
													posSearch3 = std::find(Swap_vec2.begin(),Swap_vec2.end(),iswvect2[nsearch3]);
													if(posSearch3 == Swap_vec2.end()) {
														countnosame3 += 1;
														break;
													}
												}
											}
											if(countnosame3 == 0){
												//It is! Lets record the position and move on
												nTracin2.push_back(mesearch2 + 1);
												Oldswitchvector[newnumparnts - 1][mesearch2] = stringOswvect2;
												break;
											}
										}

									}
								}
								else if(Swap_vec2.size() < iswvect2.size()){
									continue;
								}
								else {
									break;
								}

							}
							tmp2.push_back(combsparents2);
							vector <int> hvect2;
							//hvect tells us which variables are being considered always the last variable is being considered
							//e.g if ABC is our orden and we are on i equals 1 then we are looking at relationships between A and B only
							//continued: A is the only one that is either a parent or isn't a parent so hvect will be < 0 1 >
							//for A C hvect will be < 0 2 >
							for (int h2 = 0; h2 < newnumparnts; ++h2) {
								if (NewRents2[h2] == 1) {
									hvect2.push_back(h2);
								}
							}
							hvect2.push_back(newnumparnts);
							size_t shvect2 = hvect2.size();
							//Prints out hvect
							/*for (int u = 0; u < shvect; ++u) {
								std::cout << hvect[u] << " ";
							}
							std::cout << endl;*/
							//Counting the amount of values in the data that have that particular parent combination
							vector <int> Nijkovercombos2;
							for (int last2 = mincatval[CurrOrder[newnumparnts]]; last2 <= maxcatval[CurrOrder[newnumparnts]]; ++last2) {
								//(333)Creating a vector that uses the right combination
								vector <int> Test2(shvect2, last2), maxtest2;
								for (int p2 = 0; p2 < (shvect2 - 1); ++p2) {
									Test2[p2] = maxcatval[CurrOrder[hvect2[p2]]];
								}
								maxtest2 = Test2;
								for (int i2Nloopy2 = 0; i2Nloopy2 < combsparents2; ++i2Nloopy2) {
									//(444)This sets up the process for changing
									//NMpos tells me the position of the last non minimum value in the vector
									//We want to change when the position is the last position available in the vector
									int NMpos2 = 0, minpos2 = 0;
									for (NMpos2 = (shvect2 - 2); NMpos2 >= 0; --NMpos2) {
										if (Test2[NMpos2] != mincatval[CurrOrder[hvect2[NMpos2]]]) {
											break;
										}
									}
									for (minpos2 = (shvect2 - 2); minpos2 >= 0; --minpos2) {
										//minpos tells you the position of the last minimum value
										if (Test2[minpos2] == mincatval[CurrOrder[hvect2[minpos2]]]) {
											break;
										}
									}
									//Prints out Test
									/*for (int u = 0; u < shvect; ++u) {
										std::cout << Test[u] << " ";
									}
									std::cout << endl;
									std::cout << endl;
									std::cout << endl;*/
									//Count how many occurrences of the value are present in the data
									int Nijk2 = 0;
									for (int num2size2 = 0; num2size2 < tottuples; ++num2size2) {
										int countcorrect2 = 0;
										for (size_t g2 = 0; g2 < Test2.size(); ++g2) {
											//num2size cycles through tuples
											//CurrOrder[hvect[g]] represents the variable in the orden that we are considering as a parent
											if (DAT[num2size2][CurrOrder[hvect2[g2]]] == Test2[g2]) {
												countcorrect2 += 1;
											}
										}
										if (countcorrect2 == Test2.size()) {
											Nijk2 += 1;
										}
									}
									//Nijkovercombos displays data as follows
									//it starts with the smallest value for the last variable in hvect
									//and the largest values in the first n-1 variables in hvect
									//max,max-1,max-2,max-3 e.g. 2, 1, 0, 2, 1, 0
									//count,count,count,count e.g. 13, 2, 2, 3, 4, 10
									Nijkovercombos2.push_back(Nijk2);
									//(666)Now that the values have been calculated find out what the next combination of variables should be
									if ((NMpos2 == -1) & (minpos2 == (shvect2 - 2))) {
										//break when the 1st non minimum does not exist and the first minimum is found in the last position e.g. 0000
										break;
									}
									if (minpos2 < NMpos2) {
										Test2[NMpos2] = Test2[NMpos2] - 1;
									}
									else if (NMpos2 < minpos2) {
										Test2[NMpos2] = Test2[NMpos2] - 1;
										for (size_t filler2 = NMpos2 + 1; filler2 < (shvect2 - 1); ++filler2) {
											Test2[filler2] = maxtest2[filler2];
										}
									}
								}
							}
							NuNijkvector2.push_back(Nijkovercombos2);
						}
						else {
							//using this to keep track of the values that have already been used
							vector <int> iswvect3;
							for (int q3 = 0; q3 < newnumparnts; ++q3) {
								if (NewRents2[q3] == 1) {
									iswvect3.push_back(CurrOrder[q3]);
								}
							}
							string stringOswvect3;
							if(iswvect3.size() > 1) {
								stringstream swvectstring3;
								std::vector<int>::const_iterator findme3;
								for(findme3 = iswvect3.begin(); findme3 != iswvect3.end(); ++findme3) {
									if(findme3 != iswvect3.begin()) {
										swvectstring3 << " ";
									}
									swvectstring3 << *findme3;
								}
								stringOswvect3 = swvectstring3.str();
							}
							else {
								stringOswvect3 = std::to_string(iswvect3[0]);
							}
							for(size_t mesearch3 = 0; mesearch3 < Oldswitchvector[newnumparnts - 1].size(); ++mesearch3) {
								//convert from string to integer vector (new)
								vector <int> Swap_vec4;
								int svecin4;
								std::stringstream svein4(Oldswitchvector[newnumparnts - 1][mesearch3]);
								while(svein4 >> svecin4){
									Swap_vec4.push_back(svecin4);
								}
								//compare two vectors so that we can add to new tracevec the right position
								// We are only interested in checking vectors of the same size
								if( Swap_vec4.size() == iswvect3.size() ) {
									//Do they have the same elements?
									int countnosame4 = 0; //the amount that were not equal, if 0 then all are equal
									for(size_t nsearch4 = 0; nsearch4 < iswvect3.size(); ++nsearch4) {
										std::vector<int>::iterator posSearch4;
										posSearch4 = std::find(Swap_vec4.begin(),Swap_vec4.end(),iswvect3[nsearch4]);
										if(posSearch4 == Swap_vec4.end()) {
											countnosame4 += 1;
											break;
										}
									}
									if(countnosame4 == 0){
										//They do have the same elements!
										size_t mesearch4;
										//Start by making two variables to empty the currently existing results into
										double FirstVar;
										double SecondVar = OGvecparset[newnumparnts][mesearch3+1];
										//Now lets find out where the corresponding value may be found for the FirstVar
										// first convert from string to vector of integers
										vector <int> olswap_vecOG;
										int ogvecin;
										std::stringstream ogvein(Oldswitchvector[prevnumparnts - 1][posinOG2 - 1]);
										while(ogvein >> ogvecin){
											olswap_vecOG.push_back(ogvecin);
										}
										int countnosame5 = 0; //the amount that were not equal, if 0 then all are equal
										for(size_t nsearch5 = 0; nsearch5 < iswvect3.size(); ++nsearch5) {
											std::vector<int>::iterator posSearch5;
											posSearch5 = std::find(olswap_vecOG.begin(),olswap_vecOG.end(),iswvect3[nsearch5]);
											if(posSearch5 == olswap_vecOG.end()) {
												countnosame5 += 1;
												break;
											}
										}
										if ((countnosame5 == 0) && (olswap_vecOG.size() == iswvect3.size())){
											//Then both vectors have the same elements and are of the same size
											FirstVar = OGvecparset[prevnumparnts][posinOG2];
											mesearch4 = posinOG2-1;
										}
										else {
											//in this case the length of both vectors was not the same or the posinOG2 did not contain the appropriate information so lets search for the right one
											for(mesearch4 = 0; mesearch4 < Oldswitchvector[prevnumparnts - 1].size(); ++mesearch4) {
												//check if the two vectors have the same size
												// first convert from string to vector of integers
												vector <int> olswap_vecOG1;
												int ogvecin1;
												std::stringstream ogvein1(Oldswitchvector[prevnumparnts - 1][mesearch4]);
												while(ogvein1 >> ogvecin1){
													olswap_vecOG1.push_back(ogvecin1);
												}
												if(olswap_vecOG1.size() == iswvect3.size()){
													//check if the two vectors are equal
													int countnosame6 = 0;
													for(size_t nsearch6 = 0; nsearch6 < iswvect3.size(); ++nsearch6) {
														std::vector<int>::iterator posSearch6;
														posSearch6 = std::find(olswap_vecOG1.begin(),olswap_vecOG1.end(),iswvect3[nsearch6]);
														if(posSearch6 == olswap_vecOG1.end()){
															countnosame6 += 1;
														}
													}
													if(countnosame6 == 0){
														break;
													}
												}
											}
											FirstVar = OGvecparset[prevnumparnts][mesearch4 + 1];
											posinOG2 = posinOG2 - 1;
										}
										//Now lets replace the original values with the values that should be there for the new orden
										OGvecparset[prevnumparnts][mesearch4 + 1] = SecondVar;
										OGvecparset[newnumparnts][mesearch3 + 1] = FirstVar;
										Oldswitchvector[newnumparnts-1][mesearch3] = stringOswvect3;
										break;
									}
								}
								else if(Swap_vec4.size() < iswvect3.size()){
									continue;
								}
								else {
									break;
								}
							}

							posinOG2+= 1;
						}
						//Add so that we can keep track of position that was changed in The Previous vector of vectors called nuvecvarparset2
						//(666)Now that the unique values have been calculated find out what the next combination of variables should be
						if ((PosOne2 == (newnumparnts - 1)) & (SumOnes2 == it2)) {
							break;
						}
						else if ((PosOne2 == (newnumparnts - 1)) & (SumOnes2 != it2)) {
							for (NxtOne2 = (newnumparnts - 1); NxtOne2 >= 0; --NxtOne2) {
								//NxtOne tells you the position of the next closest number 1 that we would
								//like to change the position of (we will call it the important number one)
								//SumOnes3 helps keep track of the sum of all ones between now and the next important number one
								SumOnes32 += NewRents2[NxtOne2];
								if (SumOnes32 == (SumOnes22 + 1)) {
									break;
								}
							}
							if (SumOnes42 == (it2 - 1)) {
								//If all except one of the 1's are found in the last it - 1 columns
								for (int x2 = 0; x2 < newnumparnts; ++x2) {
									if (((x2 <= (NxtOne2 + SumOnes32)) & (x2 > NxtOne2)) | (x2 == (FrstOne2 + 1))) {
										//If
										NewRents2[x2] = 1;
									}
									else {
										NewRents2[x2] = 0;
									}
								}
							}
							else {
								for (int x2 = 0; x2 < newnumparnts; ++x2) {
									if (((x2 <= (NxtOne2 + SumOnes32)) & (x2 > NxtOne2)) | (x2 == FrstOne2)) {
										//If the position is that of the first 1 or it falls between the changed number one and the total
										//amount of ones that are on that side of the zero 10111
										NewRents2[x2] = 1;
									}
									else if ((x2 != FrstOne2) & (x2 != NxtOne2) & (NewRents2[x2] == 1) & (x2 < PosOne22)) {
										//If it is not the position of the first 1 and it is not the position of the 1 whose position we are interested
										// in changing
									//and the previous value at this position was 1 and the postion is below the value of the first zero spotted from the
									// right
										NewRents2[x2] = 1;
									}
									else {
										NewRents2[x2] = 0;
									}
								}
							}
						}

						else if ((PosOne2 != (newnumparnts - 1)) & (SumOnes2 != it2)) {
							for (NxtOne2 = (newnumparnts - 1); NxtOne2 >= 0; --NxtOne2) {
								//NxtOne tells you the position of the next closest number 1 that we would
								//like to change the position of (we will call it the important number one)
								//SumOnes3 helps keep track of the sum of all ones between now and the next important number one
								SumOnes32 += NewRents2[NxtOne2];
								if (SumOnes32 == 1) {
									break;
								}
							}
							if (it2 != 1) {
								for (int x2 = 0; x2 < newnumparnts; ++x2) {
									if (x2 == (NxtOne2 + 1)) {
										NewRents2[x2] = 1;
									}
									else if (x2 == NxtOne2) {
										NewRents2[x2] = 0;
									}
									else if ((NewRents2[x2] == 1) & (x2 != NxtOne2)) {
										NewRents2[x2] = 1;
									}
									else {
										NewRents2[x2] = 0;
									}
								}
							}
							else {
								for (int x2 = 0; x2 < newnumparnts; ++x2) {
									if ((x2 == (NxtOne2 + 1))) {
										NewRents2[x2] = 1;
									}
									else {
										NewRents2[x2] = 0;
									}
								}
							}
						}
					}

				}

			}
			NuRentCombos2.push_back(tmp2);
			nTracevec2.push_back(nTracin2);
			switchvector2.push_back(inswchvector2);
		}
		//Obtaining the actual score from this information
		//varinorder cycles through families2 (the amount of parent families2 that should be considered for the variable with a particular orden starting
		//the first variable in the orden)
		//keeping track of the position within the NuNijkvector2 associated with the varinorder and the qi_Uialpha
		int posinfull2 = 0;
		//finlogscore2 is the final score in natural log format
		vector< vector <double> > nuvecvarparset2;
		int icu2 = 0;
		int varinorder2 = PositionChange[1];
		//vector of all values of seclastgamma
		vector <double> vec2ndlastgamma22;
		//Uialpha cycles through all the parent sets for a particular family
		for (size_t Uialpha2 = 0; Uialpha2 < NuRentCombos2[icu2].size(); ++Uialpha2) {
			// nijkprime represents the value of 1/(ri * qi)
			double nijkprime2, nijprime2;
			double rij2 = catsi[CurrOrder[varinorder2]], PCs2 = NuRentCombos2[icu2][Uialpha2];

			nijprime2 = 1.0 / (PCs2);
			nijkprime2 = 1.0 / (rij2 * PCs2);
			//seclastgamma is the sum over all combinations for the parents in a set sum because it is logarithmic
			double seclastgamma2 = 0.0;
			//qi_Uialpha cycles through the combinations for the parents in a set
			for (int qi_Uialpha2 = 0; qi_Uialpha2 < NuRentCombos2[icu2][Uialpha2];++qi_Uialpha2) {
				double lastgamma2 = 0.0;
				double nij2 = 0.0;
				//countijk cycles through the categories of the variable with a particular orden
				//catsi is in the orden that data is input and so one must use the orden[varinorder] to first obtain the variable that we are referring to
				//and then find the categories for it
				for (int countijk2 = 0; countijk2 < catsi[CurrOrder[varinorder2]]; ++countijk2) {
					double topy2;
					//rightcol lets you find the right column/position of the value that you need for a particular category within the
					int rightcol2 = qi_Uialpha2 + (countijk2 * NuRentCombos2[icu2][Uialpha2]);
					nij2 += NuNijkvector2[posinfull2][rightcol2];
					topy2 = (nijkprime2 + NuNijkvector2[posinfull2][rightcol2]);
						//Using boost lgamma function for the product over categories and parent combinations
					lastgamma2 += boost::math::lgamma(topy2) - boost::math::lgamma(nijkprime2);
					}
				double boty2 = nij2 + nijprime2;
				seclastgamma2 += lastgamma2 + boost::math::lgamma(nijprime2) - boost::math::lgamma(boty2);
			}
			vec2ndlastgamma22.push_back(seclastgamma2);
			//Calculate sumovUialpha based on the logsumexp concept
			if ((Uialpha2 + 1) == NuRentCombos2[icu2].size()) {
				//add infor on parent set scores for each variable to this vector of vectors
				nuvecvarparset2.push_back(vec2ndlastgamma22);

			}

			posinfull2 += 1;

		}
		icu2 += 1;

		//Now make the changes in the original vector of Variable parent set scores

		//we can make these replacements using the Tracevec2 and nuvecvarparset2

		//the first row represents the first and only row variable that was affected by the change so PositionChange[1]
		for(size_t rep22 = 0; rep22 < nTracevec2[0].size(); ++rep22) {
			//rep2 will refer to the columns of Tracevec2
			int makechangecol2 = nTracevec2[0][rep22];
			int makechangerow2 = PositionChange[1];
			OGvecparset[makechangerow2][makechangecol2]=nuvecvarparset2[0][rep22];


		}

	}

	//return
	VariableScores1 = OGvecparset;
	YeOlswvector = Oldswitchvector;

}
double meswapscore(vector < vector<double> > NEWVARIABLEScore){
	double finlogscore2 = 0.0;
	for(size_t finstep1 = 0; finstep1 < NEWVARIABLEScore.size(); ++finstep1) {
		double maxseclastgamma2;
		double SumforVarScore = 0.0;
		//first find the maximum for the variable
		for(size_t finstep2 = 0; finstep2 < NEWVARIABLEScore[finstep1].size(); ++finstep2) {
			if (finstep2 == 0) {
				maxseclastgamma2 = NEWVARIABLEScore[finstep1][0];
			}
			else {
				if (maxseclastgamma2 < NEWVARIABLEScore[finstep1][finstep2]) {
					maxseclastgamma2 = NEWVARIABLEScore[finstep1][finstep2];
				}
			}
		}
		//now lets apply the logsumexp method by using the above maximum
		for(size_t finstep3 = 0; finstep3 < NEWVARIABLEScore[finstep1].size(); ++finstep3) {
			SumforVarScore += exp(NEWVARIABLEScore[finstep1][finstep3]-maxseclastgamma2);
		}
		finlogscore2 += log(SumforVarScore) + maxseclastgamma2;
	}
	return finlogscore2;
}
//Sort through Order and Score
bool wayToSort(const std::pair < vector <int>, double > & uno, const std::pair < vector <int>, double >& dos) {
	return uno.second > dos.second;
}


//Do everything else in here
int main() {

	//control time
	//std::time_t meow = std::time(0);
	//boost::random::mt19937 mygen{ static_cast<std::uint32_t>(meow) };
	//boost::random::bernoulli_distribution<> dist;
	std::srand(std::time(NULL));
	//
	// Settings
	//
	string FILENAME; // Downloads/D50C9v.txt
	int COLS; // 9
	string PRIORFILE; // Documents/Materials_for_paper/Priorset1.txt
	string PARENTFILE; // home/efraingonzalez0/Downloads/OSSTESTfile.txt;
	vector <int> max_cat; // 1 2 1 3 2 1 2 2 3
	vector <int> min_cat; // 0 0 0 0 0 0 0 0 0
	float MAXTIME; // .5
	unsigned int maxparents; // 4
	double epsilon = 0.0; // .000001
	string StartOpt; // Y
	vector <int> StartingOrder; //2 6 8 4 1 0 3 5 7
	vector<string> settings;
	vector<bool> settingsFlag;

	settings.push_back("DataFile");
	settingsFlag.push_back(false);
	settings.push_back("TotalVariables");
	settingsFlag.push_back(false);
	settings.push_back("PriorFile");
	settingsFlag.push_back(false);
	settings.push_back("ParentFile");
	settingsFlag.push_back(false);
	settings.push_back("MaximumCategory");
	settingsFlag.push_back(false);
	settings.push_back("MinimumCategory");
	settingsFlag.push_back(false);
	settings.push_back("MaximumTime");
	settingsFlag.push_back(false);
	settings.push_back("MaximumParents");
	settingsFlag.push_back(false);
	settings.push_back("EpsilonDifference");
	settingsFlag.push_back(false);
	settings.push_back("StartingOption");
	settingsFlag.push_back(false);
	settings.push_back("StartingOrder");
	settingsFlag.push_back(false);

	fstream infile;
	string ConfigFile;
	std::cout << "What is the location of the configuration file that you edited?: ";
	std::cin >> ConfigFile;
	std::cout << std::endl;
	//string ConfigFile = "/home/efraingonzalez0/Downloads/new.config";
	infile.open(ConfigFile,ios::in);
	vector<string> line;
	if(infile.is_open()) {
		std::cout << "Configuration file correctly opened" << std::endl;
        	while(!infile.eof()) {
        		string tmp;
	        	getline(infile,tmp);
	        	line.push_back(tmp);
        	}
	}
	else {
		std::cout << "Unable to open configuration file" << std::endl;
		exit(0);
	}
	infile.close();
	line.erase(line.end());
	for(size_t k = 0; k < settings.size(); ++k) {
		for(size_t i = 0; i < line.size(); ++i) {
			size_t foundit = line[i].find(settings[k]);
			if (foundit != std::string::npos) {
				//find position of "="
				unsigned int poseq = line[i].find("=");
				unsigned int posneqfir,posneqlas;
				//first position of relevant string
				for(unsigned int j = poseq + 1; j < line[i].length(); ++j) {
					if(line[i].at(j) != ' ') {
						posneqfir = j;
						break;
					}
				}
				//last position of relevant string
				for(unsigned int j2 = line[i].length() - 1; j2 >= posneqfir; ++j2) {
					if((line[i].at(j2) != ' ') && (line[i].at(j2) != '\n')){
						posneqlas = j2;
						break;
					}
				}
				if(settings[k] == "DataFile") {
					FILENAME = line[i].substr(posneqfir,posneqlas - posneqfir + 1);
				}
				else if(settings[k] == "TotalVariables") {
						COLS = stoi(line[i].substr(posneqfir,posneqlas - posneqfir + 1));
				}
				else if(settings[k] == "PriorFile") {
					PRIORFILE = line[i].substr(posneqfir,posneqlas - posneqfir + 1);
				}
				else if (settings[k] == "ParentFile") {
					PARENTFILE = line[i].substr(posneqfir, posneqlas - posneqfir + 1);
				}
				else if(settings[k] == "MaximumCategory") {//still need to edit
					string max_cat_string = line[i].substr(posneqfir,posneqlas - posneqfir + 1);
					std::stringstream macs(max_cat_string);
					int macsti;
					while( macs >> macsti) {
						max_cat.push_back(macsti);
					}
				}
				else if(settings[k] == "MinimumCategory") {//still need to edit
					string min_cat_string = line[i].substr(posneqfir,posneqlas - posneqfir + 1);
					std::stringstream mics(min_cat_string);
					int micsti;
					while( mics >> micsti) {
						min_cat.push_back(micsti);
					}
				}
				else if(settings[k] == "MaximumTime") {
					MAXTIME = stof(line[i].substr(posneqfir,posneqlas - posneqfir + 1));
				}
				else if(settings[k] == "MaximumParents") {
					maxparents = stoi(line[i].substr(posneqfir,posneqlas - posneqfir + 1));
				}
				else if(settings[k] == "EpsilonDifference") {
					epsilon = stold(line[i].substr(posneqfir,posneqlas - posneqfir + 1));
				}
				else if(settings[k] == "StartingOption") {
					StartOpt = line[i].substr(posneqfir,posneqlas - posneqfir + 1);
				}
				else if(settings[k] == "StartingOrder") {//still need to edit
					string st_ord_string = line[i].substr(posneqfir,posneqlas - posneqfir + 1);
					std::stringstream starordr(st_ord_string);
					int StOrti;
					while( starordr >> StOrti) {
						StartingOrder.push_back(StOrti);
					}
				}
				settingsFlag[k] = true;
				break;
			}
		}
	}
	for(int k=0;k<settings.size();k++) {
		if(!settingsFlag[k]) {
			std::cout << "Error: "<< settings[k] << " is not configured." << std::endl;
			exit(0);
		}
		else if((StartingOrder.size() != COLS) && (StartOpt == "Y")) {
			std::cout << "Error: Number of TotalVariables and StartingOrder size do not match" << std::endl;
			exit(0);
		}
		else if((PARENTFILE == "NONE") || (PARENTFILE == " ")){
			std::cout << "Error: In order to use this function the ParentFile must exist." << std::endl;
			exit(0);
		}
		else if(max_cat.size() != COLS) {
			std::cout << "Error: Size of MaximumCategory and the number of TotalVariables do not match" << std::endl;
			exit(0);
		}
		else if(min_cat.size() != COLS) {
			std::cout << "Error: Size of MinimumCategory and the number of TotalVariables do not match" << std::endl;
			exit(0);
		}
		else if(max_cat.size() != min_cat.size()) {
			std::cout << "Error: Size of MaximumCategory and MinimumCategory do not match" << std::endl;
			exit(0);
		}
		else if((StartOpt != "Y") && (StartOpt != "n")) {
			std::cout << "Error: You must choose between the options yes represented by the capital letter Y and no represented by the letter n." << std::endl;
			exit(0);
		}
	}

	//Choose file to analyze
	//string FILENAME;
	std::cout << "What is the name of your file?: ";
	//std::cin >> FILENAME;
	std::cout << FILENAME << std::endl;
	std::cout << std::endl;
	// Variable declarations
	fstream file;
	//int COLS;
	std::cout << "How many variables in your data file?: ";
	//std::cin >> COLS;
	std::cout << COLS << std::endl;
	std::cout << endl;
	vector < vector <int> > DAT; // 2d array as a vector of vectors
	vector <int> rowVector(COLS); // vector to add into 'array' (represents a row)
	int row = 0; // Row counter

				 // Read file
	file.open(FILENAME, ios::in); // Open file
	if (file.is_open()) { // If file has correctly opened...
						  // Output debug message
		std::cout << "Data file correctly opened" << std::endl;

		// Dynamically store data into array
		while (file.good()) { // ... and while there are no errors,
			DAT.push_back(rowVector); // add a new row,
			for (int col = 0; col<COLS; col++) {
				file >> DAT[row][col]; // fill the row with col elements
			}
			row++; // Keep track of actual row
		}
	}
	else {
		std::cout << "Error: Unable to open data file" << std::endl;
		exit(0);
	}
	file.close();
	DAT.erase(DAT.end());
	int totvars = DAT[1].size();
	int tottuples = DAT.size();


	//read in prior information
	//string PRIORFILE;
	std::cout << "What is the name of the file containing prior information? (If none put NONE) ";
	//std::cin >> PRIORFILE;
	std::cout << PRIORFILE << std::endl;
	std::cout << std::endl;
	//Priors with Rest feature type 3 includes Dr. Yoo's specifications
	vector < std::pair < vector <string>, double > > priorset;
	if (PRIORFILE != "NONE") {
		fstream prfile;
		vector < vector <string> > PRIORinput; // 2d array as a vector of vectors
		vector <string> rowVector2(3); // vector to add into 'array' (represents a row)
		int row = 0; // Row counter

					 // Read file
		prfile.open(PRIORFILE, ios::in); // Open file
		if (prfile.is_open()) { // If file has correctly opened...
							  // Output debug message
			std::cout << "File correctly opened" << std::endl;

			// Dynamically store data into array
			while (prfile.good()) { // ... and while there are no errors,
				PRIORinput.push_back(rowVector2); // add a new row,
				for (int col = 0; col < 3; col++) {
					prfile >> PRIORinput[row][col]; // fill the row with col elements
				}
				row++; // Keep track of actual row
			}
		}
		else std::cout << "Unable to open file" << std::endl;
		prfile.close();
		//print PRIORinput
		for (size_t i = 0; i < PRIORinput.size();++i) {
			for (size_t j = 0; j < PRIORinput[i].size();++j) {
				std::cout << PRIORinput[i][j] << "   ";
			}
			std::cout << std::endl;
		}
		std::cout << std::endl;
		if (PRIORinput[PRIORinput.size()-1][0] == "") {
			PRIORinput.erase(PRIORinput.end());
		}
		//std::make_pair(orden, firstcore)

		//Extract Priorsets
		//vector < std::pair < vector <string>, string > > priorsettemp;
		for (size_t i = 0;i < PRIORinput.size();++i) {
			for (int j = 0; j < 2;++j) {
				int pos = 0;
				string delimiter = ",";
				string token;
				vector <string> frstset;
				while ((pos = PRIORinput[i][j].find(delimiter)) != std::string::npos) {
					token = PRIORinput[i][j].substr(0, pos);
					frstset.push_back(token);
					PRIORinput[i][j].erase(0, pos + delimiter.length());
				}
				frstset.push_back(PRIORinput[i][j]);
				if (j == 0) {
					priorset.push_back(std::make_pair(frstset, stod(PRIORinput[i][2])));
				}
				else {
					for (size_t checkfrst = 0; checkfrst < frstset.size(); ++checkfrst) {
						if ( frstset[checkfrst] == "REST") {
							vector <int> testvectorp;
							for (int blue = 0; blue < COLS ; ++blue) {
								testvectorp.push_back(blue);
							}
							int tvposize = testvectorp.size();
							int adjustinggallo = 0;
							for (size_t pollo = 0; pollo < priorset[(i*2)].first.size(); ++pollo) {

								for (int gallo = 0; gallo < tvposize; ++gallo) {
								//std::vector<string>::iterator itnew = gallo;
									if(gallo == std::stoi(priorset[(i*2)].first[pollo])) {
										testvectorp.erase(testvectorp.begin()+ gallo - adjustinggallo);
										adjustinggallo += 1;
									}
								}
							}
							//Adding the ability to remove certain variables from consideration in REST
							if ((frstset.size() != 1)) {

								for (size_t elpollo = 1; elpollo < frstset.size(); ++elpollo) {
									if (std::stoi(frstset[elpollo]) <= 0) {
										int tvposize2 = testvectorp.size();
										for (int elgallo = 0; elgallo < tvposize2; ++elgallo) {
											//std::vector<string>::iterator itnew = gallo;
											if(testvectorp[elgallo] == (-1 * std::stoi(frstset[elpollo]))) {
												testvectorp.erase(testvectorp.begin()+ elgallo);
												break;
											}
										}
									}
								}
							}

							frstset.clear();
							for (size_t red = 0; red < testvectorp.size(); ++red) {
								frstset.push_back(std::to_string(testvectorp[red]));
							}

						}

					}
					priorset.push_back(std::make_pair(frstset,1.0 - stod(PRIORinput[i][2])));
				}
			}
		}
		//print priorsets
		std::cout << "This is the priorset " << std::endl;
		for (size_t i = 0; i < priorset.size();++i) {
			for (size_t j = 0; j < priorset[i].first.size();++j) {
				std::cout << priorset[i].first[j] << " ";
			}
			std::cout << "     " << priorset[i].second;
			std::cout << std::endl;
		}
	}

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


	std::cout << "What is the maximum categorical value of each variable?: ";
	for (size_t i = 0; i < max_cat.size();++i) {
		std::cout << max_cat[i] << " ";
	}
	std::cout << endl;
	std::cout << endl;
	std::cout << "What is the minimum categorical value of each variable?: ";
	for (size_t i = 0; i < min_cat.size();++i) {
		std::cout << min_cat[i] << " ";
	}
	std::cout << std::endl;
	//Lets convert to counts for every variable combination which would be 2^n in the case of binary variables starting with the minimum in each category:
	//categories in i
	vector <int> catsi;
	//total combinations of variables
	for (int i = 0; i < totvars; ++i) {
		catsi.push_back((max_cat[i] - min_cat[i]) + 1);
	}
	//print out catsi
	std::cout << std::endl;
	std::cout << "Amount of Categories for each variable: ";
	for (size_t i = 0; i < catsi.size();++i) {
		std::cout << catsi[i] << " ";
	}
	std::cout << std::endl;
	std::cout << std::endl;


	//Set the name of your output file
	/*string FileOut;
	std::cout << "What is the name of your output file? (Include .txt at the end of the name): ";
	std::cin >> FileOut;*/

	//NOT TRACKING REPEATS
	//Set the amount of iterations to run through of the code
	/*int MAXITER;
	std::cout << "How many iterations would you like to run through?: ";
	std::cin >> MAXITER;*/
	std::cout << std::endl;
	//Set the maximum running time for the program in integer hours
	//float MAXTIME;
	std::cout << "How long would you like the program to run (in whole hours)?: ";
	//std::cin >> MAXTIME;
	std::cout << boost::lexical_cast<string>(MAXTIME) << std::endl;
	std::cout << std::endl;
	//Set the maximum for the amount of parents for any given variable
	//unsigned int maxparents;
	std::cout << "What is the maximum amount of parents to be considered for any given variable?: ";
	//std::cin >> maxparents;
	std::cout << maxparents << std::endl;
	std::cout << endl;
	//Set the epsilon difference value that you want to use
	//If the difference is calculated to be less than epsilon the program will use the cutdeck function
	//double epsilon;
	std::cout << "What is the epsilon difference value that you would like to use?: ";
	//std::cin >> epsilon;
	std::cout << boost::lexical_cast<string>(epsilon) << std::endl;
	std::cout << endl;
	//Give the User the option to input an order or have the starting order be random
	int ordere;
	vector <int> orden;
	std::cout << "Would you like to designate a starting order? [Y/n] " << StartOpt << std::endl;

	if (StartOpt == "Y") {
		//Starting order will be set
		/*std::cout << "Specify the starting order: ";
		for (int e = 0; e < COLS;++e) {
			std::cin >> ordere;
			orden.push_back(ordere);
		}*/
		std::cout << "The specified starting order: ";
		for (int e = 0; e < COLS;++e) {
			std::cout << StartingOrder[e] << " ";
		}
		std::cout << endl;
		std::cout << std::endl;
		orden = StartingOrder;
	}
	else {
		//Starting order will be random

		vector <int> simpord;
		for (int ruff = 0; ruff < COLS; ++ruff) {
			simpord.push_back(ruff);
		}
		std::random_shuffle(simpord.begin(), simpord.end());
		//Print out what the initial random order was
		std::cout << "The random starting order is: ";
		for (int f = 0; f < COLS;++f) {
			std::cout << simpord[f] << " ";
		}
		std::cout << std::endl;
		orden = simpord;
	}
	std::cout << std::endl;

	//read in parent information

	//std::cout << "What is the name of the file containing parent information? (If none put NONE) ";
	std::cout << std::endl;
	vector < vector < vector <int> > > parentset;
	if (PARENTFILE != "NONE") {
		fstream parentfile;
		parentfile.open(PARENTFILE,ios::in);
		vector < vector <string> > PARENTinput;
		if(parentfile.is_open()) {
			std::cout << "File with Parent combinations opened correctly"  << std::endl;
			while(!parentfile.eof()) {
				vector <string> tokens;
				string tmptemp;
				getline(parentfile,tmptemp);
				istringstream iss(tmptemp);
				std::string token;
				while(getline(iss,token,'\t')) {
					tokens.push_back(token);
				}
				PARENTinput.push_back(tokens);
			}
		}
		else {
			std::cout << "Unable to open the file with parent combinations" << std::endl;
			exit(0);
		}
		parentfile.close();
		PARENTinput.erase(PARENTinput.end());
		//by this point it is separated and the stuff in PARENTinput looks like
		//"1"	"2,0"	"1,0"
		//where the above means that 1 can be a parent the combination 2 and 0 can be a parent and the combination 1 and 0 can be a parent
		//First layer represent a variable
		//Second layer represents the parent combination
		//Third layer a value within the parent combination
		for (size_t k1 = 0; k1 < PARENTinput.size(); ++k1) {
			vector < vector <int> > mysetj;
			for (size_t j1 = 0; j1 < PARENTinput[k1].size(); ++j1) {
				vector <int> myseti;
				size_t pos1 = 0;
				string delimiter1 = ",";
				string token1;
				vector <string> nthset;

				while ((pos1 = PARENTinput[k1][j1].find(delimiter1)) != std::string::npos) {
					token1 = PARENTinput[k1][j1].substr(0, pos1);
					nthset.push_back(token1);
					PARENTinput[k1][j1].erase(0, pos1 + delimiter1.length() );
				}
				nthset.push_back(PARENTinput[k1][j1]);
				//Creating the three layers deep vector
				for (size_t nuk = 0; nuk < nthset.size(); ++nuk) {
					myseti.push_back(stoi(nthset[nuk]));
				}
				//Only two layers deep
				mysetj.push_back(myseti);
			}
			//Only one layer deep
			parentset.push_back(mysetj);
		}
	}
	std::cout << "This is the parentset." << std::endl;
	//print parentsets
	for (size_t i = 0; i < parentset.size(); ++i) {
		for (size_t j = 0; j < parentset[i].size();++j) {
			for (size_t k = 0; k < parentset[i][j].size(); ++k){
				std::cout << parentset[i][j][k] << " ";
			}
			std::cout << "       ";
		}
		std::cout << std::endl;
	}
	std::cout << std::endl;
	std::cout << std::endl;
	std::cout << std::endl;
	//This vector of vectors will hold all scores for each parent set of every variable
	vector < vector <double> > OGVarParScore;
	//This is the score obtained by the OrderScore method
	double OGScore;
	vector < vector <string> > OGLabels;
	OSS(DAT, parentset,tottuples, totvars,orden, COLS, maxparents, min_cat, max_cat,catsi, OGVarParScore, OGScore, OGLabels);
	std::cout << boost::lexical_cast<string>(OGScore) << std::endl;
	std::cout << std::endl;

	return 0;

}