Repository: pgmpy/pgmpy

Issue:
PC algorithm orientation rules throw error
Code to reproduce:
```python
In [1]: from pgmpy.estimators import PC                                                                                                                                                       
                                                                                                                                                                                    
In [2]: from pgmpy.utils import get_example_model                                                                                                                                             
                                                                                                                                                                                              
In [3]: model = get_example_model('ecoli70')                                                                                                                                                  
                                                                                                                                                                                              
In [4]: df = model.simulate(int(1e3))                                                                                                                                                         
                                                                                                                                                                                              
In [5]: import numpy as np                                                                                                                                                                    
                                                                                                                                                                                              
In [6]: cpdags = []                                                                                                                                                                           
                                                                                                                                                                                              
In [7]: for i in range(10):                                                                                                                                                                   
   ...:     est = PC(df.loc[:, np.random.choice(df.columns, 20, replace=False)]).estimate(ci_te                                                                                               
   ...: st='pearsonr')                                                                                                                                                                        
   ...:     cpdags.append(est)                                                                                                                                                                
   ...:                                                                                                                                                                                       
INFO:pgmpy: Datatype (N=numerical, C=Categorical Unordered, O=Categorical Ordered) inferred from data:                                                                                        
 {'yheI': 'N', 'nuoM': 'N', 'asnA': 'N', 'cchB': 'N', 'ftsJ': 'N', 'yceP': 'N', 'icdA': 'N', 'cspA': 'N', 'yecO': 'N', 'flgD': 'N', 'aceB': 'N', 'ibpB': 'N', 'yfaD': 'N', 'lacY': 'N', 'pspA'
: 'N', 'yedE': 'N', 'yfiA': 'N', 'lacZ': 'N', 'sucD': 'N', 'dnaJ': 'N'}                                                                                                                       
Working for n conditional variables: 5: 100%|████████████████████| 5/5 [00:11<00:00,  1.44s/it]INFO:pgmpy:Reached maximum number of allowed conditional variables. Exiting                    
Working for n conditional variables: 5: 100%|████████████████████| 5/5 [00:12<00:00,  2.44s/it]       

KeyError                                  Traceback (most recent call last)
File ~/miniconda3/envs/py312/lib/python3.12/site-packages/networkx/classes/digraph.py:848, in DiGraph.remove_edge(self, u, v)
    847 try:
--> 848     del self._succ[u][v]
    849     del self._pred[v][u]

KeyError: 'sucD'

The above exception was the direct cause of the following exception:

NetworkXError                             Traceback (most recent call last)
Cell In[7], line 2
      1 for i in range(10):
----> 2     est = PC(df.loc[:, np.random.choice(df.columns, 20, replace=False)]).estimate(ci_test='pearsonr')
      3     cpdags.append(est)

File ~/pgmpy/pgmpy/estimators/PC.py:229, in PC.estimate(self, variant, ci_test, return_type, significance_level, max_cond_vars, expert_knowledge, enforce_expert_knowledge, n_jobs, show_progress, **kwargs)
    227 if not enforce_expert_knowledge:
    228     pdag = expert_knowledge.apply_expert_knowledge(pdag)
--> 229     pdag = self.apply_orientation_rules(pdag, apply_r4=True)
    231 if self.data is not None:
    232     pdag.add_nodes_from(set(self.data.columns) - set(pdag.nodes()))

File ~/pgmpy/pgmpy/estimators/PC.py:605, in PC.apply_orientation_rules(pdag, apply_r4)
    596             for Z in (
    597                 set(pdag.successors(X))
    598                 & set(pdag.predecessors(X))
    599                 & set(pdag.predecessors(Y))
    600                 & set(pdag.successors(Y))
    601             ):
    602                 for W in (
    603                     set(pdag.successors(Y)) - set(pdag.predecessors(Y))
    604                 ) & (set(pdag.predecessors(Z)) | set(pdag.successors(Z))):
--> 605                     pdag.remove_edge(X, Z)
    607     progress = num_edges > pdag.number_of_edges()
    609 edges = set(pdag.edges())

File ~/miniconda3/envs/py312/lib/python3.12/site-packages/networkx/classes/digraph.py:851, in DiGraph.remove_edge(self, u, v)
    849     del self._pred[v][u]
    850 except KeyError as err:
--> 851     raise NetworkXError(f"The edge {u}-{v} not in graph.") from err
    852 nx._clear_cache(self)

NetworkXError: The edge dnaJ-sucD not in graph.
```

