GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-6596d49e63
centrality.c
Go to the documentation of this file.
1 /*!
2  \file lib/vector/neta/centrality.c
3 
4  \brief Network Analysis library - centrality
5 
6  Centrality measures
7 
8  (C) 2009-2010 by Daniel Bundala, and the GRASS Development Team
9 
10  This program is free software under the GNU General Public License
11  (>=v2). Read the file COPYING that comes with GRASS for details.
12 
13  \author Daniel Bundala (Google Summer of Code 2009)
14  */
15 
16 #include <stdio.h>
17 #include <stdlib.h>
18 #include <grass/gis.h>
19 #include <grass/vector.h>
20 #include <grass/glocale.h>
21 #include <grass/dgl/graph.h>
22 #include <grass/neta.h>
23 
24 /*!
25  \brief Computes degree centrality measure.
26 
27  Array degree has to be properly initialised to nnodes+1 elements
28 
29  \param graph input graph
30  \param[out] degree array of degrees
31  */
32 void NetA_degree_centrality(dglGraph_s *graph, double *degree)
33 {
34  int i;
35  double nnodes = dglGet_NodeCount(graph);
36 
37  for (i = 1; i <= nnodes; i++)
38  degree[i] =
39  dglNodeGet_OutDegree(graph, dglGetNode(graph, (dglInt32_t)i)) /
40  nnodes;
41 }
42 
43 /*!
44  \brief Computes eigenvector centrality using edge costs as weights.
45 
46  \param graph input graph
47  \param iterations number of iterations
48  \param error ?
49  \param[out] eigenvector eigen vector value
50 
51  \return 0 on success
52  \return -1 on failure
53  */
54 int NetA_eigenvector_centrality(dglGraph_s *graph, int iterations, double error,
55  double *eigenvector)
56 {
57  int i, iter, nnodes;
58  double *tmp;
59 
60  nnodes = dglGet_NodeCount(graph);
61  tmp = (double *)G_calloc(nnodes + 1, sizeof(double));
62  if (!tmp) {
63  G_fatal_error(_("Out of memory"));
64  return -1;
65  }
66 
67  error *= error;
68  for (i = 1; i <= nnodes; i++)
69  eigenvector[i] = 1;
70  for (iter = 0; iter < iterations; iter++) {
71  for (i = 1; i <= nnodes; i++)
72  tmp[i] = 0;
73  dglInt32_t *node;
76 
77  dglNode_T_Initialize(&nt, graph);
78  for (node = dglNode_T_First(&nt); node; node = dglNode_T_Next(&nt)) {
79  dglInt32_t node_id = dglNodeGet_Id(graph, node);
80  double cur_value = eigenvector[node_id];
81  dglInt32_t *edge;
82 
83  dglEdgeset_T_Initialize(&et, graph,
84  dglNodeGet_OutEdgeset(graph, node));
85  for (edge = dglEdgeset_T_First(&et); edge;
86  edge = dglEdgeset_T_Next(&et))
87  tmp[dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge))] +=
88  cur_value * dglEdgeGet_Cost(graph, edge);
89 
91  }
92  dglNode_T_Release(&nt);
93  double cum_error = 0, max_value = tmp[1];
94 
95  for (i = 2; i <= nnodes; i++)
96  if (tmp[i] > max_value)
97  max_value = tmp[i];
98  for (i = 1; i <= nnodes; i++) {
99  tmp[i] /= max_value;
100  cum_error += (tmp[i] - eigenvector[i]) * (tmp[i] - eigenvector[i]);
101  eigenvector[i] = tmp[i];
102  }
103  if (cum_error < error)
104  break;
105  }
106 
107  G_free(tmp);
108  return 0;
109 }
110 
111 /*!
112  \brief Computes betweenness and closeness centrality measure using Brandes
113  algorithm.
114 
115  Edge costs must be nonnegative. If some edge costs are negative then
116  the behaviour of this method is undefined.
117 
118  \param graph input graph
119  \param[out] betweenness betweenness values
120  \param[out] closeness cloneness values
121 
122  \return 0 on success
123  \return -1 on failure
124  */
125 int NetA_betweenness_closeness(dglGraph_s *graph, double *betweenness,
126  double *closeness)
127 {
128  int i, j, nnodes, stack_size, count;
129  dglInt32_t *dst, *node, *stack, *cnt, *delta;
132  dglHeap_s heap;
133  struct ilist **prev;
134 
135  nnodes = dglGet_NodeCount(graph);
136 
137  dst = (dglInt32_t *)G_calloc(nnodes + 1, sizeof(dglInt32_t));
138  prev = (struct ilist **)G_calloc(nnodes + 1, sizeof(struct ilist *));
139  stack = (dglInt32_t *)G_calloc(nnodes, sizeof(dglInt32_t));
140  cnt = (dglInt32_t *)G_calloc(nnodes + 1, sizeof(dglInt32_t));
141  delta = (dglInt32_t *)G_calloc(nnodes + 1, sizeof(dglInt32_t));
142 
143  if (!dst || !prev || !stack || !cnt || !delta) {
144  G_fatal_error(_("Out of memory"));
145  return -1;
146  }
147 
148  for (i = 1; i <= nnodes; i++) {
149  prev[i] = Vect_new_list();
150  if (closeness)
151  closeness[i] = 0;
152  if (betweenness)
153  betweenness[i] = 0;
154  }
155 
156  count = 0;
157  G_percent_reset();
158  dglNode_T_Initialize(&nt, graph);
159  for (node = dglNode_T_First(&nt); node; node = dglNode_T_Next(&nt)) {
160  G_percent(count++, nnodes, 1);
161  dglInt32_t s = dglNodeGet_Id(graph, node);
162  dglHeapData_u heap_data;
163  dglHeapNode_s heap_node;
164 
165  stack_size = 0;
166  for (i = 1; i <= nnodes; i++)
167  Vect_reset_list(prev[i]);
168  for (i = 1; i <= nnodes; i++) {
169  cnt[i] = 0;
170  dst[i] = -1;
171  }
172  dst[s] = 0;
173  cnt[s] = 1;
174  dglHeapInit(&heap);
175  heap_data.ul = s;
176  dglHeapInsertMin(&heap, 0, ' ', heap_data);
177  while (1) {
178  dglInt32_t v, dist;
179 
180  if (!dglHeapExtractMin(&heap, &heap_node))
181  break;
182  v = heap_node.value.ul;
183  dist = heap_node.key;
184  if (dst[v] < dist)
185  continue;
186  stack[stack_size++] = v;
187 
188  dglInt32_t *edge;
189 
191  &et, graph, dglNodeGet_OutEdgeset(graph, dglGetNode(graph, v)));
192  for (edge = dglEdgeset_T_First(&et); edge;
193  edge = dglEdgeset_T_Next(&et)) {
194  dglInt32_t *to = dglEdgeGet_Tail(graph, edge);
195  dglInt32_t to_id = dglNodeGet_Id(graph, to);
196  dglInt32_t d = dglEdgeGet_Cost(graph, edge);
197 
198  if (dst[to_id] == -1 || dst[to_id] > dist + d) {
199  dst[to_id] = dist + d;
200  Vect_reset_list(prev[to_id]);
201  heap_data.ul = to_id;
202  dglHeapInsertMin(&heap, dist + d, ' ', heap_data);
203  }
204  if (dst[to_id] == dist + d) {
205  cnt[to_id] += cnt[v];
206  Vect_list_append(prev[to_id], v);
207  }
208  }
209 
211  }
212  dglHeapFree(&heap, NULL);
213  for (i = 1; i <= nnodes; i++)
214  delta[i] = 0;
215  for (i = stack_size - 1; i >= 0; i--) {
216  dglInt32_t w = stack[i];
217 
218  if (closeness)
219  closeness[s] += dst[w];
220 
221  for (j = 0; j < prev[w]->n_values; j++) {
222  dglInt32_t v = prev[w]->value[j];
223 
224  delta[v] += (cnt[v] / (double)cnt[w]) * (1.0 + delta[w]);
225  }
226  if (w != s && betweenness)
227  betweenness[w] += delta[w];
228  }
229  if (closeness)
230  closeness[s] /= (double)stack_size;
231  }
232  dglNode_T_Release(&nt);
233 
234  for (i = 1; i <= nnodes; i++)
235  Vect_destroy_list(prev[i]);
236  G_free(delta);
237  G_free(cnt);
238  G_free(stack);
239  G_free(prev);
240  G_free(dst);
241 
242  return 0;
243 }
#define NULL
Definition: ccmath.h:32
int NetA_eigenvector_centrality(dglGraph_s *graph, int iterations, double error, double *eigenvector)
Computes eigenvector centrality using edge costs as weights.
Definition: centrality.c:54
void NetA_degree_centrality(dglGraph_s *graph, double *degree)
Computes degree centrality measure.
Definition: centrality.c:32
int NetA_betweenness_closeness(dglGraph_s *graph, double *betweenness, double *closeness)
Computes betweenness and closeness centrality measure using Brandes algorithm.
Definition: centrality.c:125
void G_percent(long, long, int)
Print percent complete messages.
Definition: percent.c:61
void G_percent_reset(void)
Reset G_percent() to 0%; do not add newline.
Definition: percent.c:117
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:150
#define G_calloc(m, n)
Definition: defs/gis.h:95
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void Vect_destroy_list(struct ilist *)
Frees all memory associated with a struct ilist, including the struct itself.
int Vect_list_append(struct ilist *, int)
Append new item to the end of list if not yet present.
struct ilist * Vect_new_list(void)
Creates and initializes a struct ilist.
int Vect_reset_list(struct ilist *)
Reset ilist structure.
#define _(str)
Definition: glocale.h:10
void dglHeapInit(dglHeap_s *pheap)
Definition: heap.c:27
int dglHeapInsertMin(dglHeap_s *pheap, long key, unsigned char flags, dglHeapData_u value)
Definition: heap.c:50
void dglHeapFree(dglHeap_s *pheap, dglHeapCancelItem_fn pfnCancelItem)
Definition: heap.c:35
int dglHeapExtractMin(dglHeap_s *pheap, dglHeapNode_s *pnoderet)
Definition: heap.c:76
int count
long key
Definition: heap.h:35
dglHeapData_u value
Definition: heap.h:36
Definition: heap.h:41
List of integers.
Definition: gis.h:708
int n_values
Number of values in the list.
Definition: gis.h:716
int * value
Array of values.
Definition: gis.h:712
long dglInt32_t
Definition: type.h:36
unsigned long ul
Definition: heap.h:30
dglInt32_t * dglNode_T_First(dglNodeTraverser_s *pT)
dglInt32_t * dglNodeGet_OutEdgeset(dglGraph_s *pGraph, dglInt32_t *pnNode)
int dglEdgeset_T_Initialize(dglEdgesetTraverser_s *pT, dglGraph_s *pGraph, dglInt32_t *pnEdgeset)
int dglNode_T_Initialize(dglNodeTraverser_s *pT, dglGraph_s *pGraph)
dglInt32_t * dglEdgeset_T_Next(dglEdgesetTraverser_s *pT)
dglInt32_t * dglEdgeGet_Tail(dglGraph_s *pGraph, dglInt32_t *pnEdge)
int dglNodeGet_OutDegree(dglGraph_s *pGraph, dglInt32_t *pnNode)
dglInt32_t * dglNode_T_Next(dglNodeTraverser_s *pT)
int dglGet_NodeCount(dglGraph_s *pgraph)
void dglEdgeset_T_Release(dglEdgesetTraverser_s *pT UNUSED)
dglInt32_t dglEdgeGet_Cost(dglGraph_s *pGraph, dglInt32_t *pnEdge)
dglInt32_t * dglEdgeset_T_First(dglEdgesetTraverser_s *pT)
void dglNode_T_Release(dglNodeTraverser_s *pT)
dglInt32_t dglNodeGet_Id(dglGraph_s *pGraph, dglInt32_t *pnNode)
dglInt32_t * dglGetNode(dglGraph_s *pGraph, dglInt32_t nNodeId)