Make all graphs on up to n vertices; Number of non-isomorphic classes of such graphs

#####################################################
is_non_increasing <- function(v){
  if(length(v) == 1){return(TRUE)}
  diff <- v[1:(length(v) - 1)] - v[2:length(v)]
  all(diff >= 0)
  }

#####################################################
reduce_by_isomorphism_class <- function(adj.mats){
  # Assume undirected, ie symmetrical matrices
  # TODO - Find other simply calculated graph params
  #          from adjacency matrix to reduce # of calls
  #          to is_isomorphic_to
  require(Matrix) # doesn't rank the zero matrix correctly
  require(igraph)
  
  degree.vecs <- t(sapply(adj.mats, colSums))
  ranks       <- sapply(adj.mats, rankMatrix)
  eigs        <- t(sapply(adj.mats, function(x) eigen(x)$values))
  eigs.trunc  <- round( eigs, 3 )
  
  matrix.params <- cbind(degree.vecs, ranks, eigs.trunc)
  
  # Convert the vector of (nonincreasing) degrees,
  #   the matrix rank, and the rounded eigenvalues
  #   into a single factor to permit quick matching prior to
  #   isomorphism tests
  match.params <- data.frame(
    mat = 1:length(adj.mats),
    params = factor(
      apply(matrix.params, 1, 
      function(row){paste(row, collapse = ' ')})
      )
    )
  
  iso.classes <- lapply(levels(match.params$params),
    function(par.str){
      rows <- which(match.params$params == par.str)
      temp.classes <- list()
      temp.classes <- append(temp.classes, adj.mats[rows][1])
      if(length(rows) > 1){
        # drop adj.mats[current.row] if its is isomorphic
        #   to a graph in temp.classes
        for(r in rows){
          test.mat <- adj.mats[[r]]
          test.g <- graph.adjacency(
            test.mat,
            'undirected'
            )
          iso.to.any <- sapply(temp.classes, function(storedmat){
            stored.g <- graph.adjacency(
              storedmat,
              'undirected'
              )
            is_isomorphic_to(test.g, stored.g)
            })
          iso.to.any <- any(iso.to.any)
          if(! iso.to.any){
            temp.classes <- append(temp.classes, adj.mats[r])
            }
          }
        }
      temp.classes
    })
  do.call('c', iso.classes)
  }

#####################################################

## all graphs on n vertices
make_all_adjMats <- function(n){
  # not yet finished - works for n=1..8, but a bit slow
  maxN = 10
  stopifnot(n %in% seq(1, maxN))
  
  if(n == 1){
    return(
      list(matrix(0, nrow = n, ncol = n))
      )
    }
  adj.mats.lo <- lapply(
    # append an extra vertex to each one-smaller adjacency matrix
    make_all_adjMats(n-1),
    function(m){
      M <- matrix(0, nrow = n, ncol = n)
      M[1:(n-1), 1:(n-1)] <- m
      M
      })
  
  # all setwise combinations of 1:n-1, incl the empty set
  all.combs <- lapply(0:(n-1), function(x) combn((n-1), x))
  
  adj.mats <- lapply(adj.mats.lo, function(m){
    sublist <- list()
    # add edges for the new vertex, vertex n
    for(combmat in all.combs){
      if(nrow(combmat) == 0){
        sublist <- append(sublist, list(m))
        } else {
        for(col in 1:ncol(combmat)){
          comb <- combmat[, col]
          temp.m <- m
          temp.m[n, comb] <- 1
          temp.m[comb, n] <- 1
          # ensure deg(i)>=deg(i+1) for i=1..n-1
          if(is_non_increasing(colSums(temp.m))){
            sublist <- append(sublist, list(temp.m))
            }
          }
        }
      }
    sublist
    })
  adj.mats <- do.call('c', adj.mats)
  adj.mats <- adj.mats[order(sapply(adj.mats, sum))]
  # Becomes slow for n = 7
  # Methods to filter by isomorphism class?
  # a - ensure deg(i) >= deg(i+1) - done
  # b - group by degree vector and rank(ie, # connComp)
  #       then use is_isomorphic_to
  adj.mats <- reduce_by_isomorphism_class(adj.mats)
  adj.mats
  }

####################################################
##