Discussion:

Solution patch:
diff --git a/pgmpy/base/DAG.py b/pgmpy/base/DAG.py
index 676c471a..eee7807e 100644
--- a/pgmpy/base/DAG.py
+++ b/pgmpy/base/DAG.py
@@ -405,7 +405,7 @@ class DAG(nx.DiGraph):
         >>> from pgmpy.base import DAG
         >>> chain = DAG([('X', 'Y'), ('Y', 'Z')])
         >>> chain.get_independencies()
-        (X \u27C2 Z | Y)
+        (X \u27c2 Z | Y)
         """
         nodes = set(self.nodes())
         if not include_latents:
@@ -445,7 +445,7 @@ class DAG(nx.DiGraph):
         >>>                         ('grade', 'letter'), ('intel', 'SAT')])
         >>> ind = student.local_independencies('grade')
         >>> ind
-        (grade \u27C2 SAT | diff, intel)
+        (grade \u27c2 SAT | diff, intel)
         """
 
         independencies = Independencies()
diff --git a/pgmpy/estimators/CITests.py b/pgmpy/estimators/CITests.py
index 7891577f..dd8f7e52 100644
--- a/pgmpy/estimators/CITests.py
+++ b/pgmpy/estimators/CITests.py
@@ -46,19 +46,19 @@ def get_ci_test(test, full=False, data=None, independencies=None):
 
 def independence_match(X, Y, Z, independencies, **kwargs):
     """
-    Checks if `X \u27C2 Y | Z` is in `independencies`. This method is implemented to
+    Checks if `X \u27c2 Y | Z` is in `independencies`. This method is implemented to
     have an uniform API when the independencies are provided instead of data.
 
     Parameters
     ----------
     X: str
-        The first variable for testing the independence condition X \u27C2 Y | Z
+        The first variable for testing the independence condition X \u27c2 Y | Z
 
     Y: str
-        The second variable for testing the independence condition X \u27C2 Y | Z
+        The second variable for testing the independence condition X \u27c2 Y | Z
 
     Z: list/array-like
-        A list of conditional variable for testing the condition X \u27C2 Y | Z
+        A list of conditional variable for testing the condition X \u27c2 Y | Z
 
     data: pandas.DataFrame The dataset in which to test the indepenedence condition.
 
@@ -108,7 +108,7 @@ def chi_square(X, Y, Z, data, boolean=True, **kwargs):
         If boolean = False, Returns a tuple (chi, p_value, dof). `chi` is the
         chi-squared test statistic. The `p_value` for the test, i.e. the
         probability of observing the computed chi-square statistic (or an even
-        higher value), given the null hypothesis that X \u27C2 Y | Zs is True.
+        higher value), given the null hypothesis that X \u27c2 Y | Zs is True.
         If boolean = True, returns True if the p_value of the test is greater
         than `significance_level` else returns False.
 
@@ -168,7 +168,7 @@ def g_sq(X, Y, Z, data, boolean=True, **kwargs):
         If boolean = False, Returns a tuple (chi, p_value, dof). `chi` is the
         chi-squared test statistic. The `p_value` for the test, i.e. the
         probability of observing the computed chi-square statistic (or an even
-        higher value), given the null hypothesis that X \u27C2 Y | Zs is True.
+        higher value), given the null hypothesis that X \u27c2 Y | Zs is True.
         If boolean = True, returns True if the p_value of the test is greater
         than `significance_level` else returns False.
 
@@ -228,7 +228,7 @@ def log_likelihood(X, Y, Z, data, boolean=True, **kwargs):
         If boolean = False, Returns a tuple (chi, p_value, dof). `chi` is the
         chi-squared test statistic. The `p_value` for the test, i.e. the
         probability of observing the computed chi-square statistic (or an even
-        higher value), given the null hypothesis that X \u27C2 Y | Zs is True.
+        higher value), given the null hypothesis that X \u27c2 Y | Zs is True.
         If boolean = True, returns True if the p_value of the test is greater
         than `significance_level` else returns False.
 
@@ -287,7 +287,7 @@ def modified_log_likelihood(X, Y, Z, data, boolean=True, **kwargs):
         If boolean = False, Returns a tuple (chi, p_value, dof). `chi` is the
         chi-squared test statistic. The `p_value` for the test, i.e. the
         probability of observing the computed chi-square statistic (or an even
-        higher value), given the null hypothesis that X \u27C2 Y | Zs is True.
+        higher value), given the null hypothesis that X \u27c2 Y | Zs is True.
         If boolean = True, returns True if the p_value of the test is greater
         than `significance_level` else returns False.
 
@@ -361,7 +361,7 @@ def power_divergence(X, Y, Z, data, boolean=True, lambda_="cressie-read", **kwar
         If boolean = False, Returns a tuple (chi, p_value, dof). `chi` is the
         chi-squared test statistic. The `p_value` for the test, i.e. the
         probability of observing the computed chi-square statistic (or an even
-        higher value), given the null hypothesis that X \u27C2 Y | Zs is True.
+        higher value), given the null hypothesis that X \u27c2 Y | Zs is True.
         If boolean = True, returns True if the p_value of the test is greater
         than `significance_level` else returns False.
 
@@ -417,14 +417,14 @@ def power_divergence(X, Y, Z, data, boolean=True, lambda_="cressie-read", **kwar
             if any(contingency.sum(axis=0) == 0) or any(contingency.sum(axis=1) == 0):
                 if isinstance(z_state, str):
                     logger.info(
-                        f"Skipping the test {X} \u27C2 {Y} | {Z[0]}={z_state}. Not enough samples"
+                        f"Skipping the test {X} \u27c2 {Y} | {Z[0]}={z_state}. Not enough samples"
                     )
                 else:
                     z_str = ", ".join(
                         [f"{var}={state}" for var, state in zip(Z, z_state)]
                     )
                     logger.info(
-                        f"Skipping the test {X} \u27C2 {Y} | {z_str}. Not enough samples"
+                        f"Skipping the test {X} \u27c2 {Y} | {z_str}. Not enough samples"
                     )
             else:
                 c, _, d, _ = stats.chi2_contingency(contingency, lambda_=lambda_)
@@ -448,13 +448,13 @@ def pearsonr(X, Y, Z, data, boolean=True, **kwargs):
     Parameters
     ----------
     X: str
-        The first variable for testing the independence condition X \u27C2 Y | Z
+        The first variable for testing the independence condition X \u27c2 Y | Z
 
     Y: str
-        The second variable for testing the independence condition X \u27C2 Y | Z
+        The second variable for testing the independence condition X \u27c2 Y | Z
 
     Z: list/array-like
-        A list of conditional variable for testing the condition X \u27C2 Y | Z
+        A list of conditional variable for testing the condition X \u27c2 Y | Z
 
     data: pandas.DataFrame
         The dataset in which to test the indepenedence condition.
@@ -586,13 +586,13 @@ def pillai_trace(X, Y, Z, data, boolean=True, **kwargs):
     Parameters
     ----------
     X: str
-        The first variable for testing the independence condition X \u27C2 Y | Z
+        The first variable for testing the independence condition X \u27c2 Y | Z
 
     Y: str
-        The second variable for testing the independence condition X \u27C2 Y | Z
+        The second variable for testing the independence condition X \u27c2 Y | Z
 
     Z: list/array-like
-        A list of conditional variable for testing the condition X \u27C2 Y | Z
+        A list of conditional variable for testing the condition X \u27c2 Y | Z
 
     data: pandas.DataFrame
         The dataset in which to test the indepenedence condition.
diff --git a/pgmpy/estimators/PC.py b/pgmpy/estimators/PC.py
index 6759cd80..0debb2c5 100644
--- a/pgmpy/estimators/PC.py
+++ b/pgmpy/estimators/PC.py
@@ -394,6 +394,24 @@ class PC(StructureEstimator):
             pbar.close()
         return graph, separating_sets
 
+    @staticmethod
+    def _check_incoming_edges(pdag, u, v):
+        "Used for checking whether a new v-structure is getting formed"
+        for predecessor in pdag.predecessors(v):
+            if (
+                not pdag.has_edge(
+                    v, predecessor
+                )  # this ignores bidirected edges of 'v'
+                and not pdag.has_edge(
+                    predecessor, u
+                )  # prevent the case (by returning true) when a new unshielded
+                and not pdag.has_edge(
+                    u, predecessor
+                )  # collider may form at 'v'  i.e. predecessor--> v <--u
+            ):
+                return True
+        return False
+
     @staticmethod
     def orient_colliders(skeleton, separating_sets):
         """
@@ -518,11 +536,21 @@ class PC(StructureEstimator):
             # (Explanation in Koller & Friedman PGM, page 88)
             for pair in node_pairs:
                 X, Y = pair
-                if not pdag.has_edge(X, Y):
+                if not pdag.has_edge(X, Y) and not pdag.has_edge(Y, X):
                     for Z in (set(pdag.successors(X)) - set(pdag.predecessors(X))) & (
                         set(pdag.successors(Y)) & set(pdag.predecessors(Y))
                     ):
-                        pdag.remove_edge(Y, Z)
+                        if not PC._check_incoming_edges(pdag, Z, Y):
+                            any_directed = False
+                            for path in nx.all_simple_paths(pdag, Y, Z):
+                                is_directed = True
+                                for src, dst in list(zip(path, path[1:])):
+                                    if pdag.has_edge(dst, src):
+                                        is_directed = False
+                                if is_directed:
+                                    any_directed = True
+                            if not any_directed:
+                                pdag.remove_edge(Y, Z)
 
             # 2) for each X-Y with a directed path from X to Y, orient edges to X->Y
             for pair in node_pairs:
@@ -551,6 +579,7 @@ class PC(StructureEstimator):
                         & (set(pdag.successors(Y)) - set(pdag.predecessors(Y)))
                         & (set(pdag.successors(Z)) & set(pdag.predecessors(Z)))
                     ):
+                        # if not PC.check_incoming_edges(pdag, W):
                         pdag.remove_edge(W, Z)
 
             # This rule (rule 4 in Meek's rules) is only used in the case of a
@@ -570,8 +599,10 @@ class PC(StructureEstimator):
                         & set(pdag.successors(Y))
                     ):
                         for W in (
-                            set(pdag.successors(Y)) - set(pdag.predecessors(Y))
-                        ) & (set(pdag.predecessors(Z)) | set(pdag.successors(Z))):
+                            (set(pdag.successors(Y)) - set(pdag.predecessors(Y)))
+                            & (set(pdag.predecessors(Z)) | set(pdag.successors(Z)))
+                            & set(pdag.predecessors(X))
+                        ):
                             pdag.remove_edge(X, Z)
 
             progress = num_edges > pdag.number_of_edges()
diff --git a/pgmpy/factors/discrete/JointProbabilityDistribution.py b/pgmpy/factors/discrete/JointProbabilityDistribution.py
index 615c6d4d..af1f9f2d 100644
--- a/pgmpy/factors/discrete/JointProbabilityDistribution.py
+++ b/pgmpy/factors/discrete/JointProbabilityDistribution.py
@@ -227,9 +227,9 @@ class JointProbabilityDistribution(DiscreteFactor):
         >>> from pgmpy.factors.discrete import JointProbabilityDistribution
         >>> prob = JointProbabilityDistribution(['x1', 'x2', 'x3'], [2, 3, 2], np.ones(12)/12)
         >>> prob.get_independencies()
-        (x1 \u27C2 x2)
-        (x1 \u27C2 x3)
-        (x2 \u27C2 x3)
+        (x1 \u27c2 x2)
+        (x1 \u27c2 x3)
+        (x2 \u27c2 x3)
         """
         JPD = self.copy()
         if condition:
diff --git a/pgmpy/independencies/Independencies.py b/pgmpy/independencies/Independencies.py
index ca0303b8..27e6512a 100644
--- a/pgmpy/independencies/Independencies.py
+++ b/pgmpy/independencies/Independencies.py
@@ -163,23 +163,23 @@ class Independencies(object):
         >>> from pgmpy.independencies import Independencies
         >>> ind1 = Independencies(('A', ['B', 'C'], 'D'))
         >>> ind1.closure()
-        (A \u27C2 B | D, C)
-        (A \u27C2 B, C | D)
-        (A \u27C2 B | D)
-        (A \u27C2 C | D, B)
-        (A \u27C2 C | D)
+        (A \u27c2 B | D, C)
+        (A \u27c2 B, C | D)
+        (A \u27c2 B | D)
+        (A \u27c2 C | D, B)
+        (A \u27c2 C | D)
 
         >>> ind2 = Independencies(('W', ['X', 'Y', 'Z']))
         >>> ind2.closure()
-        (W \u27C2 Y)
-        (W \u27C2 Y | X)
-        (W \u27C2 Z | Y)
-        (W \u27C2 Z, X, Y)
-        (W \u27C2 Z)
-        (W \u27C2 Z, X)
-        (W \u27C2 X, Y)
-        (W \u27C2 Z | X)
-        (W \u27C2 Z, Y | X)
+        (W \u27c2 Y)
+        (W \u27c2 Y | X)
+        (W \u27c2 Z | Y)
+        (W \u27c2 Z, X, Y)
+        (W \u27c2 Z)
+        (W \u27c2 Z, X)
+        (W \u27c2 X, Y)
+        (W \u27c2 Z | X)
+        (W \u27c2 Z, Y | X)
         [..]
         """
 
@@ -464,13 +464,13 @@ class IndependenceAssertion(object):
 
     def __str__(self):
         if self.event3:
-            return "({event1} \u27C2 {event2} | {event3})".format(
+            return "({event1} \u27c2 {event2} | {event3})".format(
                 event1=", ".join([str(e) for e in self.event1]),
                 event2=", ".join([str(e) for e in self.event2]),
                 event3=", ".join([str(e) for e in self.event3]),
             )
         else:
-            return "({event1} \u27C2 {event2})".format(
+            return "({event1} \u27c2 {event2})".format(
                 event1=", ".join([str(e) for e in self.event1]),
                 event2=", ".join([str(e) for e in self.event2]),
             )
diff --git a/pgmpy/models/NaiveBayes.py b/pgmpy/models/NaiveBayes.py
index 2e71b4cb..59239e10 100644
--- a/pgmpy/models/NaiveBayes.py
+++ b/pgmpy/models/NaiveBayes.py
@@ -159,7 +159,7 @@ class NaiveBayes(BayesianNetwork):
         >>> model.add_edges_from([('a', 'b'), ('a', 'c'), ('a', 'd')])
         >>> ind = model.local_independencies('b')
         >>> ind
-        (b \u27C2 d, c | a)
+        (b \u27c2 d, c | a)
         """
         independencies = Independencies()
         for variable in [variables] if isinstance(variables, str) else variables:


Test patch:
diff --git a/pgmpy/tests/test_estimators/test_PC.py b/pgmpy/tests/test_estimators/test_PC.py
index cfebc4d2..ba7dd411 100644
--- a/pgmpy/tests/test_estimators/test_PC.py
+++ b/pgmpy/tests/test_estimators/test_PC.py
@@ -29,11 +29,11 @@ class TestPCFakeCITest(unittest.TestCase):
         """
         A mock CI testing function which gives False for every condition
         except for the following:
-            1. B \u27C2 C
-            2. B \u27C2 D
-            3. C \u27C2 D
-            4. A \u27C2 B | C
-            5. A \u27C2 C | B
+            1. B \u27c2 C
+            2. B \u27c2 D
+            3. C \u27c2 D
+            4. A \u27c2 B | C
+            5. A \u27c2 C | B
         """
         Z = list(Z)
         if X == "B":
@@ -91,7 +91,7 @@ class TestPCFakeCITest(unittest.TestCase):
             self.assertTrue(((u, v) in expected_edges) or ((v, u) in expected_edges))
 
 
-class TestPCEstimatorFromIndependencies(unittest.TestCase):
+class TestPCEstimatorFromIndependences(unittest.TestCase):
     def test_build_skeleton_from_ind(self):
         # Specify a set of independencies
         for variant in ["orig", "stable", "parallel"]:
@@ -156,6 +156,7 @@ class TestPCEstimatorFromIndependencies(unittest.TestCase):
             )
 
     def test_skeleton_to_pdag(self):
+        # D - A - C - B  ==> D - A -> C <- B
         skel = nx.Graph([("A", "D"), ("A", "C"), ("B", "C")])
         sep_sets = {
             frozenset({"D", "C"}): ("A",),
@@ -168,6 +169,7 @@ class TestPCEstimatorFromIndependencies(unittest.TestCase):
             set(pdag.edges()), set([("B", "C"), ("A", "D"), ("A", "C"), ("D", "A")])
         )
 
+        # C - A - B  ==> C -> A <- B
         skel = nx.Graph([("A", "B"), ("A", "C")])
         sep_sets = {frozenset({"B", "C"}): ()}
         pdag = PC.orient_colliders(skeleton=skel, separating_sets=sep_sets)
@@ -177,6 +179,7 @@ class TestPCEstimatorFromIndependencies(unittest.TestCase):
             set([("B", "A"), ("C", "A")]),
         )
 
+        # C - A - B ==> C - A - B
         sep_sets = {frozenset({"B", "C"}): ("A",)}
         pdag = PC.orient_colliders(skeleton=skel, separating_sets=sep_sets)
         pdag = PC.apply_orientation_rules(pdag)
@@ -185,6 +188,7 @@ class TestPCEstimatorFromIndependencies(unittest.TestCase):
             set([("A", "B"), ("B", "A"), ("A", "C"), ("C", "A")]),
         )
 
+        # {A, B} - C - D ==> {A, B} -> C -> D
         skel = nx.Graph([("A", "C"), ("B", "C"), ("C", "D")])
         sep_sets = {
             frozenset({"A", "B"}): tuple(),
@@ -197,6 +201,7 @@ class TestPCEstimatorFromIndependencies(unittest.TestCase):
             set(pdag.edges()), set([("A", "C"), ("B", "C"), ("C", "D")])
         )
 
+        # C - A - B - {C, D} ==> C <- A -> B <- D; B -> C
         skel = nx.Graph([("A", "B"), ("A", "C"), ("B", "C"), ("B", "D")])
         sep_sets = {frozenset({"A", "D"}): tuple(), frozenset({"C", "D"}): ("A", "B")}
         pdag = PC.orient_colliders(skeleton=skel, separating_sets=sep_sets)
@@ -224,6 +229,143 @@ class TestPCEstimatorFromIndependencies(unittest.TestCase):
             ),
         )
 
+    def test_pdag_to_cpdag(self):
+        pdag = PDAG(directed_ebunch=[("A", "B")], undirected_ebunch=[("B", "C")])
+        cpdag = PC.apply_orientation_rules(pdag, apply_r4=True)
+        self.assertSetEqual(set(cpdag.edges()), {("A", "B"), ("B", "C")})
+
+        pdag = PDAG(
+            directed_ebunch=[("A", "B")], undirected_ebunch=[("B", "C"), ("C", "D")]
+        )
+        cpdag = PC.apply_orientation_rules(pdag, apply_r4=True)
+        self.assertSetEqual(set(cpdag.edges()), {("A", "B"), ("B", "C"), ("C", "D")})
+
+        pdag = PDAG(
+            directed_ebunch=[("A", "B"), ("D", "C")], undirected_ebunch=[("B", "C")]
+        )
+        cpdag = PC.apply_orientation_rules(pdag, apply_r4=True)
+        self.assertSetEqual(
+            set(cpdag.edges()), {("A", "B"), ("D", "C"), ("B", "C"), ("C", "B")}
+        )
+
+        pdag = PDAG(
+            directed_ebunch=[("A", "B"), ("D", "C"), ("D", "B")],
+            undirected_ebunch=[("B", "C")],
+        )
+        cpdag = PC.apply_orientation_rules(pdag, apply_r4=True)
+        self.assertSetEqual(
+            set(cpdag.edges()), {("A", "B"), ("D", "C"), ("D", "B"), ("B", "C")}
+        )
+
+        pdag = PDAG(
+            directed_ebunch=[("A", "B"), ("B", "C")], undirected_ebunch=[("A", "C")]
+        )
+        cpdag = PC.apply_orientation_rules(pdag, apply_r4=True)
+        self.assertSetEqual(set(cpdag.edges()), {("A", "B"), ("B", "C"), ("A", "C")})
+
+        pdag = PDAG(
+            directed_ebunch=[("A", "B"), ("B", "C"), ("D", "C")],
+            undirected_ebunch=[("A", "C")],
+        )
+        cpdag = PC.apply_orientation_rules(pdag, apply_r4=True)
+        self.assertSetEqual(
+            set(cpdag.edges()), {("A", "B"), ("B", "C"), ("A", "C"), ("D", "C")}
+        )
+
+        # Examples taken from Perkovi\`c 2017.
+        pdag = PDAG(
+            directed_ebunch=[("V1", "X")],
+            undirected_ebunch=[("X", "V2"), ("V2", "Y"), ("X", "Y")],
+        )
+        cpdag = PC.apply_orientation_rules(pdag, apply_r4=True)
+        self.assertEqual(
+            set(cpdag.edges()),
+            {("V1", "X"), ("X", "V2"), ("X", "Y"), ("V2", "Y"), ("Y", "V2")},
+        )
+
+        pdag = PDAG(
+            directed_ebunch=[("Y", "X")],
+            undirected_ebunch=[("V1", "X"), ("X", "V2"), ("V2", "Y")],
+        )
+        cpdag = PC.apply_orientation_rules(pdag, apply_r4=True)
+        self.assertEqual(
+            set(cpdag.edges()),
+            {
+                ("X", "V1"),
+                ("Y", "X"),
+                ("X", "V2"),
+                ("V2", "X"),
+                ("V2", "Y"),
+                ("Y", "V2"),
+            },
+        )
+
+        # Exmaples from Bang 2024
+        pdag = PDAG(
+            directed_ebunch=[("B", "D"), ("C", "D")],
+            undirected_ebunch=[("A", "D"), ("A", "C")],
+        )
+        cpdag = PC.apply_orientation_rules(pdag, apply_r4=True)
+        self.assertEqual(
+            set(cpdag.edges()), {("B", "D"), ("D", "A"), ("C", "A"), ("C", "D")}
+        )
+
+        pdag = PDAG(
+            directed_ebunch=[("A", "B"), ("C", "B")],
+            undirected_ebunch=[("D", "B"), ("D", "A"), ("D", "C")],
+        )
+        cpdag = PC.apply_orientation_rules(pdag, apply_r4=True)
+        self.assertSetEqual(
+            set(cpdag.edges()),
+            {
+                ("A", "B"),
+                ("C", "B"),
+                ("D", "B"),
+                ("D", "A"),
+                ("A", "D"),
+                ("D", "C"),
+                ("C", "D"),
+            },
+        )
+
+        undirected_edges = [("A", "C"), ("B", "C"), ("D", "C")]
+        directed_edges = [("B", "D"), ("D", "A")]
+
+        pdag = PDAG(directed_ebunch=directed_edges, undirected_ebunch=undirected_edges)
+        mpdag = PC.apply_orientation_rules(pdag, apply_r4=True)
+        self.assertSetEqual(
+            set(mpdag.edges()),
+            set(
+                [
+                    ("C", "A"),
+                    ("C", "B"),
+                    ("B", "C"),
+                    ("B", "D"),
+                    ("D", "A"),
+                    ("D", "C"),
+                    ("C", "D"),
+                ]
+            ),
+        )
+
+        pdag = PDAG(directed_ebunch=directed_edges, undirected_ebunch=undirected_edges)
+        pdag = PC.apply_orientation_rules(pdag)
+        self.assertSetEqual(
+            set(pdag.edges()),
+            set(
+                [
+                    ("A", "C"),
+                    ("C", "A"),
+                    ("C", "B"),
+                    ("B", "C"),
+                    ("B", "D"),
+                    ("D", "A"),
+                    ("D", "C"),
+                    ("C", "D"),
+                ]
+            ),
+        )
+
     def test_estimate_dag(self):
         for variant in ["orig", "stable", "parallel"]:
             ind = Independencies(["B", "C"], ["A", ["B", "C"], "D"])
diff --git a/pgmpy/tests/test_factors/test_discrete/test_Factor.py b/pgmpy/tests/test_factors/test_discrete/test_Factor.py
index 5382dda2..ddf52bb0 100644
--- a/pgmpy/tests/test_factors/test_discrete/test_Factor.py
+++ b/pgmpy/tests/test_factors/test_discrete/test_Factor.py
@@ -509,7 +509,7 @@ class TestFactorMethods(unittest.TestCase):
 
     def test_reduce_scopeerror(self):
         self.assertRaises(ValueError, self.phi1.reduce, [("x4", 1)])
-        self.assertRaises(ValueError, self.phi5.reduce, [((("x1", 0.1), 0))])
+        self.assertRaises(ValueError, self.phi5.reduce, [(("x1", 0.1), 0)])
 
     def test_reduce_sizeerror(self):
         self.assertRaises(IndexError, self.phi1.reduce, [("x3", 5)])
diff --git a/pgmpy/tests/test_factors/test_discrete/test_FactorTorch.py b/pgmpy/tests/test_factors/test_discrete/test_FactorTorch.py
index 2489c745..6bb426d9 100644
--- a/pgmpy/tests/test_factors/test_discrete/test_FactorTorch.py
+++ b/pgmpy/tests/test_factors/test_discrete/test_FactorTorch.py
@@ -556,7 +556,7 @@ class TestFactorMethodsTorch(unittest.TestCase):
 
     def test_reduce_scopeerror(self):
         self.assertRaises(ValueError, self.phi1.reduce, [("x4", 1)])
-        self.assertRaises(ValueError, self.phi5.reduce, [((("x1", 0.1), 0))])
+        self.assertRaises(ValueError, self.phi5.reduce, [(("x1", 0.1), 0)])
 
     def test_reduce_sizeerror(self):
         self.assertRaises(IndexError, self.phi1.reduce, [("x3", 5)])