docs: def fit(self, X, y, headers):
          """
          Generates `num_feature_subset` sets of features each of size `size_feature_subset`.
          Thereafter, uses the input `core_algorithm` to determine feature importance scores
          for each subset. The global feature score is determined by the max score for that feature
          from all its occurences in the subsets generated. Once the final feature scores are obtained,
          the top `n_features_to_select` features can be selected.
  
          Parameters
          ----------
          X: array-like {n_samples, n_features}
              Training instances to compute the feature importance scores from
          y: array-like {n_samples}
              Training labels
          headers: array-like {n_features}
              Feature names
          Returns
          -------
          Copy of the VLSRelief instance
          """
  
          self.X_mat = X
          self._y = y
          self.headers = headers
  
          if self.core_algorithm.lower() == "multisurf":
              core = MultiSURF(n_features_to_select=self.n_features_to_select, discrete_threshold=self.discrete_threshold, verbose=self.verbose, n_jobs=self.n_jobs)
  
          elif self.core_algorithm.lower() == "multisurfstar":
              core = MultiSURFstar(n_features_to_select=self.n_features_to_select, discrete_threshold=self.discrete_threshold, verbose=self.verbose, n_jobs=self.n_jobs)
  
          elif self.core_algorithm.lower() == "surf":
              core = SURF(n_features_to_select=self.n_features_to_select, discrete_threshold=self.discrete_threshold, verbose=self.verbose, n_jobs=self.n_jobs)
  
          elif self.core_algorithm.lower() == "surfstar":
              core = SURFstar(n_features_to_select=self.n_features_to_select, discrete_threshold=self.discrete_threshold, verbose=self.verbose, n_jobs=self.n_jobs)
  
          elif self.core_algorithm.lower() == "relieff":
              core = ReliefF(n_features_to_select=self.n_features_to_select, n_neighbors=self.n_neighbors, discrete_threshold=self.discrete_threshold, verbose=self.verbose, n_jobs=self.n_jobs)
  
          total_num_features = X.shape[1]
          num_features = self.size_feature_subset
          features_scores_iter = []
          headers_iter = []
          features_selected = []
  
          for iteration in range(self.num_feature_subset):
              features_selected_id = np.random.choice(
                  range(total_num_features), num_features, replace=False)
              self.X_train = self.X_mat[:, features_selected_id]
  
              core_fit = core.fit(self.X_train, self._y)
  
              features_scores_iter.append(core_fit.feature_importances_)
              features_selected.append(features_selected_id)
              # headers_iter.append(self.headers[features_selected_id])
  
          self.features_scores_iter = features_scores_iter
          self.features_selected = features_selected
  
          zip_feat_score = [list(zip(features_selected[i], features_scores_iter[i]))
                            for i in range(len(features_selected))]
          feat_score = sorted([item for sublist in zip_feat_score for item in sublist])
          feat_score_df = pd.DataFrame(feat_score)
          feat_score_df.columns = ['feature', 'score']
          feat_score_df = feat_score_df.groupby('feature').max().reset_index()
  
          feature_scores = feat_score_df.values
  
          feature_scores = [[int(i[0]), i[1]] for i in feature_scores]
  
          self.feat_score = feature_scores
  
          head_idx = [i[0] for i in feature_scores]
          self.headers_model = list(np.array(self.headers)[head_idx])
  
          self.feature_importances_ = [i[1] for i in feature_scores]
          self.top_features_ = np.argsort(self.feature_importances_)[::-1]
          self.header_top_features_ = [self.headers_model[i] for i in self.top_features_]
  
          return self