/*Implementation of PStep_Algorithm::ColumnGeneration.

*/

#include "PStep_Algorithm.h"

ILOSTLBEGIN

//Column generation
//Solve the LP relaxation of the problem corresponding to AlgorithmMembers.CurrentTreeNode.
//New p-steps are added to AlgorithmMembers.PSteps, the solution is stored at
//AlgorithmMembers.CurrentSolution and statistics are stored in AlgorithmMembers.CurrentCGStatistics.
void PStep_Algorithm::ColumnGeneration(){

	//Initialization:
	//------------------------------------------------------------------------------------

	//Initialize support variables
	timePoint StartTime, tic, toc;
	short int NumberOfCustomers = AlgorithmMembers.Instance.NumberOfCustomers;
	short int NumberOfLocations = NumberOfCustomers + 1;
	short int NumberOfEdges = (NumberOfLocations - 1)*NumberOfLocations / 2;
	int CurrentxIndex = 0;
	//Create references to objects for easy reference
	PStep_ColumnGenerator& ColumnGenerator = *(CurrentColumnGenerator);
	PStep_ColumnGenerationStatistics& CurrentCGStatistics = AlgorithmMembers.CurrentCGStatistics;
	PStep_Solution& CurrentSolution = AlgorithmMembers.CurrentSolution;
	std::list<std::list<PStep>::iterator>& CurrentPSteps = AlgorithmMembers.CurrentPSteps;
	std::list<int>& CurrentPStepsxIndices = AlgorithmMembers.CurrentPStepsxIndices;

	//Obtain the CurrentTreeNode for easy reference
	std::list<PStep_TreeNode>::iterator CurrentTreeNode = AlgorithmMembers.CurrentTreeNode;
	//Start Algorithm timer
	StartTime = HomeMadeclock::now();

	//Initialize p-steps for the current node
	SelectCurrentPSteps();

	//Select the available arcs stored in ColumnGenerator for the current node
	SetAvailableArcs();

	//Reinitialize the decremental state-space information, if required
	if (AlgorithmMembers.DSSRUsed && AlgorithmMembers.AlgorithmSettings.DSSRReset == 2){
		ColumnGenerator.ResetDSSRNeighbourhoods();
	}

	//Reinitialize the current solution
	CurrentSolution.Reinitialize();

	//Reinitialize the current CGStatistics
	CurrentCGStatistics.EmptyStatistics();
	CurrentCGStatistics.StartNumberOfPSteps = AlgorithmMembers.NumberOfCurrentPSteps;
	//StartNumerOfCC, StartNumberOf2EH, StartNumberOfSComb, StartNumberOfFCC and StartNumberOfHMC
	//are set later in this function.

	//Generate ILOG model and CPLEX object
	//------------------------------------------------------------------------------------

	//Initialize Ilog environment
	IloEnv env;

	//Initialize objective function
	IloObjective Objective = IloMinimize(env);

	//Initialize constraints:

	//Initialize the customer assignment constraints
	//CustomerAssignmentConstraint[i] is the constraint of customer i.
	//The first customer is i=0.
	IloRangeArray CustomerAssignmentConstraints(env, NumberOfCustomers, 2, 2);
	//Initialize the degree constraints
	//DegreeConstraints[i] is the constraint of customer i.
	//The first customer is i=0.
	IloRangeArray DegreeConstraints(env, NumberOfCustomers, 0, 0);
	//Initialize the demand flow constraints.
	//DemandFlowConstraints[i] is the constraint of customer i.
	//The first customer is i=0.
	IloRangeArray DemandFlowConstraints(env, NumberOfCustomers, 0, IloInfinity);
	//Initialize the edge assignment constraints
	//EdgeFlowConstraints[(j-1)*j/2+i] is the constraint of edge (i,j) for i<j.
	//Here, i=0 represents the depot.
	IloRangeArray EdgeFlowConstraints(env, NumberOfEdges, 0, 0);

	//Capacity cuts:
	//Initialize the capacity cuts
	IloRangeArray CapacityCuts(env);

	//2 edge hypotour cuts:
	//Initialize the 2 edge hypotour cuts
	IloRangeArray TwoEdgeHypotourCuts(env);

	//Strengthened comb cuts:
	//Initialize the strengthened comb cuts
	IloRangeArray StrengthenedCombCuts(env);

	//Framed capacity cuts:
	//Initialize the framedcapacity cuts
	IloRangeArray FramedCapacityCuts(env);

	//Homogeneous multistar cuts:
	//Initialize the homogeneous multistar cuts
	IloRangeArray HomogeneousMultistarCuts(env);

	//Branching cuts:
	//Initialize the cuts that result from branching decisions
	//NOTE: If the number of vehicles is fixed, this constraint is added as a BranchingCut.
	IloRangeArray BranchingCuts(env);

	//Initialize the p-step variables, one for each p-step in CurrentPSteps
	IloNumVarArray x(env);
	//Add x's
	for (auto CurrentPStep = CurrentPSteps.begin(); CurrentPStep != CurrentPSteps.end(); ++CurrentPStep){
		//Initialize column for the new x
		IloNumColumn TempxColumn(env);
		//Create column
		CreatePStepColumn(**CurrentPStep, Objective, CustomerAssignmentConstraints, DegreeConstraints, DemandFlowConstraints, EdgeFlowConstraints, TempxColumn);
		//Add the column to the model
		x.add(IloNumVar(TempxColumn, 0, IloInfinity));
		//Update the index of the current p-step to correspond to this x variable.
		CurrentPStepsxIndices.push_back(CurrentxIndex);
		++CurrentxIndex;
		//clean
		TempxColumn.end();
	}

	//Initialize edge flow variables
	//theta[(j-1)j/2+i] provides the edge flow on (i,j) for i<j.
	IloNumVarArray theta(env);
	//Initialize counter indiciating the current theta edge index.
	short int Counter = 0;
	//Go over all customers (as end location of the edge)
	for (short int j = 1; j < NumberOfLocations; ++j){
		
		//Add the edge (depot, customer j)
		//Initialize column for the new to be added arc flow variable
		IloNumColumn TempThetaDepotColumn(env);
		//Add arc flow variable to the arcassignment constraints
		TempThetaDepotColumn += EdgeFlowConstraints[Counter](-1);
		//Update Counter
		++Counter;

		//Add variable
		theta.add(IloNumVar(TempThetaDepotColumn, 0, 2));
		TempThetaDepotColumn.end();

		//Go over all customer locations
		for (short int i = 1; i < j; ++i){

			//Initialize column for the new to be added edge flow variable
			IloNumColumn TempThetaColumn(env);
			//Add edge flow variable to the edge assignment constraints
			TempThetaColumn += EdgeFlowConstraints[Counter](-1);
			//Update Counter
			++Counter;

			//Add variable
			theta.add(IloNumVar(TempThetaColumn, 0, 1));
			TempThetaColumn.end();
		}

	}

	//Construct a variable that has a nonzero value if the solution space corresponding to the current branch
	//is void of a feasible integer solution.
	IloNumVarArray FeasibilityDummy(env);
	//Initialize column corresponding to the FeasibilityDummy
	IloNumColumn TempFeasibilityDummyColumn(env);
	//Add dummy to the objective function
	TempFeasibilityDummyColumn += Objective(AlgorithmMembers.FeasibilityDummyObjectiveCoeficient);
	//Add dummy to the assignment constraints
	for (short int i = 0; i<NumberOfCustomers; i++){
		TempFeasibilityDummyColumn += CustomerAssignmentConstraints[i](2);
	}
	//Add dummy
	FeasibilityDummy.add(IloNumVar(TempFeasibilityDummyColumn, 0, IloInfinity));
	TempFeasibilityDummyColumn.end();
	
	//Create capacity cuts
	//Obtain counter, indicating the number of CC's added so far.
	int &CurrentTotalNumberOfCC = AlgorithmMembers.CurrentCGStatistics.TotalNumberOfCC;
	CurrentTotalNumberOfCC = 0;
	//Go over all capacity cuts
	for (auto CurrentCut = AlgorithmMembers.CapacityCuts.begin(); CurrentCut != AlgorithmMembers.CapacityCuts.end(); ++CurrentCut){
		//Initialize cut
		IloExpr TempCut(env);
		//Create cut expression
		CreateCutExpression_ByEdges(*CurrentCut, theta, TempCut);
		//Add dummy
		TempCut += (CurrentCut->Bound)*FeasibilityDummy[0];
		//Add cut
		CapacityCuts.add(TempCut <= CurrentCut->Bound);
		//Update CCCounter
		++CurrentTotalNumberOfCC;
		//clean
		TempCut.end();
	}
	//Set StartNumberOfCC
	CurrentCGStatistics.StartNumberOfCC = CurrentTotalNumberOfCC;

	//Create 2 edge hypotour cuts
	//Obtain counter, indicating the number of 2EH's added so far.
	int &CurrentTotalNumberOf2EH = AlgorithmMembers.CurrentCGStatistics.TotalNumberOf2EH;
	CurrentTotalNumberOf2EH = 0;
	//Go over all 2 edge hypotour cuts
	for (auto CurrentCut = AlgorithmMembers.TwoEdgeHypotourCuts.begin(); CurrentCut != AlgorithmMembers.TwoEdgeHypotourCuts.end(); ++CurrentCut){
		//Initialize cut
		IloExpr TempCut(env);
		//Create cut expression
		CreateCutExpression_ByEdgesAndCoefficients(*CurrentCut, theta, TempCut);
		//Add dummy
		TempCut += (CurrentCut->Bound)*FeasibilityDummy[0];
		//Add cut
		TwoEdgeHypotourCuts.add(TempCut <= CurrentCut->Bound);
		//Update 2EHCounter
		++CurrentTotalNumberOf2EH;
		//clean
		TempCut.end();
	}
	//Set StartNumberOf2EH
	CurrentCGStatistics.StartNumberOf2EH = CurrentTotalNumberOf2EH;

	//Create strengthened comb cuts
	//Obtain counter, indicating the number of SComb's added so far.
	int &CurrentTotalNumberOfSComb = AlgorithmMembers.CurrentCGStatistics.TotalNumberOfSComb;
	CurrentTotalNumberOfSComb = 0;
	//Go over all strengthened comb cuts
	for (auto CurrentCut = AlgorithmMembers.StrengthenedCombCuts.begin(); CurrentCut != AlgorithmMembers.StrengthenedCombCuts.end(); ++CurrentCut){
		//Initialize cut
		IloExpr TempCut(env);
		//Create cut expression
		CreateCutExpression_ByEdgesAndCoefficients(*CurrentCut, theta, TempCut);
		//Add dummy
		TempCut += (CurrentCut->Bound)*FeasibilityDummy[0];
		//Add cut
		StrengthenedCombCuts.add(TempCut >= CurrentCut->Bound);
		//Update SCombCounter
		++CurrentTotalNumberOfSComb;
		//clean
		TempCut.end();
	}
	//Set StartNumberOfSComb
	CurrentCGStatistics.StartNumberOfSComb = CurrentTotalNumberOfSComb;

	//Create framed capacity cuts
	//Obtain counter, indicating the number of FCCs added so far.
	int &CurrentTotalNumberOfFCC = AlgorithmMembers.CurrentCGStatistics.TotalNumberOfFCC;
	CurrentTotalNumberOfFCC = 0;
	//Go over all framed capacity cuts
	for (auto CurrentCut = AlgorithmMembers.FramedCapacityCuts.begin(); CurrentCut != AlgorithmMembers.FramedCapacityCuts.end(); ++CurrentCut){
		//Initialize cut
		IloExpr TempCut(env);
		//Create cut expression
		CreateCutExpression_ByEdges(*CurrentCut, theta, TempCut);
		//Add dummy
		TempCut += (CurrentCut->Bound)*FeasibilityDummy[0];
		//Add cut
		FramedCapacityCuts.add(TempCut >= CurrentCut->Bound);
		//Update FCCCounter
		++CurrentTotalNumberOfFCC;
		//clean
		TempCut.end();
	}
	//Set StartNumberOfFCC
	CurrentCGStatistics.StartNumberOfFCC = CurrentTotalNumberOfFCC;

	//Create homogeneous multistar cuts
	//Obtain counter, indicating the number of HMC's added so far.
	int &CurrentTotalNumberOfHMC = AlgorithmMembers.CurrentCGStatistics.TotalNumberOfHMC;
	CurrentTotalNumberOfHMC = 0;
	//Go over all homogeneous multistar cuts
	for (auto CurrentCut = AlgorithmMembers.HomogeneousMultistarCuts.begin(); CurrentCut != AlgorithmMembers.HomogeneousMultistarCuts.end(); ++CurrentCut){
		//Initialize cut
		IloExpr TempCut(env);
		//Create cut expression
		CreateCutExpression_ByEdgesAndCoefficients(*CurrentCut, theta, TempCut);
		//Add dummy
		TempCut += (CurrentCut->Bound)*FeasibilityDummy[0];
		//Add cut
		HomogeneousMultistarCuts.add(TempCut >= CurrentCut->Bound);
		//Update HMC counter
		++CurrentTotalNumberOfHMC;
		//clean
		TempCut.end();
	}
	//Set StartNumberOfHMC
	CurrentCGStatistics.StartNumberOfHMC = CurrentTotalNumberOfHMC;

	//Add the various types of branching cuts:
	//Add branching cuts type 1: By nodes, lower bound
	for (auto CurrentCut = CurrentTreeNode->EdgeFlowCuts_ByNodes_LB.begin(); CurrentCut != CurrentTreeNode->EdgeFlowCuts_ByNodes_LB.end(); ++CurrentCut){
		//Initialize cut
		IloExpr TempCut(env);
		//Create cut expression
		CreateCutExpression_ByNodes(*CurrentCut, theta, TempCut);
		//Add dummy
		TempCut += (CurrentCut->Bound)*FeasibilityDummy[0];
		//Add cut
		BranchingCuts.add(TempCut >= CurrentCut->Bound);
		//clean
		TempCut.end();
	}
	//Add branching cuts type 2: By nodes, equality
	for (auto CurrentCut = CurrentTreeNode->EdgeFlowCuts_ByNodes_EQ.begin(); CurrentCut != CurrentTreeNode->EdgeFlowCuts_ByNodes_EQ.end(); ++CurrentCut){
		//Initialize cut
		IloExpr TempCut(env);
		//Create cut expression
		CreateCutExpression_ByNodes(*CurrentCut, theta, TempCut);
		//Add dummy
		TempCut += (CurrentCut->Bound)*FeasibilityDummy[0];
		//Add cut
		BranchingCuts.add(TempCut == CurrentCut->Bound);
		//clean
		TempCut.end();
	}
	//Add branching cuts type 3: By nodes, upper bound
	for (auto CurrentCut = CurrentTreeNode->EdgeFlowCuts_ByNodes_UB.begin(); CurrentCut != CurrentTreeNode->EdgeFlowCuts_ByNodes_UB.end(); ++CurrentCut){
		//Initialize cut
		IloExpr TempCut(env);
		//Create cut expression
		CreateCutExpression_ByNodes(*CurrentCut, theta, TempCut);
		//Add dummy
		TempCut += (CurrentCut->Bound)*FeasibilityDummy[0];
		//Add cut
		BranchingCuts.add(TempCut <= CurrentCut->Bound);
		//clean
		TempCut.end();
	}
	//Add the various types of branching cuts:
	//Add branching cuts type 4: By edges, lower bound
	for (auto CurrentCut = CurrentTreeNode->EdgeFlowCuts_ByEdges_LB.begin(); CurrentCut != CurrentTreeNode->EdgeFlowCuts_ByEdges_LB.end(); ++CurrentCut){
		//Initialize cut
		IloExpr TempCut(env);
		//Create cut expression
		CreateCutExpression_ByEdges(*CurrentCut, theta, TempCut);
		//Add dummy
		TempCut += (CurrentCut->Bound)*FeasibilityDummy[0];
		//Add cut
		BranchingCuts.add(TempCut >= CurrentCut->Bound);
		//clean
		TempCut.end();
	}
	//Add branching cuts type 5: By edges, equality
	for (auto CurrentCut = CurrentTreeNode->EdgeFlowCuts_ByEdges_EQ.begin(); CurrentCut != CurrentTreeNode->EdgeFlowCuts_ByEdges_EQ.end(); ++CurrentCut){
		//Initialize cut
		IloExpr TempCut(env);
		//Create cut expression
		CreateCutExpression_ByEdges(*CurrentCut, theta, TempCut);
		//Add dummy
		TempCut += (CurrentCut->Bound)*FeasibilityDummy[0];
		//Add cut
		BranchingCuts.add(TempCut == CurrentCut->Bound);
		//clean
		TempCut.end();
	}
	//Add branching cuts type 6: By edges, upper bound
	for (auto CurrentCut = CurrentTreeNode->EdgeFlowCuts_ByEdges_UB.begin(); CurrentCut != CurrentTreeNode->EdgeFlowCuts_ByEdges_UB.end(); ++CurrentCut){
		//Initialize cut
		IloExpr TempCut(env);
		//Create cut expression
		CreateCutExpression_ByEdges(*CurrentCut, theta, TempCut);
		//Add dummy
		TempCut += (CurrentCut->Bound)*FeasibilityDummy[0];
		//Add cut
		BranchingCuts.add(TempCut <= CurrentCut->Bound);
		//clean
		TempCut.end();
	}

	//Construct model
	IloModel PStepModel(env);
	PStepModel.add(Objective);
	PStepModel.add(CustomerAssignmentConstraints);
	PStepModel.add(DegreeConstraints);
	PStepModel.add(DemandFlowConstraints);
	PStepModel.add(EdgeFlowConstraints);
	PStepModel.add(CapacityCuts);
	PStepModel.add(TwoEdgeHypotourCuts);
	PStepModel.add(StrengthenedCombCuts);
	PStepModel.add(FramedCapacityCuts);
	PStepModel.add(HomogeneousMultistarCuts);
	PStepModel.add(BranchingCuts);

	//Construct CPLEX object to solve the p-step formulation of the CVRP
	IloCplex PStepLP(PStepModel);
	//Disallow cplex to print data on screen.
	PStepLP.setOut(env.getNullStream());

	//Solve the LP relaxation of CVRP
	//------------------------------------------------------------------------------------

	//Initialize indicator, indicating whether the algorithm should terminate
	//0: Do not terminate
	//1: Optimum found
	//2: Time exceeded
	short int TerminationFlag = 0;

	//Initialize iterations
	CurrentCGStatistics.TotalNumberOfIterations = 1;

	//Iterate untill some termination criterion is met
	while (TerminationFlag == 0){

		//Initialize current iteration statistics
		PStep_SingleCGIterationStatistics CurrentIterationStatistics;
		CurrentIterationStatistics.NumberOfIncludedPSteps = AlgorithmMembers.NumberOfCurrentPSteps;
		CurrentIterationStatistics.NumberOfPStepsInTheShortList = AlgorithmMembers.NumberOfPStepsInShortList;

		//Solve master problem:
		//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

		//Get LP starttime
		tic = HomeMadeclock::now();

		//Solve the LP relaxation
		PStepLP.solve();

		//Get LP Stop time
		toc = HomeMadeclock::now();

		//Process solution
		//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

		//Update the current solution:
		//----------------------------
		CurrentSolution.Reinitialize();

		//Update solution construction time
		CurrentSolution.ConstructionTime=HomeMadeclock::now();

		//Update route selection:		
		//NOTE: Below the values are taken one by one, it is more efficient
		//		to do this in one go:
		//		IloNumArray tempx(env);
		//		PStepLP.getValues(tempx, x);
		//		However, when columns are removed by the ColumnManager, this cannot be done anymore.

		//Initialize list of x variables that are still in use before column generation
		IloNumVarArray xInUseBeforeCG(env);
		//Go over all x variables that are in use and add them to the list
		for (auto CurrentPStepxIndex = CurrentPStepsxIndices.begin(); CurrentPStepxIndex != CurrentPStepsxIndices.end(); ++CurrentPStepxIndex){
			//Add the current x variable to the list
			xInUseBeforeCG.add(x[*CurrentPStepxIndex]);
		}
		//Obtain values of the x variables that are in use
		IloNumArray tempx(env);
		PStepLP.getValues(tempx, xInUseBeforeCG);

		//Initialize counter
		int Counter = 0;
		//Go over all current p-steps
		for (auto CurrentPStep = CurrentPSteps.begin(); CurrentPStep != CurrentPSteps.end(); ++CurrentPStep){
			//Check whether the current p-step was used
			if (tempx[Counter] > AlgorithmMembers.Tolerance){
				//Add the route to the solution
				CurrentSolution.PSteps.push_back(**CurrentPStep);
				//Add the fraction with which the solution was selected to the current solution
				CurrentSolution.PStepSelection.push_back(tempx[Counter]);
			}
			//Update counter
			++Counter;
		}

		//Obtain the edge flow variables
		IloNumArray tempTheta(env);
		PStepLP.getValues(tempTheta, theta);
		//Update the edge flow in the current solution
		for (short int i = 0; i < NumberOfEdges; ++i){
			CurrentSolution.EdgeFraction[i] = 0;
			CurrentSolution.EdgeFraction[i] = tempTheta[i];
		}

		//Obtain solution value
		CurrentSolution.SolutionValue = PStepLP.getObjValue();

		//clean
		tempx.end();
		tempTheta.end();

		//Update CG statistics:
		//---------------------

		//Update time
		CurrentIterationStatistics.LPTime = HomeMadeclock::seconds(tic, toc);
		CurrentCGStatistics.TotalLPTime += CurrentIterationStatistics.LPTime;
		//Obtain solution value
		CurrentIterationStatistics.Value=CurrentSolution.SolutionValue;

		//Print progress
		if (AlgorithmMembers.AlgorithmSettings.PrintCGProgress){
			std::cout << CurrentCGStatistics.TotalNumberOfIterations << ": " << AlgorithmMembers.NumberOfCurrentPSteps << " p-steps, LB " << CurrentSolution.SolutionValue << std::endl;
		}

		//Update duals:
		//-------------
		//Obtain multipliers from PStepLP
		IloNumArray tempLambda(env);
		IloNumArray tempMu(env);
		IloNumArray tempNu(env);
		IloNumArray tempPi(env);
		PStepLP.getDuals(tempLambda, CustomerAssignmentConstraints);
		PStepLP.getDuals(tempMu, DegreeConstraints);
		PStepLP.getDuals(tempNu, DemandFlowConstraints);
		PStepLP.getDuals(tempPi, EdgeFlowConstraints);
		//Update duals in ColumnGenerator:
		//Go over all customers
		for (short int i = 0; i < NumberOfCustomers; ++i){
			ColumnGenerator.Lambda[i] = tempLambda[i];
			ColumnGenerator.Mu[i] = tempMu[i];
			ColumnGenerator.Nu[i] = tempNu[i];
		}
		//Go over all edges
		for (short int i = 0; i < NumberOfEdges; ++i){
			ColumnGenerator.Pi[i] = tempPi[i];
		}
		//clean:
		tempLambda.end();
		tempMu.end();
		tempNu.end();
		tempPi.end();

		//Obtain reduced costs (for column management)
		IloNumArray xReducedCosts(env,AlgorithmMembers.NumberOfCurrentPSteps);
		PStepLP.getReducedCosts(xReducedCosts, xInUseBeforeCG);

		//Solve pricing problem
		//****************************************************************************

		//Obtain the number of p-steps currently added
		int NumberOfCurrentPStepsBeforeColumnGeneration = AlgorithmMembers.NumberOfCurrentPSteps;
		//Obtain iterator pointing to the last current p-step that was added.
		auto LastAddedCurrentPStepBeforeCG = CurrentPSteps.end();
		--LastAddedCurrentPStepBeforeCG;

		//Get pricing starttime
		tic=HomeMadeclock::now();

		//Solve the pricing problem:
		//First check shortlist
		if (!(CurrentColumnManager->SearchShortList((*CurrentColumnGenerator)))){
			//No negative reduced cost columns were found in the shortlist

			//Solve the pricing problems to generate new p-steps
			CurrentColumnGenerator->GenerateColumns();

		}

		//Get Pricing Stop time
		toc = HomeMadeclock::now();

		//Process pricing problem solution
		//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

		//Update the CG statistics
		CurrentIterationStatistics.PricingTime = HomeMadeclock::seconds(tic, toc);
		CurrentCGStatistics.TotalPricingTime += CurrentIterationStatistics.PricingTime;
		CurrentIterationStatistics.NumberOfPStepsGenerated = AlgorithmMembers.NumberOfCurrentPSteps - NumberOfCurrentPStepsBeforeColumnGeneration;
		//Initialize the seperation statistics:
		CurrentIterationStatistics.SeparateCCTime=0;
		CurrentIterationStatistics.NumberOfCC = 0;
		CurrentIterationStatistics.Separate2EHTime = 0;
		CurrentIterationStatistics.NumberOf2EH = 0;
		CurrentIterationStatistics.SeparateSCombTime = 0;
		CurrentIterationStatistics.NumberOfSComb = 0;
		CurrentIterationStatistics.SeparateFCCTime = 0;
		CurrentIterationStatistics.NumberOfFCC = 0;
		CurrentIterationStatistics.SeparateHMCTime = 0;
		CurrentIterationStatistics.NumberOfHMC = 0;

		//Check if p-steps have been added
		if (NumberOfCurrentPStepsBeforeColumnGeneration < AlgorithmMembers.NumberOfCurrentPSteps){

			//Add x's
			//Initialize CurrentPStep
			auto CurrentPStep = LastAddedCurrentPStepBeforeCG;
			++CurrentPStep;
			//Go over all new p-steps. Note that they are added to the front of the list.
			for (; CurrentPStep != CurrentPSteps.end(); ++CurrentPStep){

				//Initialize column for the new x
				IloNumColumn TempxColumn(env);
				//Create column
				CreatePStepColumn(**CurrentPStep, Objective, CustomerAssignmentConstraints, DegreeConstraints, DemandFlowConstraints, EdgeFlowConstraints, TempxColumn);
				//Add the column to the model
				x.add(IloNumVar(TempxColumn, 0, IloInfinity));
				//Update the index of the current p-step to correspond to this x variable.
				CurrentPStepsxIndices.push_back(CurrentxIndex);
				++CurrentxIndex;
				//clean
				TempxColumn.end();

			}

		}
		else{
			//No new p-steps have been found

			//Valid inequalities
			//****************************************************************************

			//Update time and value before adding inequalities
			if (CurrentCGStatistics.ValueBeforeInequalities == 0){
				CurrentCGStatistics.ValueBeforeInequalities = CurrentSolution.SolutionValue;
				toc = HomeMadeclock::now();
				CurrentCGStatistics.TimeBeforeInequalities = HomeMadeclock::seconds(StartTime,toc);
			}

			//Initialize indicator indicating whether new valid inequalities have been added
			bool NewValidInequalitiesAdded = false;

			//Obtain iterator of the last capacity cut.
			//It is initialized if CCCounter>0.
			std::list<PStep_EdgeFlowCut>::iterator LastCapacityCut;
			if (CurrentTotalNumberOfCC > 0){
				LastCapacityCut = --AlgorithmMembers.CapacityCuts.end();
			}
			//Obtain the number of cc's and the total CC separation time before attempting to separate again.
			int NumberOfCCBeforeSeparation = AlgorithmMembers.NumberOfCapacityCuts;
			double LastTotalCCSeparationTime = AlgorithmMembers.CurrentCGStatistics.TotalSeparateCCTime;
			
			//Obtain iterator of the last 2 edge hypotour cut.
			//It is initialized if CurrentTotalNumberOf2EH>0.
			std::list<PStep_EdgeAndCoefficientCut>::iterator Last2EHCut;
			if (CurrentTotalNumberOf2EH > 0){
				Last2EHCut = --AlgorithmMembers.TwoEdgeHypotourCuts.end();
			}
			//Obtain the total number of 2EH's and the 2EH separation time before attempting to separate again.
			int NumberOf2EHBeforeSeparation=AlgorithmMembers.NumberOf2EHCuts;
			double LastTotal2EHSeparationTime = AlgorithmMembers.CurrentCGStatistics.TotalSeparate2EHTime;

			//Obtain iterator of the last strengthened comb cut.
			//It is initialized if CurrentTotalNumberOfSComb>0.
			std::list<PStep_EdgeAndCoefficientCut>::iterator LastSCombCut;
			if (CurrentTotalNumberOfSComb > 0){
				LastSCombCut = --AlgorithmMembers.StrengthenedCombCuts.end();
			}
			//Obtain the total number of SComb's and the SComb separation time before attempting to separate again.
			int NumberOfSCombBeforeSeparation = AlgorithmMembers.NumberOfSCombCuts;
			double LastTotalSCombSeparationTime = AlgorithmMembers.CurrentCGStatistics.TotalSeparateSCombTime;

			//Obtain iterator of the last framed capacity cut.
			//It is initialized if FCCCounter>0.
			std::list<PStep_EdgeFlowCut>::iterator LastFramedCapacityCut;
			if (CurrentTotalNumberOfFCC > 0){
				LastFramedCapacityCut = --AlgorithmMembers.FramedCapacityCuts.end();
			}
			//Obtain the number of FCC's and the total FCC separation time before attempting to separate again.
			int NumberOfFCCBeforeSeparation = AlgorithmMembers.NumberOfFramedCapacityCuts;
			double LastTotalFCCSeparationTime = AlgorithmMembers.CurrentCGStatistics.TotalSeparateFCCTime;

			//Obtain iterator of the last homogeneous multistar cut.
			//It is initialized if CurrentTotalNumberOfHMC>0.
			std::list<PStep_EdgeAndCoefficientCut>::iterator LastHomogeneousMultistarCut;
			if (CurrentTotalNumberOfHMC > 0){
				LastHomogeneousMultistarCut = --AlgorithmMembers.HomogeneousMultistarCuts.end();
			}
			//Obtain the total number of HMC's and the HMC separation time before attempting to separate again.
			int NumberOfHMCBeforeSeparation = AlgorithmMembers.NumberOfHomogeneousMultistarCuts;
			double LastTotalHMCSeparationTime = AlgorithmMembers.CurrentCGStatistics.TotalSeparateHMCTime;

			//Separate
			//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
			
			//Separate
			CurrentSeparationAlgorithm->Separate();

			//Process output of the separation routine:
			//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
			
			//Check whether violated capacity cuts have been identified
			if (AlgorithmMembers.NumberOfCapacityCuts>NumberOfCCBeforeSeparation){
				//Violated capacity cuts have been identified, add them
				
				//Obtain indicator of the potentially new first added capacity cut
				std::list<PStep_EdgeFlowCut>::iterator FirstNewCapacityCut;
				//Check whether capacity cuts were previously added
				if (NumberOfCCBeforeSeparation == 0){
					//No cuts were previously added.
					FirstNewCapacityCut = AlgorithmMembers.CapacityCuts.begin();
				}
				else{
					//Cuts have been added before.
					FirstNewCapacityCut = LastCapacityCut;
					++FirstNewCapacityCut;
				}

				//Go over the new capacity cuts
				for (auto CurrentCut = FirstNewCapacityCut; CurrentCut != AlgorithmMembers.CapacityCuts.end(); ++CurrentCut){
					//Initialize cut
					IloExpr TempCut(env);
					//Create cut expression
					CreateCutExpression_ByEdges(*CurrentCut, theta, TempCut);
					//Add dummy
					TempCut += (CurrentCut->Bound)*FeasibilityDummy[0];
					//Add cut
					CapacityCuts.add(TempCut <= CurrentCut->Bound);
					PStepModel.add(CapacityCuts[CurrentTotalNumberOfCC]);
					//Update the number of capacity cuts added
					++CurrentTotalNumberOfCC;
					//Update the number of capacity cuts added during this iteration.
					++CurrentIterationStatistics.NumberOfCC;
					//clean
					TempCut.end();
				}

				//Update separation time
				CurrentIterationStatistics.SeparateCCTime = CurrentCGStatistics.TotalSeparateCCTime - LastTotalCCSeparationTime;

				//Update indicator stating that new violatedInequalities are added
				NewValidInequalitiesAdded = true;

			}

			//Check whether violated 2 edge hypotour cuts have been identified
			if (AlgorithmMembers.NumberOf2EHCuts > NumberOf2EHBeforeSeparation){
				//Violated two edge hypotour cuts have been identified, add them

				//Obtain indicator of the potentially new first added 2 edge hypotour cut
				std::list<PStep_EdgeAndCoefficientCut>::iterator FirstNew2EHCut;
				//Check whether 2 edge hypotour cuts were previously added
				if (NumberOf2EHBeforeSeparation == 0){
					//No cuts were previously added.
					FirstNew2EHCut = AlgorithmMembers.TwoEdgeHypotourCuts.begin();
				}
				else{
					//Cuts have been added before.
					FirstNew2EHCut = Last2EHCut;
					++FirstNew2EHCut;
				}

				//Go over the new capacity cuts
				for (auto CurrentCut = FirstNew2EHCut; CurrentCut != AlgorithmMembers.TwoEdgeHypotourCuts.end(); ++CurrentCut){
					//Initialize cut
					IloExpr TempCut(env);
					//Create cut expression
					CreateCutExpression_ByEdgesAndCoefficients(*CurrentCut, theta, TempCut);
					//Add dummy
					TempCut += (CurrentCut->Bound)*FeasibilityDummy[0];
					//Add cut
					TwoEdgeHypotourCuts.add(TempCut <= CurrentCut->Bound);
					PStepModel.add(TwoEdgeHypotourCuts[CurrentTotalNumberOf2EH]);
					//Update the number of two edge hypotour cuts added
					++CurrentTotalNumberOf2EH;
					//Update the number of two edge hypotour cuts added during this iteration.
					++CurrentIterationStatistics.NumberOf2EH;
					//clean
					TempCut.end();
				}

				//Update separation time
				CurrentIterationStatistics.Separate2EHTime = CurrentCGStatistics.TotalSeparate2EHTime - LastTotal2EHSeparationTime;

				//Update indicator stating that new violatedInequalities are added
				NewValidInequalitiesAdded = true;

			}

			//Check whether violated strengthened comb cuts have been identified
			if (AlgorithmMembers.NumberOfSCombCuts > NumberOfSCombBeforeSeparation){
				//Violated strengthened comb cuts have been identified, add them

				//Obtain indicator of the potentially new first added strengthened comb cut
				std::list<PStep_EdgeAndCoefficientCut>::iterator FirstNewSCombCut;
				//Check whether strengthened comb cuts were previously added
				if (NumberOfSCombBeforeSeparation == 0){
					//No cuts were previously added.
					FirstNewSCombCut = AlgorithmMembers.StrengthenedCombCuts.begin();
				}
				else{
					//Cuts have been added before.
					FirstNewSCombCut = LastSCombCut;
					++FirstNewSCombCut;
				}

				//Go over the new capacity cuts
				for (auto CurrentCut = FirstNewSCombCut; CurrentCut != AlgorithmMembers.StrengthenedCombCuts.end(); ++CurrentCut){
					//Initialize cut
					IloExpr TempCut(env);
					//Create cut expression
					CreateCutExpression_ByEdgesAndCoefficients(*CurrentCut, theta, TempCut);
					//Add dummy
					TempCut += (CurrentCut->Bound)*FeasibilityDummy[0];
					//Add cut
					StrengthenedCombCuts.add(TempCut >= CurrentCut->Bound);
					PStepModel.add(StrengthenedCombCuts[CurrentTotalNumberOfSComb]);
					//Update the number of strengthened comb cuts added
					++CurrentTotalNumberOfSComb;
					//Update the number of strengthened comb cuts added during this iteration.
					++CurrentIterationStatistics.NumberOfSComb;
					//clean
					TempCut.end();
				}

				//Update separation time
				CurrentIterationStatistics.SeparateSCombTime = CurrentCGStatistics.TotalSeparateSCombTime - LastTotalSCombSeparationTime;

				//Update indicator stating that new violatedInequalities are added
				NewValidInequalitiesAdded = true;

			}

			//Check whether violated framed capacity cuts have been identified
			if (AlgorithmMembers.NumberOfFramedCapacityCuts > NumberOfFCCBeforeSeparation){
				//Violated framed capacity cuts have been identified, add them

				//Obtain indicator of the potentially new first added framed capacity cut
				std::list<PStep_EdgeFlowCut>::iterator FirstNewFCCCut;
				//Check whether framed capacity cuts were previously added
				if (NumberOfFCCBeforeSeparation == 0){
					//No cuts were previously added.
					FirstNewFCCCut = AlgorithmMembers.FramedCapacityCuts.begin();
				}
				else{
					//Cuts have been added before.
					FirstNewFCCCut = LastFramedCapacityCut;
					++FirstNewFCCCut;
				}

				//Go over the new framed capacity cuts
				for (auto CurrentCut = FirstNewFCCCut; CurrentCut != AlgorithmMembers.FramedCapacityCuts.end(); ++CurrentCut){
					//Initialize cut
					IloExpr TempCut(env);
					//Create cut expression
					CreateCutExpression_ByEdges(*CurrentCut, theta, TempCut);
					//Add dummy
					TempCut += (CurrentCut->Bound)*FeasibilityDummy[0];
					//Add cut
					FramedCapacityCuts.add(TempCut >= CurrentCut->Bound);
					PStepModel.add(FramedCapacityCuts[CurrentTotalNumberOfFCC]);
					//Update the number of framed capacity cuts added
					++CurrentTotalNumberOfFCC;
					//Update the number of framed capacity cuts added during this iteration.
					++CurrentIterationStatistics.NumberOfFCC;
					//clean
					TempCut.end();
				}

				//Update separation time
				CurrentIterationStatistics.SeparateFCCTime = CurrentCGStatistics.TotalSeparateFCCTime - LastTotalFCCSeparationTime;

				//Update indicator stating that new violatedInequalities are added
				NewValidInequalitiesAdded = true;

			}

			//Check whether violated homogeneous multistar cuts have been identified
			if (AlgorithmMembers.NumberOfHomogeneousMultistarCuts > NumberOfHMCBeforeSeparation){
				//Violated homogeneous multistar cuts have been identified, add them

				//Obtain indicator of the potentially new first added homogeneous multistar cut
				std::list<PStep_EdgeAndCoefficientCut>::iterator FirstNewHomogeneousMultistarCut;
				//Check whether strengthened comb cuts were previously added
				if (NumberOfHMCBeforeSeparation == 0){
					//No cuts were previously added.
					FirstNewHomogeneousMultistarCut = AlgorithmMembers.HomogeneousMultistarCuts.begin();
				}
				else{
					//Cuts have been added before.
					FirstNewHomogeneousMultistarCut = LastHomogeneousMultistarCut;
					++FirstNewHomogeneousMultistarCut;
				}

				//Go over the new capacity cuts
				for (auto CurrentCut = FirstNewHomogeneousMultistarCut; CurrentCut != AlgorithmMembers.HomogeneousMultistarCuts.end(); ++CurrentCut){
					//Initialize cut
					IloExpr TempCut(env);
					//Create cut expression
					CreateCutExpression_ByEdgesAndCoefficients(*CurrentCut, theta, TempCut);
					//Add dummy
					TempCut += (CurrentCut->Bound)*FeasibilityDummy[0];
					//Add cut
					HomogeneousMultistarCuts.add(TempCut >= CurrentCut->Bound);
					PStepModel.add(HomogeneousMultistarCuts[CurrentTotalNumberOfHMC]);
					//Update the number of homogeneous multistar cuts added
					++CurrentTotalNumberOfHMC;
					//Update the number of homogeneous multistar cuts added during this iteration.
					++CurrentIterationStatistics.NumberOfHMC;
					//clean
					TempCut.end();
				}

				//Update separation time
				CurrentIterationStatistics.SeparateHMCTime = CurrentCGStatistics.TotalSeparateHMCTime - LastTotalHMCSeparationTime;

				//Update indicator stating that new violatedInequalities are added
				NewValidInequalitiesAdded = true;

			}

			//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
			
			//Check whether new valid inequalities are added
			if (!NewValidInequalitiesAdded){
				//No routes are found and no valid inequalities are added

				//Perform phase two column generation:
				//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

				//Get pricing starttime
				tic = HomeMadeclock::now();
				
				//Generate columns
				CurrentColumnGenerator->GenerateColumnsPhaseTwo();

				//Get Pricing Stop time
				toc = HomeMadeclock::now();

				//Process pricing problem solution:

				//Update the CG statistics
				CurrentIterationStatistics.PricingTime += HomeMadeclock::seconds(tic, toc);
				CurrentCGStatistics.TotalPricingTime += HomeMadeclock::seconds(tic, toc);
				CurrentIterationStatistics.NumberOfPStepsGenerated = AlgorithmMembers.NumberOfCurrentPSteps - NumberOfCurrentPStepsBeforeColumnGeneration;

				//Check if p-steps have been added
				if (NumberOfCurrentPStepsBeforeColumnGeneration < AlgorithmMembers.NumberOfCurrentPSteps){

					//Add x's
					//Initialize CurrentPStep
					auto CurrentPStep = LastAddedCurrentPStepBeforeCG;
					++CurrentPStep;
					//Go over all new p-steps. Note that they are added to the front of the list.
					for (; CurrentPStep != CurrentPSteps.end(); ++CurrentPStep){

						//Initialize column for the new x
						IloNumColumn TempxColumn(env);
						//Create column
						CreatePStepColumn(**CurrentPStep, Objective, CustomerAssignmentConstraints, DegreeConstraints, DemandFlowConstraints, EdgeFlowConstraints, TempxColumn);
						//Add the column to the model
						x.add(IloNumVar(TempxColumn, 0, IloInfinity));
						//Update the index of the current p-step to correspond to this x variable.
						CurrentPStepsxIndices.push_back(CurrentxIndex);
						++CurrentxIndex;
						//clean
						TempxColumn.end();

					}

				}
				else{
					//The current solution is optimal.
					TerminationFlag = 1;
				}
			}

		}

		//End of iteration:
		//****************************************************************************
		
		//Update termination flag if this does not already indicate termination
		if (TerminationFlag == 0){
			
			//Remove columns from the formulation:
			//Obtain the end of the CurrentPSteps that should be considered for removal.
			auto ConsideredCurrentPStepEnd = LastAddedCurrentPStepBeforeCG;
			++ConsideredCurrentPStepEnd;
			//Remove
			CurrentColumnManager->Remove(CurrentPSteps.begin(), ConsideredCurrentPStepEnd, x, xReducedCosts, env);

			//Get time
			toc = HomeMadeclock::now();

			//Check time:
			//Calculate time:
			double TotalComputationTime = HomeMadeclock::seconds(StartTime,toc);
			//Check whether the time limit is exceed
			if (TotalComputationTime>AlgorithmMembers.AlgorithmSettings.TimeLimit){
				//Set Flag
				TerminationFlag = 2;
			}
			else{
				//Update support parameters
				CurrentCGStatistics.TotalNumberOfIterations++;
			}

		}

		//Update CurrentCGStatistics by appending the current iteration statistics
		CurrentCGStatistics.IterationStatistics.push_back(CurrentIterationStatistics);

		//clean
		xInUseBeforeCG.end();
		xReducedCosts.end();

	}//end solve while termination criterion is not met

	//Process solution
	//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	
	//Construct flag
	short int Flag;
	if (TerminationFlag == 1){
		//Optimum found, check whether the problem is feasible
		if (PStepLP.getValue(FeasibilityDummy[0])>AlgorithmMembers.Tolerance){
			//The Problem is infeasible
			Flag = 1;
		}
		else{
			//The problem is solved to optimality
			Flag = 0;
		}
	}
	else if (TerminationFlag == 2){
		//Out of time
		Flag = 2;
	}

	//Update CurrentCGStatistics
	CurrentCGStatistics.Flag = Flag;
	CurrentCGStatistics.TotalNumberOfPSteps = AlgorithmMembers.NumberOfCurrentPSteps+AlgorithmMembers.NumberOfPStepsInShortList;
	CurrentCGStatistics.FinalValue = CurrentSolution.SolutionValue;
	//Get LP Stop time
	toc = HomeMadeclock::now();
	CurrentCGStatistics.TotalComputationTime = HomeMadeclock::seconds(StartTime,toc);
	
	//clean
	env.end();

}