<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_step(self, stars, epsf=None):
 """ A single iteration of improving an ePSF. Parameters stars : `EPSFStars` object The stars used to build the ePSF. epsf : `EPSFModel` object, optional The initial ePSF model. If not input, then the ePSF will be built from scratch. Returns ------- epsf : `EPSFModel` object The updated ePSF. """


        if len(stars) < 1:
            raise ValueError('stars must contain at least one EPSFStar or '
                             'LinkedEPSFStar object.')

        if epsf is None:
            # create an initial ePSF (array of zeros)
            epsf = self._create_initial_epsf(stars)
        else:
            # improve the input ePSF
            epsf = copy.deepcopy(epsf)

        # compute a 3D stack of 2D residual images
        residuals = self._resample_residuals(stars, epsf)

        self._residuals.append(residuals)

        # compute the sigma-clipped median along the 3D stack
        with warnings.catch_warnings():
            warnings.simplefilter('ignore', category=RuntimeWarning)
            warnings.simplefilter('ignore', category=AstropyUserWarning)
            residuals = self.sigclip(residuals, axis=0, masked=False,
                                     return_bounds=False)
            if HAS_BOTTLENECK:
                residuals = bottleneck.nanmedian(residuals, axis=0)
            else:
                residuals = np.nanmedian(residuals, axis=0)

        self._residuals_sigclip.append(residuals)

        # interpolate any missing data (np.nan)
        mask = ~np.isfinite(residuals)
        if np.any(mask):
            residuals = _interpolate_missing_data(residuals, mask,
                                                  method='cubic')

            # fill any remaining nans (outer points) with zeros
            residuals[~np.isfinite(residuals)] = 0.

        self._residuals_interp.append(residuals)

        # add the residuals to the previous ePSF image
        new_epsf = epsf.normalized_data + residuals

        # smooth the ePSF
        new_epsf = self._smooth_epsf(new_epsf)

        # recenter the ePSF
        new_epsf = self._recenter_epsf(new_epsf, epsf,
                                       centroid_func=self.recentering_func,
                                       box_size=self.recentering_boxsize,
                                       maxiters=self.recentering_maxiters,
                                       center_accuracy=1.0e-4)

        # normalize the ePSF data
        new_epsf /= np.sum(new_epsf, dtype=np.float64)

        # return the new ePSF object
        xcenter = (new_epsf.shape[1] - 1) / 2.
        ycenter = (new_epsf.shape[0] - 1) / 2.

        epsf_new = EPSFModel(data=new_epsf, origin=(xcenter, ycenter),
                             normalize=False, oversampling=epsf.oversampling)

        return epsf_new