// ================= Lagrangean B&B by Filip Cakl (2010) ========================

#include "stdafx.h"

#ifdef _MSC_VER
#pragma hdrstop
#endif

using namespace std;

#include "../../Kernel/kernelincludes.h"
#include "math.h"

#include "CovMatrix.h"
#include "LagBaBCover.h"
#include "defaults.h"

CoverOpts::CoverOpts(void)
{
	lagrange_total_iter = DEFAULT_LAG_TOTALITER;
	lagrange_iter = DEFAULT_LAG_ITER;
	lagrange_iter_dual = DEFAULT_LAG_ITER_DUAL;
	lagrange_heur_iter = DEFAULT_LAG_HEURITER;
	lagrange_init_coef = DEFAULT_LAG_KOEFINIT;
	lagrange_min_coef = (float)DEFAULT_LAG_KOEFMIN;
	lagrange_init_coef_dual = DEFAULT_LAG_KOEFINIT_DUAL;
	lagrange_min_coef_dual = (float)DEFAULT_LAG_KOEFMIN_DUAL;
	lagrange_simple_dual = DEFAULT_LAG_SDUAL;

	run_lag = false;
	use_contribAdd=false;
	use_servit=false;

}

SCPView::SCPView(CovMatrix* _matrix,CoveringMatrix* matrix) {

	int i=0;

	CLink *ColRowLink=NULL;

	row_count=matrix->countOfRows;
	column_count=matrix->countOfCols;

	// Assign problem
	matrix_all = _matrix;

	// Create initial mapping
	column_map=new int[column_count];
	row_map=new int[row_count];

	// Mapping of rows
	ColRowLink= matrix->ColLink;

	//cout << "view columns: " << column_count << endl;
	for(i=0;i<column_count;i++) {
		//cout << i << ":" <<ColRowLink->ID << endl;
		column_map[i]=ColRowLink->ID;
		ColRowLink=ColRowLink->Next;
	}
	
	// Mapping of columns
	ColRowLink = matrix->RowLink;
	//cout << endl << "view rows: " << row_count << endl;
	for(i=0;i<row_count;i++) {
		//cout << i << ":" << ColRowLink->ID << endl;
		row_map[i]=ColRowLink->ID;
		ColRowLink=ColRowLink->Next;
	}
}

SCPView::~SCPView() {
	
	// Deallocate memory if possible
	if (column_map!=NULL) delete[] column_map;
	if (row_map!=NULL) delete[] row_map;
}

unsigned char SCPView::get(int row,int column) {

	int xr = row_map[row];
	int xc = column_map[column];

	// Get the value
	return matrix_all->GetValue(xr,xc);
}


int SCPView::get_cost(int column) {

	return matrix_all->GetCost(column_map[column]);
}

void SCPView::print() {

	int i=0;
	int j=0;

	// Print general information
	cout << "SCPView {" << "R:" << row_count << ",C:" << column_count << "}" << endl;
	
	// Print matrix
	for(i=0;i<row_count;i++) {
		for(j=0;j<column_count;j++) {
			cout << get(i,j);
		}
		cout << endl;
	}
	
	// Print costs
	cout << "Costs:";
	for(j=0;j<column_count;j++) {

		cout << get_cost(j) << " ";
	}
	cout << endl;
}

LagrangeSolver::LagrangeSolver(SCPView* pv) {

	int rows = pv->row_count;
	int columns = pv->column_count;

	// Allocate space
	t=new float[rows];
	n=new float[columns];
	primal_delta=new unsigned char[rows];
	dual_delta=new unsigned char[columns];
	x=new unsigned char[columns];
	x_heur=new unsigned char[columns];
	x_best=new unsigned char[columns];
	y=new float[rows];
		
	help_indexes = new int[columns];
	help_costs=new float[columns];
	help_number_of_ones=new int[columns];

	row_is_covered=new bool[rows];
	column_is_recounted=new unsigned char[columns];
	row_is_recounted=new unsigned char[rows];

	costs=new float[columns];
	reduced_costs=new float[columns];
	modified_reduced_costs=new float[columns];

	min_costs=new float[rows];
	dual_reduced_costs=new float[rows];
	
	primal_slack=new int[rows];
	dual_slack=new float[columns];

	problemView=pv;

	rows_covered_by_column = new int*[columns];
	rows_covered_by_column_count = new int[columns];
	columns_covered_by_row = new int*[rows];
	columns_covered_by_row_count = new int[rows];

}

LagrangeSolver::~LagrangeSolver() {

	int i;

	// Deallocate memory if possible
	if (t!=NULL) delete[] t;
	if (n!=NULL) delete[] n;
	if (primal_delta!=NULL) delete[] primal_delta;
	if (dual_delta!=NULL) delete[] dual_delta;
	if (x!=NULL) delete[] x;
	if (y!=NULL) delete[] y;
	if (x_heur!=NULL) delete[] x_heur;
	if (x_best!=NULL) delete[] x_best;
	if (help_indexes!=NULL) delete[] help_indexes;
	if (help_costs!=NULL) delete[] help_costs;
	if (help_costs!=NULL) delete[] help_number_of_ones;
	if (row_is_covered!=NULL) delete[] row_is_covered;
	if (column_is_recounted!=NULL) delete[] column_is_recounted;
	if (row_is_recounted!=NULL) delete[] row_is_recounted;
	if (costs!=NULL) delete[] costs;
	if (reduced_costs!=NULL) delete[] reduced_costs;
	if (min_costs!=NULL) delete[] min_costs;
	if (dual_reduced_costs!=NULL) delete[] dual_reduced_costs;
	if (modified_reduced_costs!=NULL) delete[] modified_reduced_costs;
	if (primal_slack!=NULL) delete[] primal_slack;
	if (dual_slack!=NULL) delete[] dual_slack;

	for(i=0; i<problemView->column_count; i++) {
		if (rows_covered_by_column[i]!=NULL) 
			delete[]rows_covered_by_column[i];
	}
	if (rows_covered_by_column!=NULL)  delete[] rows_covered_by_column;
	if (rows_covered_by_column_count!=NULL) delete[]rows_covered_by_column_count;

	for(i=0; i<problemView->row_count; i++) {
		if (columns_covered_by_row[i]!=NULL) 
			delete[]columns_covered_by_row[i];
	}
	if (columns_covered_by_row!=NULL) delete[]columns_covered_by_row;
	if (columns_covered_by_row_count!=NULL) delete[]columns_covered_by_row_count;

}

void LagrangeSolver::solve(int iterations,int randomizeCosts, int heuristic, int heurIterations, float stepCoef, float stepCoefMin, int stepCoefIter,int maxLowerBound) {

	int i,j;
	int hi;
	int removeType;
	float costs;
	float minReducedCost=0.;

	if(heuristic==0 || heuristic == 2) {
		removeType = 0;
	} else if (heuristic == 1 || heuristic == 3) {
		removeType = 1; 
	} else {
		removeType = 2;
	}

	// Initialize
	
	//cout << "initializing " << endl;
	initialize(randomizeCosts,stepCoef,stepCoefMin,stepCoefIter);
	//cout << "initialized " << endl;

	hi=0;
	for(i=0;i<iterations;i++) {


		// Solve LP with current set of multipliers, update ZLP,ZLB and x
		solve_lp();

		if(hi>=heurIterations) {
			// Solve SCP with some heuristic, update ZUB and best solution found
			if(heuristic<2) {
				ha(removeType);
			} else {
				//cout << "HA_modified" << endl;
				ha_modified(removeType);
			}

			hi=0;
		} else {
			hi++;
		}

		// If (ceiled lower bound = cost of heuristic solution) then stop - found solution is optimal
		if (lbound == best_heur_value) {
			//cout << "breaking - UB reached" << endl;
			break;
		};
		
		// If ceiled lower bound == minLB then stop - cannot find solution usable by B&B alghoritm
		if (lbound >= maxLowerBound) {
			//cout <<"breaking - LB reached" << endl;
			break;
		};

		// Calculate subgradient and update Lagrangean multipliers
		if(compute_subgradient()==0){
			 //cout << "breaking - subgradient too small" << endl;
			break;
		}

		hi++;
	}


	// Find column with lowest reduced cost
	j=0;

	//cout << "going to select recommended_column " << endl;
	for(i=0;i<problemView->column_count;i++) {

		if(x_best[i]==1) {
			if(j==0) {
				//cout << "recommended_column is now " <<  problemView->column_map[i] << endl;
				recommended_column = problemView->column_map[i];
				minReducedCost=get_primal_reduced_cost(i);
				j=1;
			}

			costs = get_primal_reduced_cost(i);

			if(costs<minReducedCost) {
				minReducedCost=costs;
				//cout << "recommended_column is now " <<  problemView->column_map[i] << endl;
				recommended_column = problemView->column_map[i];
			}
		}
	}

}

void LagrangeSolver::solve_dual(int iterations,int heuristic, int heurIterations,
								float stepCoef, float stepCoefMin,int stepCoefIter, 
								int maxLowerBound,
								float stepCoefDual,  float stepCoefMinDual,int stepCoefIterDual) {

	int i,j,hi;
	int removeType;
	float costs;
	float minReducedCost=0.;

	if(heuristic==0 || heuristic == 2) {
		removeType = 0;
	} else if (heuristic == 1 || heuristic == 3) {
		removeType = 1; 
	} else {
		removeType = 2;
	}

	// Initialize
	
	initialize(1,stepCoef,stepCoefMin,stepCoefIter);

	// First iteration of alghoritm
	// Compute (primal and dual) slacks, costs and multipliers
	// Also update ZLB and ZUB
	solve_dlp();

	//cout << "ZLP :" << setw(8) <<  ZLP <<  " LB: " << setw(8) << ZLB ;

	if(heuristic<2) {
		ha(removeType);
	} else {
		ha_modified(removeType);
	}
	
	//cout << " heur :" << setw(8) <<  heur_value <<  " bestHeur: " << setw(8) << (float)best_heur_value << " UB: " <<  setw(8) << ZUB << endl  ;

	//cout << "solved dlp " << endl;
	
	/*if(covered_rows < problemView->row_count) {

		if(heuristic<2) {
			ha(removeType);
		} else {
			ha_modified(removeType);
		}

	} else {

		// All rows are covered just by loving LP
		// Update upper bound and best solution found
		costs=0.0f;
		for(j=0;j<problemView->column_count;j++) {
			if (x[j]==1) {
				costs+=problemView->get_cost(j);
			}
		}
		heur_value = costs;
		update_best_solution(costs);
	}*/

	// If (ceiled lower bound = cost of heuristic solution) then stop - found solution is optimal
	if (lbound == best_heur_value) {
		//cout << "reached LB ceiling" << endl;
		return;
	};

	// If ceiled lower bound == minLB then stop - cannot find solution usable by B&B alghoritm
	if (lbound >= maxLowerBound && maxLowerBound>0 ) {
		//cout << "reached min LB" << endl;
		return;
	};

	//cout << "lbound: " << setw(8) << lbound << " bheurval: " << setw(8) << best_heur_value << endl;

	hi=heurIterations;
	for(i=1;i<iterations;i++) {

		// Next iterations of alghoritm
		// Update (primal and dual) slacks, costs and multipliers
		// Also update ZLB and ZUB
		//cout << "solving DLP next" << endl;
		//cout << "ZLP :" << setw(8) <<  ZLP <<  " LB: " << setw(8) << ZLB ;

		solve_dlp_next();

		if(hi>=heurIterations) {
			// Solve SCP with some heuristic, update ZUB and best solution found
			if(heuristic<2) {
				ha(removeType);
			} else {
				ha_modified(removeType);
			}

			hi=0;
		}

		//cout << " heur :" << setw(8) <<  heur_value <<  " bestHeur: " << setw(8) << (float)best_heur_value << " UB: " <<  setw(8) << ZUB  << endl  ;
		
		/*
		if(covered_rows < problemView->row_count) {

			// Some rows are still uncovered
			if(hi>=heurIterations) {
				// Solve SCP with some heuristic, update ZUB and best solution found
				if(heuristic<2) {
					ha(removeType);
				} else {
					ha_modified(removeType);
				}

				hi=0;
			}
		} else {

			// All rows are covered just by loving LP
			// Update upper bound and best solution found
			costs=0.0f;
			for(j=0;j<problemView->column_count;j++) {
				if (x[j]==1) {
					costs+=problemView->get_cost(j);
				}
			}
			heur_value = costs;
			update_best_solution(costs);
		}
		*/
	
		// If both primal and dual step coefficients are too small, stop the alghoritm
		if(dual_step_change_stop==1 && primal_step_change_stop==1){
			//cout << "both steps are too small" << endl;
			break;
		}

		// If (ceiled lower bound = cost of heuristic solution) then stop - found solution is optimal
		if (lbound == best_heur_value) {
			//cout << "reached LB ceiling" << endl;
			break;
		};

		// If ceiled lower bound == minLB then stop - cannot find solution usable by B&B alghoritm
		if (lbound >= maxLowerBound) {
			//cout << "reached min LB" << endl;
			break;
		};

		hi++;
	}
	
	// Find column with lowest reduced cost
	j=0;
	for(i=0;i<problemView->column_count;i++) {

		if(x_best[i]==1) {
			if(j==0) {
				recommended_column = problemView->column_map[i];
				minReducedCost=get_primal_reduced_cost(i);
				j=1;
			}

			costs = get_primal_reduced_cost(i);

			if(costs<minReducedCost) {
				minReducedCost=costs;
				recommended_column = problemView->column_map[i];
			}
		}
	}

	//cout << "ending solving dual relaxation" << endl;
}

