/*  
 
MAXIMAL INDEPENDENT SETS 
 
26 SEPTEMBER 2001 
 
PRECONDITIONS:  
 
0. ILOG Solver 5.0. 
1. A hypergraph read in in the following format: 
 
************************************* 
number_of_vertices 
arity_1 vertex_1 vertex_2 ... 
arity_2 vertex_1 vertex_8 ... 
... 
-1 
************************************* 
 
Here number_of_vertices is the overall  
number of vertices in the hypergraph. 
 
Each following line (except the last one containing -1) 
corresponds to an hyperedge, the first number is its arity 
(i.e., how many vertices it has), the following numbers are  
the vertices of this edge. 
 
2. output file. 
3. Size k of the desired (maximal) independent set (i.e., 
it should be known how many vertices make up the desired set).
4. Whether or not to consider maximality <0/1>;
5. The index of encoding <1..6>:
 Encodings 1 to 3 : independence is implemented using the IloSum function;
                    they have no difference if maximality is not considered;
 Encodings 4 to 6 : independence is implemented using the IloDistribute function;
		    they have no difference if maximality is not considered;

 If maximality is considered:
 Encoding 1 or 4, called max1 : maximality is encoded using P <-> Q;
 Encoding 2 or 5, called max2 : maximality is encoded using (P -> Q) ^ (Q -> P);
 Encoding 3 or 6, called max3 : maximality is encoded using (P ^ Q) V (not(P) ^ not(Q));
 
 Here P is v_i = 1 (v_i is the Boolean variable corresponding to a vertex of
 the hypergraph, it is equal to 1 iff the vertex is in the desired set), 
 whereas Q is the statement that 
 for the given hyperedge e the sum of the variables v_k (k = [1..arity(e)]
 such that k!=i) is less than the arity of this hyperedge, i.e. arity(e).

POSTCONDITION: The first solution, a (maximal) independent set of size k
and some statistics (number of fails, number of nodes, etc.)

*/ 
 
#include <ilsolver/ilosolver.h> 
#include <time.h> 
 
ILOSTLBEGIN  

//max number of edges in the hypergraph 
const nbEdgesMax = 1000; 
//number of backtracks and number of nodes in the search tree
unsigned long int NbFails, NbNodes;

const IloInt max1 = 1;
const IloInt max2 = 2;
const IloInt max3 = 3;


