GRASS Programmer's Manual
6.4.2(2012)
|
00001 00016 #include <stdio.h> 00017 #include <stdlib.h> 00018 #include <grass/gis.h> 00019 #include <grass/Vect.h> 00020 #include <grass/glocale.h> 00021 #include <grass/dgl/graph.h> 00022 #include <grass/neta.h> 00023 00032 void NetA_degree_centrality(dglGraph_s * graph, double *degree) 00033 { 00034 int i; 00035 double nnodes = dglGet_NodeCount(graph); 00036 00037 for (i = 1; i <= nnodes; i++) 00038 degree[i] = 00039 dglNodeGet_OutDegree(graph, 00040 dglGetNode(graph, (dglInt32_t) i)) / nnodes; 00041 } 00042 00054 int NetA_eigenvector_centrality(dglGraph_s * graph, int iterations, 00055 double error, double *eigenvector) 00056 { 00057 int i, iter, nnodes; 00058 double *tmp; 00059 00060 nnodes = dglGet_NodeCount(graph); 00061 tmp = (double *)G_calloc(nnodes + 1, sizeof(double)); 00062 if (!tmp) { 00063 G_fatal_error(_("Out of memory")); 00064 return -1; 00065 } 00066 00067 error *= error; 00068 for (i = 1; i <= nnodes; i++) 00069 eigenvector[i] = 1; 00070 for (iter = 0; iter < iterations; iter++) { 00071 for (i = 1; i <= nnodes; i++) 00072 tmp[i] = 0; 00073 dglInt32_t *node; 00074 dglNodeTraverser_s nt; 00075 dglEdgesetTraverser_s et; 00076 00077 dglNode_T_Initialize(&nt, graph); 00078 for (node = dglNode_T_First(&nt); node; node = dglNode_T_Next(&nt)) { 00079 dglInt32_t node_id = dglNodeGet_Id(graph, node); 00080 double cur_value = eigenvector[node_id]; 00081 dglInt32_t *edge; 00082 00083 dglEdgeset_T_Initialize(&et, graph, 00084 dglNodeGet_OutEdgeset(graph, node)); 00085 for (edge = dglEdgeset_T_First(&et); edge; 00086 edge = dglEdgeset_T_Next(&et)) 00087 tmp[dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge))] += 00088 cur_value * dglEdgeGet_Cost(graph, edge); 00089 00090 dglEdgeset_T_Release(&et); 00091 } 00092 dglNode_T_Release(&nt); 00093 double cum_error = 0, max_value = tmp[1]; 00094 00095 for (i = 2; i <= nnodes; i++) 00096 if (tmp[i] > max_value) 00097 max_value = tmp[i]; 00098 for (i = 1; i <= nnodes; i++) { 00099 tmp[i] /= max_value; 00100 cum_error += 00101 (tmp[i] - eigenvector[i]) * (tmp[i] - eigenvector[i]); 00102 eigenvector[i] = tmp[i]; 00103 } 00104 if (cum_error < error) 00105 break; 00106 00107 } 00108 00109 G_free(tmp); 00110 return 0; 00111 } 00112 00126 int NetA_betweenness_closeness(dglGraph_s * graph, double *betweenness, 00127 double *closeness) 00128 { 00129 int i, j, nnodes, stack_size, count; 00130 dglInt32_t *dst, *node, *stack, *cnt, *delta; 00131 dglNodeTraverser_s nt; 00132 dglEdgesetTraverser_s et; 00133 dglHeap_s heap; 00134 struct ilist **prev; 00135 00136 nnodes = dglGet_NodeCount(graph); 00137 00138 dst = (dglInt32_t *) G_calloc(nnodes + 1, sizeof(dglInt32_t)); 00139 prev = (struct ilist **)G_calloc(nnodes + 1, sizeof(struct ilist *)); 00140 stack = (dglInt32_t *) G_calloc(nnodes, sizeof(dglInt32_t)); 00141 cnt = (dglInt32_t *) G_calloc(nnodes + 1, sizeof(dglInt32_t)); 00142 delta = (dglInt32_t *) G_calloc(nnodes + 1, sizeof(dglInt32_t)); 00143 00144 if (!dst || !prev || !stack || !cnt || !delta) { 00145 G_fatal_error(_("Out of memory")); 00146 return -1; 00147 } 00148 00149 00150 for (i = 1; i <= nnodes; i++) { 00151 prev[i] = Vect_new_list(); 00152 if (closeness) 00153 closeness[i] = 0; 00154 if (betweenness) 00155 betweenness[i] = 0; 00156 } 00157 00158 count = 0; 00159 G_percent_reset(); 00160 dglNode_T_Initialize(&nt, graph); 00161 for (node = dglNode_T_First(&nt); node; node = dglNode_T_Next(&nt)) { 00162 G_percent(count++, nnodes, 1); 00163 dglInt32_t s = dglNodeGet_Id(graph, node); 00164 dglHeapData_u heap_data; 00165 dglHeapNode_s heap_node; 00166 00167 stack_size = 0; 00168 for (i = 1; i <= nnodes; i++) 00169 Vect_reset_list(prev[i]); 00170 for (i = 1; i <= nnodes; i++) { 00171 cnt[i] = 0; 00172 dst[i] = -1; 00173 } 00174 dst[s] = 0; 00175 cnt[s] = 1; 00176 dglHeapInit(&heap); 00177 heap_data.ul = s; 00178 dglHeapInsertMin(&heap, 0, ' ', heap_data); 00179 while (1) { 00180 dglInt32_t v, dist; 00181 00182 if (!dglHeapExtractMin(&heap, &heap_node)) 00183 break; 00184 v = heap_node.value.ul; 00185 dist = heap_node.key; 00186 if (dst[v] < dist) 00187 continue; 00188 stack[stack_size++] = v; 00189 00190 dglInt32_t *edge; 00191 00192 dglEdgeset_T_Initialize(&et, graph, 00193 dglNodeGet_OutEdgeset(graph, 00194 dglGetNode(graph, 00195 v))); 00196 for (edge = dglEdgeset_T_First(&et); edge; 00197 edge = dglEdgeset_T_Next(&et)) { 00198 dglInt32_t *to = dglEdgeGet_Tail(graph, edge); 00199 dglInt32_t to_id = dglNodeGet_Id(graph, to); 00200 dglInt32_t d = dglEdgeGet_Cost(graph, edge); 00201 00202 if (dst[to_id] == -1 || dst[to_id] > dist + d) { 00203 dst[to_id] = dist + d; 00204 Vect_reset_list(prev[to_id]); 00205 heap_data.ul = to_id; 00206 dglHeapInsertMin(&heap, dist + d, ' ', heap_data); 00207 } 00208 if (dst[to_id] == dist + d) { 00209 cnt[to_id] += cnt[v]; 00210 Vect_list_append(prev[to_id], v); 00211 } 00212 } 00213 00214 dglEdgeset_T_Release(&et); 00215 } 00216 dglHeapFree(&heap, NULL); 00217 for (i = 1; i <= nnodes; i++) 00218 delta[i] = 0; 00219 for (i = stack_size - 1; i >= 0; i--) { 00220 dglInt32_t w = stack[i]; 00221 00222 if (closeness) 00223 closeness[s] += dst[w]; 00224 00225 for (j = 0; j < prev[w]->n_values; j++) { 00226 dglInt32_t v = prev[w]->value[j]; 00227 00228 delta[v] += (cnt[v] / (double)cnt[w]) * (1.0 + delta[w]); 00229 } 00230 if (w != s && betweenness) 00231 betweenness[w] += delta[w]; 00232 00233 } 00234 if (closeness) 00235 closeness[s] /= (double)stack_size; 00236 } 00237 dglNode_T_Release(&nt); 00238 00239 for (i = 1; i <= nnodes; i++) 00240 Vect_destroy_list(prev[i]); 00241 G_free(delta); 00242 G_free(cnt); 00243 G_free(stack); 00244 G_free(prev); 00245 G_free(dst); 00246 00247 return 0; 00248 };