docs: def estimate(self, start=None, tabu_length=0, max_indegree=None):
          """
          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.
  
          Parameters
          ----------
          start: DAG instance
              The starting point for the local search. By default a completely disconnected network is used.
          tabu_length: int
              If provided, the last `tabu_length` graph modifications cannot be reversed
              during the search procedure. This serves to enforce a wider exploration
              of the search space. Default value: 100.
          max_indegree: int or None
              If provided and unequal None, the procedure only searches among models
              where all nodes have at most `max_indegree` parents. Defaults to None.
  
          Returns
          -------
          model: `DAG` instance
              A `DAG` at a (local) score maximum.
  
          Examples
          --------
          >>> import pandas as pd
          >>> import numpy as np
          >>> from pgmpy.estimators import HillClimbSearch, BicScore
          >>> # create data sample with 9 random variables:
          ... data = pd.DataFrame(np.random.randint(0, 5, size=(5000, 9)), columns=list('ABCDEFGHI'))
          >>> # add 10th dependent variable
          ... data['J'] = data['A'] * data['B']
          >>> est = HillClimbSearch(data, scoring_method=BicScore(data))
          >>> best_model = est.estimate()
          >>> sorted(best_model.nodes())
          ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J']
          >>> best_model.edges()
          [('B', 'J'), ('A', 'J')]
          >>> # search a model with restriction on the number of parents:
          >>> est.estimate(max_indegree=1).edges()
          [('J', 'A'), ('B', 'J')]
          """
          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