Update to reflect some recent revisions to BCT MATLAB.

igormis · Jan 13, 2012 · 49b8e4d · 49b8e4d
1 parent 229ec44
commit 49b8e4d
Show file tree

Hide file tree

Showing 11 changed files with 96 additions and 42 deletions.
diff --git a/breadthdist.cpp b/breadthdist.cpp
@@ -28,7 +28,7 @@ MATRIX_T* BCT_NAMESPACE::breadthdist(const MATRIX_T* CIJ, MATRIX_T** D) {
 	logical_index_assign(_D, D_eq_0, GSL_POSINF);
 	MATRIX_ID(free)(D_eq_0);
 
-	// R = FP_T(D~=Inf);
+	// R = double(D~=Inf);
 	MATRIX_T* R = compare_elements(_D, fp_not_equal, GSL_POSINF);
 
 	if (D != NULL) *D = _D; else MATRIX_ID(free)(_D);

diff --git a/degrees_dir.cpp b/degrees_dir.cpp
@@ -7,7 +7,7 @@
 VECTOR_T* BCT_NAMESPACE::degrees_dir(const MATRIX_T* CIJ, VECTOR_T** id, VECTOR_T** od) {
 	if (safe_mode) check_status(CIJ, SQUARE | DIRECTED, "degrees_dir");
 
-	// CIJ = FP_T(CIJ~=0);
+	// CIJ = double(CIJ~=0);
 	// id = sum(CIJ,1);
 	VECTOR_T* _id = VECTOR_ID(alloc)(CIJ->size2);
 	for (int i = 0; i < (int)CIJ->size2; i++) {

diff --git a/degrees_und.cpp b/degrees_und.cpp
@@ -6,7 +6,7 @@
 VECTOR_T* BCT_NAMESPACE::degrees_und(const MATRIX_T* CIJ) {
 	if (safe_mode) check_status(CIJ, SQUARE | UNDIRECTED, "degrees_und");
 
-	// CIJ = FP_T(CIJ~=0);
+	// CIJ = double(CIJ~=0);
 	// deg = sum(CIJ);
 	VECTOR_T* deg = VECTOR_ID(alloc)(CIJ->size2);
 

diff --git a/efficiency.cpp b/efficiency.cpp
@@ -59,7 +59,7 @@ VECTOR_T* BCT_NAMESPACE::efficiency_local(const MATRIX_T* G) {
 				MATRIX_T* e = distance_inv(G_idx, NULL);
 				MATRIX_ID(free)(G_idx);
 
-				// E(u)=sum(e(:))./(k.^2-k);
+				// E(u)=sum(e(:))./(k^2-k);
 				VECTOR_T* e_v = to_vector(e);
 				MATRIX_ID(free)(e);
 				FP_T sum_e = sum(e_v);

diff --git a/findpaths.cpp b/findpaths.cpp
@@ -14,7 +14,7 @@
 std::vector<MATRIX_T*> BCT_NAMESPACE::findpaths(const MATRIX_T* CIJ, const VECTOR_T* sources, int qmax, VECTOR_T** plq, int* qstop, MATRIX_T** allpths, MATRIX_T** util) {
 	if (safe_mode) check_status(CIJ, SQUARE, "findpaths");
 
-	// CIJ = FP_T(CIJ~=0);
+	// CIJ = double(CIJ~=0);
 	MATRIX_T* _CIJ = compare_elements(CIJ, fp_not_equal, 0.0);
 
 	// N = size(CIJ,1);

diff --git a/findwalks.cpp b/findwalks.cpp
@@ -8,7 +8,7 @@
 std::vector<MATRIX_T*> BCT_NAMESPACE::findwalks(const MATRIX_T* CIJ, VECTOR_T** wlq) {
 	if (safe_mode) check_status(CIJ, SQUARE, "findwalks");
 
-	// CIJ = FP_T(CIJ~=0);
+	// CIJ = double(CIJ~=0);
 	MATRIX_T* _CIJ = compare_elements(CIJ, fp_not_equal, 0.0);
 
 	// N = size(CIJ,1);

diff --git a/jdegree.cpp b/jdegree.cpp
@@ -8,7 +8,7 @@
 MATRIX_T* BCT_NAMESPACE::jdegree(const MATRIX_T* CIJ) {
 	if (safe_mode) check_status(CIJ, SQUARE, "jdegree");
 
-	// CIJ = FP_T(CIJ~=0);
+	// CIJ = double(CIJ~=0);
 	MATRIX_T* _CIJ = compare_elements(CIJ, fp_not_equal, 0.0);
 
 	// N = size(CIJ,1);

diff --git a/matlab/matlab.h b/matlab/matlab.h
@@ -159,7 +159,7 @@ namespace MATLAB_NAMESPACE {
 	void ordinal_index_assign(VECTOR_T* v, const VECTOR_T* indices, FP_T value);
 	void ordinal_index_assign(VECTOR_T* v, const VECTOR_T* indices, const VECTOR_T* values);
 	VECTOR_T* logical_index(const VECTOR_T* v, const VECTOR_T* logical_v);
-	void logical_index_assign(VECTOR_T* v, const VECTOR_T* logical_v, FP_T values);
+	void logical_index_assign(VECTOR_T* v, const VECTOR_T* logical_v, FP_T value);
 	void logical_index_assign(VECTOR_T* v, const VECTOR_T* logical_v, const VECTOR_T* values);
 
 	// Matrix-by-integer indexing

diff --git a/modularity_newman.cpp b/modularity_newman.cpp
@@ -14,14 +14,23 @@ VECTOR_T* modularity(const MATRIX_T* B, int N, FP_T m);
 FP_T BCT_NAMESPACE::modularity_dir(const MATRIX_T* A, VECTOR_T** Ci) {
 	if (safe_mode) check_status(A, SQUARE | DIRECTED, "modularity_dir");
 
+	// N=length(A);
+	int N = length(A);
+
+	// n_perm = randperm(N);
+	gsl_permutation* n_perm = randperm(N);
+
+	// A = A(n_perm,n_perm);
+	VECTOR_T* n_perm_v to_vector(n_perm);
+	gsl_permutation_free(n_perm);
+	MATRIX_T* A_perm = MATRIX_ID(alloc)(N, N);
+	ordinal_index_assign(A_perm, n_perm_v, n_perm_v, A);
+
 	// Ki=sum(A,1);
-	VECTOR_T* Ki = sum(A, 1);
+	VECTOR_T* Ki = sum(A_perm, 1);
 
 	// Ko=sum(A,2);
-	VECTOR_T* Ko = sum(A, 2);
-
-	// N=length(A);
-	int N = length(A);
+	VECTOR_T* Ko = sum(A_perm, 2);
 
 	// m=sum(Ki);
 	FP_T m = sum(Ki);
@@ -36,7 +45,8 @@ FP_T BCT_NAMESPACE::modularity_dir(const MATRIX_T* A, VECTOR_T** Ci) {
 	MATRIX_ID(free)(Ki_m);
 	MATRIX_ID(transpose)(Ko_mul_Ki_transpose);
 	MATRIX_ID(scale)(Ko_mul_Ki_transpose, 1.0 / m);
-	MATRIX_T* b = copy(A);
+	MATRIX_T* b = copy(A_perm);
+	MATRIX_ID(free)(A_perm);
 	MATRIX_ID(sub)(b, Ko_mul_Ki_transpose);
 	MATRIX_ID(free)(Ko_mul_Ki_transpose);
 
@@ -72,6 +82,21 @@ FP_T BCT_NAMESPACE::modularity_dir(const MATRIX_T* A, VECTOR_T** Ci) {
 	FP_T Q = sum(sum_Q_m);
 	VECTOR_ID(free)(sum_Q_m);
 
+	// Ci_corrected = zeros(N,1);
+	VECTOR_T* _Ci_corrected = VECTOR_ID(alloc)(N);
+
+	// Ci_corrected(n_perm) = Ci;
+	for (int i = 0; i < N; i++) {
+		int index = (int)VECTOR_ID(get)(n_perm_v, i);
+		FP_T value = VECTOR_ID(get)(_Ci, i);
+		VECTOR_ID(set)(_Ci_corrected, index, value);
+	}
+	VECTOR_ID(free)(n_perm_v);
+
+	// Ci = Ci_corrected;
+	VECTOR_ID(free)(_Ci);
+	_Ci = _Ci_corrected;
+
 	if (Ci != NULL) *Ci = _Ci; else VECTOR_ID(free)(_Ci);
 	return Q;
 }
@@ -84,12 +109,21 @@ FP_T BCT_NAMESPACE::modularity_dir(const MATRIX_T* A, VECTOR_T** Ci) {
 FP_T BCT_NAMESPACE::modularity_und(const MATRIX_T* A, VECTOR_T** Ci) {
 	if (safe_mode) check_status(A, SQUARE | UNDIRECTED, "modularity_und");
 
-	// K=sum(A);
-	VECTOR_T* K = sum(A);
-
 	// N=length(A);
 	int N = length(A);
 
+	// n_perm = randperm(N);
+	gsl_permutation* n_perm = randperm(N);
+
+	// A = A(n_perm,n_perm);
+	VECTOR_T* n_perm_v to_vector(n_perm);
+	gsl_permutation_free(n_perm);
+	MATRIX_T* A_perm = MATRIX_ID(alloc)(N, N);
+	ordinal_index_assign(A_perm, n_perm_v, n_perm_v, A);
+
+	// K=sum(A);
+	VECTOR_T* K = sum(A_perm);
+
 	// m=sum(K);
 	FP_T m = sum(K);
 
@@ -101,7 +135,8 @@ FP_T BCT_NAMESPACE::modularity_und(const MATRIX_T* A, VECTOR_T** Ci) {
 	MATRIX_ID(free)(K_m_transpose);
 	MATRIX_ID(free)(K_m);
 	MATRIX_ID(scale)(K_m_transpose_mul_K_m, 1.0 / m);
-	MATRIX_T* B = copy(A);
+	MATRIX_T* B = copy(A_perm);
+	MATRIX_ID(free)(A_perm);
 	MATRIX_ID(sub)(B, K_m_transpose_mul_K_m);
 	MATRIX_ID(free)(K_m_transpose_mul_K_m);
 
@@ -130,6 +165,21 @@ FP_T BCT_NAMESPACE::modularity_und(const MATRIX_T* A, VECTOR_T** Ci) {
 	FP_T Q = sum(sum_Q_m);
 	VECTOR_ID(free)(sum_Q_m);
 
+	// Ci_corrected = zeros(N,1);
+	VECTOR_T* _Ci_corrected = VECTOR_ID(alloc)(N);
+
+	// Ci_corrected(n_perm) = Ci;
+	for (int i = 0; i < N; i++) {
+		int index = (int)VECTOR_ID(get)(n_perm_v, i);
+		FP_T value = VECTOR_ID(get)(_Ci, i);
+		VECTOR_ID(set)(_Ci_corrected, index, value);
+	}
+	VECTOR_ID(free)(n_perm_v);
+
+	// Ci = Ci_corrected;
+	VECTOR_ID(free)(_Ci);
+	_Ci = _Ci_corrected;
+
 	if (Ci != NULL) *Ci = _Ci; else VECTOR_ID(free)(_Ci);
 	return Q;
 }

diff --git a/participation_coef.cpp b/participation_coef.cpp
@@ -3,27 +3,23 @@
 #include "bct.h"
 
 /*
- * Computes nodal participation coefficient for a binary graph and its
- * corresponding community structure.  For a directed graph, computes "out-
- * neighbor" participation coefficient.
+ * Computes nodal participation coefficient and community structure.  For a
+ * directed graph, computes "out-neighbor" participation coefficient.
  */
-VECTOR_T* BCT_NAMESPACE::participation_coef(const MATRIX_T* A, const VECTOR_T* Ci) {
-	if (safe_mode) check_status(A, SQUARE | BINARY, "participation_coef");
+VECTOR_T* BCT_NAMESPACE::participation_coef(const MATRIX_T* W, const VECTOR_T* Ci) {
+	if (safe_mode) check_status(W, SQUARE | BINARY, "participation_coef");
 
-	// n=length(A);
-	int n = length(A);
+	// n=length(W);
+	int n = length(W);
 
-	// Ko=sum(A,2);
-	VECTOR_T* Ko = sum(A, 2);
+	// Ko=sum(W,2);
+	VECTOR_T* Ko = sum(W, 2);
 
-	// Ko(~Ko)=inf;
-	VECTOR_T* not_Ko = logical_not(Ko);
-	logical_index_assign(Ko, not_Ko, GSL_POSINF);
-	VECTOR_ID(free)(not_Ko);
-
-	// Gc=A*diag(Ci);
+	// Gc=(W~=0)*diag(Ci);
+	MATRIX_T* W_neq_0 = compare_elements(W, fp_not_equal, 0.0);
 	MATRIX_T* diag_Ci = diag(Ci);
-	MATRIX_T* Gc = mul(A, diag_Ci);
+	MATRIX_T* Gc = mul(W_neq_0, diag_Ci);
+	MATRIX_ID(free)(W_neq_0);
 	MATRIX_ID(free)(diag_Ci);
 
 	// Kc2=zeros(n,1);
@@ -32,26 +28,34 @@ VECTOR_T* BCT_NAMESPACE::participation_coef(const MATRIX_T* A, const VECTOR_T* C
 	// for i=1:max(Ci);
 	for (int i = 1; i <= (int)max(Ci); i++) {
 
-		// Kc2=Kc2+(sum(Gc==i,2).^2);
+		// Kc2=Kc2+(sum(W.*(Gc==i),2).^2);
 		MATRIX_T* Gc_eq_i = compare_elements(Gc, fp_equal, (FP_T)i);
-		VECTOR_T* sum_Gc_eq_i = sum(Gc_eq_i, 2);
+		MATRIX_T* W_mul_Gc_eq_i = copy(W);
+		MATRIX_ID(mul_elements)(W_mul_Gc_eq_i, Gc_eq_i);
 		MATRIX_ID(free)(Gc_eq_i);
-		VECTOR_T* sum_Gc_eq_i_pow_2 = pow_elements(sum_Gc_eq_i, 2.0);
-		VECTOR_ID(free)(sum_Gc_eq_i);
-		VECTOR_ID(add)(Kc2, sum_Gc_eq_i_pow_2);
-		VECTOR_ID(free)(sum_Gc_eq_i_pow_2);
+		VECTOR_T* sum_W_mul_Gc_eq_i = sum(W_mul_Gc_eq_i, 2);
+		MATRIX_ID(free)(W_mul_Gc_eq_i);
+		VECTOR_T* sum_W_mul_Gc_eq_i_pow_2 = pow_elements(sum_W_mul_Gc_eq_i, 2.0);
+		VECTOR_ID(free)(sum_W_mul_Gc_eq_i);
+		VECTOR_ID(add)(Kc2, sum_W_mul_Gc_eq_i_pow_2);
+		VECTOR_ID(free)(sum_W_mul_Gc_eq_i_pow_2);
 	}
 
 	MATRIX_ID(free)(Gc);
 
 	// P=ones(n,1)-Kc2./(Ko.^2);
 	VECTOR_T* Ko_pow_2 = pow_elements(Ko, 2.0);
-	VECTOR_ID(free)(Ko);
 	VECTOR_ID(div)(Kc2, Ko_pow_2);
 	VECTOR_ID(free)(Ko_pow_2);
 	VECTOR_T* P = ones_vector(n);
 	VECTOR_ID(sub)(P, Kc2);
 	VECTOR_ID(free)(Kc2);
 
+	// P(~Ko)=0;
+	VECTOR_T* not_Ko = logical_not(Ko);
+	VECTOR_ID(free)(Ko);
+	logical_index_assign(P, not_Ko, 0.0);
+	VECTOR_ID(free)(not_Ko);
+
 	return P;
 }
diff --git a/reachdist.cpp b/reachdist.cpp
@@ -105,7 +105,7 @@ void reachdist2(const MATRIX_T* CIJ, MATRIX_T** CIJpwr, MATRIX_T** R, MATRIX_T**
 	MATRIX_ID(free)(*CIJpwr);
 	*CIJpwr = temp;
 
-	// R = FP_T(R | ((CIJpwr)~=0));
+	// R = double(R | ((CIJpwr)~=0));
 	MATRIX_T* CIJpwr_neq_0 = compare_elements(*CIJpwr, fp_not_equal, 0.0);
 	MATRIX_T* R_or_CIJpwr_neq_0 = logical_or(*R, CIJpwr_neq_0);
 	MATRIX_ID(free)(CIJpwr_neq_0);