void LagrangeSolver::initialize(int randomizeMultipliers, float stepCoef, float stepCoefMin, int stepCoefIter) {

	int j=0;
	int i=0;
	int cols=0;
	int rows=0;

	int* helpMap;
	int rowsCovered;
	int colsCovered;

	int minCost = PS_INFINITY;
	int maxCost = NG_INFINITY;
	int cost;

	cols=problemView->column_count;
	rows=problemView->row_count;

	// Iteration counters
	iterations=0;
	sup_iterations=0;
	sup_dual_iterations=0;

	// ZLB and ZLP
	ZLB=0.;
	sup_ZLB=0.;
	ZLP=0.0f;

	// Zero solution vectors
	memset(x,0,cols);
	memset(x_heur,0,cols);
	memset(x_best,0,cols);

	// Set number of covered rows to 0
	covered_rows =0;

	// Mark all rows as not covered
	memset(row_is_covered,0,rows);

	// Set reduced costs to zero
	memset(reduced_costs,0,cols*sizeof(float));

	// Set intial step coefficients
	step_coefficient = stepCoef;
	dual_step_coefficient = stepCoef;	

	// Set minimal step coefficients
	step_coefficient_min = stepCoefMin;
	step_coefficient_iter = stepCoefIter;

	// Create mapping of rows covered by columns

	//cout << "map 1" << endl;
	helpMap = new int[rows];
	
	for(i=0;i<cols;i++) {
		
		
		//cout <<"column" << i << " :" ;
		// Find minimal cost of columns covering this row
		cost = problemView->get_cost(i);
		costs[i]=(float)cost;
		if(cost<minCost) {
			minCost = cost;
		} else if (cost > maxCost) {
			maxCost = cost;
		}

		memset(helpMap,0,rows);
		rowsCovered=0;
		
		// Find all rows covered by this column
		// Store their indexes into help array
		for(j=0;j<rows;j++) {
			if (problemView->get(j,i)==1) {
				helpMap[rowsCovered]=j;
				rowsCovered++;
				//cout << j << ",";
			}
		}
		
		// Copy the help array into mapping
		rows_covered_by_column[i] = new int[rowsCovered];
		memcpy(rows_covered_by_column[i],helpMap,rowsCovered*sizeof(int));
		rows_covered_by_column_count[i] = rowsCovered;
		//cout << endl;
	}
	
	if (helpMap!=NULL) 
		delete[] helpMap;

	//cout << "map 2" << endl;

	// Create mapping of columns covered by rows
	helpMap = new int[cols];
	

	float minColumnCost=0;
	for(i=0;i<rows;i++) {
		memset(helpMap,0,cols);
		colsCovered=0;
		//cout <<"row" << i << "(" << problemView->row_map[i] << "):";

		// Find all columns covering this row
		// Store their indexes into help array
		for(j=0;j<cols;j++) {
			if (problemView->get(i,j)==1) {

				helpMap[colsCovered]=j;
				colsCovered++;
				//cout << j << ",";
				if(minCost==0.)
					minColumnCost = costs[j];
				else
					if(costs[j]<minColumnCost) {
					minColumnCost=costs[j];
					//cout << "found lower cost";
					}
			}
		}

		// Copy the help array into mapping
		columns_covered_by_row[i] = new int[colsCovered];
		memcpy(columns_covered_by_row[i],helpMap,colsCovered*sizeof(int));
		columns_covered_by_row_count[i] = colsCovered;
		min_costs[i]=minColumnCost;
		//cout << endl;

	}
	

	if (helpMap!=NULL) 
		delete[] helpMap;

	
	//cout << "randomizing" << endl;

	// Initialize multipliers
	if(randomizeMultipliers) {
		randomize_multipliers(minCost,maxCost);
	
	} else {

		//cout << "getting first costs: " << endl;
		for(i=0;i<rows;i++) {
			
			//cout << "row " << i << endl;
			// First column covering given row
			j=columns_covered_by_row[i][0];
			// Cost of this column is lowest - columns are already sorted
			t[i]=costs[j];
		}
	}

	//cout << "UB" << endl;

	// ZUB
	best_heur_value = find_heur_ub();
	ZUB=(float)best_heur_value;
	sup_ZUB=(float)ZUB;


	primal_step_change_stop = 0;
	dual_step_change_stop=0;
	
	//cout << "" << endl;

	recommended_column=-1;
	// Algorithm is now initialized
}

int LagrangeSolver::find_heur_ub() {


	int cost= 0;
	int rows = problemView->row_count;
	int i,j;
	int column,row;


	// Mark all rows as not covered
	memset(row_is_covered,0,rows);

	//cout << "greedy" << endl;

	for(i=0;i<rows;i++) {

		//cout << i << " of " << rows << endl;
		
		// Select a row which is not covered
		if(row_is_covered[i]==false) {

			// Get first column covering this row (it's cost is lowest)
			column = columns_covered_by_row[i][0];
			cost += problemView->get_cost(column);

			// Add column to the solution
			x_best[column]=1;

			// Mark all affected rows as covered
			for(j=0;j<rows_covered_by_column_count[column];j++) {
				row = rows_covered_by_column[column][j];
				if(row_is_covered[row] == false) {
					row_is_covered[row] = true;
				}
			}
		}
	}

	//cout << "found simple UB and solution cost " << cost << endl;
	return cost;
}

void LagrangeSolver::randomize_multipliers(int minCost, int maxCost){

	int i;
	int range;


	range = maxCost - minCost;

	//cout << "Randomizing lagrangean multiplicators ...";

	for(i=0;i<problemView->row_count;i++) {
		/* fcost of this column */
		//t[i]= (float) (rand() % (range+1)) + minCost;
		t[i]= 1.;
	}

	for(i=0;i<problemView->column_count;i++) {
		// fcost of this column
		// n[i]= (float) (rand() % (range+1)) + minCost;
		n[i]= 1.;
	}

	//cout << "... done" << endl;


}

void LagrangeSolver::solve_lp() {

	int rows = problemView->row_count;
	int cols = problemView->column_count;

	int i,j;
	int row;
	float totalcost=0.0f;
	float cred=0.0f;

	// Set number of covered rows and selected columns to 0
	covered_rows =0;
	selected_columns = 0;

	// Mark all rows as not covered
	memset(row_is_covered,0,rows);

	for(i=0;i<cols;i++)  {

		cred=get_primal_reduced_cost(i);
		reduced_costs[i]=cred;

		// If reduced cost is negative:
		if (cred<0.0f) {

			// Column is added to the solution
			x[i]=1;

			// Column's cost is added to the cost of solution
			totalcost+=cred;

			/* Column's index is preserved for the purpose of column removal
			 used by the heuristics */
			help_indexes[selected_columns]=i;
			selected_columns+=1;

			// All rows covered by this column are marked as covered
			for(j=0;j<rows_covered_by_column_count[i];j++){
				row = rows_covered_by_column[i][j];
				if(row_is_covered[row]==false) {
					row_is_covered[row]=true;
					covered_rows++;
				}
			}
		}
		else {
			x[i]=0;
		}
	}
	
	// Add all lagrangean multipliers to the lagrangean cost
	for(i=0;i<rows;i++)
		totalcost+=t[i];

	// Update ZLB
	if (totalcost > ZLB) {
		ZLB=totalcost;
		lbound = (int)ceil(ZLB);
	} 
	ZLP=totalcost;

}

void LagrangeSolver::solve_dlp() {

	int rows = problemView->row_count;
	int cols = problemView->column_count;

	int i,j;
	int isum;
	int row;
	float sum;
	float z_LB=0.0f;
	float z_UB=0.0f;
	float cred=0.0f;
	float betterUB=0.0f;

	// Set number of covered rows and selected columns to 0
	covered_rows =0;
	selected_columns = 0;

	// Mark all rows as not covered
	memset(row_is_covered,0,rows);
	memset(y,0,rows*sizeof(float));

	// Compute reduced costs
	for(i=0;i<cols;i++)  {

		cred=get_primal_reduced_cost(i);
		reduced_costs[i]=cred;

		//If reduced cost is negative
		if (cred<0.0f) {

			// the column is added to the solution 
			x[i]=1;

			// Column's cost is added to the cost of lagrangean problem
			z_LB+=cred;
			
			/* Column's index is preserved for the purpose of column removal
			 used by the heuristics */
			help_indexes[selected_columns]=i;
			selected_columns+=1;
			
			// All rows covered by this column are marked as covered
			for(j=0;j<rows_covered_by_column_count[i];j++){
				row = rows_covered_by_column[i][j];
				if(row_is_covered[row]==false) {
					row_is_covered[row]=true;
					covered_rows++;
				}
			}

		}
		else {
			x[i]=0;
		}
	}

	// Compute dual reduced costs
	for(i=0;i<rows;i++)  {

		// Minimal cost of columns covering this row
		min_costs[i]=(float)problemView->get_cost(columns_covered_by_row[i][0]);

		cred=get_dual_reduced_cost(i);
		dual_reduced_costs[i]=cred;

		/// If dual reduced cost is positive, set value in y vector
		y[i]=0;
		if (cred>0.) {
			// set variable of dual vector
			y[i]= min_costs[i];
			z_UB+= y[i] * cred;
		} 

	}


	// Slacks and subgradients
	primal_subgradient=0;
	for(i=0;i<rows;i++) {
		isum=0;
		for(j=0;j<cols;j++) {
			isum += x[j] * problemView->get(i,j);
		}
		primal_slack[i]=1-isum;
		primal_subgradient+= (1-isum)*(1-isum);
	}

	// Dual slacks and subgradients
	dual_subgradient=0.;
	for(j=0;j<cols;j++) {
		sum=0;
		for(i=0;i<rows;i++) {
			sum += y[i] * problemView->get(i,j);
		}
		dual_slack[j]=sum - costs[j];
		dual_subgradient+=(sum-costs[j])*(sum-costs[j]);
	}

	// z_LB
	for(i=0;i<rows;i++)
		z_LB+=t[i];

	ZLP=z_LB;

	// z_UB
	for(i=0;i<cols;i++)
		z_UB+=costs[i]*n[i];

	// Minimum of z_UB, z_hUB
	if(z_UB<ZUB)
		ZUB=z_UB;

	// Maximum of z_LB, z_hLB
	if(z_LB>ZLB) {
		ZLB=z_LB;
		lbound = (int)ceil(ZLB);
	}

	// Primal step
	if(primal_subgradient>0) {
		
		if(ZUB<best_heur_value)
			betterUB = ZUB;
		else
			betterUB = (float)best_heur_value;

		primal_step = step_coefficient * (float)(betterUB - z_LB)/ (float)primal_subgradient;

	} else
		primal_step = 0.;

	// Dual step
	if(dual_subgradient>0) {
		dual_step = dual_step_coefficient  * (z_UB - ZLB)/ dual_subgradient;

	} else
		dual_step = 0.;


	// Update primal multipliers and delta vector
	for(i=0;i<rows;i++) {
		if(primal_slack[i]>0) 
			primal_delta[i]=1;
		else
			primal_delta[i]=0;
		t[i]= t[i] + primal_step * primal_slack[i];
		if(t[i]<0)
			t[i]=0;
	}

	// Update dual multipliers and delta vector
	for(j=0;j<cols;j++) {
		if(dual_slack[j]!=0) {
			dual_delta[j]=1;
		} else
			dual_delta[j]=0;
		n[j]= n[j] + dual_step * dual_slack[j];
		if(n[j]<0)
			n[j]=0;
	}

	if (sup_ZLB!=ZLB) {
		sup_iterations=0;
		sup_ZLB=ZLB;
	} else {
		sup_iterations++;
	}

	if (sup_ZUB!=ZUB) {
		sup_dual_iterations=0;
		sup_ZUB=ZUB;
	} else {
		sup_dual_iterations++;
	}

}

void LagrangeSolver::solve_dlp_next() {

	int rows = problemView->row_count;
	int cols = problemView->column_count;

	int i,j,k;
	float z_LB=0.0f;
	float z_UB=0.0f;
	float cred=0.0f;
	int column,row;
	float betterUB = 0.0f;

	// Set number of covered rows and selected columns to 0
	covered_rows =0;
	selected_columns = 0;

	// Mark all rows as not covered
	memset(row_is_covered,0,rows);
	memset(column_is_recounted,0,cols);
	memset(row_is_recounted,0,rows);


		// Change primal reduced costs
		for(i=0;i<rows;i++) {

			if(primal_delta[i]!=0) {

				for(j=0;j<columns_covered_by_row_count[i];j++) {

					column = columns_covered_by_row[i][j];

					if(column_is_recounted[column]==0) {

						cred=get_primal_reduced_cost(column);

						// Has change in reduced cost occured ?
						if( (cred<0. && reduced_costs[column] >= 0.) ||
							(cred>=0. && reduced_costs[column] < 0.)) {
								

								// If yes, update slacks for rows covered by this column
								for(k=0;k<rows_covered_by_column_count[column];k++) {
									row = rows_covered_by_column[column][k];

									if(cred>=0.) {
										primal_slack[row] +=1;
									} else {
										primal_slack[row] -=1;
									}
								}
						}
						reduced_costs[column]=cred;
						column_is_recounted[column]=1;
					}

				}
			}
		}

		for(i=0;i<cols;i++)  {

			cred=reduced_costs[i];
			// If reduced cost is negative
			if (cred<0.0f) {
				
				// Column is added to the solution
				x[i]=1;

				// Column's cost is added to the cost of lagrangean problem
				z_LB+=cred;
				
				// Column's index is preserved for the purpose of column removal used by the heuristics */
				help_indexes[selected_columns]=i;
				selected_columns+=1;
			
				// All rows covered by this column are marked as covered
				for(j=0;j<rows_covered_by_column_count[i];j++){
					row = rows_covered_by_column[i][j];
					if(row_is_covered[row]==false) {
						row_is_covered[row]=true;
						covered_rows++;
					}
				}
			}
			else {
				x[i]=0;
			}
		}

		// slacks and subgradients
		primal_subgradient=0;
		for(i=0;i<rows;i++) {
			primal_subgradient+= primal_slack[i]*primal_slack[i];
		}

		// z_LB 
		for(i=0;i<rows;i++)
			z_LB+=t[i];

		ZLP=z_LB;

		// maximum of z_LB, z_hLB
		if(z_LB>ZLB){
			ZLB=z_LB;
			lbound = (int)ceil(ZLB);
		}


		// Change dual reduced costs
		for(i=0;i<cols;i++) {

			if(dual_delta[i]!=0) {

				for(j=0;j<rows_covered_by_column_count[i];j++) {

					row = rows_covered_by_column[i][j];

					if(row_is_recounted[row]==0) {

						cred=get_dual_reduced_cost(row);

						if( (cred<=0. && dual_reduced_costs[row] > 0.) ||
							(cred>0. && dual_reduced_costs[row] <= 0.)) {

								for(k=0;k<columns_covered_by_row_count[row];k++) {
									column = columns_covered_by_row[row][k];

									if(cred>=0.) {
										dual_slack[column] +=min_costs[row];
									} else {
										dual_slack[column] -=min_costs[row];
									}
								}
						}
						dual_reduced_costs[row]=cred;
						row_is_recounted[row]=1;
					}
				}
			}
		}

		for(i=0;i<rows;i++)  {

			//If dual reduced cost is positive
			if ( dual_reduced_costs[i]>0.) {

				// Set variable of dual vector
				y[i]=min_costs[i];
				z_UB+= y[i] * dual_reduced_costs[i];
			}
			else {
				y[i]=0.;
			}
		}

		// Dual slacks and subgradients
		dual_subgradient=0.;
		for(j=0;j<cols;j++) {
			dual_subgradient+=(dual_slack[j])*(dual_slack[j]);
		}

		// z_UB
		for(i=0;i<cols;i++)
			z_UB+=costs[i]*n[i];

		// minimum of z_UB, z_hUB
		if(z_UB<ZUB)
			ZUB=z_UB;

		// Primal step
		if(primal_subgradient>0) {

			if(ZUB<best_heur_value)
				betterUB = ZUB;
			else
				betterUB = (float)best_heur_value;

			primal_step = step_coefficient * (float)(betterUB - z_LB)/ (float)primal_subgradient;

		} else
			primal_step = 0.;

		// Update multipliers
		for(i=0;i<rows;i++) {
			if(primal_slack[i]!=0) {
				primal_delta[i]=1;
			} else
				primal_delta[i]=0;
			t[i]= t[i] + primal_step * primal_slack[i];
			if(t[i]<0)
				t[i]=0;
		}


	if(primal_step_change_stop==0) {

		// If needed, update the step koeficient
		if (sup_ZLB!=ZLB) {
			sup_iterations=0;
			sup_ZLB=ZLB;
		} else {
			sup_iterations++;
			if (sup_iterations>=step_coefficient_iter) {
				sup_iterations=0;
				step_coefficient/=2;
			}
			if(step_coefficient <= step_coefficient_min) {
				primal_step_change_stop = 1;
			}
		}

	}

	

		// Dual step
		if(dual_subgradient>0) {
			dual_step = dual_step_coefficient  * (z_UB - ZLB)/ dual_subgradient;

		} else
			dual_step = 0.;

		// Update dual multipliers
		for(j=0;j<cols;j++) {
			if(dual_slack[j]!=0) {
				dual_delta[j]=1;
			} else
				dual_delta[j]=0;
			n[j]= n[j] + dual_step * dual_slack[j];
			if(n[j]<0)
				n[j]=0;
		}
		
	if(dual_step_change_stop==0) {

		// If needed, update the dual step koeficient
		if (sup_ZUB!=ZUB) {
			sup_dual_iterations=0;
			sup_ZUB=ZUB;
		} else {
			sup_dual_iterations++;
			if (sup_dual_iterations>=step_coefficient_iter) {
				sup_dual_iterations=0;
				dual_step_coefficient/=2;
			}
			if(dual_step_coefficient <= step_coefficient_min) {
				dual_step_change_stop = 1;
			}
		}

	}

}