make_all_graphs <- function(n){
  require(igraph)
  # let vertices be 1..n and assume deg(i)<=deg(i+1) for all i
  adj.mats <- make_all_adjMats(n)
  graphs <- lapply(
    adj.mats,
    function(m){
      graph_from_adjacency_matrix(
        adjmatrix = m,
        mode = 'undirected'
        )})
  graphs
  }


#######################################################
library(igraph)


# Number of non-isomorphic connected graphs on 7 vertices:

system.time(G7 <- Filter(is.connected, make_all_graphs(7)))

# Without filtering on eigenvalues
#    user  system elapsed 
#    3.35    0.00    3.37

# With filtering on eigenvalues
#    user  system elapsed 
#    2.40    0.00    2.41 

length(G7) # 524 

G7.eigen.adj <- t(sapply(G7, function(g){
  eigen(get.adjacency(g))$values
  }))
length(which(duplicated(G7.eigen.adj))) # 1

length(which(duplicated(round(G7.eigen.adj, 10)))) # 21

G7.eigen.lap <- t(sapply(G7, function(g){
  eigen(graph.laplacian(g))$values
  }))
length(which(duplicated(G7.eigen.lap))) # 2
length(which(duplicated(round(G7.eigen.lap, 10)))) # 14

r <- intersect(

  which(duplicated(round(G7.eigen.adj, 10))),

  which(duplicated(round(G7.eigen.lap, 10)))

  ) # only 441

lapply(r, function(x){

  temp <- G7.eigen.lap - matrix(

    G7.eigen.lap[x, ],

    nrow = nrow(G7.eigen.lap), ncol = 7, byrow = TRUE

    )

  which(rowSums(temp^2) < 0.001)

  })

# 374, 375, 441 

# 374 and 375 have identical degree vec

# 374 not.iso to 375 since 374 contains 3 maximal cliques and 375 contains 2 (max cliques being of size 4)

# rank(374)=5; rank(375)=6 can be seen using pracma::rref and Matrix::rankMatrix (nb, rank(441)=5 also)

# numerical rounding made it look like they had the same spectrum

eig.vs.bet <- t(sapply(G7, function(g){

  lap <- graph.laplacian(g)

  eig <- eigen(lap)

  bet <- betweenness(g)

  which.min.loading <- which.min(abs(eig$vectors[, 6]))

  min.loading <- min(abs(eig$vectors[, 6]))

  which.max.bet <- which.max(bet)

  max.bet <- max(bet)

  cbind(which.min.loading, min.loading, which.max.bet, max.bet)

  }))

# The subset of graphs where lambda2>0 and has multiplicity >= 2

# numerical accuracy < 10^-12

G7.L2.multi <- Filter(

  function(g){

    eig <- eigen(graph.laplacian(g))$values

    abs(eig[5] - eig[6]) < 1E-12

    },

  G7

  )

# 41 graphs

# The subset of graphs where lambda2>0 and has multiplicity 1

G7.L2.single <- Filter(

  function(g){

    eig <- eigen(graph.laplacian(g))$values

    !(abs(eig[5] - eig[6]) < 1E-12)

    },

  G7

  )
# 483 graphs

# Those graphs that have 0-support on at least one vertex for

# every vector in the lambda2-eigenspace

G7.L2.nosupport <- Filter(

  function(g){

    eigvals <- eigen(graph.laplacian(g))$values

    cols <- which(abs(eigvals - eigvals[6]) < 1E-12)

    eigvecs <- eigen(graph.laplacian(g))$vectors

    eig.L2 <- matrix(eigvecs[, cols], ncol = length(cols))

    no.support <- which(rowSums(abs(eig.L2)) < 1E-12)

    length(no.support) >= 1

    },

  G7

  )

# 215 graphs

#######################################################

system.time(G8 <- Filter(is.connected, make_all_graphs(8)))

# Without filtering on eigenvalues
#    user  system elapsed
#   42.17    0.14   42.33

# With filtering on eigenvalues

#    user  system elapsed 

#   19.69    0.02   19.72 

# Recommend working with 10pt accuracy to find 'isospectral' graphs

length(G8) # 4156

G8.eigen.adj <- t(sapply(G8, function(g){
  eigen(get.adjacency(g))$values
  }))

length(which(duplicated(G8.eigen.adj))) # 1
length(which(duplicated(round(G8.eigen.adj, 10)))) # 156
adj.match <- which(duplicated(round(G8.eigen.adj, 10)))

G8.eigen.lap <- t(sapply(G8, function(g){
  eigen(graph.laplacian(g))$values
  }))

length(which(duplicated(G8.eigen.lap))) # 7

length(which(duplicated(round(G8.eigen.lap, 10)))) # 110

lap.match <- which(duplicated(round(G8.eigen.lap, 10)))

r <- intersect(adj.match, lap.match) # 6 entries

# r[1] = 173

# find graphs that have the same Laplacian and Adjacency spectra:

lapply(r, function(x){

  temp <- G8.eigen.lap - matrix(

    G8.eigen.lap[x, ],

    nrow = nrow(G8.eigen.lap), ncol = 8, byrow = TRUE

    )

  which(rowSums(temp^2) < 0.001)

  })

# [1] 103 173

# [1] 228 309

# [1] 1077 1375

# [1] 2103 2110 # have the same degree distribution as well

# [1] 2676 2678 # have the same degree distrubution as well

# [1] 4004 4005 4073 # 4004/5 have the same degree distribution

# The two graphs 2103 and 2110 have the same eigenvalues (to numerical accuracy) for both Laplacian and Adjacency, have the same degree vector, but have differing betweenness vector. They are non-isomorphic, since one contains K_4, but the other doesn't.

# The eigenvectors differ for the two graphs: 

#   loadings follow the degree for the leading adjacency eigenvec;

#   same values in the eigenvec for alg-connectivity but in different places

betweenness(G8[[2103]])

[1] 9.00 3.00 3.00 0.75 0.75 0.50 0.00 0.00

betweenness(G8[[2110]])

[1] 9.00 1.25 1.25 0.00 2.50 2.50 0.50 0.00

degree(G8[[2103]])

[1] 6 4 4 3 3 3 2 1

round(eigen(get.adjacency(G8[[2103]]))$values, 3)

[1]  3.699  1.363  0.618  0.333 -0.706 -1.618 -1.689

[8] -2.000

round(eigen(graph.laplacian(G8[[2103]]))$values, 4)

[1] 7.0908 5.4142 4.8550 3.6023 2.5858 1.5155 0.9365

[8] 0.0000

# lambda_2 eigenvecs (for 2nd smallest eigenval of Laplacian)

round(eigen(graph.laplacian(G8[[2103]]))$vectors[, 7], 3)

[1]  0.057 -0.182 -0.182 -0.090 -0.090 -0.060 -0.343

[8]  0.891

round(eigen(graph.laplacian(G8[[2110]]))$vectors[, 7], 3)

[1] -0.057  0.090  0.090  0.060  0.182  0.182  0.343

[8] -0.891

# mu_1 eigenvecs (for largest eigenval of adjacency mat)

round(eigen(get.adjacency(G8[[2103]]))$vectors[, 1], 3)

[1] -0.526 -0.401 -0.401 -0.339 -0.339 -0.325 -0.217

[8] -0.142

round(eigen(get.adjacency(G8[[2110]]))$vectors[, 1], 3)

[1] -0.523 -0.431 -0.431 -0.369 -0.296 -0.296 -0.158

[8] -0.139

#######################################################################

# Finding non-isomorphic, connected graphs with the same betweenness distribution is easy:

G8.bet <- t(sapply(G8, function(g) sort(betweenness(g))))

head(which(duplicated(G8.bet)))

[1] 323 547 565 579 632 670

G8.bet[323, ]

with(

  list(x = G8.bet - matrix(G8.bet[323, ], nrow = nrow(G8.bet), ncol = 8, byrow = TRUE)),

  which(rowSums(x^2) == 0)

  )

# 115, 323

###########################################################

# TODO
# Find the smallest pair of connected graphs which have the same:

# graph params:
# a) degree distribution # G5[[4]] and G5[[5]]
# b) diameter         # G4[[2]] and G4[[3]]
# c) clique size      # G4[[1]] and G4[[2]]
# d) coclique size    # G4[[1]] and G4[[2]]
# e) betweenness      # G6[[27]] and G6[[34]]

# matrix params:
# f) rank
# g) adjacency eigenvals
# h) laplacian eigenvals

Graphs that are non-isomorphic but have identical betweenness / degree / spectrum etc..

Code to identify a smallest pair of non-isomorphic, connected, undirected graphs that share certain properties:
i)   Share the same degree multiset
ii)  Share the same adjacency spectrum
iii) Share the same laplacian spectrum
iv) Share the same betweenness multiset
v)  Share the same (degree, betweenness) multiset
vi) Share the same adjacency and laplacian spectrum

We consider a graph g1 smaller than g2 if |V(g1)| < |V(g2)| OR |E(g1)| < |E(g2)|

#####################################################
library(igraph)
library(Matrix)

# See other blogposts for make_all_graphs(n) function

####################################################

get_matching_graph_pairs <- function(

    graphset,

    matchset = c('deg', 'bet', 'adjspec', 'lapspec'),

    digits = 10

    ){

  vertex.matchset <- intersect(matchset, c('deg', 'bet'))

  spectrum.matchset  <- intersect(matchset, c('adjspec', 'lapspec'))

  graph.statistics <- NULL

  

  # vertex level matchings
  if(length(vertex.matchset) > 0){

    vertex.statistics <- lapply(graphset, function(g){

      stats.list <- lapply(vertex.matchset, function(vm){

        if(vm == 'deg'){return(data.frame(deg = degree(g)))}

        if(vm == 'bet'){return(data.frame(bet = betweenness(g)))}

        })

      stats.df <- do.call('cbind', stats.list)

      # order the rows by the arbitrary columns:

      stats.df[do.call(order, stats.df), ]

      return(unlist(stats.df))

      })

    vertex.statistics <- do.call('rbind', vertex.statistics)

    graph.statistics <- vertex.statistics
    }

  
  # spectrum matchings

  # Uses Charpoly function as defined in a previous blog
  if(length(spectrum.matchset) > 0){
    spec.statistics <- lapply(graphset, function(g){
      stats.list <- lapply(spectrum.matchset, function(sm){
        if(sm == 'adjspec'){
          return(data.frame(adj = Charpoly(get.adjacency(g))))
          }
        if(sm == 'lapspec'){
          return(data.frame(lap = Charpoly(graph.laplacian(g))))
          }
        })
      stats.df <- do.call('cbind', stats.list)
      return(unlist(stats.df))
      })
    spec.statistics <- do.call('rbind', spec.statistics)
    if(!is.null(graph.statistics)){
      graph.statistics <- cbind(graph.statistics, spec.statistics)
      } else {
      graph.statistics <- spec.statistics
      }
    }

  dups.rows <- which(duplicated(round(

    graph.statistics, digits = digits

    )))

  

  if(length(dups.rows) > 0){

    r <- dups.rows[1]

    rep.matrix <- matrix(

      graph.statistics[r, ],

      nrow = nrow(graph.statistics),

      ncol = ncol(graph.statistics),

      byrow = TRUE

      )

    delta.matrix <- round(graph.statistics - rep.matrix,

                          digits = digits)

    matching.rows <- which(rowSums(abs(delta.matrix)) == 0)

    stopifnot(length(matching.rows) > 1)

    return(graphset[matching.rows[1:2]])

    }

  return(NULL)

  }

####################################################

get_smallest_match <- function(

  matchset = c('deg', 'bet', 'adjspec', 'lapspec'),

  max.vertices = 8,

  digits = 10   

  ){

  # Consider only connected graphs

  # Start from n = 3 vertices, since only 1 graph when n<=2

  for (n in 3:max.vertices){

    # Get a list of all connected, non-isomorphic graphs 

    #   on n vertices

    graphset <- Filter(is_connected, make_all_graphs(n))

    

    # Obtain a list of two graphs that match for all of 

    #   the reqd statistics (if possible) and return them

    #   to the user:

    matching.graphs <- get_matching_graph_pairs(

      graphset,

      matchset,

      digits

      )

    if(!is.null(matching.graphs)){

      return(matching.graphs)

      }

    }

  # No graphs could be found that matched on the required statistics

  return(NULL)

  }

###########################################

With the above code, we can find the smallest pair of graphs that are non-isomorphic, but which have the same 

i) degree distribution:

deg.match <- get_smallest_match('deg')

> lapply(deg.match, degree)
[[1]]
[1] 3 2 2 2 1

[[2]]
[1] 3 2 2 2 1

> lapply(deg.match, betweenness)
[[1]]
[1] 3.5 1.0 1.0 0.5 0.0

[[2]]

[1] 4 0 0 3 0

ii) Betweenness distribution

