gcfnd.c

用改进的数据结构快速计算复杂网络的社区划分

/*-
 * Copyright (c) 2008, Alexandre P. Francisco <[e-mail]>
 *
 * Permission to use, copy, modify, and/or distribute this software for any
 * purpose with or without fee is hereby granted, provided that the above
 * copyright notice and this permission notice appear in all copies.
 *
 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
 */

#include <limits.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#ifdef _RAND_CMP
#include <sys/time.h>
#endif /* _RAND_CMP */

#ifdef _TIMER
#include <time.h>
#endif /* _TIMER */

#include <assert.h>

#include "gcfnd.h"
#include "unfnd.h"

#ifndef _BOOST_RELAXED_HEAP
#include "pqueue.h"
#else
#include "pqueue_boost.h"
#endif

/*
 * Adjacency list node definition. This data structure is specially
 * designed to provide cross reference between adjacency lists with
 * an access time of O(1).
 */
struct adj_node {
	int 	id;
	int	u;
	int	v;
	int	state;
	struct adj_node	*u_nxt;
	struct adj_node *u_prv;
	struct adj_node *v_prv;
	struct adj_node *v_nxt;
};

/*
 * Edge access macros.
 */
#define e_nxt(t,x) (((x)->u == (t)) ? (x)->u_nxt : (x)->v_nxt)
#define e_adj(t,x) (((x)->u == (t)) ? (x)->v : (x)->u)

/*
 * Given the adjacency list for each vertex and the list of edges, the
 * edge 'e' is added to the adjacency list of vertex 'x', if possible.
 */
static void
add_edge(struct adj_node **v_adj, struct adj_node *e_lst, int e, int x)
{
	struct adj_node *edge = NULL;
	
	edge = e_lst + e;
	edge->state = 1;
	if (edge->u == INT_MAX) {
		edge->u = x;
		edge->u_prv = NULL;
		edge->u_nxt = v_adj[x];
		if (v_adj[x] != NULL) {
			if (v_adj[x]->u == x)
				v_adj[x]->u_prv = edge;
			else
				v_adj[x]->v_prv = edge;
		}
		v_adj[x] = edge;
	} else if (edge->v == INT_MAX) {
		edge->v = x;
		edge->v_prv = NULL;
		edge->v_nxt = v_adj[x];
		if (v_adj[x] != NULL) {
			if (v_adj[x]->u == x)
				v_adj[x]->u_prv = edge;
			else
				v_adj[x]->v_prv = edge;
		}
		v_adj[x] = edge;
	} else
		abort();
}

/*
 * Given the adjacency list for each vertex and the list of edges, the
 * edge 'e' is removed from the adjacency list of vertex 'x'.
 */
static void
remove_edge(struct adj_node **v_adj, struct adj_node *e_lst, int e, int x)
{
	struct adj_node *edge = NULL;
	
	edge = e_lst + e;
	edge->state = 0;
	if (edge->u == x) {
		if (edge->u_nxt != NULL) {
			if (edge->u_nxt->u == x)
				edge->u_nxt->u_prv = edge->u_prv;
			else
				edge->u_nxt->v_prv = edge->u_prv;
		}
		if (edge->u_prv != NULL) {
			if (edge->u_prv->u == x)
				edge->u_prv->u_nxt = edge->u_nxt;
			else
				edge->u_prv->v_nxt = edge->u_nxt;
		} else
			v_adj[x] = edge->u_nxt;
		edge->u_nxt = edge->u_prv = NULL;
		edge->u = INT_MAX;
	} else if (edge->v == x) {
		if (edge->v_nxt != NULL) {
			if (edge->v_nxt->u == x)
				edge->v_nxt->u_prv = edge->v_prv;
			else
				edge->v_nxt->v_prv = edge->v_prv;
		}
		if (edge->v_prv != NULL) {
			if (edge->v_prv->u == x)
				edge->v_prv->u_nxt = edge->v_nxt;
			else
				edge->v_prv->v_nxt = edge->v_nxt;
		} else
			v_adj[x] = edge->v_nxt;
		edge->v_nxt = edge->v_prv = NULL;
		edge->v = INT_MAX;
	} else
		abort();
}

/* Comparison function. */
int
_dQcmp(double a, double b)
{
	if (a > b)
		return 1;

	if (a < b)
		return -1;

#ifdef _RAND_CMP
#define rand_p() (((double) rand()) / RAND_MAX)

	return (rand_p() > 0.5) ? 1 : -1;
#else
	return 0;
#endif /* _RAND_CMP */
}


/*
 * An implementation with ONE heap and cross-linked adjacency lists. At
 * least 2 times faster when compared with the original implementation.
 */
void
gcfnd(int nv, edge *edge_lst, int ne, int *com_id, int *com_nb, double *cQ)
{
	int u, v, x, e_id, i, *mtmp;
	size_t iter, *mask, *v_dg;
	double Q, maxQ, *dQ, *A;
	struct adj_node *e, *ex, *e_lst, **v_adj;
	pqueue *pq;

#ifndef _NO_CID
	disjoint_set ds;
#endif /* NO_CID */

#ifdef _PROGRESS_METER
	int mcount;
#endif /* _PROGRESS-METER */

#ifdef _TIMER
	clock_t e_time;
#endif /* _TIMER */

#ifdef _RAND_CMP
	struct timeval t;
	
	gettimeofday(&t, NULL);
        srand((unsigned) (t.tv_sec ^ t.tv_usec));
#endif /* _RAND_CMP */

	/* Memory allocation: O(V + E) space. */
	/* Adjacency for each vertex, i.e., adjacency list for short. */
	v_adj = (struct adj_node **) malloc(sizeof(struct adj_node *)*nv);
	v_dg = (size_t *) malloc(sizeof(size_t)*nv);

	/* Edges allocation. */
	e_lst = (struct adj_node *) malloc(sizeof(struct adj_node)*ne);

	/* Community find stuff. */
	dQ = (double *) malloc(sizeof(double)*ne);
	A = (double *) malloc(sizeof(double)*nv);

	/* Priority queue stuff. */
	pq = pqueue_new(ne, dQ, _dQcmp);

#ifndef _NO_CID
	/* Disjoint set. */
	ds = make_set(nv);
#endif /* _NO_CID */

	/* Auxiliary stuff. */
	mask = (size_t *) malloc(sizeof(size_t)*nv);
	mtmp = (int *) malloc(sizeof(int)*nv);

	/* Graph loading and initialization: O(V + E) time. */
	for (i = 0; i < nv; i++) {
		v_adj[i] = NULL;
		v_dg[i] = 0;
		A[i] = 0.0;
		
		mask[i] = 0;
	}

	for (i = 0; i < ne; i++) {
		u = edge_lst[i].u;
		v = edge_lst[i].v;

		e_lst[i].id = i;
		e_lst[i].state = 1;
		e_lst[i].u = u;
		e_lst[i].v = v;
		A[u] += 1.0 / ((double) (2 * ne));
		A[v] += 1.0 / ((double) (2 * ne));

		/* Add edge 'i' to the head of the adjacency list of 'u'. */
		e_lst[i].u_nxt = v_adj[u];
		e_lst[i].u_prv = NULL;
 		if (v_adj[u] != NULL) {
			if (v_adj[u]->u == u)
				v_adj[u]->u_prv = e_lst + i;
			else
				v_adj[u]->v_prv = e_lst + i;
		}
		v_adj[u] = e_lst + i;
		v_dg[u] ++;

		/* Add edge 'i' to the head of the adjacency list of 'v'. */
		e_lst[i].v_nxt = v_adj[v];
		e_lst[i].v_prv = NULL;
		if (v_adj[v] != NULL) {
			if (v_adj[v]->v == v)
				v_adj[v]->v_prv = e_lst + i;
			else
				v_adj[v]->u_prv = e_lst + i;
		}
		v_adj[v] = e_lst + i;
		v_dg[v] ++;
	}

#ifdef _TIMER
	e_time = clock();
#endif /* _TIMER */

	/* 'Q' and 'dQ' initialization: O(V + E) time. */
	maxQ = Q = 0.0;
	for (i = 0; i < nv; i++)
		Q -= A[i]*A[i];

	for (i = 0; i < ne; i++) {
		dQ[i] = 2 * (1.0 / ((double) (2 * ne))
		    - A[e_lst[i].u] * A[e_lst[i].v]);
		pqueue_push(pq, i);
	}

#ifdef _PROGRESS_METER
	mcount = 0;
#endif /* _PROGRESS_METER */

	/* Find communities. */
	iter = 0;
	while (! pqueue_empty(pq)) {
		iter++;

		/* Extract the edge with maximum 'dQ' from queue. */
		e_id = pqueue_top(pq);
		pqueue_pop(pq);
		
		/* Check if Q can increase. */
		if (dQ[e_id] < 0.0)
			break; /* Best community structure found. */

		if(e_lst[e_id].state == 0) /* A residual edge. */
			continue;

		e_lst[e_id].state = 0;

		/* Update modularity 'Q'. */
		Q += dQ[e_id];
		u = e_lst[e_id].u;
		v = e_lst[e_id].v;

		if (v_dg[u] > v_dg[v]) {
			x = u;
			u = v;
			v = x;
		}
		
		/* Remove edge from the adjacency lists of 'u' and 'v'. */
		remove_edge(v_adj, e_lst, e_id, u);
		remove_edge(v_adj, e_lst, e_id, v);
		v_dg[u] --;
		v_dg[v] --;

		/* Preprocess 'u' adjacency. */
		for (e = v_adj[u]; e != NULL; e = e_nxt(u,e)) {
			mask[e_adj(u,e)] = iter;
			mtmp[e_adj(u,e)] = e->id;
		}

		/* Add the adjacency of 'u' to the adjacency of 'v'. */
		for (e = v_adj[v]; e != NULL; e = e_nxt(v,e)) {
			if (mask[e_adj(v,e)] == iter) {
				/* 'k' is connected to 'u' and 'v'. */
				dQ[e->id] += dQ[mtmp[e_adj(v,e)]];
				x = e_adj(v,e);
				remove_edge(v_adj, e_lst, mtmp[x], x);
				remove_edge(v_adj, e_lst, mtmp[x], u);
				v_dg[x] --;
				v_dg[u] --;
			} else {
				/* 'k' is connected to 'v' but not to 'u'. */
				dQ[e->id] -= 2 * A[u] * A[e_adj(v,e)];
			}
			pqueue_update(pq, e->id);
			mask[e_adj(v,e)] = 0;
		}

		for (e = v_adj[u]; e != NULL;) {
			if(mask[e_adj(u,e)] != iter) {
				e = e_nxt(u,e);
				continue;
			}
			/* Next edge must be known, adjacency will change. */
			ex = e_nxt(u,e);

			/* Move 'e' from 'u' to 'v'. */
			x = e_adj(u,e);
			remove_edge(v_adj, e_lst, e->id, u);
			v_dg[u] --;
			add_edge(v_adj, e_lst, e->id, v);
			v_dg[v] ++;

			/* 'k' is connected to 'u' but not to 'v'. */
			dQ[e->id] -= 2 * A[v] * A[x];
			pqueue_update(pq, e->id);
			e = ex;
		}

		/* Update dendrogram and 'A'. */
#ifndef _NO_CID
		union_set(ds, u, v);
#endif /* _NO_CID */
		A[v] += A[u];

		if (maxQ < Q)
			maxQ = Q;

#ifdef _PROGRESS_METER
		if ((mcount++)%100 == 0)
			fprintf(stderr, "\rcfind: %6.2f%% (%6.5f)",
			    ((double) mcount)/(nv-1) * 100, maxQ);
#endif /* _PROGESS_METER */

	}

#ifdef _PROGRESS_METER
	fprintf(stderr, "\rcfind: 100.00%% (%6.5f)\n", maxQ);
#endif /* _PROGESS_METER */

#ifdef _TIMER
	fprintf(stderr, "%fs seconds\n",
	    (((double) clock()) - ((double) e_time)) / CLOCKS_PER_SEC);
#endif /* _TIMER */

#ifndef _NO_CID
	/* Process communities. */
	*com_nb = 0;
	memset(com_id, 0xFF, sizeof(int)*nv);
	for(i = 0; i < nv; i++){
		x = find_set(ds, i); 
		if(com_id[x] == -1)
			com_id[x] = (*com_nb) ++;
		com_id[i] = com_id[x];
	}
#endif /* _NO_CID */
	*cQ = maxQ;

	/* Memory deallocation. */
	free(v_dg);
	free(v_adj);
	free(e_lst);
	pqueue_free(pq);
#ifndef _NO_CID
	free(ds);
#endif /* _NO_CID */
	free(A);
	free(mask);
	free(mtmp);
	free(dQ);
}