Merge pull request #142 from achael/dev

merge dev into main 10/29/21

ehtim/__init__.py

@@ -16,6 +16,7 @@
 import ehtim.observing
 from ehtim.const_def import *
 from ehtim.imaging.imager_utils import imager_func
+from ehtim.modeling.modeling_utils import modeler_func
 import ehtim.imaging
 from ehtim.features import rex
 import ehtim.features
@@ -40,6 +41,7 @@
 import ehtim.array
 import ehtim.movie
 import ehtim.image
+import ehtim.model
 
 
 import warnings

ehtim/calibrating/network_cal.py

@@ -328,8 +328,8 @@ def errfunc(gvpar):
 
         # choose to only scale ampliltudes or phases
         if method == "phase":
-            g = g / np.abs(g)  # TODO: use exp(i*np.arg())?
-        if method == "amp":
+            g = g / np.abs(g)
+        elif method == "amp":
             g = np.abs(np.real(g))
 
         # append the default values to g for missing points
@@ -351,17 +351,23 @@ def errfunc(gvpar):
         else:
             verr = vis - g1 * g2.conj() * v_scan
 
-        nan_mask = np.array([not np.isnan(viter) for viter in verr]) * \
-            np.array([not np.isnan(viter) for viter in sigma_inv])
-        verr = verr[nan_mask * vis_mask]
-
-        chisq = np.sum((verr.real * sigma_inv[nan_mask * vis_mask])**2) + \
-            np.sum((verr.imag * sigma_inv[nan_mask * vis_mask])**2)
+        chi   = np.abs(verr) * sigma_inv
+        chisq = np.sum((chi * chi)[np.isfinite(chi) * vis_mask])
 
         # prior on the gains
         g_fracerr = gain_tol
-        chisq_g = np.sum((np.log(np.abs(g))**2 / g_fracerr**2))
-        chisq_v = np.sum((np.abs(v) / zbl_scan)**4)
+        if method == "phase":
+            chisq_g = 0 # because |g| == 1 so log(|g|) = 0
+        elif method == "amp":
+            logg    = np.log(g)
+            chisq_g = np.sum(logg * logg) / (g_fracerr * g_fracerr)
+        else:
+            logabsg = np.log(np.abs(g))
+            chisq_g = np.sum(logabsg * logabsg) / (g_fracerr * g_fracerr)
+
+        absv    = np.abs(v)
+        vv      = absv * absv
+        chisq_v = np.sum(vv * vv) / zbl_scan**4
         return chisq + chisq_g + chisq_v
 
     if np.max(g1_keys) > -1 or np.max(g2_keys) > -1:

ehtim/features/rex.py

@@ -797,7 +797,7 @@ def FindProfileSingle(imname, postprocdir,
                       rmin_search=RPRIOR_MIN, rmax_search=RPRIOR_MAX,
                       nrays_search=NRAYS_SEARCH, nrs_search=NRS_SEARCH,
                       thresh_search=THRESH, fov_search=FOVP_SEARCH, n_search=NSEARCH,
-                      flux_norm=NORMFLUX):
+                      flux_norm=NORMFLUX,center=False):
     """find the best ring profile for an image and save results
     """
 
