How to use the ehtim.scattering.MakeEpsilonScreenFromList function in ehtim

for i in range(N_frame):
            cur_len = np.sum(embed_mask_List[i])
            log_Frames[i] = embed(x[init_i:(init_i+cur_len)], embed_mask_List[i]).reshape((N_pixel, N_pixel))
            Frames[i] = np.exp(log_Frames[i])*(embed_mask_List[i].reshape((N_pixel, N_pixel)))
            init_i += cur_len

        if R_flow['alpha'] != 0.0:
            cur_len = np.sum(embed_mask_List[0]) #assumes all the priors have the same embedding
            Flow_x = embed(x[init_i:(init_i+2*cur_len-1):2],   embed_mask_List[i]).reshape((N_pixel, N_pixel))
            Flow_y = embed(x[(init_i+1):(init_i+2*cur_len):2], embed_mask_List[i]).reshape((N_pixel, N_pixel))
            Flow = np.transpose([Flow_x.ravel(),Flow_y.ravel()]).reshape((N_pixel, N_pixel,2))
            init_i += 2*cur_len

        if stochastic_optics == True:
            EpsilonList = x[init_i:(init_i + N**2-1)]
            Epsilon_Screen = so.MakeEpsilonScreenFromList(EpsilonList, N)
            im_List = [image.Image(Frames[j], Prior.psize, Prior.ra, Prior.dec, rf=Obsdata_List[j].rf, source=Prior.source, mjd=Prior.mjd) for j in range(N_frame)]
            scatt_im_List = [scattering_model.Scatter(im_List[j], Epsilon_Screen=so.MakeEpsilonScreenFromList(EpsilonList, N), ea_ker = ea_ker[j], sqrtQ=sqrtQ, Linearized_Approximation=True).imvec for j in range(N_frame)] #the list of scattered image vectors
            Epsilon_Screen = so.MakeEpsilonScreenFromList(EpsilonList, N)
            init_i += len(EpsilonList)

        s1 = s2 = 0.0

        if alpha_s1 != 0.0:
            s1 = static_regularizer_gradient(Frames, nprior_embed_List, embed_mask_List, Prior.total_flux(), Prior.psize, entropy1, norm_reg=norm_reg, beam_size=beam_size, alpha_A=alpha_A, **kwargs)*alpha_s1
        if alpha_s2 != 0.0:
            s2 = static_regularizer_gradient(Frames, nprior_embed_List, embed_mask_List, Prior.total_flux(), Prior.psize, entropy2, norm_reg=norm_reg, beam_size=beam_size, alpha_A=alpha_A, **kwargs)*alpha_s2

        s_dynamic_grad = cm_grad = flux_grad = s_dS = s_dF = 0.0
        if R_dI['alpha'] != 0.0: s_dynamic_grad += RdI_gradient(Frames,**R_dI)*R_dI['alpha']
        if R_dt['alpha'] != 0.0: s_dynamic_grad += Rdt_gradient(Frames, B_dt, **R_dt)*R_dt['alpha']

imgr = eh.imager.Imager(obs, gaussprior, prior_im=gaussprior, maxit=200, flux=total_flux, clipfloor=-1.)
imgr.make_image_I()
imgr.out_last().display()

# Now try deblurring before imaging
imgr_deblur = eh.imager.Imager(sm.Deblur_obs(obs), gaussprior, prior_im=gaussprior, maxit=200, flux=total_flux, clipfloor=-1.)
imgr_deblur.make_image_I()
imgr_deblur.out_last().display()

# Now image using stochastic optics
imgr_so = eh.imager.Imager(obs, gaussprior, prior_im=gaussprior, maxit=200, flux=total_flux, clipfloor=-1.)
imgr_so.make_image_I_stochastic_optics()
# Now look at the unscattered image, the scattered image, and replicate the scattering using the solved screen
imgr_so.out_last().display()
imgr_so.out_scattered_last().display()
imgr_so.scattering_model.Scatter(imgr_so.out_last(), Epsilon_Screen=so.MakeEpsilonScreenFromList(imgr_so.out_epsilon_last(), npix), ea_ker = imgr_so._ea_ker, sqrtQ=imgr_so._sqrtQ, Linearized_Approximation=True, DisplayImage=True)