int main(int argc, char *argv[]) { 

	NbFails = 0;
	NbNodes = 0;

	IloEnv env; 
	try { 
 
		/////////////////////////////////////////////////////////////////////////// 
		// 
		// READING DATA 
		// 
		/////////////////////////////////////////////////////////////////////////// 
 
		IloModel mdl(env); 
		int i=0,j=0,k=0; 
		ifstream infile; 
		ofstream outfile; 
		char * datafile; 
		char * output; 
		if (argc ==6) { 
			datafile = argv[1]; 
			output   = argv[2]; 
		} 
		else { 
			env.out() << "Usage: *.exe datafile resfile size maximalityOn <0/1> encoding <1..6>"<<endl; 
			return 0; 
		} 
 
		IloInt SizeOfSet = atoi (argv[3]); 
		IloBool MaximalityConstraintOn = atoi (argv[4]); 
		IloInt EncodingNumber = atoi (argv[5]);

		//a flag for outputting additional information
		IloBool isOutputOn = IloTrue;

		IloBool isSummationOn = IloFalse;
		IloInt Encoding=0;
		
		if (EncodingNumber<=3) isSummationOn=IloTrue;

		if (EncodingNumber==1||EncodingNumber==4) Encoding=max1;
		if (EncodingNumber==2||EncodingNumber==5) Encoding=max2;
 		if (EncodingNumber==3||EncodingNumber==6) Encoding=max3;

		infile.open(datafile); 
		outfile.open(output, ios::out); 
 
		if (!infile) {  
			env.out() << "File not found or not specified: " << datafile << endl; 
			env.end(); 
			return 0; 
		} 
		if (!outfile) {  
			env.out() << "Cannot open file for output: " << output << endl; 
			env.end(); 
			return 0; 
		} 
 
		int nbVertices; 
		const int GlobSize = nbEdgesMax; 
		infile >> nbVertices; 
		IloNumArray vertices(env,nbVertices); 
 
		for (i=0;i<nbVertices;i++) { 
			vertices[i] = i+1; 
		} 
 
		int nbEdges = 0; 
 
		//2D array of data. The hypergraph is read in from a file,
		//see comments above.

		int **DataArray; 
		DataArray = new int *[GlobSize]; 
		for (i=0;i<GlobSize;i++) { 
			DataArray[i]=new int[GlobSize]; 
		} 
		for (i=0;i<GlobSize;i++){ 
			for (j=0;j<GlobSize;j++){ 
			     DataArray[i][j]=0; 
			} 
		} 
		
		/////////////////////////////////////////////////////////////////////////////
		//
		// READ THE NUMBER OF EDGES     
		//
		/////////////////////////////////////////////////////////////////////////////

		int currArity, count = 0; 
		j = 0; 
 
		//arityArr is the array of arities of the hyperedges.
		//It is needed in order to avoid using the getSize() function
		//every time the size of an array of variables is required.

		IloNumArray arityArr(env, nbEdgesMax, 0, nbVertices, ILOINT); 
		for (i=0;i<nbEdgesMax;i++) { 
			arityArr[i]=0; 
		} 
 
		while (true) { 
			infile >> currArity;
			if (currArity==-1) break; 
			arityArr[j] = currArity; 
			nbEdges++; 
			for (int i=0;i<currArity;i++) { 
				infile >> DataArray[j][i]; 
			}
			j++; 
		} 
	
		if (isOutputOn)
			cout<<"Your hypergraph has "<<nbVertices<<" vertices and "
			    <<nbEdges<<" edges."<<endl;

		////////////////////////////////////////////////////////////////////////////
		//
		// 2D ARRAY FOR VERTEX -> EDGE MAPPING 
		//
		//////////////////////////////////////////////////////////////////////////// 

		//If vertex i is in hyperedge j,
		//EdgeLists[i][j]=1, otherwise 0.
		int **EdgeLists; 
		EdgeLists = new int *[nbVertices]; 
		for (i=0;i<nbVertices;i++) { 
			EdgeLists[i]=new int[nbEdges]; 
		} 
 
		for (i=0;i<nbVertices;i++){ 
			for (j=0;j<nbEdges;j++){ 
				EdgeLists[i][j]=0; 
			} 
		} 
 
		for (i=0;i<nbVertices;i++){ 
			for (j=0;j<GlobSize;j++){ 
				for (k=0;k<GlobSize;k++) { 
					if (DataArray[j][k]==i+1) { 
						EdgeLists[i][j]=1; 
					} 
				} 
			} 
		} 

		if (isOutputOn) {
			outfile<<"Edges - columns, rows - vertices:"<<endl;
			outfile<<endl;
			for (i=0;i<nbVertices;i++){ 
				for (j=0;j<nbEdges;j++){ 
					outfile<<EdgeLists[i][j]<<" "; 
				} 
				outfile<<endl;
			} 
			outfile<<endl;
			outfile<<"Arities array:"<<endl;
			outfile<<endl;
			for (j=0;j<nbEdges;j++){ 
				outfile<<arityArr[j]<<" "; 
			} 
			outfile<<endl;
		}		
		 
		/////////////////////////////////////////////////////////////////////////////
		//
		// VARIABLES    
		//
		/////////////////////////////////////////////////////////////////////////////
		
		//nbVertices Boolean variables:  
		//1 - a vertex is in the set, 0 - otherwise. 
		IloNumVarArray vertexPresent(env,nbVertices,0,1,ILOINT); 

		/////////////////////////////////////////////////////////////////////////////
		//
		// CONSTRAINTS     
		//
		/////////////////////////////////////////////////////////////////////////////

		///////////////////////////////////////////////////////////
		//
		// 1. constraint on the size of the maximal independent set 
		//
		///////////////////////////////////////////////////////////

		if (isSummationOn) {
			mdl.add(IloSum(vertexPresent)==SizeOfSet); 
		}
		else {
			IloNumArray values(env,2,0,1);
			IloNumVarArray cards(env,2,0,vertexPresent.getSize(),ILOINT);
			mdl.add(IloDistribute(env,cards,values,vertexPresent));
			mdl.add(cards[1]==SizeOfSet);
		}
 
		///////////////////////////////////////////////////////////
		//
		// 2. independence constraints 
		//
		///////////////////////////////////////////////////////////

		IloExpr expr(env);
 
		if (isSummationOn) {
			for (j=0;j<nbEdges;j++) { 
				IloNumVarArray edgeVarArray(env,arityArr[j]);

				int k=0;
				for (i=0;i<nbVertices;i++) {
					if (EdgeLists[i][j]) {
						edgeVarArray[k]=vertexPresent[i];
						k++;
					}
				}
				mdl.add(IloSum(edgeVarArray)<arityArr[j]);
			}
		}
		else {
			for (j=0;j<nbEdges;j++) { 
				IloNumArray values(env,2,0,1);
				IloNumVarArray cards(env,2,0,vertexPresent.getSize(),ILOINT);
				IloNumVarArray edgeVarArray(env,arityArr[j]);

				int k=0;
				for (i=0;i<nbVertices;i++) {
					if (EdgeLists[i][j]) {
						edgeVarArray[k]=vertexPresent[i];
						k++;
					}
				}

				mdl.add(IloDistribute(env,cards,values,edgeVarArray));
				mdl.add(cards[0]>=1);
			}
		}

		///////////////////////////////////////////////////////////
		//
		// 3. maximality constraints
		//
		///////////////////////////////////////////////////////////

		if (MaximalityConstraintOn) {
			IloConstraint constraint_Q, constraint_notQ, c, not_c; 
			IloBool flag;

			// encoding max1: p <-> q, see comments above;
			// encoding max2: (p -> q) && (q->p);
			// encoding max3: (p&&q)||(not(p)&&not(q));
						
			expr.remove(vertexPresent);

			if (Encoding==max1 || Encoding==max2) {
				for (i=0;i<nbVertices;i++) {
					flag = IloTrue;
					for (j=0;j<nbEdges;j++) {
						if (EdgeLists[i][j]) {
							for (int k=0;k<nbVertices;k++) {
								if (i!=k && EdgeLists[k][j]) expr += vertexPresent[k];
							}
							if (flag) {
								constraint_Q = (expr<arityArr[j]-1);
								flag = IloFalse;
							}
							else {
								constraint_Q = constraint_Q && (expr<arityArr[j]-1);
							}
							expr.remove(vertexPresent);
						}
					}

					if (Encoding==max1) {
						mdl.add((vertexPresent[i]==1)==constraint_Q);
					}
					if (Encoding==max2) {
						mdl.add((vertexPresent[i]==1)<=constraint_Q);
						mdl.add(constraint_Q<=(vertexPresent[i]==1));
					}
				}
			}


			if (Encoding==max3) {

				expr.remove(vertexPresent);
				for (i=0;i<nbVertices;i++) {
					flag = IloTrue;
					for (j=0;j<nbEdges;j++) {
						if (EdgeLists[i][j]) {
							for (int k=0;k<nbVertices;k++) {
								if (i!=k && EdgeLists[k][j]) expr += vertexPresent[k];
							}
							if (flag) {
								not_c = (expr==arityArr[j]-1);
								c = (expr<arityArr[j]-1);
								flag = IloFalse;
							}
							else {
								not_c = not_c || (expr==arityArr[j]-1);
								c = c && (expr<arityArr[j]-1);
							}
							expr.remove(vertexPresent);
						}
					}
					mdl.add((vertexPresent[i]==0)&&not_c || (vertexPresent[i]==1)&&c);
				}
			}
		}

		/////////////////////////////////////////////////////////////////////////////
		//
		// SOLVING THE PROBLEM             
		//
		/////////////////////////////////////////////////////////////////////////////

		IloSolver solver(mdl); 
		IloGoal goal = IloGenerate(env,vertexPresent); 
		IlcSearch search = solver.newSearch(goal); 
		int ctr = 1; 
		
		IloBool NoSolutions = IloTrue;

		cout<<endl;
		cout<<datafile<<": size: "<<SizeOfSet<<" maximalityOn: "<<MaximalityConstraintOn
			          <<" Encoding: "<<EncodingNumber<<endl;
		if (search.next()) { 
			NoSolutions = IloFalse;

			solver.out()<<"Solution "<<ctr<<" :"<<endl; 
			outfile<<endl;
			outfile<<"Solution "<<ctr<<" : "<<endl; 
			for (i=0;i<nbVertices;i++) { 
				if (solver.getIntValue(vertexPresent[i])) {  
					solver.out()<<" "<<i+1; 
					outfile<<" "<<i+1; 
				} 
			} 
			solver.out()<<endl; 
			outfile<<endl; 

			solver.printInformation(); 
			NbNodes +=solver.getNumberOfChoicePoints();
			NbFails +=solver.getNumberOfFails();
		}  

		if (NoSolutions) {
			solver.printInformation();
			NbNodes +=solver.getNumberOfChoicePoints();
			NbFails +=solver.getNumberOfFails();
		}

		/////////////////////////////////////////////////////////////////////////////
		//
		// CLEANING UP                  
		//
		/////////////////////////////////////////////////////////////////////////////

		for (i=0;i<nbVertices;i++) { 
			delete EdgeLists[i]; 
		} 
		delete EdgeLists; 
		 
		for (i=0;i<GlobSize;i++) { 
			delete DataArray[i]; 
		} 
		delete DataArray; 
		 
	} 
	catch (IloException& ex) { 
		cerr <<"Error: "<< ex << endl; 
	} 
	env.end(); 
	return 0; 
}