docs: def build_epsf(self, stars, init_epsf=None):
          """
          Iteratively build an ePSF from star cutouts.
  
          Parameters
          ----------
          stars : `EPSFStars` object
              The stars used to build the ePSF.
  
          init_epsf : `EPSFModel` object, optional
              The initial ePSF model.  If not input, then the ePSF will be
              built from scratch.
  
          Returns
          -------
          epsf : `EPSFModel` object
              The constructed ePSF.
  
          fitted_stars : `EPSFStars` object
              The input stars with updated centers and fluxes derived
              from fitting the output ``epsf``.
          """
  
          iter_num = 0
          center_dist_sq = self.center_accuracy_sq + 1.
          centers = stars.cutout_center_flat
          n_stars = stars.n_stars
          fit_failed = np.zeros(n_stars, dtype=bool)
          dx_dy = np.zeros((n_stars, 2), dtype=np.float)
          epsf = init_epsf
          dt = 0.
  
          while (iter_num < self.maxiters and
                  np.max(center_dist_sq) >= self.center_accuracy_sq and
                  not np.all(fit_failed)):
  
              t_start = time.time()
              iter_num += 1
  
              if self.progress_bar:
                  if iter_num == 1:
                      dt_str = ' [? s/iter]'
                  else:
                      dt_str = ' [{:.1f} s/iter]'.format(dt)
  
                  print('PROGRESS: iteration {0:d} (of max {1}){2}'
                        .format(iter_num, self.maxiters, dt_str), end='\r')
  
              # build/improve the ePSF
              epsf = self._build_epsf_step(stars, epsf=epsf)
  
              # fit the new ePSF to the stars to find improved centers
              # we catch fit warnings here -- stars with unsuccessful fits
              # are excluded from the ePSF build process
              with warnings.catch_warnings():
                  message = '.*The fit may be unsuccessful;.*'
                  warnings.filterwarnings('ignore', message=message,
                                          category=AstropyUserWarning)
                  stars = self.fitter(epsf, stars)
  
              # find all stars where the fit failed
              fit_failed = np.array([star._fit_error_status > 0
                                     for star in stars.all_stars])
              if np.all(fit_failed):
                  raise ValueError('The ePSF fitting failed for all stars.')
  
              # permanently exclude fitting any star where the fit fails
              # after 3 iterations
              if iter_num > 3 and np.any(fit_failed):
                  idx = fit_failed.nonzero()[0]
                  for i in idx:
                      stars.all_stars[i]._excluded_from_fit = True
  
              dx_dy = stars.cutout_center_flat - centers
              dx_dy = dx_dy[np.logical_not(fit_failed)]
              center_dist_sq = np.sum(dx_dy * dx_dy, axis=1, dtype=np.float64)
              centers = stars.cutout_center_flat
  
              self._nfit_failed.append(np.count_nonzero(fit_failed))
              self._center_dist_sq.append(center_dist_sq)
              self._max_center_dist_sq.append(np.max(center_dist_sq))
              self._epsf.append(epsf)
  
              dt = time.time() - t_start
  
          return epsf, stars