### Merge branch 'NIFTy_5' of https://gitlab.mpcdf.mpg.de/ift/nifty into NIFTy_5

parents 7999e816 890a6b13
 ... ... @@ -10,7 +10,7 @@ RUN apt-get update && apt-get install -y \ # Testing dependencies python3-pytest-cov jupyter \ # Optional NIFTy dependencies libfftw3-dev python3-mpi4py python3-matplotlib \ libfftw3-dev python3-mpi4py python3-matplotlib python3-pynfft \ # more optional NIFTy dependencies && pip3 install pyfftw \ && pip3 install git+https://gitlab.mpcdf.mpg.de/ift/pyHealpix.git \ ... ...
This diff is collapsed.
 ... ... @@ -38,3 +38,17 @@ Support for spherical harmonic transforms is added via:: MPI support is added via:: sudo apt-get install python3-mpi4py NIFTy documentation is provided by Sphinx. To build the documentation:: sudo apt-get install python3-sphinx-rtd-theme dvipng cd sh docs/generate.sh To view the documentation in firefox:: firefox docs/build/index.html (Note: Make sure that you reinstall nifty after each change since sphinx imports nifty from the Python path.)
 Discretization and Volume in NIFTy Discretisation and Volume in NIFTy ================================== .. note:: Some of this discussion is rather technical and may be skipped in a first read-through. ... ... @@ -160,15 +160,21 @@ Often, log-likelihoods contain integrals over the quantity of interest :math:s \int_\Omega \text{d}x\, s(x) \approx \sum_i s^i\int_{\Omega_i}\text{d}x\, 1 Here the domain of the integral :math:\Omega = \dot{\bigcup_q} \; \Omega_i is the disjoint union over smaller :math:\Omega_i, e.g. the pixels of the space, and :math:s_i is the discretized field value on the :math:i-th pixel. Here the domain of the integral :math:\Omega = \dot{\bigcup_q} \; \Omega_i is the disjoint union over smaller :math:\Omega_i, e.g. the pixels of the space, and :math:s_i is the discretised field value on the :math:i-th pixel. This introduces the weighting :math:V_i=\int_{\Omega_i}\text{d}x\, 1, also called the volume factor, a property of the space. NIFTy aids you in constructing your own log-likelihood by providing methods like :func:~field.Field.weight, which weights all pixels of a field with their corresponding volume. An integral over a :class:~field.Field :code:s can be performed by calling :code:s.weight(1).sum(), which is equivalent to :code:s.integrate(). Volume factors are also applied automatically in the following places: - :class:~operators.harmonic_operators.FFTOperator as well as all other harmonic operators. Here the zero mode of the transformed field is the integral over the original field, thus the whole field is weighted once. - some response operators, such as the :class:~library.los_response.LOSResponse. In this operator a line integral is descritized, so a 1-dimensional volume factor is applied. - In :class:~library.correlated_fields.CorrelatedField as well :class:~library.correlated_fields.MfCorrelatedField, the field is multiplied by the square root of the total volume in configuration space. This ensures that the same field reconstructed over a larger domain has the same variance in position space in the limit of infinite resolution. It also ensures that power spectra in NIFTy behave according to the definition of a power spectrum, namely the power of a k-mode is the expectation of the k-mode square, divided by the volume of the space. - :class:~operators.harmonic_operators.FFTOperator as well as all other harmonic operators. Here the zero mode of the transformed field is the integral over the original field, thus the whole field is weighted once. - Some response operators, such as the :class:~library.los_response.LOSResponse. In this operator a line integral is discretised, so a 1-dimensional volume factor is applied. - In :class:~library.correlated_fields.CorrelatedField as well as :class:~library.correlated_fields.MfCorrelatedField. Both describe fields with a smooth, a priori unknown correlation structure specified by a power spectrum. The field is multiplied by the square root of the total volume of it domain's harmonic counterpart. This ensures that the same power spectrum can be used regardless of the chosen resolution, provided the total volume of the space remains the same. It also guarantees that the power spectra in NIFTy behave according to their definition, i.e. the power of a mode :math:s_k is the expectation value of that mode squared, divided by the volume of its space :math:P(k) = \left\langle s_k^2 \right\rangle / V_k. Note that in contrast to some older versions of NIFTy, the dot product :code:s.vdot(t) of fields does **not** apply a volume factor, but instead just sums over the field components, ... ...
 ... ... @@ -74,7 +74,8 @@ from .minimization.metric_gaussian_kl import MetricGaussianKL from .sugar import * from .plot import Plot from .library.smooth_linear_amplitude import SLAmplitude, CepstrumOperator from .library.smooth_linear_amplitude import ( SLAmplitude, LinearSLAmplitude, CepstrumOperator) from .library.inverse_gamma_operator import InverseGammaOperator from .library.los_response import LOSResponse from .library.dynamic_operator import (dynamic_operator, ... ... @@ -85,6 +86,7 @@ from .library.wiener_filter_curvature import WienerFilterCurvature from .library.correlated_fields import CorrelatedField, MfCorrelatedField from .library.adjust_variances import (make_adjust_variances_hamiltonian, do_adjust_variances) from .library.nfft import NFFT from . import extra ... ...
 ... ... @@ -149,7 +149,7 @@ def ensure_default_distributed(arr): def absmax(arr): return np.linalg.norm(arr.rehape(-1), ord=np.inf) return np.linalg.norm(arr.reshape(-1), ord=np.inf) def norm(arr, ord=2): ... ...
 # 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) 2018-2019 Max-Planck-Society # # Resolve is being developed at the Max-Planck-Institut fuer Astrophysik. import numpy as np import nifty5 as ift class NFFT(ift.LinearOperator): """Performs a non-equidistant Fourier transform, i.e. a Fourier transform followed by a degridding operation. Parameters ---------- domain : RGSpace Domain of the operator. It has to be two-dimensional and have shape (2N, 2N). The coordinates of the lower left pixel of the dirty image are (-N,-N), and of the upper right pixel (N-1,N-1). uv : numpy.ndarray 2D numpy array of type float64 and shape (M,2), where M is the number of measurements. uv[i,0] and uv[i,1] contain the u and v coordinates of measurement #i, respectively. All coordinates must lie in the range [-0.5; 0,5[. """ def __init__(self, domain, uv): from pynfft.nfft import NFFT npix = domain.shape assert npix == domain.shape assert len(domain.shape) == 2 assert type(npix) == int, "npix must be integer" assert npix > 0 and ( npix % 2) == 0, "npix must be an even, positive integer" assert isinstance(uv, np.ndarray), "uv must be a Numpy array" assert uv.dtype == np.float64, "uv must be an array of float64" assert uv.ndim == 2, "uv must be a 2D array" assert uv.shape > 0, "at least one point needed" assert uv.shape == 2, "the second dimension of uv must be 2" assert np.all(uv >= -0.5) and np.all(uv <= 0.5),\ "all coordinates must lie between -0.5 and 0.5" self._domain = ift.DomainTuple.make(domain) self._target = ift.DomainTuple.make( ift.UnstructuredDomain(uv.shape)) self._capability = self.TIMES | self.ADJOINT_TIMES self.npt = uv.shape self.plan = NFFT(self.domain.shape, self.npt, m=6) self.plan.x = uv self.plan.precompute() def apply(self, x, mode): self._check_input(x, mode) if mode == self.TIMES: self.plan.f_hat = x.to_global_data() res = self.plan.trafo().copy() else: self.plan.f = x.to_global_data() res = self.plan.adjoint().copy() return ift.Field.from_global_data(self._tgt(mode), res)
 ... ... @@ -169,6 +169,18 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']): which returns on its target a power spectrum which consists out of a smooth and a linear part. ''' return LinearSLAmplitude(target=target, n_pix=n_pix, a=a, k0=k0, sm=sm, sv=sv, im=im, iv=iv, keys=keys).exp() def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']): '''LinearOperator for parametrizing smooth log-amplitudes (square roots of power spectra). Logarithm of SLAmplitude See documentation of SLAmplitude for more details ''' if not (isinstance(n_pix, int) and isinstance(target, PowerSpace)): raise TypeError ... ... @@ -196,4 +208,4 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']): loglog_ampl = 0.5*(smooth + linear) # Go from loglog-space to linear-linear-space return et @ loglog_ampl.exp() return et @ loglog_ampl
 ... ... @@ -45,6 +45,13 @@ def FuncConvolutionOperator(domain, func, space=None): 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. Notes ----- The operator assumes periodic boundaries in the input domain. This means for a sufficiently broad function a point source close to the boundary will blur into the opposite side of the image. Zero padding can be applied to avoid this behaviour. """ domain = DomainTuple.make(domain) space = utilities.infer_space(domain, space) ... ...
 ... ... @@ -20,6 +20,7 @@ import numpy as np from .. import utilities from ..domain_tuple import DomainTuple from ..field import Field from ..multi_field import MultiField from ..linearization import Linearization from ..sugar import makeDomain, makeOp from .linear_operator import LinearOperator ... ... @@ -121,7 +122,7 @@ class GaussianEnergy(EnergyOperator): """ def __init__(self, mean=None, covariance=None, domain=None): if mean is not None and not isinstance(mean, Field): if mean is not None and not isinstance(mean, (Field, MultiField)): raise TypeError if covariance is not None and not isinstance(covariance, LinearOperator): ... ... @@ -311,7 +312,6 @@ class StandardHamiltonian(EnergyOperator): Tells an internal :class:SamplingEnabler which convergence criterion to use to draw Gaussian samples. See also -------- Encoding prior knowledge in the structure of the likelihood, ... ...
 ... ... @@ -361,6 +361,13 @@ class HarmonicTransformOperator(LinearOperator): The index of the domain on which the operator should act If None, it is set to 0 if domain contains exactly one subdomain. domain[space] must be a harmonic domain. Notes ----- HarmonicTransformOperator uses a Hartley transformation to transform between harmonic and non-harmonic RGSpaces. This has the advantage that all field values are real in either space. If you require a true Fourier transform you should use FFTOperator instead. """ def __init__(self, domain, target=None, space=None): ... ...
 ... ... @@ -384,7 +384,6 @@ class _OpSum(Operator): v = x._val if lin else x v1 = v.extract(self._op1.domain) v2 = v.extract(self._op2.domain) res = None if not lin: return self._op1(v1).unite(self._op2(v2)) wm = x.want_metric ... ... @@ -393,7 +392,10 @@ class _OpSum(Operator): op = lin1._jac._myadd(lin2._jac, False) res = lin1.new(lin1._val.unite(lin2._val), op(x.jac)) if lin1._metric is not None and lin2._metric is not None: res = res.add_metric(lin1._metric + lin2._metric) from .sandwich_operator import SandwichOperator met = lin1._metric._myadd(lin2._metric, False) met = SandwichOperator.make(x.jac, met) res = res.add_metric(met) return res def _simplify_for_constant_input_nontrivial(self, c_inp): ... ...
 ... ... @@ -66,9 +66,12 @@ class ScalingOperator(EndomorphicOperator): return x if fct == 0.: return full(self.domain, 0.) if (mode & 10) != 0: MODES_WITH_ADJOINT = self.ADJOINT_TIMES | self.ADJOINT_INVERSE_TIMES MODES_WITH_INVERSE = self.INVERSE_TIMES | self.ADJOINT_INVERSE_TIMES if (mode & MODES_WITH_ADJOINT) != 0: fct = np.conj(fct) if (mode & 12) != 0: if (mode & MODES_WITH_INVERSE) != 0: fct = 1./fct return x*fct ... ...
 ... ... @@ -29,6 +29,7 @@ from .multi_field import MultiField from .operators.block_diagonal_operator import BlockDiagonalOperator from .operators.diagonal_operator import DiagonalOperator from .operators.distributors import PowerDistributor from .plot import Plot __all__ = ['PS_field', 'power_analyze', 'create_power_operator', 'create_harmonic_smoothing_operator', 'from_random', ... ... @@ -37,7 +38,7 @@ __all__ = ['PS_field', 'power_analyze', 'create_power_operator', 'sin', 'cos', 'tan', 'sinh', 'cosh', 'absolute', 'one_over', 'clip', 'sinc', 'conjugate', 'get_signal_variance', 'makeOp', 'domain_union', 'get_default_codomain'] 'get_default_codomain', 'single_plot'] def PS_field(pspace, func): ... ... @@ -434,3 +435,14 @@ def get_default_codomain(domainoid, space=None): ret = [dom for dom in domainoid] ret[space] = domainoid[space].get_default_codomain() return DomainTuple.make(ret) def single_plot(field, **kwargs): """Creates a single plot using Plot. Keyword arguments are passed to both Plot.add and Plot.output. """ p = Plot() p.add(field, **kwargs) if 'title' in kwargs: del(kwargs['title']) p.output(**kwargs)
 ... ... @@ -71,6 +71,20 @@ def testLinearInterpolator(): ift.extra.consistency_check(op) @pmp('sp', _h_spaces + _p_spaces + _pow_spaces) def testRealizer(sp): op = ift.Realizer(sp) ift.extra.consistency_check(op, np.complex128, np.float64, only_r_linear=True) @pmp('sp', _h_spaces + _p_spaces + _pow_spaces) def testConjugationOperator(sp): op = ift.ConjugationOperator(sp) ift.extra.consistency_check(op, np.complex128, np.complex128, only_r_linear=True) @pmp('args', [(ift.RGSpace(10, harmonic=True), 4, 0), (ift.RGSpace( (24, 31), distances=(0.4, 2.34), harmonic=True), 3, 0), (ift.LMSpace(4), 10, 0)]) ... ... @@ -279,3 +293,10 @@ def testValueInserter(sp, seed): ind.append(np.random.randint(0, ss-1)) op = ift.ValueInserter(sp, ind) ift.extra.consistency_check(op) def testNFFT(): dom = ift.RGSpace(2*(16,)) uv = np.array([[.2, .4], [-.22, .452]]) op = ift.NFFT(dom, uv) ift.extra.consistency_check(op)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment