docs: def fit(self, features, labels):
          """Fit model to data"""
  
          # setup data
          # imputation
          if self.clean:
              features = self.impute_data(features)
          # save the number of features
          self.n_features = features.shape[1]
          self.n_samples = features.shape[0]
          # set population size
          if type(self.population_size) is str:
              if 'x' in self.population_size: # set pop size prop to features
                  self.population_size = int(
                          float(self.population_size[:-1])*features.shape[1])
              else:
                  self.population_size = int(self.population_size)
  
          if self.verbosity >0: print("population size:",self.population_size)
          
  
          # re-initialize pipeline (needs to be here rather than init for GridSearchCV)
          if self.normalize:
              self.pipeline = Pipeline([('standardScaler',StandardScaler()), ('ml', self.ml)])
          else:
              self.pipeline = Pipeline([('ml',self.ml)])
  
          # set variable names if they haven't been set
          if self.names is None:
              self.names = ['x_'+str(i) for i in np.arange(features.shape[1])]
          # set variable data types if they haven't been set
          if self.dtypes is None:
              self.dtypes = ['f' for i in np.arange(features.shape[1])]
          
          # create terminal set
          for i in np.arange(self.n_features):
              self.term_set.append(node('x',loc=i,otype=self.dtypes[i])) # features
              # add ephemeral random constants if flag
              if self.erc: # ephemeral random constants
                  self.term_set.append(node('k',value=self.random_state.rand()))
  
          # edit function set if boolean
          if self.boolean or self.otype=='b': # include boolean functions
              self.func_set += [node('!'), node('&'), node('|'), node('=='),
                          node('>_f'), node('<_f'), node('>=_f'), node('<=_f'),
                          node('>_b'), node('<_b'), node('>=_b'), node('<=_b'),
                          node('xor_b'), node('xor_f')]
  
          # add mdr if specified
          if self.mdr:
              self.func_set += [node('mdr2')]
  
          # print few settings
          if self.verbosity > 1:
              for arg in self.get_params():
                  print('{}\t=\t{}'.format(arg, self.get_params()[arg]))
              print('')
          
          ######################################################### initial model
          # fit to original data
          with warnings.catch_warnings():
              warnings.simplefilter("ignore")
              if self.scoring_function == roc_auc_score:
                  self._best_score = self.roc_auc_cv(features,labels)
              else:
                  self._best_score = np.mean(
                                     [self.scoring_function(labels[test],
                                         self.pipeline.fit(features[train],labels[train]).
                                                       predict(features[test]))
                                     for train, test in KFold().split(features,
                                                                       labels)])
  
          initial_score = self._best_score
          if self.verbosity > 0:
              print("initial ML CV: {:1.3f}".format(self._best_score))
  
          
          ############################################# Create initial population
          # for now, force seed_with_ml to be off if otype is 'b', since data
          # types are assumed to be float
          if self.otype=='b':
              self.seed_with_ml = False
          self.pop = self.init_pop()
          # check that uuids are unique in population
          uuids = [p.id for p in self.pop.individuals]
          if len(uuids) != len(set(uuids)):
              pdb.set_trace()
          # Evaluate the entire population
          # X represents a matrix of the population outputs (number of samples x
          # population size)
          # single thread
          self.X = self.transform(features,self.pop.individuals,labels).transpose()
          # pdb.set_trace()
          # parallel:
          # X = np.asarray(Parallel(n_jobs=-1)(
          # delayed(out)(I,features,self.otype,labels) for I in self.pop.individuals),
          # order = 'F')
  
          # calculate fitness of individuals
          # fitnesses = list(map(lambda I: fitness(I,labels,self.pipeline),X))
          self.F = self.calc_fitness(self.X,labels,self.fit_choice,self.sel)
          #pdb.set_trace()
          #with Parallel(n_jobs=10) as parallel:
          ####################
  
          self.diversity=[]
          # progress bar
          pbar = tqdm(total=self.generations,disable = self.verbosity==0,
                      desc='Internal CV: {:1.3f}'.format(self._best_score))
          stall_count = 0
          ########################################################### main GP loop
          for g in np.arange(self.generations):
              if stall_count == self.max_stall:
                  if self.verbosity > 0: print('max stall count reached.')
                  break
  
              if self.track_diversity:
                  self.get_diversity(self.X)
  
              # mid verbosity printouts
              if self.verbosity > 1:
                  
                  print("generation", str(g))          
              
                  print("median fitness pop: %0.2f" % np.median(
                          [np.mean(f) for f in self.F]))
              
                  print("best fitness pop: %0.2f" % np.min(
                      [np.mean(f) for f in self.F]))
                  if self.track_diversity:
                      print("feature diversity: %0.2f" % self.diversity[-1])
  
              # high verbosity printouts    
              if self.verbosity > 2:
                  eqns = self.stacks_2_eqns(self.pop.individuals)
                  fs = [np.mean(f) for f in self.F] 
                  print("population:",[("%0.2f" % f, eqns[i]) for f,i in 
                                       zip(np.sort(fs), np.argsort(fs))])
                  #print("pop fitnesses:", ["%0.2f" % np.mean(f) for f in self.F])
              
              ####################################################### fit ml model
              if self.verbosity > 1: print("ml fitting...")
  
              tmp_score=0
              with warnings.catch_warnings():
                  warnings.simplefilter("ignore")
                  try:
                      if self.valid_loc():
                          if self.scoring_function == roc_auc_score:
                              tmp_score =  self.roc_auc_cv(self.X[self.valid_loc(),:].transpose(),
                                                           labels)
                          else:
                              tmp_score =  np.mean(
                                  [self.scoring_function(labels[test], self.pipeline.fit(
                                  self.X[self.valid_loc(),:].transpose()[train], labels[train]).
                                  predict(self.X[self.valid_loc(),:].transpose()[test]))
                                      for train, test in KFold().split(features,
                                                                       labels)])
  
                  except ValueError as detail:
                      print("warning: ValueError in ml fit. X.shape:",
                            self.X[:,self.valid_loc()].transpose().shape,
                            "labels shape:",labels.shape)
                      print("First ten entries X:",
                            self.X[self.valid_loc(),:].transpose()[:10])
                      print("First ten entries labels:",labels[:10])
                      print("equations:",self.stacks_2_eqns(self.pop.individuals))
                      print("FEW parameters:",self.get_params())
                      print("---\ndetailed error message:",
                                                   detail)
                      raise(detail)
                  
              if self.verbosity > 1:
                  print("current ml validation score:",tmp_score)
  
              #################################################### save best model
              if self.valid_loc() and tmp_score > self._best_score:
                  self._best_estimator = copy.deepcopy(self.pipeline)
                  self._best_score = tmp_score
                  stall_count = 0;
                  self._best_inds = copy.deepcopy(self.valid())
                  if self.verbosity > 1:
                      print("updated best internal CV:",self._best_score)
              else:
                  stall_count = stall_count + 1
  
              ########################################################## variation
              if self.verbosity > 2:
                  print("variation...")
              offspring,elite,elite_index = self.variation(self.pop.individuals)
  
              ################################################# evaluate offspring
              if self.verbosity > 2:
                  print("output...")
              X_offspring = self.transform(features,offspring).transpose()
  
              if self.verbosity > 2:
                  print("fitness...")
              F_offspring = self.calc_fitness(X_offspring,
                                              labels,self.fit_choice,self.sel)
  
              ########################################################### survival
              if self.verbosity > 2: print("survival..")
              survivors,survivor_index = self.survival(self.pop.individuals, offspring, 
                                                       elite, elite_index,
                                                       X = self.X, X_O=X_offspring,
                                                       F=self.F, F_O=F_offspring)
              # set survivors
              self.pop.individuals[:] = survivors
              self.X = np.vstack((self.X, X_offspring))[survivor_index]
              if 'lexicase' in self.sel:
                  self.F = np.asarray(
                          np.vstack((self.F, F_offspring))[survivor_index],
                          order='F')
              else:
                  self.F = np.asarray(
                          np.hstack((self.F,F_offspring))[survivor_index],
                          order='F')
  
              if self.verbosity > 2:
                  print("median fitness survivors: %0.2f" % np.median(
                          [np.mean(f) for f in self.F]))
              if self.verbosity>2:
                  print("best features:",
                        self.stacks_2_eqns(self._best_inds) if self._best_inds
                        else 'original')
              pbar.set_description('Internal CV: {:1.3f}'.format(self._best_score))
              pbar.update(1)
          # end of main GP loop
              ####################
          if self.verbosity > 0: print('finished. best internal val score:'
                                       ' {:1.3f}'.format(self._best_score))
          if self.verbosity > 0: print("final model:\n",self.print_model())
  
          if not self._best_estimator:
              # if no better model found, just return underlying method fit to the
              # training data
              with warnings.catch_warnings():
                  warnings.simplefilter("ignore")
                  self._best_estimator = self.pipeline.fit(features,labels)
          else:
              # fit final estimator to all the training data
              with warnings.catch_warnings():
                  warnings.simplefilter("ignore")
                  self._best_estimator.fit(self.transform(features),labels)
          return self