[ENH] Adding acyclification procedure (#17)

* Remove license in doc/index.rst (#14)
* Adding acyclification algorithm

Signed-off-by: Adam Li <[email protected]>

docs/api.rst

@@ -46,6 +46,7 @@ causal graph operations.
    pds
    pds_path
    uncovered_pd_path
+   acyclification
 
 Conversions between other package's causal graphs
 =================================================

docs/references.bib

@@ -46,6 +46,24 @@ @article{Meek1995
   journal = {Proceedings of Eleventh Conference on Uncertainty in Artificial Intelligence, Montreal, QU}
 }
 
+
+@InProceedings{Mooij2020cyclic,
+  title = 	 {Constraint-Based Causal Discovery using Partial Ancestral Graphs in the presence of Cycles},
+  author =       {M. Mooij, Joris and Claassen, Tom},
+  booktitle = 	 {Proceedings of the 36th Conference on Uncertainty in Artificial Intelligence (UAI)},
+  pages = 	 {1159--1168},
+  year = 	 {2020},
+  editor = 	 {Peters, Jonas and Sontag, David},
+  volume = 	 {124},
+  series = 	 {Proceedings of Machine Learning Research},
+  month = 	 {03--06 Aug},
+  publisher =    {PMLR},
+  pdf = 	 {http://proceedings.mlr.press/v124/m-mooij20a/m-mooij20a.pdf},
+  url = 	 {https://proceedings.mlr.press/v124/m-mooij20a.html},
+  abstract = 	 {While feedback loops are known to play important roles in many complex systems, their existence is ignored in a large part of the causal discovery literature, as systems are typically assumed to be acyclic from the outset. When applying causal discovery algorithms designed for the acyclic setting on data generated by a system that involves feedback, one would not expect to obtain correct results. In this work, we show that—surprisingly—the output of the Fast Causal Inference (FCI) algorithm is correct if it is applied to observational data generated by a system that involves feedback. More specifically, we prove that for observational data generated by a simple and sigma-faithful Structural Causal Model (SCM), FCI is sound and complete, and can be used to consistently estimate (i) the presence and absence of causal relations, (ii) the presence and absence of direct causal relations, (iii) the absence of confounders, and (iv) the absence of specific cycles in the causal graph of the SCM. We extend these results to constraint-based causal discovery algorithms that exploit certain forms of background knowledge, including the causally sufficient setting (e.g., the PC algorithm) and the Joint Causal Inference setting (e.g., the FCI-JCI algorithm).}
+}
+
+
 @book{Neapolitan2003,
   author    = {Neapolitan, Richard},
   year      = {2003},

docs/whats_new/v0.1.rst

@@ -30,6 +30,7 @@ Changelog
 - |Feature| Implement and test the :class:`pywhy_graphs.PAG` for PAGs, by `Adam Li`_ (:pr:`9`)
 - |Feature| Implement and test various PAG algorithms :func:`pywhy_graphs.algorithms.possible_ancestors`, :func:`pywhy_graphs.algorithms.possible_descendants`, :func:`pywhy_graphs.algorithms.discriminating_path`, :func:`pywhy_graphs.algorithms.pds`, :func:`pywhy_graphs.algorithms.pds_path`, and :func:`pywhy_graphs.algorithms.uncovered_pd_path`, by `Adam Li`_ (:pr:`10`)
 - |Feature| Implement an array API wrapper to convert between causal graphs in pywhy-graphs and causal graphs in ``causal-learn``, by `Adam Li`_ (:pr:`16`)
+- |Feature| Implement an acyclification algorithm for converting cyclic graphs to acyclic with :func:`pywhy_graphs.algorithms.acyclification`, by `Adam Li`_ (:pr:`17`)
 
 Code and Documentation Contributors
 -----------------------------------

pywhy_graphs/algorithms/__init__.py

@@ -1,2 +1,3 @@
+from .cyclic import *  # noqa: F403
 from .generic import *  # noqa: F403
 from .pag import *  # noqa: F403

pywhy_graphs/algorithms/cyclic.py

@@ -0,0 +1,100 @@
+import networkx as nx
+
+
+def acyclification(
+    G: nx.MixedEdgeGraph,
+    directed_edge_type: str = "directed",
+    bidirected_edge_type: str = "bidirected",
+    copy: bool = True,
+) -> nx.MixedEdgeGraph:
+    """Acyclify a cyclic graph.
+
+    Applies the acyclification procedure presented in :footcite:`Mooij2020cyclic`.
+    This converts to G to what is called :math:`G^{acy}` in the reference.
+
+    Parameters
+    ----------
+    G : nx.MixedEdgeGraph
+        A graph with cycles.
+    directed_edge_type : str
+        The name of the sub-graph of directed edges.
+    bidirected_edge_type : str
+        The name of the sub-graph of bidirected edges.
+    copy : bool
+        Whether to operate on the graph in place, or make a copy.
+
+    Returns
+    -------
+    G : nx.MixedEdgeGraph
+        The acyclified graph.
+
+    Notes
+    -----
+    This takes
+    This replaces all strongly connected components of G by fully connected
+    bidirected components without any directed edges. Then any node with an
+    edge pointing into the SC (i.e. a directed edge, or bidirected edge) is
+    made fully connected with the nodes of the SC either with a directed, or
+    bidirected edge.
+
+    References
+    ----------
+    .. footbibliography::
+    """
+    if copy:
+        G = G.copy()
+
+    # extract the subgraph of directed edges
+    directed_G: nx.DiGraph = G.get_graphs(directed_edge_type).copy()
+    bidirected_G: nx.Graph = G.get_graphs(bidirected_edge_type).copy()
+
+    # first detect all strongly connected components
+    scomps = nx.strongly_connected_components(directed_G)
+
+    # loop over all strongly connected components and their nodes
+    for comp in scomps:
+        if len(comp) == 1:
+            continue
+
+        # extract the parents, or c-components of any node
+        # in the strongly-connected component
+        scomp_parents = set()
+        scomp_c_components = set()
+        scomp_children = []
+
+        for node in comp:
+            # get any predecessors of SC
+            for parent in directed_G.predecessors(node):
+                if parent in comp:
+                    continue
+                scomp_parents.add(parent)
+
+            # get any bidirected edges pointing to elements of SC
+            for nbr in bidirected_G.neighbors(node):
+                if nbr in comp:
+                    continue
+                scomp_c_components.add(nbr)
+
+            # keep track of any edges pointing out of the SC
+            for child in directed_G.successors(node):
+                if child in comp:
+                    continue
+                scomp_children.append((node, child))
+
+        # first remove all nodes in the cycle
+        G.remove_nodes_from(comp)
+
+        # add them back in as a fully connected bidirected graph
+        bidirected_fc_G = nx.complete_graph(comp)
+        G.add_edges_from(bidirected_fc_G.edges, bidirected_edge_type)
+
+        # add back the children
+        G.add_edges_from(scomp_children, directed_edge_type)
+
+        # make all variables connect to the strongly connected component
+        for node in comp:
+            for parent in scomp_parents:
+                G.add_edge(parent, node, directed_edge_type)
+            for c_component in scomp_c_components:
+                G.add_edge(c_component, node, bidirected_edge_type)
+    return G

pywhy_graphs/algorithms/tests/test_cyclic.py

@@ -0,0 +1,60 @@
+import networkx as nx
+
+import pywhy_graphs
+
+
+def test_acyclification():
+    """Test acyclification procedure as specified in :footcite:`Mooij2020cyclic`.
+
+    Tests the graphs as presented in Figure 2.
+    """
+    directed_edges = nx.DiGraph(
+        [
+            ("x8", "x2"),
+            ("x9", "x2"),
+            ("x10", "x1"),
+            ("x2", "x4"),
+            ("x4", "x6"),  # start of cycle
+            ("x6", "x5"),
+            ("x5", "x3"),
+            ("x3", "x4"),  # end of cycle
+            ("x6", "x7"),
+        ]
+    )
+    bidirected_edges = nx.Graph([("x1", "x3")])
+    G = nx.MixedEdgeGraph([directed_edges, bidirected_edges], ["directed", "bidirected"])
+    acyclic_G = pywhy_graphs.acyclification(G)
+
+    directed_edges = nx.DiGraph(
+        [
+            ("x8", "x2"),
+            ("x9", "x2"),
+            ("x10", "x1"),
+            ("x2", "x4"),
+            ("x6", "x7"),
+            ("x2", "x3"),
+            ("x2", "x5"),
+            ("x2", "x4"),
+            ("x2", "x6"),
+        ]
+    )
+    bidirected_edges = nx.Graph(
+        [
+            ("x1", "x3"),
+            ("x4", "x6"),
+            ("x6", "x5"),
+            ("x5", "x3"),
+            ("x3", "x4"),
+            ("x4", "x5"),
+            ("x3", "x6"),
+            ("x1", "x3"),
+            ("x1", "x5"),
+            ("x1", "x4"),
+            ("x1", "x6"),
+        ]
+    )
+    expected_G = nx.MixedEdgeGraph([directed_edges, bidirected_edges], ["directed", "bidirected"])
+
+    for edge_type, graph in acyclic_G.get_graphs().items():
+        expected_graph = expected_G.get_graphs(edge_type)
+        assert nx.is_isomorphic(graph, expected_graph)

pywhy_graphs/algorithms/tests/test_pag.py

@@ -177,6 +177,7 @@ def test_discriminating_path():
 
 
 def test_uncovered_pd_path():
+    """Test basic uncovered partially directed path."""
     # If A o-> C and there is an undirected pd path
     # from A to C through u, where u and C are not adjacent
     # then orient A o-> C as A -> C
@@ -221,6 +222,7 @@ def test_uncovered_pd_path():
 
 
 def test_uncovered_pd_path_intersecting():
+    """Test basic uncovered partially directed path with intersecting paths."""
     G = pywhy_graphs.PAG()
 
     # make A o-> C