docs: def _legal_operations(self, model, tabu_list=[], max_indegree=None):
          """Generates a list of legal (= not in tabu_list) graph modifications
          for a given model, together with their score changes. Possible graph modifications:
          (1) add, (2) remove, or (3) flip a single edge. For details on scoring
          see Koller & Fridman, Probabilistic Graphical Models, Section 18.4.3.3 (page 818).
          If a number `max_indegree` is provided, only modifications that keep the number
          of parents for each node below `max_indegree` are considered."""
  
          local_score = self.scoring_method.local_score
          nodes = self.state_names.keys()
          potential_new_edges = (set(permutations(nodes, 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)])):
                  operation = ('+', (X, Y))
                  if operation not in tabu_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)
                          yield(operation, score_delta)
  
          for (X, Y) in model.edges():  # (2) remove single edge
              operation = ('-', (X, Y))
              if operation not in tabu_list:
                  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)
                  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)):
                  operation = ('flip', (X, Y))
                  if operation not in tabu_list and ('flip', (Y, X)) not in tabu_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:
                          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))
                          yield(operation, score_delta)