How to use the ehtim.imager.Imager function in ehtim

obs = scatt.observe(eht, tint_sec, tadv_sec, tstart_hr, tstop_hr, bw_hz, sgrscat=False, ampcal=True, phasecal=True)

# Generate an image prior
npix = 55 #This must be odd
prior_fwhm = 300*eh.RADPERUAS # Gaussian size in microarcssec
gaussprior = eh.image.Image(np.zeros((npix,npix)), fov/npix, im.ra, im.dec, rf=im.rf, source=im.source, mjd=im.mjd)
gaussprior = gaussprior.add_gauss(total_flux,(prior_fwhm, prior_fwhm, 0, 0, 0))

# Now try imaging
# First image with no scattering mitigation
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)

res = obs.res() # nominal array resolution, 1/longest baseline
print("Clean beam parameters: " , beamparams)
print("Nominal Resolution: " ,res)

# Generate an image prior
npix = 128
fov = 1*im.fovx()
zbl = im.total_flux() # total flux
prior_fwhm = 200*eh.RADPERUAS # Gaussian size in microarcssec
emptyprior = eh.image.make_square(obs, npix, fov)
flatprior = emptyprior.add_flat(zbl)
gaussprior = emptyprior.add_gauss(zbl, (prior_fwhm, prior_fwhm, 0, 0, 0))

# Image total flux with bispectrum
flux = zbl
imgr  = eh.imager.Imager(obs, gaussprior, gaussprior, flux,
                          data_term={'bs':100}, show_updates=False,
                          reg_term={'simple':1,'flux':100,'cm':50},
                          maxit=200, ttype='nfft')
imgr.make_image_I()

# Blur the image with a circular beam and image again to help convergance
out = imgr.out_last()
imgr.init_next = out.blur_circ(res)
imgr.prior_next = imgr.init_next
imgr.dat_term_next = {'amp':50, 'cphase':100}
imgr.reg_term_next = {'tv2':100, 'flux':1,'cm':1}
imgr.make_image_I()

out=imgr.out_last().threshold(0.01)

# Image polarization with the polarimetric ratio

method='both', ttype=ttype, caltable=True)

# Spply the calibration solution to the full (and potentially tapered) dataset
obs_sc = caltab.applycal(obs_sc, interp='nearest', extrapolate=True)

#-------------------------------------------------------------------------------
# Reconstruct an image
#-------------------------------------------------------------------------------


# First  Round of Imaging
#-------------------------
print("Round 1: Imaging with visibility amplitudes and closure quantities...")

# Initialize imaging with a Gaussian image
imgr = eh.imager.Imager(obs_sc, gaussprior, prior_im=gaussprior,
                        flux=zbl, data_term=data_term, maxit=maxit,
                        norm_reg=True, systematic_noise=systematic_noise,
                        reg_term=reg_term, ttype=ttype, cp_uv_min=uv_zblcut, stop=stop)

# Imaging
imgr.make_image_I(show_updates=False)
converge()

# Self-calibrate to the previous model (phase-only);
# The solution_interval is 0 to align phases from high and low bands if needed
obs_sc = eh.selfcal(obs_sc, imgr.out_last(), method='phase', ttype=ttype, solution_interval=0.0)


# Second  Round of Imaging
#-------------------------
print("Round 2: Imaging with visibilities and closure quantities...")

res345 = obs345.res() # nominal array resolution, 1/longest baseline
print("Nominal Resolution: " ,res230,res345)

# Generate an image prior
npix = 64
fov = 1*im.fovx()
zbl = im.total_flux() # total flux
prior_fwhm = 200*eh.RADPERUAS # Gaussian size in microarcssec
emptyprior = eh.image.make_square(obs230, npix, fov)
flatprior = emptyprior.add_flat(zbl)
gaussprior = emptyprior.add_gauss(zbl, (prior_fwhm, prior_fwhm, 0, 0, 0))


## Image both frequencies independently with complex visibilities
flux230 = zbl
imgr230  = eh.imager.Imager(obs230, gaussprior, gaussprior, flux230,
                            data_term={'vis':100}, show_updates=False,
                            reg_term={'tv2':1.e4,'l1':5},
                            maxit=200, ttype='nfft')
imgr230.make_image_I(show_updates=True)
out230 = imgr230.out_last()