void LagrangeSolver::ha(int removeType) {

	int cols = problemView->column_count;
	int rows = problemView->row_count;

	int i,j,k;
	int row,column,selectedColumn;
	float costs;

	float minReducedCost;
	int candidateColumn;

	bool selected_row_is_covered;
	bool can_be_removed;

	// Copy the vector of solution of LP
	memcpy(x_heur,x,cols);

	if(covered_rows < problemView->row_count) {

		for(i=0;i<rows;i++) {

			// Select a row which is not covered
			if(row_is_covered[i]==false) {

				minReducedCost = (float)POS_INFINITY;
				candidateColumn= -1;

				// Search all relevant columns and get their reduced costs
				for(j=0;j<columns_covered_by_row_count[i];j++) {
					column = columns_covered_by_row[i][j];

					if(x_heur[column]==0) {
						if(reduced_costs[column]<minReducedCost) {
							minReducedCost = reduced_costs[column];
							candidateColumn = column ;
						}
					}
				}

				// Get column with minimal reduced cost and add it to solution
				if(candidateColumn>=0) {
					x_heur[candidateColumn] = 1;
					help_indexes[selected_columns]=candidateColumn;
					selected_columns+=1;

					// Mark all rows covered by this column as covered
					for(j=0;j<rows_covered_by_column_count[candidateColumn];j++) {
						row = rows_covered_by_column[candidateColumn][j];

						if(row_is_covered[row] == false) {
							row_is_covered[row] =true;
							covered_rows++;
						}
					}
				} else {
					cout << endl << "ALERT - candidate column not found" << endl;
				}

			}	
		}

	}

	for(i=0;i<rows;i++) {
		if(row_is_covered[i]==false)
			cout << endl << "ALERT - not covered row found" << endl;
	}

	if(removeType==0) {

		// Test, whether given column can be removed
		for(i=cols-1;i>=0;i--) {
			if(x_heur[i]==1){

				// Temporarily remove column from solution
				x_heur[i]=0;
				can_be_removed=1;

				// For every row covered by given column
				for(j=0;j<rows_covered_by_column_count[i];j++) {
					row = rows_covered_by_column[i][j];
					selected_row_is_covered=false;

					// Test, whether given row is covered by any column
					for(k=0;k<columns_covered_by_row_count[row];k++){
						column = columns_covered_by_row[row][k];
						if(x_heur[column]==1) {
							selected_row_is_covered = true;
							break;
						}
					}

					// If not, column cannot be removed from solution
					if(selected_row_is_covered==false) {
						can_be_removed =0;
						break;
					}

				}

				// Add column back to solution if it cannot be removed
				if(can_be_removed==0) {
					x_heur[i]=1;
				}
			}
		}


	} else {
	
		// Sort all columns that are part of solution according to their reduced costs
		for(i=0;i<selected_columns;i++) {
			help_costs[i]= reduced_costs[help_indexes[i]];
			help_number_of_ones[i]= rows_covered_by_column_count[help_indexes[i]];
		}
		quick_sort_three_arrays(help_indexes,help_costs,help_number_of_ones,0,selected_columns-1);

		// Test, whether given column can be removed
		for(i=0;i<selected_columns;i++) {

			// Temporarily remove column from solution
			selectedColumn = help_indexes[i];
			x_heur[selectedColumn]=0;
			can_be_removed=1;

			// For every row covered by given column
			for(j=0;j<rows_covered_by_column_count[selectedColumn];j++) {
				row = rows_covered_by_column[selectedColumn][j];
				selected_row_is_covered=false;

				// Test, whether given row is covered by any column
				for(k=0;k<columns_covered_by_row_count[row];k++){
					column = columns_covered_by_row[row][k];
					if(x_heur[column]==1) {
						selected_row_is_covered = true;
						break;
					}
				}

				// If not, column cannot be removed from solution
				if(selected_row_is_covered==false) {
					can_be_removed =false;
					break;
				}

			}

			// Add column back to solution if it cannot be removed
			if(can_be_removed==false) {
				x_heur[selectedColumn]=1;
			}
		}
	}

	// Update best heuristic solution found
	costs=0.0f;
	for(i=0;i<cols;i++) {
		if (x_heur[i]==1) {
			costs+=problemView->get_cost(i);
		}
	}

	heur_value = costs;
	update_best_solution(costs);
}


void LagrangeSolver::ha_modified(int removeType) {

	int cols = problemView->column_count;
	int rows = problemView->row_count;

	int i,j,k;
	int row,column,selectedColumn;
	float costs;

	float minReducedCost;
	int candidateColumn;

	bool selected_row_is_covered;
	bool can_be_removed;

	// Copy the vector of solution of LP
	memcpy(x_heur,x,cols);
	memcpy(modified_reduced_costs,reduced_costs,cols*sizeof(float));

	
	for(i=0;i<rows;i++) {

		// Select a row which is not covered 
		if(row_is_covered[i]==false) {

			minReducedCost = (float)POS_INFINITY;
			candidateColumn= -1;

			// Search all relevant columns and get their modified reduced costs
			for(j=0;j<columns_covered_by_row_count[i];j++) {
				column = columns_covered_by_row[i][j];

				if(x_heur[column]==0) {
					if(modified_reduced_costs[column]<minReducedCost) {
						minReducedCost = modified_reduced_costs[column];
						candidateColumn = column ;
					}
				}
			}

			// Get column with minimal modified reduced cost and add it to solution
			if(candidateColumn>=0) {
				x_heur[candidateColumn] = 1;
				help_indexes[selected_columns]=candidateColumn;
				selected_columns+=1;

				// Mark rows covered by this column as covered
				for(j=0;j<rows_covered_by_column_count[candidateColumn];j++) {
					row = rows_covered_by_column[candidateColumn][j];

					if(row_is_covered[row] == false) {
						row_is_covered[row] = true;
						covered_rows++;

						// Modify reduced cost for every column covering newly covered row
						for(k=0;k<columns_covered_by_row_count[row];k++) {
							column = columns_covered_by_row[row][k];
							modified_reduced_costs[column] += t[row];
						}

					}
				}
			} else {
				cout << endl << "ALERT - candidate column not found" << endl;
			}

		}	
	}

	for(i=0;i<rows;i++) {
		if(row_is_covered[i]==false)
			cout << endl << "ALERT - not covered row found" << endl;
	}


	if(removeType==0) {

		// Test, whether given column can be removed
		for(i=cols-1;i>=0;i--) {
			if(x_heur[i]==1){

				// Temporarily remove column from solution
				x_heur[i]=0;
				can_be_removed=true;

				// For every row covered by given column
				for(j=0;j<rows_covered_by_column_count[i];j++) {
					row = rows_covered_by_column[i][j];
					selected_row_is_covered=false;

					// Test, whether given row is covered by any column
					for(k=0;k<columns_covered_by_row_count[row];k++){
						column = columns_covered_by_row[row][k];
						if(x_heur[column]==1) {
							selected_row_is_covered = true;
							break;
						}
					}

					// If not, column cannot be removed from solution
					if(selected_row_is_covered==false) {
						can_be_removed =false;
						break;
					}

				}

				// Add column back to solution if it cannot be removed
				if(can_be_removed==false) {
					x_heur[i]=1;
				}
			}
		}


	} else {

		if(removeType==1 ) {
			
			// Sort all columns that are part of solution according to their reduced costs
			for(i=0;i<selected_columns;i++) {
				help_costs[i]= reduced_costs[help_indexes[i]];
				help_number_of_ones[i]= rows_covered_by_column_count[help_indexes[i]];
			}
			quick_sort_three_arrays(help_indexes,help_costs,help_number_of_ones,0,selected_columns-1);

		} else {
			
			// Sort all columns that are part of solution according to their modified reduced costs
			for(i=0;i<selected_columns;i++) {
				help_costs[i]= modified_reduced_costs[help_indexes[i]];
				help_number_of_ones[i]= rows_covered_by_column_count[help_indexes[i]];
			}
			quick_sort_three_arrays(help_indexes,help_costs,help_number_of_ones,0,selected_columns-1);

		}

		// Test, whether given column can be removed
		for(i=0;i<selected_columns;i++) {

			// Temporarily remove column from solution
			selectedColumn = help_indexes[i];
			x_heur[selectedColumn]=0;
			can_be_removed=1;

			// For every row covered by given column
			for(j=0;j<rows_covered_by_column_count[selectedColumn];j++) {
				row = rows_covered_by_column[selectedColumn][j];
				selected_row_is_covered=false;

				// Test, whether given row is covered by any column
				for(k=0;k<columns_covered_by_row_count[row];k++){
					column = columns_covered_by_row[row][k];
					if(x_heur[column]==1) {
						selected_row_is_covered = true;
						break;
					}
				}

				// If not, column cannot be removed from solution
				if(selected_row_is_covered==false) {
					can_be_removed =0;
					break;
				}

			}

			// Add column back to solution if it cannot be removed 
			if(can_be_removed==0) {
				x_heur[selectedColumn]=1;
			}
		}
	}

	// Update best heuristic solution found
	costs=0.0f;
	for(i=0;i<cols;i++) {
		if (x_heur[i]==1) {
			costs+=problemView->get_cost(i);
		}
	}

	heur_value = costs;
	update_best_solution(costs);

}

int LagrangeSolver::compute_subgradient() {

	int i,j;
	int column;

	int row_count=problemView->row_count;
	//int column_count=problemView->column_count;

	int sum=0;

	// Slacks
	for(i=0;i<row_count;i++) {
		sum=0;
		for(j=0;j<columns_covered_by_row_count[i];j++) {
			column = columns_covered_by_row[i][j];
			sum += x[column];
		}
		if(t[i]==0 && (1-sum)<0 ){
			primal_slack[i]=0;
		} else {
			primal_slack[i]=1-sum;
		}
	}

	// Subgradient
	primal_subgradient=0;
	for(i=0;i<row_count;i++) {
		primal_subgradient+=primal_slack[i]*primal_slack[i];
	}

	// Should we stop ?
	if(primal_subgradient==0)
		return 0;

	// If needed, update the step koeficient
	if (sup_ZLB!=ZLB) {
		sup_iterations=0;
		sup_ZLB=ZLB;
	} else {
		sup_iterations++;
		if (sup_iterations>=step_coefficient_iter) {
			sup_iterations=0;
			step_coefficient/=2;
		}
		if(step_coefficient <= step_coefficient_min) {
			return 0;
		}
	}

	// Calculate the step
	primal_step=(float) (step_coefficient*(best_heur_value - ZLP)/primal_subgradient);

	// Update lagrange multipliers
	for(i=0;i<row_count;i++) {
		t[i]= t[i] + primal_step * primal_slack[i];
		if(t[i]<0)
			t[i]=0.;
	}

	return 1;

}


float LagrangeSolver::get_primal_reduced_cost(int columnIndex) {

	int i=0;
	int row;
	float sum=0.0f;
	

	// Cost of column
	float cred=costs[columnIndex];

	// Plus all lagrangean mutlipliers of rows covered by column
	for(i=0;i<rows_covered_by_column_count[columnIndex];i++) {
		row = rows_covered_by_column[columnIndex][i];
		sum+= t[row];
	}

	cred-=sum;

	return cred;

}

float LagrangeSolver::get_dual_reduced_cost(int rowIndex) {

	int i=0;
	int column;
	float sum=0.0f;

	float cred=1.;

	// sum of dual multipliers of all columns covering given row
	for(i=0;i<columns_covered_by_row_count[rowIndex];i++) {
		column = columns_covered_by_row[rowIndex][i];
		sum+= n[column];
	}

	cred-=sum;

	return cred;

}

void LagrangeSolver::update_best_solution(float candCost) {

	// If proposed solution cost is lower than original
	if (candCost<best_heur_value) {
		best_heur_value = (int) candCost;
		
		// Update best solution found
		memcpy(x_best,x_heur,problemView->column_count);
	}
}

int compare_three_arrays(float costs[],int number_of_ones[],int c1, int c2) {


	// Compare costs of columns
	if (costs[c1] < costs[c2]) 
		return -1; 
	else 
		if (costs[c1] > costs[c2] ) 
			return +1;

	// If costs are equal, compare number of rows covered
	if (number_of_ones[c1] > number_of_ones[c2]) 
		return -1; 
	else 
		if (number_of_ones[c1] < number_of_ones[c2]) 
			return +1;

	return 0;
}

void swap_in_arrays(int columns[], float costs[],int number_of_ones[],int left, int right){

	int index_temp;
	int	number_temp;
	float cost_temp;

	index_temp = columns[left];
	number_temp = number_of_ones[left];
	cost_temp = costs[left];

	columns[left] = columns [right];
	number_of_ones[left] = number_of_ones[right];
	costs[left] = costs[right];

	columns [right] = index_temp;
	number_of_ones[right] = number_temp;
	costs[right] = cost_temp;
}


void quick_sort_three_arrays(int columns[], float costs[],int number_of_ones[],int l, int r) {

	int i=l;
	int j=r;

	/*Pivot index*/
	int x=(l+r)/2;

	do {
		while( compare_three_arrays(costs,number_of_ones,i,x)<0 ) i++;
		while( compare_three_arrays(costs,number_of_ones,x,j)<0 ) j--;

		if (i<=j) {
			swap_in_arrays(columns,costs,number_of_ones,i,j);
			i++;
			j--;
		}

	}
	while(i<=j);

	if (l<j) quick_sort_three_arrays(columns,costs,number_of_ones,l,j);
	if (i<r) quick_sort_three_arrays(columns,costs,number_of_ones,i,r);
}

