Commit e64541c3 by Lukas Platz

### Merge branch 'NIFTy_5' into mfcorrelatedfield_withzeromodeprior

 ... @@ -37,7 +37,7 @@ NIFTy comes with reimplemented MAP and VI estimators. ... @@ -37,7 +37,7 @@ NIFTy comes with reimplemented MAP and VI estimators. Free Theory & Implicit Operators Free Theory & Implicit Operators -------------------------------- -------------------------------- A free IFT appears when the signal field :math:{s} and the noise :math:{n} of the data :math:{d} are independent, zero-centered Gaussian processes of kown covariances :math:{S} and :math:{N}, respectively, A free IFT appears when the signal field :math:{s} and the noise :math:{n} of the data :math:{d} are independent, zero-centered Gaussian processes of known covariances :math:{S} and :math:{N}, respectively, .. math:: .. math:: ... @@ -49,7 +49,7 @@ and the measurement equation is linear in both signal and noise, ... @@ -49,7 +49,7 @@ and the measurement equation is linear in both signal and noise, d= R\, s + n, d= R\, s + n, with :math:{R} being the measurement response, which maps the continous signal field into the discrete data space. with :math:{R} being the measurement response, which maps the continuous signal field into the discrete data space. This is called a free theory, as the information Hamiltonian This is called a free theory, as the information Hamiltonian ... @@ -96,8 +96,8 @@ These implicit operators can be combined into new operators, e.g. to :math:{D^{ ... @@ -96,8 +96,8 @@ These implicit operators can be combined into new operators, e.g. to :math:{D^{ The invocation of an inverse operator applied to a vector might trigger the execution of a numerical linear algebra solver. The invocation of an inverse operator applied to a vector might trigger the execution of a numerical linear algebra solver. Thus, when NIFTy calculates :math:{m = D\, j}, it actually solves :math:{D^{-1} m = j} for :math:{m} behind the scenes. Thus, when NIFTy calculates :math:{m = D\, j}, it actually solves :math:{D^{-1} m = j} for :math:{m} behind the scenes. The advantage of implicit operators to explicit matrices is the reduced memory requirements. The advantage of implicit operators compared to explicit matrices is the reduced memory consumption; The reconstruction of only a Megapixel image would otherwithe require the storage and processing of matrices with sizes of several Terabytes. for the reconstruction of just a Megapixel image the latter would already require several Terabytes. Larger images could not be dealt with due to the quadratic memory requirements of explicit operator representations. Larger images could not be dealt with due to the quadratic memory requirements of explicit operator representations. The demo codes demos/getting_started_1.py and demos/Wiener_Filter.ipynb illustrate this. The demo codes demos/getting_started_1.py and demos/Wiener_Filter.ipynb illustrate this. ... @@ -106,7 +106,7 @@ The demo codes demos/getting_started_1.py and demos/Wiener_Filter.ipynb illu ... @@ -106,7 +106,7 @@ The demo codes demos/getting_started_1.py and demos/Wiener_Filter.ipynb illu Generative Models Generative Models ----------------- ----------------- For more sophisticated measurement situations, involving non-linear measuremnts, unknown covariances, calibration constants and the like, it is recommended to formulate those as generative models for which NIFTy provides powerful inference algorithms. For more sophisticated measurement situations (involving non-linear measurements, unknown covariances, calibration constants and the like) it is recommended to formulate those as generative models for which NIFTy provides powerful inference algorithms. In a generative model, all known or unknown quantities are described as the results of generative processes, which start with simple probability distributions, like the uniform, the i.i.d. Gaussian, or the delta distribution. In a generative model, all known or unknown quantities are described as the results of generative processes, which start with simple probability distributions, like the uniform, the i.i.d. Gaussian, or the delta distribution. ... @@ -129,7 +129,7 @@ NIFTy takes advantage of this formulation in several ways: ... @@ -129,7 +129,7 @@ NIFTy takes advantage of this formulation in several ways: 1) All prior degrees of freedom have unit covariance, which improves the condition number of operators that need to be inverted. 1) All prior degrees of freedom have unit covariance, which improves the condition number of operators that need to be inverted. 2) The amplitude operator can be regarded as part of the response, :math:{R'=R\,A}. 2) The amplitude operator can be regarded as part of the response, :math:{R'=R\,A}. In general, more sophisticated responses can be constructed out of the composition of simpler operators. In general, more sophisticated responses can be obtained by combining simpler operators. 3) The response can be non-linear, e.g. :math:{R'(s)=R \exp(A\,\xi)}, see demos/getting_started_2.py. 3) The response can be non-linear, e.g. :math:{R'(s)=R \exp(A\,\xi)}, see demos/getting_started_2.py. ... @@ -168,7 +168,7 @@ Maximum a Posteriori ... @@ -168,7 +168,7 @@ Maximum a Posteriori One popular field estimation method is Maximum a Posteriori (MAP). One popular field estimation method is Maximum a Posteriori (MAP). It only requires to minimize the information Hamiltonian, e.g by a gradient descent method that stops when It only requires minimizing the information Hamiltonian, e.g. by a gradient descent method that stops when .. math:: .. math:: ... ...
 ... @@ -103,6 +103,36 @@ class LMSpace(StructuredDomain): ... @@ -103,6 +103,36 @@ class LMSpace(StructuredDomain): def get_fft_smoothing_kernel_function(self, sigma): def get_fft_smoothing_kernel_function(self, sigma): return lambda x: self._kernel(x, sigma) return lambda x: self._kernel(x, sigma) def get_conv_kernel_from_func(self, func): """Creates a convolution kernel defined by a function. Parameters ---------- func: function This function needs to take exactly one argument, which is colatitude in radians, and return the kernel amplitude at that colatitude. Assumes the function to be radially symmetric, e.g. only dependant on theta in radians""" from .gl_space import GLSpace from ..operators.harmonic_operators import HarmonicTransformOperator import pyHealpix # define azimuthally symmetric spaces for kernel transform gl = GLSpace(self.lmax + 1, 1) lm0 = gl.get_default_codomain() theta = pyHealpix.GL_thetas(gl.nlat) # evaluate the kernel function at the required thetas kernel_sphere = Field.from_global_data(gl, func(theta)) # normalize the kernel such that the integral over the sphere is 4pi kernel_sphere = kernel_sphere * (4 * np.pi / kernel_sphere.integrate()) # compute the spherical harmonic coefficients of the kernel op = HarmonicTransformOperator(lm0, gl) kernel_lm = op.adjoint_times(kernel_sphere.weight(1)).to_global_data() # evaluate the k lengths of the harmonic space k_lengths = self.get_k_length_array().to_global_data().astype(np.int) return Field.from_global_data(self, kernel_lm[k_lengths]) @property @property def lmax(self): def lmax(self): """int : maximum allowed :math:l """int : maximum allowed :math:l ... ...
 ... @@ -90,9 +90,7 @@ class RGSpace(StructuredDomain): ... @@ -90,9 +90,7 @@ class RGSpace(StructuredDomain): def scalar_dvol(self): def scalar_dvol(self): return self._dvol return self._dvol def get_k_length_array(self): def _get_dist_array(self): if (not self.harmonic): raise NotImplementedError ibegin = dobj.ibegin_from_shape(self._shape) ibegin = dobj.ibegin_from_shape(self._shape) res = np.arange(self.local_shape[0], dtype=np.float64) + ibegin[0] res = np.arange(self.local_shape[0], dtype=np.float64) + ibegin[0] res = np.minimum(res, self.shape[0]-res)*self.distances[0] res = np.minimum(res, self.shape[0]-res)*self.distances[0] ... @@ -106,6 +104,11 @@ class RGSpace(StructuredDomain): ... @@ -106,6 +104,11 @@ class RGSpace(StructuredDomain): res = np.add.outer(res, tmp) res = np.add.outer(res, tmp) return Field.from_local_data(self, np.sqrt(res)) return Field.from_local_data(self, np.sqrt(res)) def get_k_length_array(self): if (not self.harmonic): raise NotImplementedError return self._get_dist_array() def get_unique_k_lengths(self): def get_unique_k_lengths(self): if (not self.harmonic): if (not self.harmonic): raise NotImplementedError raise NotImplementedError ... @@ -146,6 +149,27 @@ class RGSpace(StructuredDomain): ... @@ -146,6 +149,27 @@ class RGSpace(StructuredDomain): raise NotImplementedError raise NotImplementedError return lambda x: self._kernel(x, sigma) return lambda x: self._kernel(x, sigma) def get_conv_kernel_from_func(self, func): """Creates a convolution kernel defined by a function. Parameters ---------- func: function This function needs to take exactly one argument, which is distance from center (in the same units as the RGSpace distances), and return the kernel amplitude at that distance. Assumes the function to be radially symmetric, e.g. only dependant on distance""" from ..operators.harmonic_operators import HarmonicTransformOperator if (not self.harmonic): raise NotImplementedError op = HarmonicTransformOperator(self, self.get_default_codomain()) dist = op.target[0]._get_dist_array() kernel = Field.from_local_data(op.target, func(dist.local_data)) kernel = kernel / kernel.integrate() return op.adjoint_times(kernel.weight(1)) def get_default_codomain(self): def get_default_codomain(self): """Returns a :class:RGSpace object representing the (position or """Returns a :class:RGSpace object representing the (position or harmonic) partner domain of self, depending on self.harmonic. harmonic) partner domain of self, depending on self.harmonic. ... ...
 # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . # # Copyright(C) 2013-2019 Max-Planck-Society # # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. import numpy as np from ..domains.rg_space import RGSpace from ..domains.lm_space import LMSpace from ..domains.hp_space import HPSpace from ..domains.gl_space import GLSpace from .endomorphic_operator import EndomorphicOperator from .harmonic_operators import HarmonicTransformOperator from .diagonal_operator import DiagonalOperator from .simple_linear_operators import WeightApplier from ..domain_tuple import DomainTuple from ..field import Field from .. import utilities def FuncConvolutionOperator(domain, func, space=None): """Convolves input with a radially symmetric kernel defined by func Parameters ---------- domain: DomainTuple Domain of the operator. func: function This function needs to take exactly one argument, which is colatitude in radians, and return the kernel amplitude at that colatitude. space: int, optional The index of the subdomain on which the operator should act If None, it is set to 0 if domain contains exactly one space. domain[space] must be of type RGSpace, HPSpace, or GLSpace. """ domain = DomainTuple.make(domain) space = utilities.infer_space(domain, space) if not isinstance(domain[space], (RGSpace, HPSpace, GLSpace)): raise TypeError("unsupported domain") codomain = domain[space].get_default_codomain() kernel = codomain.get_conv_kernel_from_func(func) return _ConvolutionOperator(domain, kernel, space) def _ConvolutionOperator(domain, kernel, space=None): domain = DomainTuple.make(domain) space = utilities.infer_space(domain, space) if len(kernel.domain) != 1: raise ValueError("kernel needs exactly one domain") if not isinstance(domain[space], (HPSpace, GLSpace, RGSpace)): raise TypeError("need RGSpace, HPSpace, or GLSpace") lm = [d for d in domain] lm[space] = lm[space].get_default_codomain() lm = DomainTuple.make(lm) if lm[space] != kernel.domain[0]: raise ValueError("Input domain and kernel are incompatible") HT = HarmonicTransformOperator(lm, domain[space], space) diag = DiagonalOperator(kernel*domain[space].total_volume, lm, (space,)) wgt = WeightApplier(domain, space, 1) return HT(diag(HT.adjoint(wgt)))