flux345 = im345.total_flux()
imgr345  = eh.imager.Imager(obs345, gaussprior, gaussprior, flux345,
                            data_term={'vis':100}, show_updates=False,
                            reg_term={'tv2':1},
                            maxit=200, ttype='nfft')
imgr345.make_image_I(show_updates=False)
out345 = imgr345.out_last()

## Image both frequencies together without multifrequency imaging (no spectral index)

reg_term={'tv2':1.e4,'l1':5},
                            maxit=200, ttype='nfft')
imgr230.make_image_I(show_updates=True)
out230 = imgr230.out_last()


flux345 = im345.total_flux()
imgr345  = eh.imager.Imager(obs345, gaussprior, gaussprior, flux345,
                            data_term={'vis':100}, show_updates=False,
                            reg_term={'tv2':1},
                            maxit=200, ttype='nfft')
imgr345.make_image_I(show_updates=False)
out345 = imgr345.out_last()

## Image both frequencies together without multifrequency imaging (no spectral index)
imgr  = eh.imager.Imager(obslist_nospec, gaussprior, gaussprior, zbl,
                            data_term={'vis':100}, show_updates=False,
                            reg_term={'tv2':1.e4,'l1':5},
                            maxit=200, ttype='nfft')
imgr.make_image_I(show_updates=False)
out_nospec = imgr.out_last()

# Image both frequencies together with spectral index
plt.close('all')
gaussprior = gaussprior.add_const_mf(3)
imgr  = eh.imager.Imager(obslist, gaussprior, gaussprior, zbl,
                            data_term={'vis':100}, show_updates=False,
                            reg_term={'tv2':1.e4,'l1':5},
                            maxit=100, ttype='nfft')
imgr.make_image_I(mf=True,show_updates=False)
out = imgr.out_last()

maxit=200, ttype='nfft')
imgr345.make_image_I(show_updates=False)
out345 = imgr345.out_last()

## Image both frequencies together without multifrequency imaging (no spectral index)
imgr  = eh.imager.Imager(obslist_nospec, gaussprior, gaussprior, zbl,
                            data_term={'vis':100}, show_updates=False,
                            reg_term={'tv2':1.e4,'l1':5},
                            maxit=200, ttype='nfft')
imgr.make_image_I(show_updates=False)
out_nospec = imgr.out_last()

# Image both frequencies together with spectral index
plt.close('all')
gaussprior = gaussprior.add_const_mf(3)
imgr  = eh.imager.Imager(obslist, gaussprior, gaussprior, zbl,
                            data_term={'vis':100}, show_updates=False,
                            reg_term={'tv2':1.e4,'l1':5},
                            maxit=100, ttype='nfft')
imgr.make_image_I(mf=True,show_updates=False)
out = imgr.out_last()

for i in range(3): # blur and reimage
    out = out.blur_circ(15*eh.RADPERUAS)
    imgr.maxit_next=500
    imgr.init_next = out
    imgr.make_image_I(mf=True,show_updates=False)
    out = imgr.out_last()

# look at results
out230_mf = out.get_image_mf(230.e9)
out345_mf = out.get_image_mf(345.e9)

gaussprior = eh.image.Image(np.zeros((npix,npix)), fov/npix, im.ra, im.dec, rf=im.rf, source=im.source, mjd=im.mjd)
gaussprior = gaussprior.add_gauss(total_flux,(prior_fwhm, prior_fwhm, 0, 0, 0))

# Now try imaging
# First image with no scattering mitigation
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()

tadv_sec = 600
tstart_hr = 0
tstop_hr = 24
bw_hz = 0.5e9
#create the observation
obs = scatt.observe(eht, tint_sec, tadv_sec, tstart_hr, tstop_hr, bw_hz, sgrscat=False, ampcal=True, phasecal=True)

# Generate an image prior
npix = 55 #This must be odd
prior_fwhm = 300*eh.RADPERUAS # Gaussian size in microarcssec
gaussprior = eh.image.Image(np.zeros((npix,npix)), fov/npix, im.ra, im.dec, rf=im.rf, source=im.source, mjd=im.mjd)
gaussprior = gaussprior.add_gauss(total_flux,(prior_fwhm, prior_fwhm, 0, 0, 0))

# Now try imaging
# First image with no scattering mitigation
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)