feat: testing M symmetry error

feat: testing M

Author: mariajmellado

frank/fourier2d.py

@@ -36,7 +36,7 @@ def __init__(self, Rmax, N, nu=0):
     def get_collocation_points(self):        
         return np.array([self.Xn, self.Yn]), np.array([self.Un, self.Vn])
 
-    def coefficients(self, u = None, v = None, direction="forward"):
+    def coefficients(self, u = None, v = None, x = None, y = None, direction="forward"):
         #start_time = time.time()
         if direction == 'forward':
             ## Normalization is dx*dy since we the DFT to be an approximation
@@ -53,14 +53,14 @@ def coefficients(self, u = None, v = None, direction="forward"):
             norm = 1 / (4*self.Xmax*self.Ymax)
             factor = 2j*np.pi
 
-            X, Y = self.Un, self.Vn
+            X, Y = u, v
             if u is None:
-                u = self.Xn
-                v = self.Yn
+                X, Y = self.Un, self.Vn
+            u = self.Xn
+            v = self.Yn
         else:
             raise AttributeError("direction must be one of {}"
                                  "".format(['forward', 'backward']))
-        
         H = norm * np.exp(factor*(np.outer(u, X) + np.outer(v, Y)))
         #print("--- %s minutes to calculate 2D-DFT coefficients---" % (time.time()/60 - start_time/60))
         return H

frank/statistical_models.py

@@ -209,22 +209,32 @@ def map_visibilities(self, u, v, V, weights, frequencies=None, geometry=None):
                 Vs = Vi[start:end]
 
                 X = self._get_mapping_coefficients(qs, ks, us, vs)
-                #print("X shape", X.shape)
-                X_CT = X.conj().T
-                #print("X_CT shape", X_CT.shape)
-                #print("w shape", np.diag(ws**(-1)).shape)
-                wXT = np.matmul(X_CT, np.diag(ws**(-1)))
-                #print("wXT shape", wXT.shape)
-                #print("V shape", Vs.shape)
-                Ms[i] += np.matmul(wXT, X)
-                js[i] += np.matmul(wXT, Vs)
+                X_CT = self._get_mapping_coefficients(qs, ks, us, vs, inverse=True)
+                wXT = np.matmul(X_CT, np.diag(ws)) # this line is the same as below.
+                #wXT = np.matmul(np.transpose(np.conjugate(X)), np.diag(ws), dtype = "complex64")
+                Ms[i] += np.matmul(wXT, X, dtype="complex128")
+                js[i] += np.matmul(wXT, Vs, dtype="complex128")
 
                 start = end
                 end = min(Ndata, end + Nstep)
-        
-        from scipy.linalg import issymmetric
-        print(Ms[0].dtype)
-        print(issymmetric(Ms[0]), "that M is a Symmetric Matrix")
+         
+        print("M type", Ms[0].dtype)
+        print("j type", js[0].dtype)
+        print("M imag-> ", " max: ", np.max(Ms[0].imag), ", min: ", np.min(Ms[0].imag) , ", mean: ", np.mean(Ms[0].imag) ,", median: ", np.median(Ms[0].imag),  ", std: ", np.std(Ms[0].imag))
+        print("M real-> ", " max: ", np.max(Ms[0].real), ", min: ", np.min(Ms[0].real) , ", mean: ", np.mean(Ms[0].real) ,", median: ", np.median(Ms[0].real),  ", std: ", np.std(Ms[0].real))
+        tol = 1e-10
+        print("with tol: ", tol, ", M is symmetric?",  np.allclose(Ms[0].real, Ms[0].T.real, atol=tol))
+        # Ms[0] = np.abs(Ms[0])
+        # from scipy.linalg import issymmetric
+        # print(issymmetric(Ms[0]), "that M is a Symmetric Matrix")
+
+
+
+        import matplotlib.pyplot as plt
+        plt.matshow(Ms[0].real, cmap="magma", vmax = np.max(Ms[0].real), vmin = np.mean(Ms[0].real))
+        plt.colorbar()
+        plt.title("M matrix, real part")
+        plt.show()
         #import sys
         #sys.exit()
 
@@ -483,7 +493,6 @@ def interpolate(self, f, r, space='Real'):
 
     def _get_mapping_coefficients(self, qs, ks, u, v, geometry=None, inverse=False):
         """Get :math:`H(q)`, such that :math:`V(q) = H(q) I_\nu`"""
-        """
         if self._vis_model == 'opt_thick':
             # Optically thick & geometrically thin
             if geometry is None:
@@ -497,16 +506,14 @@ def _get_mapping_coefficients(self, qs, ks, u, v, geometry=None, inverse=False):
             scale = np.exp(-np.outer(ks*ks, self._H2))
         else:
             raise ValueError("model not supported. Should never occur.")
-        """
         if inverse:
-            #scale = np.atleast_1d(1/scale).reshape(1,-1)
-            #qs = qs / rad_to_arcsec
+            scale = np.atleast_1d(1/scale).reshape(1,-1)
+            qs = qs / rad_to_arcsec
             direction='backward'
         else:
             direction='forward'
 
-        #H = self._DHT.coefficients(qs, direction=direction) * scale
-        H = self._DFT.coefficients(u, v, direction=direction) #* scale
+        H = self._DFT.coefficients(u, v, direction=direction)*scale
 
         return H
 
@@ -751,6 +758,7 @@ def __init__(self, DHT, M, j, p=None, scale=None, guess=None,
         S_real_inv = np.linalg.inv(S_real)
         print("--- %s minutes to calculate S_real_inv---" % (time.time()/60 - start_time/60))
         self._Sinv =  S_real_inv
+        print(self._Sinv.dtype, " Sinv dtype")
         start_time = time.time()
         self._fit()
         print("--- %s minutes to fit---" % (time.time()/60 - start_time/60))
@@ -781,8 +789,16 @@ def calculate_S_real(self, u, v, l, m, c):
         S_fspace = self.true_squared_exponential_kernel(u, v, l, m, c)
         print("--- %s minutes to calculate S---" % (time.time()/60 - start_time/60))
         start_time = time.time()
-        S_real = np.matmul(self.Ykm, np.matmul(S_fspace, self.Ykm_f))
+        S_real = np.matmul(self.Ykm, np.matmul(S_fspace, self.Ykm_f), dtype = "complex64")
         print("--- %s minutes to calculate S_real---" % (time.time()/60 - start_time/60))
+
+        print(S_real.dtype, " S_real dtype")
+        import matplotlib.pyplot as plt
+        plt.matshow(S_fspace, cmap="magma", vmin = 0, vmax = 1e-3)
+        plt.colorbar()
+        plt.title("S real matrix, real part ")
+        plt.show()
+
         return S_real
 
     def calculate_mu_cholesky(self, Dinv):
@@ -798,8 +814,8 @@ def calculate_mu_cholesky(self, Dinv):
         return mu
 
     def calculate_mu_gc(self, Dinv):
-        from scipy.sparse.linalg import bicg, bicgstab, gmres
-        method = "BiConjugate Gradient Stabilized Method"
+        from scipy.sparse.linalg import cg, bicg, bicgstab, gmres
+        method = "BiConjugate Gradient Method"
         #from scipy.sparse import csr_matrix, issparse
         #is_sparse = issparse(Dinv)
         #print("Is Dinv sparse?: ", is_sparse)