#Note that only the scattered image will fit the measured visibilities!
eh.comp_plots.plotall_obs_im_compare(obs, imgr_so.out_last(), 'uvdist', 'amp')
eh.comp_plots.plotall_obs_im_compare(obs, imgr_so.out_scattered_last(), 'uvdist', 'amp')

# Decrease the scattering regularization slightly and re-image (desired max |Epsilon| is ~2.5)
imgr_so.alpha_phi_next /= 2.0
imgr_so.init_next = imgr_so.out_last().blur_circ(obs.res())
imgr_so.epsilon_list_next = imgr_so.out_epsilon_last()
imgr_so.make_image_I_stochastic_optics()
imgr_so.out_last().display()

res = opt.minimize(self.objfunc_scattering, self._xinit, method='L-BFGS-B',
                               options=optdict, callback=self.plotcur_scattering)
        tstop = time.time()

        # Format output
        out = res.x[:N**2]
        if self.transform_next == 'log': out = np.exp(out)
        if np.any(np.invert(self._embed_mask)):
            raise Exception("Embedding is not currently implemented!")
            out = embed(out, self._embed_mask)

        outim = image.Image(out.reshape(N, N), self.prior_next.psize,
                            self.prior_next.ra, self.prior_next.dec, self.prior_next.pa,
                            rf=self.prior_next.rf, source=self.prior_next.source, mjd=self.prior_next.mjd, pulse=self.prior_next.pulse)
        outep = res.x[N**2:]
        outscatt = self.scattering_model.Scatter(outim, Epsilon_Screen=so.MakeEpsilonScreenFromList(outep, N),
                                                 ea_ker = self._ea_ker, sqrtQ=self._sqrtQ,
                                                 Linearized_Approximation=True)

        # Preserving image complex polarization fractions
        if len(self.prior_next.qvec):
            qvec = self.prior_next.qvec * out / self.prior_next.imvec
            uvec = self.prior_next.uvec * out / self.prior_next.imvec
            outim.add_qu(qvec.reshape(N, N),
                         uvec.reshape(N, N))

        # Print stats
        print("time: %f s" % (tstop - tstart))
        print("J: %f" % res.fun)
        print(res.message)

        # Append to history

def plotcur_scattering(self, minvec):
        if self._show_updates:
            if self._nit % self._update_interval == 0:
                N = self.prior_next.xdim

                imvec       = minvec[:N**2]
                EpsilonList = minvec[N**2:]
                if self.transform_next == 'log':
                    imvec = np.exp(imvec)

                IM = ehtim.image.Image(imvec.reshape(N,N), self.prior_next.psize,
                                       self.prior_next.ra, self.prior_next.dec, self.prior_next.pa,
                                       rf=self.obs_next.rf, source=self.prior_next.source, mjd=self.prior_next.mjd)
                #the scattered image vector
                scatt_im = self.scattering_model.Scatter(IM, Epsilon_Screen=so.MakeEpsilonScreenFromList(EpsilonList, N),
                                                         ea_ker = self._ea_ker, sqrtQ=self._sqrtQ, Linearized_Approximation=True).imvec

                # Calculate the chi^2 using the scattered image
                datterm = 0.
                chi2_term_dict = self.make_chisq_dict(scatt_im)
                for dname in sorted(self.dat_term_next.keys()):
                    datterm += self.dat_term_next[dname] * (chi2_term_dict[dname] - 1.)

                # Calculate the entropy using the unscattered image
                regterm = 0
                reg_term_dict = self.make_reg_dict(imvec)
                for regname in sorted(self.reg_term_next.keys()):
                    regterm += self.reg_term_next[regname] * reg_term_dict[regname]

                # Scattering screen regularization term
                chisq_epsilon = sum(EpsilonList*EpsilonList)/((N*N-1.0)/2.0)

