docs: def skeleton_to_pdag(skel, separating_sets):
          """Orients the edges of a graph skeleton based on information from
          `separating_sets` to form a DAG pattern (DAG).
  
          Parameters
          ----------
          skel: UndirectedGraph
              An undirected graph skeleton as e.g. produced by the
              estimate_skeleton method.
  
          separating_sets: dict
              A dict containing for each pair of not directly connected nodes a
              separating set ("witnessing set") of variables that makes then
              conditionally independent. (needed for edge orientation)
  
          Returns
          -------
          pdag: DAG
              An estimate for the DAG pattern of the BN underlying the data. The
              graph might contain some nodes with both-way edges (X->Y and Y->X).
              Any completion by (removing one of the both-way edges for each such
              pair) results in a I-equivalent Bayesian network DAG.
  
          Reference
          ---------
          Neapolitan, Learning Bayesian Networks, Section 10.1.2, Algorithm 10.2 (page 550)
          http://www.cs.technion.ac.il/~dang/books/Learning%20Bayesian%20Networks(Neapolitan,%20Richard).pdf
  
  
          Examples
          --------
          >>> import pandas as pd
          >>> import numpy as np
          >>> from pgmpy.estimators import ConstraintBasedEstimator
          >>> data = pd.DataFrame(np.random.randint(0, 4, size=(5000, 3)), columns=list('ABD'))
          >>> data['C'] = data['A'] - data['B']
          >>> data['D'] += data['A']
          >>> c = ConstraintBasedEstimator(data)
          >>> pdag = c.skeleton_to_pdag(*c.estimate_skeleton())
          >>> pdag.edges() # edges: A->C, B->C, A--D (not directed)
          [('B', 'C'), ('A', 'C'), ('A', 'D'), ('D', 'A')]
          """
  
          pdag = skel.to_directed()
          node_pairs = combinations(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:
              if not skel.has_edge(X, Y):
                  for Z in set(skel.neighbors(X)) & set(skel.neighbors(Y)):
                      if Z not in separating_sets[frozenset((X, Y))]:
                          pdag.remove_edges_from([(Z, X), (Z, Y)])
  
          progress = True
          while progress:  # as long as edges can be oriented (removed)
              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)
  
              # 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
  
              # 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 Z in (set(pdag.successors(X)) & set(pdag.predecessors(X)) &
                            set(pdag.successors(Y)) & set(pdag.predecessors(Y))):
                      for W in ((set(pdag.successors(X)) - set(pdag.predecessors(X))) &
                                (set(pdag.successors(Y)) - set(pdag.predecessors(Y))) &
                                (set(pdag.successors(Z)) & set(pdag.predecessors(Z)))):
                          pdag.remove(W, Z)
  
              progress = num_edges > pdag.number_of_edges()
  
          return pdag