/*Implementation of PStep_Algorithm.h

It contains the initialization functions of PStep_Algorithm

*/

#include "PStep_Algorithm.h"

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

//Initialize the algorithm corresponding to an instance of the CVRP.
//Initialize the branch-and-bound tree, i.e., AlgorithmMembers and member functions.
void PStep_Algorithm::Initialize(const CVRPInstance &NewInstance, const PStep_AlgorithmSettings &NewAlgorithmSettings){

	//Delete old members:
	//Initialize support
	short int NumberOfLocations = AlgorithmMembers.Instance.NumberOfCustomers + 1;
	//Delete current pricing problems
	if (AlgorithmMembers.PricingProblems != 0){
		for (short int i = 0; i < NumberOfLocations; ++i){
			delete[] AlgorithmMembers.PricingProblems[i];
		}
		delete[]AlgorithmMembers.PricingProblems;
	}

	//Reinitialize support
	NumberOfLocations = NewInstance.NumberOfCustomers + 1;

	//Initialize general members:
	//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	//Set instance
	AlgorithmMembers.Instance = NewInstance;

	//Set the objective coefficient of the feasibility dummy
	AlgorithmMembers.FeasibilityDummyObjectiveCoeficient = 1;
	for (short int j = 1; j < NumberOfLocations; ++j){
		AlgorithmMembers.FeasibilityDummyObjectiveCoeficient += NewInstance.TravelCosts[0][j] + NewInstance.TravelCosts[j][0];
	}

	//Set algorithm settings
	AlgorithmMembers.AlgorithmSettings = NewAlgorithmSettings;

	//Copy Tolerance
	AlgorithmMembers.Tolerance = NewAlgorithmSettings.Tolerance;

	//Initialize Column generation functions
	//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

	InitializePStep_Adders();
	InitializePStep_PricingProblemSolvers();
	InitializePStep_ColumnGenerators();
	InitializePStep_ColumnManagers();

	//Initialize network members:
	//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	
	//Initialize arcs
	InitializeArcs();

	//Initialize column generator members:
	//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

	//Initialize neighbourhoods as used for ng-path relaxation
	//NOTE: it is assumed that OutArcs are sorted per location from cheap to expensive.
	InitializeNeighbourhoods();

	//Initialize pricing problems
	InitializePricingProblems();

	//Initialize decremental state-space relaxation, if required
	//(can only be done after pricing problems are set).
	if (AlgorithmMembers.DSSRUsed && AlgorithmMembers.AlgorithmSettings.DSSRReset == 1){
		CurrentColumnGenerator->ResetDSSRNeighbourhoods();
	}

	//Initialize valid inequality members
	//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	
	//Set counters
	AlgorithmMembers.NumberOfCapacityCuts = 0;
	AlgorithmMembers.NumberOf2EHCuts = 0;
	AlgorithmMembers.NumberOfSCombCuts = 0;
	AlgorithmMembers.NumberOfFramedCapacityCuts = 0;
	AlgorithmMembers.NumberOfHomogeneousMultistarCuts = 0;

	//Initialize LysgaardCuts
	CMGR_CreateCMgr(&(AlgorithmMembers.LysgaardCuts), AlgorithmMembers.AlgorithmSettings.MaxCCPerIteration);

	//Initialize Separation functions
	//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

	InitializePStep_SeparationAlgorithms();

	//Initialize solution members
	//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

	//Initialize p-steps
	InitializePSteps();

	//Set indicator indicating whether an integer solution is found.
	AlgorithmMembers.IntegerSolutionFound = false;

	//Initialize CurrentSolution
	PStep_Solution TempSolution(NumberOfLocations);
	AlgorithmMembers.CurrentSolution = TempSolution;

	//Initialize branching members
	//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

	//Reset nodes
	AlgorithmMembers.Nodes.clear();
	//Create root node
	PStep_TreeNode InitialNode(1, -1, 1, 0, NumberOfLocations);
	//Potentially add constraints on the number of vehicles as branching cuts
	CreateFixedNumberOfVehiclesCut(InitialNode);
	//Add root node
	AlgorithmMembers.Nodes.push_back(InitialNode);
	AlgorithmMembers.CurrentNumberOfNodes = 1;
	AlgorithmMembers.LastUsedNodeID = 1;

}

//Initialize arcs
void PStep_Algorithm::InitializeArcs(){

	//Obtain number of locations
	short int NumberOfLocations = AlgorithmMembers.Instance.NumberOfCustomers + 1;
	//Obtain travel costs
	double**TravelCosts = AlgorithmMembers.Instance.TravelCosts;
	//Obtain arcs
	std::vector<std::vector<PStep_Arc>>& Arcs = AlgorithmMembers.Arcs;
	std::vector<std::vector<PStep_Arc>>& OutArcs = AlgorithmMembers.OutArcs;
	std::vector<std::vector<PStep_Arc>>& InArcs = AlgorithmMembers.InArcs;

	//Empty Arcs
	Arcs.clear();
	OutArcs.clear();
	InArcs.clear();

	//Construct Arcs:
	//---------------
	//Go over all start locations
	for (short int i = 0; i < NumberOfLocations; ++i){
		//Initialize arcs starting at i.
		std::vector<PStep_Arc> CurrentArcs;
		PStep_Arc CurrentArc;
		CurrentArc.StartNode = i;
		//Go over all end locations with index smaller than i
		for (short int j = 0; j < i; j++){
			//Create arc
			CurrentArc.EndNode = j;
			CurrentArc.Index = (i - 1)*i / 2 + j;
			//Add arc
			CurrentArcs.push_back(CurrentArc);
		}
		//Create dummy arc i-i
		CurrentArc.EndNode = i;
		CurrentArc.Index = -1;
		//Add dummy arc i-i
		CurrentArcs.push_back(CurrentArc);
		//Go over all end locations with index greater than i
		for (short int j = i + 1; j < NumberOfLocations; j++){
			//Create arc
			CurrentArc.EndNode = j;
			CurrentArc.Index = (j - 1)*j / 2 + i;
			//Add arc
			CurrentArcs.push_back(CurrentArc);
		}
		//Store CurrentArcs
		AlgorithmMembers.Arcs.push_back(CurrentArcs);
	}

	//Construct OutArcs:
	//------------------
	//Go over all start locations
	for (short int i = 0; i < NumberOfLocations; ++i){
		
		//Initialize list of arcs
		std::vector<PStep_Arc> CurrentOutArcs;

		//Go over all end locations
		for (short int j = 0; j < NumberOfLocations; ++j){
			
			//Check wheter the start and end location are not the same
			if (i != j){

				//Initialize arc
				PStep_Arc CurrentArc;
				CurrentArc.StartNode = i;
				CurrentArc.EndNode = j;
				CurrentArc.Index = (i - 1)*i / 2 + j;

				//Add the new arc at the correct position in the list of current out arcs:

				//Obtain current cost
				double CurrentArcCost = TravelCosts[i][j];

				//Initialize indicator indicating whether the CurrentArc has not yet been added
				bool CurrentArcHasNotYetBeenAdded = true;
				//Initialize iterator
				auto Arc = CurrentOutArcs.begin();;
				//Go over all current arcs in the list
				while (CurrentArcHasNotYetBeenAdded && Arc != CurrentOutArcs.end()){

					//Check whether CurrentOutArc is cheaper than Arc
					if (CurrentArcCost < TravelCosts[Arc->StartNode][Arc->EndNode]){
						//The cost of CurrentOutArc is lower than of Arc,
						//add CurrentOutArc before Arc in the list.
						CurrentOutArcs.insert(Arc, CurrentArc);
						//Update indicator
						CurrentArcHasNotYetBeenAdded = false;
					}
					else{
						//Update iterator
						++Arc;
					}

				}//end go over all current arcs

				//Check whether the current arc has not yet been added
				if (CurrentArcHasNotYetBeenAdded){
					//The current arc is the most expensive arc out of i, add it at the end of the list
					CurrentOutArcs.push_back(CurrentArc);
				}

			}//end check whether start and end location are the same

		}//end go over all end locations

		//Add CurrentOutArcs
		OutArcs.push_back(CurrentOutArcs);

	}//end go over all start locations

	//Construct InArcs:
	//------------------
	//Go over all end locations
	for (short int j = 0; j < NumberOfLocations; ++j){

		//Initialize list of arcs
		std::vector<PStep_Arc> CurrentInArcs;

		//Go over all start locations
		for (short int i = 0; i < NumberOfLocations; ++i){

			//Check wheter the end and start location are not the same
			if (j != i){

				//Initialize arc
				PStep_Arc CurrentArc;
				CurrentArc.StartNode = i;
				CurrentArc.EndNode = j;
				CurrentArc.Index = (i - 1)*i / 2 + j;

				//Add the new arc at the correct position in the list of current out arcs:

				//Obtain current cost
				double CurrentArcCost = TravelCosts[i][j];

				//Initialize indicator indicating whether the CurrentArc has not yet been added
				bool CurrentArcHasNotYetBeenAdded = true;
				//Initialize iterator
				auto Arc = CurrentInArcs.begin();;
				//Go over all current arcs in the list
				while (CurrentArcHasNotYetBeenAdded && Arc != CurrentInArcs.end()){

					//Check whether CurrentOutArc is cheaper than Arc
					if (CurrentArcCost < TravelCosts[Arc->StartNode][Arc->EndNode]){
						//The cost of CurrentOutArc is lower than of Arc,
						//add CurrentOutArc before Arc in the list.
						CurrentInArcs.insert(Arc, CurrentArc);
						//Update indicator
						CurrentArcHasNotYetBeenAdded = false;
					}
					else{
						//Update iterator
						++Arc;
					}

				}//end go over all current arcs

				//Check whether the current arc has not yet been added
				if (CurrentArcHasNotYetBeenAdded){
					//The current arc is the most expensive arc out of i, add it at the end of the list
					CurrentInArcs.push_back(CurrentArc);
				}

			}//end check whether end and start location are the same

		}//end go over all start locations

		//Add CurrentOutArcs
		InArcs.push_back(CurrentInArcs);

	}//end go over all end locations

}

//Initialize p-steps
void PStep_Algorithm::InitializePSteps(){

	//Make sure the current list is empty
	AlgorithmMembers.PSteps.clear();

	//Check type of initialization
	if (AlgorithmMembers.AlgorithmSettings.InitialPSteps == 1){
		AddAllReturnTrips();
	}
	else if (AlgorithmMembers.AlgorithmSettings.InitialPSteps == 2){
		AddGreedyInsertionRoutes();
	}
	else if (AlgorithmMembers.AlgorithmSettings.InitialPSteps == 3){
		AddGreedyInsertionRoutes();
		AddGreedyInsertionPSteps();
	}

}

//Add all return trips (0,i,0) to AlgorithmMembers.PSteps
void PStep_Algorithm::AddAllReturnTrips(){

	//Go over all customers
	for (short int i = 0; i < AlgorithmMembers.Instance.NumberOfCustomers; ++i){
		
		//Create p-step
		PStep NewPStep;
		NewPStep.Locations.push_back(0);
		NewPStep.Locations.push_back(i+1);
		NewPStep.Locations.push_back(0);
		NewPStep.StartingLoad = 0;
		NewPStep.TotalLoad = AlgorithmMembers.Instance.Demand[i];
		NewPStep.Cost = 2 * AlgorithmMembers.Instance.TravelCosts[0][i + 1];

		//Add p-step
		AlgorithmMembers.PSteps.push_back(NewPStep);

	}//end go over all customers

}

//Create a solution by applying a greedy insertion heuristic.
//In case the number of vehicles is fixed, this is satisfied by the solution.
//The solution is decomposed into p-steps which are added to AlgorithmMembers.PSteps.
//NOTE: in case the number of vehicles is fixed, we assume the number of customers
//      is at least the number of vehicles.
//		When no feasible initial routes are found, no p-steps are initialized.
void PStep_Algorithm::AddGreedyInsertionRoutes(){

	//Obtain number of customers and locations
	short int NumberOfCustomers = AlgorithmMembers.Instance.NumberOfCustomers;
	short int NumberOfLocations = NumberOfCustomers + 1;
	//Obtain vehicle capacity
	double VehicleCapacityPlusTolerance = AlgorithmMembers.Instance.VehicleCapacity+AlgorithmMembers.Tolerance;
	//Obtain costs
	double** TravelCosts = AlgorithmMembers.Instance.TravelCosts;
	//Obtain demand
	double* Demand = AlgorithmMembers.Instance.Demand;
	//Obtain indicator indicating whether the number of vehicles is not fixed
	bool NumberOfVehiclesIsNotFixed = !(AlgorithmMembers.Instance.NumberOfVehiclesIsFixed);
	
	//Initialize vector of customers that need to be scheduled:
	//The first customer has index 1.
	std::vector<short int> RemainingCustomers;
	//Go over all customer indices and add them to remaining customers
	for (short int i = 1; i < NumberOfLocations; ++i){
		RemainingCustomers.push_back(i);
	}

	//Initialize routes:
	//------------------------------------------------------------------------------------
	//Initialize Routes
	//Routes[i][j] is the j'th location on route i.
	std::list<std::list<short int>> Routes;
	std::list<double> Loads;

	//We will initialize InitialNumberOfRoutes.
	//Initialize InitialNumberOfRoutes
	short int InitialNumberOfRoutes=1;
	//Check whether the number of vehicles is fixed
	if (!NumberOfVehiclesIsNotFixed){
		//Set the InitialNumberOfRoutes to the fixed number of vehicles
		InitialNumberOfRoutes = AlgorithmMembers.Instance.FixedNumberOfVehicles;
	}

	//Initialize routes by selecting the closest customer to the depot for each of the initial routes.

	//Go over the initial routes
	for (short int Route = 0; Route < InitialNumberOfRoutes; ++Route){
		//Find the closest customer to the depot:

		//Initialize the first remaining customer as the current best
		std::vector<short int>::iterator CurrentBestCustomer = RemainingCustomers.begin();
		double CurrentBestCost = TravelCosts[0][RemainingCustomers[0]];
		//Go over all other remaining customers
		for (auto i = ++(RemainingCustomers.begin()); i != RemainingCustomers.end(); ++i){
			//Obtain customer and cost
			double CurrentCost = TravelCosts[0][*i];
			//Check whether the current customer is closer to the depot than the best
			if (CurrentCost<CurrentBestCost){
				//The current customer is closer, update the best
				CurrentBestCustomer = i;
				CurrentBestCost = CurrentCost;
			}
		}

		//Initialize a new route with the closest remaining customer.
		std::list<short int> NewRoute;
		NewRoute.push_back(*CurrentBestCustomer);
		//Add the new route to the list of Routes
		Routes.push_back(NewRoute);
		//Updat the route load
		Loads.push_back(Demand[*CurrentBestCustomer - 1]);
		//Remove the current customer from the list of remaining customers.
		RemainingCustomers.erase(CurrentBestCustomer);

	}//end go over all vehicles

	//Initialize indicator indicating whether the algorithm has not failed yet.
	bool AlgorithmDidNotFailYet = true;

	//Insert remaining customers
	//------------------------------------------------------------------------------------
	//Insert customers until no more remain or insertion fails.
	while (AlgorithmDidNotFailYet && RemainingCustomers.size() > 0){

		//Initialize best insertion
		std::vector<short int>::iterator CurrentBestCustomer;
		std::list<std::list<short int>>::iterator CurrentBestRoute;
		std::list<short int>::iterator CurrentBestPosition;
		std::list<double>::iterator CurrentBestLoad;
		double CurrentBestCost;
		double CurrentBestDemand;
		//Initialize indicator indicating whether a current best is not yet found
		bool CurrentBestNotYetFound = true;
		//Initialize indicator indicating whether the current best is a new route
		bool CurrentBestIsNewRoute=false;
		
		//Go over all remaining customers
		for (auto CurrentCustomer = RemainingCustomers.begin(); CurrentCustomer != RemainingCustomers.end(); ++CurrentCustomer){
			
			//Obtain current demand
			double CurrentDemand = Demand[(*CurrentCustomer) - 1];

			//Initialize vehicle load
			auto CurrentLoad = Loads.begin();
			//Go over all routes
			for (auto CurrentRoute = Routes.begin(); CurrentRoute != Routes.end(); ++CurrentRoute){
				//Check whether the current customer can be added to the route
				if ((*CurrentLoad) + CurrentDemand <= VehicleCapacityPlusTolerance){
					//Consider every position on the route for insertion:

					//Initialize previous location and insertion costs
					short int PreviousLocation = 0;
					double InsertionCost;

					//Go over all customers on the route
					for (auto NextCustomer = CurrentRoute->begin(); NextCustomer != CurrentRoute->end(); ++NextCustomer){
						//Compute costs of insertion
						InsertionCost = TravelCosts[PreviousLocation][*CurrentCustomer] + TravelCosts[*CurrentCustomer][*NextCustomer] - TravelCosts[PreviousLocation][*NextCustomer];
						//Check if the insertion costs is the best found so far
						if (CurrentBestNotYetFound || InsertionCost < CurrentBestCost){
							//The insertion cost is the best found so far, update the best
							CurrentBestCustomer = CurrentCustomer;
							CurrentBestRoute = CurrentRoute;
							CurrentBestPosition = NextCustomer;
							CurrentBestLoad = CurrentLoad;
							CurrentBestCost = InsertionCost;
							CurrentBestDemand = CurrentDemand;
							CurrentBestNotYetFound = false;
							CurrentBestIsNewRoute = false;
						}
						//Update the previous location
						PreviousLocation = *NextCustomer;
					}//end go over all customers on the route

					//Consider adding the current customer as the last location on the route:
					//Compute the cost of inserting the current customer at the end
					InsertionCost = TravelCosts[PreviousLocation][*CurrentCustomer] + TravelCosts[*CurrentCustomer][0] - TravelCosts[PreviousLocation][0];
					//Check if the insertion costs is the best found so far
					if (CurrentBestNotYetFound || InsertionCost < CurrentBestCost){
						//The insertion cost is the best found so far, update the best
						CurrentBestCustomer = CurrentCustomer;
						CurrentBestRoute = CurrentRoute;
						CurrentBestPosition = CurrentRoute->end();
						CurrentBestLoad = CurrentLoad;
						CurrentBestCost = InsertionCost;
						CurrentBestDemand = CurrentDemand;
						CurrentBestNotYetFound = false;
						CurrentBestIsNewRoute = false;
					}
				}//end check whether the demand of CurrentCustomer fits in the current route

				//Update current load
				++CurrentLoad;

			}//end go over all routes

			//Check whether the number of vehicles is fixed
			if (NumberOfVehiclesIsNotFixed){
				//The number of vehicles is not fixed, consider starting a new route with the current customer

				//Compute the cost of a new route
				double InsertionCost = TravelCosts[0][*CurrentCustomer] + TravelCosts[*CurrentCustomer][0];
				//Check if the insertion costs is the best found so far
				if (CurrentBestNotYetFound || InsertionCost < CurrentBestCost){

					//The insertion cost for the new route is the best found so far, update the best
					CurrentBestCustomer = CurrentCustomer;
					CurrentBestCost = InsertionCost;
					CurrentBestDemand = CurrentDemand;
					CurrentBestNotYetFound = false;
					CurrentBestIsNewRoute = true;

				}

			}

		}//end go over all remaining customers

		//Check whether a current best is found
		if (CurrentBestNotYetFound){
			//No current best is found, gthe algorithm failed to find a feasible insertion
			AlgorithmDidNotFailYet = false;
		}
		else if (CurrentBestIsNewRoute){
			//The current best is a new route, add this new route
			
			//Initialize a new route with the closest remaining customer.
			std::list<short int> NewRoute;
			NewRoute.push_back(*CurrentBestCustomer);
			//Add the new route to the list of Routes
			Routes.push_back(NewRoute);
			//Updat the route load
			Loads.push_back(Demand[*CurrentBestCustomer - 1]);
			//Update the remaining customers
			RemainingCustomers.erase(CurrentBestCustomer);

		}
		else{
			//A current best is found, perform the insert

			//Add the current best customer at the insert position of the appropriate route
			CurrentBestRoute->insert(CurrentBestPosition, *CurrentBestCustomer);
			//Update the corresponding load
			(*CurrentBestLoad) += CurrentBestDemand;
			//Update the remaining customers
			RemainingCustomers.erase(CurrentBestCustomer);

		}

	}//end while possible insert customers

	//Check whether a feasible solution has been found
	if (AlgorithmDidNotFailYet){
		//A feasible solution has been found, use this solution to initialize p-steps.

		//Obtain p
		short int p = AlgorithmMembers.AlgorithmSettings.p;

		//Go over all routes
		for (auto CurrentRoute = Routes.begin(); CurrentRoute != Routes.end(); ++CurrentRoute){

			//Add the depot to the route representation
			CurrentRoute->push_front(0);
			CurrentRoute->push_back(0);

			//Compute the number of p-steps of exactly p arcs
			int NumberOfFullPSteps= int(std::floor(((*CurrentRoute).size()-1)/p));
			//Compute the number of arcs ont the first p-step
			short int NumberOfArcsOnFirstPStep = (*CurrentRoute).size() - 1 - NumberOfFullPSteps*p;

			//Initialize Previous and Current Location
			auto PreviousLocation = CurrentRoute->begin();
			auto NextLocation = PreviousLocation;
			++NextLocation;
			//Initialize load
			double Load = 0;

			//Check whether a first p-step has to be made which does not traverse p arcs.
			if(NumberOfArcsOnFirstPStep>0){
				//Initialize p-step
				PStep NewPStep;
				NewPStep.Locations.push_back(0);
				NewPStep.Cost = 0;
				NewPStep.StartingLoad = 0;

				//Go over all arcs on the p-step
				for (short int i = 0; i < NumberOfArcsOnFirstPStep; ++i){
					//Update p-step
					NewPStep.Locations.push_back(*NextLocation);
					//If the next location is not the depot, update demand/
					if ((*NextLocation) > 0){
						Load += AlgorithmMembers.Instance.Demand[(*NextLocation)-1];
					}
					//Update costs
					NewPStep.Cost += AlgorithmMembers.Instance.TravelCosts[*PreviousLocation][*NextLocation];
					//Update locations
					PreviousLocation = NextLocation;
					++NextLocation;
				}

				//Update p-step
				NewPStep.TotalLoad = Load;

				//Add p-step
				AlgorithmMembers.PSteps.push_back(NewPStep);

			}

			//Create remaining p-steps
			for (int FullPStep = 0; FullPStep < NumberOfFullPSteps; ++FullPStep){

				//Initialize p-step
				PStep NewPStep;
				NewPStep.Locations.push_back(*PreviousLocation);
				NewPStep.Cost = 0;
				NewPStep.StartingLoad = Load;

				//Go over all arcs on the p-step
				for (short int i = 0; i < p; ++i){
					//Update p-step
					NewPStep.Locations.push_back(*NextLocation);
					//If the next location is not the depot, update demand/
					if ((*NextLocation) > 0){
						Load += AlgorithmMembers.Instance.Demand[(*NextLocation)-1];
					}
					//Update costs
					NewPStep.Cost += AlgorithmMembers.Instance.TravelCosts[*PreviousLocation][*NextLocation];
					//Update locations
					PreviousLocation = NextLocation;
					++NextLocation;
				}

				//Update p-step
				NewPStep.TotalLoad = Load;

				//Add p-step
				AlgorithmMembers.PSteps.push_back(NewPStep);

			}
			

		}//end go over all routes

	}//end check whether a feasible solution has been found

}

