How to use ehtim - 10 common examples

datatable[i]['uvis'] = float(line[7][:-1]) * np.exp(1j*float(line[8][:-1])*DEGREE)
                datatable[i]['usigma'] = float(line[10])
                i += 1

    if not vfile==0:
        f = open(vfile)
        i = 0
        for line in f:
            line = line.split()
            if not (line[0] in ['UV', 'Scan','\n']):
                datatable[i]['vvis'] = float(line[7][:-1]) * np.exp(1j*float(line[8][:-1])*DEGREE)
                datatable[i]['vsigma'] = float(line[10])
                i += 1

    # Return the data object
    return ehtim.obsdata.Obsdata(ra, dec, rf, bw, datatable, tdata, source=src, mjd=mjd, polrep='stokes')

rf = float(file.readline().split()[2]) * 1e9
    xdim = file.readline().split()
    xdim_p = int(xdim[2])
    psize_x = float(xdim[4])*RADPERAS/xdim_p
    ydim = file.readline().split()
    ydim_p = int(ydim[2])
    psize_y = float(ydim[4])*RADPERAS/ydim_p
    file.close()

    if psize_x != psize_y:
        raise Exception("Pixel dimensions in x and y are inconsistent!")

    # Load the data, convert to list format, make object
    datatable = np.loadtxt(filename, dtype=float)
    image = datatable[:,2].reshape(ydim_p, xdim_p)
    outim = ehtim.image.Image(image, psize_x, ra, dec, rf=rf, source=src, mjd=mjd, time=time, pulse=pulse,
                              polrep='stokes',pol_prim='I')

    # Look for Stokes Q and U
    qimage = uimage = vimage = np.zeros(image.shape)
    if datatable.shape[1] == 6:
        qimage = datatable[:,3].reshape(ydim_p, xdim_p)
        uimage = datatable[:,4].reshape(ydim_p, xdim_p)
        vimage = datatable[:,5].reshape(ydim_p, xdim_p)
    elif datatable.shape[1] == 5:
        qimage = datatable[:,3].reshape(ydim_p, xdim_p)
        uimage = datatable[:,4].reshape(ydim_p, xdim_p)

    if np.any((qimage != 0) + (uimage != 0)) and np.any((vimage != 0)):
        #print('Loaded Stokes I, Q, U, and V Images')
        outim.add_qu(qimage, uimage)
        outim.add_v(vimage)

#    print('Ehtim reloaded...')
    #    #except: pass
    import ehtim as eh
    path_eh = os.path.dirname(eh.__file__)
    print('Using eht-imaging library from ', path_eh)

    if force_singlepol=='LL':
        force_singlepol='L'
    if force_singlepol=='RR':
        force_singlepol='R'

    if force_singlepol=='no':
        print('reading data without singlepol, using polrep= ',polrep)
        filen = pathf.split('/')[-1]
        if polrep in ['circ','stokes']:
            obsXX = eh.io.load.load_obs_uvfits(pathf,polrep=polrep)
            print('Polrep is ', obsXX.polrep)
        else: 
            obsXX = eh.io.load.load_obs_uvfits(pathf)
            print('Polrep unspecified')
        dfXX = obsdata_2_df(obsXX)
        if 'RR' in filen:
            dfXX['polarization'] = 'RR'    
        elif 'LL' in filen:
            dfXX['polarization'] = 'LL' 
        else:
            dfXX['polarization'] = 'WTF' 
        df = dfXX.copy()
        df['band'] = band

        #Scale sigma
        df['sigma']=scale_sigma*df['sigma']

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

#skip headline
    TableMOD = np.genfromtxt(nameF, skip_header=4)
    ScaleR = 1.
    FluxConst = 1.
    Flux = FluxConst*TableMOD[:,0]
    xPS = ScaleR*TableMOD[:,1]*np.cos(np.pi/2.-(np.pi/180.)*TableMOD[:,2])*(1.e3)*RADPERUAS #to radians
    yPS = ScaleR*TableMOD[:,1]*np.sin(np.pi/2.-(np.pi/180.)*TableMOD[:,2])*(1.e3)*RADPERUAS #to radians
    NumbPoints = np.shape(yPS)[0]

    #set image parameters
    if fov==0:
        MaxR = np.amax(TableMOD[:,1]) #in mas
        fov = 1.*MaxR*(1.e3)*RADPERUAS

    image0 = np.zeros((int(npix),int(npix)))
    im = image.Image(image0, fov/npix, 0., 0., rf=86e9)

    beamMaj = beamPar[0]
    if beamMaj==0:
        beamMaj = 4.*fov/npix

    beamMin = beamPar[1]
    if beamMin==0:
        beamMin = 4.*fov/npix

    beamTh = beamPar[2]

    sigma_maj = beamMaj / (2. * np.sqrt(2. * np.log(2.)))
    sigma_min = beamMin / (2. * np.sqrt(2. * np.log(2.)))
    cth = np.cos(beamTh)
    sth = np.sin(beamTh)
    xfov = im.xdim * im.psize

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