docs: def split(self):
          """
          Split the leaf node.
          
          """
          if self.split_mode == 'random':
              # Split randomly between min and max of node's points on split dimension
              split_dim_data = self.get_data_x()[self.idxs, self.split_dim] # data on split dim
              split_min = min(split_dim_data)
              split_max = max(split_dim_data)
              split_value = split_min + np.random.rand() * (split_max - split_min)
              
          elif self.split_mode == 'median':
              # Split on median (which fall on the middle of two points for even max_points_per_region) 
              # of node's points on split dimension
              split_dim_data = self.get_data_x()[self.idxs, self.split_dim] # data on split dim
              split_value = np.median(split_dim_data)
              
          elif self.split_mode == 'middle':
              # Split on the middle of the region: might cause empty leaf
              split_dim_data = self.get_data_x()[self.idxs, self.split_dim] # data on split dim
              split_value = (self.bounds_x[0, self.split_dim] + self.bounds_x[1, self.split_dim]) / 2
              
          elif self.split_mode == 'best_interest_diff': 
              # See Baranes2012: Active Learning of Inverse Models with Intrinsically Motivated Goal Exploration in Robots
              # if strictly more than self.max_points_per_region points: chooses between self.max_points_per_region points random split values
              # the one that maximizes card(lower)*card(greater)* progress difference between the two
              # if equal or lower than self.max_points_per_region points: chooses between splits at the middle of each pair of consecutive points, 
              # the one that maximizes card(lower)*card(greater)* progress difference between the two
              split_dim_data = self.get_data_x()[self.idxs, self.split_dim] # data on split dim
              split_min = min(split_dim_data)
              split_max = max(split_dim_data)
                          
              if len(self.idxs) > self.max_points_per_region:
                  m = self.max_points_per_region # Constant that might be tuned: number of random split values to choose between
                  rand_splits = split_min + np.random.rand(m) * (split_max - split_min)
                  splits_fitness = np.zeros(m)
                  for i in range(m):
                      lower_idx = list(np.array(self.idxs)[np.nonzero(split_dim_data <= rand_splits[i])[0]])
                      greater_idx = list(np.array(self.idxs)[np.nonzero(split_dim_data > rand_splits[i])[0]])
                      splits_fitness[i] = len(lower_idx) * len(greater_idx) * abs(self.progress_idxs(lower_idx) - 
                                                                                 self.progress_idxs(greater_idx))
                  split_value = rand_splits[np.argmax(splits_fitness)]
                  
              else:
                  m = self.max_points_per_region - 1
                  splits = (np.sort(split_dim_data)[0:-1] + np.sort(split_dim_data)[1:]) / 2
                  splits_fitness = np.zeros(m)
                  for i in range(m):
                      lower_idx = list(np.array(self.idxs)[np.nonzero(split_dim_data <= splits[i])[0]])
                      greater_idx = list(np.array(self.idxs)[np.nonzero(split_dim_data > splits[i])[0]])
                      splits_fitness[i] = len(lower_idx) * len(greater_idx) * abs(self.progress_idxs(lower_idx) - 
                                                                                 self.progress_idxs(greater_idx))
                  split_value = splits[np.argmax(splits_fitness)]
              
          else:
              raise NotImplementedError
      
          lower_idx = list(np.array(self.idxs)[np.nonzero(split_dim_data <= split_value)[0]])
          greater_idx = list(np.array(self.idxs)[np.nonzero(split_dim_data > split_value)[0]])
  
          self.leafnode = False
          self.split_value = split_value
          
          split_dim = np.mod(self.split_dim + 1, np.shape(self.get_data_x())[1])
          
          l_bounds_x = np.array(self.bounds_x)
          l_bounds_x[1, self.split_dim] = split_value
          
          g_bounds_x = np.array(self.bounds_x)
          g_bounds_x[0, self.split_dim] = split_value
          
          self.lower = Tree(self.get_data_x, 
                           l_bounds_x, 
                           self.get_data_c, 
                           self.max_points_per_region, 
                           self.max_depth - 1,
                           self.split_mode, 
                           self.progress_win_size, 
                           self.progress_measure, 
                           self.sampling_mode, 
                           idxs = lower_idx, 
                           split_dim = split_dim)
          
          self.greater = Tree(self.get_data_x, 
                              g_bounds_x, 
                              self.get_data_c, 
                              self.max_points_per_region, 
                              self.max_depth - 1,
                              self.split_mode, 
                              self.progress_win_size, 
                              self.progress_measure, 
                              self.sampling_mode, 
                              idxs = greater_idx, 
                              split_dim = split_dim)