Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
Returns:
out (Image): The ensemble-average scattered image.
"""
# Inputs an unscattered image and an ensemble-average blurring kernel (2D array); returns the ensemble-average image
# The pre-computed kernel can optionally be specified (ker)
if wavelength_cm == None:
wavelength_cm = C/im.rf*100.0 #Observing wavelength [cm]
if ker is None:
ker = self.Ensemble_Average_Kernel(im, wavelength_cm, use_approximate_form)
Iim = Wrapped_Convolve((im.imvec).reshape(im.ydim, im.xdim), ker)
out = image.Image(Iim, im.psize, im.ra, im.dec, rf=C/(wavelength_cm/100.0), source=im.source, mjd=im.mjd, pulse=im.pulse)
if len(im.qvec):
Qim = Wrapped_Convolve((im.qvec).reshape(im.ydim, im.xdim), ker)
Uim = Wrapped_Convolve((im.uvec).reshape(im.ydim, im.xdim), ker)
out.add_qu(Qim, Uim)
if len(im.vvec):
Vim = Wrapped_Convolve((im.vvec).reshape(im.ydim, im.xdim), ker)
out.add_v(Vim)
return out
outcv = unpack_poltuple(res.x, xtuple, nimage, pol_solve)
if pol_prim == "amp_phase":
outcut = mcv(outcv) #change of variables
else:
outcut = outcv
if np.any(np.invert(embed_mask)):
out = embed_pol(out, embed_mask) #embed
else:
out = outcut
iimage = out[0]
qimage = make_q_image(out, pol_prim)
uimage = make_u_image(out, pol_prim)
outim = image.Image(iimage.reshape(Prior.ydim, Prior.xdim), Prior.psize,
Prior.ra, Prior.dec, rf=Prior.rf, source=Prior.source,
mjd=Prior.mjd, pulse=Prior.pulse)
outim.add_qu(qimage.reshape(Prior.ydim, Prior.xdim), uimage.reshape(Prior.ydim, Prior.xdim))
# Print stats
print("time: %f s" % (tstop - tstart))
print("J: %f" % res.fun)
print("Final Chi^2_1: %f Chi^2_2: %f" % (chisq1(outcut), chisq2(outcut)))
print(res.message)
# Return Image object
return outim
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']
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
Args:
n (int): the frame number
Returns:
(Image): the Image object of the nth frame
"""
if n<0 or n>=len(self.frames):
raise Exception("n must be in the range 0 - %i"% self.nframes)
time = self.start_hr + n * self.framedur/3600
# Make the primary image
imarr = self.frames[n].reshape(self.ydim, self.xdim)
outim = ehtim.image.Image(imarr, self.psize, self.ra, self.dec, self.pa,
polrep=self.polrep, pol_prim=self.pol_prim, time=time,
rf=self.rf, source=self.source, mjd=self.mjd, pulse=self.pulse)
# Copy over the rest of the polarizations
for pol in list(self._movdict.keys()):
if pol==self.pol_prim: continue
polframes = self._movdict[pol]
if len(polframes):
polvec = polframes[n]
polarr = polvec.reshape(self.ydim, self.xdim).copy()
outim.add_pol_image(polarr, pol)
return outim
plotcur(res.x, final=True)
# Print stats
print ("time: %f s" % (tstop - tstart))
print ("J: %f" % res.fun)
print (res.message)
#Note: the global variables are *not* released to avoid recalculation
if processes != -1:
pool.close()
#Return Frames
outim = [image.Image(Frames[i].reshape(Prior.ydim, Prior.xdim), Prior.psize,
Prior.ra, Prior.dec, rf=Obsdata_List[i].rf, source=Prior.source,
mjd=Prior.mjd, pulse=Prior.pulse) for i in range(N_frame)]
if R_flow['alpha'] == 0.0 and stochastic_optics == False:
return outim
else:
return {'Frames':outim, 'Flow':Flow, 'EpsilonList':EpsilonList }
def applyAppxWarp(im, theta, centerTheta, init_x, init_y, flowbasis_x, flowbasis_y, initTheta, method1='phase', method2='phase'):
centerIm, dImg_dTheta = calAppxWarpTerms(im, centerTheta, init_x, init_y, flowbasis_x, flowbasis_y, initTheta, method1=method1, method2=method2)
out = centerIm.imvec + np.dot(dImg_dTheta, theta - centerTheta)
outim = image.Image(np.reshape(out, (im.ydim, im.xdim)), im.psize, im.ra, im.dec, rf=im.rf, source=im.source, mjd=im.mjd, pulse=im.pulse)
return outim
if len(self.qvec)==0 or len(self.uvec)==0:
rlvec = []
lrvec = []
else:
rlvec = (self.qvec + 1j*self.uvec)
lrvec = (self.qvec - 1j*self.uvec)
imdict = {'RR':rrvec,'LL':llvec,'RL':rlvec,'LR':lrvec}
# Assemble the new image
imvec = imdict[pol_prim_out]
if len(imvec)==0:
raise Exception("for switch_polrep to %s with pol_prim_out=%s, \n"%(polrep_out,pol_prim_out) +
"output image is not defined")
newim = Image(imvec.reshape(self.ydim,self.xdim), self.psize, self.ra, self.dec, self.pa,
polrep=polrep_out, pol_prim=pol_prim_out, time=self.time,
rf=self.rf, source=self.source, mjd=self.mjd, pulse=self.pulse)
# Add in any other polarizations
for pol in list(imdict.keys()):
if pol==newim.pol_prim: continue
polvec = imdict[pol]
if len(polvec):
newim.add_pol_image(polvec.reshape(self.ydim,self.xdim), pol)
return newim
sigma = fwhm_i / (2. * np.sqrt(2. * np.log(2.)))
sigmap = sigma / self.psize
# Define a convolution function
def blur(imarr, sigma):
if np.any(np.imag(imarr)!=0):
return blur(np.real(imarr),sigma) + 1j*blur(np.imag(imarr),sigma)
imarr_blur = filt.gaussian_filter(imarr, (sigma, sigma))
return imarr_blur
# Blur the primary image
imarr = self.imvec.reshape(self.ydim, self.xdim)
imarr = blur(imarr,sigmap)
outim = Image(imarr, self.psize, self.ra, self.dec, self.pa,
polrep=self.polrep, pol_prim=self.pol_prim, time=self.time,
rf=self.rf, source=self.source, mjd=self.mjd, pulse=self.pulse)
# Blur all polarizations and copy over
for pol in list(self._imdict.keys()):
if pol==self.pol_prim: continue
polvec = self._imdict[pol]
if len(polvec):
polarr = polvec.reshape(self.ydim, self.xdim)
if fwhm_pol:
print ("Blurring polarization")
sigma = fwhm_pol / (2. * np.sqrt(2. * np.log(2.)))
sigmap = sigma / self.psize
polarr = blur(polarr, sigmap)
outim.add_pol_image(polarr, pol)
def im_new(imvec):
imarr_new = np.array([[im_new_val(imvec, x_idx , y_idx)
for x_idx in np.arange(0, -xdim_new, -1)]
for y_idx in np.arange(0, -ydim_new, -1)])
return imarr_new
# Compute new primary image vector
imarr_new = im_new(self.imvec)
# Normalize
scaling = np.sum(self.imvec) / np.sum(imarr_new)
imarr_new *= scaling
# Make new image
outim = Image(imarr_new, psize_new, self.ra, self.dec, self.pa,
polrep=self.polrep, pol_prim=self.pol_prim, time=self.time,
rf=self.rf, source=self.source, mjd=self.mjd, pulse=self.pulse)
# Interpolate all polarizations and copy over
for pol in list(self._imdict.keys()):
if pol==self.pol_prim: continue
polvec = self._imdict[pol]
if len(polvec):
polarr_new = im_new(polvec)
polarr_new *= scaling
outim.add_pol_image(polarr_new, pol)
return outim