diff --git a/frank/fourier2d.py b/frank/fourier2d.py index 8aeb084a..b0c73308 100644 --- a/frank/fourier2d.py +++ b/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 diff --git a/frank/statistical_models.py b/frank/statistical_models.py index 4552bb07..8d4ce21f 100644 --- a/frank/statistical_models.py +++ b/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)