modmaxlp.cpp

用线性规划求取复杂网络模块度算法的关键实现文件。

/*
 ****************************
 * Modularit Maximization LP
 * Implementation of the LP for modularity maximization problem.
 * Tested with CPLEX 9.0, gcc 3.2
 * Input	: Modularity Matrix
 * Output	: Distance Matrix
 * 
 * IMPORTANT: Assumes the symmetry in the input matrix.
 * 
 * Author	: Gaurav Agarwal
 * Email	: [e-mail]
 * Date		: 31 March, 2007
 ****************************
 */
 
#include <iostream>
#include <fstream>
#include <string>
#include <ilcplex/ilocplex.h>
using std::cerr;
using std::cout;
using std::endl;
using std::ifstream;
using std::ofstream;
ILOSTLBEGIN

#define m_modularity_matrix(i,j) m_modularity_matrix[i*m_total_nodes +j]
#define m_dist_matrix(i,j) m_dist_matrix[i*m_total_nodes +j]
#define m_neg_dist_matrix(i,j) m_neg_dist_matrix[i*m_total_nodes +j]
#define m_pos_dist_matrix(i,j) m_pos_dist_matrix[i*m_total_nodes +j]

typedef IloArray<IloNumVarArray> NumVarMatrix;
typedef IloArray<IloNumArray> FloatMatrix;

class ModMaxLP
{
	private:
		int m_total_nodes;
		double* m_mod_matrix;
		double* m_dist_matrix;
		double* m_pos_dist_matrix;
		double* m_neg_dist_matrix;
	
		void read_modularity_matrix(char* distFile);
		void build_model(IloModel model,NumVarMatrix X);
		void dump2File(char* opFileName,FloatMatrix dist_mat);
		void swap(int* i,int* j)
		{
			int temp=*i;
			*i=*j;
			*j=temp;
		}
		void addTriangleConstraint(IloModel model,NumVarMatrix X,int i,int j,int k);
	public:
		void solve_lp(char* simFile, char* distFile);
};

/*
 * This method will read the similarity file in a double pointer. 
 * Will allocate the necessary space before doing anything.
 */
void ModMaxLP::read_modularity_matrix(char* modFile)
{
	cout << "\tGoing to read modularity file: " << modFile << endl;
	ifstream indata;
	indata.open(modFile);
	 if(!indata) 
	 { 
      	  cerr << "\tError: file " << modFile << "could not be opened" << endl;
      	  exit(1);
	 }
	 indata >> m_total_nodes;
	 cout<< "\t\tTotal nodes: "<< m_total_nodes << endl;
	 m_modularity_matrix = new double [m_total_nodes * m_total_nodes];
	 m_dist_matrix = new double [m_total_nodes * m_total_nodes];
 	 m_pos_dist_matrix = new double [m_total_nodes * m_total_nodes];
 	 m_neg_dist_matrix = new double [m_total_nodes * m_total_nodes];
	 //read file now
	 for(int i=0;i<m_total_nodes;i++)
	 {
	 	for(int j=0;j<m_total_nodes;j++)
	 	{
	 		indata>>m_mod_matrix(i,j);
	 	}
	 }
	 indata.close();
	 
	 for(int i=0;i<m_total_nodes;i++)
	 {
	 	for(int j=0;j<m_total_nodes;j++)
	 	{
	 		if(m_mod_matrix(i,j) < 0)
	 		{
	 			m_neg_dist_matrix(i,j) = abs(m_mod_matrix(i,j));
	 			m_pos_dist_matrix(i,j) = 0;
	 		}
	 		else
	 		{
	 			m_pos_dist_matrix(i,j) = m_mod_matrix(i,j);
	 			m_neg_dist_matrix(i,j) = 0;
	 		}
	 	}
	 }
	 delete(m_mod_matrix);
	 cout<< "\tFile read successfully.."<< endl<<endl;
}
/*
 * Populates the CPLEX model for the mindisagree LP
 */
void ModMaxLP::build_model(IloModel model,NumVarMatrix X)
{
	cout<<"\t++Building Model++"<<endl;
	IloEnv env = model.getEnv();;
	cout<<"\t\tCreating variables.."<<endl;
	FloatMatrix Neg_coeff(env,m_total_nodes);
	FloatMatrix Pos_coeff(env,m_total_nodes);
	
	for(int i=0;i<m_total_nodes;i++)
	{
		X[i] = IloNumVarArray(env,i+1,0,1,ILOFLOAT); //allocating space for lower triangle only
		Neg_coeff[i] = IloNumArray(env,m_total_nodes);
		Pos_coeff[i] = IloNumArray(env,m_total_nodes);
	}
	
	cout<<"\t\tVaraibles created.."<<endl<<endl;
	
	cout <<"\t\tPopulating mod matrix (+-).."<<endl;
	for(int i=0;i<m_total_nodes;i++)
	{
		for(int j=0;j<m_total_nodes;j++)
		{
			Neg_coeff[i][j]=m_neg_dist_matrix(i,j);
			Pos_coeff[i][j]=m_pos_dist_matrix(i,j);
		}
	}
	cout <<"\t\tDone Populating mod matrix.." << endl<<endl;
	
	cout <<"\t\tBuilding objective.." << endl;
	IloExpr expr(env);
	for (int i=0;i<m_total_nodes;i++)
	{
		for(int j=0;j<=i;j++)//consider lower triangle only
		{
		  expr+=Pos_coeff[i][j]*X[i][j];
		  expr+=Neg_coeff[i][j]*(1-X[i][j]);
		}
	}
	model.add(IloMinimize(env, expr));
	expr.end();
	cout <<"\t\tDone building objective.." << endl<<endl;
	
	cout <<"\t\tBuilding constraints.."<<endl;
	int count=0;
	for(int i=0;i<m_total_nodes-1;i++)
	{
		for(int j=i+1;j<m_total_nodes-1;j++)
		{
			for(int k=j+1;k<m_total_nodes;k++)
			{
  		            //combination 1
			    addTriangleConstraint(model,X,i,j,k);	
				
			    //combination 2
			    addTriangleConstraint(model,X,j,k,i);

			    //combination 3
			    addTriangleConstraint(model,X,k,i,j);
			    count+=3;
			}
		}
	}
	for(int i=0;i<m_total_nodes;i++)
	  {
	    // diagonal == 0 constraint
	    model.add(X[i][i]==0);
	  }
	count+=m_total_nodes;
	  
	cout << "\t\tDone building constraints.."<<endl;
	cout << "\t\t"<<count<<" constraints added.."<<endl;
	cout << "\t++Done building the model++"<<endl<<endl;
}

