/*Implementation of PStep_SeparationAlgorithm.h

It contains the member functions of PStep_SeparationAlgorithm3

*/

#include "PStep_SeparationAlgorithm.h"

//PStep_SeparationAlgorithm3:
//>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

//Constructor functions
//------------------------------------------------------------------------------------

//Basic constructor
PStep_SeparationAlgorithm3::PStep_SeparationAlgorithm3(PStep_AlgorithmMembers& InitAlgorithmMembers)
	:PStep_SeparationAlgorithm2(InitAlgorithmMembers)
{
	//Obtain number of customers
	short int NumberOfCustomers = AlgorithmMembers.Instance.NumberOfCustomers;
	//Obtain Demand
	double* Demand = AlgorithmMembers.Instance.Demand;
	
	//Compute QMin:
	//-------------

	//Obtain upper bound on the number of vehicles
	short int NumOfVehiclesUB;
	//Check whether the number of vehicles is fixed
	if (AlgorithmMembers.Instance.NumberOfVehiclesIsFixed){
		//Use the fixed number of vehicles as UB.
		NumOfVehiclesUB = AlgorithmMembers.Instance.FixedNumberOfVehicles;
	}
	else{
		//Use the number of customers as upper bound on the number of vehicles.
		NumOfVehiclesUB = AlgorithmMembers.Instance.NumberOfCustomers;
	}

	//Compute total  integer demand
	int TotalDemand = 0;
	for (short int i = 0; i<NumberOfCustomers; i++){
		TotalDemand += int(Demand[i]);
	}

	//Compute QMin
	QMin = TotalDemand - (NumOfVehiclesUB - 1)*int(AlgorithmMembers.Instance.VehicleCapacity);

}

//Separate
//------------------------------------------------------------------------------------

//Separate capacity cuts using the algorithm of Lysgaard.
//These cuts are stored in AlgorithmMembers.CapacityCuts, new cuts are appended to the list.
//If it fails, separate strengthened comb cuts using the algorithm of Lysgaard.
//These cuts are stored in AlgorithmMembers.StrengthenedCombCuts, new cuts are appended to the list.
//If that also fails, separate 2 edge hypotour cuts using the algorithm of Lysgaard.
//These cuts are stored in AlgorithmMembers.TwoEdgeHypotourCuts, new cuts are appended to the list.
void PStep_SeparationAlgorithm3::Separate(){

	//Prepare the selected edges
	SetCurrentEdges();

	//Check whether there are any positively selected edges
	if (NumberOfEdgesPlusOne > 1){

		//Obtain current number of capacity cuts
		int NumberOfCCBeforeSeparation = AlgorithmMembers.NumberOfCapacityCuts;

		//Separate capacity cuts
		SeparateCCLysgaard();

		//Check whether new capacity cuts were added
		if (NumberOfCCBeforeSeparation == AlgorithmMembers.NumberOfCapacityCuts){
			//No violated capacity cuts were identified, separate strengthened comb cuts

			//Obtain current number of strentghened comb cuts
			int NumberOfSCombBeforeSeparation = AlgorithmMembers.NumberOfSCombCuts;

			//Separate strengthened comb cuts
			SeparateSCombLysgaard();

			//Check whether new strengthened comb cuts were added
			if (NumberOfSCombBeforeSeparation == AlgorithmMembers.NumberOfSCombCuts){
				//No violated strengthened comb cuts were identified, separate 2 edge hypotour cuts

				//Separate 2 edge hypotour cuts
				Separate2EHLysgaard();

			}

		}

	}

}

//Strengthened comb separation functions:
//------------------------------------------------------------------------------------