//Generate a p - step for every relevant pricing problem by applying a greedy insertion heuristic.
//For the pricing problems from the depot to the depot, add p-steps with a number
//of arcs from 2 to p with starting load zero.
//For the pricing problems from the depot to a customer, add p-steps with a number
//of arcs from 1 to p with starting load zero.
//For the pricing problems from a customer to a customer, add both a p-step with
//zero starting load and maximal starting load.
//For the pricing problems from a customer to the depot, add a p-step with maximal
//starting load.
void PStep_Algorithm::AddGreedyInsertionPSteps(){

	//Obtain support variables
	short int NumberOfCustomers = AlgorithmMembers.Instance.NumberOfCustomers;
	short int NumberOfLocations = NumberOfCustomers + 1;
	double** TravelCosts = AlgorithmMembers.Instance.TravelCosts;
	PStep_PricingProblem** (&PricingProblems) = AlgorithmMembers.PricingProblems;
	short int p = AlgorithmMembers.AlgorithmSettings.p;
	double* Demand = AlgorithmMembers.Instance.Demand;

	//Initialize list of all customers
	std::list<short int> CustomerLocations;
	//Go over all customers
	for (short int i = 1; i < NumberOfLocations; ++i){
		//Add customer
		CustomerLocations.push_back(i);
	}
	//Initialize updated list, which will be reused throughout
	std::list<short int> UpdatedCustomerLocations;

	//Initialize InsertionCosts array
	//------------------------------------------------------------------------------------

	//InsertionCosts[i][k][j] provides the cost increase of inserting location k in between i and j.
	double*** InsertionCosts=new double**[NumberOfLocations];
	//Go over all first locations
	for (short int i = 0; i < NumberOfLocations; ++i){
		//Initialize InsertionCosts for the first location i
		InsertionCosts[i] = new double*[NumberOfLocations];
		//Go over all insert locations
		for (short int k = 0; k < NumberOfLocations; ++k){
			//Initialize InsertionCosts for the insert location k
			InsertionCosts[i][k] = new double[NumberOfLocations];
			//Obtain costs from i to k
			double CostsFromiTok = TravelCosts[i][k];
			//Go over all final locations
			for (short int j = 0; j < NumberOfLocations; ++j){
				//Set costs
				InsertionCosts[i][k][j] = CostsFromiTok + TravelCosts[k][j] - TravelCosts[i][j];
			}//end go over all final locations
		}//end go over all insert locations
	}//end go over all first locations

	//Initialize p-steps
	//------------------------------------------------------------------------------------

	//Check whether the depot-depot pricing problem is relevant
	if (p > 1){
		//Obtain pricing problem
		PStep_PricingProblem &DepotDepotPricingProblem = PricingProblems[0][0];
		//Initialize p-step
		PStep NewPStep;
		NewPStep.Locations.push_back(0);
		NewPStep.Locations.push_back(0);
		NewPStep.StartingLoad = 0;
		NewPStep.TotalLoad = 0;
		NewPStep.Cost = 0;
		//Initialize a counter of the number of arcs on the p-step
		short int CurrentNumberOfArcs = 1;
		//Update list
		RemoveUnreachableLocationsAndSourceAndSink(CustomerLocations, UpdatedCustomerLocations, DepotDepotPricingProblem);
		//Initialize indicator indicating whether the greedy insert has not yet failed 
		bool NotYetFailed = true;
		//Go over all potential number of arcs of the p-step
		for (short int TargetNumberOfArcs = 2; TargetNumberOfArcs <= DepotDepotPricingProblem.MaxNumberOfArcs && NotYetFailed; ++TargetNumberOfArcs){
			NotYetFailed = GreedyInsertFillPStepAndAdd(NewPStep, CurrentNumberOfArcs, TargetNumberOfArcs, DepotDepotPricingProblem, InsertionCosts, UpdatedCustomerLocations);
		}//end go over all potential number of arcs of the p-step
	}//end check whether the depot-depot pricing problem is relevant

	//Obtain depot-customer pricing problems
	std::list<PStep_PricingProblem*>& RelevantDepotCustomerPP = AlgorithmMembers.RelevantDepotCustomerPricingProblems;
	//Go over all depot-customer pricing problems
	for (auto CurrentPP = RelevantDepotCustomerPP.begin(); CurrentPP != RelevantDepotCustomerPP.end(); ++CurrentPP){
		//Obtain current pricing problem
		PStep_PricingProblem &CurrentPricingProblem= **CurrentPP;
		//Obtain sink
		short int Sink = CurrentPricingProblem.Sink;
		//Initialize p-step
		PStep NewPStep;
		NewPStep.Locations.push_back(0);
		NewPStep.Locations.push_back(Sink);
		NewPStep.StartingLoad = 0;
		NewPStep.TotalLoad = Demand[Sink-1];
		NewPStep.Cost = TravelCosts[0][Sink];
		//Add the new p-step
		AlgorithmMembers.PSteps.push_back(NewPStep);
		//Initialize a counter of the number of arcs on the p-step
		short int CurrentNumberOfArcs = 1;
		//Update list
		RemoveUnreachableLocationsAndSourceAndSink(CustomerLocations, UpdatedCustomerLocations, CurrentPricingProblem);
		//Initialize indicator indicating whether the greedy insert has not yet failed 
		bool NotYetFailed = true;
		//Go over all potential number of arcs of the p-step
		for (short int TargetNumberOfArcs = 2; TargetNumberOfArcs <= CurrentPricingProblem.MaxNumberOfArcs && NotYetFailed; ++TargetNumberOfArcs){
			NotYetFailed = GreedyInsertFillPStepAndAdd(NewPStep, CurrentNumberOfArcs, TargetNumberOfArcs, CurrentPricingProblem, InsertionCosts, UpdatedCustomerLocations);
		}//end go over all potential number of arcs of the p-step
	}

	//Obtain customer-depot pricing problems
	std::list<PStep_PricingProblem*>& RelevantCustomerDepotPP = AlgorithmMembers.RelevantCustomerDepotPricingProblems;
	//Go over all customer-depot pricing problems
	for (auto CurrentPP = RelevantCustomerDepotPP.begin(); CurrentPP != RelevantCustomerDepotPP.end(); ++CurrentPP){
		//Obtain current pricing problem
		PStep_PricingProblem &CurrentPricingProblem = **CurrentPP;
		//Obtain source
		short int Source = CurrentPricingProblem.Source;
		//Initialize p-step
		PStep NewPStep;
		NewPStep.Locations.push_back(Source);
		NewPStep.Locations.push_back(0);
		NewPStep.StartingLoad = Demand[Source - 1];
		NewPStep.TotalLoad = NewPStep.StartingLoad;
		NewPStep.Cost = TravelCosts[Source][0];
		//Initialize a counter of the number of arcs on the p-step
		short int CurrentNumberOfArcs = 1;
		//Update list
		RemoveUnreachableLocationsAndSourceAndSink(CustomerLocations, UpdatedCustomerLocations, CurrentPricingProblem);
		//Check whether p=1
		if (p == 1){
			//Update the starting load to correspond to maximal load.
			NewPStep.StartingLoad = CurrentPricingProblem.Capacity;
			NewPStep.TotalLoad = CurrentPricingProblem.Capacity;
			//Add p-step
			AlgorithmMembers.PSteps.push_back(NewPStep);
		}
		else{
			//Construct p-step
			GreedyInsertFillPStepAndAdd(NewPStep, CurrentNumberOfArcs, p, CurrentPricingProblem, InsertionCosts, UpdatedCustomerLocations);
		}
	}

	//Obtain customer-customer pricing problems
	std::list<PStep_PricingProblem*>& RelevantCustomerCustomerPP = AlgorithmMembers.RelevantCustomerCustomerPricingProblems;
	//Go over all customer-customer pricing problems
	for (auto CurrentPP = RelevantCustomerCustomerPP.begin(); CurrentPP != RelevantCustomerCustomerPP.end(); ++CurrentPP){

		//Obtain current pricing problem
		PStep_PricingProblem &CurrentPricingProblem = **CurrentPP;
		//Obtain source and sink
		short int Source = CurrentPricingProblem.Source;
		short int Sink = CurrentPricingProblem.Sink;
		//Initialize p-step
		PStep NewPStep;
		NewPStep.Locations.push_back(Source);
		NewPStep.Locations.push_back(Sink);
		NewPStep.StartingLoad = Demand[Source - 1];
		NewPStep.TotalLoad = NewPStep.StartingLoad+Demand[Sink-1];
		NewPStep.Cost = TravelCosts[Source][Sink];
		//Initialize a counter of the number of arcs on the p-step
		short int CurrentNumberOfArcs = 1;
		//Update list
		RemoveUnreachableLocationsAndSourceAndSink(CustomerLocations, UpdatedCustomerLocations, CurrentPricingProblem);
		//Check whether p=1
		if (p == 1){
			//Add p-step with zero starting load
			AlgorithmMembers.PSteps.push_back(NewPStep);
			//Update the starting load to correspond to maximal load.
			NewPStep.StartingLoad = CurrentPricingProblem.Capacity - Demand[Sink - 1];
			NewPStep.TotalLoad = CurrentPricingProblem.Capacity;
			//Add p-step with maximal starting load
			AlgorithmMembers.PSteps.push_back(NewPStep);
		}
		else{
			//Construct p-step
			GreedyInsertFillPStepAndAdd(NewPStep, CurrentNumberOfArcs, p, CurrentPricingProblem, InsertionCosts, UpdatedCustomerLocations);
		}

		//Also consider the reverse pricing problem:
		//Obtain current pricing problem
		PStep_PricingProblem &ReversePricingProblem = PricingProblems[Sink][Source];
		//Obtain source and sink
		Source = Sink;
		Sink = ReversePricingProblem.Sink;
		//Initialize p-step
		PStep ReverseNewPStep;
		ReverseNewPStep.Locations.push_back(Source);
		ReverseNewPStep.Locations.push_back(Sink);
		ReverseNewPStep.StartingLoad = Demand[Source - 1];
		ReverseNewPStep.TotalLoad = ReverseNewPStep.StartingLoad + Demand[Sink - 1];
		ReverseNewPStep.Cost = TravelCosts[Source][Sink];
		//Initialize a counter of the number of arcs on the p-step
		CurrentNumberOfArcs = 1;
		//Update list
		RemoveUnreachableLocationsAndSourceAndSink(CustomerLocations, UpdatedCustomerLocations, ReversePricingProblem);
		//Check whether p=1
		if (p == 1){
			//Add p-step with zero starting load
			AlgorithmMembers.PSteps.push_back(ReverseNewPStep);
			//Update the starting load to correspond to maximal load.
			ReverseNewPStep.StartingLoad = CurrentPricingProblem.Capacity - Demand[Sink - 1];
			ReverseNewPStep.TotalLoad = CurrentPricingProblem.Capacity;
			//Add p-step with maximal starting load
			AlgorithmMembers.PSteps.push_back(ReverseNewPStep);
		}
		else{
			//Construct p-step
			GreedyInsertFillPStepAndAdd(ReverseNewPStep, CurrentNumberOfArcs, p, ReversePricingProblem, InsertionCosts, UpdatedCustomerLocations);
		}

	}

	//clean:
	for (short int i = 0; i < NumberOfLocations; ++i){
		for (short int k = 0; k < NumberOfLocations; ++k){
			delete[]InsertionCosts[i][k];
		}
		delete[]InsertionCosts[i];
	}
	delete[]InsertionCosts;

}

//Copy CustomerLocations to UpdatedCustomerLocations, while not copying the list
//that are not reachable or are the source or sink of the current pricing problem.
//It first empties UpdatedCustomerLocations.
//The first customer is coded as 1.
void PStep_Algorithm::RemoveUnreachableLocationsAndSourceAndSink(std::list<short int>& CustomerLocations, std::list<short int>& UpdatedCustomerLocations, PStep_PricingProblem& CurrentPP){

	//Empty
	UpdatedCustomerLocations.clear();

	//Obtain source, sink and reachability
	short int Source = CurrentPP.Source;
	short int Sink = CurrentPP.Sink;
	bool* Reachability = CurrentPP.Reachable;

	//Go over all customers in the list
	for (auto CurrentCustomer = CustomerLocations.begin(); CurrentCustomer != CustomerLocations.end(); ++CurrentCustomer){
		//Check whether the current customer should be removed
		if (*CurrentCustomer != Sink && *CurrentCustomer != Source && Reachability[(*CurrentCustomer)-1]){
			//The current customer should be copied
			UpdatedCustomerLocations.push_back(*CurrentCustomer);
		}
	}//end go over all customers

}

//Given the CurrentPStep consisting of CurrentNumberOfArcs, extend it by greedy inserts
//until it contains exactly TargetNumberOfArcs. If succesful, the output is true and 
//the identified p-step is added to AlgorithmMembers.PSteps. Otherwise the output is false.
//The pricing problem CurrentPP should correspond to the p-step. InsertionCosts[i][k][j]
//should be the costs of inserting customer k in between location i and j, where 0 is the depot.
//InsertableCustomers should contain all customer that might be inserted.
//NOTE: It is assumed that the members starting load and total load of the current p-step only include
//		the source demand and total demand respectively, such that the real starting load component
//		is 0. This will be set accordingly by this function.
bool PStep_Algorithm::GreedyInsertFillPStepAndAdd(PStep& CurrentPStep, short int& CurrentNumberOfArcs, short int TargetNumberOfArcs, PStep_PricingProblem& CurrentPP, double*** InsertionCosts, std::list<short int>& InsertableCustomers){

	//Initialize support variables
	double* Demand = CurrentPP.Demand;
	std::list<short int>& Locations = CurrentPStep.Locations;
	bool** ArcisRelevant = CurrentPP.ArcIsRelevant;
	double& TotalDemand = CurrentPStep.TotalLoad;
	double& Costs = CurrentPStep.Cost;
	double Capacity = CurrentPP.Capacity;
	double CapacityPlusTolerance = Capacity + AlgorithmMembers.Tolerance;

	//Initialize indicator indicating whether the greedy insert has not yet failed
	bool NotYetFailed = true;

	//Continue filling the p-step until the target number of locations is reached, or a faillure occurs
	while (CurrentNumberOfArcs < TargetNumberOfArcs && NotYetFailed){

		//Search for best insert:
		//------------------------------------------------------------------------------------

		//Initialize indicator indicating whether a best insert has been found
		bool BestNotYetFound = true;
		//Initialize best insertion costs
		double BestInsertionCosts;
		//Initialize best insertable location
		std::list<short int>::iterator BestInsertableCustomer;
		//Initialize best inser position (insert before this position)
		std::list<short int>::iterator BestInsertPosition;

		//Go over all locations that could be inserted
		for (std::list<short int>::iterator InsertCustomer = InsertableCustomers.begin(); InsertCustomer != InsertableCustomers.end(); ++InsertCustomer){
			
			//Obtain the insert customer
			short int k = *InsertCustomer;

			//Check whether the demand of customer k would fit the p-step
			if (TotalDemand + Demand[k - 1] <= CapacityPlusTolerance){

				//Initialize previous and next location
				std::list<short int>::iterator PreviousLocation = Locations.begin();
				short int i = *PreviousLocation;
				std::list<short int>::iterator NextLocation = PreviousLocation;
				++NextLocation;
				short int j;
				//Go over all potential insertion positions
				for (; NextLocation != Locations.end(); ++NextLocation){

					//Obtain next location index
					j = *NextLocation;

					//Check whether the current insertion is feasible and has lower costs than the current best
					if (ArcisRelevant[i][k] && ArcisRelevant[k][j]){
						//The insert is feasible

						//Check whether a best insert has already been found
						if (BestNotYetFound || InsertionCosts[i][k][j] < BestInsertionCosts){
							//The new insert is the current best
							BestNotYetFound = false;
							BestInsertionCosts = InsertionCosts[i][k][j];
							BestInsertableCustomer = InsertCustomer;
							BestInsertPosition = NextLocation;
						}

					}//end check whether inserying k in between i and j is feasible

					//Update previous location
					i = j;

				}//end go over all potential insertion positions

			}//end check whether the capacity allows including customer k

		}//end go over all locations that could be inserted

		//Execute best insert:
		//------------------------------------------------------------------------------------

		//Check whether a best was found
		if (!BestNotYetFound){
			//Execute the best insert:
			//Obtain best insertable customer
			short int BestCustomer = *BestInsertableCustomer;
			//Insert the customer
			Locations.insert(BestInsertPosition, BestCustomer);
			TotalDemand += Demand[BestCustomer - 1];
			Costs += BestInsertionCosts;
			//Remove the inserted customer from the list of insertable customers
			InsertableCustomers.erase(BestInsertableCustomer);
			//Update the number of included arcs
			++CurrentNumberOfArcs;
		}
		else{
			//No feasible insert was found, while the target number of arcs has not yet been reached.
			//We conclude that the greedy insert has failed
			NotYetFailed = false;
		}

	}//end continue filling


	//Update the p-step to account for a starting load and add:
	//------------------------------------------------------------------------------------
	
	//Check whether a propper p-step was identified
	if (NotYetFailed){
		
		//Check whether the source is the depot or the sink is not the depot
		if (CurrentPP.Source == 0 || CurrentPP.Sink != 0){
			//No update of the starting load is required.
			//Add p-step
			AlgorithmMembers.PSteps.push_back(CurrentPStep);
		}
		//Check whether the source is not the depot
		if (CurrentPP.Source!=0){
			//Update the starting load to correspond to maximal load.
			CurrentPStep.StartingLoad += Capacity - TotalDemand;
			CurrentPStep.TotalLoad = Capacity;
			//Add p-step
			AlgorithmMembers.PSteps.push_back(CurrentPStep);
		}

	}

	//Output
	return NotYetFailed;

}

//Initialize Neighbourhoods
//NOTE: it is assumed that OutArcs are sorted per location from cheap to expensive.
void PStep_Algorithm::InitializeNeighbourhoods(){

	//Check whether the algorithm makes use of neighbourhoods
	if (AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation){
		GenerateNeighbourhoods(AlgorithmMembers.AlgorithmSettings.NeighbourhoodSize);
	}

}

//Determine the n closest neighbours of each customer with respect to a certain
//distance matrix. The inward distance is considered.
//It sets AlgorithmMembers.Neighbourhoods and AlgorithmMembers.NotInNeighbourhood.
//The distances should be nonnegative.
void PStep_Algorithm::GenerateNeighbourhoods(short int n, double **Distance){

	//Obtain number of customers
	short int NumberOfCustomers = AlgorithmMembers.Instance.NumberOfCustomers;
	//Obtain neighbourhoods
	std::vector<std::vector<short int>>& Neighbourhoods = AlgorithmMembers.Neighbourhoods;

	//Empty current neighbourhoods
	Neighbourhoods.clear();
	if (AlgorithmMembers.NotInNeighbourhood != 0){
		for (short int i = 0; i < NumberOfCustomers; ++i){
			delete[] AlgorithmMembers.NotInNeighbourhood[i];
		}
		delete[]AlgorithmMembers.NotInNeighbourhood;
	}

	//Initialize template neighbourhood
	std::vector<short int> NeighbourhoodTemplate;
	//Initialize NotInNeighborhood
	AlgorithmMembers.NotInNeighbourhood = new bool*[NumberOfCustomers];
	//Obtain not in Neighborhoods
	bool** (&NotInNeighbourhood) = AlgorithmMembers.NotInNeighbourhood;
	//As a default, set NotInNeighbourhood as true for all customers
	for (short int i = 0; i < NumberOfCustomers; ++i){
		NotInNeighbourhood[i] = new bool[NumberOfCustomers];
		for (short int j = 0; j < NumberOfCustomers; ++j){
			NotInNeighbourhood[i][j] = true;
		}
	}

	//Go over all customers
	for (short int i = 0; i<NumberOfCustomers; i++){

		//Initialize the neighbourhood of the current customer.
		Neighbourhoods.push_back(NeighbourhoodTemplate);

		//Add locations to the neighbourhood while there are not yet n added or there are none left.
		while (short int(Neighbourhoods[i].size())<n && short int(Neighbourhoods[i].size())<NumberOfCustomers){

			//Find the location with the shortest distance to location i
			//Initialize the index of the nearest not selected location,
			//-1 indicates no location selected yet.
			int NearestLocation = -1;
			//Initialize the distance to the nearest location
			double NearestDistance = -1;
			//Go over all locations
			for (short int j = 0; j<NumberOfCustomers; j++){
				//Check whether j has not already been selected
				if (NotInNeighbourhood[i][j] && i!=j){
					//Check whether the current location is closer than the current closest
					if (Distance[j + 1][i + 1]<NearestDistance || NearestLocation == -1){
						//Selecte nearest location
						NearestLocation = j;
						NearestDistance = Distance[j + 1][i + 1];
					}
				}
			}
			//Add the seleceted location to the neighbourhood
			Neighbourhoods[i].push_back(NearestLocation);
			//Update NotInNeighbourhood
			NotInNeighbourhood[i][NearestLocation] = false;

		}//end while more customers should be added to the neighbourhood

	}//end go over all customers

}

//The following version of GenerateNeighbourhoods, uses AlgorithmMembers.InArcs, 
//under the assumption they are sorted.
void PStep_Algorithm::GenerateNeighbourhoods(short int n){

	//Obtain number of customers
	short int NumberOfCustomers = AlgorithmMembers.Instance.NumberOfCustomers;
	//Obtain InArcs
	std::vector<std::vector<PStep_Arc>>& InArcs = AlgorithmMembers.InArcs;
	//Obtain neighbourhoods
	std::vector<std::vector<short int >>& Neighbourhoods = AlgorithmMembers.Neighbourhoods;

	//Empty current neighbourhoods
	Neighbourhoods.clear();
	if (AlgorithmMembers.NotInNeighbourhood != 0){
		for (short int i = 0; i < NumberOfCustomers; ++i){
			delete[] AlgorithmMembers.NotInNeighbourhood[i];
		}
		delete[]AlgorithmMembers.NotInNeighbourhood;
	}

	//Initialize template neighbourhood
	std::vector<short int > NeighbourhoodTemplate;
	//Initialize NotInNeighborhood
	AlgorithmMembers.NotInNeighbourhood = new bool*[NumberOfCustomers];
	//Obtain not in Neighborhoods
	bool** (&NotInNeighbourhood) = AlgorithmMembers.NotInNeighbourhood;
	//As a default, set NotInNeighbourhood as true for all customers
	for (short int i = 0; i < NumberOfCustomers; ++i){
		NotInNeighbourhood[i] = new bool[NumberOfCustomers];
		for (short int j = 0; j < NumberOfCustomers; ++j){
			NotInNeighbourhood[i][j] = true;
		}
	}

	//Go over all customers
	for (short int i = 0; i<NumberOfCustomers; i++){

		//Initialize the neighbourhood of the current customer.
		Neighbourhoods.push_back(NeighbourhoodTemplate);

		//Initialize counter which keeps track of how many locations has been added to the current neighhbourhood.
		short int LocationsAdded = 0;
		//Add locations to the neighbourhood while there are not yet n added or there are none left.
		for (auto CurrentInArc = InArcs[i + 1].begin(); CurrentInArc != InArcs[i + 1].end() && LocationsAdded < n; ++CurrentInArc){

			//Obtain the start location
			short int StartLocation = CurrentInArc->StartNode;
			//Check whether the start location is not the depot
			if (StartLocation > 0){
				//The start location is not the depot.

				//Add the start location of CurrentInArc
				Neighbourhoods[i].push_back(StartLocation-1);
				//Update NotInNeighbourhood
				NotInNeighbourhood[i][StartLocation-1] = false;
				//Update LocationsAdded;
				++LocationsAdded;

			}

		}//end adding more customers to the neighbourhood

	}//end go over all customers

}

