code: def estimate(self, start=None, tabu_length=0, max_indegree=None):
          epsilon = 1e-8
          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):
              raise ValueError("'start' should be a DAG with the same variables as the data set, or 'None'.")
          tabu_list = []
          current_model = start
          while True:
              best_score_delta = 0
              best_operation = None
              for operation, score_delta in self._legal_operations(current_model, tabu_list, max_indegree):
                  if score_delta > best_score_delta:
                      best_operation = operation
                      best_score_delta = score_delta
              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]
              elif best_operation[0] == '-':
                  current_model.remove_edge(*best_operation[1])
                  tabu_list = ([('+', best_operation[1])] + tabu_list)[:tabu_length]
              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]
          return current_model