@@ -817,7 +817,10 @@ def FindProfileSingle(imname, postprocdir,
             im_raw = im_raw.blur_circ(blur*ehc.RADPERUAS, blur*ehc.RADPERUAS)
 
         # center image and regrid to uniform pixel size and fox
-        im = di.center_core(im_raw)
+        if center:
+            im = di.center_core(im_raw) # TODO -- why isn't this working? 
+        else:
+            im = im_raw
 
         im_search = im.regrid_image(imsize, npix)
         im = im.regrid_image(imsize, npix)

ehtim/image.py

@@ -1395,7 +1395,7 @@ def blur(imarr, gauss):
             return imarr_blur
 
         # Convolve the primary image
-        imarr = (self.imvec).reshape(self.ydim, self.xdim)
+        imarr = (self.imvec).reshape(self.ydim, self.xdim).astype('float64')
         imarr_blur = blur(imarr, gauss)
 
         # Make new image object
@@ -1409,7 +1409,7 @@ def blur(imarr, gauss):
                 continue
             polvec = self._imdict[pol]
             if len(polvec):
-                polarr = polvec.reshape(self.ydim, self.xdim)
+                polarr = polvec.reshape(self.ydim, self.xdim).astype('float64')
                 if frac_pol:
                     polarr = blur(polarr, gausspol)
                 outim.add_pol_image(polarr, pol)
@@ -1418,7 +1418,7 @@ def blur(imarr, gauss):
         mflist_out = []
         for mfvec in self._mflist:
             if len(mfvec):
-                mfarr = mfvec.reshape(self.ydim, self.xdim)
+                mfarr = mfvec.reshape(self.ydim, self.xdim).astype('float64')
                 mfarr = blur(mfarr, gauss)
                 mfvec_out = mfarr.flatten()
             else:
@@ -1503,8 +1503,10 @@ def gradim(imvec):
 
             imarr = imvec.reshape(self.ydim, self.xdim)
 
-            sx = ndi.sobel(imarr, axis=0, mode='constant')
-            sy = ndi.sobel(imarr, axis=1, mode='constant')
+            #sx = ndi.sobel(imarr, axis=0, mode='constant')
+            #sy = ndi.sobel(imarr, axis=1, mode='constant')
+            sx = ndi.sobel(imarr, axis=0, mode='nearest')
+            sy = ndi.sobel(imarr, axis=1, mode='nearest')
 
             # TODO: are these in the right order??
             if gradtype == 'x':

ehtim/imager.py

@@ -316,10 +316,15 @@ def callback_func(xcur):
                     dname_key = dname + ('_%i' % i)
                 outstr += "chi2_%s : %0.2f " % (dname_key, chi2_term_dict[dname_key])
 
-        print("time: %f s" % (tstop - tstart))
-        print("J: %f" % res.fun)
-        print(outstr)
-        print(res.message.decode())
+        try:
+            print("time: %f s" % (tstop - tstart))
+            print("J: %f" % res.fun)
+            print(outstr)
+            if isinstance(res.message,str): print(res.message)
+            else: print(res.message.decode())
+        except: # TODO -- issues for some users with res.message
+            pass
+
         print("==============================")
 
         # Embed image

ehtim/imaging/dynamical_imaging.py

@@ -305,7 +305,7 @@ def merge_obs(obs_List):
     #The important things to merge are the mjd and the data
     data_merge = np.hstack([obs.data for obs in obs_List])
 
-    return obsdata.Obsdata(obs_List[0].ra, obs_List[0].dec, obs_List[0].rf, obs_List[0].bw, data_merge, obs_List[0].tarr, source=obs_List[0].source, mjd=obs_List[0].mjd, ampcal=obs_List[0].ampcal, phasecal=obs_List[0].phasecal)
+    return obsdata.Obsdata(obs_List[0].ra, obs_List[0].dec, obs_List[0].rf, obs_List[0].bw, data_merge, obs_List[0].tarr, polrep=obs_List[0].polrep, source=obs_List[0].source, mjd=obs_List[0].mjd, ampcal=obs_List[0].ampcal, phasecal=obs_List[0].phasecal)
 
 def average_im_list(im_List):
     """Return the average of a list of images

ehtim/imaging/starwarps.py

@@ -155,7 +155,7 @@ def forwardUpdates_apxImgs(mu, Lambda_orig, obs_List, A_orig, Q_orig, init_image
 
         for k in range(0,numLinIters):
 	    # F is the derivative of the Forward model with respect to the unknown parameters
-            meas, idealmeas, F, measCov, valid = getMeasurementTerms(obs_List[t], z_List_lin[t], measurement=measurement, mask=mask, normalize=normalize)
+            meas, idealmeas, F, measCov, valid = getMeasurementTerms(obs_List[t], z_List_lin[t], measurement=measurement, tot_flux=lightcurve[t], mask=mask, normalize=normalize)
             if valid:
                 z_List_t_t[t].imvec[mask], P_List_t_t[t] = prodGaussiansLem2(F, measCov, meas, z_star_List_t_tm1[t].imvec[mask], P_star_List_t_tm1[t])
 
@@ -239,7 +239,7 @@ def backwardUpdates(mu, Lambda_orig, obs_List, A_orig, Q_orig, measurement={'vis
 
         # update
 
-        meas, idealmeas, F, measCov, valid = getMeasurementTerms(obs_List[t], apxImgs[t], measurement=measurement, mask=mask, normalize=normalize)
+        meas, idealmeas, F, measCov, valid = getMeasurementTerms(obs_List[t], apxImgs[t], measurement=measurement, tot_flux=lightcurve[t], mask=mask, normalize=normalize)
 
         if valid:
             z_t_t[t].imvec[mask], P_t_t[t] = prodGaussiansLem2(F, measCov, meas, z_star_t_tp1[t].imvec[mask], P_star_t_tp1[t])
@@ -284,6 +284,9 @@ def computeSuffStatistics(mu, Lambda, obs_List, Upsilon, theta, init_x, init_y,
 
     if mu[0].xdim!=mu[0].ydim:
         error('Error: This has only been checked thus far on square images!')
+
+    if lightcurve == None  and 'flux' in measurement.keys(): #KATIE ADDED FEB 1 2021
+        error('Error: if you are using a flux constraint you must specify a lightcurve')
 
     if list(measurement.keys())==1 and measurement.keys()[0]=='vis':
         numLinIters = 1
@@ -611,7 +614,7 @@ def prodGaussiansLem2(A, Sigma, y, mu, Q):
     return (mean, covariance)
 
 
-def getMeasurementTerms(obs, im, measurement={'vis': 1}, mask=[], normalize=False):
+def getMeasurementTerms(obs, im, measurement={'vis': 1}, tot_flux=None, mask=[], normalize=False):
     if not np.sum(mask)==len(mask):
         raise ValueError('The code doenst currently work with a mask!')
 
@@ -632,7 +635,13 @@ def getMeasurementTerms(obs, im, measurement={'vis': 1}, mask=[], normalize=Fals
 
         # check to see if you have data in the current obs
         try:
-            data, sigma, A = chisqdata(obs, im, mask, dtype=dname, ttype='direct')
+            if dname=='flux':
+                if tot_flux == None:
+                    error('Error: if you are using a flux constraint you must specify a total flux (via the lightcurve)')
+                data = np.array([tot_flux])
+                sigma = np.array([1])
+            else:
+                data, sigma, A = chisqdata(obs, im, mask, dtype=dname, ttype='direct')
             count = count + 1
         except:
             continue
@@ -653,6 +662,9 @@ def getMeasurementTerms(obs, im, measurement={'vis': 1}, mask=[], normalize=Fals
         elif dname == 'logcamp':
             F = grad_logcamp(im.imvec, A)
             ideal = logcamp(im.imvec,A)
+        elif dname == 'flux':
+            F = grad_flux(im.imvec)
+            ideal = flux(im.imvec)
 
         #turn complex matrices to real
         if not np.allclose(data.imag,0):
@@ -693,6 +705,12 @@ def grad_bs(imvec, Amatrices):
     out = pt1[:,None] * Amatrices[0] + pt2[:,None] * Amatrices[1] + pt3[:,None] * Amatrices[2]
     return out
 
+def flux(imvec):
+    return np.sum(imvec)
+
+def grad_flux(imvec):
+    return np.ones((1, len(imvec)))
+
 def cphase(imvec, Amatrices):
     """the closure phase"""
     i1 = np.dot(Amatrices[0], imvec)