Allow float- and long double-precision compilation.

Makefile

@@ -1,7 +1,4 @@
 CXXFLAGS                 = -Wall
-#To enable parallel algorithms, use following (assuming gcc)
-#CXXFLAGS                 = -Wall -fopenmp
-
 # OUTPUT_OPTION            = -o .obj/$@
 swig_flags               = -Wall -c++ -python -outputtuple
 objects                  = assortativity.o \
@@ -26,7 +23,7 @@ objects                  = assortativity.o \
                            distance_bin.o \
                            distance_wei.o \
                            efficiency.o \
-			   eigenvector_centrality.o \
+                           eigenvector_centrality.o \
                            erange.o \
                            find_motif34.o \
                            findpaths.o \
@@ -50,10 +47,12 @@ objects                  = assortativity.o \
                            makeringlatticeCIJ.o \
                            maketoeplitzCIJ.o \
                            matching_ind.o \
-                           matlab/matlab.o \
-                           matlab/matlab_double.o \
-                           matlab/matlab_float.o \
-                           matlab/matlab_long_double.o \
+                           matlab/compare.o \
+                           matlab/convert.o \
+                           matlab/functions.o \
+                           matlab/index.o \
+                           matlab/operators.o \
+                           matlab/utility.o \
                            modularity_louvain.o \
                            modularity_newman.o \
                            module_degree_zscore.o \
@@ -82,9 +81,8 @@ objects                  = assortativity.o \
 include Makefile.vars
 
 # To make a multi-architecture library:
-# Compile once for each architecture
-# Add -arch to CXXFLAGS and LDFLAGS during compilation
-# lipo -create <architecture-specific libraries> -output <universal library>
+# Compile once for each architecture by adding -arch to CXXFLAGS and LDFLAGS during compilation
+# Then execute "lipo -create <architecture-specific libraries> -output <universal library>"
 
 .PHONY: all clean install swig uninstall
 
@@ -108,13 +106,17 @@ install: libbct.a
 		mkdir $(install_dir)/include/bct; \
 		mkdir $(install_dir)/include/bct/matlab; \
 	fi
-	cp matlab/macros.h $(install_dir)/include/bct/matlab
+	if [[ "$(CXXFLAGS)" == *GSL_FLOAT* ]]; then \
+		cat bct_float.h bct.h > _bct.h; \
+	elif [[ "$(CXXFLAGS)" == *GSL_DOUBLE* ]]; then \
+		cat bct_double.h bct.h > _bct.h; \
+	elif [[ "$(CXXFLAGS)" == *GSL_LONG_DOUBLE* ]]; then \
+		cat bct_long_double.h bct.h > _bct.h; \
+	fi
+	mv _bct.h $(install_dir)/include/bct/bct.h
 	cp matlab/matlab.h $(install_dir)/include/bct/matlab
-	cp matlab/matlab_double.h $(install_dir)/include/bct/matlab
-	cp matlab/matlab_float.h $(install_dir)/include/bct/matlab
-	cp matlab/matlab_long_double.h $(install_dir)/include/bct/matlab
 	cp matlab/sort.h $(install_dir)/include/bct/matlab
-	cp bct.h $(install_dir)/include/bct
+	cp precision.h $(install_dir)/include/bct
 	cp libbct.a $(install_dir)/lib
 
 uninstall:

Makefile.vars-template

@@ -1,31 +1,35 @@
 # Arguments to be sent to the C++ compiler
 # Some arguments may already be specified in Makefile
-# To enable parallel processing with OpenMP add "-fopenmp" (assuming gcc)
-# To build for a particular architecture add "-arch <arch>",
-# where <arch> is the specific architecture; e.g., "-arch i386"
-CXXFLAGS                +=
+# To compile with float precision, replace "-DGSL_DOUBLE" with "-DGSL_FLOAT"
+# To compile with long double precision, replace "-DGSL_DOUBLE" with "-DGSL_LONG_DOUBLE"
+# WARNING: SWIG wrappers are only set up to work with double precision
+# To enable parallel processing with OpenMP, add "-fopenmp"
+# To build for a particular architecture, add "-arch <arch>" (e.g., "-arch i386")
+# On Linux, "-fPIC" may be required?
+CXXFLAGS                += -DGSL_DOUBLE
 
 # Installation directory
 install_dir              = /usr/local
 
 # The following variables are only needed for SWIG
 # If you aren't generating Python bindings, you don't need to worry about them
 
-# A typical value for the Python header file directory is
-python_dir_apple         = /Library/Frameworks/Python.framework/Versions/Current/include/python2.6
+# Typical values for python_dir in different environments
+python_dir_apple         = /Library/Frameworks/Python.framework/Versions/2.6/include/python2.6
+python_dir_linux         = /usr/include/python2.6
 
 # Python header file directory
-# This points to the C header files required to compile the SWIG bindings
+# This is the location of the C++ header files required to build Python bindings
 # You may be able to use one of the previously defined variables
 # E.g., python_dir = $(python_dir_apple)
 python_dir               =
 
 # Typical values for swig_lib_flags in different environments
-# You probably don't need to change these
 swig_lib_flags_apple     = -bundle -flat_namespace -undefined suppress
 swig_lib_flags_linux     = -shared
 
 # Arguments for generating a shared library from SWIG wrappers
-# You can probably use one of the previously defined variables
+# You may be able to use one of the previously defined variables
 # E.g., swig_lib_flags = $(swig_lib_flags_apple)
+# On Linux, "-lpython2.6" may be required?
 swig_lib_flags           =

assortativity.cpp

@@ -1,59 +1,57 @@
 #include "bct.h"
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_vector.h>
 
-double assortativity(const gsl_vector*, const gsl_matrix*);
+FP_T assortativity(const VECTOR_T*, const MATRIX_T*);
 
 /*
  * Computes assortativity for a directed graph.  Connection weights are ignored.
  */
-double bct::assortativity_dir(const gsl_matrix* CIJ) {
+FP_T bct::assortativity_dir(const MATRIX_T* CIJ) {
 	if (safe_mode) check_status(CIJ, SQUARE | DIRECTED, "assortativity_dir");
 
 	// [id,od,deg] = degrees_dir(CIJ);
-	gsl_vector* deg = degrees_dir(CIJ);
+	VECTOR_T* deg = degrees_dir(CIJ);
 
 	// [i,j] = find(CIJ>0);
-	gsl_matrix* CIJ_gt_0 = compare_elements(CIJ, fp_greater, 0.0);
-	gsl_matrix* CIJ_gt_0_ij = find_ij(CIJ_gt_0);
-	gsl_matrix_free(CIJ_gt_0);
+	MATRIX_T* CIJ_gt_0 = compare_elements(CIJ, fp_greater, 0.0);
+	MATRIX_T* CIJ_gt_0_ij = find_ij(CIJ_gt_0);
+	MATRIX_ID(free)(CIJ_gt_0);
 
-	double ret = assortativity(deg, CIJ_gt_0_ij);
-	gsl_vector_free(deg);
-	gsl_matrix_free(CIJ_gt_0_ij);
+	FP_T ret = assortativity(deg, CIJ_gt_0_ij);
+	VECTOR_ID(free)(deg);
+	MATRIX_ID(free)(CIJ_gt_0_ij);
 	return ret;
 }
 
 /*
  * Computes assortativity for an undirected graph.  Connection weights are
  * ignored.
  */
-double bct::assortativity_und(const gsl_matrix* CIJ) {
+FP_T bct::assortativity_und(const MATRIX_T* CIJ) {
 	if (safe_mode) check_status(CIJ, SQUARE | UNDIRECTED, "assortativity_und");
 
 	// [deg] = degrees_und(m);
-	gsl_vector* deg = degrees_und(CIJ);
+	VECTOR_T* deg = degrees_und(CIJ);
 
 	// [i,j] = find(triu(CIJ,1)>0);
-	gsl_matrix* triu_CIJ = triu(CIJ, 1);
-	gsl_matrix* triu_CIJ_gt_0 = compare_elements(triu_CIJ, fp_greater, 0.0);
-	gsl_matrix_free(triu_CIJ);
-	gsl_matrix* triu_CIJ_gt_0_ij = find_ij(triu_CIJ_gt_0);
-	gsl_matrix_free(triu_CIJ_gt_0);
+	MATRIX_T* triu_CIJ = triu(CIJ, 1);
+	MATRIX_T* triu_CIJ_gt_0 = compare_elements(triu_CIJ, fp_greater, 0.0);
+	MATRIX_ID(free)(triu_CIJ);
+	MATRIX_T* triu_CIJ_gt_0_ij = find_ij(triu_CIJ_gt_0);
+	MATRIX_ID(free)(triu_CIJ_gt_0);
 
-	double ret = assortativity(deg, triu_CIJ_gt_0_ij);
-	gsl_vector_free(deg);
-	gsl_matrix_free(triu_CIJ_gt_0_ij);
+	FP_T ret = assortativity(deg, triu_CIJ_gt_0_ij);
+	VECTOR_ID(free)(deg);
+	MATRIX_ID(free)(triu_CIJ_gt_0_ij);
 	return ret;
 }
 
-double assortativity(const gsl_vector* deg, const gsl_matrix* ij) {
+FP_T assortativity(const VECTOR_T* deg, const MATRIX_T* ij) {
 	using namespace bct;
 
-	gsl_vector_const_view i = gsl_matrix_const_column(ij, 0);
-	gsl_vector_const_view j = gsl_matrix_const_column(ij, 1);
-	gsl_vector* degi = gsl_vector_alloc(ij->size1);
-	gsl_vector* degj = gsl_vector_alloc(ij->size1);		
+	VECTOR_ID(const_view) i = MATRIX_ID(const_column)(ij, 0);
+	VECTOR_ID(const_view) j = MATRIX_ID(const_column)(ij, 1);
+	VECTOR_T* degi = VECTOR_ID(alloc)(ij->size1);
+	VECTOR_T* degj = VECTOR_ID(alloc)(ij->size1);		
 
 	// K = length(i);
 	int K = length(&i.vector);
@@ -62,36 +60,36 @@ double assortativity(const gsl_vector* deg, const gsl_matrix* ij) {
 	for (int k = 0; k < K; k++) {
 
 		// degi(k) = deg(i(k));
-		int i_k = (int)gsl_vector_get(&i.vector, k);
-		gsl_vector_set(degi, k, gsl_vector_get(deg, i_k));
+		int i_k = (int)VECTOR_ID(get)(&i.vector, k);
+		VECTOR_ID(set)(degi, k, VECTOR_ID(get)(deg, i_k));
 
 		// degj(k) = deg(j(k));
-		int j_k = (int)gsl_vector_get(&j.vector, k);
-		gsl_vector_set(degj, k, gsl_vector_get(deg, j_k));
+		int j_k = (int)VECTOR_ID(get)(&j.vector, k);
+		VECTOR_ID(set)(degj, k, VECTOR_ID(get)(deg, j_k));
 	}
 
 	// r = (sum(degi.*degj)/K - (sum(0.5*(degi+degj))/K)^2)/(sum(0.5*(degi.^2+degj.^2))/K - (sum(0.5*(degi+degj))/K)^2);
 
-	gsl_vector* degi_mul_degj = copy(degi);
-	gsl_vector_mul(degi_mul_degj, degj);
-	double r1 = sum(degi_mul_degj) / (double)K;
-	gsl_vector_free(degi_mul_degj);
+	VECTOR_T* degi_mul_degj = copy(degi);
+	VECTOR_ID(mul)(degi_mul_degj, degj);
+	FP_T r1 = sum(degi_mul_degj) / (FP_T)K;
+	VECTOR_ID(free)(degi_mul_degj);
 
-	gsl_vector* degi_add_degj = copy(degi);
-	gsl_vector_add(degi_add_degj, degj);
-	double r2 = 0.5 * sum(degi_add_degj) / (double)K;
-	gsl_vector_free(degi_add_degj);
+	VECTOR_T* degi_add_degj = copy(degi);
+	VECTOR_ID(add)(degi_add_degj, degj);
+	FP_T r2 = 0.5 * sum(degi_add_degj) / (FP_T)K;
+	VECTOR_ID(free)(degi_add_degj);
 	r2 *= r2;
 
-	gsl_vector* degi_pow_2_add_degj_pow_2 = pow_elements(degi, 2);
-	gsl_vector* degj_pow_2 = pow_elements(degj, 2);
-	gsl_vector_add(degi_pow_2_add_degj_pow_2, degj_pow_2);
-	gsl_vector_free(degj_pow_2);
-	double r3 = 0.5 * sum(degi_pow_2_add_degj_pow_2) / (double)K;
-	gsl_vector_free(degi_pow_2_add_degj_pow_2);
+	VECTOR_T* degi_pow_2_add_degj_pow_2 = pow_elements(degi, 2);
+	VECTOR_T* degj_pow_2 = pow_elements(degj, 2);
+	VECTOR_ID(add)(degi_pow_2_add_degj_pow_2, degj_pow_2);
+	VECTOR_ID(free)(degj_pow_2);
+	FP_T r3 = 0.5 * sum(degi_pow_2_add_degj_pow_2) / (FP_T)K;
+	VECTOR_ID(free)(degi_pow_2_add_degj_pow_2);
 
-	gsl_vector_free(degi);
-	gsl_vector_free(degj);
+	VECTOR_ID(free)(degi);
+	VECTOR_ID(free)(degj);
 
 	return (r1 - r2) / (r3 - r2);
 }