bct-cpp/charpath.cpp at master · devuci/bct-cpp · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
#include<gsl/gsl_math.h>

#include"bct.h"

/*
 * WARNING: BCT_NAMESPACE::charpath_lambda takes a distance matrix, but
 * BCT_NAMESPACE::capped_charpath_lambda takes a connection matrix. Both should be
 * lengths, not weights (called distances in CalcMetric).
*/

/*
 * Given a distance matrix, computes characteristic path length.
*/
FP_TBCT_NAMESPACE::charpath_lambda(constMATRIX_T* D) {
if (safe_mode) check_status(D, SQUARE, "charpath_lambda");

// lambda = sum(sum(D(D~=Inf)))/length(nonzeros(D~=Inf));
MATRIX_T* D_neq_inf = compare_elements(D, fp_not_equal, GSL_POSINF);
VECTOR_T* D_idx = logical_index(D, D_neq_inf);
FP_T sum_D_idx = sum(D_idx);
VECTOR_ID(free)(D_idx);
FP_T ret = sum_D_idx / (FP_T)nnz(D_neq_inf);
MATRIX_ID(free)(D_neq_inf);
return ret;
}

/*
 * Given a connection matrix, computes capped characteristic path length.
*/
FP_TBCT_NAMESPACE::capped_charpath_lambda(constMATRIX_T* L) {
if (safe_mode) check_status(L, SQUARE, "capped_charpath_lambda");
int N = L->size1;
int nonzeros = 0;
FP_T lmean = 0.0;
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
if (i == j) {
continue;
 }
FP_T l = MATRIX_ID(get)(L, i, j);
if (fp_nonzero(l)) {
 nonzeros++;
 lmean += l;
 }
 }
	}
	lmean /= nonzeros;
MATRIX_T* D = distance_wei(L);
FP_T dmax = (FP_T)N * lmean;
FP_T dmean = 0.0;
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
if (i == j) {
continue;
 }
FP_T d = MATRIX_ID(get)(D, i, j);
 dmean += (d < dmax) ? d : dmax;
 }
	}
	dmean /= N * (N - 1);
MATRIX_ID(free)(D);
return dmean;
}

/*
 * Given a distance matrix, computes eccentricity, radius, and diameter.
*/
VECTOR_T* BCT_NAMESPACE::charpath_ecc(constMATRIX_T* D, FP_T* radius, FP_T* diameter) {
if (safe_mode) check_status(D, SQUARE, "charpath_ecc");

// ecc = max(D.*(D~=Inf),[],2);
MATRIX_T* D_finite = copy(D);
MATRIX_T* D_eq_inf = compare_elements(D, fp_equal, GSL_POSINF);
logical_index_assign(D_finite, D_eq_inf, 0.0);
MATRIX_ID(free)(D_eq_inf);
VECTOR_T* ecc = max(D_finite, 2);
MATRIX_ID(free)(D_finite);

// radius = min(ecc);
if (radius != NULL) {
 *radius = min(ecc);
	}

// diameter = max(ecc);
if (diameter != NULL) {
 *diameter = max(ecc);
	}

return ecc;
}