//Initialize pricing problems
void PStep_Algorithm::InitializePricingProblems(){

	//Obtain number of customers and locations, and tolerance
	short int NumberOfCustomers = AlgorithmMembers.Instance.NumberOfCustomers;
	short int NumberOfLocations = AlgorithmMembers.Instance.NumberOfCustomers + 1;
	double Tolerance = AlgorithmMembers.Tolerance;

	//Initialize pricing problems
	AlgorithmMembers.PricingProblems = new PStep_PricingProblem*[NumberOfLocations];
	//Obtain pricing problems
	PStep_PricingProblem** PricingProblems = AlgorithmMembers.PricingProblems;
	for (short int i = 0; i < NumberOfLocations; ++i){
		PricingProblems[i] = new PStep_PricingProblem[NumberOfLocations];
	}

	//Obtain setting detailing whether adjusted capacity is used
	bool UseAdjustedCapacity = AlgorithmMembers.AlgorithmSettings.UseAdjustedCapacity;

	//Sort demand for pre-processing purposes:
	//****************************************************************************

	//Obtain p
	short int p = AlgorithmMembers.AlgorithmSettings.p;
	//Obtain Demand
	double *Demand = AlgorithmMembers.Instance.Demand;
	//Obtain Vehicle capacity
	double VehicleCapacity = AlgorithmMembers.Instance.VehicleCapacity;
	double VehicleCapacityPlusTolerance = VehicleCapacity + Tolerance;

	//Initialize list with the 2p smallest demands ordered from small to large,
	//and corresponding locations.
	double* SmallestDemand = new double[2 * p];
	double* SmallestDemandLocations = new double[2 * p];

	//Sort demand:
	//------------

	//Create vector of candidates
	std::vector<short int> Candidates;
	for (short int i = 0; i < NumberOfCustomers; ++i){
		Candidates.push_back(i);
	}

	//Initialize the amount of smallest demands to be found
	short int NumberOfSmallest = 2 * p;
	if (NumberOfCustomers < 2 * p){
		NumberOfSmallest = NumberOfCustomers;
	}
	//Find smallest demands untill NumberOfSmallest are found
	for (short int Counter = 0; Counter < NumberOfSmallest; ++Counter){

		//Initialize current smallest demand and location
		short int CurrentSmallestDemandLocation = Candidates.front();
		double CurrentSmallestDemand = Demand[CurrentSmallestDemandLocation];
		size_t CurrentSmallestIndex = 0;

		//Go over all candidates and find the current smallest demand
		for (size_t Index = 1; Index < Candidates.size(); ++Index){

			//Obtain current candidate
			short int CurrentCandidate = Candidates[Index];

			//Check whether the current candidate has smaller demand
			if (Demand[CurrentCandidate] < CurrentSmallestDemand){
				//Set the current candidate as smallest
				CurrentSmallestDemandLocation = CurrentCandidate;
				CurrentSmallestDemand = Demand[CurrentCandidate];
				CurrentSmallestIndex = Index;
			}

		}

		//Process the found smallest demand
		SmallestDemand[Counter] = CurrentSmallestDemand;
		SmallestDemandLocations[Counter] = CurrentSmallestDemandLocation;
		Candidates.erase(Candidates.begin() + CurrentSmallestIndex);

	}//end find smallest demands until NumberOfSmallest are found

	//Fill the remainder of the smallest locations with bogus capacity that will never be used.
	for (short int Counter = NumberOfSmallest; Counter < 2 * p; ++Counter){
		SmallestDemand[Counter] = 0;
		SmallestDemandLocations[Counter] = -1;
	}

	//Create depot-depot pricing problem
	//****************************************************************************
	
	//Check whether p is not 1, beacuse for p=1, this pricing problem is not relevant
	if (p != 1){
		//Initialize depot-depot pricing problem
		PStep_PricingProblem& DepotDepotPricingProblem = PricingProblems[0][0];

		//Set number of locations and customers
		DepotDepotPricingProblem.NumberOfLocations = NumberOfLocations;
		DepotDepotPricingProblem.NumberOfCustomers = NumberOfCustomers;
		//Set source and sink
		DepotDepotPricingProblem.Source = 0;
		DepotDepotPricingProblem.Sink = 0;
		//Copy demand
		DepotDepotPricingProblem.Demand = new double[NumberOfCustomers];
		for (short int i = 0; i < NumberOfCustomers; ++i){
			DepotDepotPricingProblem.Demand[i] = Demand[i];
		}
		//Set indicator indicating that the Sink may be reached before p arcs are traversed.
		DepotDepotPricingProblem.AllowedToReachSinkEarly = true;
		//Initialize reduced costs
		DepotDepotPricingProblem.ReducedArcCosts = new double*[NumberOfLocations];
		for (short int i = 0; i < NumberOfLocations; ++i){
			DepotDepotPricingProblem.ReducedArcCosts[i] = new double[NumberOfLocations];
		}

		//Set the maximum number of arcs used in the depot-depot pricing problem:
		//Initialize the maximum number of arcs
		short int& MaxNumberOfArcs = DepotDepotPricingProblem.MaxNumberOfArcs;
		MaxNumberOfArcs = 1;
		//Check how many customers fit at most in one vehicle, for each customer add one arc:
		//Initialize load
		double CurrentLoad = 0;
		//Initialize customer index
		short int CustomerIndex = 0;
		//Compute the maximum number of customers on a p-step from depot to the depot
		short int MaxNumberOfCustomersDepotToDepot = p - 1;
		if (NumberOfCustomers < MaxNumberOfCustomersDepotToDepot){
			MaxNumberOfCustomersDepotToDepot = NumberOfCustomers;
		}
		//Add customers until the vehicle capacity is reached
		while (CustomerIndex < MaxNumberOfCustomersDepotToDepot && CurrentLoad + SmallestDemand[CustomerIndex] <= VehicleCapacityPlusTolerance){
			//Update
			CurrentLoad += SmallestDemand[CustomerIndex];
			++CustomerIndex;
			++MaxNumberOfArcs;
		}

		//Set capacity
		DepotDepotPricingProblem.Capacity = VehicleCapacity;
		
		//Initialize and set reachability
		DepotDepotPricingProblem.Reachable = new bool[NumberOfCustomers];
		for (short int i = 0; i < NumberOfCustomers; ++i){
			DepotDepotPricingProblem.Reachable[i] = true;
		}

		//Initialize and set ArcIsRelevant:
		DepotDepotPricingProblem.ArcIsRelevant = new bool*[NumberOfLocations];
		for (short int i = 0; i < NumberOfLocations; ++i){
			DepotDepotPricingProblem.ArcIsRelevant[i] = new bool[NumberOfLocations];
			for (short int j = 0; j < NumberOfLocations; ++j){
				DepotDepotPricingProblem.ArcIsRelevant[i][j] = false;
			}
		}

		//OLD CODE:
		//This does not work nicely for ng path relaxation...
		/*
		//Initialize and set LayeredCapacity
		DepotDepotPricingProblem.LayeredCapacity = new double[MaxNumberOfArcs + 1];
		DepotDepotPricingProblem.LayeredCapacityPlusTolerance = new double[MaxNumberOfArcs + 1];
		double* (&LayeredCapacity) = DepotDepotPricingProblem.LayeredCapacity;
		double* (&LayeredCapacityPlusTolerance) = DepotDepotPricingProblem.LayeredCapacityPlusTolerance;
		LayeredCapacity[0] = 0;
		LayeredCapacityPlusTolerance[0] = 0;
		for (short int i = 1; i <= MaxNumberOfArcs; ++i){
			LayeredCapacity[i] = VehicleCapacity;
			LayeredCapacityPlusTolerance[i] = VehicleCapacityPlusTolerance;
		}
		*/

		//Set OutArcs:
		//Obtain out arcs
		std::vector<std::vector<PStep_Arc>>& OutArcs = AlgorithmMembers.OutArcs;
		std::list<std::list<PStep_Arc>>& PPOutArcs=DepotDepotPricingProblem.OutArcs;
		//Copy the arcs out of the depot:
		//Obtain the arcs out of the depot
		auto CurrentOutArcs = OutArcs.begin();
		//Initialize list
		PPOutArcs.emplace_front();
		std::list<PStep_Arc>& PPArcsOutOfDepot=PPOutArcs.front();
		//Go over all arcs out of the depot
		for (auto CurrentArc =CurrentOutArcs->begin(); CurrentArc != CurrentOutArcs->end(); ++CurrentArc){
			//Copy the current arc out of the depot
			PPArcsOutOfDepot.push_back(*CurrentArc);
			DepotDepotPricingProblem.ArcIsRelevant[0][CurrentArc->EndNode] = true;
		}
		//Update CurrentOutArcs to point to the out arcs out of the first customer
		++CurrentOutArcs;
		//Go over all customers, and all arcs out of those customers.
		for (; CurrentOutArcs != OutArcs.end(); ++CurrentOutArcs){
			//Initialize new list
			PPOutArcs.emplace_back();
			std::list<PStep_Arc>& CurrentPPOutArcs = PPOutArcs.back();
			//Go over all arcs out of the depot
			for (auto CurrentArc = CurrentOutArcs->begin(); CurrentArc != CurrentOutArcs->end(); ++CurrentArc){
				//Check whether the arc is into the depot or the demand of the two nodes does not exceed the capacity
				if (CurrentArc->EndNode == 0 || (p > 2 && Demand[(CurrentArc->StartNode) - 1] + Demand[(CurrentArc->EndNode) - 1] <= VehicleCapacityPlusTolerance)){
					//Copy the current arc 
					CurrentPPOutArcs.push_back(*CurrentArc);
					DepotDepotPricingProblem.ArcIsRelevant[CurrentArc->StartNode][CurrentArc->EndNode] = true;
				}
			}
		}

		//Set InArcs:
		//Obtain in arcs
		std::vector<std::vector<PStep_Arc>>& InArcs = AlgorithmMembers.InArcs;
		std::list<std::list<PStep_Arc>>& PPInArcs = DepotDepotPricingProblem.InArcs;
		//Copy the arcs into the depot:
		//Obtain the arcs into the depot
		auto CurrentInArcs = InArcs.begin();
		//Initialize list
		PPInArcs.emplace_front();
		std::list<PStep_Arc>& PPArcsIntoDepot = PPInArcs.front();
		//Go over all arcs into the depot
		for (auto CurrentArc = CurrentInArcs->begin(); CurrentArc != CurrentInArcs->end(); ++CurrentArc){
			//Copy the current arc into the depot
			PPArcsIntoDepot.push_back(*CurrentArc);
		}
		//Update CurrentInArcs to point to the arcs into the first customer
		++CurrentInArcs;
		//Go over all customers, and all arcs into those customers.
		for (; CurrentInArcs != InArcs.end(); ++CurrentInArcs){
			//Initialize new list
			PPInArcs.emplace_back();
			std::list<PStep_Arc>& CurrentPPInArcs = PPInArcs.back();
			//Go over all arcs into depot
			for (auto CurrentArc = CurrentInArcs->begin(); CurrentArc != CurrentInArcs->end(); ++CurrentArc){
				//Check whether the arc is out of the depot or the demand of the two nodes does not exceed the capacity
				if (CurrentArc->StartNode == 0 || (p > 2 && Demand[(CurrentArc->StartNode) - 1] + Demand[(CurrentArc->EndNode) - 1] <= VehicleCapacityPlusTolerance)){
					//Copy the current arc
					CurrentPPInArcs.push_back(*CurrentArc);
				}
			}
		}

		//Initialize exact column generator members (mainly label containers):
		if (AlgorithmMembers.TypeOfExactColumnGenerator == 1){
			//Initialize LayeredLabels:

			//Obtain LayeredLabels
			DepotDepotPricingProblem.LayeredLabels = new std::list<Label>*[MaxNumberOfArcs];
			std::list<Label>** LayeredLabels = DepotDepotPricingProblem.LayeredLabels;
			//Initialize source label
			LayeredLabels[0] = new std::list<Label>[1];
			//Initalize labels for layers 1 to MaxNumberOfArcs
			for (short int Layer = 1; Layer < MaxNumberOfArcs; ++Layer){
				LayeredLabels[Layer] = new std::list<Label>[NumberOfCustomers];
			}
		}
		else if (AlgorithmMembers.TypeOfExactColumnGenerator==2){
			//Initialize DoubleLayeredLabels:
			//Obtain integer capacity
			int IntegerCapacity = int(VehicleCapacity);
			//Initialize DoubleLayeredLabels
			DepotDepotPricingProblem.DoubleLayeredLabels = new std::list<Label>**[IntegerCapacity+1];
			std::list<Label>*** DoubleLayeredLabels = DepotDepotPricingProblem.DoubleLayeredLabels;
			//Initialize source label
			DoubleLayeredLabels[0] = new std::list<Label>*[1];
			DoubleLayeredLabels[0][0] = new std::list<Label>[1];
			//Initialize layers to 0 for all loads
			for (int q = 1; q <= IntegerCapacity; ++q){
				DoubleLayeredLabels[q] = 0;
			}
		}
		else if (AlgorithmMembers.TypeOfExactColumnGenerator == 3){
			
			//Compute the maximum number of forward and backward labels
			short int MaxNumberOfBackwardArcs = (MaxNumberOfArcs - 1) / 2;
			short int MaxNumberOfForwardArcs = MaxNumberOfArcs - 1 - MaxNumberOfBackwardArcs;
			//Set maximum number of forward and backward labels.
			DepotDepotPricingProblem.MaxNumberOfForwardArcs = MaxNumberOfForwardArcs;
			DepotDepotPricingProblem.MaxNumberOfBackwardArcs = MaxNumberOfBackwardArcs;
			//Set the number of forward and backward layers
			DepotDepotPricingProblem.NumberOfForwardLayers = MaxNumberOfForwardArcs + 1;
			DepotDepotPricingProblem.NumberOfBackwardLayers = MaxNumberOfBackwardArcs + 1;

			//Initialize bi-directional labels:
			
			//Obtain LayeredForwardLabels
			DepotDepotPricingProblem.LayeredForwardLabels = new std::list<ForwardLabel>*[MaxNumberOfForwardArcs+1];
			std::list<ForwardLabel>** LayeredForwardLabels = DepotDepotPricingProblem.LayeredForwardLabels;
			//Initialize source label
			LayeredForwardLabels[0] = new std::list<ForwardLabel>[1];
			//Initalize labels for layers 1 to MaxNumberOfForwardArcs
			for (short int Layer = 1; Layer <= MaxNumberOfForwardArcs; ++Layer){
				LayeredForwardLabels[Layer] = new std::list<ForwardLabel>[NumberOfCustomers];
			}
			
			//Obtain LayeredBackwardLabels
			DepotDepotPricingProblem.LayeredBackwardLabels = new std::list<BackwardLabel>*[MaxNumberOfBackwardArcs+1];
			std::list<BackwardLabel>** LayeredBackwardLabels = DepotDepotPricingProblem.LayeredBackwardLabels;
			//Initialize source label
			LayeredBackwardLabels[0] = new std::list<BackwardLabel>[1];
			//Initalize labels for layers 1 to MaxNumberOfBackwardArcs
			for (short int Layer = 1; Layer <= MaxNumberOfBackwardArcs; ++Layer){
				LayeredBackwardLabels[Layer] = new std::list<BackwardLabel>[NumberOfCustomers];
			}
			
		}
		else if (AlgorithmMembers.TypeOfExactColumnGenerator == 4){

			//Compute the maximum number of forward and backward labels
			short int MaxNumberOfBackwardArcs = (MaxNumberOfArcs - 1) / 2;
			short int MaxNumberOfForwardArcs = MaxNumberOfArcs - 1 - MaxNumberOfBackwardArcs;
			//Set maximum number of forward and backward labels.
			DepotDepotPricingProblem.MaxNumberOfForwardArcs = MaxNumberOfForwardArcs;
			DepotDepotPricingProblem.MaxNumberOfBackwardArcs = MaxNumberOfBackwardArcs;
			//Set the number of forward and backward layers
			DepotDepotPricingProblem.NumberOfForwardLayers = MaxNumberOfArcs;
			DepotDepotPricingProblem.NumberOfBackwardLayers = MaxNumberOfArcs;

			//Initialize bi-directional labels and completion bounds:

			//Obtain LayeredForwardLabels
			DepotDepotPricingProblem.LayeredForwardLabels = new std::list<ForwardLabel>*[MaxNumberOfArcs];
			std::list<ForwardLabel>** LayeredForwardLabels = DepotDepotPricingProblem.LayeredForwardLabels;
			//Initialize source label
			LayeredForwardLabels[0] = new std::list<ForwardLabel>[1];
			//Obtain ForwardCompletionBounds
			DepotDepotPricingProblem.ForwardCompletionBounds = new double*[MaxNumberOfArcs];
			double** ForwardCompletionBounds = DepotDepotPricingProblem.ForwardCompletionBounds;
			//Initalize labels and bounds for layers 1 to MaxNumberOfArcs
			for (short int Layer = 1; Layer < MaxNumberOfArcs; ++Layer){
				LayeredForwardLabels[Layer] = new std::list<ForwardLabel>[NumberOfCustomers];
				ForwardCompletionBounds[Layer] = new double[NumberOfCustomers];
			}

			//Obtain LayeredBackwardLabels
			DepotDepotPricingProblem.LayeredBackwardLabels = new std::list<BackwardLabel>*[MaxNumberOfArcs];
			std::list<BackwardLabel>** LayeredBackwardLabels = DepotDepotPricingProblem.LayeredBackwardLabels;
			//Initialize source label
			LayeredBackwardLabels[0] = new std::list<BackwardLabel>[1];
			//Obtain BackwardCompletionBounds
			DepotDepotPricingProblem.BackwardCompletionBounds = new double*[MaxNumberOfArcs];
			double** BackwardCompletionBounds = DepotDepotPricingProblem.BackwardCompletionBounds;
			//Initalize labels and bounds for layers 1 to MaxNumberOfArcs
			for (short int Layer = 1; Layer < MaxNumberOfArcs; ++Layer){
				LayeredBackwardLabels[Layer] = new std::list<BackwardLabel>[NumberOfCustomers];
				BackwardCompletionBounds[Layer] = new double[NumberOfCustomers];
			}

		}

	}//end check whether p is not 1

	//Create depot-customer pricing problem
	//****************************************************************************

	//Go over all customers
	for (short int j = 1; j < NumberOfLocations; ++j){
		
		//Initialize support
		short int jMinusOne = j - 1;

		//Determine the capacity:
		//Initialize excluded customers
		std::vector<short int> ExcludedCustomers;
		ExcludedCustomers.push_back(jMinusOne);
		//Initialize Index indicating which customer is the first next smallest customer to consider
		short int SmallestIndex = 0;

		//Compute the remaining capacity after visiting p-1 customers on the on a potential next p-step starting at j and ending at the depot.
		double RemainingCapacity = ComputeCapacity(VehicleCapacity, SmallestDemand, SmallestDemandLocations, p - 1, ExcludedCustomers, NumberOfCustomers, SmallestIndex);

		//Initialize current capacity, it will represent adjusted or vehicle capacity depending on UseAdjustedCapacity.
		double CurrentCapacity;
		//Check whether adjusted capacity is used
		if (UseAdjustedCapacity){
			//Compute adjusted capacity
			CurrentCapacity = RemainingCapacity;
		}
		else{
			//Set the current capacity to vehicle capacity
			CurrentCapacity = VehicleCapacity;
		}

		//Obtain demand perturbed by tolerance:
		double PerturbedDemand = Demand[jMinusOne] - AlgorithmMembers.Tolerance;

		//Check whether the pricing problem is relevant
		if (PerturbedDemand <= RemainingCapacity && p < NumberOfLocations) {
			//The pricing problem to create a p-step from the depot to customer j is relevant.

			//Add the pricing problem to the list of relevant pricing problems
			AlgorithmMembers.RelevantDepotCustomerPricingProblems.push_back(&(PricingProblems[0][j]));

			//Obtain the current pricing problem
			PStep_PricingProblem& CurrentPricingProblem = PricingProblems[0][j];

			//Create depot-customer pricing problem:
			//--------------------------------------

			//Set number of locations and customers
			CurrentPricingProblem.NumberOfLocations = NumberOfLocations;
			CurrentPricingProblem.NumberOfCustomers = NumberOfCustomers;
			//Set source and sink
			CurrentPricingProblem.Source = 0;
			CurrentPricingProblem.Sink = j;
			//Copy demand
			CurrentPricingProblem.Demand = new double[NumberOfCustomers];
			for (short int i = 0; i < NumberOfCustomers; ++i){
				CurrentPricingProblem.Demand[i] = Demand[i];
			}
			//Set indicator indicating that the Sink may be reached before p arcs are traversed.
			CurrentPricingProblem.AllowedToReachSinkEarly = true;
			//Initialize reduced costs
			CurrentPricingProblem.ReducedArcCosts = new double*[NumberOfLocations];
			for (short int i = 0; i < NumberOfLocations; ++i){
				CurrentPricingProblem.ReducedArcCosts[i] = new double[NumberOfLocations];
			}

			//Determine how many arcs can be traversed on the depot customer p-step:
			//Initialize the maximum number of arcs
			short int& MaxNumberOfArcs = CurrentPricingProblem.MaxNumberOfArcs;
			MaxNumberOfArcs = 0;
			//Initialize temporary capacity
			double TempCapacity = CurrentCapacity;
			//Go over all potential additional visits.
			while (MaxNumberOfArcs < p && PerturbedDemand <= TempCapacity){
				//Update visits
				++MaxNumberOfArcs;
				//Update capacity
				TempCapacity = ComputeCapacity(TempCapacity, SmallestDemand, SmallestDemandLocations, 1, ExcludedCustomers, NumberOfCustomers, SmallestIndex);
			}

			//Set capacity
			CurrentPricingProblem.Capacity = CurrentCapacity;

			//Initialize and set reachability
			CurrentPricingProblem.Reachable = new bool[NumberOfCustomers];
			bool* Reachable = CurrentPricingProblem.Reachable;
			//Check whether MaxNumberOfArcs=1
			if (MaxNumberOfArcs == 1){
				//Go over all customers and set reachability
				for (short int i = 0; i < NumberOfCustomers; ++i){
					Reachable[i] = false;
				}
				Reachable[jMinusOne] = true;
			}
			else{
				//Go over all customers and set reachability
				for (short int i = 0; i < NumberOfCustomers; ++i){
					//Check whether i=j
					if (i == jMinusOne){
						//The sink is reachable
						Reachable[jMinusOne] = true;
					}
					else{
						//Check whether a 0-j p-step exists in which also i is visited.

						//Initialize current excluded customers
						std::vector<short int> CurrentExcludedCustomers = ExcludedCustomers;
						CurrentExcludedCustomers.push_back(i);
						//Compute capacity taking into account that i is visited
						//Check whether adjusted capacity is used
						if (UseAdjustedCapacity){
							//Compute adjusted capacity
							TempCapacity = ComputeCapacity(VehicleCapacity, SmallestDemand, SmallestDemandLocations, p - 1, CurrentExcludedCustomers, NumberOfCustomers);
						}
						else{
							TempCapacity = VehicleCapacity;
						}
						//Check whether the demand of j and k is less than the vehicle capacity
						if (PerturbedDemand + Demand[i] <= TempCapacity){
							//Customer i could be visited
							Reachable[i] = true;
						}
						else{
							//Customer i cannot be reached
							Reachable[i] = false;
						}
					}
				}
			}

			//Initialize ArcIsRelevant
			CurrentPricingProblem.ArcIsRelevant = new bool*[NumberOfLocations];
			bool** ArcIsRelevant = CurrentPricingProblem.ArcIsRelevant;
			for (short int k = 0; k < NumberOfLocations; ++k){
				ArcIsRelevant[k] = new bool[NumberOfLocations];
				for (short int l = 0; l < NumberOfLocations; ++l){
					ArcIsRelevant[k][l] = false;
				}
			}

			//OLD CODE:
			//This does not work nicely for ng path relaxation...
			/*
			//Initialize and set LayeredCapacity
			CurrentPricingProblem.LayeredCapacity = new double[MaxNumberOfArcs + 1];
			CurrentPricingProblem.LayeredCapacityPlusTolerance = new double[MaxNumberOfArcs + 1];
			double* (&LayeredCapacity) = CurrentPricingProblem.LayeredCapacity;
			double* (&LayeredCapacityPlusTolerance) = CurrentPricingProblem.LayeredCapacityPlusTolerance;
			LayeredCapacity[0] = 0;
			LayeredCapacityPlusTolerance[0] = 0;
			double CurrentCapacityPlusTolerance = CurrentCapacity + Tolerance;
			for (short int i = 1; i <= MaxNumberOfArcs; ++i){
				LayeredCapacity[i] = CurrentCapacity;
				LayeredCapacityPlusTolerance[i] = CurrentCapacityPlusTolerance;
			}
			*/

			//Set OutArcs and ArcIsRelevant:
			//------------------------------

			//Obtain out arcs
			std::vector<std::vector<PStep_Arc>>& OutArcs = AlgorithmMembers.OutArcs;
			std::list<std::list<PStep_Arc>>& PPOutArcs = CurrentPricingProblem.OutArcs;

			//Set OutArcs and ArcIsRelevant for the arcs out of the depot
			auto CurrentOutArcs = OutArcs.begin();
			//Initialize arcs out of the depot
			PPOutArcs.emplace_front();
			//Obtain arcs out of the depot
			std::list<PStep_Arc>& PPArcsOutOfDepot = PPOutArcs.front();
			//Go over all arcs out of the depot
			for (auto CurrentArc = CurrentOutArcs->begin(); CurrentArc != CurrentOutArcs->end(); ++CurrentArc){
				//Obtain end node
				short int EndNode = CurrentArc->EndNode;
				//Check whether the end node can be reached
				if (Reachable[EndNode - 1]){
					//The arc (0,EndNode) can be used, add it and make it relevant
					PPArcsOutOfDepot.push_back(*CurrentArc);
					ArcIsRelevant[0][EndNode] = true;
				}
			}//end go over all arcs out of the depot

			//Go over all locations as start nodes of the arc:
			//Initialize from arcs out of location 1.
			short int StartNode = 1;
			++CurrentOutArcs;
			//Check whether MaxNumberOfArcs=1
			if (MaxNumberOfArcs == 1){
				//Only the arc (0,i) is feasible, add empty lists of outarcs
				for (; CurrentOutArcs != OutArcs.end(); ++CurrentOutArcs){
					//Initialize new list
					PPOutArcs.emplace_back();
				}
			}
			else{
				//MaxNumberOfArcs>1
				for (; CurrentOutArcs != OutArcs.end(); ++CurrentOutArcs){

					//Initialize new list
					PPOutArcs.emplace_back();

					//Check whether the current start node is reachable, and it is not the sink
					if (Reachable[StartNode - 1] && StartNode != j){
						//StartNode is reachable, set OutArcs.

						//Obtain new list
						std::list<PStep_Arc>& CurrentPPOutArcs = PPOutArcs.back();

						//Check whether MaxNumberOfArcs=2
						if (MaxNumberOfArcs == 2){
							//Only the arc (StartNode,j) could be used.

							//Go over all arcs out of StartNode
							for (auto CurrentArc = CurrentOutArcs->begin(); CurrentArc != CurrentOutArcs->end(); ++CurrentArc){

								//Obtain end node
								short int EndNode = CurrentArc->EndNode;
								//Check whether the end node is the sink, in which case the arc can immediately be added
								if (EndNode == j){
									//The arc (0,EndNode) can be used, add it and make it relevant
									CurrentPPOutArcs.push_back(*CurrentArc);
									ArcIsRelevant[StartNode][EndNode] = true;
								}
							}
						}
						else{
							//MaxNumberOfArcs>2

							//Initialize current excluded customers
							std::vector<short int> CurrentExcludedCustomers = ExcludedCustomers;
							CurrentExcludedCustomers.push_back(StartNode - 1);
							//Initialize the perturbed demand of j and the startNode
							double PerturbedDemandIncludingStartNode = PerturbedDemand + Demand[StartNode - 1];

							//Go over all arcs out of StartNode
							for (auto CurrentArc = CurrentOutArcs->begin(); CurrentArc != CurrentOutArcs->end(); ++CurrentArc){

								//Obtain end node
								short int EndNode = CurrentArc->EndNode;
								//Check whether the end node is the sink, in which case the arc can immediately be added
								if (EndNode == j){
									//The arc (0,EndNode) can be used, add it and make it relevant
									CurrentPPOutArcs.push_back(*CurrentArc);
									ArcIsRelevant[StartNode][EndNode] = true;
								}
								//Otherwise check whether the end node is not the depot
								else if (EndNode != 0){
									//StartNode, EndNode and j are all different, check whether they fit on one p-step.

									//Initialize current excluded customers
									std::vector<short int> CurrentExcludedCustomersForArc = CurrentExcludedCustomers;
									CurrentExcludedCustomersForArc.push_back(EndNode - 1);
									//Check whether adjusted capacity is used
									if (UseAdjustedCapacity){
										//Compute adjusted capacity
										TempCapacity = ComputeCapacity(VehicleCapacity, SmallestDemand, SmallestDemandLocations, p - 1, CurrentExcludedCustomersForArc, NumberOfCustomers);
									}
									else{
										TempCapacity = VehicleCapacity;
									}
									//Check capacity
									if (PerturbedDemandIncludingStartNode + Demand[EndNode - 1] <= TempCapacity){
										//The arc (StartNode,EndNode) can be used, add it and make it relevant
										CurrentPPOutArcs.push_back(*CurrentArc);
										ArcIsRelevant[StartNode][EndNode] = true;
									}

								}

							}//end go over all arcs out of StartNode
						}//end else MaxNumberOfArcs>2

					}//end check whether the current start node is reachable

					//Update StartNode
					++StartNode;

				}//end go over all locations as start node of the arc
			}//end else MaxNumberOfArcs>1

			//Set InArcs:
			//------------------------------

			//Obtain in arcs
			std::vector<std::vector<PStep_Arc>>& InArcs = AlgorithmMembers.InArcs;
			std::list<std::list<PStep_Arc>>& PPInArcs = CurrentPricingProblem.InArcs;

			//Set InArcs and ArcIsRelevant for the arcs into the depot
			auto CurrentInArcs = InArcs.begin();
			//Initialize arcs into the depot and leave it empty
			PPInArcs.emplace_front();

			//Go over all locations as end nodes of the arc:
			//Initialize from arcs into location 1.
			short int EndNode = 1;
			++CurrentInArcs;
			for (; CurrentInArcs != InArcs.end(); ++CurrentInArcs){

				//Initialize new list
				PPInArcs.emplace_back();

				//Check whether the end node is the sink
				if (EndNode == j){

					//Obtain new list
					std::list<PStep_Arc>& CurrentPPInArcs = PPInArcs.back();

					//Check whether MaxNumberOfArcs=1
					if (MaxNumberOfArcs == 1){
						//Only the arc (0,j) is feasible
						
						//Go over all InArcs
						for (auto CurrentArc = CurrentInArcs->begin(); CurrentArc != CurrentInArcs->end(); ++CurrentArc){
							//Check whether the startnode is the depot
							if (CurrentArc->StartNode == 0){
								//Add the in arc
								CurrentPPInArcs.push_back(*CurrentArc);
							}
						}
					}
					else{
						//MaxNumberOfArcs>1

						//The end node is the sink, all arcs into j from a reachable node and the depot are feasible.

						//Go over all InArcs
						for (auto CurrentArc = CurrentInArcs->begin(); CurrentArc != CurrentInArcs->end(); ++CurrentArc){
							//Obtain the start node
							StartNode = CurrentArc->StartNode;
							//Check whether the startnode is the depot or reachable
							if (StartNode == 0 || Reachable[StartNode - 1]){
								//Add the in arc
								CurrentPPInArcs.push_back(*CurrentArc);
							}
						}

					}

				}
				//Check whether the current end node, which is not the sink, is reachable
				else if (Reachable[EndNode - 1] && MaxNumberOfArcs > 1){
					//EndNode is reachable and not the sink, set InArcs.

					//Obtain new list
					std::list<PStep_Arc>& CurrentPPInArcs = PPInArcs.back();

					//Check whether MaxNumberOfArcs>2
					if (MaxNumberOfArcs == 2){
						//Only the arc (0,EndNode) is feasible
						//Go over all arcs into EndNode
						for (auto CurrentArc = CurrentInArcs->begin(); CurrentArc != CurrentInArcs->end(); ++CurrentArc){
							//Obtain start node
							StartNode = CurrentArc->StartNode;
							//Check whether the start node is the depot, in which case the arc can immediately be added
							if (StartNode == 0){
								//The arc (0,EndNode) can be used, add it.
								CurrentPPInArcs.push_back(*CurrentArc);
							}
						}
					}
					else{
						//MaxNumberOfArcs>2

						//Initialize current excluded customers
						std::vector<short int> CurrentExcludedCustomers = ExcludedCustomers;
						CurrentExcludedCustomers.push_back(EndNode - 1);
						//Initialize the perturbed demand of j and the end node
						double PerturbedDemandIncludingEndNode = PerturbedDemand + Demand[EndNode - 1];

						//Go over all arcs into EndNode
						for (auto CurrentArc = CurrentInArcs->begin(); CurrentArc != CurrentInArcs->end(); ++CurrentArc){

							//Obtain start node
							StartNode = CurrentArc->StartNode;
							//Check whether the start node is the depot, in which case the arc can immediately be added
							if (StartNode == 0){
								//The arc (0,EndNode) can be used, add it.
								CurrentPPInArcs.push_back(*CurrentArc);
							}
							//Otherwise check whether the start node is not the sink
							else if (StartNode != j){
								//StartNode, EndNode and j are all different, check whether they fit on one p-step.

								//Initialize current excluded customers
								std::vector<short int> CurrentExcludedCustomersForArc = CurrentExcludedCustomers;
								CurrentExcludedCustomersForArc.push_back(StartNode - 1);
								//Check whether adjusted capacity is used
								if (UseAdjustedCapacity){
									//Compute adjusted capacity
									TempCapacity = ComputeCapacity(VehicleCapacity, SmallestDemand, SmallestDemandLocations, p - 1, CurrentExcludedCustomersForArc, NumberOfCustomers);
								}
								else{
									TempCapacity = VehicleCapacity;
								}
								//Check capacity
								if (PerturbedDemandIncludingEndNode + Demand[StartNode - 1] <= TempCapacity){
									//The arc (StartNode,EndNode) can be used, add it.
									CurrentPPInArcs.push_back(*CurrentArc);
								}

							}

						}//end go over all arcs into EndNode
					}//end else MaxNumberOfArcs>2

				}//end check whether the current end node is reachable and MaxNumberOfArcs>1

				//Update EndNode
				++EndNode;

			}//end go over all locations as end node of the arc

			//Initialize Labels:
			if (AlgorithmMembers.TypeOfExactColumnGenerator == 1){
				//Initialize LayeredLabels:

				//Obtain LayeredLabels
				CurrentPricingProblem.LayeredLabels = new std::list<Label>*[MaxNumberOfArcs];
				std::list<Label>** LayeredLabels = CurrentPricingProblem.LayeredLabels;
				//Initialize source label
				LayeredLabels[0] = new std::list<Label>[1];
				//Initalize labels for layers 1 to MaxNumberOfArcs
				for (short int Layer = 1; Layer < MaxNumberOfArcs; ++Layer){
					LayeredLabels[Layer] = new std::list<Label>[NumberOfCustomers];
				}
			}
			else if (AlgorithmMembers.TypeOfExactColumnGenerator == 2){
				//Initialize DoubleLayeredLabels:
				
				//Obtain integer capacity
				int IntegerCapacity = int(CurrentCapacity);
				//Initialize DoubleLayeredLabels
				CurrentPricingProblem.DoubleLayeredLabels = new std::list<Label>**[IntegerCapacity+1];
				std::list<Label>*** DoubleLayeredLabels = CurrentPricingProblem.DoubleLayeredLabels;
				//Initialize source label
				DoubleLayeredLabels[0] = new std::list<Label>*[1];
				DoubleLayeredLabels[0][0] = new std::list<Label>[1];
				//Initialize layers to 0 for all loads
				for (int q = 1; q <= IntegerCapacity; ++q){
					DoubleLayeredLabels[q] = 0;
				}
			}
			else if (AlgorithmMembers.TypeOfExactColumnGenerator == 3){
				
				//Compute the maximum number of forward and backward labels
				short int MaxNumberOfBackwardArcs = (MaxNumberOfArcs - 1) / 2;
				short int MaxNumberOfForwardArcs = MaxNumberOfArcs - 1 - MaxNumberOfBackwardArcs;
				//Set maximum number of forward and backward labels.
				CurrentPricingProblem.MaxNumberOfForwardArcs = MaxNumberOfForwardArcs;
				CurrentPricingProblem.MaxNumberOfBackwardArcs = MaxNumberOfBackwardArcs;
				//Set the number of forward and backward layers
				CurrentPricingProblem.NumberOfForwardLayers = MaxNumberOfForwardArcs + 1;
				CurrentPricingProblem.NumberOfBackwardLayers = MaxNumberOfBackwardArcs + 1;
				
				//Initialize bi-directional labels:
				
				//Obtain LayeredForwardLabels
				CurrentPricingProblem.LayeredForwardLabels = new std::list<ForwardLabel>*[MaxNumberOfForwardArcs+1];
				std::list<ForwardLabel>** LayeredForwardLabels = CurrentPricingProblem.LayeredForwardLabels;
				//Initialize source label
				LayeredForwardLabels[0] = new std::list<ForwardLabel>[1];
				//Initalize labels for layers 1 to MaxNumberOfForwardArcs
				for (short int Layer = 1; Layer <= MaxNumberOfForwardArcs; ++Layer){
					LayeredForwardLabels[Layer] = new std::list<ForwardLabel>[NumberOfCustomers];
				}

				//Obtain LayeredBackwardLabels
				CurrentPricingProblem.LayeredBackwardLabels = new std::list<BackwardLabel>*[MaxNumberOfBackwardArcs+1];
				std::list<BackwardLabel>** LayeredBackwardLabels = CurrentPricingProblem.LayeredBackwardLabels;
				//Initialize source label
				LayeredBackwardLabels[0] = new std::list<BackwardLabel>[1];
				//Initalize labels for layers 1 to MaxNumberOfBackwardArcs
				for (short int Layer = 1; Layer <= MaxNumberOfBackwardArcs; ++Layer){
					LayeredBackwardLabels[Layer] = new std::list<BackwardLabel>[NumberOfCustomers];
				}
				
			}
			else if (AlgorithmMembers.TypeOfExactColumnGenerator == 4){

				//Compute the maximum number of forward and backward labels
				short int MaxNumberOfBackwardArcs = (MaxNumberOfArcs - 1) / 2;
				short int MaxNumberOfForwardArcs = MaxNumberOfArcs - 1 - MaxNumberOfBackwardArcs;
				//Set maximum number of forward and backward labels.
				CurrentPricingProblem.MaxNumberOfForwardArcs = MaxNumberOfForwardArcs;
				CurrentPricingProblem.MaxNumberOfBackwardArcs = MaxNumberOfBackwardArcs;
				//Set the number of forward and backward layers
				CurrentPricingProblem.NumberOfForwardLayers = MaxNumberOfArcs;
				CurrentPricingProblem.NumberOfBackwardLayers = MaxNumberOfArcs;

				//Initialize bi-directional labels and completion bounds:

				//Obtain LayeredForwardLabels
				CurrentPricingProblem.LayeredForwardLabels = new std::list<ForwardLabel>*[MaxNumberOfArcs];
				std::list<ForwardLabel>** LayeredForwardLabels = CurrentPricingProblem.LayeredForwardLabels;
				//Initialize source label
				LayeredForwardLabels[0] = new std::list<ForwardLabel>[1];
				//Obtain ForwardCompletionBounds
				CurrentPricingProblem.ForwardCompletionBounds = new double*[MaxNumberOfArcs];
				double** ForwardCompletionBounds = CurrentPricingProblem.ForwardCompletionBounds;
				//Initalize labels and completion bounds for layers 1 to MaxNumberOfArcs
				for (short int Layer = 1; Layer < MaxNumberOfArcs; ++Layer){
					LayeredForwardLabels[Layer] = new std::list<ForwardLabel>[NumberOfCustomers];
					ForwardCompletionBounds[Layer] = new double[NumberOfCustomers];
				}

				//Obtain LayeredBackwardLabels
				CurrentPricingProblem.LayeredBackwardLabels = new std::list<BackwardLabel>*[MaxNumberOfArcs];
				std::list<BackwardLabel>** LayeredBackwardLabels = CurrentPricingProblem.LayeredBackwardLabels;
				//Initialize source label
				LayeredBackwardLabels[0] = new std::list<BackwardLabel>[1];
				//Obtain BackwardCompletionBounds
				CurrentPricingProblem.BackwardCompletionBounds = new double*[MaxNumberOfArcs];
				double** BackwardCompletionBounds = CurrentPricingProblem.BackwardCompletionBounds;
				//Initalize labels and bounds for layers 1 to MaxNumberOfArcs
				for (short int Layer = 1; Layer < MaxNumberOfArcs; ++Layer){
					LayeredBackwardLabels[Layer] = new std::list<BackwardLabel>[NumberOfCustomers];
					BackwardCompletionBounds[Layer] = new double[NumberOfCustomers];
				}

			}


		}//end check whether the current pricing problem is relevant.

	}//end go over all customers

	//Create customer-depot pricing problem
	//****************************************************************************
	
	//Go over all customers
	for (short int i = 1; i < NumberOfLocations; ++i){

		//Initialize support
		short int iMinusOne = i - 1;

		//Determine the capacity:
		//Initialize excluded customers
		std::vector<short int> ExcludedCustomers;
		ExcludedCustomers.push_back(iMinusOne);

		//Compute the remaining capacity after subtracting the smallest p-1 customers, not i.
		double RemainingCapacity = ComputeCapacity(VehicleCapacity, SmallestDemand, SmallestDemandLocations, p - 1, ExcludedCustomers, NumberOfCustomers);

		//Obtain demand perturbed by tolerance:
		double PerturbedDemand = Demand[iMinusOne] - Tolerance;

		//Check whether the pricing problem is relevant
		if (PerturbedDemand <= RemainingCapacity && p < NumberOfLocations){
			//The pricing problem to create a p-step from the depot to customer j is relevant.

			//Add the pricing problem to the list of relevant pricing problems
			AlgorithmMembers.RelevantCustomerDepotPricingProblems.push_back(&(PricingProblems[i][0]));

			//Obtain the current pricing problem
			PStep_PricingProblem& CurrentPricingProblem = PricingProblems[i][0];

			//Create customer-depot pricing problem:
			//--------------------------------------

			//Set number of locations and customers
			CurrentPricingProblem.NumberOfLocations = NumberOfLocations;
			CurrentPricingProblem.NumberOfCustomers = NumberOfCustomers;
			//Set source and sink
			CurrentPricingProblem.Source = i;
			CurrentPricingProblem.Sink = 0;
			//Copy demand
			CurrentPricingProblem.Demand = new double[NumberOfCustomers];
			for (short int j = 0; j < NumberOfCustomers; ++j){
				CurrentPricingProblem.Demand[j] = Demand[j];
			}
			//Set indicator indicating that the Sink may not be reached before p arcs are traversed.
			CurrentPricingProblem.AllowedToReachSinkEarly = false;
			//Initialize reduced costs
			CurrentPricingProblem.ReducedArcCosts = new double*[NumberOfLocations];
			for (short int j = 0; j < NumberOfLocations; ++j){
				CurrentPricingProblem.ReducedArcCosts[j] = new double[NumberOfLocations];
			}

			//Set the maximum number of arcs that can be traversed on the customer-depot p-step:
			//Initialize the maximum number of arcs
			CurrentPricingProblem.MaxNumberOfArcs=p;

			//Set capacity
			CurrentPricingProblem.Capacity = VehicleCapacity;

			//Initialize and set reachability
			CurrentPricingProblem.Reachable = new bool[NumberOfCustomers];
			bool* Reachable = CurrentPricingProblem.Reachable;
			//Check whether p=1
			if (p == 1){
				//Go over all customers and set reachability
				for (short int j = 0; j < NumberOfCustomers; ++j){
					Reachable[j] = false;
				}
			}
			else{
				//p is larger than 1.
				//Go over all customers and set reachability
				for (short int j = 0; j < NumberOfCustomers; ++j){
					//Check whether i=j
					if (j == iMinusOne){
						//The source is not reachable
						Reachable[iMinusOne] = false;
					}
					else{
						//Check whether an i-0 p-step exists in which also j is visited.

						//Initialize current excluded customers
						std::vector<short int> CurrentExcludedCustomers = ExcludedCustomers;
						CurrentExcludedCustomers.push_back(j);
						//Compute the remaining capacity after subtracting the smallest p-2 customers, not i and j.
						double RemainingCapacity = ComputeCapacity(VehicleCapacity, SmallestDemand, SmallestDemandLocations, p - 2, ExcludedCustomers, NumberOfCustomers);

						//Check whether the demand of j and k is less than the vehicle capacity
						if (PerturbedDemand + Demand[j] <= RemainingCapacity){
							//Customer j could be visited
							Reachable[j] = true;
						}
						else{
							//Customer j cannot be reached
							Reachable[j] = false;
						}
					}
				}
			}

			//Initialize ArcIsRelevant
			CurrentPricingProblem.ArcIsRelevant = new bool*[NumberOfLocations];
			bool** ArcIsRelevant = CurrentPricingProblem.ArcIsRelevant;
			for (short int k = 0; k < NumberOfLocations; ++k){
				ArcIsRelevant[k] = new bool[NumberOfLocations];
				for (short int l = 0; l < NumberOfLocations; ++l){
					ArcIsRelevant[k][l] = false;
				}
			}

			//OLD CODE:
			//This does not work nicely for ng path relaxation...
			/*
			//Initialize and set LayeredCapacity
			CurrentPricingProblem.LayeredCapacity = new double[p + 1];
			CurrentPricingProblem.LayeredCapacityPlusTolerance = new double[p + 1];
			double* (&LayeredCapacity) = CurrentPricingProblem.LayeredCapacity;
			double* (&LayeredCapacityPlusTolerance) = CurrentPricingProblem.LayeredCapacityPlusTolerance;
			//Set bogus capacity for Layer 0
			LayeredCapacity[0] = 0;
			LayeredCapacityPlusTolerance[0] = 0;
			//Set the final LayeredCapacities
			LayeredCapacity[p - 1] = VehicleCapacity;
			LayeredCapacityPlusTolerance[p - 1] = VehicleCapacityPlusTolerance;
			LayeredCapacity[p] = VehicleCapacity;
			LayeredCapacityPlusTolerance[p] = VehicleCapacityPlusTolerance;
			//Initialize smallest index
			short int SmallestIndex = 0;
			//Go over all remaining layers
			for (short int Layer = p - 2; 0 < Layer; --Layer){
				//Compute and set the LayeredCapacity
				LayeredCapacity[Layer] = ComputeCapacity(LayeredCapacity[Layer+1], SmallestDemand, SmallestDemandLocations, 1, ExcludedCustomers, NumberOfCustomers, SmallestIndex);
				LayeredCapacityPlusTolerance[Layer] = LayeredCapacity[Layer] + Tolerance;
			}
			*/

			//Set OutArcs and ArcIsRelevant:
			//------------------------------

			//Obtain out arcs
			std::vector<std::vector<PStep_Arc>>& OutArcs = AlgorithmMembers.OutArcs;
			std::list<std::list<PStep_Arc>>& PPOutArcs = CurrentPricingProblem.OutArcs;

			//Set OutArcs and ArcIsRelevant for the arcs out of the depot
			auto CurrentOutArcs = OutArcs.begin();
			//Initialize arcs out of the depot and leave empty
			PPOutArcs.emplace_front();

			//Go over all locations as start nodes of the arc:
			//Initialize from arcs out of location 1.
			short int StartNode = 1;
			++CurrentOutArcs;
			for (; CurrentOutArcs != OutArcs.end(); ++CurrentOutArcs){

				//Initialize new list
				PPOutArcs.emplace_back();

				//Check whether the current start node is the source
				if (StartNode==i){

					//Obtain new list
					std::list<PStep_Arc>& CurrentPPOutArcs = PPOutArcs.back();

					//Check if p=1
					if (p == 1){
						//Only the arc (i,0) is feasible
						
						//Go over all OutArcs and check whether (i,0) is present.
						for (auto CurrentArc = CurrentOutArcs->begin(); CurrentArc != CurrentOutArcs->end(); ++CurrentArc){
							//Obtain the end node
							short int EndNode = CurrentArc->EndNode;
							//Check whether the Endnode is the depot
							if (EndNode == 0){
								//Add the out arc
								CurrentPPOutArcs.push_back(*CurrentArc);
								ArcIsRelevant[i][EndNode] = true;
							}
						}
					}
					else{
						//The start node is the source, all arcs out of i into a reachable node, except the depot, are feasible.

						//Go over all OutArcs
						for (auto CurrentArc = CurrentOutArcs->begin(); CurrentArc != CurrentOutArcs->end(); ++CurrentArc){
							//Obtain the end node
							short int EndNode = CurrentArc->EndNode;
							//Check whether the Endnode is not the depot and reachable
							if (EndNode != 0 && Reachable[EndNode - 1]){
								//Add the out arc
								CurrentPPOutArcs.push_back(*CurrentArc);
								ArcIsRelevant[i][EndNode] = true;
							}
						}
					}
				}
				//Check whether the current start node is reachable
				else if (Reachable[StartNode - 1] && p > 1){
					//StartNode is reachable and it is not the source, set OutArcs.

					//Obtain new list
					std::list<PStep_Arc>& CurrentPPOutArcs = PPOutArcs.back();

					//Check whether p==2
					if (p == 2){
						//If p==2, only arcs (StartNode,0) are feasible

						//Go over all arcs out of StartNode
						for (auto CurrentArc = CurrentOutArcs->begin(); CurrentArc != CurrentOutArcs->end(); ++CurrentArc){
							//Obtain end node
							short int EndNode = CurrentArc->EndNode;
							//Check whether the end node is the sink, in which case the arc can immediately be added
							if (EndNode == 0){
								//The arc (StartNode,0) can be used, add it and make it relevant
								CurrentPPOutArcs.push_back(*CurrentArc);
								ArcIsRelevant[StartNode][EndNode] = true;
							}
						}

					}
					else{
						//p>2

						//Initialize current excluded customers
						std::vector<short int> CurrentExcludedCustomers = ExcludedCustomers;
						CurrentExcludedCustomers.push_back(StartNode - 1);
						//Initialize the perturbed demand of j and the startNode
						double PerturbedDemandIncludingStartNode = PerturbedDemand + Demand[StartNode - 1];

						//Go over all arcs out of StartNode
						for (auto CurrentArc = CurrentOutArcs->begin(); CurrentArc != CurrentOutArcs->end(); ++CurrentArc){

							//Obtain end node
							short int EndNode = CurrentArc->EndNode;
							//Check whether the end node is the sink, in which case the arc can immediately be added
							if (EndNode == 0){
								//The arc (StartNode,0) can be used, add it and make it relevant
								CurrentPPOutArcs.push_back(*CurrentArc);
								ArcIsRelevant[StartNode][EndNode] = true;
							}
							//Otherwise check whether the end node is not the source
							else if (EndNode != i){
								//StartNode, EndNode and i are all different, check whether they fit on one p-step.

								//Initialize current excluded customers
								std::vector<short int> CurrentExcludedCustomersForArc = CurrentExcludedCustomers;
								CurrentExcludedCustomersForArc.push_back(EndNode - 1);
								//Compute the remaining capacity after subtracting the p-3 smallest customers not i, StartNode and EndNode
								RemainingCapacity = ComputeCapacity(VehicleCapacity, SmallestDemand, SmallestDemandLocations, p - 3, CurrentExcludedCustomersForArc, NumberOfCustomers);

								//Check capacity
								if (PerturbedDemandIncludingStartNode + Demand[EndNode - 1] <= RemainingCapacity){
									//The arc (StartNode,EndNode) can be used, add it and make it relevant
									CurrentPPOutArcs.push_back(*CurrentArc);
									ArcIsRelevant[StartNode][EndNode] = true;
								}

							}

						}//end go over all arcs out of StartNode

					}//end else p>2

				}//end check whether the current start node is reachable

				//Update StartNode
				++StartNode;

			}//end go over all locations as start node of the arc

			//Set InArcs:
			//------------------------------

			//Obtain in arcs
			std::vector<std::vector<PStep_Arc>>& InArcs = AlgorithmMembers.InArcs;
			std::list<std::list<PStep_Arc>>& PPInArcs = CurrentPricingProblem.InArcs;

			//Set InArcs and ArcIsRelevant for the arcs into the depot
			auto CurrentInArcs = InArcs.begin();
			//Initialize arcs into the depot and leave it empty
			PPInArcs.emplace_front();
			//Obtain arcs into the depot
			std::list<PStep_Arc>& PPArcsIntoDepot = PPInArcs.front();
			//Check if p==1
			if (p == 1){
				//Only the arc (i,0) is relevant.
				//Go over all InArcs and check whether (i,0) is present.
				for (auto CurrentArc = CurrentInArcs->begin(); CurrentArc != CurrentInArcs->end(); ++CurrentArc){
					//Obtain the end node
					short int StartNode = CurrentArc->StartNode;
					//Check whether the StartNode is the depot
					if (StartNode == i){
						//Add the in arc
						PPArcsIntoDepot.push_back(*CurrentArc);
					}
				}

				//Add empty InArcs
				++CurrentInArcs;
				for (; CurrentInArcs != InArcs.end(); ++CurrentInArcs){
					//Initialize new list
					PPInArcs.emplace_back();
				}

			}
			else{
				//Go over all arcs into the depot
				for (auto CurrentArc = CurrentInArcs->begin(); CurrentArc != CurrentInArcs->end(); ++CurrentArc){
					//Obtain start node
					short int StartNode = CurrentArc->StartNode;
					//Check whether the start node can be reached
					if (Reachable[StartNode - 1]){
						//The arc (StartNode,0) can be used, add it and make it relevant
						PPArcsIntoDepot.push_back(*CurrentArc);
					}
				}//end go over all arcs into the depot

				//Go over all locations as end nodes of the arc:
				//Initialize from arcs into location 1.
				short int EndNode = 1;
				++CurrentInArcs;
				for (; CurrentInArcs != InArcs.end(); ++CurrentInArcs){

					//Initialize new list
					PPInArcs.emplace_back();

					//Check whether the current end node is reachable and not the source
					if (EndNode != i && Reachable[EndNode - 1]){
						//EndNode is reachable and not the source, set InArcs.

						//Obtain new list
						std::list<PStep_Arc>& CurrentPPInArcs = PPInArcs.back();

						//Check whether p==2
						if (p == 2){
							//If p==2, only arcs (i,EndNode) are feasible

							//Go over all arcs into EndNode
							for (auto CurrentArc = CurrentInArcs->begin(); CurrentArc != CurrentInArcs->end(); ++CurrentArc){
								//Obtain start node
								StartNode = CurrentArc->StartNode;
								//Check whether the start node is the source, in which case the arc can immediately be added
								if (StartNode == i){
									//The arc (i,EndNode) can be used, add it.
									CurrentPPInArcs.push_back(*CurrentArc);
								}
							}

						}
						else{
							//p>2

							//Initialize current excluded customers
							std::vector<short int> CurrentExcludedCustomers = ExcludedCustomers;
							CurrentExcludedCustomers.push_back(EndNode - 1);
							//Initialize the perturbed demand of j and the end node
							double PerturbedDemandIncludingEndNode = PerturbedDemand + Demand[EndNode - 1];

							//Go over all arcs into EndNode
							for (auto CurrentArc = CurrentInArcs->begin(); CurrentArc != CurrentInArcs->end(); ++CurrentArc){

								//Obtain start node
								StartNode = CurrentArc->StartNode;
								//Check whether the start node is the source, in which case the arc can immediately be added
								if (StartNode == i){
									//The arc (i,EndNode) can be used, add it.
									CurrentPPInArcs.push_back(*CurrentArc);
								}
								//Otherwise check whether the start node is not the depot
								else if (StartNode != 0){
									//StartNode, EndNode and i are all different, check whether they fit on one p-step.

									//Initialize current excluded customers
									std::vector<short int> CurrentExcludedCustomersForArc = CurrentExcludedCustomers;
									CurrentExcludedCustomersForArc.push_back(StartNode - 1);
									//Compute the remaining capacity after subtracting the p-3 smallest customers not i, StartNode and EndNode
									RemainingCapacity = ComputeCapacity(VehicleCapacity, SmallestDemand, SmallestDemandLocations, p - 3, CurrentExcludedCustomersForArc, NumberOfCustomers);

									//Check capacity
									if (PerturbedDemandIncludingEndNode + Demand[StartNode - 1] <= RemainingCapacity){
										//The arc (StartNode,EndNode) can be used, add it.
										CurrentPPInArcs.push_back(*CurrentArc);
									}

								}

							}//end go over all arcs into EndNode

						}//end else p>2

					}//end check whether the current end node is reachable

					//Update EndNode
					++EndNode;

				}//end go over all locations as end node of the arc

			}//end else p!=1

			//Initialize Labels:
			if (AlgorithmMembers.TypeOfExactColumnGenerator == 1){
				//Initialize LayeredLabels:

				//Obtain LayeredLabels
				CurrentPricingProblem.LayeredLabels = new std::list<Label>*[p];
				std::list<Label>** LayeredLabels = CurrentPricingProblem.LayeredLabels;
				//Initialize source label
				LayeredLabels[0] = new std::list<Label>[1];
				//Initalize labels for layers 1 to p
				for (short int Layer = 1; Layer < p; ++Layer){
					LayeredLabels[Layer] = new std::list<Label>[NumberOfCustomers];
				}
			}
			else if (AlgorithmMembers.TypeOfExactColumnGenerator == 2){
				//Initialize DoubleLayeredLabels:

				//Obtain integer capacity
				int IntegerCapacity = int(VehicleCapacity);
				//Initialize DoubleLayeredLabels
				CurrentPricingProblem.DoubleLayeredLabels = new std::list<Label>**[IntegerCapacity+1];
				std::list<Label>*** DoubleLayeredLabels = CurrentPricingProblem.DoubleLayeredLabels;
				//Initialize source label
				DoubleLayeredLabels[0] = new std::list<Label>*[1];
				DoubleLayeredLabels[0][0] = new std::list<Label>[1];
				//Initialize layers to 0 for all loads
				for (int q = 1; q <= IntegerCapacity; ++q){
					DoubleLayeredLabels[q] = 0;
				}
			}
			else if (AlgorithmMembers.TypeOfExactColumnGenerator == 3){
				
				//Compute the maximum number of forward and backward labels
				short int MaxNumberOfBackwardArcs = (p - 1) / 2;
				short int MaxNumberOfForwardArcs = p - 1 - MaxNumberOfBackwardArcs;
				//Set maximum number of forward and backward labels.
				CurrentPricingProblem.MaxNumberOfForwardArcs = MaxNumberOfForwardArcs;
				CurrentPricingProblem.MaxNumberOfBackwardArcs = MaxNumberOfBackwardArcs;
				//Set the number of forward and backward layers
				CurrentPricingProblem.NumberOfForwardLayers = MaxNumberOfForwardArcs + 1;
				CurrentPricingProblem.NumberOfBackwardLayers = MaxNumberOfBackwardArcs + 1;
				
				//Initialize bi-directional labels:
				
				//Obtain LayeredForwardLabels
				CurrentPricingProblem.LayeredForwardLabels = new std::list<ForwardLabel>*[MaxNumberOfForwardArcs+1];
				std::list<ForwardLabel>** LayeredForwardLabels = CurrentPricingProblem.LayeredForwardLabels;
				//Initialize source label
				LayeredForwardLabels[0] = new std::list<ForwardLabel>[1];
				//Initalize labels for layers 1 to MaxNumberOfForwardArcs
				for (short int Layer = 1; Layer <= MaxNumberOfForwardArcs; ++Layer){
					LayeredForwardLabels[Layer] = new std::list<ForwardLabel>[NumberOfCustomers];
				}

				//Obtain LayeredBackwardLabels
				CurrentPricingProblem.LayeredBackwardLabels = new std::list<BackwardLabel>*[MaxNumberOfBackwardArcs+1];
				std::list<BackwardLabel>** LayeredBackwardLabels = CurrentPricingProblem.LayeredBackwardLabels;
				//Initialize source label
				LayeredBackwardLabels[0] = new std::list<BackwardLabel>[1];
				//Initalize labels for layers 1 to MaxNumberOfBackwardArcs
				for (short int Layer = 1; Layer <= MaxNumberOfBackwardArcs; ++Layer){
					LayeredBackwardLabels[Layer] = new std::list<BackwardLabel>[NumberOfCustomers];
				}
				
			}
			else if (AlgorithmMembers.TypeOfExactColumnGenerator == 4){

				//Compute the maximum number of forward and backward labels
				short int MaxNumberOfBackwardArcs = (p - 1) / 2;
				short int MaxNumberOfForwardArcs = p - 1 - MaxNumberOfBackwardArcs;
				//Set maximum number of forward and backward labels.
				CurrentPricingProblem.MaxNumberOfForwardArcs = MaxNumberOfForwardArcs;
				CurrentPricingProblem.MaxNumberOfBackwardArcs = MaxNumberOfBackwardArcs;
				//Set the number of forward and backward layers
				CurrentPricingProblem.NumberOfForwardLayers = p;
				CurrentPricingProblem.NumberOfBackwardLayers = p;

				//Initialize bi-directional labels and completion bounds:

				//Obtain LayeredForwardLabels
				CurrentPricingProblem.LayeredForwardLabels = new std::list<ForwardLabel>*[p];
				std::list<ForwardLabel>** LayeredForwardLabels = CurrentPricingProblem.LayeredForwardLabels;
				//Initialize source label
				LayeredForwardLabels[0] = new std::list<ForwardLabel>[1];
				//Obtain ForwardCompletionBounds
				CurrentPricingProblem.ForwardCompletionBounds = new double*[p];
				double** ForwardCompletionBounds = CurrentPricingProblem.ForwardCompletionBounds;
				//Initalize labels and bounds for layers 1 to p
				for (short int Layer = 1; Layer < p; ++Layer){
					LayeredForwardLabels[Layer] = new std::list<ForwardLabel>[NumberOfCustomers];
					ForwardCompletionBounds[Layer] = new double[NumberOfCustomers];
				}

				//Obtain LayeredBackwardLabels
				CurrentPricingProblem.LayeredBackwardLabels = new std::list<BackwardLabel>*[p];
				std::list<BackwardLabel>** LayeredBackwardLabels = CurrentPricingProblem.LayeredBackwardLabels;
				//Initialize source label
				LayeredBackwardLabels[0] = new std::list<BackwardLabel>[1];
				//Obtain BackwardCompletionBounds
				CurrentPricingProblem.BackwardCompletionBounds = new double*[p];
				double** BackwardCompletionBounds = CurrentPricingProblem.BackwardCompletionBounds;
				//Initalize labels and bounds for layers 1 to p
				for (short int Layer = 1; Layer < p; ++Layer){
					LayeredBackwardLabels[Layer] = new std::list<BackwardLabel>[NumberOfCustomers];
					BackwardCompletionBounds[Layer] = new double[NumberOfCustomers];
				}

			}

		}//end check whether the current pricing problem is relevant.

	}//end go over all customers

	//Create customer-customer pricing problem
	//****************************************************************************

	//Go over all pairs of customers
	for (short int i = 1; i < NumberOfCustomers; ++i){
		for (short int j = i + 1; j < NumberOfLocations; ++j){

			//Initialize support
			short int iMinusOne = i - 1;
			short int jMinusOne = j - 1;

			//Determine the capacity:
			//Initialize excluded customers
			std::vector<short int> ExcludedCustomers;
			ExcludedCustomers.push_back(iMinusOne);
			ExcludedCustomers.push_back(jMinusOne);
			//Initialize Index indicating which customer is the first next smallest customer to consider
			short int SmallestIndex = 0;
			//Initialize current capacity, it will represent adjusted or vehicle capacity depending on UseAdjustedCapacity.
			double CurrentCapacity;
			//Check whether adjusted capacity is used
			if (UseAdjustedCapacity){
				//Compute adjusted capacity
				CurrentCapacity = ComputeCapacity(VehicleCapacity, SmallestDemand, SmallestDemandLocations, p - 1, ExcludedCustomers, NumberOfCustomers, SmallestIndex);
			}
			else{
				//Set the current capacity to vehicle capacity
				CurrentCapacity = VehicleCapacity;
			}

			//OLD CODE:
			//This does not work nicely for ng path relaxation...
			/*
			//Store the value of SmallestIndex for use when setting LayeredCapacity
			short int LayeredCapacitySmallestIndex = SmallestIndex;
			*/

			//Compute the remaining capacity after visiting p-1 customers on the i-j p-step.
			double RemainingCapacity = ComputeCapacity(CurrentCapacity, SmallestDemand, SmallestDemandLocations, p - 1, ExcludedCustomers, NumberOfCustomers, SmallestIndex);

			//Obtain demand perturbed by tolerance:
			double PerturbedDemand = Demand[iMinusOne] + Demand[jMinusOne] - Tolerance;

			//Check whether the pricing problem is relevant
			if (PerturbedDemand <= RemainingCapacity && 2*p < NumberOfLocations){
				//The pricing problem to create a p-step from the depot to customer j is relevant.

				//Add the pricing problem to the list of relevant pricing problems
				//NOTE: j-i is also relevant, but not added to the list.
				AlgorithmMembers.RelevantCustomerCustomerPricingProblems.push_back(&(PricingProblems[i][j]));

				//Obtain the current pricing problems
				PStep_PricingProblem& ijPricingProblem = PricingProblems[i][j];
				PStep_PricingProblem& jiPricingProblem = PricingProblems[j][i];

				//Create depot-customer pricing problem:
				//--------------------------------------

				//Set number of locations and customers
				ijPricingProblem.NumberOfLocations = NumberOfLocations;
				ijPricingProblem.NumberOfCustomers = NumberOfCustomers;
				jiPricingProblem.NumberOfLocations = NumberOfLocations;
				jiPricingProblem.NumberOfCustomers = NumberOfCustomers;
				//Set source and sink
				ijPricingProblem.Source = i;
				ijPricingProblem.Sink = j;
				jiPricingProblem.Source = j;
				jiPricingProblem.Sink = i;
				//Copy demand
				ijPricingProblem.Demand = new double[NumberOfCustomers];
				jiPricingProblem.Demand = new double[NumberOfCustomers];
				for (short int k = 0; k < NumberOfCustomers; ++k){
					ijPricingProblem.Demand[k] = Demand[k];
					jiPricingProblem.Demand[k] = Demand[k];
				}
				//Set indicator indicating that the Sink may be reached before p arcs are traversed.
				ijPricingProblem.AllowedToReachSinkEarly = false;
				jiPricingProblem.AllowedToReachSinkEarly = false;
				//Initialize reduced costs
				ijPricingProblem.ReducedArcCosts = new double*[NumberOfLocations];
				jiPricingProblem.ReducedArcCosts = new double*[NumberOfLocations];
				for (short int k = 0; k < NumberOfLocations; ++k){
					ijPricingProblem.ReducedArcCosts[k] = new double[NumberOfLocations];
					jiPricingProblem.ReducedArcCosts[k] = new double[NumberOfLocations];
				}

				//Set the maximum number of arcs
				ijPricingProblem.MaxNumberOfArcs = p;
				jiPricingProblem.MaxNumberOfArcs = p;

				//Set capacity
				ijPricingProblem.Capacity = CurrentCapacity;
				jiPricingProblem.Capacity = CurrentCapacity;

				//Initialize and set reachability
				ijPricingProblem.Reachable = new bool[NumberOfCustomers];
				jiPricingProblem.Reachable = new bool[NumberOfCustomers];
				bool* ijReachable = ijPricingProblem.Reachable;
				bool* jiReachable = jiPricingProblem.Reachable;
				//Check whether MaxNumberOfArcs=1
				if (p == 1){
					//Go over all customers and set reachability
					for (short int k = 0; k < NumberOfCustomers; ++k){
						ijReachable[k] = false;
						jiReachable[k] = false;
					}
					ijReachable[jMinusOne] = true;
					jiReachable[iMinusOne] = true;
				}
				else{
					//Go over all customers and set reachability
					for (short int k = 0; k < NumberOfCustomers; ++k){
						//Check whether k=j
						if (k == jMinusOne){
							//The sink is reachable and the source is not
							ijReachable[jMinusOne] = true;
							jiReachable[jMinusOne] = false;
						}
						else if (k == iMinusOne){
							//The sink is reachable and the source is not
							ijReachable[iMinusOne] = false;
							jiReachable[iMinusOne] = true;
						}
						else{
							//Check whether an i-j p-step exists in which also k is visited.

							//Initialize current excluded customers
							std::vector<short int> CurrentExcludedCustomers = ExcludedCustomers;
							CurrentExcludedCustomers.push_back(k);
							//Compute capacity taking into account that i is visited
							//Check whether adjusted capacity is used
							if (UseAdjustedCapacity){
								//Compute adjusted capacity
								RemainingCapacity = ComputeCapacity(VehicleCapacity, SmallestDemand, SmallestDemandLocations, 2*p - 3, CurrentExcludedCustomers, NumberOfCustomers);
							}
							else{
								//Compute adjusted capacity
								RemainingCapacity = ComputeCapacity(VehicleCapacity, SmallestDemand, SmallestDemandLocations, p-2, CurrentExcludedCustomers, NumberOfCustomers);
							}
							//Check whether the demand is less than the vehicle capacity
							if (PerturbedDemand + Demand[k] <= RemainingCapacity){
								//Customer k could be visited
								ijReachable[k] = true;
								jiReachable[k] = true;
							}
							else{
								//Customer k cannot be reached
								ijReachable[k] = false;
								jiReachable[k] = false;
							}
						}
					}
				}

				//Initialize ArcIsRelevant
				ijPricingProblem.ArcIsRelevant = new bool*[NumberOfLocations];
				jiPricingProblem.ArcIsRelevant = new bool*[NumberOfLocations];
				bool** ijArcIsRelevant = ijPricingProblem.ArcIsRelevant;
				bool** jiArcIsRelevant = jiPricingProblem.ArcIsRelevant;
				for (short int k = 0; k < NumberOfLocations; ++k){
					ijArcIsRelevant[k] = new bool[NumberOfLocations];
					jiArcIsRelevant[k] = new bool[NumberOfLocations];
					for (short int l = 0; l < NumberOfLocations; ++l){
						ijArcIsRelevant[k][l] = false;
						jiArcIsRelevant[k][l] = false;
					}
				}

				//OLD CODE:
				//This does not work nicely for ng path relaxation...
				/*
				//Initialize and set LayeredCapacity
				ijPricingProblem.LayeredCapacity = new double[p + 1];
				jiPricingProblem.LayeredCapacity = new double[p + 1];
				ijPricingProblem.LayeredCapacityPlusTolerance = new double[p + 1];
				jiPricingProblem.LayeredCapacityPlusTolerance = new double[p + 1];
				double* (&ijLayeredCapacity) = ijPricingProblem.LayeredCapacity;
				double* (&jiLayeredCapacity) = jiPricingProblem.LayeredCapacity;
				double* (&ijLayeredCapacityPlusTolerance) = ijPricingProblem.LayeredCapacityPlusTolerance;
				double* (&jiLayeredCapacityPlusTolerance) = jiPricingProblem.LayeredCapacityPlusTolerance;
				//Set bogus capacity for Layer 0
				ijLayeredCapacity[0] = 0;
				jiLayeredCapacity[0] = 0;
				ijLayeredCapacityPlusTolerance[0] = 0;
				jiLayeredCapacityPlusTolerance[0] = 0;
				//Set the final LayeredCapacity
				ijLayeredCapacity[p - 1] = CurrentCapacity;
				jiLayeredCapacity[p - 1] = CurrentCapacity;
				ijLayeredCapacityPlusTolerance[p - 1] = CurrentCapacity+Tolerance;
				jiLayeredCapacityPlusTolerance[p - 1] = ijLayeredCapacityPlusTolerance[p - 1];
				ijLayeredCapacity[p] = CurrentCapacity;
				jiLayeredCapacity[p] = CurrentCapacity;
				ijLayeredCapacityPlusTolerance[p] = ijLayeredCapacityPlusTolerance[p - 1];
				jiLayeredCapacityPlusTolerance[p] = jiLayeredCapacityPlusTolerance[p - 1];
				//Go over all remaining layers
				for (short int Layer = p - 2; 0 < Layer; --Layer){
					//Compute and set the LayeredCapacity
					ijLayeredCapacity[Layer] = ComputeCapacity(ijLayeredCapacity[Layer + 1], SmallestDemand, SmallestDemandLocations, 1, ExcludedCustomers, NumberOfCustomers, LayeredCapacitySmallestIndex);
					jiLayeredCapacity[Layer] = ijLayeredCapacity[Layer];
					ijLayeredCapacityPlusTolerance[Layer] = ijLayeredCapacity[Layer] + Tolerance;
					jiLayeredCapacityPlusTolerance[Layer] = ijLayeredCapacityPlusTolerance[Layer];
				}
				*/

				//Set OutArcs and ArcIsRelevant:
				//------------------------------

				//Obtain out arcs
				std::vector<std::vector<PStep_Arc>>& OutArcs = AlgorithmMembers.OutArcs;
				std::list<std::list<PStep_Arc>>& ijOutArcs = ijPricingProblem.OutArcs;
				std::list<std::list<PStep_Arc>>& jiOutArcs = jiPricingProblem.OutArcs;

				//Set OutArcs and ArcIsRelevant for the arcs out of the depot
				auto CurrentOutArcs = OutArcs.begin();
				//Initialize arcs out of the depot and leave empty
				ijOutArcs.emplace_front();
				jiOutArcs.emplace_front();

				//Go over all locations as start nodes of the arc:
				//Initialize from arcs out of location 1.
				short int StartNode = 1;
				++CurrentOutArcs;
				for (; CurrentOutArcs != OutArcs.end(); ++CurrentOutArcs){

					//Initialize new list
					ijOutArcs.emplace_back();
					jiOutArcs.emplace_back();

					//Check whether StartNode is i
					if (StartNode == i){
						//When i is the sink no arc out of i is feasible.

						//Obtain new list
						std::list<PStep_Arc>& ijCurrentOutArcs = ijOutArcs.back();

						//Check whether p=1
						if (p == 1){
							//When i is the source, only (i,j) is feasible.
							
							//Go over all arcs out of StartNode
							for (auto CurrentArc = CurrentOutArcs->begin(); CurrentArc != CurrentOutArcs->end(); ++CurrentArc){
								//Obtain end node
								short int EndNode = CurrentArc->EndNode;
								//Check whether the end node is j, in which case the arc can immediately be added for ijPricingProblem
								if (EndNode == j){
									//The arc (i,j) can be used, add it and make it relevant
									ijCurrentOutArcs.push_back(*CurrentArc);
									ijArcIsRelevant[StartNode][EndNode] = true;
								}
							}

						}
						else{
							//p>1

							//When i is the source, all reachable locations except the depot and sink are feasible

							//Go over all arcs out of StartNode
							for (auto CurrentArc = CurrentOutArcs->begin(); CurrentArc != CurrentOutArcs->end(); ++CurrentArc){
								//Obtain end node
								short int EndNode = CurrentArc->EndNode;
								//Check whether the end node is not the depot, the sink and reachable,
								//in which case the arc can immediately be added for ijPricingProblem
								if (EndNode != 0 && EndNode != j && ijReachable[EndNode - 1]){
									//The arc (i,EndNode) can be used, add it and make it relevant
									ijCurrentOutArcs.push_back(*CurrentArc);
									ijArcIsRelevant[StartNode][EndNode] = true;
								}
							}
						}

					}//end check whether StartNode is i
					//Check whether the StartNode is j
					else if (StartNode == j){
						//When j is the sink no arc out of j is feasible.

						//Obtain new list
						std::list<PStep_Arc>& jiCurrentOutArcs = jiOutArcs.back();

						//Check whether p=1
						if (p == 1){
							//When j is the source, only (j,i) is feasible.

							//Go over all arcs out of StartNode
							for (auto CurrentArc = CurrentOutArcs->begin(); CurrentArc != CurrentOutArcs->end(); ++CurrentArc){
								//Obtain end node
								short int EndNode = CurrentArc->EndNode;
								//Check whether the end node is j, in which case the arc can immediately be added for jiPricingProblem
								if (EndNode == i){
									//The arc (j,i) can be used, add it and make it relevant
									jiCurrentOutArcs.push_back(*CurrentArc);
									jiArcIsRelevant[StartNode][EndNode] = true;
								}
							}

						}
						else{
							//p>1

							//When j is the source, all reachable locations except the depot and sink are feasible

							//Go over all arcs out of StartNode
							for (auto CurrentArc = CurrentOutArcs->begin(); CurrentArc != CurrentOutArcs->end(); ++CurrentArc){
								//Obtain end node
								short int EndNode = CurrentArc->EndNode;
								//Check whether the end node is not the depot, the sink and reachable,
								//in which case the arc can immediately be added for jiPricingProblem
								if (EndNode != 0 && EndNode != i && jiReachable[EndNode - 1]){
									//The arc (j,EndNode) can be used, add it and make it relevant
									jiCurrentOutArcs.push_back(*CurrentArc);
									jiArcIsRelevant[StartNode][EndNode] = true;
								}
							}
						}

					}//end check whether StartNode is j
					//Check whether the current start node is reachable,
					//(ijReachable and jiReachable should be equal for StartNode not i or j)
					//and p>1
					else if (ijReachable[StartNode - 1] && p>1){
						//StartNode is reachable, set OutArcs.

						//Obtain new lists
						std::list<PStep_Arc>& ijCurrentOutArcs = ijOutArcs.back();
						std::list<PStep_Arc>& jiCurrentOutArcs = jiOutArcs.back();

						//Check whether p=2
						if (p == 2){
							//Only the arc (StartNode,Sink) could be used.

							//Go over all arcs out of StartNode
							for (auto CurrentArc = CurrentOutArcs->begin(); CurrentArc != CurrentOutArcs->end(); ++CurrentArc){

								//Obtain end node
								short int EndNode = CurrentArc->EndNode;
								//Check whether the end node is the sink j, in which case the arc can immediately be added to ijPricingProblem
								if (EndNode == j){
									//The arc (StartNode,j) can be used, add it and make it relevant
									ijCurrentOutArcs.push_back(*CurrentArc);
									ijArcIsRelevant[StartNode][EndNode] = true;
								}
								//Check whether the end node is the sink i, in which case the arc can immediately be added to jiPricingProblem
								if (EndNode == i){
									//The arc (StartNode,j) can be used, add it and make it relevant
									jiCurrentOutArcs.push_back(*CurrentArc);
									jiArcIsRelevant[StartNode][EndNode] = true;
								}
							}
						}
						else{
							//p>2

							//Initialize current excluded customers
							std::vector<short int> CurrentExcludedCustomers = ExcludedCustomers;
							CurrentExcludedCustomers.push_back(StartNode - 1);
							//Initialize the perturbed demand of j and the startNode
							double PerturbedDemandIncludingStartNode = PerturbedDemand + Demand[StartNode - 1];

							//Go over all arcs out of StartNode
							for (auto CurrentArc = CurrentOutArcs->begin(); CurrentArc != CurrentOutArcs->end(); ++CurrentArc){

								//Obtain end node
								short int EndNode = CurrentArc->EndNode;
								//Check whether the end node is the sink j, in which case the arc can immediately be added to ijPricingProblem
								if (EndNode == j){
									//The arc (StartNode,EndNode) can be used in ijPricingProblem, add it and make it relevant
									ijCurrentOutArcs.push_back(*CurrentArc);
									ijArcIsRelevant[StartNode][EndNode] = true;
									//The arc (StartNode,EndNode) cannot be used in jiPricingProblem.
								}
								//Check whether the end node is the sink i, in which case the arc can immediately be added to jiPricingProblem
								else if (EndNode == i){
									//The arc (StartNode,EndNode) can be used in jiPricingProblem, add it and make it relevant
									jiCurrentOutArcs.push_back(*CurrentArc);
									jiArcIsRelevant[StartNode][EndNode] = true;
									//The arc (StartNode,EndNode) cannot be used in ijPricingProblem.
								}
								//Otherwise check whether the end node is not the depot
								else if (EndNode != 0){
									//StartNode, EndNode, i and j are all different, check whether they fit on one p-step.

									//Initialize current excluded customers
									std::vector<short int> CurrentExcludedCustomersForArc = CurrentExcludedCustomers;
									CurrentExcludedCustomersForArc.push_back(EndNode - 1);
									//Check whether adjusted capacity is used
									if (UseAdjustedCapacity){
										//Compute remaining capacity
										RemainingCapacity = ComputeCapacity(VehicleCapacity, SmallestDemand, SmallestDemandLocations, 2 * p - 4, CurrentExcludedCustomersForArc, NumberOfCustomers);
									}
									else{
										//Compute remaining capacity
										RemainingCapacity = ComputeCapacity(VehicleCapacity, SmallestDemand, SmallestDemandLocations, p - 3, CurrentExcludedCustomersForArc, NumberOfCustomers);
									}
									//Check capacity
									if (PerturbedDemandIncludingStartNode + Demand[EndNode - 1] <= RemainingCapacity){
										//The arc (StartNode,EndNode) can be used, add it and make it relevant
										ijCurrentOutArcs.push_back(*CurrentArc);
										ijArcIsRelevant[StartNode][EndNode] = true;
										jiCurrentOutArcs.push_back(*CurrentArc);
										jiArcIsRelevant[StartNode][EndNode] = true;
									}

								}

							}//end go over all arcs out of StartNode
						}//end else p>2

					}//end check whether the current start node is reachable

					//Update StartNode
					++StartNode;

				}//end go over all locations as start node of the arc

				//Set InArcs:
				//------------------------------

				//Obtain in arcs
				std::vector<std::vector<PStep_Arc>>& InArcs = AlgorithmMembers.InArcs;
				std::list<std::list<PStep_Arc>>& ijInArcs = ijPricingProblem.InArcs;
				std::list<std::list<PStep_Arc>>& jiInArcs = jiPricingProblem.InArcs;

				//Set InArcs and ArcIsRelevant for the arcs into the depot
				auto CurrentInArcs = InArcs.begin();
				//Initialize arcs into the depot and leave it empty
				ijInArcs.emplace_front();
				jiInArcs.emplace_front();

				//Go over all locations as end nodes of the arc:
				//Initialize from arcs into location 1.
				short int EndNode = 1;
				++CurrentInArcs;
				for (; CurrentInArcs != InArcs.end(); ++CurrentInArcs){

					//Initialize new list
					ijInArcs.emplace_back();
					jiInArcs.emplace_back();

					//Check whether the end node is j
					if (EndNode == j){
						//When j is the source no arc into j is feasible.

						//Obtain new list
						std::list<PStep_Arc>& ijCurrentInArcs = ijInArcs.back();

						//Check whether p=1
						if (p == 1){
							//When j is the sink, only (i,j) is feasible.

							//Go over all InArcs
							for (auto CurrentArc = CurrentInArcs->begin(); CurrentArc != CurrentInArcs->end(); ++CurrentArc){
								//Check whether the startnode is the source
								if (CurrentArc->StartNode == i){
									//Add the in arc
									ijCurrentInArcs.push_back(*CurrentArc);
								}
							}
						}
						else{
							//p>1

							//When j is the sink, all reachable locations except the depot and source are feasible

							//Go over all InArcs
							for (auto CurrentArc = CurrentInArcs->begin(); CurrentArc != CurrentInArcs->end(); ++CurrentArc){
								//Obtain the start node
								StartNode = CurrentArc->StartNode;
								//Check whether the startnode is not the depot, the source and reachable
								if (StartNode != 0 && StartNode !=i && ijReachable[StartNode - 1]){
									//Add the in arc
									ijCurrentInArcs.push_back(*CurrentArc);
								}
							}

						}

					}//end check whether the end node is j
					//Check whether the end node is i
					else if (EndNode == i){
						//When i is the source no arc into i is feasible.

						//Obtain new list
						std::list<PStep_Arc>& jiCurrentInArcs = jiInArcs.back();

						//Check whether p=1
						if (p == 1){
							//When i is the sink, only (j,i) is feasible.

							//Go over all InArcs
							for (auto CurrentArc = CurrentInArcs->begin(); CurrentArc != CurrentInArcs->end(); ++CurrentArc){
								//Check whether the startnode is the source
								if (CurrentArc->StartNode == j){
									//Add the in arc
									jiCurrentInArcs.push_back(*CurrentArc);
								}
							}
						}
						else{
							//p>1

							//When i is the sink, all reachable locations except the depot and source are feasible

							//Go over all InArcs
							for (auto CurrentArc = CurrentInArcs->begin(); CurrentArc != CurrentInArcs->end(); ++CurrentArc){
								//Obtain the start node
								StartNode = CurrentArc->StartNode;
								//Check whether the startnode is not the depot, the source and reachable
								if (StartNode != 0 && StartNode != j && ijReachable[StartNode - 1]){
									//Add the in arc
									jiCurrentInArcs.push_back(*CurrentArc);
								}
							}

						}

					}//end check whether the end node is j
					//Check whether the current start node is reachable,
					//(ijReachable and jiReachable should be equal for StartNode not i or j)
					//and p>1
					else if (ijReachable[EndNode - 1] && p > 1){
						//EndNode is reachable and not the sink, set InArcs.

						//Obtain new list
						std::list<PStep_Arc>& ijCurrentInArcs = ijInArcs.back();
						std::list<PStep_Arc>& jiCurrentInArcs = jiInArcs.back();

						//Check whether p=2
						if (p == 2){
							//Only the arc (Source,EndNode) could be used.

							//Go over all arcs into EndNode
							for (auto CurrentArc = CurrentInArcs->begin(); CurrentArc != CurrentInArcs->end(); ++CurrentArc){
								//Obtain start node
								StartNode = CurrentArc->StartNode;
								//Check whether the start node is the source i, in which case the arc can immediately be added to ijPricingProblem
								if (StartNode == i){
									//The arc (i,EndNode) can be used, add it.
									ijCurrentInArcs.push_back(*CurrentArc);
								}
								//Check whether the start node is the source j, in which case the arc can immediately be added to jiPricingProblem
								if (StartNode == j){
									//The arc (j,EndNode) can be used, add it.
									jiCurrentInArcs.push_back(*CurrentArc);
								}
							}

						}
						else{
							//p>2

							//Initialize current excluded customers
							std::vector<short int> CurrentExcludedCustomers = ExcludedCustomers;
							CurrentExcludedCustomers.push_back(EndNode - 1);
							//Initialize the perturbed demand of j and the end node
							double PerturbedDemandIncludingEndNode = PerturbedDemand + Demand[EndNode - 1];

							//Go over all arcs into EndNode
							for (auto CurrentArc = CurrentInArcs->begin(); CurrentArc != CurrentInArcs->end(); ++CurrentArc){

								//Obtain start node
								StartNode = CurrentArc->StartNode;
								//Check whether the start node is the source i, in which case the arc can immediately be added to ijPricingProblem
								if (StartNode == i){
									//The arc (i,EndNode) can be used, add it.
									ijCurrentInArcs.push_back(*CurrentArc);
								}
								//Check whether the start node is the source j, in which case the arc can immediately be added to jiPricingProblem
								else if (StartNode == j){
									//The arc (j,EndNode) can be used, add it.
									jiCurrentInArcs.push_back(*CurrentArc);
								}
								//Otherwise check whether the start node is not the depot
								else if (StartNode != 0){
									//StartNode, EndNode, i and j are all different, check whether they fit on one p-step.

									//Initialize current excluded customers
									std::vector<short int> CurrentExcludedCustomersForArc = CurrentExcludedCustomers;
									CurrentExcludedCustomersForArc.push_back(StartNode - 1);
									//Check whether adjusted capacity is used
									if (UseAdjustedCapacity){
										//Compute remaining capacity capacity
										RemainingCapacity = ComputeCapacity(VehicleCapacity, SmallestDemand, SmallestDemandLocations, 2*p - 4, CurrentExcludedCustomersForArc, NumberOfCustomers);
									}
									else{
										//Compute remaining capacity capacity
										RemainingCapacity = ComputeCapacity(VehicleCapacity, SmallestDemand, SmallestDemandLocations, p - 3, CurrentExcludedCustomersForArc, NumberOfCustomers);
									}
									//Check capacity
									if (PerturbedDemandIncludingEndNode + Demand[StartNode - 1] <= RemainingCapacity){
										//The arc (StartNode,EndNode) can be used, add it.
										ijCurrentInArcs.push_back(*CurrentArc);
										jiCurrentInArcs.push_back(*CurrentArc);
									}

								}

							}//end go over all arcs into EndNode
						}//end else MaxNumberOfArcs>2

					}//end check whether the current end node is reachable and MaxNumberOfArcs>1

					//Update EndNode
					++EndNode;

				}//end go over all locations as end node of the arc

				//Initialize Labels:
				if (AlgorithmMembers.TypeOfExactColumnGenerator == 1){
					//Initialize LayeredLabels:

					//Obtain LayeredLabels
					ijPricingProblem.LayeredLabels = new std::list<Label>*[p];
					jiPricingProblem.LayeredLabels = new std::list<Label>*[p];
					std::list<Label>** ijLayeredLabels = ijPricingProblem.LayeredLabels;
					std::list<Label>** jiLayeredLabels = jiPricingProblem.LayeredLabels;
					//Initialize source label
					ijLayeredLabels[0] = new std::list<Label>[1];
					jiLayeredLabels[0] = new std::list<Label>[1];
					//Initalize labels for layers 1 to p
					for (short int Layer = 1; Layer < p; ++Layer){
						ijLayeredLabels[Layer] = new std::list<Label>[NumberOfCustomers];
						jiLayeredLabels[Layer] = new std::list<Label>[NumberOfCustomers];
					}
				}
				else if (AlgorithmMembers.TypeOfExactColumnGenerator == 2){
					//Initialize DoubleLayeredLabels:

					//Obtain integer capacity
					int IntegerCapacity = int(CurrentCapacity);
					//Initialize DoubleLayeredLabels
					ijPricingProblem.DoubleLayeredLabels = new std::list<Label>**[IntegerCapacity+1];
					jiPricingProblem.DoubleLayeredLabels = new std::list<Label>**[IntegerCapacity+1];
					std::list<Label>*** ijDoubleLayeredLabels = ijPricingProblem.DoubleLayeredLabels;
					std::list<Label>*** jiDoubleLayeredLabels = jiPricingProblem.DoubleLayeredLabels;
					//Initialize source label
					ijDoubleLayeredLabels[0] = new std::list<Label>*[1];
					jiDoubleLayeredLabels[0] = new std::list<Label>*[1];
					ijDoubleLayeredLabels[0][0] = new std::list<Label>[1];
					jiDoubleLayeredLabels[0][0] = new std::list<Label>[1];
					//Initialize layers to 0 for all loads
					for (int q = 1; q <= IntegerCapacity; ++q){
						ijDoubleLayeredLabels[q] = 0;
						jiDoubleLayeredLabels[q] = 0;
					}
				}
				else if (AlgorithmMembers.TypeOfExactColumnGenerator == 3){
					
					//Compute the maximum number of forward and backward labels
					short int MaxNumberOfBackwardArcs = (p - 1) / 2;
					short int MaxNumberOfForwardArcs = p - 1 - MaxNumberOfBackwardArcs;
					//Set maximum number of forward and backward labels.
					ijPricingProblem.MaxNumberOfForwardArcs = MaxNumberOfForwardArcs;
					jiPricingProblem.MaxNumberOfForwardArcs = MaxNumberOfForwardArcs;
					ijPricingProblem.MaxNumberOfBackwardArcs = MaxNumberOfBackwardArcs;
					jiPricingProblem.MaxNumberOfBackwardArcs = MaxNumberOfBackwardArcs;
					//Set the number of forward and backward layers
					ijPricingProblem.NumberOfForwardLayers = MaxNumberOfForwardArcs + 1;
					jiPricingProblem.NumberOfForwardLayers = MaxNumberOfForwardArcs + 1;
					ijPricingProblem.NumberOfBackwardLayers = MaxNumberOfBackwardArcs + 1;
					jiPricingProblem.NumberOfBackwardLayers = MaxNumberOfBackwardArcs + 1;
					
					//Initialize bi-directional labels:
					
					//Obtain LayeredForwardLabels
					ijPricingProblem.LayeredForwardLabels = new std::list<ForwardLabel>*[MaxNumberOfForwardArcs+1];
					jiPricingProblem.LayeredForwardLabels = new std::list<ForwardLabel>*[MaxNumberOfForwardArcs+1];
					std::list<ForwardLabel>** ijLayeredForwardLabels = ijPricingProblem.LayeredForwardLabels;
					std::list<ForwardLabel>** jiLayeredForwardLabels = jiPricingProblem.LayeredForwardLabels;
					//Initialize source label
					ijLayeredForwardLabels[0] = new std::list<ForwardLabel>[1];
					jiLayeredForwardLabels[0] = new std::list<ForwardLabel>[1];
					//Initalize labels for layers 1 to MaxNumberOfForwardArcs
					for (short int Layer = 1; Layer <= MaxNumberOfForwardArcs; ++Layer){
						ijLayeredForwardLabels[Layer] = new std::list<ForwardLabel>[NumberOfCustomers];
						jiLayeredForwardLabels[Layer] = new std::list<ForwardLabel>[NumberOfCustomers];
					}

					//Obtain LayeredBackwardLabels
					ijPricingProblem.LayeredBackwardLabels = new std::list<BackwardLabel>*[MaxNumberOfBackwardArcs+1];
					jiPricingProblem.LayeredBackwardLabels = new std::list<BackwardLabel>*[MaxNumberOfBackwardArcs+1];
					std::list<BackwardLabel>** ijLayeredBackwardLabels = ijPricingProblem.LayeredBackwardLabels;
					std::list<BackwardLabel>** jiLayeredBackwardLabels = jiPricingProblem.LayeredBackwardLabels;
					//Initialize source label
					ijLayeredBackwardLabels[0] = new std::list<BackwardLabel>[1];
					jiLayeredBackwardLabels[0] = new std::list<BackwardLabel>[1];
					//Initalize labels for layers 1 to MaxNumberOfBackwardArcs
					for (short int Layer = 1; Layer <= MaxNumberOfBackwardArcs; ++Layer){
						ijLayeredBackwardLabels[Layer] = new std::list<BackwardLabel>[NumberOfCustomers];
						jiLayeredBackwardLabels[Layer] = new std::list<BackwardLabel>[NumberOfCustomers];
					}
					
				}
				else if (AlgorithmMembers.TypeOfExactColumnGenerator == 4){

					//Compute the maximum number of forward and backward labels
					short int MaxNumberOfBackwardArcs = (p - 1) / 2;
					short int MaxNumberOfForwardArcs = p - 1 - MaxNumberOfBackwardArcs;
					//Set maximum number of forward and backward labels.
					ijPricingProblem.MaxNumberOfForwardArcs = MaxNumberOfForwardArcs;
					jiPricingProblem.MaxNumberOfForwardArcs = MaxNumberOfForwardArcs;
					ijPricingProblem.MaxNumberOfBackwardArcs = MaxNumberOfBackwardArcs;
					jiPricingProblem.MaxNumberOfBackwardArcs = MaxNumberOfBackwardArcs;
					//Set the number of forward and backward layers
					ijPricingProblem.NumberOfForwardLayers = p;
					jiPricingProblem.NumberOfForwardLayers = p;
					ijPricingProblem.NumberOfBackwardLayers = p;
					jiPricingProblem.NumberOfBackwardLayers = p;

					//Initialize bi-directional labels and completion bounds:

					//Obtain LayeredForwardLabels
					ijPricingProblem.LayeredForwardLabels = new std::list<ForwardLabel>*[p];
					jiPricingProblem.LayeredForwardLabels = new std::list<ForwardLabel>*[p];
					std::list<ForwardLabel>** ijLayeredForwardLabels = ijPricingProblem.LayeredForwardLabels;
					std::list<ForwardLabel>** jiLayeredForwardLabels = jiPricingProblem.LayeredForwardLabels;
					//Initialize source label
					ijLayeredForwardLabels[0] = new std::list<ForwardLabel>[1];
					jiLayeredForwardLabels[0] = new std::list<ForwardLabel>[1];
					//Obtain ForwardCompletionBounds
					ijPricingProblem.ForwardCompletionBounds = new double*[p];
					jiPricingProblem.ForwardCompletionBounds = new double*[p];
					double** ijForwardCompletionBounds = ijPricingProblem.ForwardCompletionBounds;
					double** jiForwardCompletionBounds = jiPricingProblem.ForwardCompletionBounds;
					//Initalize labels and bounds for layers 1 to p
					for (short int Layer = 1; Layer < p; ++Layer){
						ijLayeredForwardLabels[Layer] = new std::list<ForwardLabel>[NumberOfCustomers];
						jiLayeredForwardLabels[Layer] = new std::list<ForwardLabel>[NumberOfCustomers];
						ijForwardCompletionBounds[Layer] = new double[NumberOfCustomers];
						jiForwardCompletionBounds[Layer] = new double[NumberOfCustomers];
					}

					//Obtain LayeredBackwardLabels
					ijPricingProblem.LayeredBackwardLabels = new std::list<BackwardLabel>*[p];
					jiPricingProblem.LayeredBackwardLabels = new std::list<BackwardLabel>*[p];
					std::list<BackwardLabel>** ijLayeredBackwardLabels = ijPricingProblem.LayeredBackwardLabels;
					std::list<BackwardLabel>** jiLayeredBackwardLabels = jiPricingProblem.LayeredBackwardLabels;
					//Initialize source label
					ijLayeredBackwardLabels[0] = new std::list<BackwardLabel>[1];
					jiLayeredBackwardLabels[0] = new std::list<BackwardLabel>[1];
					//Obtain BackwardCompletionBounds
					ijPricingProblem.BackwardCompletionBounds = new double*[p];
					jiPricingProblem.BackwardCompletionBounds = new double*[p];
					double** ijBackwardCompletionBounds = ijPricingProblem.BackwardCompletionBounds;
					double** jiBackwardCompletionBounds = jiPricingProblem.BackwardCompletionBounds;
					//Initalize labels and bounds for layers 1 to p
					for (short int Layer = 1; Layer < p; ++Layer){
						ijLayeredBackwardLabels[Layer] = new std::list<BackwardLabel>[NumberOfCustomers];
						jiLayeredBackwardLabels[Layer] = new std::list<BackwardLabel>[NumberOfCustomers];
						ijBackwardCompletionBounds[Layer] = new double[NumberOfCustomers];
						jiBackwardCompletionBounds[Layer] = new double[NumberOfCustomers];
					}

				}

			}//end check whether the current pricing problem is relevant.

		}//end go over all customers j
	}//end go over all customers i

	//clean:
	//------
	delete[]SmallestDemand;
	delete[]SmallestDemandLocations;

}