void swap_in_links(CLink* links[],int left, int right){

	CLink* link_temp;

	link_temp = links[left];
	links[left] = links[right];
	links[right] = link_temp;
}

int compare_column_links(CLink* links[],int c1, int c2){

	if (links[c1]->Cost < links[c2]->Cost) 
		return -1; 
	else 
		if (links[c1]->Cost > links[c2]->Cost) 
			return +1;


	if (links[c1]->NumberOf1>links[c2]->NumberOf1) 
		return -1; 
	else 
		if (links[c1]->NumberOf1<links[c2]->NumberOf1) 
			return +1;

	return 0;
}

int compare_row_links(CLink* links[],int r1, int r2){

	if (links[r1]->NumberOf1 < links[r2]->NumberOf1) 
		return -1; 
	else 
		if (links[r1]->NumberOf1 > links[r2]->NumberOf1) 
			return +1;

	return 0;
}


void quick_sort_column_links(CLink* links[], int l, int r) {

	int i=l;
	int j=r;

	/*Pivot index*/
	int x=(l+r)/2;

	do {
		while(compare_column_links(links,i,x)<0) i++;
		while(compare_column_links(links,x,j)<0) j--;

		if (i<=j) {
			swap_in_links(links,i,j);
			i++;
			j--;
		}

	}
	while(i<=j);

	if (l<j) quick_sort_column_links(links,l,j);
	if (i<r) quick_sort_column_links(links,i,r);
}

void quick_sort_row_links(CLink* links[], int l, int r) {

	int i=l;
	int j=r;

	/*Pivot index*/
	int x=(l+r)/2;

	do {
		while(compare_row_links(links,i,x)<0) i++;
		while(compare_row_links(links,x,j)<0) j--;

		if (i<=j) {
			swap_in_links(links,i,j);
			i++;
			j--;
		}

	}
	while(i<=j);

	if (l<j) quick_sort_row_links(links,l,j);
	if (i<r) quick_sort_row_links(links,i,r);
}

void LagBaBCover::release (struct CoveringMatrix * CM) {
	
	struct CLink *Aux;
	struct CLink *Aux2;
	if ( CM->ColLink !=0 ) {	
		
		// Delete all rows
		Aux=CM->ColLink->Next;
		delete (CM->ColLink);
		while ( Aux!=0 )
		{
			Aux2=Aux->Next;
			delete (Aux);
			Aux=Aux2;
		}
	}
	if ( CM->RowLink !=0 ) {

		// Delete all columns
		Aux=CM->RowLink->Next;
		delete (CM->RowLink);
		while ( Aux!=0 )
		{
			Aux2=Aux->Next;
			delete (Aux);
			Aux=Aux2;
		}
	}
	
	// Delete header of structure
	delete (CM);
}

void LagBaBCover::release ( struct header_of_solution_UCP *SolUCP ) { 
	
	CLink *Aux;
	CLink *Aux2;

	if ( SolUCP->SetC !=0 ) {	

		// Delete all column links
		Aux=SolUCP->SetC->Next;
		delete ( SolUCP->SetC );
		while ( Aux!=0 )
		{
			Aux2=Aux->Next;
			delete (Aux);
			Aux=Aux2;
		}
	}

	// Delete header
	delete ( SolUCP );
}

struct header_of_solution_UCP *LagBaBCover::create_New_Empty_Solution ( void ){

	struct header_of_solution_UCP *Solution;

	Solution=new header_of_solution_UCP;
	
	if ( Solution==NULL) 
		throw RuntimeException( "Aura: Not enough memory" );
	
	// Set all antries to zero
	Solution->Cost=0;
	Solution->SetC=0;

	return Solution;
}

void LagBaBCover::append_to_solution (struct header_of_solution_UCP *Source,struct header_of_solution_UCP *Target ) {
	
	struct CLink *Aux;

	// If there are no cols in solution of UCP
	if ( Target->SetC==0 ) 
		Target->SetC=Source->SetC;

	// else connect to the end of link
	else {
		
		Aux=Target->SetC;
		while ( Aux->Next!=0 )
			Aux=Aux->Next;
		Aux->Next=Source->SetC;
	}

	Target->Cost+=Source->Cost;

	// Delete header of partial solution,it is no inter needed
	delete (Source);
}


void LagBaBCover::add_to_solution ( struct header_of_solution_UCP *SolutionUCP,struct CLink *Column ){
	
	struct CLink *Solution;
	struct CLink *Aux;

	Solution=new CLink;
	if ( Solution==NULL) throw RuntimeException( "Aura: Not enough memory" );
	Solution->ID=Column->ID;
	Solution->Cost=Column->Cost;
	//Solution->Offset=Column->Offset;
	Solution->NumberOf1=Column->NumberOf1;
	Solution->Next=0;
	
	// If there are no entries in solution
	if ( SolutionUCP->SetC==0 ) {
		SolutionUCP->SetC=Solution;
		SolutionUCP->Cost=Column->Cost;

	// else append it to the end od links
	} else {
		
		SolutionUCP->Cost+=Column->Cost;
		Aux=SolutionUCP->SetC;
		while ( Aux->Next!=0 )
			Aux=Aux->Next;
		Aux->Next=Solution;
	}

}

void LagBaBCover::add_to_solution_masked( struct header_of_solution_UCP *TargetSolution,struct CoveringMatrix* Matrix,unsigned char* mask, int len) {
	
	struct CLink *Solution = NULL;
	struct CLink *ColPointer = NULL;
	struct CLink *TargetPointer = NULL;
	int columnCount = len;
	int i;

	// set pointers
	ColPointer = Matrix->ColLink;

	if ( TargetSolution->SetC!=0 ){ 
		TargetPointer=TargetSolution->SetC;
		while ( TargetPointer->Next!=0 )
			TargetPointer=TargetPointer->Next;
	}

	ColPointer = Matrix->ColLink;
	while(ColPointer->Next != 0) {
		ColPointer=ColPointer->Next;
	}

	// Iterate over mask
	for(i=columnCount-1;i>=0;i--) {
		
		// if given column is part of solution
		if(mask[i]==1) {

			// add column link to solution
			Solution= new CLink();

			// find column link in matrix
			ColPointer = Matrix->ColLink;
			while(ColPointer->ID != i) {
				ColPointer = ColPointer->Next;
				if(ColPointer==0)
					throw BoomException( "LagBaBCover::add_to_solution_masked: could not find column link" );
			}

			// copy content of column link
			Solution->ID=ColPointer->ID;
			Solution->Cost=ColPointer->Cost;
			Solution->NumberOf1=ColPointer->NumberOf1;
			Solution->Next=0;

			// add new column poiinter to the end of solution
			if ( TargetSolution->SetC==0 ){ 
				TargetSolution->SetC=Solution;
				TargetSolution->Cost=Solution->Cost;
				TargetPointer=Solution;
			} else {
				TargetSolution->Cost+=Solution->Cost;
				TargetPointer->Next=Solution;
				TargetPointer=Solution;
			}
		}
	}
}

void LagBaBCover::remove_col_from_matrix ( struct CoveringMatrix *CM, struct CLink *colToRemove ) {
	
	struct CLink *Entry;

	Entry=CM->ColLink;

	if ( Entry !=0 ) {

		// found the right entry
		while ( Entry->ID!=colToRemove->ID )
			Entry=Entry->Next;

		// entry is in the middle
		if ( Entry->Next != 0&&Entry->Previous!=0 ) {
			Entry->Next->Previous=Entry->Previous;
			Entry->Previous->Next=Entry->Next;
		} else
			// entry is at the beginning
			if ( Entry->Next!=0&&Entry->Previous==0	) {
				Entry->Next->Previous=0;
				CM->ColLink=Entry->Next;
			} else
				// entry is at the end
				if ( Entry->Next == 0&&Entry->Previous!=0 ) {
					Entry->Previous->Next=0;
				} else
					// entry is the only column
					if ( Entry->Next == 0&&Entry->Previous==0 ) {
						
						// if it is the only column in the matrix,all rows must be deleted
						// to set the matrix empty
						CM->ColLink=0;
					}

					// decrease dimension of matrix
					CM->countOfCols--;
					// delete column
					delete (Entry);
	}
}

void LagBaBCover::remove_col_from_matrix ( struct CoveringMatrix *CM, int ID ) {

	struct CLink *Entry;

	Entry=CM->ColLink;

	if ( Entry !=0 ) {

		// found the right entry
		while ( Entry->ID!=ID )
			Entry=Entry->Next;
		
		// entry is in the middle
		if ( Entry->Next != 0&&Entry->Previous!=0 ) {
			Entry->Next->Previous=Entry->Previous;
			Entry->Previous->Next=Entry->Next;
		} else
			// entry is at the beginning
			if ( Entry->Next!=0&&Entry->Previous==0	) {
				Entry->Next->Previous=0;
				CM->ColLink=Entry->Next;
			} else
				// entry is at the end
				if ( Entry->Next == 0&&Entry->Previous!=0 ) {
					Entry->Previous->Next=0;
				} else
					// entry is the only column
					if ( Entry->Next == 0&&Entry->Previous==0 ) {
						
						// if it is the only column in the matrix,all rows must be deleted
						// to set the matrix empty
						CM->ColLink=0;
					}

					// decrease dimension of matrix
					CM->countOfCols--;

					// delete column
					delete (Entry);
	}
}

void LagBaBCover::remove_row_from_matrix ( struct CoveringMatrix *CM, struct CLink *rowToRemove ) {
	
	struct CLink *Entry;
	Entry=CM->RowLink;

	if ( Entry !=0 ) {


		// found the right entry
		while ( Entry->ID!=rowToRemove->ID )
			Entry=Entry->Next;
		
		// entry is in the middle
		if ( Entry->Next != 0 &&Entry->Previous!=0 ) {
			Entry->Next->Previous=Entry->Previous;
			Entry->Previous->Next=Entry->Next;
		} else
			// entry is at the beginning
			if ( Entry->Next!=0 &&Entry->Previous==0){
				Entry->Next->Previous=0;
				CM->RowLink=Entry->Next;
			} else
				// entry is at the end
				if(Entry->Next == 0 &&Entry->Previous!=0) {
					Entry->Previous->Next=0;
				} else
					// entry is the only row
					if ( Entry->Next == 0&&Entry->Previous==0 ) {
						CM->RowLink=0;
					}

					//delete row
					delete (Entry);

					//decrease dimension of matrix
					CM->countOfRows--;
	}
}

void LagBaBCover::remove_col_with_rows ( struct CLink * Column, struct CoveringMatrix *CM ) {
	
	struct CLink *RowLink;
	struct CLink *Aux;

	RowLink = CM->RowLink;
	while ( RowLink!=0) {
		if ( CM->Array[Column->ID][RowLink->ID]==1 ) {
			// if column intersect row, that row will be deleted
			Aux=RowLink->Next;
			remove_row_from_matrix(CM,RowLink );
		} else 
			Aux=RowLink->Next;
		RowLink=Aux;
	}

	if ( CM->ColLink!=0 ) 
		remove_col_from_matrix (CM,Column);
}

