fastcommunity_mh.cc

大型复杂网络社区划分的快速算法

		// ---------------------------------
		// Find largest dQ
		if (ioparm.textFlag > 0) { h->printHeapTop10(); cout << endl; }
		dQmax = h->popMaximum();					// select maximum dQ_ij // convention: insert i into j
		if (dQmax.m < -4000000000.0) { break; }		// no more joins possible
		cout << "Q["<<t-1<<"] = "<<Q[t-1];
		
		// ---------------------------------
		// Merge the chosen communities
		cout << "\tdQ = " << dQmax.m << "\t  |H| = " << h->heapSize() << "\n";
		if (dq[dQmax.i].v == NULL || dq[dQmax.j].v == NULL) {
			cout << "WARNING: invalid join (" << dQmax.i << " " << dQmax.j << ") found at top of heap\n"; cin >> pauseme;
		}
		isupport = dq[dQmax.i].v->returnNodecount();
		jsupport = dq[dQmax.j].v->returnNodecount();
		if (isupport < jsupport) {
			cout << "  join: " << dQmax.i << " -> " << dQmax.j << "\t";
			cout << "(" << isupport << " -> " << jsupport << ")\n";
			mergeCommunities(dQmax.i, dQmax.j);	// merge community i into community j
			joins[t].x = dQmax.i;				// record merge of i(x) into j(y)
			joins[t].y = dQmax.j;				// 
		} else {								// 
			cout << "  join: " << dQmax.i << " <- " << dQmax.j << "\t";
			cout << "(" << isupport << " <- " << jsupport << ")\n";
			dq[dQmax.i].heap_ptr = dq[dQmax.j].heap_ptr; // take community j's heap pointer
			dq[dQmax.i].heap_ptr->i = dQmax.i;			//   mark it as i's
			dq[dQmax.i].heap_ptr->j = dQmax.j;			//   mark it as i's
			mergeCommunities(dQmax.j, dQmax.i);	// merge community j into community i
			joins[t].x = dQmax.j;				// record merge of j(x) into i(y)
			joins[t].y = dQmax.i;				// 
		}									// 
		Q[t] = dQmax.m + Q[t-1];					// record Q(t)
		
		// ---------------------------------
		// Record join to file
		ofstream fjoins(ioparm.f_joins.c_str(), ios::app);   // open file for writing the next join
		fjoins << joins[t].x-1 << "\t" << joins[t].y-1 << "\t";	// convert to external format
		if ((Q[t] > 0.0 && Q[t] < 0.0000000000001) || (Q[t] < 0.0 && Q[t] > -0.0000000000001))
			{ fjoins << 0.0; } else { fjoins << Q[t]; }
		fjoins << "\t" << t << "\n";
		fjoins.close();
		// Note that it is the .joins file which contains both the dendrogram and the corresponding
		// Q values. The file format is tab-delimited columns of data, where the columns are:
		// 1. the community which grows
		// 2. the community which was absorbed
		// 3. the modularity value Q after the join
		// 4. the time step value
		
		// ---------------------------------
		// If cutstep valid, then do some work
		if (t <= ioparm.cutstep) { groupListsUpdate(joins[t].x, joins[t].y); }
		if (t == ioparm.cutstep) { recordNetwork(); recordGroupLists(); groupListsStats(); }

		// ---------------------------------
		// Record the support data to file
		if (ioparm.suppFlag) {
			dqSupport();
			ofstream fsupp(ioparm.f_support.c_str(), ios::app);
			// time   remaining support   mean support   support_i --   support_j
			fsupp << t << "\t" << supportTot << "\t" << supportAve << "\t" << isupport;
			if (isupport < jsupport) { fsupp  << "\t->\t"; }
			else { fsupp << "\t<-\t"; }
			fsupp << jsupport << "\n";
			fsupp.close();
		}
		if (Q[t] > Qmax.y) { Qmax.y = Q[t]; Qmax.x = t; }
		
		t++;									// increment time
	} // ------------- end community merging loop
	cout << "Q["<<t-1<<"] = "<<Q[t-1] << endl;
	
	// ----------------------------------------------------------------------
	// Record some results
	t1 = time(&t1);
	ofstream fout(ioparm.f_parm.c_str(), ios::app);
	fout << "---MODULARITY---\n";
	fout << "MAXQ------:\t" << Qmax.y  << "\n";
	fout << "STEP------:\t" << Qmax.x  << "\n";
	fout << "EXIT------:\t" << asctime(localtime(&t1));
	fout.close();

	cout << "exited safely" << endl;
	return 1;
}


// ------------------------------------------------------------------------------------ //
// ------------------------------------------------------------------------------------ //
// ------------------------------------------------------------------------------------ //
// FUNCTION DEFINITIONS --------------------------------------------------------------- //