bet.match <- get_smallest_match('bet')

> lapply(bet.match, degree)
[[1]]
[1] 6 2 2 2 2 2 2

[[2]]
[1] 6 3 3 3 1 1 1

> lapply(bet.match, betweenness)
[[1]]
[1] 12  0  0  0  0  0  0

[[2]]
[1] 12  0  0  0  0  0  0

iii) Distribution of (degree, betweenness) ordered pairs
Note the adjacency eigenvalues differ for these two graphs

db.match  <- get_smallest_match(c('deg', 'bet'))

> lapply(db.match, degree)
[[1]]
[1] 4 4 4 4 4 4 4 4

[[2]]
[1] 4 4 4 4 4 4 4 4

> lapply(db.match, betweenness)
[[1]]
[1] 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5

[[2]]
[1] 1.5 1.5 1.5 1.5 1.5 1.5 1.5 1.5

> lapply(db.match, function(x) round(eigen(get.adjacency(x))$values, 10))
[[1]]
[1]  4  2  0  0  0 -2 -2 -2

[[2]]
[1]  4.000000  1.414214  1.414214  0.000000 -1.414214 -1.414214 -2.000000
[8] -2.000000

iv) Adjacency spectrum

adj.match <- get_smallest_match('adjspec')

lapply(adj.match, plot)

v) Laplacian spectrum

lap.match <- get_smallest_match('lapspec')

lapply(lap.match, plot)

vi) Adjacency spectrum and Degree:

da.match <- get_smallest_match(c('deg', 'adjspec'))

vii) Laplacian spectrum and Degree:

dl.match <- get_smallest_match(c('lapspec', 'deg'))

Note that I could not find a pair of graphs on <= 9 vertices (code is pretty slow for > 8 vertices though) that had 
- both identical Laplacian spectrum and identical adjacency spectrum
- both identical adjacency spectrum and betweenness distribution
- both identical Laplacian spectrum and betweenness distribution 

Characteristic polynomials

In R, we can obtain the characteristic polynomial for a symmetric matrix (such as the adjacency matrices of undirected graphs: see previous posts) using the pracma::charpoly function. This returns the coefficients of the polynomial det(xI - A) in decreasing powers of x.

library(pracma)
library(igraph)
library(Matrix)

# Characteristic polynomial of the 3x3 identity matrix:

diag(c(1,1,1))
#      [,1] [,2] [,3]
# [1,]    1    0    0
# [2,]    0    1    0
# [3,]    0    0    1

charpoly(diag(c(1,1,1)))

[1]  1 -3  3 -1

# x^3 -3x^2 +3x -1

# adjacency matrix for the Petersen graph:
pet <- get.adjacency(make_graph('petersen'))

pet
#10 x 10 sparse Matrix of class "dgCMatrix"
#                         
# [1,] . 1 . . 1 1 . . . .
# [2,] 1 . 1 . . . 1 . . .
# [3,] . 1 . 1 . . . 1 . .
# [4,] . . 1 . 1 . . . 1 .
# [5,] 1 . . 1 . . . . . 1
# [6,] 1 . . . . . . 1 1 .
# [7,] . 1 . . . . . . 1 1
# [8,] . . 1 . . 1 . . . 1
# [9,] . . . 1 . 1 1 . . .
#[10,] . . . . 1 . 1 1 . .

charpoly(pet)

# Error: is.numeric(a) is not TRUE

igraph adjacency matrices are of 'Matrix' class and can't be used directly in charpoly (which expects base R numeric matrices as input). We define a function to return the coefficients of the characteristic polynomial for a graph, or a 'Matrix' or a matrix:

Charpoly <- function(graph.or.matrix){

  require(pracma)

  require(igraph)

  require(Matrix)

  if(is(graph.or.matrix, 'igraph')){

    return(Charpoly(get.adjacency(graph.or.matrix)))

    }

  if(is(graph.or.matrix, 'Matrix')){

    return(Charpoly(as.matrix(graph.or.matrix)))

    }

  charpoly(graph.or.matrix)

  }

Charpoly(make_graph('petersen'))

 [1]    1    0  -15    0   75  -24 -165  120  120 -160   48

That is, the characteristic polynomial is 

x^10 - 15x^8 + 75x^6 -24x^5 ....

Charpoly(pet)

 [1]    1    0  -15    0   75  -24 -165  120  120 -160   48