//Initialize PStep_Adders
void PStep_Algorithm::InitializePStep_Adders(){
	
	//Check which p-step adder needs to be used
	if (AlgorithmMembers.AlgorithmSettings.AddPStepsSetting == 1){
		ListOfAllPStepAdders.emplace_front(AlgorithmMembers);
		CurrentPStepAdder = &ListOfAllPStepAdders.front();
	}

}

//Initialize PStep_PricingProblemSolvers
void PStep_Algorithm::InitializePStep_PricingProblemSolvers(){

	//Check which solver needs to be used
	if (AlgorithmMembers.AlgorithmSettings.PricingProblemSelection == 1){
		ListOfAllPricingProblemSequentialSolvers.emplace_front();
		CurrentPricingProblemSolver = &ListOfAllPricingProblemSequentialSolvers.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.PricingProblemSelection == 2){
		ListOfAllPricingProblemSequentialSolverSavingsTimers.emplace_front();
		CurrentPricingProblemSolver = &ListOfAllPricingProblemSequentialSolverSavingsTimers.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.PricingProblemSelection == 3){
		ListOfAllPricingProblemParallelSolvers.emplace_front();
		CurrentPricingProblemSolver = &ListOfAllPricingProblemParallelSolvers.front();
	}

}

//Initialize PStep_ColumnGenerators
void PStep_Algorithm::InitializePStep_ColumnGenerators(){

	//Initialize generic settings
	AlgorithmMembers.LimitedArcsUsed = false;
	AlgorithmMembers.DSSRUsed = false;

	//Check which column generator needs to be used.
	if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 1){
		AlgorithmMembers.TypeOfExactColumnGenerator = 1;
		ListOfColumnGenerator1s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator1s.front();
		
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 2){
		AlgorithmMembers.TypeOfExactColumnGenerator = 1;
		ListOfColumnGenerator2s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator2s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 3){
		AlgorithmMembers.TypeOfExactColumnGenerator = 1;
		ListOfColumnGenerator3s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator3s.front();
		AlgorithmMembers.LimitedArcsUsed = true;
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 4){
		AlgorithmMembers.TypeOfExactColumnGenerator = 1;
		ListOfColumnGenerator4s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator4s.front();
		AlgorithmMembers.LimitedArcsUsed = true;
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 5){
		AlgorithmMembers.TypeOfExactColumnGenerator = 1;
		ListOfColumnGenerator5s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator5s.front();
		AlgorithmMembers.DSSRUsed = true;
		//Check whether ng-path relaxation is not applied
		if (!(AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation)){
			//Set neighborhood size to n and apply ng-path relaxation
			//This is required to apply decremental state-space relaxation.
			AlgorithmMembers.AlgorithmSettings.NeighbourhoodSize = AlgorithmMembers.Instance.NumberOfCustomers - 1;
			AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation = true;
		}
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 6){
		AlgorithmMembers.TypeOfExactColumnGenerator = 1;
		ListOfColumnGenerator6s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator6s.front();
		AlgorithmMembers.DSSRUsed = true;
		//Check whether ng-path relaxation is not applied
		if (!(AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation)){
			//Set neighborhood size to n and apply ng-path relaxation
			//This is required to apply decremental state-space relaxation.
			AlgorithmMembers.AlgorithmSettings.NeighbourhoodSize = AlgorithmMembers.Instance.NumberOfCustomers - 1;
			AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation = true;
		}
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 7){
		AlgorithmMembers.TypeOfExactColumnGenerator = 1;
		ListOfColumnGenerator7s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator7s.front();
		AlgorithmMembers.LimitedArcsUsed = true;
		AlgorithmMembers.DSSRUsed = true;
		//Check whether ng-path relaxation is not applied
		if (!(AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation)){
			//Set neighborhood size to n and apply ng-path relaxation
			//This is required to apply decremental state-space relaxation.
			AlgorithmMembers.AlgorithmSettings.NeighbourhoodSize = AlgorithmMembers.Instance.NumberOfCustomers - 1;
			AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation = true;
		}
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 8){
		AlgorithmMembers.TypeOfExactColumnGenerator = 1;
		ListOfColumnGenerator8s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator8s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 9){
		AlgorithmMembers.TypeOfExactColumnGenerator = 1;
		ListOfColumnGenerator9s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator9s.front();
		AlgorithmMembers.LimitedArcsUsed = true;
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 10){
		AlgorithmMembers.TypeOfExactColumnGenerator = 1;
		ListOfColumnGenerator10s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator10s.front();
		AlgorithmMembers.LimitedArcsUsed = true;
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 11){
		AlgorithmMembers.TypeOfExactColumnGenerator = 1;
		ListOfColumnGenerator11s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator11s.front();
		AlgorithmMembers.LimitedArcsUsed = true;
		AlgorithmMembers.DSSRUsed = true;
		//Check whether ng-path relaxation is not applied
		if (!(AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation)){
			//Set neighborhood size to n and apply ng-path relaxation
			//This is required to apply decremental state-space relaxation.
			AlgorithmMembers.AlgorithmSettings.NeighbourhoodSize = AlgorithmMembers.Instance.NumberOfCustomers - 1;
			AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation = true;
		}
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 12){
		AlgorithmMembers.TypeOfExactColumnGenerator = 2;
		ListOfColumnGenerator12s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator12s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 13){
		AlgorithmMembers.TypeOfExactColumnGenerator = 2;
		ListOfColumnGenerator13s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator13s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 14){
		AlgorithmMembers.TypeOfExactColumnGenerator = 3;
		ListOfColumnGenerator14s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator14s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 15){
		AlgorithmMembers.TypeOfExactColumnGenerator = 3;
		ListOfColumnGenerator15s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator15s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 16){
		AlgorithmMembers.TypeOfExactColumnGenerator = 3;
		ListOfColumnGenerator16s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator16s.front();
		AlgorithmMembers.LimitedArcsUsed = true;
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 17){
		AlgorithmMembers.TypeOfExactColumnGenerator = 3;
		ListOfColumnGenerator17s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator17s.front();
		AlgorithmMembers.LimitedArcsUsed = true;
		AlgorithmMembers.DSSRUsed = true;
		//Check whether ng-path relaxation is not applied
		if (!(AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation)){
			//Set neighborhood size to n and apply ng-path relaxation
			//This is required to apply decremental state-space relaxation.
			AlgorithmMembers.AlgorithmSettings.NeighbourhoodSize = AlgorithmMembers.Instance.NumberOfCustomers - 1;
			AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation = true;
		}
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 18){
		AlgorithmMembers.TypeOfExactColumnGenerator = 3;
		ListOfColumnGenerator18s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator18s.front();
		AlgorithmMembers.LimitedArcsUsed = true;
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 19){
		AlgorithmMembers.TypeOfExactColumnGenerator = 4;
		ListOfColumnGenerator19s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator19s.front();
		AlgorithmMembers.DSSRUsed = true;
		//Check whether ng-path relaxation is not applied
		if (!(AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation)){
			//Set neighborhood size to n and apply ng-path relaxation
			//This is required to apply decremental state-space relaxation.
			AlgorithmMembers.AlgorithmSettings.NeighbourhoodSize = AlgorithmMembers.Instance.NumberOfCustomers - 1;
			AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation = true;
		}
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 20){
		AlgorithmMembers.TypeOfExactColumnGenerator = 4;
		ListOfColumnGenerator20s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator20s.front();
		AlgorithmMembers.LimitedArcsUsed = true;
		AlgorithmMembers.DSSRUsed = true;
		//Check whether ng-path relaxation is not applied
		if (!(AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation)){
			//Set neighborhood size to n and apply ng-path relaxation
			//This is required to apply decremental state-space relaxation.
			AlgorithmMembers.AlgorithmSettings.NeighbourhoodSize = AlgorithmMembers.Instance.NumberOfCustomers - 1;
			AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation = true;
		}
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnGenerator == 21){
		AlgorithmMembers.TypeOfExactColumnGenerator = 4;
		ListOfColumnGenerator21s.emplace_front(AlgorithmMembers, *CurrentPStepAdder, *CurrentPricingProblemSolver);
		CurrentColumnGenerator = &ListOfColumnGenerator21s.front();
		AlgorithmMembers.LimitedArcsUsed = true;
		AlgorithmMembers.DSSRUsed = true;
		//Check whether ng-path relaxation is not applied
		if (!(AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation)){
			//Set neighborhood size to n and apply ng-path relaxation
			//This is required to apply decremental state-space relaxation.
			AlgorithmMembers.AlgorithmSettings.NeighbourhoodSize = AlgorithmMembers.Instance.NumberOfCustomers - 1;
			AlgorithmMembers.AlgorithmSettings.UseNGPathRelaxation = true;
		}
	}

}

//Initialize PStep_ColumnManagers
void PStep_Algorithm::InitializePStep_ColumnManagers(){

	//Check which column generator needs to be used.
	if (AlgorithmMembers.AlgorithmSettings.ColumnManager == 0){
		ListOfColumnManagers.emplace_front();
		CurrentColumnManager = &ListOfColumnManagers.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnManager == 1){
		ListOfColumnManager1s.emplace_front(AlgorithmMembers);
		CurrentColumnManager = &ListOfColumnManager1s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnManager == 2){
		ListOfColumnManager2s.emplace_front(AlgorithmMembers);
		CurrentColumnManager = &ListOfColumnManager2s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.ColumnManager == 3){
		ListOfColumnManager3s.emplace_front(AlgorithmMembers);
		CurrentColumnManager = &ListOfColumnManager3s.front();
	}

}

//Initialize PStep_SeparationAlgorithms
void PStep_Algorithm::InitializePStep_SeparationAlgorithms(){
	
	//Check which manager should be used
	if (AlgorithmMembers.AlgorithmSettings.SeparationStrategy == 0){
		ListOfSeperationManager0s.emplace_front(AlgorithmMembers,CurrentSeparationAlgorithm);
		CurrentSeparationManager = &ListOfSeperationManager0s.front();		
	}
	else if (AlgorithmMembers.AlgorithmSettings.SeparationStrategy == 1){
		ListOfSeperationManager1s.emplace_front(AlgorithmMembers, CurrentSeparationAlgorithm);
		CurrentSeparationManager = &ListOfSeperationManager1s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.SeparationStrategy == 2){
		ListOfSeperationManager2s.emplace_front(AlgorithmMembers, CurrentSeparationAlgorithm);
		CurrentSeparationManager = &ListOfSeperationManager2s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.SeparationStrategy == 3){
		ListOfSeperationManager3s.emplace_front(AlgorithmMembers, CurrentSeparationAlgorithm);
		CurrentSeparationManager = &ListOfSeperationManager3s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.SeparationStrategy == 4){
		ListOfSeperationManager4s.emplace_front(AlgorithmMembers, CurrentSeparationAlgorithm);
		CurrentSeparationManager = &ListOfSeperationManager4s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.SeparationStrategy == 5){
		ListOfSeperationManager5s.emplace_front(AlgorithmMembers, CurrentSeparationAlgorithm);
		CurrentSeparationManager = &ListOfSeperationManager5s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.SeparationStrategy == 6){
		ListOfSeperationManager6s.emplace_front(AlgorithmMembers, CurrentSeparationAlgorithm);
		CurrentSeparationManager = &ListOfSeperationManager6s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.SeparationStrategy == 7){
		ListOfSeperationManager7s.emplace_front(AlgorithmMembers, CurrentSeparationAlgorithm);
		CurrentSeparationManager = &ListOfSeperationManager7s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.SeparationStrategy == 8){
		ListOfSeperationManager8s.emplace_front(AlgorithmMembers, CurrentSeparationAlgorithm);
		CurrentSeparationManager = &ListOfSeperationManager8s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.SeparationStrategy == 9){
		ListOfSeperationManager9s.emplace_front(AlgorithmMembers, CurrentSeparationAlgorithm);
		CurrentSeparationManager = &ListOfSeperationManager9s.front();
	}
	else if (AlgorithmMembers.AlgorithmSettings.SeparationStrategy == 10){
		ListOfSeperationManager10s.emplace_front(AlgorithmMembers, CurrentSeparationAlgorithm);
		CurrentSeparationManager = &ListOfSeperationManager10s.front();
	}
	
}

