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")
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment