Repository: pgmpy/pgmpy

Issue:
An empty whitelist produces different results than no whitelist for HillClimbSearch
@ankurankan --- Posting as a new issue since the other issue got closed.  I pulled dev/latest and ran some tests.  Something is off with the whitelist functionality in the HillClimbSearch for Bayesian Networks.  I have a fix locally, happy to discuss first and then submit a PR.  I'm also fairly certain that this same bug makes it possible for whitelisted edges to get removed by a flip operation.

## Description
For a HillClimbSearch, specifying an empty whitelist does not reuslt in the same behavior as specifying no whitelist.  More generally, I *think* what's going on is that specifying an empty whitelist will cause HillClimbSearch to never do any flip operations.

I found this issue while attempting to show that the whitelist logic for the flip operation was still broken in dev.  Since this bug originates with the same line of code, I figured it's sufficient to post this test.

I'd be happy to submit a PR, but I'd like to discuss how the _legal_operations algorithm guarantees edges in the whitelist will be included (esp. if max_iters < len(whitelist) ) 

pgmpy 0.1.10dev
Python 3.8.0
Windows 10

## Expected Behavior
The networks should be identical whether you pass an empty whitelist or no whitelist at all.

## Actual Behavior
The networks are not identical



```
import pgmpy
import pandas as pd
from pgmpy.estimators.HillClimbSearch import HillClimbSearch
from pgmpy.estimators import K2Score
from pgmpy.models import BayesianModel
import numpy as np

X = np.asarray(["A", "B", "A", "B", "A",
    "B", "C", "C", "B", "A",
    "C", "B", "A", "C", "B",
    "D", "E", "D", "E", "D",
    "A", "B", "C", "D", "E"]).reshape((5,-1))

# synthesize data
X = pd.DataFrame({"A": X[:,0], "B": X[:,1], "C": X[:,2], "D": X[:,3], "E": X[:,4]})


# Initialize the starting model and hill climb search
model = BayesianModel([("A", "B"), ("C", "A"), ("E", "A"), ("D", "B"), ("D", "C"), ("D", "E")])
model.add_nodes_from(list(X.columns.values))
search = HillClimbSearch(data = X, scoring_method=K2Score(X))

# Kick off the hill climb search
m1 = search.estimate(start = model, max_indegree = 2, max_iter = 100, white_list = [])
print((m1.edges()))
m2 = search.estimate(start = model, max_indegree = 2, max_iter = 100)
print((m2.edges()))

assert set(m1.edges()) == set(m2.edges()) # Fails
```

Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/pgmpy/pgmpy/pull/1312?src=pr&el=h1) Report
> Merging [#1312](https://codecov.io/gh/pgmpy/pgmpy/pull/1312?src=pr&el=desc) into [dev](https://codecov.io/gh/pgmpy/pgmpy/commit/45a769eb8338eec6740781ab46e74376d680c515&el=desc) will **decrease** coverage by `0.02%`.
> The diff coverage is `92.85%`.

[![Impacted file tree graph](https://codecov.io/gh/pgmpy/pgmpy/pull/1312/graphs/tree.svg?width=650&height=150&src=pr&token=UaJMCdHaEF)](https://codecov.io/gh/pgmpy/pgmpy/pull/1312?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##              dev    #1312      +/-   ##
==========================================
- Coverage   93.08%   93.05%   -0.03%     
==========================================
  Files         132      132              
  Lines       13132    13152      +20     
==========================================
+ Hits        12224    12239      +15     
- Misses        908      913       +5     
```


| [Impacted Files](https://codecov.io/gh/pgmpy/pgmpy/pull/1312?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [pgmpy/estimators/HillClimbSearch.py](https://codecov.io/gh/pgmpy/pgmpy/pull/1312/diff?src=pr&el=tree#diff-cGdtcHkvZXN0aW1hdG9ycy9IaWxsQ2xpbWJTZWFyY2gucHk=) | `92.13% <88.63%> (-4.21%)` | :arrow_down: |
| [pgmpy/estimators/ExhaustiveSearch.py](https://codecov.io/gh/pgmpy/pgmpy/pull/1312/diff?src=pr&el=tree#diff-cGdtcHkvZXN0aW1hdG9ycy9FeGhhdXN0aXZlU2VhcmNoLnB5) | `92.10% <100.00%> (-0.21%)` | :arrow_down: |
| [...gmpy/tests/test\_estimators/test\_HillClimbSearch.py](https://codecov.io/gh/pgmpy/pgmpy/pull/1312/diff?src=pr&el=tree#diff-cGdtcHkvdGVzdHMvdGVzdF9lc3RpbWF0b3JzL3Rlc3RfSGlsbENsaW1iU2VhcmNoLnB5) | `98.57% <100.00%> (-1.43%)` | :arrow_down: |
| [pgmpy/tests/test\_estimators/test\_PC.py](https://codecov.io/gh/pgmpy/pgmpy/pull/1312/diff?src=pr&el=tree#diff-cGdtcHkvdGVzdHMvdGVzdF9lc3RpbWF0b3JzL3Rlc3RfUEMucHk=) | `95.97% <100.00%> (+0.14%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/pgmpy/pgmpy/pull/1312?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/pgmpy/pgmpy/pull/1312?src=pr&el=footer). Last update [45a769e...969688c](https://codecov.io/gh/pgmpy/pgmpy/pull/1312?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2503d0ff..e186fddc 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -9,6 +9,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 1. PC estimator with original, stable, and parallel variants.
 2. PDAG class to represent partially directed DAGs.
 3. `pgmpy.utils.get_example_model` function to fetch models from bnlearn repository.
+4. Refactor HillClimbSearch with a new feature to specify fixed edges in the model.
+5. Adds a global `SHOW_PROGRESS` variable.
+6. Adds Chow-Liu structure learning algorithm.
+7. Add `pgmpy.utils.get_example_model` to fetch models from bnlearn's repository.
 
 ### Changed
 1. Refactors ConstraintBasedEstimators into PC with a lot of general improvements.
diff --git a/pgmpy/estimators/ExhaustiveSearch.py b/pgmpy/estimators/ExhaustiveSearch.py
index 369e6113..fb303b19 100644
--- a/pgmpy/estimators/ExhaustiveSearch.py
+++ b/pgmpy/estimators/ExhaustiveSearch.py
@@ -104,9 +104,8 @@ class ExhaustiveSearch(StructureEstimator):
         all_graphs = powerset(edges)  # 2^(n*(n-1)) graphs
 
         for graph_edges in all_graphs:
-            graph = nx.DiGraph()
+            graph = nx.DiGraph(graph_edges)
             graph.add_nodes_from(nodes)
-            graph.add_edges_from(graph_edges)
             if nx.is_directed_acyclic_graph(graph):
                 yield graph
 
diff --git a/pgmpy/estimators/HillClimbSearch.py b/pgmpy/estimators/HillClimbSearch.py
index 4235eb39..46ae82a5 100644
--- a/pgmpy/estimators/HillClimbSearch.py
+++ b/pgmpy/estimators/HillClimbSearch.py
@@ -1,10 +1,18 @@
 #!/usr/bin/env python
 from itertools import permutations
+from collections import deque
 
 import networkx as nx
 from tqdm import trange
 
-from pgmpy.estimators import StructureEstimator, K2Score, ScoreCache
+from pgmpy.estimators import (
+    StructureScore,
+    StructureEstimator,
+    K2Score,
+    ScoreCache,
+    BDeuScore,
+    BicScore,
+)
 from pgmpy.base import DAG
 from pgmpy.global_vars import SHOW_PROGRESS
 
@@ -22,10 +30,6 @@ class HillClimbSearch(StructureEstimator):
             (If some values in the data are missing the data cells should be set to `numpy.NaN`.
             Note that pandas converts each column containing `numpy.NaN`s to dtype `float`.)
 
-        scoring_method: Instance of a `StructureScore`-subclass (`K2Score` is used as default)
-            An instance of `K2Score`, `BDeuScore`, or `BicScore`.
-            This score is optimized during structure estimation by the `estimate`-method.
-
         state_names: dict (optional)
             A dict indicating, for each variable, the discrete set of states (or values)
             that the variable can take. If unspecified, the observed values in the data set
@@ -47,18 +51,12 @@ class HillClimbSearch(StructureEstimator):
         Koller & Friedman, Probabilistic Graphical Models - Principles and Techniques, 2009
         Section 18.4.3 (page 811ff)
         """
-        if scoring_method is not None:
-            if use_cache:
-                self.scoring_method = ScoreCache.ScoreCache(scoring_method, data)
-            else:
-                self.scoring_method = scoring_method
-        else:
-            self.scoring_method = ScoreCache.ScoreCache(K2Score(data, **kwargs), data)
+        self.use_cache = use_cache
 
         super(HillClimbSearch, self).__init__(data, **kwargs)
 
     def _legal_operations(
-        self, model, tabu_list=[], max_indegree=None, black_list=None, white_list=None
+        self, model, score, tabu_list, max_indegree, black_list, white_list, fixed_edges
     ):
         """Generates a list of legal (= not in tabu_list) graph modifications
         for a given model, together with their score changes. Possible graph modifications:
@@ -70,68 +68,73 @@ class HillClimbSearch(StructureEstimator):
         edges or to limit the search.
         """
 
-        local_score = self.scoring_method.local_score
-        nodes = self.state_names.keys()
+        tabu_list = set(tabu_list)
+
+        # Step 1: Get all legal operations for adding edges.
         potential_new_edges = (
-            set(permutations(nodes, 2))
+            set(permutations(self.variables, 2))
             - set(model.edges())
             - set([(Y, X) for (X, Y) in model.edges()])
         )
 
-        for (X, Y) in potential_new_edges:  # (1) add single edge
-            if nx.is_directed_acyclic_graph(nx.DiGraph(list(model.edges()) + [(X, Y)])):
+        for (X, Y) in potential_new_edges:
+            # Check if adding (X, Y) will create a cycle.
+            if not nx.has_path(model, Y, X):
                 operation = ("+", (X, Y))
                 if (
-                    operation not in tabu_list
-                    and (black_list is None or (X, Y) not in black_list)
-                    and (white_list is None or (X, Y) in white_list)
+                    (operation not in tabu_list)
+                    and ((X, Y) not in black_list)
+                    and ((X, Y) in white_list)
                 ):
                     old_parents = model.get_parents(Y)
                     new_parents = old_parents + [X]
-                    if max_indegree is None or len(new_parents) <= max_indegree:
-                        score_delta = local_score(Y, new_parents) - local_score(
-                            Y, old_parents
-                        )
+                    if len(new_parents) <= max_indegree:
+                        score_delta = score(Y, new_parents) - score(Y, old_parents)
                         yield (operation, score_delta)
 
-        for (X, Y) in model.edges():  # (2) remove single edge
+        # Step 2: Get all legal operations for removing edges
+        for (X, Y) in model.edges():
             operation = ("-", (X, Y))
-            if operation not in tabu_list:
+            if (operation not in tabu_list) and ((X, Y) not in fixed_edges):
                 old_parents = model.get_parents(Y)
                 new_parents = old_parents[:]
                 new_parents.remove(X)
-                score_delta = local_score(Y, new_parents) - local_score(Y, old_parents)
+                score_delta = score(Y, new_parents) - score(Y, old_parents)
                 yield (operation, score_delta)
 
-        for (X, Y) in model.edges():  # (3) flip single edge
-            new_edges = list(model.edges()) + [(Y, X)]
-            new_edges.remove((X, Y))
-            if nx.is_directed_acyclic_graph(nx.DiGraph(new_edges)):
+        # Step 3: Get all legal operations for flipping edges
+        for (X, Y) in model.edges():
+            # Check if flipping creates any cycles
+            if not any(
+                map(lambda path: len(path) > 2, nx.all_simple_paths(model, X, Y))
+            ):
                 operation = ("flip", (X, Y))
                 if (
-                    operation not in tabu_list
-                    and ("flip", (Y, X)) not in tabu_list
-                    and (black_list is None or (Y, X) not in black_list)
-                    and (white_list is None or (Y, X) in white_list)
+                    ((operation not in tabu_list) and ("flip", (Y, X)) not in tabu_list)
+                    and ((X, Y) not in fixed_edges)
+                    and ((Y, X) not in black_list)
+                    and ((Y, X) in white_list)
                 ):
                     old_X_parents = model.get_parents(X)
                     old_Y_parents = model.get_parents(Y)
                     new_X_parents = old_X_parents + [Y]
                     new_Y_parents = old_Y_parents[:]
                     new_Y_parents.remove(X)
-                    if max_indegree is None or len(new_X_parents) <= max_indegree:
+                    if len(new_X_parents) <= max_indegree:
                         score_delta = (
-                            local_score(X, new_X_parents)
-                            + local_score(Y, new_Y_parents)
-                            - local_score(X, old_X_parents)
-                            - local_score(Y, old_Y_parents)
+                            score(X, new_X_parents)
+                            + score(Y, new_Y_parents)
+                            - score(X, old_X_parents)
+                            - score(Y, old_Y_parents)
                         )
                         yield (operation, score_delta)
 
     def estimate(
         self,
-        start=None,
-        tabu_length=0,
+        scoring_method="k2score",
+        start_dag=None,
+        fixed_edges=set(),
+        tabu_length=100,
         max_indegree=None,
         black_list=None,
         white_list=None,
@@ -140,15 +143,27 @@ class HillClimbSearch(StructureEstimator):
         show_progress=True,
     ):
         """
-        Performs local hill climb search to estimates the `DAG` structure
-        that has optimal score, according to the scoring method supplied in the constructor.
-        Starts at model `start` and proceeds by step-by-step network modifications
-        until a local maximum is reached. Only estimates network structure, no parametrization.
+        Performs local hill climb search to estimates the `DAG` structure that
+        has optimal score, according to the scoring method supplied. Starts at
+        model `start_dag` and proceeds by step-by-step network modifications
+        until a local maximum is reached. Only estimates network structure, no
+        parametrization.
 
         Parameters
         ----------
-        start: DAG instance
-            The starting point for the local search. By default a completely disconnected network is used.
+        scoring_method: str or StructureScore instance
+            The score to be optimized during structure estimation.  Supported
+            structure scores: k2score, bdeuscore, bicscore. Also accepts a
+            custom score but it should be an instance of `StructureScore`.
+
+        start_dag: DAG instance
+            The starting point for the local search. By default a completely
+            disconnected network is used.
+
+        fixed_edges: iterable
+            A list of edges that will always be there in the final learned model.
+            The algorithm will add these edges at the start of the algorithm and
+            will never change it.
 
         tabu_length: int
             If provided, the last `tabu_length` graph modifications cannot be reversed
@@ -198,46 +213,102 @@ class HillClimbSearch(StructureEstimator):
         >>> est.estimate(max_indegree=1).edges()
         [('J', 'A'), ('B', 'J')]
         """
-        nodes = self.state_names.keys()
-        if start is None:
-            start = DAG()
-            start.add_nodes_from(nodes)
-        elif not isinstance(start, DAG) or not set(start.nodes()) == set(nodes):
+
+        # Step 1: Initial checks and setup for arguments
+        # Step 1.1: Check scoring_method
+        supported_methods = {
+            "k2score": K2Score,
+            "bdeuscore": BDeuScore,
+            "bicscore": BicScore,
+        }
+        if (scoring_method.lower() not in supported_methods) and (
+            not isinstance(scoring_method, StructureScore)
+        ):
             raise ValueError(
-                "'start' should be a DAG with the same variables as the data set, or 'None'."
+                "scoring_method should either be one of k2score, bdeuscore, bicscore, or an instance of StructureScore"
             )
 
-        tabu_list = []
-        current_model = start
+        if isinstance(scoring_method, str):
+            score = supported_methods[scoring_method.lower()](data=self.data)
+        else:
+            score = scoring_method(data=self.data)
+
+        if self.use_cache:
+            score_fn = ScoreCache.ScoreCache(score, self.data).local_score
+        else:
+            score_fn = score.local_score
+
+        # Step 1.2: Check the start_dag
+        if start_dag is None:
+            start_dag = DAG()
+            start_dag.add_nodes_from(self.variables)
+        elif not isinstance(start_dag, DAG) or not set(start_dag.nodes()) == set(
+            self.variables
+        ):
+            raise ValueError(
+                "'start_dag' should be a DAG with the same variables as the data set, or 'None'."
+            )
+
+        # Step 1.3: Check fixed_edges
+        if not hasattr(fixed_edges, "__iter__"):
+            raise ValueError("fixed_edges must be an iterable")
+        else:
+            fixed_edges = set(fixed_edges)
+            start_dag.add_edges_from(fixed_edges)
+            if not nx.is_directed_acyclic_graph(start_dag):
+                raise ValueError(
+                    "fixed_edges creates a cycle in start_dag. Please modify either fixed_edges or start_dag."
+                )
+
+        # Step 1.4: Check black list and white list
+        black_list = set() if black_list is None else set(black_list)
+        white_list = (
+            set([(u, v) for u in self.variables for v in self.variables])
+            if white_list is None
+            else set(white_list)
+        )
+
+        # Step 1.5: Initialize max_indegree, tabu_list, and progress bar
+        if max_indegree is None:
+            max_indegree = float("inf")
+
+        tabu_list = deque(maxlen=tabu_length)
+        current_model = start_dag
 
         if show_progress and SHOW_PROGRESS:
             iteration = trange(int(max_iter))
         else:
             iteration = range(int(max_iter))
 
-        for iter_no in iteration:
-            best_score_delta = 0
-            best_operation = None
-
-            for operation, score_delta in self._legal_operations(
-                current_model, tabu_list, max_indegree, black_list, white_list
-            ):
-                if score_delta > best_score_delta:
-                    best_operation = operation
-                    best_score_delta = score_delta
+        # Step 2: For each iteration, find the best scoring operation and
+        #         do that to the current model.
+        for _ in iteration:
+            best_operation, best_score_delta = max(
+                self._legal_operations(
+                    current_model,
+                    score_fn,
+                    tabu_list,
+                    max_indegree,
+                    black_list,
+                    white_list,
+                    fixed_edges,
+                ),
+                key=lambda t: t[1],
+            )
 
             if best_operation is None or best_score_delta < epsilon:
                 break
             elif best_operation[0] == "+":
                 current_model.add_edge(*best_operation[1])
-                tabu_list = ([("-", best_operation[1])] + tabu_list)[:tabu_length]
+                tabu_list.append(("-", best_operation[1]))
             elif best_operation[0] == "-":
                 current_model.remove_edge(*best_operation[1])
-                tabu_list = ([("+", best_operation[1])] + tabu_list)[:tabu_length]
+                tabu_list.append(("+", best_operation[1]))
             elif best_operation[0] == "flip":
                 X, Y = best_operation[1]
                 current_model.remove_edge(X, Y)
                 current_model.add_edge(Y, X)
-                tabu_list = ([best_operation] + tabu_list)[:tabu_length]
+                tabu_list.append(best_operation)
 
+        # Step 3: Return if no more improvements or maximum iterations reached.
         return current_model


Test patch:
diff --git a/pgmpy/tests/test_estimators/test_HillClimbSearch.py b/pgmpy/tests/test_estimators/test_HillClimbSearch.py
index f3313f73..a1c147b9 100644
--- a/pgmpy/tests/test_estimators/test_HillClimbSearch.py
+++ b/pgmpy/tests/test_estimators/test_HillClimbSearch.py
@@ -10,16 +10,23 @@ from pgmpy.models import BayesianModel
 class TestHillClimbEstimator(unittest.TestCase):
     def setUp(self):
         self.rand_data = pd.DataFrame(
-            np.random.randint(0, 5, size=(5000, 2)), columns=list("AB")
+            np.random.randint(0, 5, size=(int(1e4), 2)), columns=list("AB")
         )
         self.rand_data["C"] = self.rand_data["B"]
-        self.est_rand = HillClimbSearch(
-            self.rand_data, scoring_method=K2Score(self.rand_data)
-        )
+        self.est_rand = HillClimbSearch(self.rand_data)
+        self.score_rand = K2Score(self.rand_data).local_score
+
         self.model1 = BayesianModel()
         self.model1.add_nodes_from(["A", "B", "C"])
+        self.model1_possible_edges = set(
+            [(u, v) for u in self.model1.nodes() for v in self.model1.nodes()]
+        )
+
         self.model2 = self.model1.copy()
         self.model2.add_edge("A", "B")
+        self.model2_possible_edges = set(
+            [(u, v) for u in self.model2.nodes() for v in self.model2.nodes()]
+        )
 
         # link to dataset: "https://www.kaggle.com/c/titanic/download/train.csv"
         self.titanic_data = pd.read_csv(
@@ -28,12 +35,25 @@ class TestHillClimbEstimator(unittest.TestCase):
         self.titanic_data1 = self.titanic_data[
             ["Survived", "Sex", "Pclass", "Age", "Embarked"]
         ]
-        self.titanic_data2 = self.titanic_data[["Survived", "Sex", "Pclass"]]
         self.est_titanic1 = HillClimbSearch(self.titanic_data1)
+        self.score_titanic1 = K2Score(self.titanic_data1).local_score
+
+        self.titanic_data2 = self.titanic_data[["Survived", "Sex", "Pclass"]]
         self.est_titanic2 = HillClimbSearch(self.titanic_data2)
+        self.score_titanic2 = K2Score(self.titanic_data2).local_score
 
     def test_legal_operations(self):
-        model2_legal_ops = list(self.est_rand._legal_operations(self.model2))
+        model2_legal_ops = list(
+            self.est_rand._legal_operations(
+                model=self.model2,
+                score=self.score_rand,
+                tabu_list=set(),
+                max_indegree=float("inf"),
+                black_list=set(),
+                white_list=self.model2_possible_edges,
+                fixed_edges=set(),
+            )
+        )
         model2_legal_ops_ref = [
             (("+", ("C", "A")), -28.15602208305154),
             (("+", ("A", "C")), -28.155467430966382),
@@ -50,7 +70,13 @@ class TestHillClimbEstimator(unittest.TestCase):
     def test_legal_operations_blacklist_whitelist(self):
         model2_legal_ops_bl = list(
             self.est_rand._legal_operations(
-                self.model2, black_list=[("A", "B"), ("A", "C"), ("C", "A"), ("C", "B")]
+                model=self.model2,
+                score=self.score_rand,
+                tabu_list=set(),
+                max_indegree=float("inf"),
+                black_list=set([("A", "B"), ("A", "C"), ("C", "A"), ("C", "B")]),
+                white_list=self.model2_possible_edges,
+                fixed_edges=set(),
             )
         )
         model2_legal_ops_bl_ref = [
@@ -64,7 +90,13 @@ class TestHillClimbEstimator(unittest.TestCase):
 
         model2_legal_ops_wl = list(
             self.est_rand._legal_operations(
-                self.model2, white_list=[("A", "B"), ("A", "C"), ("C", "A"), ("A", "B")]
+                model=self.model2,
+                score=self.score_rand,
+                tabu_list=set(),
+                max_indegree=float("inf"),
+                black_list=set(),
+                white_list=set([("A", "B"), ("A", "C"), ("C", "A"), ("A", "B")]),
+                fixed_edges=set(),
             )
         )
         model2_legal_ops_wl_ref = [
@@ -77,12 +109,21 @@ class TestHillClimbEstimator(unittest.TestCase):
         )
 
     def test_legal_operations_titanic(self):
-        est = self.est_titanic1
         start_model = BayesianModel(
             [("Survived", "Sex"), ("Pclass", "Age"), ("Pclass", "Embarked")]
         )
-
-        legal_ops = est._legal_operations(start_model)
+        all_possible_edges = set(
+            [(u, v) for u in start_model.nodes() for v in start_model.nodes()]
+        )
+        legal_ops = self.est_titanic1._legal_operations(
+            model=start_model,
+            score=self.score_titanic1,
+            tabu_list=[],
+            max_indegree=float("inf"),
+            black_list=set(),
+            white_list=all_possible_edges,
+            fixed_edges=set(),
+        )
         self.assertEqual(len(list(legal_ops)), 20)
 
         tabu_list = [
@@ -90,15 +131,38 @@ class TestHillClimbEstimator(unittest.TestCase):
             ("-", ("Survived", "Pclass")),
             ("flip", ("Age", "Pclass")),
         ]
-        legal_ops_tabu = est._legal_operations(start_model, tabu_list=tabu_list)
+        legal_ops_tabu = self.est_titanic1._legal_operations(
+            model=start_model,
+            score=self.score_titanic1,
+            tabu_list=tabu_list,
+            max_indegree=float("inf"),
+            black_list=set(),
+            white_list=all_possible_edges,
+            fixed_edges=set(),
+        )
         self.assertEqual(len(list(legal_ops_tabu)), 18)
 
-        legal_ops_indegree = est._legal_operations(start_model, max_indegree=1)
+        legal_ops_indegree = self.est_titanic1._legal_operations(
+            model=start_model,
+            score=self.score_titanic1,
+            tabu_list=[],
+            max_indegree=1,
+            black_list=set(),
+            white_list=all_possible_edges,
+            fixed_edges=set(),
+        )
         self.assertEqual(len(list(legal_ops_indegree)), 11)
 
-        legal_ops_both = list(
-            est._legal_operations(start_model, tabu_list=tabu_list, max_indegree=1)
+        legal_ops_both = self.est_titanic1._legal_operations(
+            model=start_model,
+            score=self.score_titanic1,
+            tabu_list=tabu_list,
+            max_indegree=1,
+            black_list=set(),
+            white_list=all_possible_edges,
+            fixed_edges=set(),
         )
+
         legal_ops_both_ref = {
             ("+", ("Embarked", "Survived")): 10.050632580087495,
             ("+", ("Survived", "Pclass")): 41.8886804654893,
@@ -123,17 +187,25 @@ class TestHillClimbEstimator(unittest.TestCase):
             list(est1.edges()) == [("B", "C")] or list(est1.edges()) == [("C", "B")]
         )
 
-        est2 = self.est_rand.estimate(start=BayesianModel([("A", "B"), ("A", "C")]))
+        est2 = self.est_rand.estimate(start_dag=BayesianModel([("A", "B"), ("A", "C")]))
         self.assertTrue(
             list(est2.edges()) == [("B", "C")] or list(est2.edges()) == [("C", "B")]
         )
 
+        est3 = self.est_rand.estimate(fixed_edges=[("B", "C")])
+        self.assertTrue([("B", "C")] == list(est3.edges()))
+
     def test_estimate_titanic(self):
         self.assertSetEqual(
             set(self.est_titanic2.estimate().edges()),
             set([("Survived", "Pclass"), ("Sex", "Pclass"), ("Sex", "Survived")]),
         )
 
+        self.assertTrue(
+            ("Pclass", "Survived")
+            in self.est_titanic2.estimate(fixed_edges=[("Pclass", "Survived")]).edges()
+        )
+
     def tearDown(self):
         del self.rand_data
         del self.est_rand
diff --git a/pgmpy/tests/test_estimators/test_PC.py b/pgmpy/tests/test_estimators/test_PC.py
index fac68722..84cd5ab6 100644
--- a/pgmpy/tests/test_estimators/test_PC.py
+++ b/pgmpy/tests/test_estimators/test_PC.py
@@ -8,13 +8,12 @@ from pgmpy.estimators import PC
 from pgmpy.models import BayesianModel
 from pgmpy.independencies import Independencies
 
-np.random.seed(42)
-
 # This class tests examples from: Le, Thuc, et al. "A fast PC algorithm for
 # high dimensional causal discovery with multi-core PCs." IEEE/ACM transactions
 # on computational biology and bioinformatics (2016).
 class TestPCFakeCITest(unittest.TestCase):
     def setUp(self):
+        np.random.seed(42)
         self.fake_data = pd.DataFrame(
             np.random.random((1000, 4)), columns=["A", "B", "C", "D"]
         )
@@ -193,13 +192,17 @@ class TestPCEstimatorFromDiscreteData(unittest.TestCase):
     def test_build_skeleton(self):
         for variant in ["orig", "stable", "parallel"]:
             # Fake dataset no: 1
+            np.random.seed(42)
             data = pd.DataFrame(
                 np.random.randint(0, 2, size=(10000, 5)), columns=list("ABCDE")
             )
             data["F"] = data["A"] + data["B"] + data["C"]
             est = PC(data=data)
             skel, sep_sets = est.estimate(
-                variant=variant, ci_test="chi_square", return_type="skeleton"
+                variant=variant,
+                ci_test="chi_square",
+                return_type="skeleton",
+                significance_level=0.005,
             )
             expected_edges = {("A", "F"), ("B", "F"), ("C", "F")}
             expected_sepsets = {
@@ -231,6 +234,7 @@ class TestPCEstimatorFromDiscreteData(unittest.TestCase):
                 else:
                     return False
 
+            np.random.seed(42)
             fake_data = pd.DataFrame(
                 np.random.randint(low=0, high=2, size=(10000, 3)),
                 columns=["X", "Y", "Z"],
@@ -249,13 +253,18 @@ class TestPCEstimatorFromDiscreteData(unittest.TestCase):
 
     def test_build_dag(self):
         for variant in ["orig", "stable", "parallel"]:
+            np.random.seed(42)
             data = pd.DataFrame(
                 np.random.randint(0, 3, size=(10000, 3)), columns=list("XYZ")
             )
             data["sum"] = data.sum(axis=1)
             est = PC(data=data)
-            dag = est.estimate(variant=variant, ci_test="chi_square", return_type="dag")
-
+            dag = est.estimate(
+                variant=variant,
+                ci_test="chi_square",
+                return_type="dag",
+                significance_level=0.001,
+            )
             expected_edges = {("Z", "sum"), ("X", "sum"), ("Y", "sum")}
             self.assertEqual(set(dag.edges()), expected_edges)
 
@@ -264,6 +273,7 @@ class TestPCEstimatorFromContinuousData(unittest.TestCase):
     def test_build_skeleton(self):
         for variant in ["orig", "stable", "parallel"]:
             # Fake dataset no: 1
+            np.random.seed(42)
             data = pd.DataFrame(np.random.randn(10000, 5), columns=list("ABCDE"))
             data["F"] = data["A"] + data["B"] + data["C"]
             est = PC(data=data)
@@ -306,6 +316,7 @@ class TestPCEstimatorFromContinuousData(unittest.TestCase):
                 else:
                     return False
 
+            np.random.seed(42)
             data = pd.DataFrame(np.random.randn(10000, 3), columns=list("XYZ"))
             est = PC(data=data)
             skel, sep_sets = est.estimate(
@@ -322,6 +333,7 @@ class TestPCEstimatorFromContinuousData(unittest.TestCase):
 
     def test_build_dag(self):
         for variant in ["orig", "stable", "parallel"]:
+            np.random.seed(42)
             data = pd.DataFrame(np.random.randn(10000, 3), columns=list("XYZ"))
             data["sum"] = data.sum(axis=1)
             est = PC(data=data)