Geografic-Information-System/vector/v.net.centrality/main.c

/****************************************************************
 *
 * MODULE:     v.net.centrality
 *
 * AUTHOR(S):  Daniel Bundala
 *
 * PURPOSE:    This module computes various centrality measures
 *
 * COPYRIGHT:  (C) 2002-2005 by the GRASS Development Team
 *
 *             This program is free software under the
 *             GNU General Public License (>=v2).
 *             Read the file COPYING that comes with GRASS
 *             for details.
 *
 ****************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <grass/gis.h>
#include <grass/vector.h>
#include <grass/glocale.h>
#include <grass/dbmi.h>
#include <grass/neta.h>

/*Global variables */
struct Option *deg_opt, *close_opt, *betw_opt, *eigen_opt;
double *deg, *closeness, *betw, *eigen;

/* Attribute table */
dbString sql, tmp;
dbDriver *driver;
struct field_info *Fi;

void append_string(dbString * string, char *s)
{
    db_append_string(string, ", ");
    db_append_string(string, s);
    db_append_string(string, " double precision");
}

void append_double(dbString * string, double d)
{
    char buf[50];

    sprintf(buf, ",%f", d);
    db_append_string(string, buf);
}

void process_node(int node, int cat)
{
    char buf[2000];

    sprintf(buf, "INSERT INTO %s VALUES(%d", Fi->table, cat);
    db_set_string(&sql, buf);
    if (deg_opt->answer)
	append_double(&sql, deg[node]);
    if (close_opt->answer)
	append_double(&sql, closeness[node]);
    if (betw_opt->answer)
	append_double(&sql, betw[node]);
    if (eigen_opt->answer)
	append_double(&sql, eigen[node]);

    db_append_string(&sql, ")");
    if (db_execute_immediate(driver, &sql) != DB_OK) {
	db_close_database_shutdown_driver(driver);
	G_fatal_error(_("Cannot insert new record: %s"), db_get_string(&sql));
    }
}

