For prime graph decomposition, one repeatedly looks for cliques that separate the graph.
Please refer to Appendix B.2.2 (page 124) of https://web.math.ku.dk/~lauritzen/papers/gmnotes.pdf for a formal definition of prime graph decomposition.
I am wondering if it has already been implemented in igraph? If not, can I contribute to it? I have already written a function in R that takes an igraph object and output a list of prime components and a list of separators.
Thanks!
Contributions are welcome!
I’m sympathetic towards the idea of including this in igraph, but keep in mind that the igraph team is quite busy, so it’ll take some time to get there.
If we incorporate it into igraph, it should be implemented in C. An R prototype is a good first step. Would you mind sharing your R implementation? Are you comfortable programming in C?
Let’s use this thread to discuss the details.
Thank you! At present I have the following code. I just run through all cliques until I found one that separates the graph. I do not have much experience in C but I don’t mind it.
prime_decom <- function(graph) {
prime_components <- list()
min_sep <- list()
subgraph_stack <- groups(components(graph))
while (length(subgraph_stack) >= 1) {
found = FALSE # a flag
# pop the last graph from stack
curr_graph_v <- subgraph_stack[[length(subgraph_stack)]] # get
subgraph_stack <- subgraph_stack[-length(subgraph_stack)] # remove
curr_graph_v <- sort(curr_graph_v)
curr_graph <- induced_subgraph(graph, curr_graph_v, impl = c("auto"))
all_cliques <- cliques(curr_graph)
for (cli in all_cliques) {
# obtain the induced subgraph by removing cli
remaining_v <- setdiff(V(curr_graph), cli)
removed <- induced_subgraph(curr_graph, remaining_v, impl = c("auto"))
# get connected components of the induced subgraph
clu <- components(removed)
if (clu$no > 1) {
found = TRUE
connected_components <- groups(clu)
# append subgraph_stack
length(subgraph_stack) <- length(subgraph_stack) + clu$no
for (i in 1:clu$no) {
vertices <- union(remaining_v[connected_components[[i]]], cli)
subgraph_stack[[length(subgraph_stack) - i + 1]] <- curr_graph_v[vertices]
}
# append min_sep
length(min_sep) <- length(min_sep) + (clu$no - 1)
for (i in 1:(clu$no - 1)) {
min_sep[[length(min_sep) - i + 1]] <- curr_graph_v[cli]
}
# break the for loop
break
}
}
if (!found) {
# it is a prime component
prime_components <- c(prime_components, list(curr_graph_v))
}
}
return (list(components = prime_components, separators = min_sep))
}
Are you aware of any work on efficient algorithms for decomposition?
Would it work (better) to look for minimal separators and check if they’re cliques?
Have a look at:
https://igraph.org/c/html/latest/igraph-Separators.html#igraph_all_minimal_st_separators
https://igraph.org/c/html/latest/igraph-Separators.html#igraph_minimum_size_separators
Looking for minimal separators is indeed faster. Here I have updated the code.
Since I cannot find the function is_clique in R, I created a subgraph and check the number of edges. I see that this function is available in C.
library(igraph)
prime_decom2 <- function(graph) {
prime_components <- list()
min_sep <- list()
subgraph_stack <- groups(components(graph))
while (length(subgraph_stack) >= 1) {
found = FALSE # a flag
# pop the last graph from stack
curr_graph_v <- subgraph_stack[[length(subgraph_stack)]] # get
subgraph_stack <- subgraph_stack[-length(subgraph_stack)] # remove
curr_graph <- induced_subgraph(graph, curr_graph_v, impl = c("auto"))
## if curr_graph is complete, then prime component
num_edges <- length(curr_graph_v) * (length(curr_graph_v) - 1) / 2
if (ecount(curr_graph) != num_edges) {
all_separators <- min_separators(curr_graph)
for (sep in all_separators) {
induced_sep <- induced_subgraph(curr_graph, sep, impl = c("auto"))
num_edges <- length(sep) * (length(sep) - 1) / 2
if (ecount(induced_sep) == num_edges) {
found = TRUE
curr_graph_v <- sort(curr_graph_v)
# obtain the induced subgraph by removing sep
remaining_v <- setdiff(V(curr_graph), sep)
removed <- induced_subgraph(curr_graph, remaining_v, impl = c("auto"))
# get connected components of the induced subgraph
clu <- components(removed)
connected_components <- groups(clu)
# append subgraph_stack
length(subgraph_stack) <- length(subgraph_stack) + clu$no
for (i in 1:clu$no) {
vertices <- union(remaining_v[connected_components[[i]]], sep)
subgraph_stack[[length(subgraph_stack) - i + 1]] <- curr_graph_v[vertices]
}
# append min_sep
length(min_sep) <- length(min_sep) + (clu$no - 1)
for (i in 1:(clu$no - 1)) {
min_sep[[length(min_sep) - i + 1]] <- curr_graph_v[sep]
}
# break the for loop
break
}
}
}
if (!found) {
# it is a prime component
prime_components <- c(prime_components, list(curr_graph_v))
}
}
return (list(components = prime_components, separators = min_sep))
}
Looks like there is an O(mn) algorithm for this decomposition. I will update my code here once I have it written.
Great find! I would suggest doing a bit of additional literature search before embarking on the implementation. It is useful to find out:
Instead of going for the best known complexity algorithm, at first, I would go for algorithms which seem easier to implement and have a clear description available.
A good next step in literature search is to see what papers cite this one. I haven’t looked at the contents, but the title of this one is promising:
Thank you szhorvat.
I implemented an algorithm presented in a recent paper https://epubs.siam.org/doi/10.1137/16M1059837.
The algorithm goes like this:
(1) Given a graph G, get a minimal chordal completion H.
(2) Get a clique tree of H (so we need the cliques and separators of H, see chordal_factor() below).
[Fact: the minimal separators of H which forms a clique in G are exactly the minimal separators of G.]
(3) Check the minimal separators of H one by one. If the separator does not form a clique in G, then contract that edge in the clique tree.
(4) The resulting clique tree becomes the atom tree, which gives the prime components.
In (1), I rely on the assumption that the output newgraph of the igraph function is_chordal() gives a minimal chordal completion. Could you please verify this?
The algorithm theoretically takes O(n^alpha log n) time, where alpha < 2.3729, using fast matrix multiplication to compute the matrix C (middle of prime_decom3()). Here I simply use %*% in R.
###### input graph must be chordal
chordal_factor <- function(graph) {
n <- vcount(graph)
if (n == 0) {
return(list("cliques" = list(), "separators" = list()))
}
PEO <- as.vector(max_cardinality(graph)$alpham1)
cliques <- list()
kdeg <- vector()
representatives <- vector()
perfectCliques <- list()
separators <- list()
# all possible cliques
for (i in 1:n) {
v <- PEO[i]
neighbours <- intersect(neighbors(graph, v), PEO[i:n])
cliques[[i]] <- c(v, neighbours)
kdeg[i] <- length(cliques[[i]]) - 1
}
# perfect sequence of cliques
perfectCliques <- list(cliques[[1]])
representatives[1] <- PEO[1]
if (n > 1) {
for (i in 2:n) {
if (kdeg[[i]] >= kdeg[[i-1]]) {
representatives <- c(representatives, PEO[i])
perfectCliques <- append(perfectCliques, list(cliques[[i]]))
}
}
perfectCliques <- rev(perfectCliques)
}
# separators
tempUnion <- perfectCliques[[1]]
num_maximal_cliques <- length(perfectCliques)
if (num_maximal_cliques > 1) {
for (i in 2:num_maximal_cliques) {
separators <- append(separators, list(intersect(tempUnion, perfectCliques[[i]])))
tempUnion <- union(tempUnion, perfectCliques[[i]])
}
}
return(list("cliques" = perfectCliques, "separators" = separators))
}
###### input graph is an undirected graph
# assume is_chordal() gives minimal chordal completion
prime_decom3 <- function(graph) {
n <- vcount(graph)
# get a minimal chordal completion of graph
chordal <- is_chordal(graph, newgraph = TRUE)
if (chordal$chordal) {
return (chordal_factor(graph))
} else {
chordal <- chordal$newgraph
factor <- chordal_factor(chordal)
cliques <- factor$cliques
separators <- factor$separators
}
min_sep <- list()
# define clique_matrix and C
clique_matrix <- matrix(0, nrow = n, ncol = length(separators))
for (i in seq_along(separators)) {
clique_matrix[separators[[i]], i] <- 1
}
C <- as_adjacency_matrix(graph, "both") %*% clique_matrix # use fast matrix multiplication
# check
remove <- vector()
for (i in seq_along(separators)) {
sep <- separators[[i]]
if (all(C[as.vector(sep), i] == length(sep) - 1)) {
min_sep <- append(min_sep, list(sep))
} else {
remove <- c(remove, i)
for (j in 1:i) {
if (all(sep %in% cliques[[j]])) {
cliques[[j]] <- union(cliques[[i]], cliques[[j]])
cliques[[i]] <- vector()
break
}
}
}
}
return (list(components = cliques[-remove], separators = min_sep))
}
If is_chordal() does not give minimal chordal completion, I can also implement the algorithm described in Section 4 of the review paper you posted. It’s self-contained.
It turns out that is_chordal() does not give a minimal chordal completion. I found a counter example:
set.seed(124)
graph <- sample_grg(20, 0.4)
chordal <- is_chordal(graph, newgraph = TRUE)$newgraph
print(ecount(chordal)) # 59
print(ecount(graph)) # 49
print(is_chordal(delete_edges(chordal, 54))) # true
I think it is a good idea to modify the implementation of is_chordal() as well (or add a new function for minimal chordal completion if the current chordal completion has other nice properties).
I’m sorry, I forgot to respond here. I suspected that it doesn’t give a minimal one, I didn’t have time to find a reference or example yesterday.
The current algorithm is quite efficient. I’m in favour of writing a new function for a minimal chordal completion.
Took a bit of time to get used to C and the igraph library. Here I have implemented the 2 algorithms described in An Introduction to Clique Minimal Separator Decomposition. I mostly followed the pseudocode in the paper.
This draft of code seems to be working. I have not, for example, handled the errors, destroyed variables.
#include <igraph.h>
void MCS_M_plus(const igraph_t *g,
igraph_vector_int_t *alpha,
igraph_vector_int_t *min_sep_gen,
igraph_t *minimal_chordal);
void Atoms(const igraph_t *g,
const igraph_vector_int_t *alpha,
const igraph_vector_int_t *min_sep_gen,
const igraph_t *minimal_chordal,
igraph_vector_int_list_t *A,
igraph_vector_int_list_t *S_H,
igraph_vector_int_list_t *S_C);
int main() {
int n = 20;
int failed_cases = 0;
float radius = 0.5;
int num_cases = 10;
for (int i = 0; i < num_cases; i++) {
printf("======================================\n");
igraph_t g;
srand(i);
igraph_grg_game(&g, n, radius, 0, NULL, NULL);
igraph_t minimal_chordal;
igraph_vector_int_t alpha, min_sep_gen;
igraph_vector_int_init(&alpha, 0);
igraph_vector_int_init(&min_sep_gen, 0);
MCS_M_plus(&g, &alpha, &min_sep_gen, &minimal_chordal);
printf("alpha ");
igraph_vector_int_print(&alpha);
// igraph_is_chordal() on g
igraph_bool_t chordal;
igraph_vector_int_t fill_in;
igraph_vector_int_init(&fill_in, 0);
igraph_t newgraph;
igraph_is_chordal(&g, NULL, NULL, &chordal, &fill_in, &newgraph);
// compare number of edges
printf("g: %lld, minimal_chordal: %lld, newgraph: %lld, difference: %lld \n", igraph_ecount(&g), igraph_ecount(&minimal_chordal), igraph_ecount(&newgraph), igraph_ecount(&minimal_chordal) - igraph_ecount(&newgraph));
// check minimal_chordal is chordal
igraph_bool_t chordal_m;
igraph_vector_int_t fill_in_m;
igraph_vector_int_init(&fill_in_m, 0);
igraph_t newgraph_m;
igraph_is_chordal(&minimal_chordal, NULL, NULL, &chordal_m, &fill_in_m, &newgraph_m);
if (!chordal_m) {
printf("minimal_chordal is not chordal!");
failed_cases ++;
}
// check minimal_chordal is minimal
igraph_t g_small;
igraph_t newgraph_small;
igraph_es_t es;
igraph_vector_int_t fill_in_small;
igraph_bool_t chordal_small;
for (int j = 0; j < igraph_ecount(&minimal_chordal) - igraph_ecount(&g); j++) {
igraph_copy(&g_small, &minimal_chordal);
igraph_es_1(&es, igraph_ecount(&g)+j);
igraph_delete_edges(&g_small, es);
igraph_vector_int_init(&fill_in_small, n);
igraph_is_chordal(&g_small, NULL, NULL, &chordal_small, &fill_in_small, &newgraph_small);
if (chordal_small == 1) {
printf("minimal_chordal is not minimal! \n");
printf("i is %d \n", i);
printf("j is %d \n", j);
failed_cases ++;
}
}
igraph_vector_int_list_t A;
igraph_vector_int_list_t S_H;
igraph_vector_int_list_t S_C;
igraph_vector_int_list_init(&A, 0);
igraph_vector_int_list_init(&S_H, 0);
igraph_vector_int_list_init(&S_C, 0);
Atoms(&g, &alpha, &min_sep_gen, &minimal_chordal, &A, &S_H, &S_C);
printf("The set of atoms: \n");
igraph_integer_t size_A = igraph_vector_int_list_size(&A);
for (int j = 0; j < size_A; j++) {
igraph_vector_int_print(igraph_vector_int_list_get_ptr(&A, j));
}
printf("The separators: \n");
igraph_integer_t size_S_C = igraph_vector_int_list_size(&S_C);
for (int j = 0; j < size_S_C; j++) {
igraph_vector_int_print(igraph_vector_int_list_get_ptr(&S_C, j));
}
}
printf("======================================\n");
printf("Number of failed cases is %d. \n", failed_cases);
return 0;
}
void MCS_M_plus(const igraph_t *g,
igraph_vector_int_t *alpha,
igraph_vector_int_t *min_sep_gen,
igraph_t *minimal_chordal) {
// the number of vertices in the graph
const igraph_integer_t n = igraph_vcount(g);
// init
igraph_copy(minimal_chordal, g);
igraph_vector_int_init(alpha, n);
igraph_vector_int_init(min_sep_gen, 0);
igraph_vector_int_t F; // to store the fill_in edges
igraph_vector_int_init(&F, 0);
igraph_integer_t s = -1;
// initialise the labels of all vertices as 0
// label = -1 if the vertex has been chosen
igraph_vector_int_t label;
igraph_vector_int_init(&label, n);
// for loop
igraph_integer_t x, y, z, min_label, size_Y, size_Z, index;
igraph_vector_int_t Y, Z, reached;
igraph_vector_int_list_t reach;
igraph_bool_t x_y_adj;
for (int i = n-1; i >= 0; i--) {
// choose a vertex x of g of maximal label
igraph_vector_int_which_minmax(&label, &min_label, &x);
// Y <- N_{g}(x)
igraph_vector_int_init(&Y, 0);
igraph_neighbors(g, &Y, x, IGRAPH_ALL);
// remove from Y the vertices with label -1
size_Y = igraph_vector_int_size(&Y);
index = 0;
for (int j = 0; j < size_Y; j++) {
if (VECTOR(label)[VECTOR(Y)[j]] > -1) {
VECTOR(Y)[index] = VECTOR(Y)[j];
index ++;
}
}
igraph_vector_int_resize(&Y, index);
// if loop
if (VECTOR(label)[x] <= s) {
igraph_vector_int_push_back(min_sep_gen, x);
}
// s <- label(x)
s = VECTOR(label)[x];
// mark x reached and mark all other vertices of g unreached
igraph_vector_int_init(&reached, n);
VECTOR(reached)[x] = 1;
// for loop, initialise reach as a list of vectors
igraph_vector_int_list_init(&reach, n);
for (int j = 0; j < n; j++) {
igraph_vector_int_init(igraph_vector_int_list_get_ptr(&reach, j), 0);
}
// for loop, neighbours of x
size_Y = igraph_vector_int_size(&Y);
for (int j = 0; j < size_Y; j++) {
y = VECTOR(Y)[j];
// make y reached
VECTOR(reached)[y] = 1;
// add y to reach(label(y))
igraph_vector_int_push_back(igraph_vector_int_list_get_ptr(&reach, VECTOR(label)[y]), y);
}
// for loop, other vertices
for (int j = 0; j < n; j++) {
// while reach(j) \neq empty
while (igraph_vector_int_size(igraph_vector_int_list_get_ptr(&reach, j)) > 0) {
// remove a vertex y from reach(j)
y = igraph_vector_int_pop_back(igraph_vector_int_list_get_ptr(&reach, j));
// for loop, neighbours of y, Z is the neighbours of y in g
igraph_vector_int_init(&Z, 0);
igraph_neighbors(g, &Z, y, IGRAPH_ALL);
size_Z = igraph_vector_int_size(&Z);
for (int k = 0; k < size_Z; k++) {
z = VECTOR(Z)[k];
// if z is unreached
if (VECTOR(label)[z] > -1 && VECTOR(reached)[z] == 0) {
// mark z reached
VECTOR(reached)[z] = 1;
// if label(z) > j
if (VECTOR(label)[z] > j) {
// Y <- Y + {z}
igraph_vector_int_push_back(&Y, z);
// add z to reach(label(z))
igraph_vector_int_push_back(igraph_vector_int_list_get_ptr(&reach, VECTOR(label)[z]), z);
} else {
// add z to reach(j)
igraph_vector_int_push_back(igraph_vector_int_list_get_ptr(&reach, j), z);
}
}
}
}
}
// for loop
size_Y = igraph_vector_int_size(&Y);
for (int j = 0; j < size_Y; j++) {
y = VECTOR(Y)[j];
// F <- F + {xy} if {x,y} is not an edge in g, here one can also add everything then use igraph_simplify() to remove multiple edges
igraph_are_adjacent(g, x, y, &x_y_adj);
if (!x_y_adj) {
igraph_vector_int_push_back(&F, x);
igraph_vector_int_push_back(&F, y);
}
// label(y) <- label(y) + 1
VECTOR(label)[y] ++;
}
// alpha(i) <- x
VECTOR(*alpha)[i] = x;
// remove x from g, i.e. set the label -1
VECTOR(label)[x] = -1;
}
// add the edges of F to minimal_chordal
igraph_add_edges(minimal_chordal, &F, NULL);
}
void Atoms(const igraph_t *g,
const igraph_vector_int_t *alpha,
const igraph_vector_int_t *min_sep_gen,
const igraph_t *minimal_chordal,
igraph_vector_int_list_t *A,
igraph_vector_int_list_t *S_H,
igraph_vector_int_list_t *S_C) {
igraph_vector_int_list_init(A, 0);
igraph_vector_int_list_init(S_H, 0);
igraph_vector_int_list_init(S_C, 0);
// the number of vertices in the graph
const igraph_integer_t n = igraph_vcount(g);
// label = 1 if the vertex has been chosen, 0 otherwise (for minimal_chordal)
igraph_vector_int_t label;
igraph_vector_int_init(&label, n);
// g1 is G'
igraph_t g1;
igraph_copy(&g1, g);
// g1_v is the vertices of g1
igraph_vector_int_t g1_v;
igraph_vector_int_init_range(&g1_v, 0, n);
igraph_vector_int_t map, invmap;
igraph_vector_int_init(&map, 0);
igraph_vector_int_init(&invmap, 0);
// for loop
igraph_integer_t x, size_S, size_C, index;
igraph_vector_int_t S,g1_minus_S_v, C, g1_minus_C_v;
igraph_bool_t S_is_clique;
igraph_vs_t S_vids, g1_minus_S_vids, g1_minus_C_vids;
igraph_t g1_minus_S, g2;
for (int i = 0; i < n; i++) {
// x <- alpha(i)
x = VECTOR(*alpha)[i];
// if x in X
if (igraph_vector_int_contains(min_sep_gen, x)) {
// S <- N_{H'}(x)
//// get the neighbours of x in H
igraph_vector_int_init(&S, 0);
igraph_neighbors(minimal_chordal, &S, x, IGRAPH_ALL);
//// remove those neighbours not in H', i.e., with label != 0
size_S = igraph_vector_int_size(&S);
index = 0;
for (int j = 0; j < size_S; j++) {
if (VECTOR(label)[VECTOR(S)[j]] == 0) {
VECTOR(S)[index] = VECTOR(S)[j];
index ++;
}
}
igraph_vector_int_resize(&S, index);
//// make a copy of S
igraph_vector_int_t S_copy;
igraph_vector_int_init_copy(&S_copy, &S);
// S_H append S
if (igraph_vector_int_size(&S_copy) > 0) {
igraph_vector_int_list_push_back(S_H, &S_copy);
}
// if S is a clique in G
igraph_vs_vector(&S_vids, &S_copy);
igraph_is_clique(g, S_vids, 0, &S_is_clique);
if (S_is_clique == 1) {
// S_C append S
if (igraph_vector_int_size(&S_copy) > 0) {
igraph_vector_int_list_push_back(S_C, &S_copy);
}
// C <- the connected component of G' - S containing x
//// g1_minus_S is the graph G'-S
igraph_vector_int_init(&g1_minus_S_v, 0);
igraph_vector_int_difference_sorted(&g1_v, &S, &g1_minus_S_v);
igraph_vs_vector(&g1_minus_S_vids, &g1_minus_S_v);
igraph_induced_subgraph_map(g, &g1_minus_S, g1_minus_S_vids, IGRAPH_SUBGRAPH_AUTO, &map, &invmap);
//// C is the component containing x
igraph_vector_int_init(&C, 0);
igraph_subcomponent(&g1_minus_S, &C, VECTOR(map)[x] - 1, IGRAPH_ALL);
size_C = igraph_vector_int_size(&C);
for (int k = 0; k < size_C; k++) {
VECTOR(C)[k] = VECTOR(invmap)[VECTOR(C)[k]];
}
igraph_vector_int_sort(&C);
// A append S union C
igraph_vector_int_append(&S, &C);
igraph_vector_int_list_push_back(A, &S);
// G' <- G' - C
igraph_vector_int_init (&g1_minus_C_v, 0);
igraph_vector_int_difference_sorted(&g1_v, &C, &g1_minus_C_v);
igraph_vs_vector(&g1_minus_C_vids, &g1_minus_C_v);
igraph_induced_subgraph_map(g, &g2, g1_minus_C_vids, IGRAPH_SUBGRAPH_AUTO, &map, &invmap);
g1_v = g1_minus_C_v;
g1 = g2;
}
}
//remove x from H', i.e., set the label to 1
VECTOR(label)[x] = 1;
}
// A append g1_v
igraph_vector_int_list_push_back(A, &g1_v);
}
The output looks good to me. But I am not sure how to compare the results with the code I had in R. Do you use Rcpp or .C()?
======================================
alpha 18 17 3 19 15 16 14 13 12 11 10 9 7 6 8 5 4 2 1 0
g: 83, minimal_chordal: 95, newgraph: 90, difference: 5
The set of atoms:
10 13 17 18
0 2 4 5 3
9 11 12 13 14 15 16 19
1 2 4 5 8 6 7 9 10 11 12 13 14 15 16
0 1 2 4 5 8
The separators:
10 13
0 2 4 5
9 11 12 13 14 15 16
1 2 4 5 8
======================================
alpha 14 19 16 15 18 17 13 12 11 2 10 9 6 8 7 5 4 3 1 0
g: 86, minimal_chordal: 91, newgraph: 93, difference: -2
The set of atoms:
10 19 14
10 12 13 17 18 19
8 9 11 12 13 17 18 15 16
1 4 5 6 2
0 1 3 4 5 6 7 8 9 10 11 12 13 17 18
The separators:
10 19
10 12 13 17 18
8 9 11 12 13 17 18
1 4 5 6
======================================
alpha 15 18 13 19 17 16 14 10 8 7 6 9 12 11 5 4 3 2 1 0
g: 86, minimal_chordal: 116, newgraph: 108, difference: 8
The set of atoms:
11 12 14 16 17 19 15
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 16 17 18 19
The separators:
11 12 14 16 17 19
======================================
alpha 19 17 3 16 18 11 10 8 15 14 13 12 9 6 7 5 4 2 1 0
g: 105, minimal_chordal: 111, newgraph: 108, difference: 3
The set of atoms:
8 9 10 11 12 16 18 17 19
5 6 7 8 9 10 11 12 3
5 6 7 9 12 13 8 10 11 14 15 16 18
4 5 6 7 9 12 13
0 4 5 7 6
0 1 4 5 7
0 1 2 4 5
The separators:
8 9 10 11 12 16 18
5 6 7 8 9 10 11 12
5 6 7 9 12 13
4 5 6 7
0 4 5 7
0 1 4 5
======================================
alpha 15 17 19 18 14 16 13 12 11 9 10 8 7 2 6 5 4 3 1 0
g: 80, minimal_chordal: 86, newgraph: 84, difference: 2
The set of atoms:
8 11 13 15
9 12 14 16 18 19 17
9 10 12 16 14 18 19
8 9 10 12 16
6 8 10 11 13
7 8 9 10 12
1 3 4 5 6 2 7 8 9 10 11
0 1 3 4 5 6
The separators:
8 11 13
9 12 14 16 18 19
9 10 12 16
8 9 10 12
6 8 10 11
7 8 9 10
1 3 4 5 6
======================================
alpha 19 18 17 16 13 12 8 7 6 14 15 11 5 2 10 9 4 3 1 0
g: 107, minimal_chordal: 114, newgraph: 116, difference: -2
The set of atoms:
8 12 13 14 15 16 17 18
2 5 6 11 14 15 7 8 12 13 16 17
2 4 5 11 14 15 6
0 3 4 9 10 2 5 11 14 15 19
0 1 3 4 9 10
The separators:
8 12 13 14 15 16 17
2 5 6 11 14 15
2 4 5 11 14 15
0 3 4 9 10
======================================
alpha 19 14 18 17 16 15 12 11 13 6 5 4 2 10 9 8 7 3 1 0
g: 86, minimal_chordal: 96, newgraph: 93, difference: 3
The set of atoms:
12 14 19
11 12 14
13 15 16 17 18
10 13 15 16 17
7 9 10 13 15
0 3 7 8 9 10 2 4 5 6 11 12 13
0 1 3 7 8 9 10
The separators:
12 14
11 12
13 15 16 17
10 13 15
7 9 10 13
0 3 7 8 9 10
======================================
alpha 17 14 19 18 1 16 12 7 5 15 13 11 10 6 9 3 8 4 2 0
g: 98, minimal_chordal: 108, newgraph: 109, difference: -1
The set of atoms:
8 10 13 14 15 17
8 10 11 13 15 14
3 4 5 6 7 9 1
0 4 8 3 5 6 7 9 10 11 12 13 15 16 18 19
0 2 4 8
The separators:
8 10 13 14 15
8 10 11 13 15
3 4 5 6 7 9
0 4 8
======================================
alpha 16 15 10 7 19 18 17 14 1 13 12 11 9 6 8 5 4 3 2 0
g: 87, minimal_chordal: 96, newgraph: 97, difference: -1
The set of atoms:
7 10 14 15 16
3 4 5 6 8 9 11 1 7 10 12 13 14 17 18 19
2 3 4 5 8 6 9 11
0 2 3 4 5 8
The separators:
7 10 14
3 4 5 6 8 9 11
2 3 4 5 8
======================================
alpha 1 3 18 17 15 19 16 14 11 9 13 12 10 6 5 2 8 7 4 0
g: 89, minimal_chordal: 103, newgraph: 106, difference: -3
The set of atoms:
2 3 5 6 9 1
2 5 6 7 9 3
11 14 15 16 17 19 18
4 7 8 2 5 6 9 10 11 12 13 14 15 16 17 19
0 4 7 8
The separators:
2 3 5 6 9
2 5 6 7 9
11 14 15 16 17 19
4 7 8
======================================
Number of failed cases is 0.
Notice that in some cases igraph_is_chordal() gives a chordal completion with fewer edges. So I agree that we should keep is_chordal() as is. The function MCS_M_plus() always gives a minimal chordal completion, which could also be useful.
Could you please open a pull request at GitHub - igraph/igraph: Library for the analysis of networks ?
That will make it easier me me to have a look at it.
You don’t need to worry about code quality or proper integration into igraph for now. All I’m looking for is an easier way to test and review, which a barebones PR can do.
I am extremely overwhelmed at the moment so it will take some timebefore I can take a look, probably on the week of December 16th.
Thanks szhorvat. I have opened a pull request.
Please ignore the above R code. Here is an updated one written in R, by brute force. Both atom functions should return the same atoms and separators, unless I made some mistakes.
library(igraph)
atoms <- function(graph) {
prime_components <- list()
min_sep <- list()
subgraph_stack <- groups(components(graph))
while (length(subgraph_stack) >= 1) {
found = FALSE # a flag
# pop the last graph from stack
curr_graph_v <- subgraph_stack[[length(subgraph_stack)]] # get
subgraph_stack <- subgraph_stack[-length(subgraph_stack)] # remove
curr_graph_v <- sort(curr_graph_v)
curr_graph <- induced_subgraph(graph, curr_graph_v, impl = c("auto"))
all_cliques <- cliques(curr_graph)
# sort the cliques by length to make sure we look at minimal clique
all_cliques <- all_cliques[order(sapply(all_cliques, length))]
for (cli in all_cliques) {
# obtain the induced subgraph by removing cli
remaining_v <- setdiff(V(curr_graph), cli)
removed <- induced_subgraph(curr_graph, remaining_v, impl = c("auto"))
# get connected components of the induced subgraph
clu <- components(removed)
if (clu$no > 1) {
found = TRUE
connected_components <- groups(clu)
# append subgraph_stack
length(subgraph_stack) <- length(subgraph_stack) + clu$no
for (i in 1:clu$no) {
vertices <- union(remaining_v[connected_components[[i]]], cli)
subgraph_stack[[length(subgraph_stack) - i + 1]] <- curr_graph_v[vertices]
}
# append min_sep
length(min_sep) <- length(min_sep) + (clu$no - 1)
for (i in 1:(clu$no - 1)) {
min_sep[[length(min_sep) - i + 1]] <- curr_graph_v[cli]
}
# break the for loop
break
}
}
if (!found) {
# it is a prime component
prime_components <- c(prime_components, list(curr_graph_v))
}
}
return (list(components = prime_components, separators = min_sep))
}