centrality.c

Go to the documentation of this file.

#include <stdio.h>
#include <stdlib.h>
#include <grass/gis.h>
#include <grass/Vect.h>
#include <grass/glocale.h>
#include <grass/dgl/graph.h>
#include <grass/neta.h>

void NetA_degree_centrality(dglGraph_s * graph, double *degree)
{
    int i;
    double nnodes = dglGet_NodeCount(graph);

    for (i = 1; i <= nnodes; i++)
        degree[i] =
            dglNodeGet_OutDegree(graph,
                                 dglGetNode(graph, (dglInt32_t) i)) / nnodes;
}

int NetA_eigenvector_centrality(dglGraph_s * graph, int iterations,
                                double error, double *eigenvector)
{
    int i, iter, nnodes;
    double *tmp;

    nnodes = dglGet_NodeCount(graph);
    tmp = (double *)G_calloc(nnodes + 1, sizeof(double));
    if (!tmp) {
        G_fatal_error(_("Out of memory"));
        return -1;
    }

    error *= error;
    for (i = 1; i <= nnodes; i++)
        eigenvector[i] = 1;
    for (iter = 0; iter < iterations; iter++) {
        for (i = 1; i <= nnodes; i++)
            tmp[i] = 0;
        dglInt32_t *node;
        dglNodeTraverser_s nt;
        dglEdgesetTraverser_s et;

        dglNode_T_Initialize(&nt, graph);
        for (node = dglNode_T_First(&nt); node; node = dglNode_T_Next(&nt)) {
            dglInt32_t node_id = dglNodeGet_Id(graph, node);
            double cur_value = eigenvector[node_id];
            dglInt32_t *edge;

            dglEdgeset_T_Initialize(&et, graph,
                                    dglNodeGet_OutEdgeset(graph, node));
            for (edge = dglEdgeset_T_First(&et); edge;
                 edge = dglEdgeset_T_Next(&et))
                tmp[dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge))] +=
                    cur_value * dglEdgeGet_Cost(graph, edge);

            dglEdgeset_T_Release(&et);
        }
        dglNode_T_Release(&nt);
        double cum_error = 0, max_value = tmp[1];

        for (i = 2; i <= nnodes; i++)
            if (tmp[i] > max_value)
                max_value = tmp[i];
        for (i = 1; i <= nnodes; i++) {
            tmp[i] /= max_value;
            cum_error +=
                (tmp[i] - eigenvector[i]) * (tmp[i] - eigenvector[i]);
            eigenvector[i] = tmp[i];
        }
        if (cum_error < error)
            break;

    }

    G_free(tmp);
    return 0;
}

int NetA_betweenness_closeness(dglGraph_s * graph, double *betweenness,
                               double *closeness)
{
    int i, j, nnodes, stack_size, count;
    dglInt32_t *dst, *node, *stack, *cnt, *delta;
    dglNodeTraverser_s nt;
    dglEdgesetTraverser_s et;
    dglHeap_s heap;
    struct ilist **prev;

    nnodes = dglGet_NodeCount(graph);

    dst = (dglInt32_t *) G_calloc(nnodes + 1, sizeof(dglInt32_t));
    prev = (struct ilist **)G_calloc(nnodes + 1, sizeof(struct ilist *));
    stack = (dglInt32_t *) G_calloc(nnodes, sizeof(dglInt32_t));
    cnt = (dglInt32_t *) G_calloc(nnodes + 1, sizeof(dglInt32_t));
    delta = (dglInt32_t *) G_calloc(nnodes + 1, sizeof(dglInt32_t));

    if (!dst || !prev || !stack || !cnt || !delta) {
        G_fatal_error(_("Out of memory"));
        return -1;
    }


    for (i = 1; i <= nnodes; i++) {
        prev[i] = Vect_new_list();
        if (closeness)
            closeness[i] = 0;
        if (betweenness)
            betweenness[i] = 0;
    }

    count = 0;
    G_percent_reset();
    dglNode_T_Initialize(&nt, graph);
    for (node = dglNode_T_First(&nt); node; node = dglNode_T_Next(&nt)) {
        G_percent(count++, nnodes, 1);
        dglInt32_t s = dglNodeGet_Id(graph, node);
        dglHeapData_u heap_data;
        dglHeapNode_s heap_node;

        stack_size = 0;
        for (i = 1; i <= nnodes; i++)
            Vect_reset_list(prev[i]);
        for (i = 1; i <= nnodes; i++) {
            cnt[i] = 0;
            dst[i] = -1;
        }
        dst[s] = 0;
        cnt[s] = 1;
        dglHeapInit(&heap);
        heap_data.ul = s;
        dglHeapInsertMin(&heap, 0, ' ', heap_data);
        while (1) {
            dglInt32_t v, dist;

            if (!dglHeapExtractMin(&heap, &heap_node))
                break;
            v = heap_node.value.ul;
            dist = heap_node.key;
            if (dst[v] < dist)
                continue;
            stack[stack_size++] = v;

            dglInt32_t *edge;

            dglEdgeset_T_Initialize(&et, graph,
                                    dglNodeGet_OutEdgeset(graph,
                                                          dglGetNode(graph,
                                                                     v)));
            for (edge = dglEdgeset_T_First(&et); edge;
                 edge = dglEdgeset_T_Next(&et)) {
                dglInt32_t *to = dglEdgeGet_Tail(graph, edge);
                dglInt32_t to_id = dglNodeGet_Id(graph, to);
                dglInt32_t d = dglEdgeGet_Cost(graph, edge);

                if (dst[to_id] == -1 || dst[to_id] > dist + d) {
                    dst[to_id] = dist + d;
                    Vect_reset_list(prev[to_id]);
                    heap_data.ul = to_id;
                    dglHeapInsertMin(&heap, dist + d, ' ', heap_data);
                }
                if (dst[to_id] == dist + d) {
                    cnt[to_id] += cnt[v];
                    Vect_list_append(prev[to_id], v);
                }
            }

            dglEdgeset_T_Release(&et);
        }
        dglHeapFree(&heap, NULL);
        for (i = 1; i <= nnodes; i++)
            delta[i] = 0;
        for (i = stack_size - 1; i >= 0; i--) {
            dglInt32_t w = stack[i];

            if (closeness)
                closeness[s] += dst[w];

            for (j = 0; j < prev[w]->n_values; j++) {
                dglInt32_t v = prev[w]->value[j];

                delta[v] += (cnt[v] / (double)cnt[w]) * (1.0 + delta[w]);
            }
            if (w != s && betweenness)
                betweenness[w] += delta[w];

        }
        if (closeness)
            closeness[s] /= (double)stack_size;
    }
    dglNode_T_Release(&nt);

    for (i = 1; i <= nnodes; i++)
        Vect_destroy_list(prev[i]);
    G_free(delta);
    G_free(cnt);
    G_free(stack);
    G_free(prev);
    G_free(dst);

    return 0;
};