int main(int argc, char *argv[])
{
    struct Map_info In, Out;
    static struct line_pnts *Points;
    struct line_cats *Cats;
    struct GModule *module;	/* GRASS module for parsing arguments */
    struct Option *map_in, *map_out;
    struct Option *cat_opt, *where_opt, *afield_opt, *nfield_opt, *abcol,
                  *afcol, *ncol;
    struct Option *iter_opt, *error_opt;
    struct Flag *geo_f, *add_f;
    int chcat, with_z;
    int afield, nfield, mask_type;
    struct varray *varray;
    dglGraph_s *graph;
    int i, geo, nnodes, nlines, j, max_cat;
    char buf[2000], *covered;

    /* initialize GIS environment */
    G_gisinit(argv[0]);		/* reads grass env, stores program name to G_program_name() */

    /* initialize module */
    module = G_define_module();
    G_add_keyword(_("vector"));
    G_add_keyword(_("network"));
    G_add_keyword(_("centrality measures"));
    module->description =
	_("Computes degree, centrality, betweeness, closeness and eigenvector "
	 "centrality measures in the network.");

    /* Define the different options as defined in gis.h */
    map_in = G_define_standard_option(G_OPT_V_INPUT);

    afield_opt = G_define_standard_option(G_OPT_V_FIELD);
    afield_opt->key = "arc_layer";
    afield_opt->answer = "1";
    afield_opt->label = _("Arc layer");
    afield_opt->guisection = _("Cost");

    nfield_opt = G_define_standard_option(G_OPT_V_FIELD);
    nfield_opt->key = "node_layer";
    nfield_opt->answer = "2";
    nfield_opt->label = _("Node layer");
    nfield_opt->guisection = _("Cost");

    map_out = G_define_standard_option(G_OPT_V_OUTPUT);

    cat_opt = G_define_standard_option(G_OPT_V_CATS);
    cat_opt->guisection = _("Selection");
    where_opt = G_define_standard_option(G_OPT_DB_WHERE);
    where_opt->guisection = _("Selection");

    afcol = G_define_standard_option(G_OPT_DB_COLUMN);
    afcol->key = "arc_column";
    afcol->required = NO;
    afcol->description =
	_("Arc forward/both direction(s) cost column (number)");
    afcol->guisection = _("Cost");

    abcol = G_define_standard_option(G_OPT_DB_COLUMN);
    abcol->key = "arc_backward_column";
    abcol->required = NO;
    abcol->description = _("Arc backward direction cost column (number)");
    abcol->guisection = _("Cost");

    ncol = G_define_option();
    ncol->key = "node_column";
    ncol->type = TYPE_STRING;
    ncol->required = NO;
    ncol->description = _("Node cost column (number)");
    ncol->guisection = _("Cost");

    deg_opt = G_define_standard_option(G_OPT_DB_COLUMN);
    deg_opt->key = "degree";
    deg_opt->required = NO;
    deg_opt->description = _("Name of degree centrality column");
    deg_opt->guisection = _("Columns");

    close_opt = G_define_standard_option(G_OPT_DB_COLUMN);
    close_opt->key = "closeness";
    close_opt->required = NO;
    close_opt->description = _("Name of closeness centrality column");
    close_opt->guisection = _("Columns");

    betw_opt = G_define_standard_option(G_OPT_DB_COLUMN);
    betw_opt->key = "betweenness";
    betw_opt->required = NO;
    betw_opt->description = _("Name of betweenness centrality column");
    betw_opt->guisection = _("Columns");

    eigen_opt = G_define_standard_option(G_OPT_DB_COLUMN);
    eigen_opt->key = "eigenvector";
    eigen_opt->required = NO;
    eigen_opt->description = _("Name of eigenvector centrality column");
    eigen_opt->guisection = _("Columns");

    iter_opt = G_define_option();
    iter_opt->key = "iterations";
    iter_opt->answer = "1000";
    iter_opt->type = TYPE_INTEGER;
    iter_opt->required = NO;
    iter_opt->description =
	_("Maximum number of iterations to compute eigenvector centrality");

    error_opt = G_define_option();
    error_opt->key = "error";
    error_opt->answer = "0.1";
    error_opt->type = TYPE_DOUBLE;
    error_opt->required = NO;
    error_opt->description =
	_("Cumulative error tolerance for eigenvector centrality");

    geo_f = G_define_flag();
    geo_f->key = 'g';
    geo_f->description =
	_("Use geodesic calculation for longitude-latitude locations");

    add_f = G_define_flag();
    add_f->key = 'a';
    add_f->description = _("Add points on nodes");

    /* options and flags parser */
    if (G_parser(argc, argv))
	exit(EXIT_FAILURE);
    /* TODO: make an option for this */
    mask_type = GV_LINE | GV_BOUNDARY;

    Points = Vect_new_line_struct();
    Cats = Vect_new_cats_struct();

    Vect_check_input_output_name(map_in->answer, map_out->answer,
				 G_FATAL_EXIT);

    Vect_set_open_level(2);

    if (1 > Vect_open_old(&In, map_in->answer, ""))
	G_fatal_error(_("Unable to open vector map <%s>"), map_in->answer);

    with_z = Vect_is_3d(&In);

    if (0 > Vect_open_new(&Out, map_out->answer, with_z)) {
	Vect_close(&In);
	G_fatal_error(_("Unable to create vector map <%s>"), map_out->answer);
    }


    if (geo_f->answer) {
	geo = 1;
	if (G_projection() != PROJECTION_LL)
	    G_warning(_("The current projection is not longitude-latitude"));
    }
    else
	geo = 0;

    /* parse filter option and select appropriate lines */
    afield = Vect_get_field_number(&In, afield_opt->answer);
    nfield = Vect_get_field_number(&In, nfield_opt->answer);

    if (where_opt->answer || cat_opt->answer) {
	chcat = (NetA_initialise_varray(&In, afield, GV_POINT,
	                                where_opt->answer,
					cat_opt->answer, &varray) > 0);
    }
    else
	chcat = 0;

    /* Create table */
    Fi = Vect_default_field_info(&Out, 1, NULL, GV_1TABLE);
    Vect_map_add_dblink(&Out, 1, NULL, Fi->table, GV_KEY_COLUMN, Fi->database,
			Fi->driver);
    db_init_string(&sql);
    driver = db_start_driver_open_database(Fi->driver, Fi->database);
    if (driver == NULL)
	G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
		      Fi->database, Fi->driver);
    db_set_error_handler_driver(driver);

    db_init_string(&tmp);
    if (deg_opt->answer)
	append_string(&tmp, deg_opt->answer);
    if (close_opt->answer)
	append_string(&tmp, close_opt->answer);
    if (betw_opt->answer)
	append_string(&tmp, betw_opt->answer);
    if (eigen_opt->answer)
	append_string(&tmp, eigen_opt->answer);
    sprintf(buf,
	    "create table %s(cat integer%s)", Fi->table, db_get_string(&tmp));

    db_set_string(&sql, buf);
    G_debug(2, "%s", db_get_string(&sql));

    if (db_execute_immediate(driver, &sql) != DB_OK) {
	G_fatal_error(_("Unable to create table: '%s'"), db_get_string(&sql));
    }

    if (db_create_index2(driver, Fi->table, GV_KEY_COLUMN) != DB_OK)
	G_warning(_("Cannot create index"));

    if (db_grant_on_table
	(driver, Fi->table, DB_PRIV_SELECT, DB_GROUP | DB_PUBLIC) != DB_OK)
	G_fatal_error(_("Cannot grant privileges on table <%s>"), Fi->table);

    db_begin_transaction(driver);

    Vect_copy_head_data(&In, &Out);
    Vect_hist_copy(&In, &Out);
    Vect_hist_command(&Out);

    if (0 != Vect_net_build_graph(&In, mask_type, afield, nfield, afcol->answer,
                                  abcol->answer, ncol->answer, geo, 0))
        G_fatal_error(_("Unable to build graph for vector map <%s>"), Vect_get_full_name(&In));
    
    graph = Vect_net_get_graph(&In);
    nnodes = dglGet_NodeCount(graph);

    deg = closeness = betw = eigen = NULL;

    covered = (char *)G_calloc(nnodes + 1, sizeof(char));
    if (!covered)
	G_fatal_error(_("Out of memory"));

    if (deg_opt->answer) {
	deg = (double *)G_calloc(nnodes + 1, sizeof(double));
	if (!deg)
	    G_fatal_error(_("Out of memory"));
    }

    if (close_opt->answer) {
	closeness = (double *)G_calloc(nnodes + 1, sizeof(double));
	if (!closeness)
	    G_fatal_error(_("Out of memory"));
    }

    if (betw_opt->answer) {
	betw = (double *)G_calloc(nnodes + 1, sizeof(double));
	if (!betw)
	    G_fatal_error(_("Out of memory"));
    }

    if (eigen_opt->answer) {
	eigen = (double *)G_calloc(nnodes + 1, sizeof(double));
	if (!eigen)
	    G_fatal_error(_("Out of memory"));
    }


    if (deg_opt->answer) {
	G_message(_("Computing degree centrality measure"));
	NetA_degree_centrality(graph, deg);
    }
    if (betw_opt->answer || close_opt->answer) {
	G_message(_("Computing betweenness and/or closeness centrality measure"));
	NetA_betweenness_closeness(graph, betw, closeness);
	if (closeness)
	    for (i = 1; i <= nnodes; i++)
		closeness[i] /= (double)In.dgraph.cost_multip;
    }
    if (eigen_opt->answer) {
	G_message(_("Computing eigenvector centrality measure"));
	NetA_eigenvector_centrality(graph, atoi(iter_opt->answer),
				    atof(error_opt->answer), eigen);
    }


    nlines = Vect_get_num_lines(&In);
    G_message(_("Writing data into the table..."));
    G_percent_reset();
    for (i = 1; i <= nlines; i++) {
	G_percent(i, nlines, 1);
	int type = Vect_read_line(&In, Points, Cats, i);

	if (type == GV_POINT && (!chcat || varray->c[i])) {
	    int cat, node;

	    if (!Vect_cat_get(Cats, nfield, &cat))
		continue;
	    Vect_reset_cats(Cats);
	    Vect_cat_set(Cats, 1, cat);
	    Vect_write_line(&Out, type, Points, Cats);
	    node = Vect_find_node(&In, Points->x[0], Points->y[0], Points->z[0], 0, 0);
	    process_node(node, cat);
	    covered[node] = 1;
	}
    }

    if (add_f->answer && !chcat) {
	max_cat = 0;
	for (i = 1; i <= nlines; i++) {
	    Vect_read_line(&In, NULL, Cats, i);
	    for (j = 0; j < Cats->n_cats; j++)
		if (Cats->cat[j] > max_cat)
		    max_cat = Cats->cat[j];
	}
	max_cat++;
	for (i = 1; i <= nnodes; i++)
	    if (!covered[i]) {
		Vect_reset_cats(Cats);
		Vect_cat_set(Cats, 1, max_cat);
		NetA_add_point_on_node(&In, &Out, i, Cats);
		process_node(i, max_cat);
		max_cat++;
	    }

    }

    db_commit_transaction(driver);
    db_close_database_shutdown_driver(driver);

    G_free(covered);
    if (deg)
	G_free(deg);
    if (closeness)
	G_free(closeness);
    if (betw)
	G_free(betw);
    if (eigen)
	G_free(eigen);
    Vect_build(&Out);

    Vect_close(&In);
    Vect_close(&Out);

    exit(EXIT_SUCCESS);
}