void buildDeltaQMatrix() {
	
	// Given that we've now populated a sparse (unordered) adjacency matrix e (e), 
	// we now need to construct the intial dQ matrix according to the definition of dQ
	// which may be derived from the definition of modularity Q:
	//    Q(t) = \sum_{i} (e_{ii} - a_{i}^2) = Tr(e) - ||e^2||
	// thus dQ is
	//    dQ_{i,j} = 2* ( e_{i,j} - a_{i}a_{j} )
	//    where a_{i} = \sum_{j} e_{i,j} (i.e., the sum over the ith row)
	// To create dQ, we must insert each value of dQ_{i,j} into a binary search tree,
	// for the jth column. That is, dQ is simply an array of such binary search trees,
	// each of which represents the dQ_{x,j} adjacency vector. Having created dQ as
	// such, we may happily delete the matrix e in order to free up more memory.
	// The next step is to create a max-heap data structure, which contains the entries
	// of the following form (value, s, t), where the heap-key is 'value'. Accessing the
	// root of the heap gives us the next dQ value, and the indices (s,t) of the vectors
	// in dQ which need to be updated as a result of the merge.
	
	
	// First we compute e_{i,j}, and the compute+store the a_{i} values. These will be used
	// shortly when we compute each dQ_{i,j}.
	edge   *current;
	double  eij = (double)(0.5/gparm.m);				// intially each e_{i,j} = 1/m
	for (int i=1; i<gparm.maxid; i++) {				// for each row
		a[i] = 0.0;								// 
		if (e[i].so != 0) {							//    ensure it exists
			current = &e[i];						//    grab first edge
			a[i] = eij;							//    initialize a[i]
			while (current->next != NULL) {			//    loop through remaining edges
				a[i] += eij;						//       add another eij
				current = current->next;				//
			}
			Q[0] += -1.0*a[i]*a[i];					// calculate initial value of Q
		}
	}

	// now we create an empty (ordered) sparse matrix dq[]
	dq = new nodenub [gparm.maxid];						// initialize dq matrix
	for (int i=0; i<gparm.maxid; i++) {					// 
		dq[i].heap_ptr = NULL;							// no pointer in the heap at first
		if (e[i].so != 0) { dq[i].v = new vektor(2+(int)floor(gparm.m*a[i])); }
		else {			dq[i].v = NULL; }
	}
	h = new maxheap(gparm.n);						// allocate max-heap of size = number of nodes
	
	// Now we do all the work, which happens as we compute and insert each dQ_{i,j} into 
	// the corresponding (ordered) sparse vector dq[i]. While computing each dQ for a
	// row i, we track the maximum dQmax of the row and its (row,col) indices (i,j). Upon
	// finishing all dQ's for a row, we insert the tuple into the max-heap hQmax. That
	// insertion returns the itemaddress, which we then store in the nodenub heap_ptr for 
	// that row's vector.
	double    dQ;
	tuple	dQmax;										// for heaping the row maxes
	tuple*    itemaddress;									// stores address of item in maxheap

	for (int i=1; i<gparm.maxid; i++) {
		if (e[i].so != 0) {
			current = &e[i];								// grab first edge
			dQ      = 2.0*(eij-(a[current->so]*a[current->si]));   // compute its dQ
			dQmax.m = dQ;									// assume it is maximum so far
			dQmax.i = current->so;							// store its (row,col)
			dQmax.j = current->si;							// 
			dq[i].v->insertItem(current->si, dQ);				// insert its dQ
			while (current->next != NULL) {					// 
				current = current->next;						// step to next edge
				dQ = 2.0*(eij-(a[current->so]*a[current->si]));	// compute new dQ
				if (dQ > dQmax.m) {							// if dQ larger than current max
					dQmax.m = dQ;							//    replace it as maximum so far
					dQmax.j = current->si;					//    and store its (col)
				}
				dq[i].v->insertItem(current->si, dQ);			// insert it into vector[i]
			}
			dq[i].heap_ptr = h->insertItem(dQmax);				// store the pointer to its loc in heap
		}
	}

	delete [] elist;								// free-up adjacency matrix memory in two shots
	delete [] e;									// 
	return;
}

// ------------------------------------------------------------------------------------ //
// ------------------------------------------------------------------------------------ //
// ------------------------------------------------------------------------------------ //
// ------------------------------------------------------------------------------------ //
// ------------------------------------------------------------------------------------ //