log_Frames[i] = embed(x[init_i:(init_i+cur_len)], embed_mask_List[i]).reshape((N_pixel, N_pixel))
                Frames[i] = np.exp(log_Frames[i])*(embed_mask_List[i].reshape((N_pixel, N_pixel)))
                init_i += cur_len

            if R_flow['alpha'] != 0.0:
                cur_len = np.sum(embed_mask_List[0]) #assumes all the priors have the same embedding
                Flow_x = embed(x[init_i:(init_i+2*cur_len-1):2],   embed_mask_List[i]).reshape((N_pixel, N_pixel))
                Flow_y = embed(x[(init_i+1):(init_i+2*cur_len):2], embed_mask_List[i]).reshape((N_pixel, N_pixel))
                Flow = np.transpose([Flow_x.ravel(),Flow_y.ravel()]).reshape((N_pixel, N_pixel,2))
                init_i += 2*cur_len

            if stochastic_optics == True:
                EpsilonList = x[init_i:(init_i + N**2-1)]
                im_List = [image.Image(Frames[j], Prior.psize, Prior.ra, Prior.dec, rf=Obsdata_List[j].rf, source=Prior.source, mjd=Prior.mjd) for j in range(N_frame)]

                scatt_im_List = [scattering_model.Scatter(im_List[j], Epsilon_Screen=so.MakeEpsilonScreenFromList(EpsilonList, N), ea_ker = ea_ker[j], sqrtQ=sqrtQ, Linearized_Approximation=True).imvec
                                 for j in range(N_frame)] #the list of scattered image vectors

            s1 = s2 = 0.0

            if alpha_s1 != 0.0:

                s1 = static_regularizer(Frames, nprior_embed_List, embed_mask_List, Prior.total_flux(),
                                        Prior.psize, entropy1, norm_reg=norm_reg, beam_size=beam_size, alpha_A=alpha_A, **kwargs)*alpha_s1
            if alpha_s2 != 0.0:
                s2 = static_regularizer(Frames, nprior_embed_List, embed_mask_List, Prior.total_flux(),
                                        Prior.psize, entropy2, norm_reg=norm_reg, beam_size=beam_size, alpha_A=alpha_A, **kwargs)*alpha_s2

            s_dynamic = cm = s_dS = s_dF = 0.0

            if R_dI['alpha'] != 0.0: s_dynamic += RdI(Frames, **R_dI)*R_dI['alpha']
            if R_dt['alpha'] != 0.0: s_dynamic += Rdt(Frames, B_dt, **R_dt)*R_dt['alpha']

wavelengthbar = wavelength/(2.0*np.pi) #lambda/(2pi) [cm]
        N = self.prior_next.xdim
        #Field of view, in cm, at the scattering screen
        FOV = self.prior_next.psize * N * self.scattering_model.observer_screen_distance
        rF = self.scattering_model.rF(wavelength)

        imvec       = minvec[:N**2]
        EpsilonList = minvec[N**2:]
        if self.transform_next == 'log':
            imvec = np.exp(imvec)

        IM = ehtim.image.Image(imvec.reshape(N,N), self.prior_next.psize,
                               self.prior_next.ra, self.prior_next.dec, self.prior_next.pa,
                               rf=self.obs_next.rf, source=self.prior_next.source, mjd=self.prior_next.mjd)
        #the scattered image vector
        scatt_im = self.scattering_model.Scatter(IM, Epsilon_Screen=so.MakeEpsilonScreenFromList(EpsilonList, N),
                                                 ea_ker = self._ea_ker, sqrtQ=self._sqrtQ,
                                                 Linearized_Approximation=True).imvec

        EA_Image = self.scattering_model.Ensemble_Average_Blur(IM, ker = self._ea_ker)
        EA_Gradient = so.Wrapped_Gradient((EA_Image.imvec/(FOV/N)).reshape(N, N))
        #The gradient signs don't actually matter, but let's make them match intuition (i.e., right to left, bottom to top)
        EA_Gradient_x = -EA_Gradient[1]
        EA_Gradient_y = -EA_Gradient[0]

        Epsilon_Screen = so.MakeEpsilonScreenFromList(EpsilonList, N)
        phi = self.scattering_model.MakePhaseScreen(Epsilon_Screen, IM, obs_frequency_Hz=self.obs_next.rf,sqrtQ_init=self._sqrtQ).imvec.reshape((N, N))
        phi_Gradient = so.Wrapped_Gradient(phi/(FOV/N))
        phi_Gradient_x = -phi_Gradient[1]
        phi_Gradient_y = -phi_Gradient[0]

        #Entropy gradient; wrt unscattered image so unchanged by scattering