//Separate strengthened comb cuts using the algorithm of Lysgaard.
//These cuts are stored in AlgorithmMembers.StrengthenedCombCuts, new cuts are appended to the list.
void PStep_SeparationAlgorithm3::SeparateSCombLysgaard(){

	//Obtain start time of the separation algorithm
	timePoint tic = HomeMadeclock::now();

	//Initialize
	//-------------------------------------------------

	//Initialize support variables
	short int NumberOfCustomers = AlgorithmMembers.Instance.NumberOfCustomers;
	short int NumberOfLocations = NumberOfCustomers + 1;
	short int TotalNumberOfEdges = NumberOfLocations*(NumberOfLocations - 1) / 2;
	double Tolerance = AlgorithmMembers.Tolerance;
	std::vector < std::vector < PStep_Arc >> &Arcs = AlgorithmMembers.Arcs;
	double MinimumViolation = AlgorithmMembers.AlgorithmSettings.MinimumSCombViolation;

	//Obtain the maximum number of strengthened comb cuts added per iteration
	int MaxCuts = AlgorithmMembers.AlgorithmSettings.MaxSCombPerIteration;

	//Initalize max violation
	double MaxViolation;
	//Initialize CnstrMgrPointer to contain the generated cuts
	CnstrMgrPointer MyCutsCMP;
	CMGR_CreateCMgr(&MyCutsCMP, MaxCuts);

	//Use Lysgaard's algorithm:
	//-------------------------------------------------
	COMBSEP_SeparateCombs(NumberOfCustomers, LysgaardDemand, int(AlgorithmMembers.Instance.VehicleCapacity), QMin, NumberOfEdgesPlusOne - 1, EdgeTail, EdgeHead, EdgeX, MaxCuts, &MaxViolation, MyCutsCMP);
	
	//Process new cuts
	//-------------------------------------------------

	//Initialize indicators, indicating whether an edge is in use.
	//EdgeInUse[(j-1)*j/2+i] is corresponds to edge (i,j) for i<j.
	bool* EdgeInUse = new bool[TotalNumberOfEdges];
	//Initialize coefficients
	double* Coefficient = new double[TotalNumberOfEdges];

	//Initialize indicators, indicating whether a location is in a tooth
	bool* LocationInTooth = new bool[NumberOfLocations];

	//Go over all new cuts
	for (int i = 0; i<MyCutsCMP->Size; i++){

		//Initialize strengthened comb cut
		AlgorithmMembers.StrengthenedCombCuts.emplace_back();
		//Obtain new cut
		PStep_EdgeAndCoefficientCut &NewSCombCut = AlgorithmMembers.StrengthenedCombCuts.back();
		//Update counter
		++AlgorithmMembers.NumberOfSCombCuts;

		//Reset EdgeInUse and 
		for (short int j = 0; j < TotalNumberOfEdges; ++j){
			EdgeInUse[j] = false;
		}

		//Obtain number of teeth
		short int NoOfTeeth = MyCutsCMP->CPL[i]->Key;
		
		//Go over all teeth (t=0 is the handle)
		for (short int t = 0; t <= NoOfTeeth; t++){

			//Reset LocationInTooth
			for (short int j = 0; j < NumberOfLocations; ++j){
				LocationInTooth[j] = false;
			}

			//Check whether we are reading the handle
			if (t == 0){
				//Read the handle

				//Go over all locations in the handle
				for (short int k = 1; k <= MyCutsCMP->CPL[i]->IntListSize; k++){
					//Obtain current location
					short int j = MyCutsCMP->CPL[i]->IntList[k];
					if (j == NumberOfLocations){
						//Recode the depot
						j = 0;
					}
					//Update Location In Tooth
					LocationInTooth[j] = true;
				}

			}
			else{
				//Read a tooth

				//Do some Lysgaard stuff to get the proper indices...
				short int MinIdx = MyCutsCMP->CPL[i]->ExtList[t];
				short int MaxIdx;
				if (t == NoOfTeeth){
					MaxIdx = MyCutsCMP->CPL[i]->ExtListSize;
				}
				else{
					MaxIdx = MyCutsCMP->CPL[i]->ExtList[t + 1] - 1;
				}
				//Go over all locations in the tooth
				for (short int k = MinIdx; k <= MaxIdx; k++){
					//Obtain current location
					short int j = MyCutsCMP->CPL[i]->ExtList[k];
					//Check whether the current location is the depot.
					if (j == NumberOfLocations){
						//Recode the depot
						j = 0;
					}
					//Update Location In Tooth
					LocationInTooth[j] = true;
				}

			}

			//Create a list of locations in and out of the tooth
			std::list<short int> LocationsInTooth;
			std::list<short int> LocationsNotInTooth;
			//Go over all locations and put them in the right list
			for (short int j = 0; j < NumberOfLocations; ++j){
				//Check whether the location is in the tooth
				if (LocationInTooth[j]){
					//The location is in the tooth, add it to the list
					LocationsInTooth.push_back(j);
				}
				else{
					//The location is not in the tooth, add it to the list
					LocationsNotInTooth.push_back(j);
				}
			}

			//Set edge in use and coefficients of the edges out of the current tooth:
			//Go over all locations in the tooth
			for (auto k = LocationsInTooth.begin(); k != LocationsInTooth.end(); ++k){
				//Go over all locations not in the tooth
				for (auto j = LocationsNotInTooth.begin(); j != LocationsNotInTooth.end(); ++j){					
					//Obtain the edge index of the corresponding edge
					short int Index = Arcs[*k][*j].Index;
					//Check whether this edge is already in use
					if (EdgeInUse[Index]){
						//The edge is already in use, update the coefficient
						Coefficient[Index] += 1;
					}
					else{
						//The edge is not yet in use, put it to use
						EdgeInUse[Index] = true;
						Coefficient[Index] = 1;
					}
				}
			}

		}//end go over all teeth
		
		//Add all edges and coefficients to the cut
		for (short int j = 0; j < TotalNumberOfEdges; ++j){
			//Check whether the edge is in use
			if (EdgeInUse[j]){
				//The edge is in use, add it
				NewSCombCut.EdgeAndCoefficients.emplace_back();
				NewSCombCut.EdgeAndCoefficients.back().EdgeIndex = j;
				NewSCombCut.EdgeAndCoefficients.back().Coefficient = Coefficient[j];
			}
		}

		//Add bound to the strengthened comb cut
		NewSCombCut.Bound = MyCutsCMP->CPL[i]->RHS;

		//Check whether the violation of the cut is large enough
		if (ComputeViolation(NewSCombCut, false) >= MinimumViolation){
			//The violation is large enough, so add the cut:
			
			//Move the current cut to LysgaardCuts
			CMGR_MoveCnstr(MyCutsCMP, AlgorithmMembers.LysgaardCuts, i, 0);
		}
		else{
			//Remove cut again
			AlgorithmMembers.StrengthenedCombCuts.pop_back();
			--AlgorithmMembers.NumberOfSCombCuts;
		}

	}

	//clean
	delete[]EdgeInUse;
	delete[]Coefficient;
	delete[]LocationInTooth;
	CMGR_FreeMemCMgr(&MyCutsCMP);

	//Terminate:
	//-------------------------------------------------
	//Obtain end time of the separation algorithm
	timePoint toc = HomeMadeclock::now();
	//Update strengthened comb separation time
	AlgorithmMembers.CurrentCGStatistics.TotalSeparateSCombTime += HomeMadeclock::seconds(tic, toc);

}