void buildFilenames() {

	ioparm.f_input   = ioparm.d_in  + ioparm.filename;
	ioparm.f_parm    = ioparm.d_out + ioparm.s_scratch + "-fc_"  + ioparm.s_label + ".info";
	ioparm.f_joins   = ioparm.d_out + ioparm.s_scratch + "-fc_"  + ioparm.s_label + ".joins";
	ioparm.f_support = ioparm.d_out + ioparm.s_scratch + "-fc_"  + ioparm.s_label + ".supp";
	ioparm.f_net     = ioparm.d_out + ioparm.s_scratch + "-fc_"  + ioparm.s_label + ".wpairs";
	ioparm.f_group   = ioparm.d_out + ioparm.s_scratch + "-fc_"  + ioparm.s_label + ".groups";
	ioparm.f_gstats  = ioparm.d_out + ioparm.s_scratch + "-fc_"  + ioparm.s_label + ".hist";
	
	if (true) { ofstream flog(ioparm.f_parm.c_str(), ios::trunc); flog.close(); }
	time_t t; t = time(&t);
	ofstream flog(ioparm.f_parm.c_str(), ios::app);
	flog << "FASTCOMMUNITY_INFERENCE_ALGORITHM\n";
	flog << "START-----:\t" << asctime(localtime(&t));
	flog << "---FILES--------\n";
	flog << "DIRECTORY-:\t" << ioparm.d_out		<< "\n";
	flog << "F_IN------:\t" << ioparm.filename   << "\n";
	flog << "F_JOINS---:\t" << ioparm.s_scratch + "-fc_"  + ioparm.s_label + ".joins" << "\n";
	flog << "F_INFO----:\t" << ioparm.s_scratch + "-fc_"  + ioparm.s_label + ".info"  << "\n";
	if (ioparm.suppFlag) {
		flog << "F_SUPP----:\t" << ioparm.s_scratch + "-fc_"  + ioparm.s_label + ".supp" << "\n"; }
	if (ioparm.cutstep>0) {
		flog << "F_NET-----:\t" << ioparm.s_scratch + "-fc_"  + ioparm.s_label + ".wpairs" << "\n";
		flog << "F_GROUPS--:\t" << ioparm.s_scratch + "-fc_"  + ioparm.s_label + ".groups" << "\n";
		flog << "F_GDIST---:\t" << ioparm.s_scratch + "-fc_"  + ioparm.s_label + ".hist"   << "\n";
	}
	flog.close();
	
	return;
}

// ------------------------------------------------------------------------------------
// returns the support of the dQ[]

void dqSupport() {
	int    total = 0;
	int    count = 0;
	for (int i=0; i<gparm.maxid; i++) {
		if (dq[i].heap_ptr != NULL) { total += dq[i].v->returnNodecount(); count++; }
	}
	supportTot = total;
	supportAve = total/(double)count;
	return;
}

// ------------------------------------------------------------------------------------

void groupListsSetup() {
	
	list *newList;
	c = new stub [gparm.maxid];
	for (int i=0; i<gparm.maxid; i++) {
		if (e[i].so != 0) {								// note: internal indexing
			newList = new list;							// create new community member
			newList->index = i;							//    with index i
			c[i].members   = newList;					// point ith community at newList
			c[i].size		= 1;							// point ith community at newList
			c[i].last		= newList;					// point last[] at that element too
			c[i].valid	= true;						// mark as valid community
		}
	}
	
	return;
}

// ------------------------------------------------------------------------------------
// function for computing statistics on the list of groups

void groupListsStats() {

	gstats.numgroups = 0;
	gstats.maxsize   = 0;
	gstats.minsize   = gparm.maxid;
	double count     = 0.0;
	for (int i=0; i<gparm.maxid; i++) {
		if (c[i].valid) {
			gstats.numgroups++;							// count number of communities
			count += 1.0;
			if (c[i].size > gstats.maxsize) { gstats.maxsize = c[i].size; }  // find biggest community
			if (c[i].size < gstats.minsize) { gstats.minsize = c[i].size; }  // find smallest community
			// compute mean group size
			gstats.meansize = (double)(c[i].size)/count + (((double)(count-1.0)/count)*gstats.meansize);
		}
	}
	
	count = 0.0;
	gstats.sizehist = new double [gstats.maxsize+1];
	for (int i=0; i<gstats.maxsize+1; i++) { gstats.sizehist[i] = 0; }
	for (int i=0; i<gparm.maxid; i++) {
		if (c[i].valid) {
			gstats.sizehist[c[i].size] += 1.0;					// tabulate histogram of sizes
			count += 1.0;
		}
	}
	// convert histogram to pdf, and write it to disk
	for (int i=0; i<gstats.maxsize+1; i++) { gstats.sizehist[i] = gstats.sizehist[i]/count; }
	ofstream fgstat(ioparm.f_gstats.c_str(), ios::trunc);
	for (int i=gstats.minsize; i<gstats.maxsize+1; i++) {
		fgstat << i << "\t" << gstats.sizehist[i] << "\n";
	}
	fgstat.close();
	
	// record some statistics
	time_t t1; t1 = time(&t1);
	ofstream fout(ioparm.f_parm.c_str(), ios::app);
	fout << "---GROUPS-------\n";
	fout << "NUMGROUPS-:\t" << gstats.numgroups  << "\n";
	fout << "MINSIZE---:\t" << gstats.minsize    << "\n";
	fout << "MEANSIZE--:\t" << gstats.meansize   << "\n";
	fout << "MAXSIZE---:\t" << gstats.maxsize    << "\n";
	fout.close();
	
	return;