From d5aab36a77f8b80a364866054afe5d39cfce58f7 Mon Sep 17 00:00:00 2001 From: mariajmellado Date: Tue, 16 Apr 2024 00:13:52 -0400 Subject: [PATCH] feat: finding best parameters for SE Kernel --- frank/fitExample/AS209_1mm_FrankFit.py | 44 +++++++++++ frank/radial_fitters.py | 5 +- frank/statistical_models.py | 101 ++++++++++++++++++++++++- 3 files changed, 147 insertions(+), 3 deletions(-) create mode 100644 frank/fitExample/AS209_1mm_FrankFit.py diff --git a/frank/fitExample/AS209_1mm_FrankFit.py b/frank/fitExample/AS209_1mm_FrankFit.py new file mode 100644 index 00000000..c7e0e4da --- /dev/null +++ b/frank/fitExample/AS209_1mm_FrankFit.py @@ -0,0 +1,44 @@ +import frank +import numpy as np +import matplotlib.pyplot as plt + +from frank.geometry import SourceGeometry +from frank.radial_fitters import FrankFitter, FourierBesselFitter + +# Huang 2018 +inc = 34.97 +pa = 85.76 +dra = 1.9e-3 +ddec = -2.5e-3 +r_out = 1.9 + +# Frank Parameters +n_pts = 300 +alpha = 1.3 +w_smooth = 1e-1 + +# UVtable AS209 at 1mm with removed visibilities between . +dir = "./" +data_file = dir + 'AS209_continuum_prom_1chan_30s_keepflagsFalse_removed1.txt' + +# Loading data +u, v, Re, Im, Weights = np.loadtxt(data_file, unpack = True, skiprows = 1) +vis = Re + Im*1j + +geom = SourceGeometry(inc= inc, PA= pa, dRA= dra, dDec= ddec) + +" Fitting with frank " +#FF = FrankFitter(r_out, n_pts, geom, alpha = alpha, weights_smooth = w_smooth) +FF = FourierBesselFitter(r_out, n_pts, geom) +sol = FF.fit(u, v, vis, Weights) +#setattr(sol, 'positive', sol.solve_non_negative()) + +fig = plt.figure(num = 1, figsize = (10, 5)) +plt.plot(sol.r, sol.mean / 1e10) +plt.xlabel('Radius ["]', size = 15) +plt.ylabel(r'Brightness profile [$10^{10}$ Jy sr$^{-1}$]', size = 15) +plt.title('Frank Fit AS209 1mm', size = 15) +plt.xlim(0, 1.3) +#plt.savefig('FrankFit_AS209_1mm.jpg', bbox_inches='tight') +plt.show() +plt.close() \ No newline at end of file diff --git a/frank/radial_fitters.py b/frank/radial_fitters.py index 81a0e65c..dfb6f87f 100644 --- a/frank/radial_fitters.py +++ b/frank/radial_fitters.py @@ -510,6 +510,8 @@ def _build_matrices(self, mapping): self._M = mapping['M'] self._j = mapping['j'] + self._V = mapping['V'] + self._Wvalues = mapping['W'] self._H0 = mapping['null_likelihood'] @@ -574,7 +576,8 @@ def fit(self, u, v, V, weights=1): def _fit(self): """Fit step. Computes the best fit given the pre-processed data""" fit = GaussianModel(self._DHT, self._M, self._j, - noise_likelihood=self._H0) + noise_likelihood=self._H0, + Wvalues= self._Wvalues, V = self._V) self._sol = FrankGaussianFit(self._vis_map, fit, self._info, geometry=self._geometry.clone()) diff --git a/frank/statistical_models.py b/frank/statistical_models.py index 46db07b7..70c538aa 100644 --- a/frank/statistical_models.py +++ b/frank/statistical_models.py @@ -628,9 +628,12 @@ class GaussianModel: """ def __init__(self, DHT, M, j, p=None, scale=None, guess=None, - Nfields=None, noise_likelihood=0): + Nfields=None, noise_likelihood=0, + Wvalues = None, V = None): self._DHT = DHT + self._Wvalues = Wvalues + self._V = V # Correct shape of design matrix etc. if len(M.shape) == 2: @@ -687,7 +690,19 @@ def __init__(self, DHT, M, j, p=None, scale=None, guess=None, self._Sinv[sn:en, sn:en] += Sj[n] else: - self._Sinv = None + #self._Sinv = None + q_array = self._DHT.q + + def true_squared_exponential_kernel(q, p, l): + + q1, q2 = np.meshgrid(q, q) + p1, p2 = np.meshgrid(p, p) + SE_Kernel = np.sqrt(p1 * p2) * np.exp(-0.5*(q1-q2)**2 / l**2) + return SE_Kernel + + Ykm = self._DHT.coefficients(direction="backward") + # We continue after set M matrix because is needed to calculate + # best parameters for S matrix. # Compute the design matrix self._M = np.zeros([Nr*Nfields, Nr*Nfields], dtype='f8') @@ -707,7 +722,89 @@ def __init__(self, DHT, M, j, p=None, scale=None, guess=None, self._like_noise = noise_likelihood + # M is already defined, so we find best parameters for S matrix and use it. + m, c , l = self.minimizeS() + pI = np.exp(m*np.log(q_array) + c) + S_fspace = true_squared_exponential_kernel(q_array, pI, l) + S_real = np.dot(np.transpose(Ykm), np.dot(S_fspace, Ykm)) + S_real_inv = np.linalg.inv(S_real) + self._Sinv = S_real_inv + self._fit() + + def minimizeS(self): + from scipy.optimize import minimize + from scipy.special import gamma + V = self._V + + def calculate_S(m, c, l): + q_array = self._DHT.q + p_array = c*(q_array**m) + def true_squared_exponential_kernel(q, p, l): + q1, q2 = np.meshgrid(q, q) + p1, p2 = np.meshgrid(p, p) + SE_Kernel = np.sqrt(p1 * p2) * np.exp(-0.5*(q1-q2)**2 / l**2) + return SE_Kernel + + Ykm = self._DHT.coefficients(direction="backward") + S_fspace = true_squared_exponential_kernel(q_array, p_array, l) + S_real = np.dot(np.transpose(Ykm), np.dot(S_fspace, Ykm)) + return S_real + + def calculate_D(S): + S_real_inv = np.linalg.inv(S) + Dinv = self._M + S_real_inv + D = np.linalg.inv(Dinv) + return [Dinv, D] + + def calculate_mu(Dinv): + try: + Dchol = scipy.linalg.cho_factor(Dinv) + mu = scipy.linalg.cho_solve(Dchol, self._j) + + except np.linalg.LinAlgError: + U, s, V = scipy.linalg.svd(Dinv, full_matrices=False) + s1 = np.where(s > 0, 1. / s, 0) + mu = np.dot(V.T, np.multiply(np.dot(U.T, self._j), s1)) + return mu + + def likelihood(param, data): + m, c, l = param + Wvalues = self._Wvalues + N = np.diag(1/Wvalues) + + alpha = 1.3 + l0 = 1e7 + + # Create an Inverse Gamma distribution function + def inv_gamma_function(l, alpha, beta): + return ((gamma(alpha)*beta)**(-1))*((beta/l)**(alpha + 1))*np.exp(-beta/l) + + S = calculate_S(m,c, l) + [Dinv, D] = calculate_D(S) + mu = calculate_mu(Dinv) + logdetS = np.linalg.slogdet(S)[1] + logdetD = np.linalg.slogdet(D)[1] + logdetN = np.linalg.slogdet(N)[1] + factor = np.log(2*np.pi) + + log_likelihood = 2*np.log(np.abs((1/m)*(1/c))) \ + + 2*np.log(inv_gamma_function(l, alpha, l0)) \ + - 0.5*(factor + logdetN) \ + - 0.5*(factor + logdetS) \ + + 0.5*(factor + logdetD) \ + + 0.5*np.dot(np.transpose(self._j), mu) \ + - 0.5*np.dot(np.transpose(data), np.dot(np.diag(Wvalues), data)) + return -log_likelihood + + result = minimize(likelihood, x0=np.array([-5, 60, 5e4]), args=(V,), + bounds=[(-6, 6), (1, 70), (1e3, 1e6)], + method="Powell", tol=1e-6, + ) + m, c, l = result.x + print("Result: ", "m: ", m, "c: ", c, "l: ", "{:e}".format(l)) + return [m, c, l] + def _fit(self): """Compute the mean and variance"""