<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 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