struct CoveringMatrix* LagBaBCover::duplicate ( struct CoveringMatrix *CM) {
	
	struct CoveringMatrix *Copy;
	struct CLink *Rows;
	struct CLink *Columns;
	struct CLink *CMRowLink;
	struct CLink *CMColLink;

	CMRowLink=CM->RowLink;
	CMColLink=CM->ColLink;
	Copy=new CoveringMatrix;
	if ( Copy==NULL) throw RuntimeException( "Aura: Not enough memory" );
	Copy->Cost=CM->Cost;
	Copy->countOfOnes=CM->countOfOnes;
	Copy->countOfCols=CM->countOfCols;
	Copy->countOfRows=CM->countOfRows;
	//Copy->MSIR=CM->MSIR;
	//Copy->lbound=CM->lbound;
	Copy->Array=CM->Array;
	// Set all antries of created link to values of link of original matrix
	Copy->Dirty=CM->Dirty;
	Copy->ColLink=NULL;
	Copy->RowLink=NULL;
	
	// if there are another entries,add it to duplicated matrix
	if ( CM->RowLink!=NULL ) {
		Copy->RowLink=new CLink;
		if ( Copy->RowLink==NULL) throw RuntimeException( "Aura: Not enough memory" );
		Copy->RowLink->Previous=NULL;
		Copy->RowLink->Next=NULL;
		//Copy->RowLink->Offset=CM->RowLink->Offset;
		Copy->RowLink->ID=CM->RowLink->ID;
		Copy->RowLink->NumberOf1=CM->RowLink->NumberOf1;
		Copy->RowLink->Cost=CM->RowLink->Cost;
		Rows=Copy->RowLink;
		CMRowLink=CM->RowLink->Next;


		while ( CMRowLink!=NULL )
		{
			Rows->Next=new CLink;
			if ( Rows->Next==NULL ) throw RuntimeException( "Aura: Not enough memory" );

			Rows->Next->Previous=Rows;
			Rows=Rows->Next;
			//Rows->Offset=CMRowLink->Offset;
			Rows->ID=CMRowLink->ID;
			Rows->NumberOf1=CMRowLink->NumberOf1;
			Rows->Cost=CMRowLink->Cost;
			// Set all antries of created link to values of link of original matrix
			Rows->Next=NULL;
			CMRowLink=CMRowLink->Next;
		}
	}

	if ( CM->ColLink!=NULL )
	{

		Copy->ColLink = new CLink;
		if ( Copy->ColLink==NULL) throw RuntimeException( "Aura: Not enough memory" );

		Copy->ColLink->Previous=NULL;
		Copy->ColLink->Next=NULL;
		//Copy->ColLink->Offset=CM->ColLink->Offset;
		Copy->ColLink->NumberOf1=CM->ColLink->NumberOf1;
		Copy->ColLink->Cost=CM->ColLink->Cost;
		Copy->ColLink->ID=CM->ColLink->ID;
		Columns=Copy->ColLink;
		CMColLink=CM->ColLink->Next;
		// Set all antries of created link to values of link of original matrix

		while ( CMColLink!=NULL )
		{
			Columns->Next=new CLink;
			if ( Columns->Next==NULL ) throw RuntimeException( "Aura: Not enough memory" );

			Columns->Next->Previous=Columns;
			Columns=Columns->Next;

			//Columns->Offset=CMColLink->Offset;
			Columns->ID=CMColLink->ID;
			Columns->NumberOf1=CMColLink->NumberOf1;
			Columns->Cost=CMColLink->Cost;

			Columns->Next=NULL;
			CMColLink=CMColLink->Next;
			// Set all antries of created link to values of link of original matrix
		}
	}
	return Copy;
}

void LagBaBCover::create_matrix ( struct CoveringMatrix *Cm,int *Costs ) {
	
	int I;
	struct CLink *Columns;
	struct CLink *Rows;
	
	// Create row links
	if (Cm->countOfCols!=0||Cm->countOfRows!=0)
	{
		Cm->RowLink = new CLink;
		if ( Cm->RowLink==NULL ) throw RuntimeException( "Aura: Not enough memory" );
		Cm->RowLink->Previous=0;
		//Cm->RowLink->Offset=Cm->countOfCols;
		Cm->RowLink->ID=0;
		Cm->RowLink->NumberOf1=0;
		Rows=Cm->RowLink;
		if (Cm->countOfRows>1)
			for ( I=2;I<=Cm->countOfRows;I++)
			{
				Rows->Next=new CLink;
				if ( Rows->Next==NULL) throw RuntimeException( "Aura: Not enough memory" );
				Rows->Next->Previous=Rows;
				Rows=Rows->Next;
				//Rows->Offset=I*Cm->countOfCols*sizeof(char);
				Rows->ID=I-1;
				Rows->NumberOf1=0;
				Rows->Next=0;
			}
		else  Rows->Next=0;
		// Create column links
		Cm->ColLink = new CLink;
		if ( Cm->ColLink==NULL) throw RuntimeException( "Aura: Not enough memory" );
		Cm->ColLink->Previous=0;
		//Cm->ColLink->Offset=0;
		Cm->ColLink->NumberOf1=0;
		Cm->ColLink->ID=Cm->countOfCols-1;
		
		// Cost of the column is kept at the beginning of array
		Cm->ColLink->Cost=Costs[Cm->countOfCols-1];

		//upper bound is computed here
		uBound+=Cm->ColLink->Cost;
		Columns=Cm->ColLink;
		if (Cm->countOfCols>1)
			for ( I=Cm->countOfCols-2;I>=0;I--)
			{
				Columns->Next=new CLink;
				if ( Columns->Next==NULL) throw RuntimeException( "Aura: Not enough memory" );
				Columns->Next->Previous=Columns;
				Columns=Columns->Next;
				//Columns->Offset=I*sizeof(char);
				Columns->ID=I;
				Columns->NumberOf1=0;
				Columns->Cost=Costs[I];
				uBound+=Columns->Cost;
				Columns->Next=0;
			}
		else  Columns->Next=0;
		uBound++;
		MaxCost=uBound;
	}
	else {Cm->RowLink=0;Cm->ColLink=0;}
}

void LagBaBCover::sort_matrix ( struct CoveringMatrix *Cm ) {

	CLink** links;
	CLink* link;
	int i;


	find_number_of_ones(Cm);

	// create array with column links
	links = new CLink*[Cm->countOfCols];
	link = Cm->ColLink;
	for(i=0;i<Cm->countOfCols;i++) {
		links[i] = link;
		link=link->Next;

	}
	
	// sort this array according to column costs and number of rows covered
	quick_sort_column_links(links,0,i-1);

	// rebuild column links in covering matrix
	Cm->ColLink = links[0];
	
	for(i=0;i<(Cm->countOfCols-1);i++) {
		links[i]->Next=links[i+1];
	}
	for(i=1;i<Cm->countOfCols;i++) {
		links[i]->Previous=links[i-1];
	}
	links[Cm->countOfCols-1]->Next = 0;
	links[0]->Previous = 0;

	delete[] links;
	
	// create array with row links
	links = new CLink*[Cm->countOfRows];
	link = Cm->RowLink;
	for(i=0;i<Cm->countOfRows;i++) {
		links[i] = link;
		link=link->Next;

	}
	
	// sort this array
	quick_sort_row_links(links,0,i-1);


	// rebuild row links in covering matrix
	Cm->RowLink = links[0];

	for(i=0;i<(Cm->countOfRows-1);i++) {
		links[i]->Next=links[i+1];
	}
	for(i=1;i<Cm->countOfRows;i++) {
		links[i]->Previous=links[i-1];
	}
	links[Cm->countOfRows-1]->Next = 0;
	links[0]->Previous = 0;

	delete[] links;

}

void check_number_of_ones ( struct CoveringMatrix *CM) {

	struct CLink *Columns;
	struct CLink *Rows;

	Rows=CM->RowLink;
	while ( Rows!=0) {
		//cout << " row " << Rows->ID << " : " << Rows->NumberOf1 << endl;
		Rows=Rows->Next;
	}

	Columns=CM->ColLink;
	while ( Columns!=0) {
		//cout << " row " << Columns->ID << " : " << Columns->NumberOf1 << endl;
		Columns=Columns->Next;
	}
}

void LagBaBCover::find_number_of_ones ( struct CoveringMatrix *CM) {
	
	int One;
	struct CLink *Columns,*colWithoutOnes;
	struct CLink *Rows;
	
	// reset counters
	CM->countOfOnes=0;
	Rows=CM->RowLink;
	while ( Rows!=0) {
		Rows->NumberOf1=0;
		Rows=Rows->Next;
	}

	Columns=CM->ColLink;
	while ( Columns!=0) {
		Columns->NumberOf1=0;  
		One=0;
		Rows=CM->RowLink;
		
		// find all rows covered by this column
		while ( Rows!=0) {
			One = CM->Array[Columns->ID][Rows->ID];
			if (One) {
				Rows->NumberOf1+=One;
				Columns->NumberOf1+=One;
				CM->countOfOnes++;
			}
			Rows=Rows->Next;
		}


		
		// if no rows are covered by this column
		if (Columns->NumberOf1==0) {
			colWithoutOnes=Columns;
			Columns=Columns->Next;
			remove_col_from_matrix(CM,colWithoutOnes);
		} else 
			Columns=Columns->Next;
	}

	check_number_of_ones(CM);

}

void LagBaBCover::solve_row_dominance ( struct CoveringMatrix *CM,struct header_of_solution_UCP *SolutionUCP ) {
	
	int Dominance;
	struct CLink *RowLink;
	struct CLink *Aux;
	struct CLink *ColLink;
	struct CLink *EssentialColumn = NULL;
	struct CLink *ComparedRow;
	int Dirty=0;
	int NumberOf1 = 0;


	find_number_of_ones(CM);
	ComparedRow=CM->RowLink;
	// Row to compare
	while ( ComparedRow != 0 ) { 
		
		RowLink=CM->RowLink;
		
		// Other rows are compared, looking for dominance
		while ( RowLink!=0 ) {

			Dominance=0;
			if (ComparedRow->NumberOf1 > RowLink->NumberOf1) {
				
				// There is no chance of dominance
				RowLink=RowLink->Next;
			} else {
				
				ColLink=CM->ColLink;
				NumberOf1=0;
				EssentialColumn=NULL;
				while ( ColLink!=0 ) {

					if ( CM->Array[ColLink->ID][ComparedRow->ID]==1 ) {
						if ( CM->Array[ColLink->ID][RowLink->ID] ==1 ) {   
							// If ComparedRow has ones at the same position as the second, domitates
							if ( Dominance!=-1 ) 
								Dominance=1;
						} else  { 
							// Dominance failed
							Dominance=-1; 
						}

						EssentialColumn=ColLink;
						NumberOf1++;
					}
					ColLink=ColLink->Next;
				}
				
				// If compare row with itself, it cannot be removed
				if (Dominance==1&&RowLink->ID!=ComparedRow->ID) {
					
					Aux=RowLink->Next;
					remove_row_from_matrix (CM,RowLink );
					Dirty = 1;
				} else 
					Aux=RowLink->Next;
				RowLink=Aux;
			}
		}

		if ( NumberOf1==1) {	
			
			/* If row is covered by only one column,that column is essential.
			* It will be added to solution and removed from matrix with
			* all rows it covers */
			add_to_solution ( SolutionUCP,EssentialColumn);
			Aux=ComparedRow->Next;
			remove_col_with_rows(EssentialColumn,CM);
			Dirty = 1;
			ComparedRow=Aux;
		} else 
			ComparedRow=ComparedRow->Next;
	}
	CM->Dirty=Dirty;
}

