1
2context("Graphlets")
3
4sortgl <- function(x) {
5  cl <- lapply(x$cliques, sort)
6  n <- sapply(cl, length)
7  list(cliques=cl[order(n)], thresholds=x$thresholds[order(n)])
8}
9
10test_that("Graphlets work for some simple graphs", {
11  library(igraph)
12
13  g <- make_full_graph(5)
14  E(g)$weight <- 1
15  gl <- graphlet_basis(g)
16
17  expect_that(names(gl), equals(c("cliques", "thresholds")))
18  expect_that(length(gl$cliques), equals(1))
19  expect_that(sort(gl$cliques[[1]]), equals(1:vcount(g)))
20  expect_that(gl$thresholds, equals(1))
21
22  g2 <- make_full_graph(5)
23  E(g2)$weight <- 1
24  E(g2)[1%--%2]$weight <- 2
25  gl2 <- sortgl(graphlet_basis(g2))
26
27  expect_that(gl2, equals(list(cliques=list(1:2, 1:5), thresholds=c(2,1))))
28})
29
30test_that("Graphlets filtering works", {
31  library(igraph)
32  gt <- data.frame(from  =c("A", "A", "B", "B", "B", "C", "C", "D"),
33                   to    =c("B", "C", "C", "D", "E", "D", "E", "E"),
34                   weight=c( 8 ,  8 ,  8 ,  5 ,  5 ,  5 ,  5 ,  5 ))
35
36  g <- graph_from_data_frame(gt, directed=FALSE, vertices=data.frame(LETTERS[1:5]))
37  gl <- sortgl(graphlet_basis(g))
38
39  expect_that(gl$cliques, equals(list(1:3, 2:5)))
40  expect_that(gl$thresholds, equals(c(8, 5)))
41})
42
43## Naive version of graphlets
44
45unvs <- function(x) lapply(x, as.vector)
46
47threshold.net <- function(graph, level) {
48  N <- vcount(graph)
49  graph.t <- delete_edges(graph, which(E(graph)$weight < level))
50
51  clqt <- unvs(max_cliques(graph.t))
52  clqt <- lapply(clqt, sort)
53  clqt[order(sapply(clqt, length), decreasing=TRUE)]
54}
55
56graphlets.old <- function(graph) {
57
58  if (!is_weighted(graph)) { stop("Graph not weighted") }
59  if (min(E(graph)$weight) <= 0 || any(!is.finite(E(graph)$weight))) {
60    stop("Edge weights must be non-negative and finite")
61  }
62
63  ## Do all thresholds
64  cl <- lapply(sort(unique(E(graph)$weight)), function(w) {
65    threshold.net(graph, w)
66  })
67
68  ## Put the cliques in one long list
69  clv <- unlist(cl, recursive=FALSE)
70
71  ## Sort the vertices within the cliques
72  cls <- lapply(clv, sort)
73
74  ## Delete duplicate cliques
75  clu <- unique(cls)
76
77  ## Delete cliques that consist of single vertices
78  clf <- clu[sapply(clu, length) != 1]
79
80  clf
81}
82
83test_that("Graphlets work for a bigger graph", {
84  library(igraph)
85  set.seed(42)
86  g <- make_graph("zachary")
87  E(g)$weight <- sample(1:5, ecount(g), replace=TRUE)
88
89  gl <- graphlet_basis(g)
90  gl2 <- graphlets.old(g)
91
92  glo <- sort(sapply(unvs(gl$cliques), paste, collapse="-"))
93  gl2o <- sort(sapply(gl2, paste, collapse="-"))
94
95  expect_that(glo, equals(gl2o))
96})
97
98graphlets.project.old <- function(graph, cliques, iter, Mu=NULL) {
99
100  if (!is_weighted(graph)) { stop("Graph not weighted") }
101  if (min(E(graph)$weight) <= 0 || any(!is.finite(E(graph)$weight))) {
102    stop("Edge weights must be non-negative and finite")
103  }
104  if (length(iter) != 1 || !is.numeric(iter) ||
105      !is.finite(iter) || iter != as.integer(iter)) {
106    stop("`iter' must be a non-negative finite integer scalar")
107  }
108
109  clf <- cliques
110
111  ## Create vertex-clique list first
112  vcl <- vector(length=vcount(graph), mode="list")
113  for (i in 1:length(clf)) {
114    for (j in clf[[i]]) {
115      vcl[[j]] <- c(vcl[[j]], i)
116    }
117  }
118
119  ## Create edge-clique list from this, it is useful to have the edge list
120  ## of the graph at hand
121  el <- as_edgelist(graph, names=FALSE)
122  ecl <- vector(length=ecount(graph), mode="list")
123  for (i in 1:ecount(graph)) {
124    edge <- el[i,]
125    ecl[[i]] <- intersect(vcl[[edge[1]]], vcl[[edge[2]]])
126  }
127
128  ## We will also need a clique-edge list, the edges in the cliques
129  system.time({
130    cel <- vector(length=length(clf), mode="list")
131    for (i in 1:length(ecl)) {
132      for (j in ecl[[i]]) {
133        cel[[j]] <- c(cel[[j]], i)
134      }
135    }
136  })
137
138  ## OK, we are ready to do the projection now
139  if (is.null(Mu)) { Mu <- rep(1, length(clf)) }
140  origw <- E(graph)$weight
141  w <- numeric(length(ecl))
142  a <- sapply(clf, function(x) length(x) * (length(x) + 1) / 2)
143  for (i in 1:iter) {
144    for (j in 1:length(ecl)) {
145      w[j] <- sum(Mu[ ecl[[j]] ])
146    }
147    for (j in 1:length(clf)) {
148      Mu[j] <- Mu[j] * sum(origw[cel[[j]]] / (w[cel[[j]]] + .0001)) / a[j]
149    }
150  }
151
152  ## Sort the cliques according to their weights
153  Smb <- sort(Mu, decreasing=TRUE, index=TRUE)
154  list(cliques=clf[Smb$ix], Mu=Mu[Smb$ix])
155}
156
157test_that("Graphlet projection works", {
158  library(igraph)
159
160  D1 <- matrix(0, 5, 5)
161  D2 <- matrix(0, 5, 5)
162  D3 <- matrix(0, 5, 5)
163  D1[1:3, 1:3] <- 2
164  D2[3:5, 3:5] <- 3
165  D3[2:5, 2:5] <- 1
166  g <- graph_from_adjacency_matrix(D1 + D2 + D3, mode="undirected", weighted=TRUE)
167  g <- simplify(g)
168
169  gl <- graphlet_basis(g)
170  glp <- graphlets(g)
171  glp2 <- graphlets.project.old(g, cliques=gl$cliques, iter=1000)
172
173  glp$cliques <- unvs(glp$cliques)
174  expect_that(glp, equals(glp2))
175})
176