fcm.R

# 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")
}

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!