From 4d3a144129c9013ec8bd6403c99d15c326e66427 Mon Sep 17 00:00:00 2001 From: ninsbl Date: Fri, 31 Mar 2023 02:09:22 +0200 Subject: [PATCH 1/3] python network --- .../r.connectivity.network.py | 2377 +++++++---------- 1 file changed, 1011 insertions(+), 1366 deletions(-) diff --git a/src/raster/r.connectivity/r.connectivity.network/r.connectivity.network.py b/src/raster/r.connectivity/r.connectivity.network/r.connectivity.network.py index 3a552f4c30..a4823e9ccd 100755 --- a/src/raster/r.connectivity/r.connectivity.network/r.connectivity.network.py +++ b/src/raster/r.connectivity/r.connectivity.network/r.connectivity.network.py @@ -5,7 +5,7 @@ MODULE: r.connectivity.network AUTHOR(S): Stefan Blumentrath PURPOSE: Compute connectivity measures for a set of habitat patches - based on graph-theory (using the igraph-package in R). + based on graph-theory (usig the igraph-package in R). Recently, graph-theory has been characterised as an efficient and useful tool for conservation planning @@ -17,7 +17,7 @@ r.connectivity.network is the 2nd tool of the r.connectivity.* toolchain and performs the (core) network - analysis (using the igraph-package in R) on the network + analysis (usig the igraph-package in R) on the network data prepared with r.connectivity.distance. This network data is analysed on graph, edge and vertex level. @@ -26,8 +26,8 @@ largest cluster, average cluster size, diameter, density, modularity, number of communities, community size (in number of vertices) Network characteristics are visualised - in a plot showing an overview over the number of connections, - number of components and the size of the largest + in a plot showing an overview over number of connections, + number of components and and the size of the largest network component within the network with regards to cost-distance between patches. @@ -143,11 +143,11 @@ # %end # Temporarily disabled -##%flag -##% key: y -##% description: Calculate edge betweenness community (default is do not). -##% guisection: Measures -##%end +## %flag +## % key: y +## % description: Calculate edge betweenness community (default is do not). +## % guisection: Measures +## %end # %option # % key: base @@ -196,11 +196,6 @@ # % answer: 1 # %end -# %flag -# % key: i -# % description: Install required R packages in an interactive session if they are missing -# %end - # %flag # % key: r # % description: Remove indirect connections from network @@ -208,10 +203,15 @@ import atexit import os +from io import BytesIO import sys import platform import warnings import numpy as np +import matplotlib +# Required for Windows +matplotlib.use('wx') +import matplotlib.pyplot as plt import grass.script as grass import grass.script.task as task import grass.script.db as grass_db @@ -221,7 +221,6 @@ if "GISBASE" not in os.environ: sys.exit("You must be in GRASS GIS to run this program") - def cleanup(): """tmp_maps = grass.read_command("g.list", type=['vector', 'raster'], @@ -240,7 +239,7 @@ def is_number(string): """Check if a string can be converted to numeric""" try: - complex(string) # for int, long, float and complex + complex(string) # for int, long, float and complex except ValueError: return False @@ -250,162 +249,109 @@ def is_number(string): def main(): """Do the main work""" + # Lazy import try: - import rpy2 - import rpy2.rinterface - - rpy2.rinterface.set_initoptions( - (b"rpy2", b"--no-save", b"--no-restore", b"--quiet") - ) - import rpy2.robjects as robjects - - # rpy2 throws lots of warnings (that cannot be suppressed) - # when packages are loaded - warnings.filterwarnings("ignore") - import rpy2.robjects.packages as rpackages - from rpy2.robjects.vectors import StrVector - import rpy2.robjects.numpy2ri + import igraph except ImportError: - grass.fatal( - _( - "Cannot import rpy2 (https://rpy2.bitbucket.io)" - " library." - " Please install it (pip install rpy2)" - " or ensure that it is on path" - " (use PYTHONPATH variable)." - ) - ) - - import matplotlib - - matplotlib.use("wxAGG") # required by windows - import matplotlib.pyplot as plt - - # Input variables - network_map = options["input"] + gscript.fatal("Could not import igraph library! To install it run:", + "pip install python-igraph") + #Input variables + network_map = options['input'] # network_mapset = network_map.split('@')[0] # network = network_map.split('@')[1] if len(network_map.split('@')) # > 1 else None - prefix = options["prefix"] - cores = options["cores"] - convergence_treshold = options["convergence_threshold"] - euler = np.exp(1) - base = float(options["base"]) - exponent = float(options["exponent"]) - qml_style_dir = options["qml_style"] - if is_number(options["connectivity_cutoff"]): - if float(options["connectivity_cutoff"]) > 0: - connectivity_cutoff = float(options["connectivity_cutoff"]) + prefix = options['prefix'] + cores = options['cores'] + convergence_treshold = options['convergence_threshold'] + base = float(options['base']) + exponent = float(options['exponent']) + qml_style_dir = options['qml_style'] + if is_number(options['connectivity_cutoff']): + if float(options['connectivity_cutoff']) > 0: + connectivity_cutoff = float(options['connectivity_cutoff']) else: grass.fatal('"connectivity_cutoff" has to be > 0.') else: grass.fatal('Option "connectivity_cutoff" is not given as a number.') - lnbh_cutoff = options["lnbh_cutoff"] - cl_thresh = options["cl_thresh"] + lnbh_cutoff = options['lnbh_cutoff'] + cl_thresh = options['cl_thresh'] cd_cutoff = connectivity_cutoff - kernel_plot = options["kernel_plot"] - overview_plot = options["overview_plot"] + kernel_plot = options['kernel_plot'] + overview_plot = options['overview_plot'] verbose = grass.verbosity() == 3 overwrite = grass.overwrite() - edge_output = "{}_edge_measures".format(prefix) - edge_output_tmp = "{}_edge_measures_tmp".format(prefix) - vertex_output = "{}_vertex_measures".format(prefix) - vertex_output_tmp = "{}_vertex_measures_tmp".format(prefix) - network_output = "{}_network_measures".format(prefix) + edge_output = '{}_edge_measures'.format(prefix) + edge_output_tmp = '{}_edge_measures_tmp'.format(prefix) + vertex_output = '{}_vertex_measures'.format(prefix) + vertex_output_tmp = '{}_vertex_measures_tmp'.format(prefix) + network_output = '{}_network_measures'.format(prefix) # Check if input parameters are correct for table in [vertex_output, edge_output, network_output]: if grass.db.db_table_exist(table) and not overwrite: - grass.fatal( - 'Table "{}" already exists. \ - Use --o flag to overwrite'.format( - table - ) - ) + grass.fatal('Table "{}" already exists. \ + Use --o flag to overwrite'.format(table)) if qml_style_dir: if not os.path.exists(qml_style_dir): - grass.fatal( - 'QML output requested but directory "{}" \ - does not exists.'.format( - qml_style_dir - ) - ) + grass.fatal('QML output requested but directory "{}" \ + does not exists.'.format(qml_style_dir)) if not os.path.isdir(qml_style_dir): - grass.fatal( - 'QML output requested but "{}" is not a \ - directory.'.format( - qml_style_dir - ) - ) + grass.fatal('QML output requested but "{}" is not a \ + directory.'.format(qml_style_dir)) if not os.access(qml_style_dir, os.R_OK): - grass.fatal( - 'Output directory "{}" for QML files is not \ - writable.'.format( - qml_style_dir - ) - ) + grass.fatal('Output directory "{}" for QML files is not \ + writable.'.format(qml_style_dir)) for plot in [kernel_plot, overview_plot]: if plot: if not os.path.exists(os.path.dirname(plot)): - grass.fatal( - 'Directory for output "{}" does not \ - exists.'.format( - plot - ) - ) + grass.fatal('Directory for output "{}" does not \ + exists.'.format(plot)) if not os.access(os.path.dirname(plot), os.R_OK): - grass.fatal( - 'Output directory "{}" is not \ - writable.'.format( - plot - ) - ) + grass.fatal('Output directory "{}" is not \ + writable.'.format(plot)) # Visualise negative exponential decay kernel and exit - x_flag = flags["x"] + x_flag = flags['x'] # Calculate edge betweenness community (default is do not) # y_flag = flags['y'] # Install missing R packages - i_flag = flags["i"] + i_flag = flags['i'] # Remove indirect connections from network - r_flag = flags["r"] + r_flag = flags['r'] db_connection = grass_db.db_connection() # parse_command('db.connect', flags='g') - db_name = grass.read_command("db.databases").rstrip("\n") + db_name = grass.read_command('db.databases').rstrip('\n') - net_hist_str = grass.parse_command( - "v.info", map=network_map, flags="h", delimiter=":" - )["COMMAND"].split("\n")[0] + net_hist_str = grass.parse_command('v.info', map=network_map, + flags='h', delimiter=':' + )['COMMAND'].split('\n')[0] # grass.parse_command('v.info', map=network_map, flags='h' # ).split('\n')[0].split(': ')[1] - # Parsing goes wrong in some cases (extracting with -tp) as input + # Parsing goes wrong in some cases (extractig with -tp) as input dist_cmd_dict = task.cmdstring_to_tuple(net_hist_str) - command = os.environ["CMDLINE"] + command = os.environ['CMDLINE'] - dist_prefix = dist_cmd_dict[1]["prefix"] - pop_proxy = dist_cmd_dict[1]["pop_proxy"] + dist_prefix = dist_cmd_dict[1]['prefix'] + pop_proxy = dist_cmd_dict[1]['pop_proxy'] - in_vertices = "{}_vertices".format(dist_prefix) + in_vertices = '{}_vertices'.format(dist_prefix) # OS adjustment os_type = platform.system() - robjects.numpy2ri.activate() - grass.verbose("prefix is {}".format(prefix)) grass.verbose("cores is {}".format(cores)) grass.verbose("convergence_treshold is {}".format(convergence_treshold)) - grass.verbose("euler is {}".format(euler)) grass.verbose("base is {}".format(base)) grass.verbose("exponent is {}".format(exponent)) grass.verbose("cd_cutoff is {}".format(cd_cutoff)) @@ -418,31 +364,8 @@ def main(): grass.verbose("VERTICES is {}_vertices".format(prefix)) grass.verbose("cl_thresh is {}".format(cl_thresh)) - # Check if R is installed - if not grass.find_program("R"): - grass.fatal( - "R is required, but cannot be found on the system.\n \ - Please make sure that R is installed and the path \ - to R is added to the environment variables \ - (see: http://grass.osgeo.org/wiki/R_statistics#MS_Windows). \ - After that a restart of GRASS GIS is required." - ) - - # Visualise the negative exponential decay kernel and exit (if requested) - if x_flag: - rscript = """euler <- {euler} - basis <- {base} - exponent <- {exponent} - cd_cutoff <- {cd_cutoff} - matplot((0:cd_cutoff/1000), - (euler^((basis*(10^exponent))*(0:cd_cutoff))), - type=\"l\", - xlab=c(\"Cost distance (in 1000)\"), - ylab=\"Potential flow between patches\") - locator(n = 1, type = \"n\")""".format( - euler=euler, base=base, exponent=exponent, cd_cutoff=cd_cutoff - ) + #Visualise the negative exponential decay kernel and exit (if requested) if x_flag or kernel_plot: if x_flag: plt.ion() @@ -451,7 +374,8 @@ def main(): # ax = plt.axes() x_values = np.linspace(0, cd_cutoff, 1000) - plt.plot(x_values, np.exp(base * (np.power(10.0, exponent)) * x_values)) + plt.plot(x_values, np.exp(base*(np.power(10.0, + exponent))*x_values)) # Simplified version could be # plt.plot(x, np.exp(exponent*x)) @@ -464,1248 +388,969 @@ def main(): fig.savefig(kernel_plot) if cores > 1 and os_type == "Windows": - grass.warning( - "Parallel processing not yet supported on MS Windows. \ - Setting number of cores to 1." - ) + grass.warning('Parallel processing not yet supported on MS Windows. \ + Setting number of cores to 1.') cores = 1 - # Check if ghostscript is installed - if overview_plot: - if os_type == "Windows": - # grass.warning("Ghostscript is required. Assuming - # ghostscript is installed on MS Windows") - # Check for ghostscript deactivated on Windows because the - # path to gswin32 is not defined in the environment settings - # by ghostscript installation automatically - if not grass.shutil_which("gswin32") and not grass.shutil_which( - "gswin32.exe" - ): - grass.fatal( - "ghostscript is required for postscript \ - output, cannot find ghostscript \ - (gswin32.exe), please install ghostscript \ - first (or check environment path settings)" - ) - - else: - if not grass.shutil_which("gs") and not grass.shutil_which("ghostscript"): - grass.fatal( - "ghostscript is required for postscript \ - output, please install ghostscript first" - ) - - req_packages = ["igraph", "nlme", "rgrass7", "DBI"] - if cores > 1 and os_type != "Windows": - req_packages = req_packages + ["doMC", "foreach", "iterators", "codetools"] - - rscript = "unique(sort(row.names(installed.packages())))" - installed_packages = robjects.r(rscript) - - missing_packages = np.setdiff1d(req_packages, installed_packages) - - if i_flag and missing_packages: - grass.warning( - "Installing required R-packages {} on \ - request.".format( - ",".join(missing_packages) - ) - ) - utils = rpackages.importr("utils") - utils.chooseCRANmirror(ind=1) - utils.install_packages(StrVector(missing_packages)) - elif missing_packages: - grass.fatal( - "The following required R-packages {} are issing \ - on your system. Please install \ - them.".format( - ",".join(missing_packages) - ) - ) - - # Check igraph version - rscript = """if(as.double(substr( - packageDescription("igraph", - fields = c("Version")),1,3) - )<=0.6) packageDescription("igraph", - fields = c("Version"))""" - igraph_version = robjects.r(rscript) - - if not isinstance(igraph_version, rpy2.rinterface.RNULLType): - if i_flag: - grass.warning( - "The required igraph version is 0.6-2 or \ - later, the installed version is {}. \ - Installing latest version of igraph on \ - request.".format( - igraph_version - ) - ) - utils.install_packages("igraph") - else: - grass.fatal( - "The required igraph version is 0.6-2 or \ - later, the installed version is {}. \ - Please download and install at least igraph \ - version 0.6-2".format( - igraph_version - ) - ) - - rscript = """ -######################################################################## -###Preparation -######################################################################## -verbose <- {} - -if (verbose) {{ -cat("Loading variables and packages...\n") -}} - -library(igraph) -library(nlme) -library(codetools) -library(rgrass7) -library(DBI) -""".format( - "{}".format(verbose).upper() - ) - - if cores > 1: - rscript += """ -# library(multicore) -library(iterators) -library(foreach) -library(doMC) - -registerDoMC() -options(cores = {}) -""".format( - cores - ) - - rscript += """ -# Variables -convergence_threshold <- {convergence_treshold} -euler <- {euler} -basis <- {base} -exponent <- {exponent} -remove_indirect <- {r_flag} -cd_cutoff <- {cd_cutoff} -connectivity_cutoff <- {connectivity_cutoff} -lnbh_cutoff <- {lnbh_cutoff} -cl_thresh <- {cl_thresh} -db <- "{database}" -db_driver <- "{driver}" -db_schema <- "{schema}" -qml_directory <- '{qml_style_dir}' -in_edges <- "{network_map}" -in_vertices <- '{in_vertices}' -overwrite <- {overwrite} -command <- '{command}' -pop_proxy <- '{pop_proxy}' - -#Output files -#Path to kernel plot file -decay_kernel_ps <- '{kernel_plot}' - -#Path to output overview plot file -network_overview_ps <- '{overview_plot}' - -# Tables with connectivity measures -network_output <- '{network_output}' -vertex_output <- '{vertex_output}' -edge_output <- '{edge_output}' -""".format( - convergence_treshold=convergence_treshold, - euler=euler, - base=base, - exponent=exponent, - pop_proxy=pop_proxy, - r_flag="{}".format(r_flag).upper(), - cd_cutoff=cd_cutoff, - connectivity_cutoff=connectivity_cutoff, - lnbh_cutoff=lnbh_cutoff, - cl_thresh=cl_thresh, - kernel_plot=kernel_plot, - overview_plot=overview_plot, - driver=db_connection["driver"], - database=db_name, - schema=db_connection["schema"], - qml_style_dir=qml_style_dir, - in_vertices=in_vertices, - edge_output=edge_output_tmp, - vertex_output=vertex_output_tmp, - network_output=network_output, - network_map=network_map, - overwrite="{}".format(overwrite).upper(), - command=command, - ) - - rscript += """ -#network_measures <- paste(prefix, "network_measures", sep="_") -#edge_output <- paste(prefix, "edge_measures", sep="_") -#vertex_output <- paste(prefix, "vertex_measures", sep="_") -#network_output <- paste(prefix, "network_measures", sep="_") - -if (decay_kernel_ps != "") {{ -if (verbose) {{ -cat("Exporting a plot of the negative exponential decay kernel...\n") -}} - -######################################################################## -#Plot the negative exponential decay kernel -ps.options(encoding="CP1257.enc", family="Helvetica", horizontal=FALSE, fonts="Helvetica", paper="a4", pointsize=1/96, height=3.5, width=3.5) -postscript(decay_kernel_ps) -matplot((0:cd_cutoff/(10^(nchar(as.integer(cd_cutoff))-2))), (euler^((basis*(10^exponent))*(0:cd_cutoff))), type="l", xlab=paste(c("Cost distance (in ", as.character(10^(nchar(as.integer(cd_cutoff))-2)), ")"), sep="", collapse=""), ylab="Potential flow between patches") -off <- dev.off(dev.cur()) -}} - -######################################################################## -###Construct graph -if (verbose) {{ -cat("Reading and preparing input data...\n") -}} - -goutput <- execGRASS('v.db.select', map=in_vertices, columns=paste('cat', pop_proxy, sep=','), separator=',', intern=TRUE) -con <- textConnection(goutput) -v <- read.csv(con, header=TRUE) -names(v) <- c('patch_id', 'pop_proxy') - -#Read and group edge data -goutput <- execGRASS('v.db.select', map=in_edges, columns='cat,from_p,to_p,dist', separator=',', intern=TRUE) -con <- textConnection(goutput) -e_df <- read.csv(con, header=TRUE) - -#e_df <- data.frame(from_p=e_pre[,1], to_p=e_pre[,2], dist=e_pre[,3]) -#e_groups <- groupedData(dist ~ from_p | from_p/to_p, data=as.data.frame(e_df), FUN=mean) -#e_grouped <- gsummary(e_groups, mean) -con_id_u_pre <- unlist(mclapply(mc.cores={0}, 1:length(e_df$from_p), function(x) paste(sort(c(e_df$from_p[x], e_df$to_p[x]))[1], sort(c(e_df$from_p[x], e_df$to_p[x]))[2], sep="_"))) -con_id_u_pre_df <- data.frame(con_id_u_pre=unique(con_id_u_pre), con_id_u=1:length(unique(con_id_u_pre))) -e_grouped_df_pre <- merge(data.frame(con_id_u_pre=con_id_u_pre, con_id=1:length(e_df$from_p), from_p=e_df$from_p, to_p=e_df$to_p, dist=e_df$dist), con_id_u_pre_df, all.x=TRUE) -e_groups_ud <- groupedData(dist ~ 1 | con_id_u, data=e_grouped_df_pre, FUN=mean) -e_grouped_ud <- gsummary(e_groups_ud, mean) -e_grouped_ud_df <- data.frame(con_id_u=as.vector(e_grouped_ud$con_id_u), dist_ud=as.vector(e_grouped_ud$dist)) -e_grouped_df <- merge(data.frame(con_id=as.integer(e_grouped_df_pre$con_id), con_id_u=as.integer(e_grouped_df_pre$con_id_u), from_p=as.integer(e_grouped_df_pre$from_p), to_p=as.integer(e_grouped_df_pre$to_p), dist=as.double(e_grouped_df_pre$dist)), e_grouped_ud_df, all.x=TRUE) - -##Remove connections longer than cost distance threshold if requested -#if(remove_longer_cutoff == 1) {{ -#con_id_true <- e_grouped_df$con_id[grep(TRUE, e_grouped_df$dist_ud10000000000000000,10000000000000000,max(mf_u))-ifelse(min(mf_u)<0.000000000001,0.000000000001,min(mf_u)))/(max(mf_u*-1)-min(mf_u*-1))+ifelse(min(mf_u)<0.000000000001,0.000000000001,min(mf_u)) - -flow_df <- data.frame(from_p=e$from_p, mf_i=mf_i) -sum_mf_i_groups <- groupedData(mf_i ~ from_p | from_p, data=as.data.frame(flow_df), FUN=sum) -sum_mf_i_pre <- gsummary(sum_mf_i_groups, sum) -sum_mf_i <- data.frame(from_p=as.integer(as.character(sum_mf_i_pre[,1])), sum_mf_i=as.double(sum_mf_i_pre[,2])) - -tmp <- merge(merge(e, data.frame(con_id=e$con_id, distance_weight_e_ud=distance_weight_e_ud, distance_weight_e=distance_weight_e, mf_o=mf_o, mf_o_inv=mf_o_inv, mf_i=mf_i, mf_i_inv=mf_i_inv, mf_u=mf_u, mf_inv_u=mf_inv_u)), sum_mf_i) -cf <- tmp$mf_o*(tmp$mf_i/tmp$sum_mf_i) -#cf <- (cf-rep(min(cf),length(cf)))*(1/(rep(max(cf),length(cf))-rep(min(cf),length(cf)))) -#cf_inv <- cf^-1 -#igraphs limits for weights in edge.betweenness function are aproximately max 12600000000000000 and min 0.000000000001 -cf_inv <- ((cf*-1)-min(cf*-1))*(ifelse(max(cf)>10000000000000000,10000000000000000,max(cf))-ifelse(min(cf)<0.000000000001,0.000000000001,min(cf)))/(max(cf*-1)-min(cf*-1))+ifelse(min(cf)<0.000000000001,0.000000000001,min(cf)) - -e <- merge(tmp, data.frame(con_id=e$con_id, cf=cf, cf_inv=cf_inv), by="con_id", sort=TRUE) - -cf_u_df <- data.frame(con_id_u=e$con_id_u, cf=e$cf, cf_inv=e$cf_inv) -cf_u_groups <- groupedData(cf ~ 1 | con_id_u, data=cf_u_df, FUN=sum) -cf_u_grouped <- gsummary(cf_u_groups, sum) - -e <- merge(e, data.frame(con_id_u=as.vector(cf_u_grouped$con_id_u), cf_u=as.vector(cf_u_grouped$cf), cf_inv_u=((as.vector(cf_u_grouped$cf)*-1)-min(as.vector(cf_u_grouped$cf)*-1))*(ifelse(max(as.vector(cf_u_grouped$cf))>10000000000000000,10000000000000000,max(as.vector(cf_u_grouped$cf)))-ifelse(min(as.vector(cf_u_grouped$cf))<0.000000000001,0.000000000001,min(as.vector(cf_u_grouped$cf))))/(max(as.vector(cf_u_grouped$cf)*-1)-min(as.vector(cf_u_grouped$cf)*-1))+ifelse(min(as.vector(cf_u_grouped$cf))<0.000000000001,0.000000000001,min(as.vector(cf_u_grouped$cf)))), all.x=TRUE) - -#Collapse edge data for undirected graph -con_id_u_unique <- unique(e$con_id_u) -ud_groups <- groupedData(con_id ~ 1 | con_id_u, data=data.frame(con_id_u=e$con_id_u, con_id=e$con_id), FUN=mean) -ud_grouped <- gsummary(ud_groups) -ud_grouped_df <- data.frame(con_id_u=as.vector(ud_grouped$con_id_u), con_id_avg=as.vector(ud_grouped$con_id)) -e_ud_pre <- merge(e, ud_grouped_df, all.x=TRUE) -first_con_id <- data.frame(con_id=e_ud_pre$con_id[e_ud_pre$con_id_avg>e_ud_pre$con_id]) -e_ud <- merge(first_con_id, e) - -#Merge vertices and edges to graph-object -#Build directed graph -if (verbose) {{ -cat("Building the graph...\n") -}} - -g <- graph.empty() -g <- add.vertices(g, nrow(v), patch_id=as.character(v$patch_id), pop_proxy=as.numeric(v$pop_proxy)) -names <- V(g)$patch_id -ids <- 1:length(names) -names(ids) <- names -from <- as.character(e$from_p) -to <- as.character(e$to_p) -edges <- matrix(c(ids[from], ids[to]), nc=2) -g <- add.edges(g, t(edges), con_id=e$con_id, con_id_u=e$con_id_u, from_p=e$from_p, from_pop=e$from_pop, to_p=e$to_p, to_pop=e$to_pop, cost_distance=e$cost_distance, cd_u=e$cd_u, distance_weight_e=e$distance_weight_e, distance_weight_e_ud=e$distance_weight_e_ud, mf_o=e$mf_o, mf_o_inv=e$mf_o_inv, mf_i=e$mf_i, mf_i_inv=e$mf_i_inv, mf_u=e$mf_u, mf_inv_u=e$mf_inv_u, cf=e$cf, cf_inv=e$cf_inv, cf_u=e$cf_u, cf_inv_u=e$cf_inv_u) - -#Build undirected graph -g_ud <- as.undirected(graph.empty()) -g_ud <- add.vertices(g_ud, nrow(v), patch_id=as.character(v$patch_id), pop_proxy=as.numeric(v$pop_proxy)) -names <- V(g_ud)$patch_id -ids <- 1:length(names) -names(ids) <- names -from <- as.character(e_ud$from_p) -to <- as.character(e_ud$to_p) -edges <- matrix(c(ids[from], ids[to]), nc=2) -g_ud <- add.edges(g_ud, t(edges), con_id=e_ud$con_id, con_id_u=e_ud$con_id_u, from_p=e_ud$from_p, from_pop=e_ud$from_pop, to_p=e_ud$to_p, to_pop=e_ud$to_pop, cost_distance=e_ud$cost_distance, cd_u=e_ud$cd_u, distance_weight_e=e_ud$distance_weight_e, mf_o=e_ud$mf_o, mf_o_inv=e_ud$mf_o_inv, mf_i=e_ud$mf_i, mf_i_inv=e_ud$mf_i_inv, mf_u=e_ud$mf_u, mf_inv_u=e_ud$mf_inv_u, cf=e_ud$cf, cf_inv=e_ud$cf_inv, cf_u=e_ud$cf_u, cf_inv_u=e_ud$cf_inv_u) - -#Classify connection with regards to shortest paths -#con_id_u <- unlist(lapply(0:(length(E(g))-1), function(x) E(g, path=c(sort(array(c(ends(g, x)[1], ends(g, x)[2])))[1], sort(array(c(ends(g, x)[1], ends(g, x)[2])))[2])))) -#l1 <- unlist(lapply(1:length(V(g_ud)), function(x) rep(x, length(V(g_ud))))) -#l2 <- rep(1:length(V(g_ud)), length(V(g_ud))) -#fa <- function(a) get.shortest.paths(g_ud, l1[a], l2[a], weights=E(g_ud)$cd_u) -#sp <- fa(1:length(l1)) -#Set edge attributes -#E(g)$con_id_u <- con_id_u - -#E(g)$cd_u <- unlist(mclapply(mc.cores={0}, 0:(length(E(g))-1), function(x) (E(g)[x]$cost_distance+E(g, path=c(ends(g, x)[2], ends(g, x)[1]))$cost_distance)/2)) - -#################################### -###Can the following be vectorised???? - -E(g_ud)$isshort_cd <- as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud))), function(x) length(get.shortest.paths(g_ud, from=ends(g_ud, x)[1], to=ends(g_ud, x)[2], weights=E(g_ud)$cd_u)$vpath[[1]])))==2) -E(g_ud)$isshort_mf <- as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud))), function(x) length(get.shortest.paths(g_ud, from=ends(g_ud, x)[1], to=ends(g_ud, x)[2], weights=E(g_ud)$mf_inv_u)$vpath[[1]])))==2) -E(g_ud)$isshort_cf <- as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud))), function(x) length(get.shortest.paths(g_ud, from=ends(g_ud, x)[1], to=ends(g_ud, x)[2], weights=E(g_ud)$cf_inv_u)$vpath[[1]])))==2) -E(g_ud)$isshort <- ifelse((E(g_ud)$isshort_cd==0 & E(g_ud)$isshort_mf==0 & E(g_ud)$isshort_cf==0), 0, 1) -#################################### - -#Add edge attributes to directed graph (for export) -#E(g)$isshort[grep(TRUE, E(g)$con_id_u %in% E(g_ud)$con_id_u[grep(TRUE, E(g_ud)$isshort)])] <- TRUE -#E(g)$isshort[grep(TRUE, E(g)$con_id_u %in% E(g_ud)$con_id_u[grep(FALSE, E(g_ud)$isshort)])] <- FALSE - -#Remove indirect connections if requested -#if(remove_indirect == 1) {{ -g_ud_d <- delete.edges(g_ud, E(g_ud)[E(g_ud)$isshort==FALSE]) -#g_ud <- delete.edges(g_ud, E(g_ud)[E(g_ud)$isshort==FALSE]) -#}} - -#Remove connections above connectivity threshold if requested -if(connectivity_cutoff >= 0.0) {{ -g_ud_cd <- delete.edges(g_ud, E(g_ud)[grep(TRUE, E(g_ud)$cd_u>=connectivity_cutoff)]) -g_ud_d_cd <- delete.edges(g_ud_d, E(g_ud_d)[grep(TRUE, E(g_ud_d)$cd_u>=connectivity_cutoff)]) -}} - -######################################################################## -###Analysis on graph level -######################################################################## - -if (verbose) {{ -cat("Starting analysis on graph level...\n") -}} - -###graph measures -vertices_n <- length(V(g_ud_d)) -edges_n <- length(E(g_ud)) -edges_n_d <- length(E(g_ud_d)) -edges_n_cd <- length(E(g_ud_cd)) -edges_n_d_cd <- length(E(g_ud_d_cd)) - - -######################################################################## -######Calculate number of clusters -#On the undirected graph with only direct edges and on the undirected graph with only direct edges shorter cost distance threshold -cl_no_d <- clusters(g_ud_d)$no -cl_no_d_cd <- clusters(g_ud_d_cd)$no - -cls_size_d <- as.vector(gsummary(groupedData(pop_size ~ 1 | cl, data.frame(cl=clusters(g_ud_d)$membership, pop_size=V(g_ud_d)$pop_proxy), order.groups=TRUE, FUN=sum), sum)$pop_size) -cls_size_d_cd <- as.vector(gsummary(groupedData(pop_size ~ 1 | cl, data.frame(cl=clusters(g_ud_d_cd)$membership, pop_size=V(g_ud_d_cd)$pop_proxy), order.groups=TRUE, FUN=sum), sum)$pop_size) - -cl_max_size_d <- max(cls_size_d) -cl_max_size_d_cd <- max(cls_size_d_cd) - -cl_mean_size_d <- mean(cls_size_d) -cl_mean_size_d_cd <- mean(cls_size_d_cd) - -diam <- diameter(g_ud, directed=FALSE, unconnected=TRUE, weights=E(g_ud)$cd_u) -diam_d <- diameter(g_ud_d, directed=FALSE, unconnected=TRUE, weights=E(g_ud_d)$cd_u) -diam_d_cd <- diameter(g_ud_d_cd, directed=FALSE, unconnected=TRUE, weights=E(g_ud_d_cd)$cd_u) - -density <- graph.density(g, loops=FALSE) -density_u <- graph.density(g_ud, loops=FALSE) -density_ud <- graph.density(g_ud_d, loops=FALSE) -density_udc <- graph.density(g_ud_d_cd, loops=FALSE) - -if (network_overview_ps != "") {{ -######Edge removal operations -idx <- sort(E(g_ud)$cd_u, decreasing=TRUE, na.last=NA, index.return=TRUE)$ix -cl_del_count <- unlist(mclapply(mc.cores={0}, 1:length(idx), function(x) ifelse((clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$no/vertices_n)>=convergence_threshold||E(g_ud)$cd_u[idx[x]]<=(connectivity_cutoff+(connectivity_cutoff*0.25)),clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$no,NA))) -cl_del_max_size <- unlist(mclapply(mc.cores={0}, 1:length(idx), function(x) ifelse((clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$no/vertices_n)>=convergence_threshold||E(g_ud)$cd_u[idx[x]]<=(connectivity_cutoff+(connectivity_cutoff*0.25)),max(as.vector(gsummary(groupedData(pop_size ~ 1 | cl, data.frame(cl=clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$membership, pop_size=V(g_ud)$pop_proxy), order.groups=TRUE, FUN=sum), sum)$pop_size)),NA))) -cl_del_diam <- unlist(mclapply(mc.cores={0}, 1:length(idx), function(x) ifelse((clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$no/vertices_n)>=convergence_threshold||E(g_ud)$cd_u[idx[x]]<=(connectivity_cutoff+(connectivity_cutoff*0.25)),diameter(delete.edges(g_ud, E(g_ud)[idx[1:x]]), directed=FALSE, unconnected=TRUE, weights=E(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$cd_u),NA))) -extract <- sort(grep(FALSE, is.na(cl_del_count)), decreasing=TRUE) -df_edgeremoval <- data.frame(distance=E(g_ud)$cd_u[((idx[extract]))], cl_del_count=cl_del_count[extract], cl_del_max_size=cl_del_max_size[extract], cl_del_diam=cl_del_diam[extract]) -dist_inv_ix <- sort(df_edgeremoval$distance, decreasing=FALSE, na.last=NA, index.return=TRUE)$ix - -ps.options(encoding="CP1257.enc", family="Helvetica", horizontal=FALSE, fonts="Helvetica", paper="a4", pointsize=1/96, height=3.5, width=3.5) -postscript(network_overview_ps) - -###Check axis labels!!! -matplot(df_edgeremoval$distance/(10^(nchar(as.integer(max(df_edgeremoval$distance)))-2)), df_edgeremoval$cl_del_count*100/vertices_n, type="l", ylab="", xlab=c("Connectivity threshold", paste("(Cost distance between patches in ", as.character(10^(nchar(as.integer(max(df_edgeremoval$distance)))-2)), ")", sep="")), lty=1, yaxt="n", yaxs="r", ylim=c(0, 100), yaxs="i", xaxs="i") -if (connectivity_cutoff>0) abline(v=connectivity_cutoff/10^(nchar(as.integer(max(df_edgeremoval$distance)))-2), col="red", lty=3) -ifelse(connectivity_cutoff>0, legend("topleft", - #inset=c(0,-0.15), bty = "n", - xpd=TRUE, - c("Clusters (in % of maximum possible clusters)", "Size of the largest cluster (in % of total population size)", "Number of edges (in % of maximum possible number of edges)", "Diameter (in % of diameter of the entire graph)", "Connectivity threshold used in analysis"), lty=c(1, 2, 3, 4, 3), col=c("black", "black", "black", "black", "red"), inset=0.005, bty="o", box.lty=0, bg="White"), legend("topleft", c("Clusters (in % of maximum possible clusters)", "Size of the largest cluster (in % of total population size)", "Number of edges (in % of maximum possible number of edges)", "Diameter (in % of diameter of the entire graph)"), lty=c(1, 2, 3, 4), col=c("black", "black", "black", "black"), inset=0.005, bty="o", box.lty=0, bg="White")) -axis(2, seq.int(0, 100, 25), labels=c("0 %", "25 %", "50 %", "75 %", "100 %"), yaxs="r") -matplot(df_edgeremoval$distance/(10^(nchar(as.integer(max(df_edgeremoval$distance)))-2)), df_edgeremoval$cl_del_max_size*100/sum(V(g)$pop_proxy), type="l", lty=2, yaxt="n", yaxs="i", add=TRUE) -lines(df_edgeremoval$distance/(10^(nchar(as.integer(max(df_edgeremoval$distance)))-2)), dist_inv_ix*100/edges_n, type="l", lty=3) -lines(df_edgeremoval$distance/(10^(nchar(as.integer(max(df_edgeremoval$distance)))-2)), df_edgeremoval$cl_del_diam*100/diam_d, type="l", lty=4) -#Label axis 4!!! -#axis(4, seq.int(0, 100, 25), yaxs="i", labels=seq.int(0, ceiling((as.integer(edges_n)/(10^(nchar(as.integer(edges_n))-2))))*(10^(nchar(as.integer(edges_n))-2)), ceiling((as.integer(edges_n)/(10^(nchar(as.integer(edges_n))-2))))*(10^(nchar(as.integer(edges_n))-2))/4)) -off <- dev.off(dev.cur()) - -#idx <- sort(E(g_ud)$cd_u, decreasing=TRUE, na.last=NA, index.return=TRUE)$ix -#cl_del_count <- unlist(mclapply(mc.cores={0}, 1:length(idx), function(x) ifelse((clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$no/vertices_n)>=convergence_threshold,clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$no,NA))) -#cl_del_max_size <- unlist(mclapply(mc.cores={0}, 1:length(idx), function(x) ifelse((clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$no/vertices_n)>=convergence_threshold,max(as.vector(gsummary(groupedData(pop_size ~ 1 | cl, data.frame(cl=clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$membership, pop_size=V(g_ud)$pop_proxy), order.groups=TRUE, FUN=sum), sum)$pop_size)),NA))) -#cl_del_diam <- unlist(mclapply(mc.cores={0}, 1:length(idx), function(x) ifelse((clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$no/vertices_n)>=convergence_threshold,diameter(delete.edges(g_ud, E(g_ud)[idx[1:x]]), directed=FALSE, unconnected=TRUE, weights=E(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$cd_u),NA))) -#extract <- sort(grep(FALSE, is.na(cl_del_count)), decreasing=TRUE) -#df_edgeremoval <- data.frame(distance=E(g_ud)$cd_u[((idx[extract]))], cl_del_count=cl_del_count[extract], cl_del_max_size=cl_del_max_size[extract], cl_del_diam=cl_del_diam[extract]) -#dist_inv_ix <- sort(df_edgeremoval$distance, decreasing=FALSE, na.last=NA, index.return=TRUE)$ix - -#ps.options(encoding="CP1257.enc", family="Helvetica", horizontal=FALSE, fonts="Helvetica", paper="a4", pointsize=1/96, height=3.5, width=3.5) -#postscript(network_overview_ps) - -####Check axis labels!!! -#matplot(df_edgeremoval$distance/(10^(nchar(as.integer(max(df_edgeremoval$distance)))-2)), df_edgeremoval$cl_del_count*100/vertices_n, type="l", xlab=c("Grenseverdi for antatt konnektivitet mellom vernområdene", paste("basert på den funksjonale avstanden mellom dem i ", as.character(10^(nchar(as.integer(max(df_edgeremoval$distance)))-2)), sep="")), ylab="Antall vernområder / nettverkskomponenter", lty=1, yaxt="n", yaxs="r", ylim=c(0, 100), yaxs="i", xaxs="i") -#if (connectivity_cutoff>0) abline(v=connectivity_cutoff, col="red", lty=3) -#ifelse(connectivity_cutoff>0, legend("topleft", c("Andel nettverkskomponenter", "Størrelsen av den største nettverkskomponenten", "(i andel populasjons størrelse)", "Antall forbindelser i nettverket", "Andel av diameter av nettverket", "Antatt grenseverdi for konnektivitet"), lty=c(1, 2, 0, 3, 4, 3), col=c("black", "black", "black", "black", "black", "red"), inset=0.005, bty="o", box.lty=0, bg="White"), legend("topleft", c("Andel nettverkskomponenter", "Størrelsen av den største nettverkskomponenten", "(i andel populasjons størrelse)", "Antall forbindelser i nettverket", "Andel av diameter av nettverket"), lty=c(1, 2, 0, 3, 4), col=c("black", "black", "black", "black", "black"), inset=0.005, bty="o", box.lty=0, bg="White")) -#axis(2, seq.int(0, 100, 25), yaxs="r") -#matplot(df_edgeremoval$distance, df_edgeremoval$cl_del_max_size*100/sum(V(g)$pop_proxy), type="l", lty=2, yaxt="n", yaxs="i", add=TRUE) -#lines(df_edgeremoval$distance, dist_inv_ix*100/edges_n, type="l", lty=3) -#lines(df_edgeremoval$distance, df_edgeremoval$cl_del_diam*100/diam_d, type="l", lty=4) -##Label axis 4!!! -##axis(4, seq.int(0, maks_forbindelser_ud*maks_nodes/maks_forbindelser_ud, 100*maks_nodes/maks_forbindelser_ud), yaxs="i", labels=seq.int(0, maks_forbindelser_ud, 100), ylab="Antall forbindelser") -#off <- dev.off(dev.cur()) -}} - -######################################################################## -###Analysis on edge level -######################################################################## - -if (verbose) {{ -cat("Starting analysis on edge level...\n") -}} - -######################################################################## -######Calculate minimum spanning trees (MST) on the undirected graph with only direct edges -###Minimum spanning tree weighted by maximum potential flow - -E(g_ud_d)$mf_mst_ud <- as.integer(0) -E(g_ud_d)[grep(TRUE, E(g_ud_d)$con_id %in% E(minimum.spanning.tree(g_ud_d, weights=E(g_ud_d)$mf_inv_u))$con_id)]$mf_mst_ud <- 1 - -###Minimum spanning tree weighted by cost distance -E(g_ud_d)$cd_mst_ud <- as.integer(0) -E(g_ud_d)[grep(TRUE, E(g_ud_d)$con_id %in% E(minimum.spanning.tree(g_ud_d, weights=E(g_ud_d)$cd_u))$con_id)]$cd_mst_ud <- 1 - -E(g_ud_d)$cf_mst_ud <- as.integer(0) -E(g_ud_d)[grep(TRUE, E(g_ud_d)$con_id %in% E(minimum.spanning.tree(g_ud_d, weights=E(g_ud_d)$cf_inv_u))$con_id)]$cf_mst_ud <- 1 - -######Calculate minimum spanning trees (MST) on the undirected graph with only direct edges shorter cost distance threshold -###Minimum spanning tree weighted by maximum potential flow -E(g_ud_d_cd)$mf_mst_udc <- as.integer(0) -E(g_ud_d_cd)[grep(TRUE, E(g_ud_d_cd)$con_id %in% E(minimum.spanning.tree(g_ud_d_cd, weights=E(g_ud_d_cd)$mf_inv_u))$con_id)]$mf_mst_udc <- 1 - - -###Minimum spanning tree weighted by cost distance -E(g_ud_d_cd)$cd_mst_udc <- as.integer(0) -E(g_ud_d_cd)[grep(TRUE, E(g_ud_d_cd)$con_id %in% E(minimum.spanning.tree(g_ud_d_cd, weights=E(g_ud_d_cd)$cd_u))$con_id)]$cd_mst_udc <- 1 - -###Minimum spanning tree weighted by competing potential flow -E(g_ud_d_cd)$cf_mst_udc <- as.integer(0) -E(g_ud_d_cd)[grep(TRUE, E(g_ud_d_cd)$con_id %in% E(minimum.spanning.tree(g_ud_d_cd, weights=E(g_ud_d_cd)$cf_inv_u))$con_id)]$cf_mst_udc <- 1 - -######################################################################## -###Identify bridges for the undirected graph -E(g_ud)$is_br_u <- as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud))), function(x) clusters(delete.edges(g_ud, E(g_ud)[as.integer(x)]))$no-cl_no_d))) - -#Identify bridges for the undirected graph with only direct edges -E(g_ud_d)$is_br_ud <- as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud_d))), function(x) clusters(delete.edges(g_ud_d, E(g_ud_d)[as.integer(x)]))$no-cl_no_d))) - -###Identify bridges for the undirected graph with only direct edges shorter cost distance threshold -E(g_ud_d_cd)$is_br_udc <- as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud_d_cd))), function(x) clusters(delete.edges(g_ud_d_cd, E(g_ud_d_cd)[as.integer(x)]))$no-cl_no_d_cd))) - -######################################################################## - -biconnected_components_d <- biconnected.components(g_ud) - -#Identify component edges (biconnected components) for the undirected graph with only direct edges -E(g_ud)$bc_e_u <- as.integer(0) -for (c in 1:biconnected_components_d$no) {{ -E(g_ud)$bc_e_u[unlist(biconnected_components_d$component_edges[c])] <- c -}} -#Identify tree edges (biconnected components) for the undirected graph with only direct edges -E(g_ud)$bc_te_u <- as.integer(0) -for (c in 1:biconnected_components_d$no) {{ -E(g_ud)$bc_te_u[unlist(biconnected_components_d$tree_edges[c])] <- c -}} - -biconnected_components_d <- biconnected.components(g_ud_d) -#Identify component edges (biconnected components) for the undirected graph with only direct edges -E(g_ud_d)$bc_e_ud <- as.integer(0) -for (c in 1:biconnected_components_d$no) {{ -E(g_ud_d)$bc_e_ud[unlist(biconnected_components_d$component_edges[c])] <- c -}} -#Identify tree edges (biconnected components) for the undirected graph with only direct edges -E(g_ud_d)$bc_te_ud <- as.integer(0) -for (c in 1:biconnected_components_d$no) {{ -E(g_ud_d)$bc_te_ud[unlist(biconnected_components_d$tree_edges[c])] <- c -}} - -biconnected_components_d_cd <- biconnected.components(g_ud_d_cd) -#Identify component edges (biconnected components) for the undirected graph with only direct edges shorter cost distance threshold -E(g_ud_d_cd)$bc_e_udc <- as.integer(0) -for (c in 1:biconnected_components_d_cd$no) {{ -E(g_ud_d_cd)$bc_e_udc[unlist(biconnected_components_d_cd$component_edges[c])] <- c -}} -#Identify tree edges (biconnected components) for the undirected graph with only direct edges shorter cost distance threshold -E(g_ud_d_cd)$bc_te_udc <- as.integer(0) -for (c in 1:biconnected_components_d_cd$no) {{ -E(g_ud_d_cd)$bc_te_udc[unlist(biconnected_components_d_cd$tree_edges[c])] <- c -}} - -######################################################################## -###Calculate edge betweenness -##Results of the internal edge.betweenness.estimate function are slightly different from the one used here (always larger values): -##see: E(g_ud_d)$cd_leb_ud <- edge.betweenness.estimate(g_ud_d, e=E(g_ud_d), directed=FALSE, lnbh_cutoff*connectivity_cutoff, weights=E(g_ud_d)$cd_u) -##But since the internal edge.betweenness.estimate function is not reasonable to use with other weights than cost distance, the workarounds are used nevertheless -##Further investigation necessary!!! - -###Calculate edge betweenness on the entire undirected graph with only direct edges - -path_lengths <- shortest.paths(g_ud_d, v=V(g_ud_d), to=V(g_ud_d), weights=E(g_ud_d)$cd_u) - -#weighted by cost distance -E(g_ud_d)$cd_eb_ud <- edge.betweenness(g_ud_d, e=E(g_ud_d), directed=FALSE, weights=E(g_ud_d)$cd_u) -E(g_ud_d)$cd_leb_ud <- hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d)), function(x) get.shortest.paths(g_ud_d, x, V(g_ud_d)[grep(TRUE, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d)$cd_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d))), plot=FALSE)$counts/2 -#E(g_ud_d)$cd_leb_ud <- edge.betweenness.estimate(g_ud_d, e=E(g_ud_d), directed=FALSE, lnbh_cutoff*connectivity_cutoff, weights=E(g_ud_d)$cd_u) - -#weighted by maximum potential flow -E(g_ud_d)$mf_eb_ud <- edge.betweenness(g_ud_d, e=E(g_ud_d), directed=FALSE, weights=E(g_ud_d)$mf_inv_u) -E(g_ud_d)$mf_leb_ud <- hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d)), function(x) get.shortest.paths(g_ud_d, x, V(g_ud_d)[grep(TRUE, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d)$mf_inv_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d))), plot=FALSE)$counts/2 - -##weighted by competing potential flow -E(g_ud_d)$cf_eb_ud <- edge.betweenness(g_ud_d, e=E(g_ud_d), directed=FALSE, weights=E(g_ud_d)$cf_inv_u) -E(g_ud_d)$cf_leb_ud <- hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d)), function(x) get.shortest.paths(g_ud_d, x, V(g_ud_d)[grep(TRUE, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d)$cf_inv_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d))), plot=FALSE)$counts/2 - -######################################################################## -###Calculate edge betweenness on the entire undirected graph with only direct edges shorter cost distance threshold - -path_lengths_cd <- shortest.paths(g_ud_d_cd, v=V(g_ud_d_cd), to=V(g_ud_d_cd), weights=E(g_ud_d_cd)$cd_u) - -#weighted by cost distance -E(g_ud_d_cd)$cd_eb_udc <- edge.betweenness(g_ud_d_cd, e=E(g_ud_d_cd), directed=FALSE, weights=E(g_ud_d_cd)$cd_u) -E(g_ud_d_cd)$cd_leb_udc <- hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d_cd)), function(x) get.shortest.paths(g_ud_d_cd, x, V(g_ud_d_cd)[grep(TRUE, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)$cd_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d_cd))), plot=FALSE)$counts/2 - - -#weighted by maximum potential flow -E(g_ud_d_cd)$mf_eb_udc <- edge.betweenness(g_ud_d_cd, e=E(g_ud_d_cd), directed=FALSE, weights=E(g_ud_d_cd)$mf_inv_u) -E(g_ud_d_cd)$mf_leb_udc <- hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d_cd)), function(x) get.shortest.paths(g_ud_d_cd, x, V(g_ud_d_cd)[grep(TRUE, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)$mf_inv_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d_cd))), plot=FALSE)$counts/2 - -##weighted by competing potential flow -E(g_ud_d_cd)$cf_eb_udc <- edge.betweenness(g_ud_d_cd, e=E(g_ud_d_cd), directed=FALSE, weights=E(g_ud_d_cd)$cf_inv_u) -E(g_ud_d_cd)$cf_leb_udc <- hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d_cd)), function(x) get.shortest.paths(g_ud_d_cd, x, V(g_ud_d_cd)[grep(TRUE, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)$cf_inv_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d_cd))), plot=FALSE)$counts/2 - -if (cl_thresh > 0) {{ -######################################################################### -#######Calculate edge betweenness community on the undirected graph -##weighted by competing potential flow -ebc <- edge.betweenness.community(g_ud, weights=E(g_ud)$cf_inv_u, directed=FALSE, edge.betweenness=TRUE, merges=TRUE, bridges=TRUE, modularity=TRUE, membership=TRUE) - -E(g_ud)$cf_iebc_v <- ebc$edge.betweenness -#cf_iebc_m <- ebc$merges -E(g_ud)$cf_iebc_r <- ebc$removed.edges -E(g_ud)$cf_iebc_b <- unlist(mclapply(mc.cores={0}, 1:length(E(g_ud)), function(x) ifelse(x %in% ebc$bridges, 1, 0))) -E(g_ud)$cf_iebc_c <- crossing(ebc, g_ud) - -#V(g)$cf_iebc_mo <- ebc$modularity -V(g)$cf_iebc_me <- ebc$membership - -com_struct <- as.character(cutat(ebc, cl_no_d)) -for (i in (cl_no_d+1):(cl_no_d+cl_thresh)) {{ -com_struct <- paste(com_struct, cutat(ebc, i), sep=';') -}} -V(g)$cf_iebc_cs <- as.character(com_struct) -V(g)$cf_iebc_cl <- cutat(ebc, (cl_no_d+cl_thresh)) - -cf_iebc_u_mod <- modularity(ebc) -com_no_u <- length(ebc) -com_sizes_u <- sizes(ebc) -com_sizes_u_names <- paste("Size of comumity ", names(sizes(ebc)), " (at maximum modularity score (from iebc))", sep="") - -#The following workaround was written when igraph did not support weights when calculating edge.betweenness.community -#The results of this procedure differ significantly from the (new) edge.betweenness.community in igraph 0.6-2 with weights (above) -#In a first visual the results of the following procedure look more reasonable, but further investigation is necessary!!! - -n <- as.integer(1) -ebc_con_id_u <- as.integer(0) -ebc_rank <- as.integer(0) -ebc_value <- as.integer(0) -ebc_clust <- as.integer(0) -V(g_ud)$cf_ebc_cs <- as.character(clusters(g_ud)$membership) -del_eb_g <- g_ud - -for (edge in 1:length(E(del_eb_g))) {{ -##Recalculate edge betweenness for the remaining edges -E(del_eb_g)$cf_eb_ud <- edge.betweenness(del_eb_g, e=E(del_eb_g), directed=FALSE, weights=E(del_eb_g)$cf_inv_u) -##Identify edge with highest edge betweenness value -idx <- sort.int(E(del_eb_g)$cf_eb_ud, decreasing=TRUE, na.last=NA, index.return=TRUE)$ix[1] -##Save edge betweenness value for this edge -ebc_con_id_u[n] <- E(del_eb_g)$con_id_u[idx] -ebc_rank[n] <- n -ebc_clust[n] <- clusters(del_eb_g)$no -if((ebc_clust[n]-cl_no_d)<(cl_no_d+cl_thresh)) {{if(n>1) {{if(ebc_clust[n]>ebc_clust[n-1]) V(g_ud)$cf_ebc_cs <- paste(V(g_ud)$cf_ebc_cs, clusters(del_eb_g)$membership, sep=";")}}}} -ebc_value[n] <- E(del_eb_g)$cf_eb_ud[idx] -##Delete edge with highest edge betweenness value -del_eb_g <- delete.edges(del_eb_g, E(del_eb_g)[idx]) -n <- n+1 -}} - -V(g_ud)$cf_ebc_cl <- unlist(mclapply(mc.cores={0}, 1:length(V(g_ud)), function(x) strsplit(V(g_ud)$cf_ebc_cs, ";")[[x]][(cl_no_d+cl_thresh)])) - -#ebc_df <- merge(data.frame(con_id=E(g)$con_id, con_id_u=E(g)$con_id_u), data.frame(con_id_u=ebc_con_id_u, ebc_rank=ebc_rank, ebc_clust=ebc_clust, ebc_value=ebc_value), all.x=TRUE) -idx <- sort.int(ebc_con_id_u, index.return=TRUE)$ix -E(g_ud)$cf_ebc_v <- ebc_value[idx] -E(g_ud)$cf_ebc_r <- ebc_rank[idx] -E(g_ud)$cf_ebc_c <- ebc_clust[idx] - -E(g_ud)$cf_ebc_vi <- strftime(Sys.time()+E(g_ud)$cf_ebc_r) - -######################################################################### -###Calculate community connectors (based on community structure from ebc) -com_membership <- data.frame(patch_id=V(g_ud_d_cd)$patch_id, com=V(g_ud)$cf_ebc_cl) -com_pc_pre <- merge(merge(data.frame(id=1:length(E(g)), from_p=E(g)$from_p, to_p=E(g)$to_p), com_membership, by.x="from_p", by.y="patch_id"), com_membership, by.x="to_p", by.y="patch_id")[order(merge(merge(data.frame(id=1:length(E(g)), from_p=E(g)$from_p, to_p=E(g)$to_p), com_membership, by.x="from_p", by.y="patch_id"), com_membership, by.x="to_p", by.y="patch_id")$id) ,] -E(g)$cf_ebc_cc <- ifelse(com_pc_pre$com.x==com_pc_pre$com.y,0,1) - - -}} else {{ -if (verbose) {{ -cat("Skipping community algorithms...\n") -}} -}} -######################################################################## -######Calculate cluster membership -V(g_ud_d)$cl_ud <- clusters(g_ud_d)$membership -V(g_ud_d_cd)$cl_udc <- clusters(g_ud_d_cd)$membership - -######################################################################### -###Calculate potential cluster connectors (based on cost distance threshold) -cl_membership <- data.frame(patch_id=V(g_ud_d_cd)$patch_id, cl=V(g_ud_d_cd)$cl_udc) -cl_pc_pre <- merge(merge(data.frame(id=1:length(E(g)), from_p=E(g)$from_p, to_p=E(g)$to_p), cl_membership, by.x="from_p", by.y="patch_id"), cl_membership, by.x="to_p", by.y="patch_id")[order(merge(merge(data.frame(id=1:length(E(g)), from_p=E(g)$from_p, to_p=E(g)$to_p), cl_membership, by.x="from_p", by.y="patch_id"), cl_membership, by.x="to_p", by.y="patch_id")$id) ,] -E(g)$cl_pc <- ifelse(cl_pc_pre$cl.x==cl_pc_pre$cl.y,0,1) - -if (verbose) {{ -cat("Starting analysis on vertex level...\n") -}} - -######################################################################## -###Analysis on vertex level -######################################################################## - -######################################################################## -######Calculate degree centrality -#On the directed graph -#V(g)$deg_dir <- as.integer(degree(g, v=V(g), mode=c("in"))) - -#On the entire undirected graph -V(g_ud)$deg_u <- as.integer(degree(g_ud, v=V(g_ud))) - -#On the undirected graph with edges shorter cost distance threshold -V(g_ud_cd)$deg_uc <- as.integer(degree(g_ud_cd, v=V(g_ud_cd))) - -#On the undirected graph with only direct edges -V(g_ud_d)$deg_ud <- as.integer(degree(g_ud_d, v=V(g_ud_d))) - -#On the undirected graph with only direct edges shorter cost distance threshold -V(g_ud_d_cd)$deg_udc <- as.integer(degree(g_ud_d_cd, v=V(g_ud_d_cd))) - -######################################################################## -######Calculate closeness centrality -V(g_ud_d)$cd_cl_ud <- closeness(g_ud_d, vids=V(g_ud_d), weights = E(g_ud_d)$cd_u, normalized = FALSE) -V(g_ud_d)$mf_cl_ud <- closeness(g_ud_d, vids=V(g_ud_d), weights = E(g_ud_d)$mf_inv_u, normalized = FALSE) -V(g_ud_d)$cf_cl_ud <- closeness(g_ud_d, vids=V(g_ud_d), weights = E(g_ud_d)$cf_inv_u, normalized = FALSE) - -V(g_ud_d_cd)$cd_cl_udc <- closeness(g_ud_d_cd, vids=V(g_ud_d_cd), weights = E(g_ud_d_cd)$cd_u, normalized = FALSE) -V(g_ud_d_cd)$mf_cl_udc <- closeness(g_ud_d_cd, vids=V(g_ud_d_cd), weights = E(g_ud_d_cd)$mf_inv_u, normalized = FALSE) -V(g_ud_d_cd)$cf_cl_udc <- closeness(g_ud_d_cd, vids=V(g_ud_d_cd), weights = E(g_ud_d_cd)$cf_inv_u, normalized = FALSE) - -######################################################################## -######Calculate biconnected components - -##On the undirected graph with only direct edges -#Number of new clusters when a vertex is deleted -V(g_ud_d)$art_ud <- as.integer(unlist(mclapply(mc.cores={0}, 1:(length(V(g_ud_d))), function(x) ifelse((clusters(delete.vertices(g_ud_d, V(g_ud_d)[x]))$no-cl_no_d)<0,0,clusters(delete.vertices(g_ud_d, V(g_ud_d)[x]))$no-cl_no_d)))) -#Vertex is articulation point -V(g_ud_d)$art_p_ud <- as.integer(0) -V(g_ud_d)$art_p_ud[biconnected_components_d$articulation_points] <- as.integer(1) - -##On the undirected graph with only direct edges shorter cost distance threshold -#Number of new clusters when a vertex is deleted -V(g_ud_d_cd)$art_udc <- unlist(mclapply(mc.cores={0}, 1:(length(V(g_ud_d_cd))), function(x) ifelse((clusters(delete.vertices(g_ud_d_cd, V(g_ud_d_cd)[x]))$no-cl_no_d_cd)<0,0,clusters(delete.vertices(g_ud_d_cd, V(g_ud_d_cd)[x]))$no-cl_no_d_cd))) -#Vertex is articulation point -V(g_ud_d_cd)$art_p_udc <- as.integer(0) -V(g_ud_d_cd)$art_p_udc[biconnected_components_d_cd$articulation_points] <- as.integer(1) - -######################################################################## -######Calculate eigenvector centrality -#Calculate eigenvector for the entire directed graph - -#Calculate eigenvector centrality weighted by competing potential flow -cf_evc_groups <- groupedData(cf ~ 1 | to_p, data=e, FUN=sum) -cf_evc <- gsummary(cf_evc_groups, sum) -df_cf_evc_d <- data.frame(patch_id=cf_evc$to_p, cf_evc_d=cf_evc$cf) - -#Calculate eigenvector centrality weighted by maximum potential flow -mf_o_evc_groups <- groupedData(mf_o ~ 1 | to_p, data=e, FUN=sum) -mf_o_evc <- gsummary(mf_o_evc_groups, sum) -df_mf_i_evc_d <- data.frame(patch_id=mf_o_evc$to_p, mf_evc_d=mf_o_evc$mf_o) - -df_evc_d <- merge(df_cf_evc_d, df_mf_i_evc_d) - -######Calculate eigenvector centrality on the directed graph with edges shorter cost distance threshold -#Calculate eigenvector centrality weighted by competing potential flow -evc_udc_pre <- data.frame(to_p=e$to_p[grep(TRUE, e$cd_u 0) {{ -grml <- rbind(grml, -c("Density of the graph with only edges shorter cost distance threshold", density_udc), -c("Modularity (from iebc) of the entire graph (undirected) weighted by cf", cf_iebc_u_mod), -c("Number of communities (at maximum modularity score (from iebc)) of the entire (undirected) graph weighted by cf", com_no_u), -c("com_sizes_u", toString(com_sizes_u, sep=',')), -c("com_sizes_u_names", toString(com_sizes_u_names, sep=',')) -) -}} - -network_overview_measures <- data.frame( - measure=grml[,1], - value=grml[,2]) - -# Write dataframe -dbWriteTable(con, network_output, network_overview_measures, overwrite=overwrite, row.names=FALSE) - -############################################################### -#####Merge vertexmeasures in a dataframe and save it - -###Create initial export data frame - -vertex_export_df_ud <- as.data.frame(1:length(V(g_ud))) - -####Adjust vertex-attributes first - -g_ud <- remove.vertex.attribute(g_ud, "pop_proxy") - -vertex_attribute_list <- list.vertex.attributes(g_ud) -if(length(vertex_attribute_list)>0) {{ - for (vl in vertex_attribute_list) {{ - vertex_export_df_ud <- as.data.frame(cbind(vertex_export_df_ud, get.vertex.attribute(g_ud, vl))) - }} - - vertex_export_df_ud <- vertex_export_df_ud[2:length(vertex_export_df_ud)] - names(vertex_export_df_ud) <- vertex_attribute_list -}} - -###Create initial export data frame - -vertex_export_df_ud_d <- as.data.frame(1:length(V(g_ud_d))) - -####Adjust vertex-attributes first - -g_ud_d <- remove.vertex.attribute(g_ud_d, "pop_proxy") - -vertex_attribute_list <- list.vertex.attributes(g_ud_d) -if(length(vertex_attribute_list)>0) {{ - for (vl in vertex_attribute_list) {{ - vertex_export_df_ud_d <- as.data.frame(cbind(vertex_export_df_ud_d, get.vertex.attribute(g_ud_d, vl))) - }} - vertex_export_df_ud_d <- vertex_export_df_ud_d[2:length(vertex_export_df_ud_d)] - names(vertex_export_df_ud_d) <- vertex_attribute_list -}} + #Check igraph version + if igraph.__version_info__ < (0, 6, 2): + grass.fatal('The required igraph version is 0.6.2 or \ + later, the installed version is {}. \ + Please download and install at least igraph \ + version 0.6.2'.format(".".join(igraph.__version_info__))) + + driver = db_connection['driver'] + database = db_name + schema = db_connection['schema'] + command = command + + gscript.verbose("Reading and preparing input data...\n") + vertices = np.gefromtext(BytesIO(gscript.read_command('v.db.select', flags="c", map=in_vertices, columns="cat,{}".format(pop_proxy), separator=',').strip("\n")), dtype=None, names=['patch_id', 'pop_proxy'], delimiter=",") + + #Read edge data + edges = np.gefromtext(BytesIO(gscript.read_command('v.db.select', flags="c", map=in_edges, columns='cat,from_p,to_p,dist', separator=',').strip("\n")), dtype=None, names=['cat', 'from_patch', 'to_patch', 'distance'], delimiter=",") + + # THe following can be done in igraph graph.simplify(combine_edges=dict(weight="mean")) + # https://igraph.org/python/doc/igraph.GraphBase-class.html#simplify + """con_id_u_pre = np.unique(np.sort(e_df[:,[1,2]]), axis=0) + con_id_u_pre_df = data.frame(con_id_u_pre=unique(con_id_u_pre), con_id_u=1:length(unique(con_id_u_pre))) + e_grouped_df_pre = merge(data.frame(con_id_u_pre=con_id_u_pre, con_id=1:length(e_df$from_p), from_p=e_df$from_p, to_p=e_df$to_p, dist=e_df$dist), con_id_u_pre_df, all.x=True) + e_groups_ud = groupedData(dist ~ 1 | con_id_u, data=e_grouped_df_pre, FUN=mean) + e_grouped_ud = gsummary(e_groups_ud, mean) + e_grouped_ud_df = data.frame(con_id_u=as.vector(e_grouped_ud$con_id_u), dist_ud=as.vector(e_grouped_ud$dist)) + e_grouped_df = merge(data.frame(con_id=as.integer(e_grouped_df_pre$con_id), con_id_u=as.integer(e_grouped_df_pre$con_id_u), from_p=as.integer(e_grouped_df_pre$from_p), to_p=as.integer(e_grouped_df_pre$to_p), dist=as.double(e_grouped_df_pre$dist)), e_grouped_ud_df, all.x=True)""" + + # Merge vertices and edges to graph-object + # Build directed graph + gscript.verbose("Building the graph...") + + g = Graph().as_directed() + g.add_vertices(vertices[:,0].astype("i").astype("S")) + g.vs["patch_id"] = vertices[:,0].astype("i") + g.vs["pop_proxy"] = vertices[:,1] + + g.add_edges(edges[:,[1,2]].astype("i").astype("S")) + g.es["cat"] = edges[:,0].astype("i") + g.es["from_patch"] = edges[:,1].astype("i").astype("S") + g.es["to_patch"] = edges[:,2].astype("i").astype("S") + g.es["cost_distance"] = edges[:,3] + + #g = add.edges(g, t(edges), con_id=e$con_id, con_id_u=e$con_id_u, from_p=e$from_p, from_pop=e$from_pop, to_p=e$to_p, to_pop=e$to_pop, cost_distance=e$cost_distance, cd_u=e$cd_u, distance_weight_e=e$distance_weight_e, distance_weight_e_ud=e$distance_weight_e_ud, mf_o=e$mf_o, mf_o_inv=e$mf_o_inv, mf_i=e$mf_i, mf_i_inv=e$mf_i_inv, mf_u=e$mf_u, mf_inv_u=e$mf_inv_u, cf=e$cf, cf_inv=e$cf_inv, cf_u=e$cf_u, cf_inv_u=e$cf_inv_u) + + ##Remove connections longer than cost distance threshold if requested + #if(remove_longer_cutoff == 1) {{ + #con_id_True = e_grouped_df$con_id[grep(True, e_grouped_df$dist_ud10000000000000000,10000000000000000,max(g.es["maximum_flow_sum"]))-ifelse(min(g.es["maximum_flow_sum"])<0.000000000001,0.000000000001,min(g.es["maximum_flow_sum"])))/(max(g.es["maximum_flow_sum"]*-1)-min(g.es["maximum_flow_sum"]*-1))+ifelse(min(g.es["maximum_flow_sum"])<0.000000000001,0.000000000001,min(g.es["maximum_flow_sum"])) + + flow_df = data.frame(from_p=e$from_p, mf_i=mf_i) + sum_mf_i_groups = groupedData(mf_i ~ from_p | from_p, data=as.data.frame(flow_df), FUN=sum) + sum_mf_i_pre = gsummary(sum_mf_i_groups, sum) + sum_mf_i = data.frame(from_p=as.integer(as.character(sum_mf_i_pre[,1])), sum_mf_i=as.double(sum_mf_i_pre[,2])) + + tmp = merge(merge(e, data.frame(con_id=e$con_id, distance_weight_e_ud=distance_weight_e_ud, distance_weight_e=distance_weight_e, mf_o=mf_o, mf_o_inv=mf_o_inv, mf_i=mf_i, mf_i_inv=mf_i_inv, mf_u=mf_u, mf_inv_u=mf_inv_u)), sum_mf_i) + cf = tmp$mf_o*(tmp$mf_i/tmp$sum_mf_i) + # cf = (cf-rep(min(cf),length(cf)))*(1/(rep(max(cf),length(cf))-rep(min(cf),length(cf)))) + # cf_inv = cf^-1 + # igraphs limits for weights in edge.betweenness function are aproximately max 12600000000000000 and min 0.000000000001 + cf_inv = ((cf*-1)-min(cf*-1))*(ifelse(max(cf)>10000000000000000,10000000000000000,max(cf))-ifelse(min(cf)<0.000000000001,0.000000000001,min(cf)))/(max(cf*-1)-min(cf*-1))+ifelse(min(cf)<0.000000000001,0.000000000001,min(cf)) + + e = merge(tmp, data.frame(con_id=e$con_id, cf=cf, cf_inv=cf_inv), by="con_id", sort=True) + + cf_u_df = data.frame(con_id_u=e$con_id_u, cf=e$cf, cf_inv=e$cf_inv) + cf_u_groups = groupedData(cf ~ 1 | con_id_u, data=cf_u_df, FUN=sum) + cf_u_grouped = gsummary(cf_u_groups, sum) + + e = merge(e, data.frame(con_id_u=as.vector(cf_u_grouped$con_id_u), cf_u=as.vector(cf_u_grouped$cf), cf_inv_u=((as.vector(cf_u_grouped$cf)*-1)-min(as.vector(cf_u_grouped$cf)*-1))*(ifelse(max(as.vector(cf_u_grouped$cf))>10000000000000000,10000000000000000,max(as.vector(cf_u_grouped$cf)))-ifelse(min(as.vector(cf_u_grouped$cf))<0.000000000001,0.000000000001,min(as.vector(cf_u_grouped$cf))))/(max(as.vector(cf_u_grouped$cf)*-1)-min(as.vector(cf_u_grouped$cf)*-1))+ifelse(min(as.vector(cf_u_grouped$cf))<0.000000000001,0.000000000001,min(as.vector(cf_u_grouped$cf)))), all.x=True) + + # Collapse edge data for undirected graph + con_id_u_unique = np.unique(e$con_id_u) + ud_groups = groupedData(con_id ~ 1 | con_id_u, data=data.frame(con_id_u=e$con_id_u, con_id=e$con_id), FUN=mean) + ud_grouped = gsummary(ud_groups) + ud_grouped_df = data.frame(con_id_u=as.vector(ud_grouped$con_id_u), con_id_avg=as.vector(ud_grouped$con_id)) + e_ud_pre = merge(e, ud_grouped_df, all.x=True) + first_con_id = data.frame(con_id=e_ud_pre$con_id[e_ud_pre$con_id_avg>e_ud_pre$con_id]) + e_ud = merge(first_con_id, e) + + # Build undirected graph + # graph.simplify(combine_edges=dict(weight="mean")) + # https://igraph.org/python/doc/igraph.GraphBase-class.html#simplify + g_ud = deepcopy(g) + g_ud = g_ud.simplify(combine_edges=dict(mf_o="mean", mf_i="mean", )) + g_ud.es["exponential_distance_weight_undirected"] = np.exp((basis*(10^exponent))*np.array(g_ud.es["cost_distance_undirected"])) + + """ + g_ud = add.edges(g_ud, t(edges), con_id=e_ud$con_id, + con_id_u=e_ud$con_id_u, + from_p=e_ud$from_p, + from_pop=e_ud$from_pop, + to_p=e_ud$to_p, + to_pop=e_ud$to_pop, + cost_distance=e_ud$cost_distance, + cd_u=e_ud$cd_u, + distance_weight_e=e_ud$distance_weight_e, + mf_o=e_ud$mf_o, + mf_o_inv=e_ud$mf_o_inv, + mf_i=e_ud$mf_i, + mf_i_inv=e_ud$mf_i_inv, + mf_u=e_ud$mf_u, + mf_inv_u=e_ud$mf_inv_u, + cf=e_ud$cf, + cf_inv=e_ud$cf_inv, + cf_u=e_ud$cf_u, + cf_inv_u=e_ud$cf_inv_u) + """ + + # Classify connection with regards to shortest paths + # con_id_u = unlist(lapply(0:(length(E(g))-1), function(x) E(g, path=c(sort(array(c(ends(g, x)[1], ends(g, x)[2])))[1], sort(array(c(ends(g, x)[1], ends(g, x)[2])))[2])))) + # l1 = unlist(lapply(1:length(V(g_ud)), function(x) rep(x, length(V(g_ud))))) + # l2 = rep(1:length(V(g_ud)), length(V(g_ud))) + # fa = function(a) get.shortest.paths(g_ud, l1[a], l2[a], weights=g_ud.es["cd_u) + # sp = fa(1:length(l1)) + # Set edge attributes + # g.es["con_id_u = con_id_u + + # g.es["cd_u = unlist(mclapply(mc.cores={0}, 0:(length(E(g))-1), function(x) (E(g)[x]$cost_distance+E(g, path=c(ends(g, x)[2], ends(g, x)[1]))$cost_distance)/2)) + + #################################### + ### Can the following be vectorised???? + + g_ud.es["is_short_cd"] = as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud))), function(x) length(get.shortest.paths(g_ud, from=ends(g_ud, x)[1], to=ends(g_ud, x)[2], weights=g_ud.es["cd_u)$vpath[[1]])))==2) + g_ud.es["is_short_mf"] = as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud))), function(x) length(get.shortest.paths(g_ud, from=ends(g_ud, x)[1], to=ends(g_ud, x)[2], weights=g_ud.es["mf_inv_u)$vpath[[1]])))==2) + g_ud.es["is_short_cf"] = as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud))), function(x) length(get.shortest.paths(g_ud, from=ends(g_ud, x)[1], to=ends(g_ud, x)[2], weights=g_ud.es["cf_inv_u)$vpath[[1]])))==2) + # g_ud.es["is_short"] = ifelse((g_ud.es["isshort_cd==0 & g_ud.es["isshort_mf==0 & g_ud.es["isshort_cf==0), 0, 1) + g_ud.es["is_short"] = np.min(g_ud.es["is_short_cd"], g_ud.es["is_short_mf"], g_ud.es["is_short_cf"]) + + #################################### + + # Add edge attributes to directed graph (for export) + # g.es["isshort[grep(True, g.es["con_id_u %in% g_ud.es["con_id_u[grep(True, g_ud.es["isshort)])] = True + # g.es["isshort[grep(True, g.es["con_id_u %in% g_ud.es["con_id_u[grep(False, g_ud.es["isshort)])] = False + + # Remove indirect connections if requested + if remove_indirect == 1: + g_ud_d = deepcopy(g_ud) + g_ud_d.delete.edges(g_ud.es.select[is_short_eq=0]) + + # Remove connections above connectivity threshold if requested + if connectivity_cutoff >= 0.0: + g_ud_cd = deepcopy(g_ud) + g_ud_cd = g_ud_cd.delete_edges(E(g_ud)[grep(True, g_ud.es["cd_u>=connectivity_cutoff)]) + g_ud_d_cd = deepcopy(g_ud_d) + g_ud_d_cd = g_ud_d_cd.delete_edges(E(g_ud_d)[grep(True, g_ud_d.es["cd_u>=connectivity_cutoff)]) + + ######################################################################## + ### Analysis on graph level + ######################################################################## + + gscript.verbose("Starting analysis on graph level...") + + ### graph measures + vertices_n = len(g_ud.vs) + edges_n = len(g_ud.es) + edges_n_d = len(g_ud_d.es) + edges_n_cd = len(g_ud_cd.es) + edges_n_d_cd = len(g_ud_d_cd.es) + + + ######################################################################## + ###### Calculate number of clusters + # On the undirected graph with only direct edges and on the undirected graph with only direct edges shorter cost distance threshold + cl_no_d = g_ud_d.clusters()["no"] + cl_no_d_cd = g_ud_d_cd.clusters()["no"] + + cls_size_d = as.vector(gsummary(groupedData(pop_size ~ 1 | cl, data.frame(cl=clusters(g_ud_d)$membership, pop_size=g_ud_d.vs["pop_proxy), order.groups=True, FUN=sum), sum)$pop_size) + cls_size_d_cd = as.vector(gsummary(groupedData(pop_size ~ 1 | cl, data.frame(cl=clusters(g_ud_d_cd)$membership, pop_size=g_ud_d_cd.vs["pop_proxy), order.groups=True, FUN=sum), sum)$pop_size) + + cl_max_size_d = max(cls_size_d) + cl_max_size_d_cd = max(cls_size_d_cd) + + cl_mean_size_d = mean(cls_size_d) + cl_mean_size_d_cd = mean(cls_size_d_cd) + + diam = g_ud.diameter(directed=False, unconnected=True, weights=g_ud.es["cd_u) + diam_d = g_ud_d.diameter(directed=False, unconnected=True, weights=g_ud_d.es["cd_u) + diam_d_cd = g_ud_d_cd.diameter(directed=False, unconnected=True, weights=E(g_ud_d_cd)$cd_u) + + density = g.density(g, loops=False) + density_u = g_ud.density(g_ud, loops=False) + density_ud = g_ud_d.density(g_ud_d, loops=False) + density_udc = g_ud_d_cd.density(g_ud_d_cd, loops=False) + + if network_overview_ps: + ###### Edge removal operations + for edge_length in np.sort(g_ud.es["cd_u"], decreasing=True, na.last=NA): + cl_del_count = unlist(mclapply(mc.cores={0}, 1:length(idx), function(x) ifelse((clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$no/vertices_n)>=convergence_threshold||g_ud.es["cd_u[idx[x]]<=(connectivity_cutoff+(connectivity_cutoff*0.25)),clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$no,NA))) + cl_del_max_size = unlist(mclapply(mc.cores={0}, 1:length(idx), function(x) ifelse((clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$no/vertices_n)>=convergence_threshold||g_ud.es["cd_u[idx[x]]<=(connectivity_cutoff+(connectivity_cutoff*0.25)),max(as.vector(gsummary(groupedData(pop_size ~ 1 | cl, data.frame(cl=clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$membership, pop_size=g_ud.vs["pop_proxy), order.groups=True, FUN=sum), sum)$pop_size)),NA))) + cl_del_diam = unlist(mclapply(mc.cores={0}, 1:length(idx), function(x) ifelse((clusters(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$no/vertices_n)>=convergence_threshold||g_ud.es["cd_u[idx[x]]<=(connectivity_cutoff+(connectivity_cutoff*0.25)),diameter(delete.edges(g_ud, E(g_ud)[idx[1:x]]), directed=False, unconnected=True, weights=E(delete.edges(g_ud, E(g_ud)[idx[1:x]]))$cd_u),NA))) + extract = sort(grep(False, is.na(cl_del_count)), decreasing=True) + df_edgeremoval = data.frame(distance=g_ud.es["cd_u[((idx[extract]))], cl_del_count=cl_del_count[extract], cl_del_max_size=cl_del_max_size[extract], cl_del_diam=cl_del_diam[extract]) + dist_inv_ix = sort(df_edgeremoval$distance, decreasing=False, na.last=NA, index.return=True)$ix + + ps.options(encoding="CP1257.enc", family="Helvetica", horizontal=False, fonts="Helvetica", paper="a4", pointsize=1/96, height=3.5, width=3.5) + postscript(network_overview_ps) + + ### Check axis lables!!! + matplot(df_edgeremoval$distance/(10^(nchar(as.integer(max(df_edgeremoval$distance)))-2)), df_edgeremoval$cl_del_count*100/vertices_n, type="l", ylab="", xlab=c("Connectivity threshold", paste("(Cost distance between patches in ", as.character(10^(nchar(as.integer(max(df_edgeremoval$distance)))-2)), ")", sep="")), lty=1, yaxt="n", yaxs="r", ylim=c(0, 100), yaxs="i", xaxs="i") + if (connectivity_cutoff>0) abline(v=connectivity_cutoff/10^(nchar(as.integer(max(df_edgeremoval$distance)))-2), col="red", lty=3) + ifelse(connectivity_cutoff>0, legend("topleft", + #inset=c(0,-0.15), bty = "n", + xpd=True, + c("Clusters (in % of maximum possible clusters)", "Size of the largest cluster (in % of total population size)", "Number of edges (in % of maximum possible number of edges)", "Diameter (in % of diameter of the entire graph)", "Connectivity threshold used in analysis"), lty=c(1, 2, 3, 4, 3), col=c("black", "black", "black", "black", "red"), inset=0.005, bty="o", box.lty=0, bg="White"), legend("topleft", c("Clusters (in % of maximum possible clusters)", "Size of the largest cluster (in % of total population size)", "Number of edges (in % of maximum possible number of edges)", "Diameter (in % of diameter of the entire graph)"), lty=c(1, 2, 3, 4), col=c("black", "black", "black", "black"), inset=0.005, bty="o", box.lty=0, bg="White")) + axis(2, seq.int(0, 100, 25), labels=c("0 %", "25 %", "50 %", "75 %", "100 %"), yaxs="r") + matplot(df_edgeremoval$distance/(10^(nchar(as.integer(max(df_edgeremoval$distance)))-2)), df_edgeremoval$cl_del_max_size*100/sum(g.vs["pop_proxy), type="l", lty=2, yaxt="n", yaxs="i", add=True) + lines(df_edgeremoval$distance/(10^(nchar(as.integer(max(df_edgeremoval$distance)))-2)), dist_inv_ix*100/edges_n, type="l", lty=3) + lines(df_edgeremoval$distance/(10^(nchar(as.integer(max(df_edgeremoval$distance)))-2)), df_edgeremoval$cl_del_diam*100/diam_d, type="l", lty=4) + # Lable axis 4!!! + # axis(4, seq.int(0, 100, 25), yaxs="i", labels=seq.int(0, ceiling((as.integer(edges_n)/(10^(nchar(as.integer(edges_n))-2))))*(10^(nchar(as.integer(edges_n))-2)), ceiling((as.integer(edges_n)/(10^(nchar(as.integer(edges_n))-2))))*(10^(nchar(as.integer(edges_n))-2))/4)) + + ######################################################################## + ###Analysis on edge level + ######################################################################## + + gscript.verbose("Starting analysis on edge level...") + + ######################################################################## + ###### Calculate minimum spanning trees (MST) on the undirected graph with only direct edges + ### Minimum spanning tree weighted by maximum potential flow + + g_ud_d.es["mf_mst_ud"] = 0 + E(g_ud_d.es["mf_mst_ud"] = [grep(True, g_ud_d.es["con_id %in% E(minimum.spanning.tree(g_ud_d, weights=g_ud_d.es["mf_inv_u))$con_id)]$mf_mst_ud = 1 + + ### Minimum spanning tree weighted by cost distance + g_ud_d.es["cd_mst_ud"] = 0 + g_ud_d.es["cd_mst_ud"] = [grep(True, g_ud_d.es["con_id %in% E(minimum.spanning.tree(g_ud_d, weights=g_ud_d.es["cd_u))$con_id)]$cd_mst_ud = 1 + + g_ud_d.es["cf_mst_ud"] = 0 + g_ud_d.es["cf_mst_ud"] = [grep(True, g_ud_d.es["con_id %in% E(minimum.spanning.tree(g_ud_d, weights=g_ud_d.es["cf_inv_u))$con_id)]$cf_mst_ud = 1 + + ###### Calculate minimum spanning trees (MST) on the undirected graph with only direct edges shorter cost distance threshold + ### Minimum spanning tree weighted by maximum potential flow + g_ud_d_cd.es["mf_mst_udc"] = 0 + g_ud_d_cd.es["mf_mst_udc"] = [grep(True, E(g_ud_d_cd)$con_id %in% E(minimum.spanning.tree(g_ud_d_cd, weights=E(g_ud_d_cd)$mf_inv_u))$con_id)]$mf_mst_udc = 1 + + + ### Minimum spanning tree weighted by cost distance + g_ud_d_cd.es["cd_mst_udc"] = 0 + g_ud_d_cd.es["cd_mst_udc"] = [grep(True, E(g_ud_d_cd)$con_id %in% E(minimum.spanning.tree(g_ud_d_cd, weights=E(g_ud_d_cd)$cd_u))$con_id)]$cd_mst_udc = 1 + + ### Minimum spanning tree weighted by competing potential flow + g_ud_d_cd.es["cf_mst_udc"] = 0 + g_ud_d_cd.es["cf_mst_udc"] = [grep(True, E(g_ud_d_cd)$con_id %in% E(minimum.spanning.tree(g_ud_d_cd, weights=E(g_ud_d_cd)$cf_inv_u))$con_id)]$cf_mst_udc = 1 + + ######################################################################## + ### Identify bridges for the undirected graph + E(g_ud.es["is_br_u"] = as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud))), function(x) clusters(delete.edges(g_ud, E(g_ud)[as.integer(x)]))$no-cl_no_d))) + + # Identify bridges for the undirected graph with only direct edges + E(g_ud_d.es["is_br_ud"] = as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud_d))), function(x) clusters(delete.edges(g_ud_d, E(g_ud_d)[as.integer(x)]))$no-cl_no_d))) + + ### Identify bridges for the undirected graph with only direct edges shorter cost distance threshold + E(g_ud_d_cd.es["is_br_udc"] = as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud_d_cd))), function(x) clusters(delete.edges(g_ud_d_cd, E(g_ud_d_cd)[as.integer(x)]))$no-cl_no_d_cd))) + + ######################################################################## + + biconnected_components_d = biconnected.components(g_ud) + + # Identify component edges (biconnected components) for the undirected graph with only direct edges + g_ud.es["bc_e_u = 0 + for c in 1:biconnected_components_d$no: + g_ud.es["bc_e_u[unlist(biconnected_components_d$component_edges[c])] = c + + # Identify tree edges (biconnected components) for the undirected graph with only direct edges + g_ud.es["bc_te_u = 0 + for c in 1:biconnected_components_d$no: + g_ud.es["bc_te_u[unlist(biconnected_components_d$tree_edges[c])] = c + + biconnected_components_d = biconnected.components(g_ud_d) + # Identify component edges (biconnected components) for the undirected graph with only direct edges + g_ud_d.es["bc_e_ud = 0 + for c in 1:biconnected_components_d$no: + g_ud_d.es["bc_e_ud[unlist(biconnected_components_d$component_edges[c])] = c + + # Identify tree edges (biconnected components) for the undirected graph with only direct edges + g_ud_d.es["bc_te_ud = 0 + for c in 1:biconnected_components_d$no: + g_ud_d.es["bc_te_ud[unlist(biconnected_components_d$tree_edges[c])] = c + + biconnected_components_d_cd = biconnected.components(g_ud_d_cd) + # Identify component edges (biconnected components) for the undirected graph with only direct edges shorter cost distance threshold + E(g_ud_d_cd)$bc_e_udc = 0 + for c in 1:biconnected_components_d_cd$no: + E(g_ud_d_cd)$bc_e_udc[unlist(biconnected_components_d_cd$component_edges[c])] = c + + # Identify tree edges (biconnected components) for the undirected graph with only direct edges shorter cost distance threshold + E(g_ud_d_cd)$bc_te_udc = 0 + for c in 1:biconnected_components_d_cd$no: + E(g_ud_d_cd)$bc_te_udc[unlist(biconnected_components_d_cd$tree_edges[c])] = c + + ######################################################################## + ### Calculate edge betweenness + ## Results of the internal edge.betweeness.estimate function are slightly different from the one used here (always larger values): + ## see: g_ud_d.es["cd_leb_ud = edge.betweenness.estimate(g_ud_d, e=E(g_ud_d), directed=False, lnbh_cutoff*connectivity_cutoff, weights=g_ud_d.es["cd_u) + ## But since the internal edge.betweeness.estimate function is not reasonable to use with other weights than cost distance, the workarounds are used nevertheless + ## Further investigation necessary!!! + + ### Calculate edge betweenness on the entire undirected graph with only direct edges + + path_lengths = shortest.paths(g_ud_d, v=V(g_ud_d), to=V(g_ud_d), weights=g_ud_d.es["cd_u) + + # weighted by cost distance + E(g_ud_d.es.["cd_eb_ud"] = edge.betweenness(g_ud_d, e=E(g_ud_d), directed=False, weights=g_ud_d.es["cd_u) + E(g_ud_d.es.["cd_leb_ud"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d)), function(x) get.shortest.paths(g_ud_d, x, V(g_ud_d)[grep(True, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=g_ud_d.es["cd_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d))), plot=False)$counts/2 + # g_ud_d.es["cd_leb_ud = edge.betweenness.estimate(g_ud_d, e=E(g_ud_d), directed=False, lnbh_cutoff*connectivity_cutoff, weights=g_ud_d.es["cd_u) + + # weighted by maximum potential flow + E(g_ud_d.es.["mf_eb_ud"] = edge.betweenness(g_ud_d, e=E(g_ud_d), directed=False, weights=g_ud_d.es["mf_inv_u) + E(g_ud_d.es.["mf_leb_ud"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d)), function(x) get.shortest.paths(g_ud_d, x, V(g_ud_d)[grep(True, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=g_ud_d.es["mf_inv_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d))), plot=False)$counts/2 + + # weighted by competing potential flow + E(g_ud_d.es.["cf_eb_ud"] = edge.betweenness(g_ud_d, e=E(g_ud_d), directed=False, weights=g_ud_d.es["cf_inv_u) + E(g_ud_d.es.["cf_leb_ud"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d)), function(x) get.shortest.paths(g_ud_d, x, V(g_ud_d)[grep(True, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=g_ud_d.es["cf_inv_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d))), plot=False)$counts/2 + + ######################################################################## + ### Calculate edge betweenness on the entire undirected graph with only direct edges shorter cost distance threshold + + path_lengths_cd = shortest.paths(g_ud_d_cd, v=V(g_ud_d_cd), to=V(g_ud_d_cd), weights=E(g_ud_d_cd)$cd_u) + + # weighted by cost distance + g_ud_d_cd.es["cd_eb_udc"] = edge.betweenness(g_ud_d_cd, e=E(g_ud_d_cd), directed=False, weights=E(g_ud_d_cd)$cd_u) + g_ud_d_cd.es["cd_leb_udc"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d_cd)), function(x) get.shortest.paths(g_ud_d_cd, x, V(g_ud_d_cd)[grep(True, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)$cd_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d_cd))), plot=False)$counts/2 + + + # weighted by maximum potential flow + E(g_ud_d_cd.es["mf_eb_udc"] = edge.betweenness(g_ud_d_cd, e=E(g_ud_d_cd), directed=False, weights=E(g_ud_d_cd)$mf_inv_u) + E(g_ud_d_cd.es["mf_leb_udc"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d_cd)), function(x) get.shortest.paths(g_ud_d_cd, x, V(g_ud_d_cd)[grep(True, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)$mf_inv_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d_cd))), plot=False)$counts/2 + + # weighted by competing potential flow + E(g_ud_d_cd.es["cf_eb_udc"] = edge.betweenness(g_ud_d_cd, e=E(g_ud_d_cd), directed=False, weights=E(g_ud_d_cd)$cf_inv_u) + E(g_ud_d_cd.es["cf_leb_udc"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d_cd)), function(x) get.shortest.paths(g_ud_d_cd, x, V(g_ud_d_cd)[grep(True, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)$cf_inv_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d_cd))), plot=False)$counts/2 + + if cl_thresh > 0: + ######################################################################### + ####### Calculate edge betweenness community on the undirected graph + ## weighted by competing potential flow + ebc = edge.betweenness.community(g_ud, weights=g_ud.es["cf_inv_u, directed=False, edge.betweenness=True, merges=True, bridges=True, modularity=True, membership=True) + + E(g_ud.es["cf_iebc_v"] = ebc$edge.betweenness + # cf_iebc_m = ebc$merges + E(g_ud.es["cf_iebc_r"] = ebc$removed.edges + E(g_ud.es["cf_iebc_b"] = unlist(mclapply(mc.cores={0}, 1:length(E(g_ud)), function(x) ifelse(x %in% ebc$bridges, 1, 0))) + E(g_ud.es["cf_iebc_c"] = crossing(ebc, g_ud) + + # g.vs["cf_iebc_mo = ebc$modularity + V(g.vs["cf_iebc_me"] = ebc$membership + + com_struct = as.character(cutat(ebc, cl_no_d)) + for i in (cl_no_d+1):(cl_no_d+cl_thresh): + com_struct = paste(com_struct, cutat(ebc, i), sep=';') + + V(g.vs["cf_iebc_cs"] = as.character(com_struct) + V(g.vs["cf_iebc_cl"] = cutat(ebc, (cl_no_d+cl_thresh)) + + cf_iebc_u_mod = modularity(ebc) + com_no_u = length(ebc) + com_sizes_u = sizes(ebc) + com_sizes_u_names = paste("Size of comumity ", names(sizes(ebc)), " (at maximum modularity score (from iebc))", sep="") + + # The following workaround was written when igraph did not support weights when calculating edge.betweenness.community + # The results of this procedure differ significantly from the (new) edge.betweenness.community in igraph 0.6-2 with weights (above) + # In a first visual the results of the following procedure look more reasonable, but further investigation is necessary!!! + + n = 1 + ebc_con_id_u = 0 + ebc_rank = 0 + ebc_value = 0 + ebc_clust = 0 + g_ud.vs["cf_ebc_cs = as.character(clusters(g_ud)$membership) + del_eb_g = g_ud + + for edge in del_eb_g.es: + # Recalculate edge betweenness for the remaining edges + E(del_eb_g)$cf_eb_ud = edge.betweenness(del_eb_g, e=E(del_eb_g), directed=False, weights=E(del_eb_g)$cf_inv_u) + # Identify edge with highest edge betweenness value + idx = sort.int(E(del_eb_g)$cf_eb_ud, decreasing=True, na.last=NA, index.return=True)$ix[1] + # Save edge betweenness value for this edge + ebc_con_id_u[n] = E(del_eb_g)$con_id_u[idx] + ebc_rank[n] = n + ebc_clust[n] = clusters(del_eb_g)$no + if((ebc_clust[n]-cl_no_d)<(cl_no_d+cl_thresh)) {{if(n>1) {{if(ebc_clust[n]>ebc_clust[n-1]) g_ud.vs["cf_ebc_cs = paste(g_ud.vs["cf_ebc_cs, clusters(del_eb_g)$membership, sep=";")}}}} + ebc_value[n] = E(del_eb_g)$cf_eb_ud[idx] + # Delete edge with highest edge betweenness value + del_eb_g = delete.edges(del_eb_g, E(del_eb_g)[idx]) + n = n+1 + + g_ud.vs["cf_ebc_cl = unlist(mclapply(mc.cores={0}, 1:length(V(g_ud)), function(x) strsplit(g_ud.vs["cf_ebc_cs, ";")[[x]][(cl_no_d+cl_thresh)])) + + # ebc_df = merge(data.frame(con_id=g.es["con_id, con_id_u=g.es["con_id_u), data.frame(con_id_u=ebc_con_id_u, ebc_rank=ebc_rank, ebc_clust=ebc_clust, ebc_value=ebc_value), all.x=True) + idx = sort.int(ebc_con_id_u, index.return=True)$ix + g_ud.es["cf_ebc_v = ebc_value[idx] + g_ud.es["cf_ebc_r = ebc_rank[idx] + g_ud.es["cf_ebc_c = ebc_clust[idx] + + g_ud.es["cf_ebc_vi = strftime(Sys.time()+g_ud.es["cf_ebc_r) + + ######################################################################### + ### Calculate community connectors (based on community structure from ebc) + com_membership = data.frame(patch_id=g_ud_d_cd.vs["patch_id, com=g_ud.vs["cf_ebc_cl) + com_pc_pre = merge(merge(data.frame(id=1:length(E(g)), from_p=g.es["from_p, to_p=g.es["to_p), com_membership, by.x="from_p", by.y="patch_id"), com_membership, by.x="to_p", by.y="patch_id")[order(merge(merge(data.frame(id=1:length(E(g)), from_p=g.es["from_p, to_p=g.es["to_p), com_membership, by.x="from_p", by.y="patch_id"), com_membership, by.x="to_p", by.y="patch_id")$id) ,] + g.es["cf_ebc_cc = ifelse(com_pc_pre$com.x==com_pc_pre$com.y,0,1) -###Create initial export data frame - -vertex_export_df_ud_d_cd <- as.data.frame(1:length(V(g_ud_d_cd))) - -####Adjust vertex-attributes first - -g_ud_d_cd <- remove.vertex.attribute(g_ud_d_cd, "pop_proxy") - -vertex_attribute_list <- list.vertex.attributes(g_ud_d_cd) - -if(length(vertex_attribute_list)>0) {{ - for (vl in vertex_attribute_list) {{ - vertex_export_df_ud_d_cd <- as.data.frame(cbind(vertex_export_df_ud_d_cd, get.vertex.attribute(g_ud_d_cd, vl))) - }} - vertex_export_df_ud_d_cd <- vertex_export_df_ud_d_cd[2:length(vertex_export_df_ud_d_cd)] - names(vertex_export_df_ud_d_cd) <- vertex_attribute_list -}} -###Create initial export data frame - -vertex_export_df_d <- as.data.frame(1:length(V(g))) - -vertex_attribute_list <- list.vertex.attributes(g) -if(length(vertex_attribute_list)>0) {{ - - for (vl in vertex_attribute_list) {{ - vertex_export_df_d <- as.data.frame(cbind(vertex_export_df_d, get.vertex.attribute(g, vl))) - }} - vertex_export_df_d <- vertex_export_df_d[2:length(vertex_export_df_d)] - names(vertex_export_df_d) <- vertex_attribute_list -}} - -vertex_export_df <- merge(merge(merge(vertex_export_df_d, vertex_export_df_ud, by="patch_id"), vertex_export_df_ud_d, by="patch_id"), vertex_export_df_ud_d_cd, by="patch_id") - -dbWriteTable(con, vertex_output, vertex_export_df, overwrite=overwrite, row.names=FALSE) - -######################################################################### - - -###Create initial export data frame for the undirected graph - -edge_export_df_ud <- as.data.frame(E(g_ud)$con_id_u) - -####Adjust edge attributes for the undirected graphs first -g_ud <- remove.edge.attribute(g_ud, "con_id") -g_ud <- remove.edge.attribute(g_ud, "con_id_u") -g_ud <- remove.edge.attribute(g_ud, "from_p") -g_ud <- remove.edge.attribute(g_ud, "from_pop") -g_ud <- remove.edge.attribute(g_ud, "to_p") -g_ud <- remove.edge.attribute(g_ud, "to_pop") -g_ud <- remove.edge.attribute(g_ud, "cost_distance") -g_ud <- remove.edge.attribute(g_ud, "cd_u") -g_ud <- remove.edge.attribute(g_ud, "distance_weight_e") -#g_ud <- remove.edge.attribute(g_ud, "distance_weight_e_ud") -g_ud <- remove.edge.attribute(g_ud, "mf_o") -g_ud <- remove.edge.attribute(g_ud, "mf_i") -g_ud <- remove.edge.attribute(g_ud, "mf_o_inv") -g_ud <- remove.edge.attribute(g_ud, "mf_i_inv") -g_ud <- remove.edge.attribute(g_ud, "mf_u") -g_ud <- remove.edge.attribute(g_ud, "mf_inv_u") -g_ud <- remove.edge.attribute(g_ud, "cf") -g_ud <- remove.edge.attribute(g_ud, "cf_inv") -g_ud <- remove.edge.attribute(g_ud, "cf_u") -g_ud <- remove.edge.attribute(g_ud, "cf_inv_u") - -edge_attribute_list <- list.edge.attributes(g_ud) - -if(length(edge_attribute_list )>0) {{ - for (el in edge_attribute_list) {{ - edge_export_df_ud <- as.data.frame(cbind(edge_export_df_ud, get.edge.attribute(g_ud, el))) - }} - names(edge_export_df_ud) <- append("con_id_u", edge_attribute_list) -}} - -###Create initial export data frame for the undirected graph with only direct edges - -edge_export_df_ud_d <- as.data.frame(E(g_ud_d)$con_id_u) - -####Adjust edge attributes for the undirected graph with only direct edges -g_ud_d <- remove.edge.attribute(g_ud_d, "con_id") -g_ud_d <- remove.edge.attribute(g_ud_d, "con_id_u") -g_ud_d <- remove.edge.attribute(g_ud_d, "from_p") -g_ud_d <- remove.edge.attribute(g_ud_d, "from_pop") -g_ud_d <- remove.edge.attribute(g_ud_d, "to_p") -g_ud_d <- remove.edge.attribute(g_ud_d, "to_pop") -g_ud_d <- remove.edge.attribute(g_ud_d, "cost_distance") -g_ud_d <- remove.edge.attribute(g_ud_d, "cd_u") -g_ud_d <- remove.edge.attribute(g_ud_d, "distance_weight_e") -#g_ud_d <- remove.edge.attribute(g_ud_d, "distance_weight_e_ud") -g_ud_d <- remove.edge.attribute(g_ud_d, "mf_o") -g_ud_d <- remove.edge.attribute(g_ud_d, "mf_i") -g_ud_d <- remove.edge.attribute(g_ud_d, "mf_o_inv") -g_ud_d <- remove.edge.attribute(g_ud_d, "mf_i_inv") -g_ud_d <- remove.edge.attribute(g_ud_d, "mf_u") -g_ud_d <- remove.edge.attribute(g_ud_d, "mf_inv_u") -g_ud_d <- remove.edge.attribute(g_ud_d, "cf") -g_ud_d <- remove.edge.attribute(g_ud_d, "cf_inv") -g_ud_d <- remove.edge.attribute(g_ud_d, "cf_u") -g_ud_d <- remove.edge.attribute(g_ud_d, "cf_inv_u") -g_ud_d <- remove.edge.attribute(g_ud_d, "isshort") -g_ud_d <- remove.edge.attribute(g_ud_d, "isshort_cd") -g_ud_d <- remove.edge.attribute(g_ud_d, "isshort_mf") -g_ud_d <- remove.edge.attribute(g_ud_d, "isshort_cf") - -edge_attribute_list <- list.edge.attributes(g_ud_d) - -if(length(edge_attribute_list )>0) {{ - for (el in edge_attribute_list) {{ - edge_export_df_ud_d <- as.data.frame(cbind(edge_export_df_ud_d, get.edge.attribute(g_ud_d, el))) - }} - names(edge_export_df_ud_d) <- append("con_id_u", edge_attribute_list) -}} -###Create initial export data frame for the undirected graph with only direct edges - -edge_export_df_ud_d_cd <- as.data.frame(E(g_ud_d_cd)$con_id_u) - -####Adjust edge attributes for the undirected graph with only direct edges -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "con_id") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "con_id_u") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "from_p") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "from_pop") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "to_p") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "to_pop") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "cost_distance") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "cd_u") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "distance_weight_e") -#g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "distance_weight_e_ud") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "mf_o") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "mf_i") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "mf_o_inv") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "mf_i_inv") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "mf_u") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "mf_inv_u") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "cf") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "cf_inv") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "cf_u") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "cf_inv_u") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "isshort") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "isshort_cd") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "isshort_mf") -g_ud_d_cd <- remove.edge.attribute(g_ud_d_cd, "isshort_cf") - - -edge_attribute_list <- list.edge.attributes(g_ud_d_cd) - -if(length(edge_attribute_list )>0) {{ - for (el in edge_attribute_list) {{ - edge_export_df_ud_d_cd <- as.data.frame(cbind(edge_export_df_ud_d_cd, get.edge.attribute(g_ud_d_cd, el))) - }} - names(edge_export_df_ud_d_cd) <- append("con_id_u", edge_attribute_list) -}} - - -###Create initial export data frame for the directed graph -edge_export_df <- as.data.frame(1:length(E(g))) - -####Adjust edge attributes for the directed graph first - -E(g)$cd <- E(g)$cost_distance -E(g)$distk <- E(g)$distance_weight_e -E(g)$distk_u <- E(g)$distance_weight_e_ud -g <- remove.edge.attribute(g, "distance_weight_e") -g <- remove.edge.attribute(g, "distance_weight_e_ud") -g <- remove.edge.attribute(g, "cost_distance") - -edge_attribute_list <- list.edge.attributes(g) - -if(length(edge_attribute_list )>0) {{ - for (el in edge_attribute_list) {{ - edge_export_df <- as.data.frame(cbind(edge_export_df, get.edge.attribute(g, el))) - }} - names(edge_export_df) <- append("id", edge_attribute_list) -}} - -export_df_list <-c("edge_export_df_ud", "edge_export_df_ud_d", "edge_export_df_ud_d_cd") -for (df in export_df_list) {{ - if(length(names(get(df)))) {{ - edge_export_df_final <- merge(edge_export_df, get(df), all.x=TRUE, by="con_id_u", suffixes=c("_x", "_y")) - edge_export_df <- edge_export_df_final + else: + gscript.verbose("Skipping comunity algorithms...") + + ######################################################################## + ###### Calculate cluster membership + g_ud_d.vs["cl_ud = clusters(g_ud_d)$membership + g_ud_d_cd.vs["cl_udc = clusters(g_ud_d_cd)$membership + + ######################################################################### + ### Calculate potential cluster connectors (based on cost distance threshold) + cl_membership = data.frame(patch_id=g_ud_d_cd.vs["patch_id, cl=g_ud_d_cd.vs["cl_udc) + cl_pc_pre = merge(merge(data.frame(id=1:length(E(g)), from_p=g.es["from_p, to_p=g.es["to_p), cl_membership, by.x="from_p", by.y="patch_id"), cl_membership, by.x="to_p", by.y="patch_id")[order(merge(merge(data.frame(id=1:length(E(g)), from_p=g.es["from_p, to_p=g.es["to_p), cl_membership, by.x="from_p", by.y="patch_id"), cl_membership, by.x="to_p", by.y="patch_id")$id) ,] + g.es["cl_pc = ifelse(cl_pc_pre$cl.x==cl_pc_pre$cl.y,0,1) + + gscript.verbose("Starting analysis on vertex level...") + + ######################################################################## + ### Analysis on vertex level + ######################################################################## + + ######################################################################## + ###### Calculate degree centrality + # On the directed graph + # g.vs["deg_dir = as.integer(degree(g, v=V(g), mode=c("in"))) + + # On the entire undirected graph + g_ud.vs["deg_u = degree(g_ud, v=V(g_ud)) + + # On the undirected graph with edges shorter cost distance threshold + g_ud_cd.vs["deg_uc = degree(g_ud_cd, v=V(g_ud_cd)) + + # On the undirected graph with only direct edges + g_ud_d.vs["deg_ud = degree(g_ud_d, v=V(g_ud_d)) + + # On the undirected graph with only direct edges shorter cost distance threshold + g_ud_d_cd.vs["deg_udc = degree(g_ud_d_cd, v=V(g_ud_d_cd)) + + ######################################################################## + ###### Calculate closeness centrality + g_ud_d.vs["cd_cl_ud = closeness(g_ud_d, vids=V(g_ud_d), weights = g_ud_d.es["cd_u, normalized = False) + g_ud_d.vs["mf_cl_ud = closeness(g_ud_d, vids=V(g_ud_d), weights = g_ud_d.es["mf_inv_u, normalized = False) + g_ud_d.vs["cf_cl_ud = closeness(g_ud_d, vids=V(g_ud_d), weights = g_ud_d.es["cf_inv_u, normalized = False) + + g_ud_d_cd.vs["cd_cl_udc = closeness(g_ud_d_cd, vids=V(g_ud_d_cd), weights = E(g_ud_d_cd)$cd_u, normalized = False) + g_ud_d_cd.vs["mf_cl_udc = closeness(g_ud_d_cd, vids=V(g_ud_d_cd), weights = E(g_ud_d_cd)$mf_inv_u, normalized = False) + g_ud_d_cd.vs["cf_cl_udc = closeness(g_ud_d_cd, vids=V(g_ud_d_cd), weights = E(g_ud_d_cd)$cf_inv_u, normalized = False) + + ######################################################################## + ###### Calculate biconnected components + + # On the undirected graph with only direct edges + # Number of new clusters when a vertex is deleted + g_ud_d.vs["art_ud = unlist(mclapply(mc.cores={0}, 1:(length(V(g_ud_d))), function(x) ifelse((clusters(delete.vertices(g_ud_d, V(g_ud_d)[x]))$no-cl_no_d)<0,0,clusters(delete.vertices(g_ud_d, V(g_ud_d)[x]))$no-cl_no_d))) + # Vertex is articulation point + g_ud_d.vs["art_p_ud = 0 + g_ud_d.vs["art_p_ud[biconnected_components_d$articulation_points] = 1 + + # On the undirected graph with only direct edges shorter cost distance threshold + # Number of new clusters when a vertex is deleted + g_ud_d_cd.vs["art_udc = unlist(mclapply(mc.cores={0}, 1:(length(V(g_ud_d_cd))), function(x) ifelse((clusters(delete.vertices(g_ud_d_cd, V(g_ud_d_cd)[x]))$no-cl_no_d_cd)<0,0,clusters(delete.vertices(g_ud_d_cd, V(g_ud_d_cd)[x]))$no-cl_no_d_cd))) + # Vertex is articulation point + g_ud_d_cd.vs["art_p_udc = 0 + g_ud_d_cd.vs["art_p_udc[biconnected_components_d_cd$articulation_points] = 1 + + ######################################################################## + ###### Calculate eigenvector centrality + # Calculate eigenvector for the entire directed graph + + # Calculate eigenvector centrality weighted by competing potential flow + cf_evc_groups = groupedData(cf ~ 1 | to_p, data=e, FUN=sum) + cf_evc = gsummary(cf_evc_groups, sum) + df_cf_evc_d = data.frame(patch_id=cf_evc$to_p, cf_evc_d=cf_evc$cf) + + # Calculate eigenvector centrality weighted by maximum potential flow + mf_o_evc_groups = groupedData(mf_o ~ 1 | to_p, data=e, FUN=sum) + mf_o_evc = gsummary(mf_o_evc_groups, sum) + df_mf_i_evc_d = data.frame(patch_id=mf_o_evc$to_p, mf_evc_d=mf_o_evc$mf_o) + + df_evc_d = merge(df_cf_evc_d, df_mf_i_evc_d) + + ###### Calculate eigenvector centrality on the directed graph with edges shorter cost distance threshold + # Calculate eigenvector centrality weighted by competing potential flow + evc_udc_pre = data.frame(to_p=e$to_p[grep(True, e$cd_u', - ' ', - ' ', - ' ', - ' ') -colortable_bin <- c(' ') -for (attribute in names(edge_export_df_final)) {{ - -#Skip id and geom columns -if(attribute %in% c("id", "con_id", "con_id_u", "from_p", "to_p", "WKT", "cf_ebc_vi")) {{ next }} - -st_mod <- storage.mode(unlist(edge_export_df_final[grep(1, match(names(edge_export_df_final), attribute))])) -att_val <- unlist(edge_export_df_final[grep(1, match(names(edge_export_df_final), attribute))]) - -#Write header -qml <- c("") -qml <- append(qml, '') -qml <- append(qml, ' 255') - -if((max(att_val, na.rm=TRUE)-min(att_val, na.rm=TRUE)==1)) {{ - -#More header -qml <- append(qml, paste(' ', sep='')) - -#Categories -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') - -#Write symbols -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, colortable_bin) -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -}} -else {{ - -attribute_quantile <- quantile(edge_export_df_final[grep(1, match(names(edge_export_df_final), attribute))], probs=seq(0, 1, 1/5), na.rm=TRUE, type=8) - -#Write more header -qml <- append(qml, paste(' ', sep='')) - -#Write ranges -ranges <- character() -qml <- append(qml, ' ') -for (quant in 1:no_quantile) {{ -ranges <- append(ranges, paste(' ', sep='')) -}} -qml <- append(qml, ranges) -qml <- append(qml, ' ') - -#Write symbols -qml <- append(qml, ' ') - -for (quant in 1:no_quantile) {{ -qml <- append(qml, paste(' ', - sep='')) -qml <- append(qml, paste(' ', - sep='')) -qml <- append(qml, ' ') -qml <- append(qml, colortable[quant]) -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -}} -}} -#Write Footer -qml <- append(qml, ' ') - - -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, ' ') -qml <- append(qml, paste(' "', attribute, '"', sep='')) -qml <- append(qml, ' ') -qml <- append(qml, '') - -#Save qml-file -qml_output <- paste(qml_directory, paste(paste("edge_measures_", attribute, sep=''), "qml", sep='.'), sep="/") -con_qml <- file(qml_output, open="wt") -write.table(qml, con_qml, append = FALSE, quote = FALSE, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = FALSE, col.names = FALSE) -close(con_qml) -}} -}} -######################################################################## -#Close R -######################################################################## -#q() -""".format( - cores + ############################################## + # Export network overview measures + grml = rbind(c("Command", command), + c("Number of vertices", vertices_n), + c("Number of edges (undirected)", edges_n), + c("Number of direct edges (undirected)", edges_n_d), + c("Number of edges shorter than cost distance threshold (undirected)", edges_n_cd), + c("Number of direct edges shorter than cost distance threshold (undirected)", edges_n_d_cd), + c("Number of clusters of the entire graph", cl_no_d), + c("Number of clusters of the graph with only edges shorter cost distance threshold", cl_no_d_cd), + c("Size of the largest cluster of the entire graph", cl_max_size_d), + c("Size of the largest cluster of the graph with only edges shorter cost distance threshold", cl_max_size_d_cd), + c("Average size of the clusters of the entire graph", cl_mean_size_d), + c("Average size of the clusters of the graph with only edges shorter cost distance threshold", cl_mean_size_d_cd), + c("Diameter of the entire graph (undirected)", diam), + c("Diameter of the graph with only direct edges (undirected)", diam_d), + c("Diameter of the graph with only edges shorter cost distance threshold", diam_d_cd), + c("Density of the entire graph (directed)", density), + c("Density of the entire graph (undirected)", density_u), + c("Density of the graph with only direct edges (undirected)", density_ud)) + + if cl_thresh > 0: + grml = rbind(grml, + c("Density of the graph with only edges shorter cost distance threshold", density_udc), + c("Modularity (from iebc) of the entire graph (undirected) weighted by cf", cf_iebc_u_mod), + c("Number of communities (at maximum modularity score (from iebc)) of the entire (undirected) graph weighted by cf", com_no_u), + c("com_sizes_u", toString(com_sizes_u, sep=',')), + c("com_sizes_u_names", toString(com_sizes_u_names, sep=',')) ) - if cores <= 1 or os_type == "Windows": - rscript.replace("mclapply(mc.cores={}, ".format(cores), "lapply(") - - with open(os.path.join(qml_style_dir, "net_r_script.r"), "w") as rs_file: - rs_file.write(rscript) + network_overview_measures = data.frame( + measure=grml[,1], + value=grml[,2]) - robjects.r(rscript) + # Write dataframe + dbWriteTable(con, network_output, network_overview_measures, overwrite=overwrite, row.names=False) + + ############################################################### + ##### Merge vertexmeasures in a dataframe and save it + + ### Create initial export data frame + vertex_export_df_ud = as.data.frame(1:length(V(g_ud))) + + #### Adjust vertex-attributes first + + g_ud = remove.vertex.attribute(g_ud, "pop_proxy") + + vertex_attribute_list = list.vertex.attributes(g_ud) + if len(vertex_attribute_list) > 0: + for vl in vertex_attribute_list: + vertex_export_df_ud = as.data.frame(cbind(vertex_export_df_ud, get.vertex.attribute(g_ud, vl))) + + vertex_export_df_ud = vertex_export_df_ud[2:length(vertex_export_df_ud)] + names(vertex_export_df_ud) = vertex_attribute_list + + ### Create initial export data frame + + vertex_export_df_ud_d = as.data.frame(1:length(V(g_ud_d))) + + #### Adjust vertex-attributes first + g_ud_d = remove.vertex.attribute(g_ud_d, "pop_proxy") + + vertex_attribute_list = list.vertex.attributes(g_ud_d) + if length(vertex_attribute_list) > 0: + for vl in vertex_attribute_list: + vertex_export_df_ud_d = as.data.frame(cbind(vertex_export_df_ud_d, get.vertex.attribute(g_ud_d, vl))) + + vertex_export_df_ud_d = vertex_export_df_ud_d[2:length(vertex_export_df_ud_d)] + names(vertex_export_df_ud_d) = vertex_attribute_list + + ###Create initial export data frame + vertex_export_df_ud_d_cd = as.data.frame(1:length(V(g_ud_d_cd))) + + ####Adjust vertex-attributes first + g_ud_d_cd = remove.vertex.attribute(g_ud_d_cd, "pop_proxy") + vertex_attribute_list = list.vertex.attributes(g_ud_d_cd) + + if len(vertex_attribute_list)>0: + for vl in vertex_attribute_list: + vertex_export_df_ud_d_cd = as.data.frame(cbind(vertex_export_df_ud_d_cd, get.vertex.attribute(g_ud_d_cd, vl))) + + vertex_export_df_ud_d_cd = vertex_export_df_ud_d_cd[2:length(vertex_export_df_ud_d_cd)] + names(vertex_export_df_ud_d_cd) = vertex_attribute_list + + ###Create initial export data frame + vertex_export_df_d = as.data.frame(1:length(V(g))) + vertex_attribute_list = list.vertex.attributes(g) + if len(vertex_attribute_list)>0: + for vl in vertex_attribute_list: + vertex_export_df_d = as.data.frame(cbind(vertex_export_df_d, get.vertex.attribute(g, vl))) + + vertex_export_df_d = vertex_export_df_d[2:length(vertex_export_df_d)] + names(vertex_export_df_d) = vertex_attribute_list + + vertex_export_df = merge(merge(merge(vertex_export_df_d, vertex_export_df_ud, by="patch_id"), vertex_export_df_ud_d, by="patch_id"), vertex_export_df_ud_d_cd, by="patch_id") + dbWriteTable(con, vertex_output, vertex_export_df, overwrite=overwrite, row.names=False) + + ######################################################################### + ###Create initial export data frame for the undirected graph + edge_export_df_ud = as.data.frame(g_ud.es["con_id_u) + + ####Adjust edge attributes for the undirected graphs first + g_ud = remove.edge.attribute(g_ud, "con_id") + g_ud = remove.edge.attribute(g_ud, "con_id_u") + g_ud = remove.edge.attribute(g_ud, "from_p") + g_ud = remove.edge.attribute(g_ud, "from_pop") + g_ud = remove.edge.attribute(g_ud, "to_p") + g_ud = remove.edge.attribute(g_ud, "to_pop") + g_ud = remove.edge.attribute(g_ud, "cost_distance") + g_ud = remove.edge.attribute(g_ud, "cd_u") + g_ud = remove.edge.attribute(g_ud, "distance_weight_e") + #g_ud = remove.edge.attribute(g_ud, "distance_weight_e_ud") + g_ud = remove.edge.attribute(g_ud, "mf_o") + g_ud = remove.edge.attribute(g_ud, "mf_i") + g_ud = remove.edge.attribute(g_ud, "mf_o_inv") + g_ud = remove.edge.attribute(g_ud, "mf_i_inv") + g_ud = remove.edge.attribute(g_ud, "mf_u") + g_ud = remove.edge.attribute(g_ud, "mf_inv_u") + g_ud = remove.edge.attribute(g_ud, "cf") + g_ud = remove.edge.attribute(g_ud, "cf_inv") + g_ud = remove.edge.attribute(g_ud, "cf_u") + g_ud = remove.edge.attribute(g_ud, "cf_inv_u") + + edge_attribute_list = list.edge.attributes(g_ud) + + if length(edge_attribute_list )>0: + for el in edge_attribute_list: + edge_export_df_ud = as.data.frame(cbind(edge_export_df_ud, get.edge.attribute(g_ud, el))) + names(edge_export_df_ud) = append("con_id_u", edge_attribute_list) + + ###Create initial export data frame for the undirected graph with only direct edges + edge_export_df_ud_d = as.data.frame(g_ud_d.es["con_id_u) + + ####Adjust edge attributes for the undirected graph with only direct edges + g_ud_d = remove.edge.attribute(g_ud_d, "con_id") + g_ud_d = remove.edge.attribute(g_ud_d, "con_id_u") + g_ud_d = remove.edge.attribute(g_ud_d, "from_p") + g_ud_d = remove.edge.attribute(g_ud_d, "from_pop") + g_ud_d = remove.edge.attribute(g_ud_d, "to_p") + g_ud_d = remove.edge.attribute(g_ud_d, "to_pop") + g_ud_d = remove.edge.attribute(g_ud_d, "cost_distance") + g_ud_d = remove.edge.attribute(g_ud_d, "cd_u") + g_ud_d = remove.edge.attribute(g_ud_d, "distance_weight_e") + #g_ud_d = remove.edge.attribute(g_ud_d, "distance_weight_e_ud") + g_ud_d = remove.edge.attribute(g_ud_d, "mf_o") + g_ud_d = remove.edge.attribute(g_ud_d, "mf_i") + g_ud_d = remove.edge.attribute(g_ud_d, "mf_o_inv") + g_ud_d = remove.edge.attribute(g_ud_d, "mf_i_inv") + g_ud_d = remove.edge.attribute(g_ud_d, "mf_u") + g_ud_d = remove.edge.attribute(g_ud_d, "mf_inv_u") + g_ud_d = remove.edge.attribute(g_ud_d, "cf") + g_ud_d = remove.edge.attribute(g_ud_d, "cf_inv") + g_ud_d = remove.edge.attribute(g_ud_d, "cf_u") + g_ud_d = remove.edge.attribute(g_ud_d, "cf_inv_u") + g_ud_d = remove.edge.attribute(g_ud_d, "isshort") + g_ud_d = remove.edge.attribute(g_ud_d, "isshort_cd") + g_ud_d = remove.edge.attribute(g_ud_d, "isshort_mf") + g_ud_d = remove.edge.attribute(g_ud_d, "isshort_cf") + + edge_attribute_list = list.edge.attributes(g_ud_d) + + if len(edge_attribute_list )>0: + for el in edge_attribute_list: + edge_export_df_ud_d = as.data.frame(cbind(edge_export_df_ud_d, get.edge.attribute(g_ud_d, el))) + + names(edge_export_df_ud_d) = append("con_id_u", edge_attribute_list) + + ###Create initial export data frame for the undirected graph with only direct edges + edge_export_df_ud_d_cd = as.data.frame(E(g_ud_d_cd)$con_id_u) + + ####Adjust edge attributes for the undirected graph with only direct edges + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "con_id") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "con_id_u") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "from_p") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "from_pop") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "to_p") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "to_pop") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "cost_distance") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "cd_u") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "distance_weight_e") + #g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "distance_weight_e_ud") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "mf_o") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "mf_i") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "mf_o_inv") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "mf_i_inv") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "mf_u") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "mf_inv_u") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "cf") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "cf_inv") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "cf_u") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "cf_inv_u") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "isshort") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "isshort_cd") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "isshort_mf") + g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "isshort_cf") + + edge_attribute_list = list.edge.attributes(g_ud_d_cd) + + if len(edge_attribute_list )>0: + for el in edge_attribute_list: + edge_export_df_ud_d_cd = as.data.frame(cbind(edge_export_df_ud_d_cd, get.edge.attribute(g_ud_d_cd, el))) + names(edge_export_df_ud_d_cd) = append("con_id_u", edge_attribute_list) + + ###Create initial export data frame for the directed graph + edge_export_df = as.data.frame(1:length(E(g))) + + ####Adjust edge attributes for the directed graph first + g.es["cd"] = g.es["cost_distance + g.es["distk"] = g.es["distance_weight_e + g.es["distk_u"] = g.es["distance_weight_e_ud + g = remove.edge.attribute(g, "distance_weight_e") + g = remove.edge.attribute(g, "distance_weight_e_ud") + g = remove.edge.attribute(g, "cost_distance") + + edge_attribute_list = list.edge.attributes(g) + + if len(edge_attribute_list )>0: + for el in edge_attribute_list: + edge_export_df = as.data.frame(cbind(edge_export_df, get.edge.attribute(g, el))) + + names(edge_export_df) = append("id", edge_attribute_list) + + export_df_list =c("edge_export_df_ud", "edge_export_df_ud_d", "edge_export_df_ud_d_cd") + for df in export_df_list: + if len(names(get(df))): + edge_export_df_final = merge(edge_export_df, get(df), all.x=True, by="con_id_u", suffixes=c("_x", "_y")) + edge_export_df = edge_export_df_final + + dbWriteTable(con, edge_output, edge_export_df, overwrite=overwrite, row.names=False) + dbDisconnect(con) + + if qml_directory != '': + ######################### + #CREATE .qml-files for edge measures visualistion in QGIS + + no_quantile = 5 + colortable = c(' ', + ' ', + ' ', + ' ', + ' ') + colortable_bin = c(' ') + for (attribute in names(edge_export_df_final)) {{ + + #Skip id and geom columns + if(attribute %in% c("id", "con_id", "con_id_u", "from_p", "to_p", "WKT", "cf_ebc_vi")) {{ next }} + + st_mod = storage.mode(unlist(edge_export_df_final[grep(1, match(names(edge_export_df_final), attribute))])) + att_val = unlist(edge_export_df_final[grep(1, match(names(edge_export_df_final), attribute))]) + + #Write header + qml = "" + qml += '' + qml += ' 255' + + if max(att_val, na.rm=True)-min(att_val, na.rm=True)==1: + + #More header + qml += paste(' ', sep='')) + + #Categories + qml += ' ') + qml += ' ') + qml += ' ') + + #Write symbols + qml += ' ') + qml += ' ') + qml += ' ') + qml += ' ') + qml += colortable_bin) + qml += ' ') + qml += ' ') + qml += ' ') + qml += ' ') + qml += ' ') + qml += ' ') + qml += ' ') + qml += ' ') + else: - grass.run_command( - "g.copy", quiet=True, vector="{},{}".format(network_map, edge_output) - ) - grass.run_command( - "g.copy", quiet=True, vector="{},{}".format(in_vertices, vertex_output) - ) + attribute_quantile = quantile(edge_export_df_final[grep(1, match(names(edge_export_df_final), attribute))], probs=seq(0, 1, 1/5), na.rm=True, type=8) + + #Write more header + qml += paste(' ', sep='')) + + #Write ranges + ranges = character() + qml += ' ' + for quant in 1:no_quantile: + ranges = append(ranges, paste(' ', sep='') + + qml += ranges + qml += ' ' + + #Write symbols + qml += ' ' + + for quant in 1:no_quantile: + qml += paste(' ', + sep='') + qml += paste(' ', + sep='') + qml += ' ' + qml += colortable[quant] + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + + #Write Footer + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += paste(' "', attribute, '"', sep='') + qml += ' ' + qml += '' + + #Save qml-file + qml_output = paste(qml_directory, paste(paste("edge_measures_", attribute, sep=''), "qml", sep='.'), sep="/") + con_qml = file(qml_output, open="wt") + write.table(qml, con_qml, append = False, quote = False, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = False, col.names = False) + close(con_qml) + + grass.run_command('g.copy', quiet=True, + vector='{},{}'.format(network_map, + edge_output)) + grass.run_command('g.copy', quiet=True, + vector='{},{}'.format(in_vertices, + vertex_output)) # Use v.db.connect instead of v.db.join (much faster) - grass.run_command( - "v.db.join", - map=edge_output, - column="cat", - other_table=edge_output_tmp, - other_column="con_id", - quiet=True, - ) - grass.run_command( - "v.db.join", - map=vertex_output, - column="cat", - other_table=vertex_output_tmp, - other_column="patch_id", - quiet=True, - ) + grass.run_command('v.db.join', map=edge_output, column='cat', + other_table=edge_output_tmp, + other_column='con_id', quiet=True) + grass.run_command('v.db.join', map=vertex_output, column='cat', + other_table=vertex_output_tmp, + other_column='patch_id', quiet=True) - update_history = "{}\n{}".format(net_hist_str, os.environ["CMDLINE"]) + update_history = '{}\n{}'.format(net_hist_str, + os.environ['CMDLINE']) - grass.run_command( - "v.support", - flags="h", - map=vertex_output, - person=os.environ["USER"], - cmdhist=update_history, - ) + grass.run_command('v.support', flags='h', map=vertex_output, + person=os.environ['USER'], + cmdhist=update_history) - grass.run_command( - "v.support", - flags="h", - map=edge_output, - person=os.environ["USER"], - cmdhist=update_history, - ) + grass.run_command('v.support', flags='h', map=edge_output, + person=os.environ['USER'], + cmdhist=update_history) if __name__ == "__main__": options, flags = grass.parser() atexit.register(cleanup) - sys.exit(main()) + sys.exit(main()) \ No newline at end of file From 2c32389e736d1d87d1062fb7043fa189cbeb090b Mon Sep 17 00:00:00 2001 From: ninsbl Date: Sat, 28 Oct 2023 09:16:10 +0200 Subject: [PATCH 2/3] add thredds requirement --- .github/workflows/extra_requirements.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/extra_requirements.txt b/.github/workflows/extra_requirements.txt index f519e3c1ac..c6b6e7c4dc 100644 --- a/.github/workflows/extra_requirements.txt +++ b/.github/workflows/extra_requirements.txt @@ -4,3 +4,4 @@ scikit-learn GDAL==${GDAL_VERSION} pygbif cf_units +thredds_crawler From e6f0a719f677712b65adac7d2372dd403cdb8916 Mon Sep 17 00:00:00 2001 From: ninsbl Date: Wed, 19 Jun 2024 01:29:35 +0200 Subject: [PATCH 3/3] start converting to Python --- .../r.connectivity.network.py | 1021 ++++++++++------- 1 file changed, 623 insertions(+), 398 deletions(-) diff --git a/src/raster/r.connectivity/r.connectivity.network/r.connectivity.network.py b/src/raster/r.connectivity/r.connectivity.network/r.connectivity.network.py index a4823e9ccd..44e04cc92e 100755 --- a/src/raster/r.connectivity/r.connectivity.network/r.connectivity.network.py +++ b/src/raster/r.connectivity/r.connectivity.network/r.connectivity.network.py @@ -1,5 +1,4 @@ #!/usr/bin/env python -# -*- coding: UTF-8 -*- """ MODULE: r.connectivity.network @@ -73,8 +72,8 @@ - replace R with Python - add RGB columns instead of QML - Fix history assignment -- - grass.parse_command('v.support', map=edges, flags='g')[comments] -- - grass.parse_command('v.support', map=nodes, flags='g')[comments] +- - gscript.parse_command('v.support', map=edges, flags='g')[comments] +- - gscript.parse_command('v.support', map=nodes, flags='g')[comments] """ # %Module @@ -209,10 +208,11 @@ import warnings import numpy as np import matplotlib + # Required for Windows -matplotlib.use('wx') +matplotlib.use("wx") import matplotlib.pyplot as plt -import grass.script as grass +import grass.script as gscript import grass.script.task as task import grass.script.db as grass_db @@ -221,17 +221,18 @@ if "GISBASE" not in os.environ: sys.exit("You must be in GRASS GIS to run this program") + def cleanup(): - """tmp_maps = grass.read_command("g.list", + """tmp_maps = gscript.read_command("g.list", type=['vector', 'raster'], pattern='{}*'.format(TMP_prefix), separator=',') if tmp_maps: - grass.run_command("g.remove", type=['vector', 'raster'], + gscript.run_command("g.remove", type=['vector', 'raster'], pattern='{}*'.format(TMP_prefix), quiet=True, flags='f') - #grass.del_temp_region()""" + #gscript.del_temp_region()""" pass @@ -239,7 +240,7 @@ def is_number(string): """Check if a string can be converted to numeric""" try: - complex(string) # for int, long, float and complex + complex(string) # for int, long, float and complex except ValueError: return False @@ -253,119 +254,144 @@ def main(): try: import igraph except ImportError: - gscript.fatal("Could not import igraph library! To install it run:", - "pip install python-igraph") - #Input variables - network_map = options['input'] + gscript.fatal( + "Could not import igraph library! To install it run:", + "pip install python-igraph", + ) + # Input variables + network_map = options["input"] # network_mapset = network_map.split('@')[0] # network = network_map.split('@')[1] if len(network_map.split('@')) # > 1 else None - prefix = options['prefix'] - cores = options['cores'] - convergence_treshold = options['convergence_threshold'] - base = float(options['base']) - exponent = float(options['exponent']) - qml_style_dir = options['qml_style'] - if is_number(options['connectivity_cutoff']): - if float(options['connectivity_cutoff']) > 0: - connectivity_cutoff = float(options['connectivity_cutoff']) + prefix = options["prefix"] + cores = options["cores"] + convergence_treshold = options["convergence_threshold"] + base = float(options["base"]) + exponent = float(options["exponent"]) + qml_style_dir = options["qml_style"] + if is_number(options["connectivity_cutoff"]): + if float(options["connectivity_cutoff"]) > 0: + connectivity_cutoff = float(options["connectivity_cutoff"]) else: - grass.fatal('"connectivity_cutoff" has to be > 0.') + gscript.fatal('"connectivity_cutoff" has to be > 0.') else: - grass.fatal('Option "connectivity_cutoff" is not given as a number.') - lnbh_cutoff = options['lnbh_cutoff'] - cl_thresh = options['cl_thresh'] + gscript.fatal('Option "connectivity_cutoff" is not given as a number.') + lnbh_cutoff = options["lnbh_cutoff"] + cl_thresh = options["cl_thresh"] cd_cutoff = connectivity_cutoff - kernel_plot = options['kernel_plot'] - overview_plot = options['overview_plot'] + kernel_plot = options["kernel_plot"] + overview_plot = options["overview_plot"] - verbose = grass.verbosity() == 3 - overwrite = grass.overwrite() + verbose = gscript.verbosity() == 3 + overwrite = gscript.overwrite() - edge_output = '{}_edge_measures'.format(prefix) - edge_output_tmp = '{}_edge_measures_tmp'.format(prefix) - vertex_output = '{}_vertex_measures'.format(prefix) - vertex_output_tmp = '{}_vertex_measures_tmp'.format(prefix) - network_output = '{}_network_measures'.format(prefix) + edge_output = f"{prefix}_edge_measures" + edge_output_tmp = f"{prefix}_edge_measures_tmp" + vertex_output = f"{prefix}_vertex_measures" + vertex_output_tmp = f"{prefix}_vertex_measures_tmp" + network_output = f"{prefix}_network_measures" # Check if input parameters are correct for table in [vertex_output, edge_output, network_output]: - if grass.db.db_table_exist(table) and not overwrite: - grass.fatal('Table "{}" already exists. \ - Use --o flag to overwrite'.format(table)) + if gscript.db.db_table_exist(table) and not overwrite: + gscript.fatal( + _( + 'Table "{}" already exists. \ + Use --o flag to overwrite' + ).format(table) + ) if qml_style_dir: if not os.path.exists(qml_style_dir): - grass.fatal('QML output requested but directory "{}" \ - does not exists.'.format(qml_style_dir)) + gscript.fatal( + _( + 'QML output requested but directory "{}" \ + does not exists.' + ).format(qml_style_dir) + ) if not os.path.isdir(qml_style_dir): - grass.fatal('QML output requested but "{}" is not a \ - directory.'.format(qml_style_dir)) + gscript.fatal( + _( + 'QML output requested but "{}" is not a \ + directory.' + ).format(qml_style_dir) + ) if not os.access(qml_style_dir, os.R_OK): - grass.fatal('Output directory "{}" for QML files is not \ - writable.'.format(qml_style_dir)) + gscript.fatal( + _( + 'Output directory "{}" for QML files is not \ + writable.' + ).format(qml_style_dir) + ) for plot in [kernel_plot, overview_plot]: if plot: if not os.path.exists(os.path.dirname(plot)): - grass.fatal('Directory for output "{}" does not \ - exists.'.format(plot)) + gscript.fatal( + _( + 'Directory for output "{}" does not \ + exists.' + ).format(plot) + ) if not os.access(os.path.dirname(plot), os.R_OK): - grass.fatal('Output directory "{}" is not \ - writable.'.format(plot)) + gscript.fatal( + _( + 'Output directory "{}" is not \ + writable.' + ).format(plot) + ) # Visualise negative exponential decay kernel and exit - x_flag = flags['x'] + x_flag = flags["x"] # Calculate edge betweenness community (default is do not) # y_flag = flags['y'] # Install missing R packages - i_flag = flags['i'] + i_flag = flags["i"] # Remove indirect connections from network - r_flag = flags['r'] + r_flag = flags["r"] db_connection = grass_db.db_connection() # parse_command('db.connect', flags='g') - db_name = grass.read_command('db.databases').rstrip('\n') + db_name = gscript.read_command("db.databases").rstrip("\n") - net_hist_str = grass.parse_command('v.info', map=network_map, - flags='h', delimiter=':' - )['COMMAND'].split('\n')[0] - # grass.parse_command('v.info', map=network_map, flags='h' + net_hist_str = gscript.parse_command( + "v.info", map=network_map, flags="h", delimiter=":" + )["COMMAND"].split("\n")[0] + # gscript.parse_command('v.info', map=network_map, flags='h' # ).split('\n')[0].split(': ')[1] # Parsing goes wrong in some cases (extractig with -tp) as input dist_cmd_dict = task.cmdstring_to_tuple(net_hist_str) - command = os.environ['CMDLINE'] + command = os.environ["CMDLINE"] - dist_prefix = dist_cmd_dict[1]['prefix'] - pop_proxy = dist_cmd_dict[1]['pop_proxy'] + dist_prefix = dist_cmd_dict[1]["prefix"] + pop_proxy = dist_cmd_dict[1]["pop_proxy"] - in_vertices = '{}_vertices'.format(dist_prefix) + in_vertices = "{}_vertices".format(dist_prefix) # OS adjustment os_type = platform.system() - grass.verbose("prefix is {}".format(prefix)) - grass.verbose("cores is {}".format(cores)) - grass.verbose("convergence_treshold is {}".format(convergence_treshold)) - grass.verbose("base is {}".format(base)) - grass.verbose("exponent is {}".format(exponent)) - grass.verbose("cd_cutoff is {}".format(cd_cutoff)) - grass.verbose("connectivity_cutoff is {}".format(connectivity_cutoff)) - grass.verbose("lnbh_cutoff is {}".format(lnbh_cutoff)) - grass.verbose("kernel_plot is {}".format(kernel_plot)) - grass.verbose("overview_plot is {}".format(overview_plot)) - grass.verbose("dist_prefix is {}".format(dist_prefix)) - grass.verbose("EDGES is {}_edges".format(prefix)) - grass.verbose("VERTICES is {}_vertices".format(prefix)) - grass.verbose("cl_thresh is {}".format(cl_thresh)) - - - #Visualise the negative exponential decay kernel and exit (if requested) + gscript.verbose(_("prefix is {}").format(prefix)) + gscript.verbose(_("cores is {}").format(cores)) + gscript.verbose(_("convergence_treshold is {}").format(convergence_treshold)) + gscript.verbose(_("base is {}").format(base)) + gscript.verbose(_("exponent is {}").format(exponent)) + gscript.verbose(_("cd_cutoff is {}").format(cd_cutoff)) + gscript.verbose(_("connectivity_cutoff is {}").format(connectivity_cutoff)) + gscript.verbose(_("lnbh_cutoff is {}").format(lnbh_cutoff)) + gscript.verbose(_("kernel_plot is {}").format(kernel_plot)) + gscript.verbose(_("overview_plot is {}").format(overview_plot)) + gscript.verbose(_("dist_prefix is {}").format(dist_prefix)) + gscript.verbose(_("EDGES is {}_edges").format(prefix)) + gscript.verbose(_("VERTICES is {}_vertices").format(prefix)) + gscript.verbose(_("cl_thresh is {}").format(cl_thresh)) + + # Visualise the negative exponential decay kernel and exit (if requested) if x_flag or kernel_plot: if x_flag: plt.ion() @@ -374,8 +400,7 @@ def main(): # ax = plt.axes() x_values = np.linspace(0, cd_cutoff, 1000) - plt.plot(x_values, np.exp(base*(np.power(10.0, - exponent))*x_values)) + plt.plot(x_values, np.exp(base * (np.power(10.0, exponent)) * x_values)) # Simplified version could be # plt.plot(x, np.exp(exponent*x)) @@ -388,30 +413,63 @@ def main(): fig.savefig(kernel_plot) if cores > 1 and os_type == "Windows": - grass.warning('Parallel processing not yet supported on MS Windows. \ - Setting number of cores to 1.') + gscript.warning( + _( + "Parallel processing not yet supported on MS Windows. \ + Setting number of cores to 1." + ) + ) cores = 1 - - #Check igraph version + # Check igraph version if igraph.__version_info__ < (0, 6, 2): - grass.fatal('The required igraph version is 0.6.2 or \ + gscript.fatal( + _( + "The required igraph version is 0.6.2 or \ later, the installed version is {}. \ Please download and install at least igraph \ - version 0.6.2'.format(".".join(igraph.__version_info__))) + version 0.6.2" + ).format(".".join(igraph.__version_info__)) + ) - driver = db_connection['driver'] + driver = db_connection["driver"] database = db_name - schema = db_connection['schema'] + schema = db_connection["schema"] command = command - gscript.verbose("Reading and preparing input data...\n") - vertices = np.gefromtext(BytesIO(gscript.read_command('v.db.select', flags="c", map=in_vertices, columns="cat,{}".format(pop_proxy), separator=',').strip("\n")), dtype=None, names=['patch_id', 'pop_proxy'], delimiter=",") + gscript.verbose(_("Reading and preparing input data...\n")) + vertices = np.gefromtext( + BytesIO( + gscript.read_command( + "v.db.select", + flags="c", + map=in_vertices, + columns="cat,{}".format(pop_proxy), + separator=",", + ).strip("\n") + ), + dtype=None, + names=["patch_id", "pop_proxy"], + delimiter=",", + ) - #Read edge data - edges = np.gefromtext(BytesIO(gscript.read_command('v.db.select', flags="c", map=in_edges, columns='cat,from_p,to_p,dist', separator=',').strip("\n")), dtype=None, names=['cat', 'from_patch', 'to_patch', 'distance'], delimiter=",") + # Read edge data + edges = np.gefromtext( + BytesIO( + gscript.read_command( + "v.db.select", + flags="c", + map=in_edges, + columns="cat,from_p,to_p,dist", + separator=",", + ).strip("\n") + ), + dtype=None, + names=["cat", "from_patch", "to_patch", "distance"], + delimiter=",", + ) - # THe following can be done in igraph graph.simplify(combine_edges=dict(weight="mean")) + # The following can be done in igraph graph.simplify(combine_edges=dict(weight="mean")) # https://igraph.org/python/doc/igraph.GraphBase-class.html#simplify """con_id_u_pre = np.unique(np.sort(e_df[:,[1,2]]), axis=0) con_id_u_pre_df = data.frame(con_id_u_pre=unique(con_id_u_pre), con_id_u=1:length(unique(con_id_u_pre))) @@ -423,48 +481,53 @@ def main(): # Merge vertices and edges to graph-object # Build directed graph - gscript.verbose("Building the graph...") + gscript.verbose(_("Building the graph...")) g = Graph().as_directed() - g.add_vertices(vertices[:,0].astype("i").astype("S")) - g.vs["patch_id"] = vertices[:,0].astype("i") - g.vs["pop_proxy"] = vertices[:,1] - - g.add_edges(edges[:,[1,2]].astype("i").astype("S")) - g.es["cat"] = edges[:,0].astype("i") - g.es["from_patch"] = edges[:,1].astype("i").astype("S") - g.es["to_patch"] = edges[:,2].astype("i").astype("S") - g.es["cost_distance"] = edges[:,3] - - #g = add.edges(g, t(edges), con_id=e$con_id, con_id_u=e$con_id_u, from_p=e$from_p, from_pop=e$from_pop, to_p=e$to_p, to_pop=e$to_pop, cost_distance=e$cost_distance, cd_u=e$cd_u, distance_weight_e=e$distance_weight_e, distance_weight_e_ud=e$distance_weight_e_ud, mf_o=e$mf_o, mf_o_inv=e$mf_o_inv, mf_i=e$mf_i, mf_i_inv=e$mf_i_inv, mf_u=e$mf_u, mf_inv_u=e$mf_inv_u, cf=e$cf, cf_inv=e$cf_inv, cf_u=e$cf_u, cf_inv_u=e$cf_inv_u) + g.add_vertices(vertices[:, 0].astype("i").astype("S")) + g.vs["patch_id"] = vertices[:, 0].astype("i") + g.vs["pop_proxy"] = vertices[:, 1] + + g.add_edges(edges[:, [1, 2]].astype("i").astype("S")) + g.es["cat"] = edges[:, 0].astype("i") + g.es["from_patch"] = edges[:, 1].astype("i").astype("S") + g.es["to_patch"] = edges[:, 2].astype("i").astype("S") + g.es["cost_distance"] = edges[:, 3] + + # g = add.edges(g, t(edges), con_id=e$con_id, con_id_u=e$con_id_u, from_p=e$from_p, from_pop=e$from_pop, to_p=e$to_p, to_pop=e$to_pop, cost_distance=e$cost_distance, cd_u=e$cd_u, distance_weight_e=e$distance_weight_e, distance_weight_e_ud=e$distance_weight_e_ud, mf_o=e$mf_o, mf_o_inv=e$mf_o_inv, mf_i=e$mf_i, mf_i_inv=e$mf_i_inv, mf_u=e$mf_u, mf_inv_u=e$mf_inv_u, cf=e$cf, cf_inv=e$cf_inv, cf_u=e$cf_u, cf_inv_u=e$cf_inv_u) ##Remove connections longer than cost distance threshold if requested - #if(remove_longer_cutoff == 1) {{ - #con_id_True = e_grouped_df$con_id[grep(True, e_grouped_df$dist_ud10000000000000000,10000000000000000,max(g.es["maximum_flow_sum"]))-ifelse(min(g.es["maximum_flow_sum"])<0.000000000001,0.000000000001,min(g.es["maximum_flow_sum"])))/(max(g.es["maximum_flow_sum"]*-1)-min(g.es["maximum_flow_sum"]*-1))+ifelse(min(g.es["maximum_flow_sum"])<0.000000000001,0.000000000001,min(g.es["maximum_flow_sum"])) + g.es["maximum_flow_out_inverse"] = g.es["mf_o"] ** -1 + g.es["maximum_flow_in_inverse"] = g.es["mf_i"] ** -1 + # g.es["maximum_flow_sum_inverse"] = ((g.es["maximum_flow_sum"]*-1)-min(g.es["maximum_flow_sum"]*-1))*(ifelse(max(g.es["maximum_flow_sum"])>10000000000000000,10000000000000000,max(g.es["maximum_flow_sum"]))-ifelse(min(g.es["maximum_flow_sum"])<0.000000000001,0.000000000001,min(g.es["maximum_flow_sum"])))/(max(g.es["maximum_flow_sum"]*-1)-min(g.es["maximum_flow_sum"]*-1))+ifelse(min(g.es["maximum_flow_sum"])<0.000000000001,0.000000000001,min(g.es["maximum_flow_sum"])) + """ flow_df = data.frame(from_p=e$from_p, mf_i=mf_i) sum_mf_i_groups = groupedData(mf_i ~ from_p | from_p, data=as.data.frame(flow_df), FUN=sum) sum_mf_i_pre = gsummary(sum_mf_i_groups, sum) @@ -493,13 +556,21 @@ def main(): e_ud_pre = merge(e, ud_grouped_df, all.x=True) first_con_id = data.frame(con_id=e_ud_pre$con_id[e_ud_pre$con_id_avg>e_ud_pre$con_id]) e_ud = merge(first_con_id, e) + """ # Build undirected graph # graph.simplify(combine_edges=dict(weight="mean")) # https://igraph.org/python/doc/igraph.GraphBase-class.html#simplify g_ud = deepcopy(g) - g_ud = g_ud.simplify(combine_edges=dict(mf_o="mean", mf_i="mean", )) - g_ud.es["exponential_distance_weight_undirected"] = np.exp((basis*(10^exponent))*np.array(g_ud.es["cost_distance_undirected"])) + g_ud = g_ud.simplify( + combine_edges=dict( + mf_o="mean", + mf_i="mean", + ) + ) + g_ud.es["exponential_distance_weight_undirected"] = np.exp( + (basis * (10 ^ exponent)) * np.array(g_ud.es["cost_distance_undirected"]) + ) """ g_ud = add.edges(g_ud, t(edges), con_id=e_ud$con_id, @@ -534,6 +605,7 @@ def main(): # g.es["cd_u = unlist(mclapply(mc.cores={0}, 0:(length(E(g))-1), function(x) (E(g)[x]$cost_distance+E(g, path=c(ends(g, x)[2], ends(g, x)[1]))$cost_distance)/2)) + """ #################################### ### Can the following be vectorised???? @@ -541,7 +613,10 @@ def main(): g_ud.es["is_short_mf"] = as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud))), function(x) length(get.shortest.paths(g_ud, from=ends(g_ud, x)[1], to=ends(g_ud, x)[2], weights=g_ud.es["mf_inv_u)$vpath[[1]])))==2) g_ud.es["is_short_cf"] = as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud))), function(x) length(get.shortest.paths(g_ud, from=ends(g_ud, x)[1], to=ends(g_ud, x)[2], weights=g_ud.es["cf_inv_u)$vpath[[1]])))==2) # g_ud.es["is_short"] = ifelse((g_ud.es["isshort_cd==0 & g_ud.es["isshort_mf==0 & g_ud.es["isshort_cf==0), 0, 1) - g_ud.es["is_short"] = np.min(g_ud.es["is_short_cd"], g_ud.es["is_short_mf"], g_ud.es["is_short_cf"]) + """ + g_ud.es["is_short"] = np.min( + g_ud.es["is_short_cd"], g_ud.es["is_short_mf"], g_ud.es["is_short_cf"] + ) #################################### @@ -552,20 +627,21 @@ def main(): # Remove indirect connections if requested if remove_indirect == 1: g_ud_d = deepcopy(g_ud) - g_ud_d.delete.edges(g_ud.es.select[is_short_eq=0]) + g_ud_d.delete.edges(g_ud.es.select(is_short_eq=0)) # Remove connections above connectivity threshold if requested + """ if connectivity_cutoff >= 0.0: g_ud_cd = deepcopy(g_ud) g_ud_cd = g_ud_cd.delete_edges(E(g_ud)[grep(True, g_ud.es["cd_u>=connectivity_cutoff)]) g_ud_d_cd = deepcopy(g_ud_d) g_ud_d_cd = g_ud_d_cd.delete_edges(E(g_ud_d)[grep(True, g_ud_d.es["cd_u>=connectivity_cutoff)]) - + """ ######################################################################## ### Analysis on graph level ######################################################################## - gscript.verbose("Starting analysis on graph level...") + gscript.verbose(_("Starting analysis on graph level...")) ### graph measures vertices_n = len(g_ud.vs) @@ -574,15 +650,16 @@ def main(): edges_n_cd = len(g_ud_cd.es) edges_n_d_cd = len(g_ud_d_cd.es) - ######################################################################## ###### Calculate number of clusters # On the undirected graph with only direct edges and on the undirected graph with only direct edges shorter cost distance threshold cl_no_d = g_ud_d.clusters()["no"] cl_no_d_cd = g_ud_d_cd.clusters()["no"] + """ cls_size_d = as.vector(gsummary(groupedData(pop_size ~ 1 | cl, data.frame(cl=clusters(g_ud_d)$membership, pop_size=g_ud_d.vs["pop_proxy), order.groups=True, FUN=sum), sum)$pop_size) cls_size_d_cd = as.vector(gsummary(groupedData(pop_size ~ 1 | cl, data.frame(cl=clusters(g_ud_d_cd)$membership, pop_size=g_ud_d_cd.vs["pop_proxy), order.groups=True, FUN=sum), sum)$pop_size) + """ cl_max_size_d = max(cls_size_d) cl_max_size_d_cd = max(cls_size_d_cd) @@ -590,15 +667,20 @@ def main(): cl_mean_size_d = mean(cls_size_d) cl_mean_size_d_cd = mean(cls_size_d_cd) - diam = g_ud.diameter(directed=False, unconnected=True, weights=g_ud.es["cd_u) - diam_d = g_ud_d.diameter(directed=False, unconnected=True, weights=g_ud_d.es["cd_u) - diam_d_cd = g_ud_d_cd.diameter(directed=False, unconnected=True, weights=E(g_ud_d_cd)$cd_u) + diam = g_ud.diameter(directed=False, unconnected=True, weights=g_ud.es["cd_u"]) + diam_d = g_ud_d.diameter( + directed=False, unconnected=True, weights=g_ud_d.es["cd_u"] + ) + diam_d_cd = g_ud_d_cd.diameter( + directed=False, unconnected=True, weights=g_ud_d_cd.es["cd_u"] + ) density = g.density(g, loops=False) density_u = g_ud.density(g_ud, loops=False) density_ud = g_ud_d.density(g_ud_d, loops=False) density_udc = g_ud_d_cd.density(g_ud_d_cd, loops=False) + """ if network_overview_ps: ###### Edge removal operations for edge_length in np.sort(g_ud.es["cd_u"], decreasing=True, na.last=NA): @@ -625,7 +707,7 @@ def main(): lines(df_edgeremoval$distance/(10^(nchar(as.integer(max(df_edgeremoval$distance)))-2)), df_edgeremoval$cl_del_diam*100/diam_d, type="l", lty=4) # Lable axis 4!!! # axis(4, seq.int(0, 100, 25), yaxs="i", labels=seq.int(0, ceiling((as.integer(edges_n)/(10^(nchar(as.integer(edges_n))-2))))*(10^(nchar(as.integer(edges_n))-2)), ceiling((as.integer(edges_n)/(10^(nchar(as.integer(edges_n))-2))))*(10^(nchar(as.integer(edges_n))-2))/4)) - + """ ######################################################################## ###Analysis on edge level ######################################################################## @@ -637,74 +719,77 @@ def main(): ### Minimum spanning tree weighted by maximum potential flow g_ud_d.es["mf_mst_ud"] = 0 - E(g_ud_d.es["mf_mst_ud"] = [grep(True, g_ud_d.es["con_id %in% E(minimum.spanning.tree(g_ud_d, weights=g_ud_d.es["mf_inv_u))$con_id)]$mf_mst_ud = 1 + # E(g_ud_d.es["mf_mst_ud"] = [grep(True, g_ud_d.es["con_id %in% E(minimum.spanning.tree(g_ud_d, weights=g_ud_d.es["mf_inv_u))$con_id)]$mf_mst_ud = 1 ### Minimum spanning tree weighted by cost distance g_ud_d.es["cd_mst_ud"] = 0 - g_ud_d.es["cd_mst_ud"] = [grep(True, g_ud_d.es["con_id %in% E(minimum.spanning.tree(g_ud_d, weights=g_ud_d.es["cd_u))$con_id)]$cd_mst_ud = 1 + # g_ud_d.es["cd_mst_ud"] = [grep(True, g_ud_d.es["con_id %in% E(minimum.spanning.tree(g_ud_d, weights=g_ud_d.es["cd_u))$con_id)]$cd_mst_ud = 1 g_ud_d.es["cf_mst_ud"] = 0 - g_ud_d.es["cf_mst_ud"] = [grep(True, g_ud_d.es["con_id %in% E(minimum.spanning.tree(g_ud_d, weights=g_ud_d.es["cf_inv_u))$con_id)]$cf_mst_ud = 1 + # g_ud_d.es["cf_mst_ud"] = [grep(True, g_ud_d.es["con_id %in% E(minimum.spanning.tree(g_ud_d, weights=g_ud_d.es["cf_inv_u))$con_id)]$cf_mst_ud = 1 ###### Calculate minimum spanning trees (MST) on the undirected graph with only direct edges shorter cost distance threshold ### Minimum spanning tree weighted by maximum potential flow g_ud_d_cd.es["mf_mst_udc"] = 0 - g_ud_d_cd.es["mf_mst_udc"] = [grep(True, E(g_ud_d_cd)$con_id %in% E(minimum.spanning.tree(g_ud_d_cd, weights=E(g_ud_d_cd)$mf_inv_u))$con_id)]$mf_mst_udc = 1 - + # g_ud_d_cd.es["mf_mst_udc"] = [grep(True, E(g_ud_d_cd)$con_id %in% E(minimum.spanning.tree(g_ud_d_cd, weights=E(g_ud_d_cd)$mf_inv_u))$con_id)]$mf_mst_udc = 1 ### Minimum spanning tree weighted by cost distance g_ud_d_cd.es["cd_mst_udc"] = 0 - g_ud_d_cd.es["cd_mst_udc"] = [grep(True, E(g_ud_d_cd)$con_id %in% E(minimum.spanning.tree(g_ud_d_cd, weights=E(g_ud_d_cd)$cd_u))$con_id)]$cd_mst_udc = 1 + # g_ud_d_cd.es["cd_mst_udc"] = [grep(True, E(g_ud_d_cd)$con_id %in% E(minimum.spanning.tree(g_ud_d_cd, weights=E(g_ud_d_cd)$cd_u))$con_id)]$cd_mst_udc = 1 ### Minimum spanning tree weighted by competing potential flow g_ud_d_cd.es["cf_mst_udc"] = 0 - g_ud_d_cd.es["cf_mst_udc"] = [grep(True, E(g_ud_d_cd)$con_id %in% E(minimum.spanning.tree(g_ud_d_cd, weights=E(g_ud_d_cd)$cf_inv_u))$con_id)]$cf_mst_udc = 1 + # g_ud_d_cd.es["cf_mst_udc"] = [grep(True, E(g_ud_d_cd)$con_id %in% E(minimum.spanning.tree(g_ud_d_cd, weights=E(g_ud_d_cd)$cf_inv_u))$con_id)]$cf_mst_udc = 1 ######################################################################## ### Identify bridges for the undirected graph - E(g_ud.es["is_br_u"] = as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud))), function(x) clusters(delete.edges(g_ud, E(g_ud)[as.integer(x)]))$no-cl_no_d))) + # E(g_ud.es["is_br_u"] = as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud))), function(x) clusters(delete.edges(g_ud, E(g_ud)[as.integer(x)]))$no-cl_no_d))) # Identify bridges for the undirected graph with only direct edges - E(g_ud_d.es["is_br_ud"] = as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud_d))), function(x) clusters(delete.edges(g_ud_d, E(g_ud_d)[as.integer(x)]))$no-cl_no_d))) + # E(g_ud_d.es["is_br_ud"] = as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud_d))), function(x) clusters(delete.edges(g_ud_d, E(g_ud_d)[as.integer(x)]))$no-cl_no_d))) ### Identify bridges for the undirected graph with only direct edges shorter cost distance threshold - E(g_ud_d_cd.es["is_br_udc"] = as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud_d_cd))), function(x) clusters(delete.edges(g_ud_d_cd, E(g_ud_d_cd)[as.integer(x)]))$no-cl_no_d_cd))) + # E(g_ud_d_cd.es["is_br_udc"] = as.integer(unlist(mclapply(mc.cores={0}, 1:(length(E(g_ud_d_cd))), function(x) clusters(delete.edges(g_ud_d_cd, E(g_ud_d_cd)[as.integer(x)]))$no-cl_no_d_cd))) ######################################################################## biconnected_components_d = biconnected.components(g_ud) # Identify component edges (biconnected components) for the undirected graph with only direct edges - g_ud.es["bc_e_u = 0 - for c in 1:biconnected_components_d$no: - g_ud.es["bc_e_u[unlist(biconnected_components_d$component_edges[c])] = c + g_ud.es["bc_e_u"] = 0 + for c in range(biconnected_components_d["no"]): + g_ud.es["bc_e_u"][unlist(biconnected_components_d["component_edges"][c])] = c # Identify tree edges (biconnected components) for the undirected graph with only direct edges - g_ud.es["bc_te_u = 0 - for c in 1:biconnected_components_d$no: - g_ud.es["bc_te_u[unlist(biconnected_components_d$tree_edges[c])] = c + g_ud.es["bc_te_u"] = 0 + for c in range(biconnected_components_d["no"]): + g_ud.es["bc_te_u"][unlist(biconnected_components_d["tree_edges"][c])] = c biconnected_components_d = biconnected.components(g_ud_d) # Identify component edges (biconnected components) for the undirected graph with only direct edges - g_ud_d.es["bc_e_ud = 0 - for c in 1:biconnected_components_d$no: - g_ud_d.es["bc_e_ud[unlist(biconnected_components_d$component_edges[c])] = c + g_ud_d.es["bc_e_ud"] = 0 + for c in range(biconnected_components_d["no"]): + g_ud_d.es["bc_e_ud"][unlist(biconnected_components_d["component_edges"][c])] = c # Identify tree edges (biconnected components) for the undirected graph with only direct edges - g_ud_d.es["bc_te_ud = 0 - for c in 1:biconnected_components_d$no: - g_ud_d.es["bc_te_ud[unlist(biconnected_components_d$tree_edges[c])] = c + g_ud_d.es["bc_te_ud"] = 0 + for c in range(biconnected_components_d["no"]): + g_ud_d.es["bc_te_ud"][unlist(biconnected_components_d["tree_edges"][c])] = c biconnected_components_d_cd = biconnected.components(g_ud_d_cd) # Identify component edges (biconnected components) for the undirected graph with only direct edges shorter cost distance threshold - E(g_ud_d_cd)$bc_e_udc = 0 - for c in 1:biconnected_components_d_cd$no: - E(g_ud_d_cd)$bc_e_udc[unlist(biconnected_components_d_cd$component_edges[c])] = c + g_ud_d_cd.es["bc_e_udc"] = 0 + for c in range(biconnected_components_d_cd["no"]): + g_ud_d_cd.es["bc_e_udc"][ + unlist(biconnected_components_d_cd["component_edges"][c]) + ] = c # Identify tree edges (biconnected components) for the undirected graph with only direct edges shorter cost distance threshold - E(g_ud_d_cd)$bc_te_udc = 0 - for c in 1:biconnected_components_d_cd$no: - E(g_ud_d_cd)$bc_te_udc[unlist(biconnected_components_d_cd$tree_edges[c])] = c + g_ud_d_cd.es["bc_te_udc"] = 0 + for c in range(biconnected_components_d_cd["no"]): + g_ud_d_cd.es["bc_te_udc"][ + unlist(biconnected_components_d_cd["tree_edges"][c]) + ] = c ######################################################################## ### Calculate edge betweenness @@ -715,54 +800,79 @@ def main(): ### Calculate edge betweenness on the entire undirected graph with only direct edges - path_lengths = shortest.paths(g_ud_d, v=V(g_ud_d), to=V(g_ud_d), weights=g_ud_d.es["cd_u) + path_lengths = shortest.paths( + g_ud_d, v=g_ud_d.vs, to=g_ud_d.vs, weights=g_ud_d.es["cd_u"] + ) # weighted by cost distance - E(g_ud_d.es.["cd_eb_ud"] = edge.betweenness(g_ud_d, e=E(g_ud_d), directed=False, weights=g_ud_d.es["cd_u) - E(g_ud_d.es.["cd_leb_ud"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d)), function(x) get.shortest.paths(g_ud_d, x, V(g_ud_d)[grep(True, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=g_ud_d.es["cd_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d))), plot=False)$counts/2 + g_ud_d.es["cd_eb_ud"] = edge.betweenness( + g_ud_d, e=g_ud_d.es, directed=False, weights=g_ud_d.es["cd_u"] + ) + # g_ud_d.es.["cd_leb_ud"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d)), function(x) get.shortest.paths(g_ud_d, x, V(g_ud_d)[grep(True, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=g_ud_d.es["cd_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d))), plot=False)$counts/2 # g_ud_d.es["cd_leb_ud = edge.betweenness.estimate(g_ud_d, e=E(g_ud_d), directed=False, lnbh_cutoff*connectivity_cutoff, weights=g_ud_d.es["cd_u) # weighted by maximum potential flow - E(g_ud_d.es.["mf_eb_ud"] = edge.betweenness(g_ud_d, e=E(g_ud_d), directed=False, weights=g_ud_d.es["mf_inv_u) - E(g_ud_d.es.["mf_leb_ud"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d)), function(x) get.shortest.paths(g_ud_d, x, V(g_ud_d)[grep(True, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=g_ud_d.es["mf_inv_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d))), plot=False)$counts/2 + g_ud_d.es["mf_eb_ud"] = edge.betweenness( + g_ud_d, e=g_ud_d.es, directed=False, weights=g_ud_d.es["mf_inv_u"] + ) + # g_ud_d.es.["mf_leb_ud"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d)), function(x) get.shortest.paths(g_ud_d, x, V(g_ud_d)[grep(True, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=g_ud_d.es["mf_inv_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d))), plot=False)$counts/2 # weighted by competing potential flow - E(g_ud_d.es.["cf_eb_ud"] = edge.betweenness(g_ud_d, e=E(g_ud_d), directed=False, weights=g_ud_d.es["cf_inv_u) - E(g_ud_d.es.["cf_leb_ud"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d)), function(x) get.shortest.paths(g_ud_d, x, V(g_ud_d)[grep(True, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=g_ud_d.es["cf_inv_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d))), plot=False)$counts/2 + g_ud_d.es["cf_eb_ud"] = edge.betweenness( + g_ud_d, e=g_ud_d.es, directed=False, weights=g_ud_d.es["cf_inv_u"] + ) + # g_ud_d.es.["cf_leb_ud"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d)), function(x) get.shortest.paths(g_ud_d, x, V(g_ud_d)[grep(True, path_lengths[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=g_ud_d.es["cf_inv_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d))), plot=False)$counts/2 ######################################################################## ### Calculate edge betweenness on the entire undirected graph with only direct edges shorter cost distance threshold - path_lengths_cd = shortest.paths(g_ud_d_cd, v=V(g_ud_d_cd), to=V(g_ud_d_cd), weights=E(g_ud_d_cd)$cd_u) + path_lengths_cd = shortest.paths( + g_ud_d_cd, v=g_ud_d_cd.vs, to=g_ud_d_cd.vs, weights=g_ud_d_cd.es["cd_u"] + ) # weighted by cost distance - g_ud_d_cd.es["cd_eb_udc"] = edge.betweenness(g_ud_d_cd, e=E(g_ud_d_cd), directed=False, weights=E(g_ud_d_cd)$cd_u) - g_ud_d_cd.es["cd_leb_udc"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d_cd)), function(x) get.shortest.paths(g_ud_d_cd, x, V(g_ud_d_cd)[grep(True, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)$cd_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d_cd))), plot=False)$counts/2 - + g_ud_d_cd.es["cd_eb_udc"] = edge.betweenness( + g_ud_d_cd, e=g_ud_d_cd.es, directed=False, weights=g_ud_d_cd.es["cd_u"] + ) + # g_ud_d_cd.es["cd_leb_udc"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d_cd)), function(x) get.shortest.paths(g_ud_d_cd, x, V(g_ud_d_cd)[grep(True, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)$cd_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d_cd))), plot=False)$counts/2 # weighted by maximum potential flow - E(g_ud_d_cd.es["mf_eb_udc"] = edge.betweenness(g_ud_d_cd, e=E(g_ud_d_cd), directed=False, weights=E(g_ud_d_cd)$mf_inv_u) - E(g_ud_d_cd.es["mf_leb_udc"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d_cd)), function(x) get.shortest.paths(g_ud_d_cd, x, V(g_ud_d_cd)[grep(True, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)$mf_inv_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d_cd))), plot=False)$counts/2 + g_ud_d_cd.es["mf_eb_udc"] = edge.betweenness( + g_ud_d_cd, e=g_ud_d_cd.es, directed=False, weights=g_ud_d_cd.es["mf_inv_u"] + ) + # g_ud_d_cd.es["mf_leb_udc"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d_cd)), function(x) get.shortest.paths(g_ud_d_cd, x, V(g_ud_d_cd)[grep(True, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)$mf_inv_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d_cd))), plot=False)$counts/2 # weighted by competing potential flow - E(g_ud_d_cd.es["cf_eb_udc"] = edge.betweenness(g_ud_d_cd, e=E(g_ud_d_cd), directed=False, weights=E(g_ud_d_cd)$cf_inv_u) - E(g_ud_d_cd.es["cf_leb_udc"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d_cd)), function(x) get.shortest.paths(g_ud_d_cd, x, V(g_ud_d_cd)[grep(True, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)$cf_inv_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d_cd))), plot=False)$counts/2 + g_ud_d_cd.es["cf_eb_udc"] = edge.betweenness( + g_ud_d_cd, e=g_ud_d_cd.es, directed=False, weights=g_ud_d_cd.es["cf_inv_u"] + ) + # g_ud_d_cd.es["cf_leb_udc"] = hist(unlist(mclapply(mc.cores={0}, 1:length(V(g_ud_d_cd)), function(x) get.shortest.paths(g_ud_d_cd, x, V(g_ud_d_cd)[grep(True, path_lengths_cd[x,]<(lnbh_cutoff*connectivity_cutoff))], weights=E(g_ud_d_cd)$cf_inv_u, output=c("epath"))$epath)), breaks=0:(length(E(g_ud_d_cd))), plot=False)$counts/2 if cl_thresh > 0: ######################################################################### ####### Calculate edge betweenness community on the undirected graph ## weighted by competing potential flow - ebc = edge.betweenness.community(g_ud, weights=g_ud.es["cf_inv_u, directed=False, edge.betweenness=True, merges=True, bridges=True, modularity=True, membership=True) - - E(g_ud.es["cf_iebc_v"] = ebc$edge.betweenness + ebc = edge.betweenness.community( + g_ud, + weights=g_ud.es["cf_inv_u"], + directed=False, + edge_betweenness=True, + merges=True, + bridges=True, + modularity=True, + membership=True, + ) + + g_ud.es["cf_iebc_v"] = ebc["edge.betweenness"] # cf_iebc_m = ebc$merges - E(g_ud.es["cf_iebc_r"] = ebc$removed.edges - E(g_ud.es["cf_iebc_b"] = unlist(mclapply(mc.cores={0}, 1:length(E(g_ud)), function(x) ifelse(x %in% ebc$bridges, 1, 0))) - E(g_ud.es["cf_iebc_c"] = crossing(ebc, g_ud) + g_ud.es["cf_iebc_r"] = ebc["removed.edges"] + # g_ud.es["cf_iebc_b"] = unlist(mclapply(mc.cores={0}, 1:length(E(g_ud)), function(x) ifelse(x %in% ebc$bridges, 1, 0))) + g_ud.es["cf_iebc_c"] = crossing(ebc, g_ud) # g.vs["cf_iebc_mo = ebc$modularity - V(g.vs["cf_iebc_me"] = ebc$membership + g.vs["cf_iebc_me"] = ebc.es["membership"] + """ com_struct = as.character(cutat(ebc, cl_no_d)) for i in (cl_no_d+1):(cl_no_d+cl_thresh): com_struct = paste(com_struct, cutat(ebc, i), sep=';') @@ -774,6 +884,7 @@ def main(): com_no_u = length(ebc) com_sizes_u = sizes(ebc) com_sizes_u_names = paste("Size of comumity ", names(sizes(ebc)), " (at maximum modularity score (from iebc))", sep="") + """ # The following workaround was written when igraph did not support weights when calculating edge.betweenness.community # The results of this procedure differ significantly from the (new) edge.betweenness.community in igraph 0.6-2 with weights (above) @@ -784,55 +895,59 @@ def main(): ebc_rank = 0 ebc_value = 0 ebc_clust = 0 - g_ud.vs["cf_ebc_cs = as.character(clusters(g_ud)$membership) + g_ud.vs["cf_ebc_cs"] = str(clusters(g_ud)["membership"]) del_eb_g = g_ud for edge in del_eb_g.es: # Recalculate edge betweenness for the remaining edges - E(del_eb_g)$cf_eb_ud = edge.betweenness(del_eb_g, e=E(del_eb_g), directed=False, weights=E(del_eb_g)$cf_inv_u) + del_eb_g.es["cf_eb_ud"] = edge.betweenness( + del_eb_g, e=del_eb_g.es, directed=False, weights=del_eb_g.es["cf_inv_u"] + ) # Identify edge with highest edge betweenness value - idx = sort.int(E(del_eb_g)$cf_eb_ud, decreasing=True, na.last=NA, index.return=True)$ix[1] + # idx = sort.int(del_eb_g.es["cf_eb_ud"], decreasing=True, na.last=NA, index.return=True).es["ix"][1] # Save edge betweenness value for this edge - ebc_con_id_u[n] = E(del_eb_g)$con_id_u[idx] + ebc_con_id_u[n] = del_eb_g.es["con_id_u"][idx] ebc_rank[n] = n - ebc_clust[n] = clusters(del_eb_g)$no - if((ebc_clust[n]-cl_no_d)<(cl_no_d+cl_thresh)) {{if(n>1) {{if(ebc_clust[n]>ebc_clust[n-1]) g_ud.vs["cf_ebc_cs = paste(g_ud.vs["cf_ebc_cs, clusters(del_eb_g)$membership, sep=";")}}}} - ebc_value[n] = E(del_eb_g)$cf_eb_ud[idx] + ebc_clust[n] = clusters(del_eb_g)["no"] + # if((ebc_clust[n]-cl_no_d)<(cl_no_d+cl_thresh)) {{if(n>1) {{if(ebc_clust[n]>ebc_clust[n-1]) g_ud.vs["cf_ebc_cs = paste(g_ud.vs["cf_ebc_cs, clusters(del_eb_g)$membership, sep=";")}}}} + ebc_value[n] = del_eb_g.es["cf_eb_ud"][idx] # Delete edge with highest edge betweenness value - del_eb_g = delete.edges(del_eb_g, E(del_eb_g)[idx]) - n = n+1 + del_eb_g = delete.edges(del_eb_g, del_eb_g.es[idx]) + n += 1 - g_ud.vs["cf_ebc_cl = unlist(mclapply(mc.cores={0}, 1:length(V(g_ud)), function(x) strsplit(g_ud.vs["cf_ebc_cs, ";")[[x]][(cl_no_d+cl_thresh)])) + # g_ud.vs["cf_ebc_cl"] = unlist(mclapply(mc.cores={0}, 1:length(V(g_ud)), function(x) strsplit(g_ud.vs["cf_ebc_cs, ";")[[x]][(cl_no_d+cl_thresh)])) # ebc_df = merge(data.frame(con_id=g.es["con_id, con_id_u=g.es["con_id_u), data.frame(con_id_u=ebc_con_id_u, ebc_rank=ebc_rank, ebc_clust=ebc_clust, ebc_value=ebc_value), all.x=True) - idx = sort.int(ebc_con_id_u, index.return=True)$ix - g_ud.es["cf_ebc_v = ebc_value[idx] - g_ud.es["cf_ebc_r = ebc_rank[idx] - g_ud.es["cf_ebc_c = ebc_clust[idx] + # idx = sort.int(ebc_con_id_u, index.return=True)["ix"] + g_ud.es["cf_ebc_v"] = ebc_value[idx] + g_ud.es["cf_ebc_r"] = ebc_rank[idx] + g_ud.es["cf_ebc_c"] = ebc_clust[idx] - g_ud.es["cf_ebc_vi = strftime(Sys.time()+g_ud.es["cf_ebc_r) + g_ud.es["cf_ebc_vi"] = strftime(Sys.time() + g_ud.es["cf_ebc_r"]) ######################################################################### ### Calculate community connectors (based on community structure from ebc) - com_membership = data.frame(patch_id=g_ud_d_cd.vs["patch_id, com=g_ud.vs["cf_ebc_cl) + """ + com_membership = data.frame(patch_id=g_ud_d_cd.vs["patch_id"], com=g_ud.vs["cf_ebc_cl"]) com_pc_pre = merge(merge(data.frame(id=1:length(E(g)), from_p=g.es["from_p, to_p=g.es["to_p), com_membership, by.x="from_p", by.y="patch_id"), com_membership, by.x="to_p", by.y="patch_id")[order(merge(merge(data.frame(id=1:length(E(g)), from_p=g.es["from_p, to_p=g.es["to_p), com_membership, by.x="from_p", by.y="patch_id"), com_membership, by.x="to_p", by.y="patch_id")$id) ,] - g.es["cf_ebc_cc = ifelse(com_pc_pre$com.x==com_pc_pre$com.y,0,1) - + g.es["cf_ebc_cc"] = ifelse(com_pc_pre$com.x==com_pc_pre$com.y,0,1) + """ else: - gscript.verbose("Skipping comunity algorithms...") + gscript.verbose(_("Skipping comunity algorithms...")) ######################################################################## ###### Calculate cluster membership - g_ud_d.vs["cl_ud = clusters(g_ud_d)$membership - g_ud_d_cd.vs["cl_udc = clusters(g_ud_d_cd)$membership + g_ud_d.vs["cl_ud"] = clusters(g_ud_d)["membership"] + g_ud_d_cd.vs["cl_udc"] = clusters(g_ud_d_cd)["membership"] ######################################################################### ### Calculate potential cluster connectors (based on cost distance threshold) + """ cl_membership = data.frame(patch_id=g_ud_d_cd.vs["patch_id, cl=g_ud_d_cd.vs["cl_udc) cl_pc_pre = merge(merge(data.frame(id=1:length(E(g)), from_p=g.es["from_p, to_p=g.es["to_p), cl_membership, by.x="from_p", by.y="patch_id"), cl_membership, by.x="to_p", by.y="patch_id")[order(merge(merge(data.frame(id=1:length(E(g)), from_p=g.es["from_p, to_p=g.es["to_p), cl_membership, by.x="from_p", by.y="patch_id"), cl_membership, by.x="to_p", by.y="patch_id")$id) ,] - g.es["cl_pc = ifelse(cl_pc_pre$cl.x==cl_pc_pre$cl.y,0,1) - - gscript.verbose("Starting analysis on vertex level...") + g.es["cl_pc"] = ifelse(cl_pc_pre$cl.x==cl_pc_pre$cl.y,0,1) + """ + gscript.verbose(_("Starting analysis on vertex level...")) ######################################################################## ### Analysis on vertex level @@ -844,52 +959,64 @@ def main(): # g.vs["deg_dir = as.integer(degree(g, v=V(g), mode=c("in"))) # On the entire undirected graph - g_ud.vs["deg_u = degree(g_ud, v=V(g_ud)) + g_ud.vs["deg_u"] = degree(g_ud, v=g_ud.vs) # On the undirected graph with edges shorter cost distance threshold - g_ud_cd.vs["deg_uc = degree(g_ud_cd, v=V(g_ud_cd)) + g_ud_cd.vs["deg_uc"] = degree(g_ud_cd, v=g_ud_cd.vs) # On the undirected graph with only direct edges - g_ud_d.vs["deg_ud = degree(g_ud_d, v=V(g_ud_d)) + g_ud_d.vs["deg_ud"] = degree(g_ud_d, v=g_ud_d.vs) # On the undirected graph with only direct edges shorter cost distance threshold - g_ud_d_cd.vs["deg_udc = degree(g_ud_d_cd, v=V(g_ud_d_cd)) + g_ud_d_cd.vs["deg_udc"] = degree(g_ud_d_cd, v=g_ud_d_cd.vs) ######################################################################## ###### Calculate closeness centrality - g_ud_d.vs["cd_cl_ud = closeness(g_ud_d, vids=V(g_ud_d), weights = g_ud_d.es["cd_u, normalized = False) - g_ud_d.vs["mf_cl_ud = closeness(g_ud_d, vids=V(g_ud_d), weights = g_ud_d.es["mf_inv_u, normalized = False) - g_ud_d.vs["cf_cl_ud = closeness(g_ud_d, vids=V(g_ud_d), weights = g_ud_d.es["cf_inv_u, normalized = False) + g_ud_d.vs["cd_cl_ud"] = closeness( + g_ud_d, vids=g_ud_d.vs, weights=g_ud_d.es["cd_u"], normalized=False + ) + g_ud_d.vs["mf_cl_ud"] = closeness( + g_ud_d, vids=g_ud_d.vs, weights=g_ud_d.es["mf_inv_u"], normalized=False + ) + g_ud_d.vs["cf_cl_ud"] = closeness( + g_ud_d, vids=g_ud_d.vs, weights=g_ud_d.es["cf_inv_u"], normalized=False + ) - g_ud_d_cd.vs["cd_cl_udc = closeness(g_ud_d_cd, vids=V(g_ud_d_cd), weights = E(g_ud_d_cd)$cd_u, normalized = False) - g_ud_d_cd.vs["mf_cl_udc = closeness(g_ud_d_cd, vids=V(g_ud_d_cd), weights = E(g_ud_d_cd)$mf_inv_u, normalized = False) - g_ud_d_cd.vs["cf_cl_udc = closeness(g_ud_d_cd, vids=V(g_ud_d_cd), weights = E(g_ud_d_cd)$cf_inv_u, normalized = False) + g_ud_d_cd.vs["cd_cl_udc"] = closeness( + g_ud_d_cd, vids=g_ud_d_cd.vs, weights=g_ud_d_cd.es["cd_u"], normalized=False + ) + g_ud_d_cd.vs["mf_cl_udc"] = closeness( + g_ud_d_cd, vids=g_ud_d_cd.vs, weights=g_ud_d_cd["mf_inv_u"], normalized=False + ) + g_ud_d_cd.vs["cf_cl_udc"] = closeness( + g_ud_d_cd, vids=g_ud_d_cd.vs, weights=g_ud_d_cd["cf_inv_u"], normalized=False + ) ######################################################################## ###### Calculate biconnected components # On the undirected graph with only direct edges # Number of new clusters when a vertex is deleted - g_ud_d.vs["art_ud = unlist(mclapply(mc.cores={0}, 1:(length(V(g_ud_d))), function(x) ifelse((clusters(delete.vertices(g_ud_d, V(g_ud_d)[x]))$no-cl_no_d)<0,0,clusters(delete.vertices(g_ud_d, V(g_ud_d)[x]))$no-cl_no_d))) + # g_ud_d.vs["art_ud"] = unlist(mclapply(mc.cores={0}, 1:(length(V(g_ud_d))), function(x) ifelse((clusters(delete.vertices(g_ud_d, V(g_ud_d)[x]))$no-cl_no_d)<0,0,clusters(delete.vertices(g_ud_d, V(g_ud_d)[x]))$no-cl_no_d))) # Vertex is articulation point - g_ud_d.vs["art_p_ud = 0 - g_ud_d.vs["art_p_ud[biconnected_components_d$articulation_points] = 1 + g_ud_d.vs["art_p_ud"] = 0 + g_ud_d.vs[art_p_ud[biconnected_components_d["articulation_points"]]] = 1 # On the undirected graph with only direct edges shorter cost distance threshold # Number of new clusters when a vertex is deleted - g_ud_d_cd.vs["art_udc = unlist(mclapply(mc.cores={0}, 1:(length(V(g_ud_d_cd))), function(x) ifelse((clusters(delete.vertices(g_ud_d_cd, V(g_ud_d_cd)[x]))$no-cl_no_d_cd)<0,0,clusters(delete.vertices(g_ud_d_cd, V(g_ud_d_cd)[x]))$no-cl_no_d_cd))) + # g_ud_d_cd.vs["art_udc"] = unlist(mclapply(mc.cores={0}, 1:(length(V(g_ud_d_cd))), function(x) ifelse((clusters(delete.vertices(g_ud_d_cd, V(g_ud_d_cd)[x]))$no-cl_no_d_cd)<0,0,clusters(delete.vertices(g_ud_d_cd, V(g_ud_d_cd)[x]))$no-cl_no_d_cd))) # Vertex is articulation point - g_ud_d_cd.vs["art_p_udc = 0 - g_ud_d_cd.vs["art_p_udc[biconnected_components_d_cd$articulation_points] = 1 + g_ud_d_cd.vs["art_p_udc"] = 0 + g_ud_d_cd.vs["art_p_udc"][biconnected_components_d_cd["articulation_points"]] = 1 ######################################################################## ###### Calculate eigenvector centrality # Calculate eigenvector for the entire directed graph - + """ # Calculate eigenvector centrality weighted by competing potential flow cf_evc_groups = groupedData(cf ~ 1 | to_p, data=e, FUN=sum) cf_evc = gsummary(cf_evc_groups, sum) - df_cf_evc_d = data.frame(patch_id=cf_evc$to_p, cf_evc_d=cf_evc$cf) + df_cf_evc_d = data.frame(patch_id=cf_evc["to_p"], cf_evc_d=cf_evc["cf"]) # Calculate eigenvector centrality weighted by maximum potential flow mf_o_evc_groups = groupedData(mf_o ~ 1 | to_p, data=e, FUN=sum) @@ -1005,12 +1132,13 @@ def main(): # Write dataframe dbWriteTable(con, network_output, network_overview_measures, overwrite=overwrite, row.names=False) + """ ############################################################### ##### Merge vertexmeasures in a dataframe and save it ### Create initial export data frame - vertex_export_df_ud = as.data.frame(1:length(V(g_ud))) + # vertex_export_df_ud = as.data.frame(1:length(V(g_ud))) #### Adjust vertex-attributes first @@ -1019,14 +1147,15 @@ def main(): vertex_attribute_list = list.vertex.attributes(g_ud) if len(vertex_attribute_list) > 0: for vl in vertex_attribute_list: - vertex_export_df_ud = as.data.frame(cbind(vertex_export_df_ud, get.vertex.attribute(g_ud, vl))) + # vertex_export_df_ud = as.data.frame(cbind(vertex_export_df_ud, get.vertex.attribute(g_ud, vl))) + continue - vertex_export_df_ud = vertex_export_df_ud[2:length(vertex_export_df_ud)] - names(vertex_export_df_ud) = vertex_attribute_list + vertex_export_df_ud = vertex_export_df_ud[2 : length(vertex_export_df_ud)] + vertex_export_df_ud["names"] = vertex_attribute_list ### Create initial export data frame - vertex_export_df_ud_d = as.data.frame(1:length(V(g_ud_d))) + # vertex_export_df_ud_d = as.data.frame(1:length(V(g_ud_d))) #### Adjust vertex-attributes first g_ud_d = remove.vertex.attribute(g_ud_d, "pop_proxy") @@ -1034,41 +1163,44 @@ def main(): vertex_attribute_list = list.vertex.attributes(g_ud_d) if length(vertex_attribute_list) > 0: for vl in vertex_attribute_list: - vertex_export_df_ud_d = as.data.frame(cbind(vertex_export_df_ud_d, get.vertex.attribute(g_ud_d, vl))) + # vertex_export_df_ud_d = as.data.frame(cbind(vertex_export_df_ud_d, get.vertex.attribute(g_ud_d, vl))) + continue - vertex_export_df_ud_d = vertex_export_df_ud_d[2:length(vertex_export_df_ud_d)] - names(vertex_export_df_ud_d) = vertex_attribute_list + # vertex_export_df_ud_d = vertex_export_df_ud_d[2:length(vertex_export_df_ud_d)] + vertex_export_df_ud_d["names"] = vertex_attribute_list ###Create initial export data frame - vertex_export_df_ud_d_cd = as.data.frame(1:length(V(g_ud_d_cd))) + # vertex_export_df_ud_d_cd = as.data.frame(1:length(V(g_ud_d_cd))) ####Adjust vertex-attributes first g_ud_d_cd = remove.vertex.attribute(g_ud_d_cd, "pop_proxy") vertex_attribute_list = list.vertex.attributes(g_ud_d_cd) - if len(vertex_attribute_list)>0: + if len(vertex_attribute_list) > 0: for vl in vertex_attribute_list: - vertex_export_df_ud_d_cd = as.data.frame(cbind(vertex_export_df_ud_d_cd, get.vertex.attribute(g_ud_d_cd, vl))) + # vertex_export_df_ud_d_cd = as.data.frame(cbind(vertex_export_df_ud_d_cd, get.vertex.attribute(g_ud_d_cd, vl))) + continue - vertex_export_df_ud_d_cd = vertex_export_df_ud_d_cd[2:length(vertex_export_df_ud_d_cd)] - names(vertex_export_df_ud_d_cd) = vertex_attribute_list + # vertex_export_df_ud_d_cd = vertex_export_df_ud_d_cd[2:length(vertex_export_df_ud_d_cd)] + vertex_export_df_ud_d_cd["names"] = vertex_attribute_list ###Create initial export data frame - vertex_export_df_d = as.data.frame(1:length(V(g))) + # vertex_export_df_d = as.data.frame(1:length(V(g))) vertex_attribute_list = list.vertex.attributes(g) - if len(vertex_attribute_list)>0: + if len(vertex_attribute_list) > 0: for vl in vertex_attribute_list: - vertex_export_df_d = as.data.frame(cbind(vertex_export_df_d, get.vertex.attribute(g, vl))) + # vertex_export_df_d = as.data.frame(cbind(vertex_export_df_d, get.vertex.attribute(g, vl))) + continue - vertex_export_df_d = vertex_export_df_d[2:length(vertex_export_df_d)] - names(vertex_export_df_d) = vertex_attribute_list + # vertex_export_df_d = vertex_export_df_d[2:length(vertex_export_df_d)] + vertex_export_df_d["names"] = vertex_attribute_list - vertex_export_df = merge(merge(merge(vertex_export_df_d, vertex_export_df_ud, by="patch_id"), vertex_export_df_ud_d, by="patch_id"), vertex_export_df_ud_d_cd, by="patch_id") - dbWriteTable(con, vertex_output, vertex_export_df, overwrite=overwrite, row.names=False) + # vertex_export_df = merge(merge(merge(vertex_export_df_d, vertex_export_df_ud, by="patch_id"), vertex_export_df_ud_d, by="patch_id"), vertex_export_df_ud_d_cd, by="patch_id") + # dbWriteTable(con, vertex_output, vertex_export_df, overwrite=overwrite, row.names=False) ######################################################################### ###Create initial export data frame for the undirected graph - edge_export_df_ud = as.data.frame(g_ud.es["con_id_u) + # edge_export_df_ud = as.data.frame(g_ud.es["con_id_u) ####Adjust edge attributes for the undirected graphs first g_ud = remove.edge.attribute(g_ud, "con_id") @@ -1080,7 +1212,7 @@ def main(): g_ud = remove.edge.attribute(g_ud, "cost_distance") g_ud = remove.edge.attribute(g_ud, "cd_u") g_ud = remove.edge.attribute(g_ud, "distance_weight_e") - #g_ud = remove.edge.attribute(g_ud, "distance_weight_e_ud") + # g_ud = remove.edge.attribute(g_ud, "distance_weight_e_ud") g_ud = remove.edge.attribute(g_ud, "mf_o") g_ud = remove.edge.attribute(g_ud, "mf_i") g_ud = remove.edge.attribute(g_ud, "mf_o_inv") @@ -1094,13 +1226,14 @@ def main(): edge_attribute_list = list.edge.attributes(g_ud) - if length(edge_attribute_list )>0: + if length(edge_attribute_list) > 0: for el in edge_attribute_list: - edge_export_df_ud = as.data.frame(cbind(edge_export_df_ud, get.edge.attribute(g_ud, el))) - names(edge_export_df_ud) = append("con_id_u", edge_attribute_list) + # edge_export_df_ud = as.data.frame(cbind(edge_export_df_ud, get.edge.attribute(g_ud, el))) + continue + edge_export_df_ud["names"] = append("con_id_u", edge_attribute_list) ###Create initial export data frame for the undirected graph with only direct edges - edge_export_df_ud_d = as.data.frame(g_ud_d.es["con_id_u) + # edge_export_df_ud_d = as.data.frame(g_ud_d.es["con_id_u) ####Adjust edge attributes for the undirected graph with only direct edges g_ud_d = remove.edge.attribute(g_ud_d, "con_id") @@ -1112,7 +1245,7 @@ def main(): g_ud_d = remove.edge.attribute(g_ud_d, "cost_distance") g_ud_d = remove.edge.attribute(g_ud_d, "cd_u") g_ud_d = remove.edge.attribute(g_ud_d, "distance_weight_e") - #g_ud_d = remove.edge.attribute(g_ud_d, "distance_weight_e_ud") + # g_ud_d = remove.edge.attribute(g_ud_d, "distance_weight_e_ud") g_ud_d = remove.edge.attribute(g_ud_d, "mf_o") g_ud_d = remove.edge.attribute(g_ud_d, "mf_i") g_ud_d = remove.edge.attribute(g_ud_d, "mf_o_inv") @@ -1130,14 +1263,15 @@ def main(): edge_attribute_list = list.edge.attributes(g_ud_d) - if len(edge_attribute_list )>0: + if len(edge_attribute_list) > 0: for el in edge_attribute_list: - edge_export_df_ud_d = as.data.frame(cbind(edge_export_df_ud_d, get.edge.attribute(g_ud_d, el))) + # edge_export_df_ud_d = as.data.frame(cbind(edge_export_df_ud_d, get.edge.attribute(g_ud_d, el))) + pass - names(edge_export_df_ud_d) = append("con_id_u", edge_attribute_list) + edge_export_df_ud_d["names"] = append("con_id_u", edge_attribute_list) ###Create initial export data frame for the undirected graph with only direct edges - edge_export_df_ud_d_cd = as.data.frame(E(g_ud_d_cd)$con_id_u) + # edge_export_df_ud_d_cd = as.data.frame(E(g_ud_d_cd)$con_id_u) ####Adjust edge attributes for the undirected graph with only direct edges g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "con_id") @@ -1149,7 +1283,7 @@ def main(): g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "cost_distance") g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "cd_u") g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "distance_weight_e") - #g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "distance_weight_e_ud") + # g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "distance_weight_e_ud") g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "mf_o") g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "mf_i") g_ud_d_cd = remove.edge.attribute(g_ud_d_cd, "mf_o_inv") @@ -1167,190 +1301,281 @@ def main(): edge_attribute_list = list.edge.attributes(g_ud_d_cd) - if len(edge_attribute_list )>0: + if len(edge_attribute_list) > 0: for el in edge_attribute_list: - edge_export_df_ud_d_cd = as.data.frame(cbind(edge_export_df_ud_d_cd, get.edge.attribute(g_ud_d_cd, el))) - names(edge_export_df_ud_d_cd) = append("con_id_u", edge_attribute_list) + # edge_export_df_ud_d_cd = as.data.frame(cbind(edge_export_df_ud_d_cd, get.edge.attribute(g_ud_d_cd, el))) + pass + edge_export_df_ud_d_cd["names"] = append("con_id_u", edge_attribute_list) ###Create initial export data frame for the directed graph - edge_export_df = as.data.frame(1:length(E(g))) + # edge_export_df = as.data.frame(1:length(E(g))) ####Adjust edge attributes for the directed graph first - g.es["cd"] = g.es["cost_distance - g.es["distk"] = g.es["distance_weight_e - g.es["distk_u"] = g.es["distance_weight_e_ud + g.es["cd"] = g.es["cost_distance"] + g.es["distk"] = g.es["distance_weight_e"] + g.es["distk_u"] = g.es["distance_weight_e_ud"] g = remove.edge.attribute(g, "distance_weight_e") g = remove.edge.attribute(g, "distance_weight_e_ud") g = remove.edge.attribute(g, "cost_distance") edge_attribute_list = list.edge.attributes(g) - if len(edge_attribute_list )>0: + if len(edge_attribute_list) > 0: for el in edge_attribute_list: - edge_export_df = as.data.frame(cbind(edge_export_df, get.edge.attribute(g, el))) + # edge_export_df = as.data.frame(cbind(edge_export_df, get.edge.attribute(g, el))) + pass - names(edge_export_df) = append("id", edge_attribute_list) + edge_export_df["names"] = append("id", edge_attribute_list) - export_df_list =c("edge_export_df_ud", "edge_export_df_ud_d", "edge_export_df_ud_d_cd") + export_df_list = ( + "edge_export_df_ud", + "edge_export_df_ud_d", + "edge_export_df_ud_d_cd", + ) for df in export_df_list: if len(names(get(df))): - edge_export_df_final = merge(edge_export_df, get(df), all.x=True, by="con_id_u", suffixes=c("_x", "_y")) + edge_export_df_final = merge( + edge_export_df, + get(df), + all_x=True, + by="con_id_u", + suffixes=("_x", "_y"), + ) edge_export_df = edge_export_df_final - dbWriteTable(con, edge_output, edge_export_df, overwrite=overwrite, row.names=False) + dbWriteTable(con, edge_output, edge_export_df, overwrite=overwrite, row_names=False) dbDisconnect(con) - if qml_directory != '': + if qml_directory != "": ######################### - #CREATE .qml-files for edge measures visualistion in QGIS + # CREATE .qml-files for edge measures visualistion in QGIS no_quantile = 5 - colortable = c(' ', - ' ', - ' ', - ' ', - ' ') - colortable_bin = c(' ') - for (attribute in names(edge_export_df_final)) {{ - - #Skip id and geom columns - if(attribute %in% c("id", "con_id", "con_id_u", "from_p", "to_p", "WKT", "cf_ebc_vi")) {{ next }} - - st_mod = storage.mode(unlist(edge_export_df_final[grep(1, match(names(edge_export_df_final), attribute))])) - att_val = unlist(edge_export_df_final[grep(1, match(names(edge_export_df_final), attribute))]) - - #Write header + colortable = ( + ' ', + ' ', + ' ', + ' ', + ' ', + ) + colortable_bin = ' ' + for attribute in edge_export_df_final: + + # Skip id and geom columns + if attribute in ( + "id", + "con_id", + "con_id_u", + "from_p", + "to_p", + "WKT", + "cf_ebc_vi", + ): + continue + + st_mod = storage.mode( + unlist( + edge_export_df_final[ + grep(1, match(names(edge_export_df_final), attribute)) + ] + ) + ) + att_val = unlist( + edge_export_df_final[ + grep(1, match(names(edge_export_df_final), attribute)) + ] + ) + + # Write header qml = "" qml += '' - qml += ' 255' - - if max(att_val, na.rm=True)-min(att_val, na.rm=True)==1: - - #More header - qml += paste(' ', sep='')) - - #Categories - qml += ' ') - qml += ' ') - qml += ' ') - - #Write symbols - qml += ' ') - qml += ' ') - qml += ' ') - qml += ' ') - qml += colortable_bin) - qml += ' ') - qml += ' ') - qml += ' ') - qml += ' ') - qml += ' ') - qml += ' ') - qml += ' ') - qml += ' ') + qml += " 255" + + if max(att_val, na_rm=True) - min(att_val, na_rm=True) == 1: + + # More header + qml += paste( + ' ', + sep="", + ) + + # Categories + qml += " " + qml += ' ' + qml += " " + + # Write symbols + qml += " " + qml += ' ' + qml += ' ' + qml += ' ' + qml += colortable_bin + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += " " + qml += " " else: - attribute_quantile = quantile(edge_export_df_final[grep(1, match(names(edge_export_df_final), attribute))], probs=seq(0, 1, 1/5), na.rm=True, type=8) - - #Write more header - qml += paste(' ', sep='')) - - #Write ranges - ranges = character() - qml += ' ' - for quant in 1:no_quantile: - ranges = append(ranges, paste(' ', sep='') + # attribute_quantile = quantile(edge_export_df_final[grep(1, match(names(edge_export_df_final), attribute))], probs=seq(0, 1, 1/5), na.rm=True, type=8) + + # Write more header + qml += paste( + ' ', + sep="", + ) + + # Write ranges + ranges = str() + qml += " " + for quant in range(no_quantile): + ranges = append( + ranges, + paste( + ' ', + sep="", + ), + ) qml += ranges - qml += ' ' - - #Write symbols - qml += ' ' - - for quant in 1:no_quantile: - qml += paste(' ', - sep='') - qml += paste(' ', - sep='') - qml += ' ' + qml += " " + + # Write symbols + qml += " " + + for quant in range(no_quantile): + qml += paste( + ' ', + sep="", + ) + qml += paste( + ' ', + sep="", + ) + qml += ' ' qml += colortable[quant] - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' - - #Write Footer - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' - qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += " " + qml += " " + + # Write Footer + qml += " " + qml += " " + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += ' ' + qml += " " + qml += " " + qml += " " qml += ' ' qml += ' ' qml += ' ' - qml += ' ' - qml += ' ' - qml += paste(' "', attribute, '"', sep='') - qml += ' ' - qml += '' - - #Save qml-file - qml_output = paste(qml_directory, paste(paste("edge_measures_", attribute, sep=''), "qml", sep='.'), sep="/") + qml += " " + qml += " " + qml += ' "' + attribute + '"' + qml += " " + qml += "" + + # Save qml-file + qml_output = paste( + qml_directory, + paste(paste("edge_measures_", attribute, sep=""), "qml", sep="."), + sep="/", + ) con_qml = file(qml_output, open="wt") - write.table(qml, con_qml, append = False, quote = False, sep = " ", eol = "\n", na = "NA", dec = ".", row.names = False, col.names = False) + write.table( + qml, + con_qml, + append=False, + quote=False, + sep=" ", + eol="\n", + na="NA", + dec=".", + row_names=False, + col_names=False, + ) close(con_qml) - grass.run_command('g.copy', quiet=True, - vector='{},{}'.format(network_map, - edge_output)) - grass.run_command('g.copy', quiet=True, - vector='{},{}'.format(in_vertices, - vertex_output)) + gscript.run_command( + "g.copy", quiet=True, vector="{},{}".format(network_map, edge_output) + ) + gscript.run_command( + "g.copy", quiet=True, vector="{},{}".format(in_vertices, vertex_output) + ) # Use v.db.connect instead of v.db.join (much faster) - grass.run_command('v.db.join', map=edge_output, column='cat', - other_table=edge_output_tmp, - other_column='con_id', quiet=True) - grass.run_command('v.db.join', map=vertex_output, column='cat', - other_table=vertex_output_tmp, - other_column='patch_id', quiet=True) + gscript.run_command( + "v.db.join", + map=edge_output, + column="cat", + other_table=edge_output_tmp, + other_column="con_id", + quiet=True, + ) + gscript.run_command( + "v.db.join", + map=vertex_output, + column="cat", + other_table=vertex_output_tmp, + other_column="patch_id", + quiet=True, + ) - update_history = '{}\n{}'.format(net_hist_str, - os.environ['CMDLINE']) + update_history = "{}\n{}".format(net_hist_str, os.environ["CMDLINE"]) - grass.run_command('v.support', flags='h', map=vertex_output, - person=os.environ['USER'], - cmdhist=update_history) + gscript.run_command( + "v.support", + flags="h", + map=vertex_output, + person=os.environ["USER"], + cmdhist=update_history, + ) - grass.run_command('v.support', flags='h', map=edge_output, - person=os.environ['USER'], - cmdhist=update_history) + gscript.run_command( + "v.support", + flags="h", + map=edge_output, + person=os.environ["USER"], + cmdhist=update_history, + ) if __name__ == "__main__": - options, flags = grass.parser() + options, flags = gscript.parser() atexit.register(cleanup) - sys.exit(main()) \ No newline at end of file + sys.exit(main())