//Demand functions:
//------------------------------------------------------------------------------------

//Compute the vehicle capacity minus the demand of the NumberOfInclusions first other locations
//excluding all customers in the list of ExclusiveCustomers.
double PStep_Algorithm::ComputeCapacity(double Vehiclecapacity, double* DemandList, double* DemandCustomerList, short int NumberOfInclusions, std::vector<short int>& ExclusiveCustomers, short int NumberOfCustomers, short int& CustomerIndex){

	//Initialize output
	double Output = Vehiclecapacity;

	//Initialize counter to count up to NumberOfInclusions processed customers
	short int Counter = 0;
	//Until p customers are processed, consider the next smallest demand
	while (Counter < NumberOfInclusions && CustomerIndex<NumberOfCustomers){
		//Check whether the next smallest customer is j
		bool CustomerIsValid = true;
		//Check whether the current customer is included in the exclusive list.
		for (auto ExclusiveCustomer = ExclusiveCustomers.begin(); ExclusiveCustomer != ExclusiveCustomers.end() && CustomerIsValid; ++ExclusiveCustomer){
			if (DemandCustomerList[CustomerIndex] == *ExclusiveCustomer){
				CustomerIsValid = false;
			}
		}
		//Check whether the current customer was not included in the exclusive list
		if (CustomerIsValid){
			//Update tempCapacity and counter
			Output -= DemandList[CustomerIndex];
			++Counter;
		}
		++CustomerIndex;
	}

	return Output;
}


double PStep_Algorithm::ComputeCapacity(double Vehiclecapacity, double* DemandList, double* DemandCustomerList, short int NumberOfInclusions, std::vector<short int>& ExclusiveCustomers, short int NumberOfCustomers){

	short int DummyIndex = 0;
	return ComputeCapacity(Vehiclecapacity, DemandList, DemandCustomerList, NumberOfInclusions, ExclusiveCustomers, NumberOfCustomers, DummyIndex);

}