Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
# Copyright 2015 JP de Vooght <jp@vooght.de>
#
# Permission to use, copy, modify, and distribute this software for any
# purpose with or without fee is hereby granted, provided that the above
# notice and this permission notice appear in all copies.
#
# THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
# WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
# MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
# ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
# WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
# ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
# OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
suppressPackageStartupMessages(library(igraph))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(reshape2))
suppressPackageStartupMessages(library(scales))
# use this function first to check the consistency
# of the fuzzy cognitive map
# maps are provided in GML format as generated by yEd
# returns the adjacency matrix
# yEd can be found at http://www.yWorks.com/
check.map = function(filename) {
if(!file.exists(filename)) {
stop("File does not exist.")
}
# all edges must have labels
g = read.graph(filename, format="gml")
if(any(E(g)$label == "")) {
p = paste(which(E(g)$label == ""), collapse=",")
stop(paste("Missing weights in map", p))
}
# all edge labels are valid numeric values
E(g)$weight = suppressWarnings(as.numeric(E(g)$label))
if(any(is.na(E(g)$weight))) {
p = paste(which(is.na(E(g)$weight)), collapse=",")
stop(paste("Failed to read weight value", p))
}
# no edges have weights above 1
if(any(E(g)$weight > 1)) {
p = paste(which(E(g)$weight > 1), collapse=",")
stop(paste("Found weights greater than 1", p))
}
# no edges have weights below -1
if(any(E(g)$weight < -1)) {
p = paste(which(E(g)$weight < -1), collapse=",")
stop(paste("Found weights smaller than -1", p))
}
# create the adjacency matrix
mat = as.matrix(as_adjacency_matrix(g, attr="weight"))
if(any(diag(mat) > 0)) {
p = paste(which(diag(mat) > 0), collapse=",")
warning(paste("Diagonal contains non-zero elements", p))
}
return(g)
}
# computes map-level indices
# returns a data frame of indices in one row
# See Oezesmi and Oezesmi, Ecological Modelling 176 (2004) 43-64
map.indices = function(filename) {
if(!file.exists(filename)) {
stop("File does not exist.")
}
g = read.graph(filename, format="gml")
E(g)$weight = suppressWarnings(abs(as.numeric(E(g)$label)))
mat = as.matrix(as_adjacency_matrix(g, attr="weight"))
Connections = ecount(g)
Density = graph.density(g, loops=any(diag(mat) > 0))
Concepts = vcount(g)
ConnectionsPerVariable = Connections / Concepts
Transmitters = length(
intersect(which(strength(g, mode="out") > 0),
which(strength(g, mode="in") == 0)))
Receivers = length(
intersect(which(strength(g, mode="in") > 0),
which(strength(g, mode="out") == 0)))
Complexity = Receivers / Transmitters
Isolates = length(which(strength(g) == 0))
Ordinary = Concepts - Transmitters - Receivers - Isolates
SelfLoops = sum(diag(mat) != 0)
Hierarchy = var(strength(g, mode="out")) *
12/(vcount(g)*(vcount(g)-1)*(vcount(g)+1))
res = data.frame(Connections, Density, Concepts, Transmitters,
Receivers, Isolates, Ordinary, SelfLoops,
ConnectionsPerVariable, Complexity, Hierarchy)
colnames(res) = c("Connections", "Density", "Concepts",
"Transmitters", "Receivers", "Isolates",
"Ordinary", "Self-loops", "CPV", "Complexity",
"Hierarchy")
return(res)
}
# computes concept-level indices
# returns a data frame where nrows=concepts and ncols=indices
# side-effect: generates bar chart of concepts ordered by centrality
# See Oezesmi and Oezesmi, Ecological Modelling 176 (2004) 43-64
concept.indices = function(filename) {
if(!file.exists(filename)) {
stop("File does not exist.")
}
g = read.graph(filename, format="gml")
E(g)$weight = suppressWarnings(abs(as.numeric(E(g)$label)))
Outdegree = strength(g, mode="out")
Indegree = strength(g, mode="in")
Centrality = Outdegree + Indegree
Concept = vcount(g)
Transmitter = replace(rep(F, Concept),
intersect(
which(strength(g, mode="out") > 0),
which(strength(g, mode="in") == 0)),
T)
Receiver = replace(rep(F, Concept),
intersect(
which(strength(g, mode="in") > 0),
which(strength(g, mode="out") == 0)),
T)
Ordinary = replace(rep(F, Concept),
intersect(
which(strength(g, mode="out") > 0),
which(strength(g, mode="in") > 0)),
T)
Isolates = replace(rep(F, Concept), which(strength(g) == 0), T)
Index = c("Concept", "Outdegree", "Indegree", "Centrality",
"Transmitter", "Receiver", "Ordinary", "Isolates")
res = data.frame(V(g)$label, Outdegree, Indegree,
Centrality, Transmitter, Receiver,
Ordinary, Isolates)
colnames(res) = Index
dat = melt(res[,c("Concept", "Outdegree", "Indegree")], id="Concept")
print(ggplot(dat, aes(x=reorder(Concept, -value), y=value, fill=variable)) +
geom_bar(stat="identity") +
theme(legend.position="top", legend.direction="horizontal",
axis.text.x=element_text(angle=45, hjust=1)) +
labs(x="Concept", y="Centrality", fill=""))
return(res)
}
# linear activation function
fcm.lin = function(x) {
return(x)
}
# sigmoid activation function
fcm.sigmoid = function(x, lambda=1, h=0) {
return(1/(1+exp(-lambda*(x-h))))
}
# exponential activation function
fcm.exp = function(x, lambda=1) {
if(x < 0) {
return(-1+exp(lambda*x))
} else {
return(1-exp(-lambda*x))
}
}
# trivalent activation function
fcm.trival = function(x, lo=-1, hi=1) {
if(x < lo) {
return(-1)
} else if(x > hi) {
return(1)
} else {
return(x)
}
}
# FCM iteration function
# returns a table of (iterations x state vectors)
fcm.iterate = function(mat, init, peg, iter, trans, fn, ...) {
acc = matrix(0, iter, ncol(mat))
acc[1,] = init
if(!missing(peg)) {
acc[1,which(!is.na(peg))] = peg[!is.na(peg)]
}
for(i in 2:iter) {
if(trans) {
acc[i,] = sapply(t(acc[i-1,]) %*% mat, fn, ...)
} else {
acc[i,] = sapply(mat %*% acc[i-1,], fn, ...)
}
if(!missing(peg)) {
acc[i,which(!is.na(peg))] = peg[!is.na(peg)]
}
}
return(acc)
}
# computes baseline scenario (all initial values set to 1)
# arguments include number of iterations, trans flag and activation function
# the default case for dot products in iterations is the transposed form
# when activation function missing, assumes linear
# returns steady-state vector
# side-effect: generates line-chart of steady-state vector across iterations
baseline.scenario = function(filename, iter=20, trans=T, fn, ...) {
if(!file.exists(filename)) {
stop("File does not exist.")
}
g = read.graph(filename, format="gml")
E(g)$weight = suppressWarnings(as.numeric(E(g)$label))
mat = as_adjacency_matrix(g, attr="weight")
if(missing(fn)) {
fn = fcm.lin
}
acc = fcm.iterate(mat=mat, init=rep(1, vcount(g)), iter=iter, trans=trans, fn=fn, ...)
if(!isTRUE(all.equal(acc[iter,], acc[iter-1,]))) {
warning("Convergence not reached. Try increasing the number of iterations.")
}
print(ggplot(melt(acc, varnames=c("Iteration", "Concept")),
aes(x=Iteration, y=value)) +
geom_line(aes(group=Concept, colour=factor(Concept))) +
labs(x="Iteration", y="Value", colour="Concept") +
ggtitle("Baseline Scenario"))
res = data.frame(V(g)$label, acc[iter,])
colnames(res) = c("Concept", "Value")
return(res)
}
# computes alternative scenario
# in addition to baseline parameters, alt.scenario takes 2 vectors
# init is a vector of initial values for the scenario
# peg is a vector of values to be kept constant
# e.g. c(NA, NA, NA, NA, 1) would peg the 5th concept at 1 throughout computations
# returns steady-state vector
# side-effect: generates line-chart of steady-state vector across iterations
# TODO indicate which concepts are pegged in legend
alt.scenario = function(filename, init, peg, iter=20, trans=T, fn, ...) {
if(!file.exists(filename)) {
stop("File does not exist.")
}
g = read.graph(filename, format="gml")
E(g)$weight = suppressWarnings(as.numeric(E(g)$label))
mat = as_adjacency_matrix(g, attr="weight")
if(missing(init)) {
init = rep(1, vcount(g))
}
if(length(init) != vcount(g)) {
stop("Scenario vector length not equal to concepts count.")
}
if(any(init > 1) || any(init < -1)) {
stop("Scenario vector must contains numbers between -1 and 1.")
}
if(missing(peg)) {
peg = rep(NA, vcount(g))
}
if(length(peg) != vcount(g)) {
stop("Peg vector length not equal to concepts count.")
}
if(any(na.omit(peg) > 1) || any(na.omit(peg) < -1)) {
stop("Peg vector must contains numbers between -1 and 1.")
}
if(missing(fn)) {
fn = fcm.lin
}
acc = fcm.iterate(mat=mat, init=init, peg=peg, iter=iter, trans=trans, fn, ...)
if(!isTRUE(all.equal(acc[iter,], acc[iter-1,]))) {
warning("Convergence not reached. Try increasing the number of iterations.")
}
print(ggplot(melt(acc, varnames=c("Iteration", "Concept")),
aes(x=Iteration, y=value)) + ylim(-1, 1) +
geom_line(aes(group=Concept, colour=factor(Concept))) +
labs(x="Iteration", y="Value", colour="Concept") +
ggtitle(bquote(atop("Scenario", atop(.(paste(paste(init, collapse=", "))))))))
res = data.frame(V(g)$label, acc[iter,])
colnames(res) = c("Concept","Value")
return(res)
}
# computes differences with baseline scenario
# returns data frame including change and difference in percent
# side effect: generates bar chart of differences across scenarios
diff.scenarios = function(baseline, ...) {
vecs = list(...)
base = deparse(substitute(baseline))
dots = sapply(substitute(list(...))[-1], deparse)
if(length(vecs) == 0) {
stop("Missing alternative scenario(s).")
}
for(vec in vecs) {
if(paste(vec$Concept,collapse="") != paste(vec$Concept,collapse="")) {
stop("Scenarios do not share the same concepts.")
}
}
res = data.frame(baseline)
colnames(res) = c("Concept", base)
for(i in seq_along(vecs)) {
diff = vecs[[i]][,2] - baseline[,2]
chng = (diff / baseline[,2]) * 100
res[,c(dots[i],
paste("Diff",dots[i],sep="_"),
paste("Chng",dots[i],sep="_"))] = c(vecs[[i]][,2], diff, chng)
}
dat = melt(res[,c("Concept", paste("Diff",dots,sep="_"))], id="Concept")
p = ggplot(dat, aes(x=Concept, y=value, fill=variable)) +
geom_bar(stat="identity", position="dodge") +
theme(legend.position="top", legend.direction="horizontal",
axis.text.x=element_text(angle=45, hjust=1)) +
labs(x="Concept", y=paste("Difference from ", base), fill="") +
scale_fill_discrete(labels=dots)
suppressWarnings(print(p))
return(res)
}
# plot fuzzy cognitive map
# TODO improve layout
draw.map = function(filename) {
if(!file.exists(filename)) {
stop("File does not exist.")
}
radian.rescale = function(x, start=0, dir=1) {
c.rotate = function(x) (x + start) %% (2 * pi) * dir
c.rotate(rescale(x, c(0, 2 * pi), range(x)))
}
g = read.graph(filename, format="gml")
E(g)$weight = suppressWarnings(as.numeric(E(g)$label))
E(g)$color = ifelse(E(g)$weight < 0, "red", "blue")
V(g)$color = "white"
ec = rep(0, ecount(g))
ec[which(is.mutual(g))] = 0.5
plot(g, main=filename, layout=layout.circle, edge.curved=ec,
vertex.label.degree=radian.rescale(x=1:vcount(g), dir=-1, start=0),
vertex.label.dist=0.3, vertex.size=1,
vertex.color="black", vertex.frame.color="black")
}

IvansaR commented Sep 13, 2017

Hi, I have been trying to use your script, but I receive an error. The gml was created in yEd

read.graph(file.choose(),format = "gml")
Error in read.graph.gml(file, ...) :
At foreign.c:1127 : Parse error in GML file, line 1 (syntax error, unexpected LISTCLOSE, expecting KEYWORD), Parse error

May you please indicate what the solution may be?

Thanks!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment