Repository: pgmpy/pgmpy

Issue:
Some bugs in PC -> skeleton_to_pdag
### Some bugs in PC -> skeleton_to_pdag
When I debug the module `PC`, I found that `node_pairs` in `skeleton_to_pdag` is defined as a `combination` iterator. So after first iteration, it will not yield any more. Check this answer:

https://stackoverflow.com/a/25336739

In another words, steps 2, 3, and 4 in `skeleton_to_pdag` are not performed in all circumstances. In addition, all codes inside of those `for` loop (steps 2, 3, and 4) have never been tested out. Bugs in there are not reflected by the tests results.

### Your environment
* pgmpy 0.1.12
* Python 3.6.11
* Operating System Ubuntu Linux x86_64

### Steps to reproduce
Using a simple example of `PC` estimation

```
from itertools import combinations
from pgmpy.estimators import PC
from pgmpy.utils import get_example_model
from pgmpy.sampling import BayesianModelSampling
import networkx as nx
import matplotlib.pyplot as plt

if __name__ == "__main__":
    model = get_example_model('asia')
    samples = BayesianModelSampling(model).forward_sample(size=int(1e5))
    est = PC(data=samples)
    dag = est.estimate(variant='stable', return_type="dag", max_cond_vars=5, significance_level=0.01)
    nx.draw_circular(dag, with_labels=True)
    plt.show()
```

### Expected behaviour
`node_pairs` in `PC` are desired to be iterated several times.

### Actual behaviour
In the first step's loop, everything is correct.
![Screenshot from 2020-12-08 09-01-02](https://user-images.githubusercontent.com/63798002/101494012-ca4c1b00-3934-11eb-8244-e5a1b42a2a8a.png)

In the second loop, `node_loop` will not yield, thus, step 2, 3, and 4 are all skipped in all cases. Some bugs inside of loop, too.
![Screenshot from 2020-12-08 09-02-40](https://user-images.githubusercontent.com/63798002/101494202-041d2180-3935-11eb-86b1-adc8747421d0.png)

I would like to open a pull request to fix this. Thanks!



Discussion:
ankurankan: @yzh211 Thanks a lot for the PR. Could you please also add some tests covering this fix?
ankurankan: @yzh211 Also, travis is complaining about the formatting: https://travis-ci.org/github/pgmpy/pgmpy/jobs/748475706. Running black (https://github.com/psf/black) on the file would fix it.
yzh211: @ankurankan Hi, thanks for the information! My bad for missing the tests! More bugs are fixed in the `skeleton_to_pdag`:

1. In step 2, there should be no direct path from X to Y. 
2. I replaced `combimations` to `permutations`. Because in some steps, it is directional sensitive for node pair (X, Y). Using combinations cannot cover all cases (X->Y and Y->X), thus will yield wrong results.
3. Tests were added for steps 2, 3, and 4.


ankurankan: @yzh211 Thanks for adding the tests and travis is also passing now. But I have a couple of questions about the test cases. 

1. In the first test case (https://github.com/pgmpy/pgmpy/pull/1360/files#diff-f54dc3d67c8bb46a60471463019dc1fc1e7e80f0068fbfd53ec8e8db45923304R169), the first separating set ( `frozenset({"A", "B"}): ("D",)`) doesn't seem to comply with the oriented model because if we condition the pdag on `D`, it would open the collider at `C` making `A` and `B` d-connected, which contradicts the separating set. 

2. In the second test case (https://github.com/pgmpy/pgmpy/pull/1360/files#diff-f54dc3d67c8bb46a60471463019dc1fc1e7e80f0068fbfd53ec8e8db45923304R180), I am seeing a similar issue. The separating set (`frozenset({"C", "D"}): ("B",)`) is getting violated as conditioning on `B` opens the path ( C <- A -> B <- D) making `C` and `D` d-connected, which again violates the separating set.  

The third test case seems to comply with the separating sets though. Am I not understanding something correctly here?
yzh211: @ankurankan Thanks for the questions!

1. I think you are correct on this one. I have made the corresponding change, please check if it is valid now. 

2. For the second case, there is still a block `D->B->C` as conditioning on `B`. Based on my understanding, d-separation is valid if *any*  of the three rules holds. https://en.wikipedia.org/wiki/Bayesian_network

In fact, the second test case encompasses steps 2, 3 in the algorithm, so the first testing case maybe not necessary.

Please let me know your idea?
ankurankan: @yzh211 Thanks for checking again. But I think there is still some bug in the code because in the first case (with the separating sets that you had specified before) I wouldn't expect it to return the model that it did before. For the second case, yes you are right that `D->B->C` is blocked when conditioned on `B` but it also opens the other path `C<-A->B<-D`. And for two variables to be d-separated, all the paths between them need to be blocked.

Also, I haven't looked at your code yet, so I will also give it a try and see if I can find where the bug is.
yzh211: > @yzh211 Thanks for checking again. But I think there is still some bug in the code because in the first case (with the separating sets that you had specified before) I wouldn't expect it to return the model that it did before. For the second case, yes you are right that `D->B->C` is blocked when conditioned on `B` but it also opens the other path `C<-A->B<-D`. And for two variables to be d-separated, all the paths between them need to be blocked.
> 
> Also, I haven't looked at your code yet, so I will also give it a try and see if I can find where the bug is.

Thanks!

Patch:
diff --git a/README.md b/README.md
index d1a909f9..e7477581 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@
 pgmpy
 =====
-[![Build Status](https://travis-ci.org/pgmpy/pgmpy.svg?style=flat)](https://travis-ci.org/pgmpy/pgmpy)
+[![Build Status](https://api.travis-ci.com/pgmpy/pgmpy.svg?branch=dev)](https://travis-ci.com/pgmpy/pgmpy)
 [![Appveyor](https://ci.appveyor.com/api/projects/status/github/pgmpy/pgmpy?branch=dev)](https://www.appveyor.com/)
 [![codecov](https://codecov.io/gh/pgmpy/pgmpy/branch/dev/graph/badge.svg)](https://codecov.io/gh/pgmpy/pgmpy)
 [![Codacy Badge](https://api.codacy.com/project/badge/Grade/78a8256c90654c6892627f6d8bbcea14)](https://www.codacy.com/gh/pgmpy/pgmpy?utm_source=github.com&amp;utm_medium=referral&amp;utm_content=pgmpy/pgmpy&amp;utm_campaign=Badge_Grade)
diff --git a/appveyor.yml b/appveyor.yml
index 19403b19..410767fe 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -7,7 +7,7 @@ environment:
 
 install:
   - "%PYTHON%\\python.exe -m pip install pip -U"
-  - "%PYTHON%\\python.exe -m pip install torch==1.4.0+cpu -f https://download.pytorch.org/whl/torch_stable.html"
+  - "%PYTHON%\\python.exe -m pip install torch==1.7.0+cpu -f https://download.pytorch.org/whl/torch_stable.html"
   - "%PYTHON%\\python.exe -m pip install -r requirements.txt -U"
   - "%PYTHON%\\python.exe -m pip install -e ."
 
diff --git a/pgmpy/estimators/PC.py b/pgmpy/estimators/PC.py
index f44b6c2f..259e7c63 100644
--- a/pgmpy/estimators/PC.py
+++ b/pgmpy/estimators/PC.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 
 from warnings import warn
-from itertools import combinations
+from itertools import combinations, permutations
 
 import networkx as nx
 from tqdm import tqdm
@@ -422,11 +422,12 @@ class PC(StructureEstimator):
         """
 
         pdag = skeleton.to_directed()
-        node_pairs = combinations(pdag.nodes(), 2)
+        node_pairs = list(permutations(pdag.nodes(), 2))
 
         # 1) for each X-Z-Y, if Z not in the separating set of X,Y, then orient edges as X->Z<-Y
         # (Algorithm 3.4 in Koller & Friedman PGM, page 86)
-        for X, Y in node_pairs:
+        for pair in node_pairs:
+            X, Y = pair
             if not skeleton.has_edge(X, Y):
                 for Z in set(skeleton.neighbors(X)) & set(skeleton.neighbors(Y)):
                     if Z not in separating_sets[frozenset((X, Y))]:
@@ -437,25 +438,31 @@ class PC(StructureEstimator):
             num_edges = pdag.number_of_edges()
 
             # 2) for each X->Z-Y, orient edges to Z->Y
-            for X, Y in node_pairs:
-                for Z in (set(pdag.successors(X)) - set(pdag.predecessors(X))) & (
-                    set(pdag.successors(Y)) & set(pdag.predecessors(Y))
-                ):
-                    pdag.remove(Y, Z)
+            # (Explanation in Koller & Friedman PGM, page 88)
+            for pair in node_pairs:
+                X, Y = pair
+                if not pdag.has_edge(X, Y):
+                    for Z in (set(pdag.successors(X)) - set(pdag.predecessors(X))) & (
+                        set(pdag.successors(Y)) & set(pdag.predecessors(Y))
+                    ):
+                        pdag.remove_edge(Y, Z)
 
             # 3) for each X-Y with a directed path from X to Y, orient edges to X->Y
-            for X, Y in node_pairs:
-                for path in nx.all_simple_paths(pdag, X, Y):
-                    is_directed = True
-                    for src, dst in path:
-                        if pdag.has_edge(dst, src):
-                            is_directed = False
-                    if is_directed:
-                        pdag.remove(Y, X)
-                        break
+            for pair in node_pairs:
+                X, Y = pair
+                if pdag.has_edge(Y, X) and pdag.has_edge(X, Y):
+                    for path in nx.all_simple_paths(pdag, X, Y):
+                        is_directed = True
+                        for src, dst in list(zip(path, path[1:])):
+                            if pdag.has_edge(dst, src):
+                                is_directed = False
+                        if is_directed:
+                            pdag.remove_edge(Y, X)
+                            break
 
             # 4) for each X-Z-Y with X->W, Y->W, and Z-W, orient edges to Z->W
-            for X, Y in node_pairs:
+            for pair in node_pairs:
+                X, Y = pair
                 for Z in (
                     set(pdag.successors(X))
                     & set(pdag.predecessors(X))
@@ -467,7 +474,7 @@ class PC(StructureEstimator):
                         & (set(pdag.successors(Y)) - set(pdag.predecessors(Y)))
                         & (set(pdag.successors(Z)) & set(pdag.predecessors(Z)))
                     ):
-                        pdag.remove(W, Z)
+                        pdag.remove_edge(W, Z)
 
             progress = num_edges > pdag.number_of_edges()
 


Test patch:
diff --git a/pgmpy/tests/test_estimators/test_PC.py b/pgmpy/tests/test_estimators/test_PC.py
index 0bbcded1..eae1f5f7 100644
--- a/pgmpy/tests/test_estimators/test_PC.py
+++ b/pgmpy/tests/test_estimators/test_PC.py
@@ -166,6 +166,47 @@ class TestPCEstimatorFromIndependencies(unittest.TestCase):
             set([("A", "B"), ("B", "A"), ("A", "C"), ("C", "A")]),
         )
 
+        skel = nx.Graph([("A", "C"), ("B", "C"), ("C", "D")])
+        sep_sets = {
+            frozenset({"A", "B"}): tuple(),
+            frozenset({"A", "D"}): ("C",),
+            frozenset({"B", "D"}): ("C",),
+        }
+        pdag = PC.skeleton_to_pdag(skeleton=skel, separating_sets=sep_sets)
+        self.assertSetEqual(
+            set(pdag.edges()), set([("A", "C"), ("B", "C"), ("C", "D")])
+        )
+
+        skel = nx.Graph([("A", "B"), ("A", "C"), ("B", "C"), ("B", "D")])
+        sep_sets = {
+            frozenset({"A", "D"}): tuple(),
+            frozenset({"C", "D"}): ("A", "B"),
+        }
+        pdag = PC.skeleton_to_pdag(skeleton=skel, separating_sets=sep_sets)
+        self.assertSetEqual(
+            set(pdag.edges()), set([("A", "B"), ("B", "C"), ("A", "C"), ("D", "B")])
+        )
+
+        skel = nx.Graph([("A", "B"), ("B", "C"), ("A", "D"), ("B", "D"), ("C", "D")])
+        sep_sets = {
+            frozenset({"A", "C"}): ("B",),
+        }
+        pdag = PC.skeleton_to_pdag(skeleton=skel, separating_sets=sep_sets)
+        self.assertSetEqual(
+            set(pdag.edges()),
+            set(
+                [
+                    ("A", "B"),
+                    ("B", "A"),
+                    ("B", "C"),
+                    ("C", "B"),
+                    ("A", "D"),
+                    ("B", "D"),
+                    ("C", "D"),
+                ]
+            ),
+        )
+
     def test_estimate_dag(self):
         for variant in ["orig", "stable", "parallel"]:
             ind = Independencies(["B", "C"], ["A", ["B", "C"], "D"])