# Is Prime Graph Decomposition implemented?

**URL:** <https://igraph.discourse.group/t/is-prime-graph-decomposition-implemented/1895>\
**Category:** Development\
**Tags:** R\
**Created:** [15 November 2024 11:39 UTC](https://igraph.discourse.group/t/is-prime-graph-decomposition-implemented/1895 "2024-11-15T11:39:30Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![ching-wong](https://yyz2.discourse-cdn.com/free1/user_avatar/igraph.discourse.group/ching-wong/32/1130_2.png) [@ching-wong](https://igraph.discourse.group/u/ching-wong)\
**Post date:** [15 November 2024 11:39 UTC](https://igraph.discourse.group/t/is-prime-graph-decomposition-implemented/1895/1 "2024-11-15T11:39:30Z")

</div>

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](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!

---

<div class="post-metadata">

**Author:** ![szhorvat](https://yyz2.discourse-cdn.com/free1/user_avatar/igraph.discourse.group/szhorvat/32/3_2.png) [@szhorvat](https://igraph.discourse.group/u/szhorvat)\
**Post date:** [15 November 2024 13:41 UTC](https://igraph.discourse.group/t/is-prime-graph-decomposition-implemented/1895/2 "2024-11-15T13:41:29Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ching-wong](https://yyz2.discourse-cdn.com/free1/user_avatar/igraph.discourse.group/ching-wong/32/1130_2.png) [@ching-wong](https://igraph.discourse.group/u/ching-wong)\
**Post date:** [15 November 2024 14:57 UTC](https://igraph.discourse.group/t/is-prime-graph-decomposition-implemented/1895/3 "2024-11-15T14:57:07Z")

</div>

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.

```auto
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))
}

```

---

<div class="post-metadata">

**Author:** ![szhorvat](https://yyz2.discourse-cdn.com/free1/user_avatar/igraph.discourse.group/szhorvat/32/3_2.png) [@szhorvat](https://igraph.discourse.group/u/szhorvat)\
**Post date:** [15 November 2024 19:13 UTC](https://igraph.discourse.group/t/is-prime-graph-decomposition-implemented/1895/4 "2024-11-15T19:13:05Z")

</div>

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:

> **[Minimum size vertex separators — min\_separators](https://r.igraph.org/reference/min_separators.html)**
>
> Find all vertex sets of minimal size whose removal separates the graph into
> more components

> **[Minimum size vertex separators — min\_st\_separators](https://r.igraph.org/reference/min_st_separators.html)**
>
> List all vertex sets that are minimal \\((s,t)\\) separators for some
> \\(s\\) and \\(t\\), in an undirected graph.

[https://igraph.org/c/html/latest/igraph-Separators.html#igraph\_all\_minimal\_st\_separators](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](https://igraph.org/c/html/latest/igraph-Separators.html#igraph_minimum_size_separators)

---

<div class="post-metadata">

**Author:** ![ching-wong](https://yyz2.discourse-cdn.com/free1/user_avatar/igraph.discourse.group/ching-wong/32/1130_2.png) [@ching-wong](https://igraph.discourse.group/u/ching-wong)\
**Post date:** [18 November 2024 08:28 UTC](https://igraph.discourse.group/t/is-prime-graph-decomposition-implemented/1895/5 "2024-11-18T08:28:48Z")

</div>

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.

```auto
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))
}

```

---

<div class="post-metadata">

**Author:** ![ching-wong](https://yyz2.discourse-cdn.com/free1/user_avatar/igraph.discourse.group/ching-wong/32/1130_2.png) [@ching-wong](https://igraph.discourse.group/u/ching-wong)\
**Post date:** [18 November 2024 09:28 UTC](https://igraph.discourse.group/t/is-prime-graph-decomposition-implemented/1895/6 "2024-11-18T09:28:48Z")

</div>

Looks like there is an O(mn) algorithm for this decomposition. I will update my code here once I have it written.

> **[Optimal decomposition by clique separators](https://www.sciencedirect.com/science/article/pii/0012365X9390510Z?ref=pdf_download&fr=RR-2&rr=8e46ec82eeac3b57)**
>
> Decompositions of a graph by clique separators are investigated which have the additional property that they do not generate new maximal prime subgrap…

---

<div class="post-metadata">

**Author:** ![szhorvat](https://yyz2.discourse-cdn.com/free1/user_avatar/igraph.discourse.group/szhorvat/32/3_2.png) [@szhorvat](https://igraph.discourse.group/u/szhorvat)\
**Post date:** [18 November 2024 09:36 UTC](https://igraph.discourse.group/t/is-prime-graph-decomposition-implemented/1895/7 "2024-11-18T09:36:11Z")

</div>

Great find! I would suggest doing a bit of additional literature search before embarking on the implementation. It is useful to find out:

- Are there more modern algorithms?
- Are there easier to follow explanations of the algorithm, or existing implementations that we can use for testing?

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:

> **[An Introduction to Clique Minimal Separator Decomposition](https://www.mdpi.com/1999-4893/3/2/197)**
>
> This paper is a review which presents and explains the decomposition of graphs by clique minimal separators. The pace is leisurely, we give many examples and figures. Easy algorithms are provided to implement this decomposition. The historical and...

---

<div class="post-metadata">

**Author:** ![ching-wong](https://yyz2.discourse-cdn.com/free1/user_avatar/igraph.discourse.group/ching-wong/32/1130_2.png) [@ching-wong](https://igraph.discourse.group/u/ching-wong)\
**Post date:** [19 November 2024 11:26 UTC](https://igraph.discourse.group/t/is-prime-graph-decomposition-implemented/1895/8 "2024-11-19T11:26:45Z")

</div>

Thank you szhorvat.

I implemented an algorithm presented in a recent paper [https://epubs.siam.org/doi/10.1137/16M1059837](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.

```auto
###### 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))
}

```

```auto
###### 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.

---

<div class="post-metadata">

**Author:** ![ching-wong](https://yyz2.discourse-cdn.com/free1/user_avatar/igraph.discourse.group/ching-wong/32/1130_2.png) [@ching-wong](https://igraph.discourse.group/u/ching-wong)\
**Post date:** [20 November 2024 08:31 UTC](https://igraph.discourse.group/t/is-prime-graph-decomposition-implemented/1895/9 "2024-11-20T08:31:12Z")

</div>

It turns out that is\_chordal() does not give a minimal chordal completion. I found a counter example:

```auto
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).

---

<div class="post-metadata">

**Author:** ![szhorvat](https://yyz2.discourse-cdn.com/free1/user_avatar/igraph.discourse.group/szhorvat/32/3_2.png) [@szhorvat](https://igraph.discourse.group/u/szhorvat)\
**Post date:** [20 November 2024 09:22 UTC](https://igraph.discourse.group/t/is-prime-graph-decomposition-implemented/1895/10 "2024-11-20T09:22:47Z")

</div>

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.

> [@ching-wong](#):
>
> 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).

The current algorithm is quite efficient. I’m in favour of writing a new function for a _minimal_ chordal completion.

---

<div class="post-metadata">

**Author:** ![ching-wong](https://yyz2.discourse-cdn.com/free1/user_avatar/igraph.discourse.group/ching-wong/32/1130_2.png) [@ching-wong](https://igraph.discourse.group/u/ching-wong)\
**Post date:** [28 November 2024 13:05 UTC](https://igraph.discourse.group/t/is-prime-graph-decomposition-implemented/1895/11 "2024-11-28T13:05:22Z")

</div>

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.

```auto
#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()?

```auto
======================================
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.

---

<div class="post-metadata">

**Author:** ![szhorvat](https://yyz2.discourse-cdn.com/free1/user_avatar/igraph.discourse.group/szhorvat/32/3_2.png) [@szhorvat](https://igraph.discourse.group/u/szhorvat)\
**Post date:** [29 November 2024 10:40 UTC](https://igraph.discourse.group/t/is-prime-graph-decomposition-implemented/1895/12 "2024-11-29T10:40:43Z")

</div>

Could you please open a pull request at [GitHub - igraph/igraph: Library for the analysis of networks](https://github.com/igraph/igraph/) ?

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.

---

<div class="post-metadata">

**Author:** ![ching-wong](https://yyz2.discourse-cdn.com/free1/user_avatar/igraph.discourse.group/ching-wong/32/1130_2.png) [@ching-wong](https://igraph.discourse.group/u/ching-wong)\
**Post date:** [2 December 2024 11:18 UTC](https://igraph.discourse.group/t/is-prime-graph-decomposition-implemented/1895/13 "2024-12-02T11:18:06Z")

</div>

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.

```auto
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))
}

```