void LagBaBCover::solve_column_dominance ( struct CoveringMatrix *CM) {
	
	int Dominance;
	int Dirty=0;
	struct CLink *RowLink;
	struct CLink *Aux;
	struct CLink *ColLink;
	struct CLink *ComparedCol;

	find_number_of_ones(CM);

	ComparedCol=CM->ColLink;
	//Column to compare
	while ( ComparedCol != 0 ) {

		ColLink=CM->ColLink;
		while ( ColLink!=0 ) {
			
			if (ComparedCol->NumberOf1<ColLink->NumberOf1) {	
				//there is no chance dominance
				ColLink=ColLink->Next;
			} else {
				
				Aux=ColLink->Next;

				//To disable compare column with itself
				if ( ColLink->ID!=ComparedCol->ID ) {
					
					Dominance=0;

					RowLink=CM->RowLink;
					while ( RowLink!=0 ) {
						if (CM->Array[ColLink->ID][RowLink->ID] ==1 ) {
							if ( CM->Array[ComparedCol->ID][RowLink->ID] ==1 ){   
								
								// If ComparedCol has ones at the same position as second col
								// it has less ones than second column,will be removed
								Dominance=1;
							} else  {

								//dominance failed
								Dominance=0; 
								break; 
							}
						}
						RowLink=RowLink->Next;
					}

					//column will be removed here
					if (Dominance==1)  {

						//keep column with less cost than copared col in matrix
						if ( ColLink->Cost >= ComparedCol->Cost )  {
							remove_col_from_matrix (CM,ColLink );
							Dirty = 1;
						}
					}
				}
				ColLink=Aux;
			}
		}

		ComparedCol=ComparedCol->Next;
	}
	CM->Dirty=Dirty;
}

int LagBaBCover::solve_dominance ( struct CoveringMatrix *CM,struct header_of_solution_UCP *SolutionUCP ){ 
	
	int ColDom=0,return_Value=0;
	
	if ( CM->ColLink !=0 ){

		CM->Dirty=1;
		//Dominance will be solved as long as any change has been made
		
		while (CM->Dirty) {
			
			solve_row_dominance(CM,SolutionUCP);
			if (CM->Dirty) return_Value=1;
			if ( CM->Dirty==1||ColDom==0 ) {
				
				CM->Dirty=0;
				solve_column_dominance(CM);
				if (CM->Dirty) 
					return_Value =1;
				ColDom=1;
			}
		}
	}
	return return_Value;
}

void LagBaBCover::check_essentiality ( struct CoveringMatrix *CM,struct header_of_solution_UCP *SolutionUCP ){ 
	
	struct CLink *RowLink;
	struct CLink *ColLink;
	struct CLink *EssentialColumn = NULL;

	//cout << " finding number of 1s" << endl;
	find_number_of_ones(CM);

	RowLink=CM->RowLink;
	// Row to compare

	//cout << "checking rows" << endl;
	while ( RowLink != 0 ) { 

		if ( RowLink->NumberOf1==1) {	
			
			ColLink = CM->ColLink;

			while(ColLink!=0) {
				if ( CM->Array[ColLink->ID][RowLink->ID]==1 ) {
					EssentialColumn = ColLink;
					break;
				}
				ColLink=ColLink->Next;
			}

			
			/* If row is covered by only one column,that column is essential.
			* It will be added to solution and removed from matrix with
			* all rows it covers */
			add_to_solution ( SolutionUCP,EssentialColumn);
			remove_col_with_rows(EssentialColumn,CM);
			RowLink=CM->RowLink;
		} else 
			RowLink=RowLink->Next;
	}


}

float LagBaBCover::evaluate_column (struct CoveringMatrix *CM,struct CLink *ColumnLink ) {
	
	struct CLink *row_link;
	row_link=CM->RowLink;

	float SP=0;

	while (row_link!=0) {
		if ( CM->Array[ColumnLink->ID][row_link->ID]  )
			SP+=(float)1/row_link->NumberOf1;

		row_link=row_link->Next;
	}
	return SP;
}

struct CLink* LagBaBCover::select_col_contrib (struct CoveringMatrix *CM) {
	
	struct CLink *ColLink;
	struct CLink *SelectedColumn = NULL;
	float SP1,SP2;

	find_number_of_ones(CM);
	if ( CM->ColLink !=0 ) {
		ColLink=CM->ColLink;
		SelectedColumn=CM->ColLink;

		while ( ColLink!=0 ) {
			SP1=evaluate_column (CM,ColLink);
			SP1/=(float)ColLink->Cost;

			SP2=evaluate_column (CM,SelectedColumn);
			SP2/=SelectedColumn->Cost;

			if (SP1>SP2) 
				SelectedColumn=ColLink;

			ColLink=ColLink->Next;
		}	     
	}

	return SelectedColumn;
}