void ModMaxLP::addTriangleConstraint(IloModel model,NumVarMatrix X, int i,int j,int k)
 {
 	int i1=i,i2=i,j1=j,j2=j,k1=k,k2=k;
				 
	if(i1<j1) 
	{
	 	swap(&i1,&j1);
	}
	if(j2<k1)
	{
	 	swap(&j2,&k1);
	}
	if(i2<k2)
	{
	 	swap(&i2,&k2);
	}
			 
	if(i1<j1 || j2<k1 || i2<k2)
	{
	 	printf("dissastor(1):i1:%d i2:%d j1:%d j2:%d k1:%d k2:%d\n",i1,i2,j1,j2,k1,k2);
	 	printf("i:%d j:%d  k:%d\n",i,j,k); 
	}
	model.add(X[i1][j1] + X[j2][k1] >= X[i2][k2]);		
 }

/*
 * Writes the distance matrix to a file
 */
void ModMaxLP::dump2File(char* distFile,FloatMatrix dist_mat)
{
	cout<<"\tBegining to write the output to :"<<distFile<<endl;
	ofstream outdata;
	outdata.precision(10);
	outdata.open(distFile);
	if(!outdata) 
	 { 
      	cerr << "\tError: file " << distFile << "could not be opened" << endl;
      	exit(1);
	 }
	for(int i=0;i<m_total_nodes;i++)
	{
		for(int j=0;j<m_total_nodes;j++)
		{
			//, check i<j, revert
			if(i<j)
			{
				outdata<<dist_mat[j][i]<<" ";
			}
			else
			{
				outdata<<dist_mat[i][j]<<" ";
			}
		}
		outdata<<endl;
	}
	outdata.close();
	cout<<"\tDone writing the output.."<<endl;
}

/*
 * Method which does all the work thru sub-methods.
 * Reads similarity matrix.
 * Creates Model.
 * Solves LP
 * Dumps output
 */
void ModMaxLP::solve_lp(char* simFile,char* distFile)
{
	cout<<">>>>>>>>>>>>>>>>>>>Starting<<<<<<<<<<<<<<<<<<<"<<endl<<endl;
	IloEnv env;
	IloModel model(env);
	read_modularity_matrix(simFile);
	NumVarMatrix X(env,m_total_nodes);
	
	build_model(model,X);
	
	    //free some memory here..even if lil bit (we are in C-liperry territory) ;)
	free(m_pos_dist_matrix);
	free(m_neg_dist_matrix);
	IloCplex cplex(model);
	//The following flag causes to reduce the CPLEX memory requirement but produces very small errors in the order of 10^-5
	//cplex.setParam(IloCplex::FinalFactor, false);
	try
	  {
	    cplex.solve();
	  }
	catch(IloException& e)
	  {
	    cout<<e;
	    cout<<e.getMessage();
	  }
	if(cplex.getStatus() == IloAlgorithm::Infeasible)
	{
		cout << "\tInfeasible.."<<endl;
	}
	else if(cplex.getStatus() == IloAlgorithm::Optimal)
	{
		cout << "\tSolved optimally.."<<endl;
		cout << "\tOptimal value: "<<cplex.getObjValue()<<endl;
	}
	else
	{
		cout<<"\tUnknown terminating status: "<<cplex.getStatus()<<endl;
	}
	
	cout<<"\tGoing to retreive variables.."<<endl;
	FloatMatrix dist_sol(env,m_total_nodes);
	for(int i=0;i<m_total_nodes;i++)
	{
	    dist_sol[i]=IloNumArray(env,i+1);
	}
	for(int i=0;i<m_total_nodes;i++)
	{
		cplex.getValues(dist_sol[i],X[i]);
	}
	cout<<"\tVariables retreived.."<<endl;
	dump2File(distFile,dist_sol);
	cout<<endl<<">>>>>>>>>>>>>>>>>>>Completed<<<<<<<<<<<<<<<<<<<"<<endl;
}


int main(int argc, char* argv[])
{
	if(argc!=3)
	{
		cout<<"Wrong usage.."<<endl;
		cout<<"Usage:Solver <modularity_matrix_file> <output_dist_mat_file (writable)>"<<endl;
	}
	else
	{
		ModMaxLP ob;
		ob.solve_lp(argv[1],argv[2]);
	}
}