log_Frames[i] = embed(x[init_i:(init_i+cur_len)], embed_mask_List[i]).reshape((N_pixel, N_pixel))
            Frames[i] = np.exp(log_Frames[i])*(embed_mask_List[i].reshape((N_pixel, N_pixel)))
            init_i += cur_len

        if R_flow['alpha'] != 0.0:
            cur_len = np.sum(embed_mask_List[0]) #assumes all the priors have the same embedding
            Flow_x = embed(x[init_i:(init_i+2*cur_len-1):2],   embed_mask_List[i]).reshape((N_pixel, N_pixel))
            Flow_y = embed(x[(init_i+1):(init_i+2*cur_len):2], embed_mask_List[i]).reshape((N_pixel, N_pixel))
            Flow = np.transpose([Flow_x.ravel(),Flow_y.ravel()]).reshape((N_pixel, N_pixel,2))
            init_i += 2*cur_len

        if stochastic_optics == True:
            EpsilonList = x[init_i:(init_i + N**2-1)]
            Epsilon_Screen = so.MakeEpsilonScreenFromList(EpsilonList, N)
            im_List = [image.Image(Frames[j], Prior.psize, Prior.ra, Prior.dec, rf=Obsdata_List[j].rf, source=Prior.source, mjd=Prior.mjd) for j in range(N_frame)]
            scatt_im_List = [scattering_model.Scatter(im_List[j], Epsilon_Screen=so.MakeEpsilonScreenFromList(EpsilonList, N), ea_ker = ea_ker[j], sqrtQ=sqrtQ, Linearized_Approximation=True).imvec for j in range(N_frame)] #the list of scattered image vectors
            Epsilon_Screen = so.MakeEpsilonScreenFromList(EpsilonList, N)
            init_i += len(EpsilonList)

        s1 = s2 = 0.0

        if alpha_s1 != 0.0:
            s1 = static_regularizer_gradient(Frames, nprior_embed_List, embed_mask_List, Prior.total_flux(), Prior.psize, entropy1, norm_reg=norm_reg, beam_size=beam_size, alpha_A=alpha_A, **kwargs)*alpha_s1
        if alpha_s2 != 0.0:
            s2 = static_regularizer_gradient(Frames, nprior_embed_List, embed_mask_List, Prior.total_flux(), Prior.psize, entropy2, norm_reg=norm_reg, beam_size=beam_size, alpha_A=alpha_A, **kwargs)*alpha_s2

        s_dynamic_grad = cm_grad = flux_grad = s_dS = s_dF = 0.0
        if R_dI['alpha'] != 0.0: s_dynamic_grad += RdI_gradient(Frames,**R_dI)*R_dI['alpha']
        if R_dt['alpha'] != 0.0: s_dynamic_grad += Rdt_gradient(Frames, B_dt, **R_dt)*R_dt['alpha']

        if alpha_dS1 != 0.0: s_dS += RdS_gradient(Frames, nprior_embed_List, embed_mask_List, entropy1, norm_reg, beam_size=beam_size, alpha_A=alpha_A)*alpha_dS1
        if alpha_dS2 != 0.0: s_dS += RdS_gradient(Frames, nprior_embed_List, embed_mask_List, entropy2, norm_reg, beam_size=beam_size, alpha_A=alpha_A)*alpha_dS2