int LagBaBCover::branch_bound_lagrange ( struct CoveringMatrix *CM,struct header_of_solution_UCP *SolUCP,long double instance ) {

	struct CLink *SelectedColumn = NULL;
	struct CLink *SolutionColumn;
	struct CoveringMatrix *CopyCM2;
	struct CoveringMatrix *CopyCM3 = NULL;
	struct header_of_solution_UCP *SolUCP2;
	struct header_of_solution_UCP *SolUCP1;
	struct header_of_solution_UCP *SolLagrange = NULL;
	struct header_of_solution_UCP *Dominance_Solution;
	int lBound;
	int diff,Remove_Column_ID,Selected_Column_ID = -1;
	unsigned int i;
	unsigned char* solutionMask = NULL;
	unsigned char* lgSolution;
	float time;

	SCPView* lgview;
	LagrangeSolver* lgsolver;

	NumCalls++;
	

	//Solution given by essential columns
	Dominance_Solution=create_New_Empty_Solution();

	//cout << "B&B " << setw(8) << instance << ", Upper Bound: " << setw(6) << uBound;
	if(options.dominance==2 || (options.dominance==1 && instance==0)) {
		//cout << "solving dominance " << endl;
		solve_dominance(CM,Dominance_Solution);
	} else {
		//cout << "checking essentiality " << endl;
		check_essentiality(CM,Dominance_Solution);
	}
		

	//cout << "added dominance cost " << Dominance_Solution->Cost;
	CM->Cost+=Dominance_Solution->Cost;
	//cout << " ; Overall state cost: " << setw(6) << CM->Cost ;
	//cout << " ;  Solution cost (without domination): " << setw(6) << SolUCP->Cost << endl;
	append_to_solution (Dominance_Solution,SolUCP );

	//cout << "appended" << endl;
	
	if ( CM->ColLink!=0 && CM->RowLink!=0){

		//cout << "matrix is not empty" << endl;

		//needed for sorting
		//cout << "num of 1s" << endl;
		find_number_of_ones(CM);
		
		//cout << "sorting" << endl;
		sort_matrix(CM);

		
		time = GetMyCPUTime();

		// new problem view
		//cout << "new view" << endl;
		lgview = new SCPView(covMatrix,CM);

		// new lagrange structure
		
		//cout << "new solver" << endl;
		lgsolver = new LagrangeSolver(lgview);
		
		diff = (int)uBound - (int)CM->Cost;

		//solve using lagrange
		
		//cout << "solving" << endl;
	
		if(options.lagrange_simple_dual) {
		lgsolver->solve(options.lagrange_total_iter,0,2,options.lagrange_heur_iter,
			options.lagrange_init_coef,options.lagrange_min_coef,options.lagrange_iter,diff);
		//cout << "solved" << endl;
		} else {
		lgsolver->solve_dual(options.lagrange_total_iter,2,0,options.lagrange_init_coef,
		options.lagrange_min_coef,options.lagrange_iter,diff,
			options.lagrange_init_coef_dual,options.lagrange_min_coef_dual,options.lagrange_iter_dual);
		}
		
		//lBound =  lgsolver->get_lower_bound();
		//cout << "B&B " << setw(8) << instance << ", Upper Bound: " << setw(6) << uBound;
		//cout << " ; LB: " << setw(6) << lBound << ",maximum is: "; 
		//cout << setw(6) << uBound-CM->Cost << ". (found solution with cost: ";
		//cout << setw(6) << (int)lgsolver->get_lgsolution_cost() << ")" <<  endl;
		//return 0;
		
		//cout << "got it" << endl;

		Selected_Column_ID = lgsolver->get_recommended_column();

		lBound =  lgsolver->get_lower_bound();

		lagrangeTime+= GetMyCPUTime() - time;

		if(lBound >= (uBound-CM->Cost)){
			//LB is not acceptable
			//cout << "B&B " << setw(8) << instance << ", Upper Bound: " << setw(6) << uBound;
			//cout << " ; LB: " << setw(6) << lBound << " too high, maximum is: "; 
			//cout << setw(6) << uBound-CM->Cost << ". (found solution with cost: ";
			//cout << setw(6) << (int)lgsolver->get_lgsolution_cost() << ")" <<  endl;
			
			SolUCP->Cost=MaxCost;

			delete lgsolver;
			delete lgview;
			return 1;

		} else {
			// found LB is acceptable

			if (lBound == lgsolver->get_lgsolution_cost() )  {
				//found feasible solution with acceptable LB
				//this solution cannot be better

				SolLagrange = create_New_Empty_Solution();
				
				//cout << "B&B " << setw(8) << instance << ", Upper Bound: " << setw(6) << uBound;
				//cout << " ; LB: " << setw(6) << lBound << " leads to feasible solution with cost: "; 
				//cout << setw(6) << (int)lgsolver->get_lgsolution_cost();
				//cout << " Total cost is: " << setw(6) << (int)lgsolver->get_lgsolution_cost()+CM->Cost << endl;
				
				solutionMask = new unsigned char[covMatrix->GetColNum()];
				memset(solutionMask,0,covMatrix->GetColNum());

				/*printf("vypisuji column map: \n");
				for(i =0;i<problem->column_count;i++) 
				printf("%3d,",lgview->column_map[i]);
				printf("\n");*/

				lgSolution = lgsolver->get_best_solution();

				for(i =0;i<lgview->get_column_count();i++) {
					if(lgSolution[i]==1) {
						solutionMask[lgview->get_column_map(i)]=1;
					}
				}

				/*for(i =problem->get_column_count()-1;i>=0;i--)
					if(solutionMask[i]==1)
						cout << "i: " << i << " - " << solutionMask[i];*/

				add_to_solution_masked(SolLagrange,CM,solutionMask,covMatrix->GetColNum());
				append_to_solution(SolLagrange,SolUCP);

				//save new upper bound
				uBound= CM->Cost + (long)lgsolver->get_lgsolution_cost();

				//mark best solution as that in alghoritm
				bestLgSolutionUsed=0;

				bestLgState=instance;
				//cout << "B&B " << instance << "," << uBound << " - returning 0" << endl;
				delete lgsolver;
				delete lgview;
				return 0;
			} else 
				if(CM->Cost + lgsolver->get_lgsolution_cost() < uBound) {

					//found good solution, but still can be found better

					//cout << "B&B " << setw(8) << instance << ", Upper Bound: " << setw(6) << uBound;
					//cout << " ; LB: " << setw(6) << lBound << " leads to solution with cost: "; 
					//cout << setw(6) << (int)lgsolver->get_lgsolution_cost();
					//cout << " Total cost is: " << setw(6) << (int)lgsolver->get_lgsolution_cost()+CM->Cost << endl;

					//save new upper bound
					uBound= CM->Cost + (long)lgsolver->get_lgsolution_cost();
				
					// cout << "CM cost is " << CM->Cost << " and sol cost is " << SolUCP->Cost << endl;
			
					for(i=0;i<covMatrix->GetColNum();i++) {
						bestLgSolution[i]=0;
					}

					//store this solution somewhere else
					SolutionColumn=SolUCP->SetC;
					while(SolutionColumn!=NULL){
						bestLgSolution[SolutionColumn->ID]=1;
						//cout << "Column from solution " << SolutionColumn->ID << endl;
						SolutionColumn=SolutionColumn->Next;
					}

					lgSolution = lgsolver->get_best_solution();
					for(i =0;i<lgview->get_column_count();i++) {
						if(lgSolution[i]==1) {
							//cout << "Column from lagrange " << lgview->get_column_map(i) << endl;
							bestLgSolution[lgview->get_column_map(i)]=1;
						} 
					}
				
					/*for(i=0; i<problem->get_column_count();i++) {
					if(bestLgSolution[i]==1) {
					cout << "column " << i << endl;
					}
					}*/

					//mark best solution as that in bb
					bestLgSolutionCost=CM->Cost + (long)lgsolver->get_lgsolution_cost();
					bestLgSolutionUsed=1;
					bestLgState=instance;

					SolLagrange = create_New_Empty_Solution();
					solutionMask = new unsigned char[covMatrix->GetColNum()];
					memset(solutionMask,0,covMatrix->GetColNum());

					lgSolution = lgsolver->get_best_solution();

					for(i =0;i<lgview->get_column_count();i++) {
						if(lgSolution[i]==1) {
							solutionMask[lgview->get_column_map(i)]=1;
						}
					}
					CopyCM3=duplicate (CM);

				} else {
					
					//cout << "B&B " << setw(8) << instance << ", Upper Bound: " << setw(6) << uBound;
					//cout << " ; LB: " << setw(6) << lBound << " , maximum is: " << setw(6) << uBound-CM->Cost << endl;
				}
		}

		//cout << "deleting" << endl;
		delete lgsolver;
		delete lgview;

	}

	// CM is empty
	if ( CM->ColLink==0||CM->RowLink==0 )
	{

		//cout << "CM is empty" << endl;
		if ( CM->Cost<uBound ) {
			uBound=CM->Cost;
			//cout << "B&B " << setw(8) << instance << ", Upper Bound: " << setw(6) << uBound << " - found new best solution with cost " << CM->Cost << endl;
			//cout << "B&B " << instance << "," << uBound << " - returning 0" << endl;
			bestLgSolutionUsed=0;
			return 0;
		} else {
			//cout << "B&B " << setw(8) << instance << ", Upper Bound: " << setw(6) << uBound << " - found solution with cost " << CM->Cost << endl;
			//cout << "B&B " << instance << "," << uBound << " - returning 1" << endl;
			return 1;
		}
	}   else  {

		
		//cout << "going to select column" << endl;
		if(options.use_contribAdd==true) {
			
			SelectedColumn=select_col_contrib (CM);
			//cout << "selected column by ContribAdd is" << SelectedColumn->ID << endl;
			Selected_Column_ID = SelectedColumn->ID;
		}else {
			
			//cout << "selected column by lagrange is :" << Selected_Column_ID << endl;
			SelectedColumn=CM->ColLink;
			while ( SelectedColumn->ID!=Selected_Column_ID )
				SelectedColumn=SelectedColumn->Next;
			
			/*if(SelectedColumn==NULL) {
				cout << "column not found!, WTF?" << endl;
			} else {
				cout << "selected " << SelectedColumn->ID << endl;
			}*/

		}
		
		//cout << "column selected by B&B : " << Remove_Column_ID << endl;
		//cout << " selected column " << Selected_Column_ID << " with cost " << SelectedColumn->Cost << endl;

		Remove_Column_ID = Selected_Column_ID;

		CopyCM2=duplicate (CM);

		//Solution given by first call B&B
		SolUCP1=create_New_Empty_Solution();
		
		//First branch

		//Add selected column to solution		

		add_to_solution ( SolUCP1,SelectedColumn);
		CM->Cost+=SolUCP1->Cost;
		 
		//Remove selected column from matrix
		remove_col_with_rows ( SelectedColumn,CM );

		if ( branch_bound_lagrange(CM,SolUCP1,2*instance+1)==1 ) 
			SolUCP1->Cost=MaxCost;
		
		//Second branch
		
		//Remove selected column from matrix
		remove_col_from_matrix ( CopyCM2,Remove_Column_ID);
		SolUCP2=create_New_Empty_Solution();

		if ( branch_bound_lagrange(CopyCM2,SolUCP2,2*instance+2)==1 )
			SolUCP2->Cost=MaxCost;
		release (CopyCM2);   
	}


	//Solution
	if ( SolUCP1->Cost < SolUCP2->Cost ) {
		//Solution 1 was succesfull
		release (SolUCP2);
		append_to_solution ( SolUCP1,SolUCP );
		//cout << "B&B " << instance << "," << uBound << " - returning 0" << endl;
		//if(CopyCM3 != NULL ) release(CopyCM3);
		//if(SolLagrange!= NULL) release(SolLagrange);
		return 0;
	} else {
		//Solution 2 was succesfull
		if ( SolUCP2->Cost!=MaxCost ) {
			release (SolUCP1);
			append_to_solution ( SolUCP2,SolUCP );

			//if(CopyCM3 != NULL ) release(CopyCM3);
			//if(SolLagrange!= NULL) release(SolLagrange);
			//cout << "B&B " << instance << "," << uBound << " - returning 0" << endl;
			return 0;
		} else {

			if ((bestLgSolutionUsed==1) && (bestLgState==instance)) {
				release (SolUCP1);
				release (SolUCP2);
				
				add_to_solution_masked(SolLagrange,CopyCM3,solutionMask,covMatrix->GetColNum());
				append_to_solution(SolLagrange,SolUCP);

				//mark best solution as that in alghoritm
				bestLgSolutionUsed=0;
				bestLgState=instance;
				//cout << "B&B " << instance << "," << uBound << " - returning 0" << endl;

				//release(CopyCM3);
				//release(SolLagrange);


				return 0;
			}
			//No solution was succesfull
			release (SolUCP1);
			release (SolUCP2);
			//cout << "B&B " << instance << "," << uBound << " - returning 1" << endl;
			//if(CopyCM3 != NULL ) release(CopyCM3);
			//if(SolLagrange!= NULL) release(SolLagrange);
			return 1;
		}
	}
}

struct CoveringMatrix *LagBaBCover::rebuild_matrix (CovMatrix *CVM)
{
	int I,J;
	int * Costs;

	struct CoveringMatrix *Cm;

	Cm=new CoveringMatrix;
	if ( Cm==NULL ) throw RuntimeException( "Aura: Not enough memory" );
	Cm->Cost=0;
	Cm->countOfCols=0;
	Cm->countOfRows=0;
	Cm->countOfOnes=0;

	// First load dimensions of array to allocate memory
	Cm->countOfCols=CVM->GetColNum();
	Cm->countOfRows=CVM->GetRowNum();            
	Costs=new int [Cm->countOfCols];
	if ( Costs==NULL) throw RuntimeException( "Aura: Not enough memory" );

	Cm->Array=new char * [Cm->countOfCols];
	if ( Cm->Array==0 ) throw RuntimeException( "Aura: Not enough memory" );

	for (I=0;I<Cm->countOfCols;I++)
	{
		Cm->Array[I]=new char [Cm->countOfRows]; 
		if ( Cm->Array[I]==0 ) throw RuntimeException( "Aura: Not enough memory" );
	}

	for (I=0;I<Cm->countOfCols;I++)
	{
		Costs[I]=CVM->GetCost(I);

	}

	for (J=0;J<Cm->countOfRows;J++)

		for (I=0;I<Cm->countOfCols;I++)
			Cm->Array[I][J]=CVM->GetValue(J,I);

	create_matrix ( Cm,Costs);

	delete [] Costs;

	return Cm;

}


list<unsigned int> *LagBaBCover::rebuild_solution (struct header_of_solution_UCP *SolutionUCP)
{
	struct CLink *Aux;
	list<unsigned int> *Solution;
	Aux=SolutionUCP->SetC;

	Solution = new list<unsigned int>;		// !!! this is not deallocated by this class !!!

	while (Aux!=0)
	{
		Solution->push_back(Aux->ID);
		//cout << Aux->ID << endl;
		Aux=Aux->Next;

	}
	return Solution;
}

LagBaBCover::LagBaBCover (CovMatrix *_matrix )
{
	uBound=0;
	NumCalls=0;

	//Create covering matrix from original matrix structure
	covering_matrix=rebuild_matrix( _matrix );

	covMatrix = _matrix;

	//Create new header of solution
	solution_UCP=create_New_Empty_Solution();

	bestLgSolution =  new int[covMatrix->GetColNum()];
	bestLgSolutionCost = 0;
	bestLgSolutionUsed = 0;

	lagrangeTime=0;

	//isAllCovered = false;


}
LagBaBCover::~LagBaBCover ()
{
	release (covering_matrix);
	release (solution_UCP);

}

list<unsigned int> *LagBaBCover::Run()
{
	unsigned int i;
	list<unsigned int> *Solution;

	Solution = new list<unsigned int>;		// !!! this is not deallocated by this class !!!

	branch_bound_lagrange ( covering_matrix,solution_UCP,0);

	if(bestLgSolutionUsed==1) {
		//cout << "stored best solution with cost " << bestLgSolutionCost << endl;

		for(i=0;i<covMatrix->GetColNum();i++) {
			if(bestLgSolution[i]==1) {
				
				Solution->push_back(i);
				//cout << "column " << i << endl;
			}
		}

	} else {
		//cout << "got solution" << endl;
		Solution=rebuild_solution (solution_UCP);
	}
	//cout << "Number of expanded states: " << NumCalls << endl;
	//cout << "Time spend by solving lagrange " << lagrangeTime << endl;

	//isAllCovered = check_coverage(covering_matrix,Solution);

	return Solution;
}

bool LagBaBCover::check_coverage(struct CoveringMatrix *CM, list<unsigned int> solution){

	bool* rowIsCovered = new bool[CM->countOfRows];
	int i;
	list<unsigned int>::iterator sol_it;

	for(i=0;i<CM->countOfRows;i++)
			rowIsCovered[i] = false;

	for ( sol_it = solution.begin(); sol_it != solution.end(); sol_it++ ){
		
		for(i=0;i<CM->countOfRows;i++)
			if(CM->Array[*sol_it][i]==1)
				rowIsCovered[i] = true;		
		
	}


	for(i=0;i<CM->countOfRows;i++)
			if(rowIsCovered[i] == false)
				return false;
	return true;

}
