docs: def _best_split(self, y_true, X, cost_mat):
          """ Private function to calculate the split that gives the best gain.
  
          Parameters
          ----------
          y_true : array indicator matrix
              Ground truth (correct) labels.
  
          X : array-like of shape = [n_samples, n_features]
              The input samples.
  
          cost_mat : array-like of shape = [n_samples, 4]
              Cost matrix of the classification problem
              Where the columns represents the costs of: false positives, false negatives,
              true positives and true negatives, for each example.
  
          Returns
          -------
          tuple(split : tuple(j, l), gain : float, left node prediction : int,
                y_pred : int, y_prob : float)
  
          """
  
          n_samples, n_features = X.shape
          num_pct = self.num_pct
  
          cost_base, y_pred, y_prob = self._node_cost(y_true, cost_mat)
  
          # Calculate the gain of all features each split in num_pct
          gains = np.zeros((n_features, num_pct))
          pred = np.zeros((n_features, num_pct))
          splits = np.zeros((n_features, num_pct))
  
          # Selected features
          selected_features = np.arange(0, self.n_features_)
          # Add random state
          np.random.shuffle(selected_features)
          selected_features = selected_features[:self.max_features_]
          selected_features.sort()
  
          #TODO:  # Skip the CPU intensive evaluation of the impurity criterion for
                  # features that were already detected as constant (hence not suitable
                  # for good splitting) by ancestor nodes and save the information on
                  # newly discovered constant features to spare computation on descendant
                  # nodes.
  
          # For each feature test all possible splits
          for j in selected_features:
              splits[j, :] = np.percentile(X[:, j], np.arange(0, 100, 100.0 / num_pct).tolist())
  
              for l in range(num_pct):
                  # Avoid repeated values, since np.percentile may return repeated values
                  if l == 0 or (l > 0 and splits[j, l] != splits[j, l - 1]):
                      split = (j, splits[j, l])
                      gains[j, l], pred[j, l] = self._calculate_gain(cost_base, y_true, X, cost_mat, split)
  
          best_split = np.unravel_index(gains.argmax(), gains.shape)
  
          return (best_split[0], splits[best_split]), gains.max(), pred[best_split], y_pred, y_prob