Merge branch 'mfcorrelatedfield_withzeromodeprior' of...

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

demos/getting_started_1.py

@@ -93,12 +93,10 @@ if __name__ == '__main__':
 
     # Apply the mapping
     prior_correlation_structure = PD(ift.PS_field(power_space, power_spectrum))
-    print(prior_correlation_structure)
 
     # Insert the result into the diagonal of an harmonic space operator
     S = ift.DiagonalOperator(prior_correlation_structure)
     # S is the prior field covariance
-    assert False
 
     # Build instrument response consisting of a discretization, mask
     # and harmonic transformaion

demos/misc/convolution.py

+import numpy as np
+import nifty5 as ift
+
+
+def convtest(test_signal, delta, func):
+    domain = test_signal.domain
+
+    # Create Convolution Operator
+    conv_op = ift.FuncConvolutionOperator(domain, func)
+
+    # Convolve, Adjoint-Convolve
+    conv_signal = conv_op(test_signal)
+    cac_signal = conv_op.adjoint_times(conv_signal)
+
+    print(test_signal.integrate(), conv_signal.integrate(),
+          cac_signal.integrate())
+
+    # generate kernel image
+    conv_delta = conv_op(delta)
+
+    # Plot results
+    plot = ift.Plot()
+    plot.add(signal, title='Signal')
+    plot.add(conv_signal, title='Signal Convolved')
+    plot.add(cac_signal, title='Signal, Conv, Adj-Conv')
+    plot.add(conv_delta, title='Kernel')
+    plot.output()
+
+
+# Healpix test
+nside = 64
+npix = 12 * nside * nside
+
+domain = ift.HPSpace(nside)
+
+# Define test signal (some point sources)
+signal_vals = np.zeros(npix, dtype=np.float64)
+for i in range(0, npix, npix//12 + 27):
+    signal_vals[i] = 500.
+
+signal = ift.from_global_data(domain, signal_vals)
+
+delta_vals = np.zeros(npix, dtype=np.float64)
+delta_vals[0] = 1.0
+delta = ift.from_global_data(domain, delta_vals)
+
+
+# Define kernel function
+def func(theta):
+    ct = np.cos(theta)
+    return 1. * np.logical_and(ct > 0.7, ct <= 0.8)
+
+
+convtest(signal, delta, func)
+
+domain = ift.RGSpace((100, 100))
+# Define test signal (some point sources)
+signal_vals = np.zeros(domain.shape, dtype=np.float64)
+signal_vals[35, 70] = 5000.
+signal_vals[45, 8] = 5000.
+signal = ift.from_global_data(domain, signal_vals)
+
+# Define delta signal, generate kernel image
+delta_vals = np.zeros(domain.shape, dtype=np.float64)
+delta_vals[0, 0] = 1.0
+delta = ift.from_global_data(domain, delta_vals)
+
+
+# Define kernel function
+def func(dist):
+    return 1. * np.logical_and(dist > 0.1, dist <= 0.2)
+
+
+convtest(signal, delta, func)

docs/source/ift.rst

@@ -37,7 +37,7 @@ NIFTy comes with reimplemented MAP and VI estimators.
 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::
 
@@ -49,7 +49,7 @@ and the measurement equation is linear in both signal and noise,
 
     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
 
@@ -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.
 
 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 reconstruction of only a Megapixel image would otherwithe require the storage and processing of matrices with sizes of several Terabytes.
+The advantage of implicit operators compared to explicit matrices is the reduced memory consumption;
+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.
 
 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
 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.
 
@@ -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.
 
 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`.
 
@@ -168,7 +168,7 @@ Maximum a Posteriori
 
 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::
 

nifty5/__init__.py

@@ -50,6 +50,7 @@ from .operators.value_inserter import ValueInserter
 from .operators.energy_operators import (
     EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood,
     BernoulliEnergy, StandardHamiltonian, AveragedEnergy)
+from .operators.convolution_operators import FuncConvolutionOperator
 
 from .probing import probe_with_posterior_samples, probe_diagonal, \
     StatCalculator
@@ -81,7 +82,8 @@ from .library.dynamic_operator import (dynamic_operator,
 from .library.light_cone_operator import LightConeOperator
 
 from .library.wiener_filter_curvature import WienerFilterCurvature
-from .library.correlated_fields import CorrelatedField, MfCorrelatedField
+from .library.correlated_fields import (CorrelatedField, MfCorrelatedField,
+                                        MfPartiallyCorrelatedField)
 from .library.adjust_variances import (make_adjust_variances_hamiltonian,
                                        do_adjust_variances)
 
@@ -93,5 +95,10 @@ from .logger import logger
 
 from .linearization import Linearization
 
+from . import internal_config
+_scheme = internal_config.parallelization_scheme()
+if _scheme == "Samples":
+    from .minimization.metric_gaussian_kl_mpi import MetricGaussianKL_MPI
+
 # We deliberately don't set __all__ here, because we don't want people to do a
 # "from nifty5 import *"; that would swamp the global namespace.

nifty5/dobj.py

@@ -15,11 +15,19 @@
 #
 # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
 
-try:
-    from mpi4py import MPI
-    if MPI.COMM_WORLD.Get_size() == 1:
-        from .data_objects.numpy_do import *
-    else:
-        from .data_objects.distributed_do import *
-except ImportError:
+from . import internal_config
+
+
+_scheme = internal_config.parallelization_scheme()
+
+if _scheme in ("Samples", "None"):
     from .data_objects.numpy_do import *
+else:
+    try:
+        from mpi4py import MPI
+        if MPI.COMM_WORLD.Get_size() == 1:
+            from .data_objects.numpy_do import *
+        else:
+            from .data_objects.distributed_do import *
+    except ImportError:
+        from .data_objects.numpy_do import *

nifty5/domains/lm_space.py

@@ -103,6 +103,36 @@ class LMSpace(StructuredDomain):
     def get_fft_smoothing_kernel_function(self, 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
     def lmax(self):
         """int : maximum allowed :math:`l`

nifty5/domains/rg_space.py

@@ -90,9 +90,7 @@ class RGSpace(StructuredDomain):
     def scalar_dvol(self):
         return self._dvol
 
-    def get_k_length_array(self):
-        if (not self.harmonic):
-            raise NotImplementedError
+    def _get_dist_array(self):
         ibegin = dobj.ibegin_from_shape(self._shape)
         res = np.arange(self.local_shape[0], dtype=np.float64) + ibegin[0]
         res = np.minimum(res, self.shape[0]-res)*self.distances[0]
@@ -106,6 +104,11 @@ class RGSpace(StructuredDomain):
             res = np.add.outer(res, tmp)
         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):
         if (not self.harmonic):
             raise NotImplementedError
@@ -146,6 +149,27 @@ class RGSpace(StructuredDomain):
             raise NotImplementedError
         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):
         """Returns a :class:`RGSpace` object representing the (position or
         harmonic) partner domain of `self`, depending on `self.harmonic`.

nifty5/field.py

@@ -636,6 +636,8 @@ class Field(object):
         return 0.5*(1.+self.tanh())
 
     def clip(self, min=None, max=None):
+        min = min.local_data if isinstance(min, Field) else min
+        max = max.local_data if isinstance(max, Field) else max
         return Field(self._domain, dobj.clip(self._val, min, max))
 
     def one_over(self):

nifty5/internal_config.py

+# 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 <http://www.gnu.org/licenses/>.
+#
+# Copyright(C) 2019 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+
+# Internal configuration switches, typically for experimental features.
+# Leave unchanged unless you know what you are doing!
+
+def parallelization_scheme():
+    return "Standard"
+    # return "Samples"
+    # return "None"

nifty5/library/correlated_fields.py

@@ -141,3 +141,54 @@ def MfCorrelatedField(target, amplitudes, name='xi'):
     vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
 
     return ht(vol*A*ducktape(hsp, None, name))
+
+
+def MfPartiallyCorrelatedField(target, energy_amplitude, name='xi_u'):
+    """Constructs an operator which turns white Gaussian excitation fields
+    into a correlated field defined on a DomainTuple with two entries.
+    On the first domain, a white correlation structure is assumed. On
+    the second domain, the correlation structures given by energy_amplitude
+    is used.
+
+    This operator may be used as a model for multi-frequency reconstructions
+    with correlation structure only in the energy direction.
+
+    Parameters
+    ----------
+    target : Domain, DomainTuple or tuple of Domain
+        Target of the operator. Must contain exactly two spaces.
+        It is assumed that the second space is the energy domain.
+    energy_amplitude: Operator
+        amplitude operator for the energy correlation structure
+    name : string
+        :class:`MultiField` key for xi-field.
+
+    Returns
+    -------
+    Operator
+        Correlated field
+
+    Notes
+    -----
+    In NIFTy, non-harmonic RGSpaces are by definition periodic. Therefore
+    the operator constructed by this method will output a correlated field
+    with *periodic* boundary conditions. If a non-periodic field is needed,
+    one needs to combine this operator with a :class:`FieldZeroPadder` or even
+    two (one for the energy and one for the spatial subdomain)
+    """
+    tgt = DomainTuple.make(target)
+    if len(tgt) != 2:
+        raise ValueError
+
+    h_space = DomainTuple.make([dom.get_default_codomain() for dom in tgt])
+    ht1 = HarmonicTransformOperator(h_space, target=tgt[0], space=0)
+    ht2 = HarmonicTransformOperator(ht1.target, target=tgt[1], space=1)
+    ht = ht2 @ ht1
+
+    p_space = energy_amplitude.target[0]
+    power_distributor = PowerDistributor(h_space[-1], p_space)
+    A = power_distributor(energy_amplitude)
+
+    dd = ContractionOperator(h_space, 0).adjoint
+
+    return ht((dd @ A)*ducktape(h_space, None, name))

nifty5/library/smooth_linear_amplitude.py

@@ -118,7 +118,7 @@ def CepstrumOperator(target, a, k0):
     return sym @ qht @ makeOp(cepstrum.sqrt())
 
 
-def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
+def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za=None, zq=None,
                 keys=['tau', 'phi', 'zero_mode']):
     '''Operator for parametrizing smooth amplitudes (square roots of power
     spectra).
@@ -168,10 +168,10 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
         LogRGSpace (see :class:`ExpTransform`).
     iv : float
         Prior standard deviation of y-intercept of power law.
-    za : float
+    za : float, optional
         The alpha-parameter of the inverse-gamma distribution.
         Setting the a seperate prior on the zeroGmode of the amplitude model.
-    zq : float
+    zq : float, optional
         The q-parameter of the inverse-gamma distribution.
 
     Returns
@@ -186,10 +186,17 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
 
     a, k0 = float(a), float(k0)
     sm, sv, im, iv = float(sm), float(sv), float(im), float(iv)
-    za, zq = float(za), float(zq)
     if sv <= 0 or iv <= 0:
         raise ValueError
 
+    if za is not None and zq is not None:
+        separate_zero_mode_prior = True
+        za, zq = float(za), float(zq)
+    else:
+        separate_zero_mode_prior = False
+        if za is not None or zq is not None:
+            raise ValueError("za and zq have to be given together")
+
     et = ExpTransform(target, n_pix)
     dom = et.domain[0]
 
@@ -208,14 +215,18 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, za, zq,
     # Combine linear and smooth component
     loglog_ampl = 0.5*(smooth + linear)
 
-    zero_mode = ValueInserter(et.target, (0,)*len(et.target.shape))
-
-    mask = np.ones(et.target.shape)
-    mask[(0,)*len(et.target.shape)] = 0.
-    mask = from_global_data(et.target, mask)
-    mask = DiagonalOperator(mask)
-    zero_mode = zero_mode @ InverseGammaOperator(
-        zero_mode.domain, za, zq) @ FieldAdapter(
-            zero_mode.domain, keys[2])
-    # Go from loglog-space to linear-linear-space
-    return (mask @ (et @ loglog_ampl.exp()) + zero_mode)
+    if not separate_zero_mode_prior:
+        # Go from loglog-space to linear-linear-space
+        return et @ loglog_ampl.exp()
+    else:
+        zero_mode = ValueInserter(et.target, (0,)*len(et.target.shape))
+
+        mask = np.ones(et.target.shape)
+        mask[(0,)*len(et.target.shape)] = 0.
+        mask = from_global_data(et.target, mask)
+        mask = DiagonalOperator(mask)
+        zero_mode = zero_mode @ InverseGammaOperator(
+            zero_mode.domain, za, zq) @ FieldAdapter(
+                zero_mode.domain, keys[2])
+
+        return mask @ (et @ loglog_ampl.exp()) + zero_mode

nifty5/minimization/metric_gaussian_kl_mpi.py

+# 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 <http://www.gnu.org/licenses/>.
+#
+# Copyright(C) 2013-2019 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
+
+from .. import utilities
+from ..linearization import Linearization
+from ..operators.energy_operators import StandardHamiltonian
+from .energy import Energy
+from mpi4py import MPI
+import numpy as np
+from ..field import Field
+from ..multi_field import MultiField
+
+
+_comm = MPI.COMM_WORLD
+ntask = _comm.Get_size()
+rank = _comm.Get_rank()
+master = (rank == 0)
+
+
+def _shareRange(nwork, nshares, myshare):
+    nbase = nwork//nshares
+    additional = nwork % nshares
+    lo = myshare*nbase + min(myshare, additional)
+    hi = lo + nbase + int(myshare < additional)
+    return lo, hi
+
+
+def np_allreduce_sum(arr):
+    arr = np.array(arr)
+    res = np.empty_like(arr)
+    _comm.Allreduce(arr, res, MPI.SUM)
+    return res
+
+
+def allreduce_sum_field(fld):
+    if isinstance(fld, Field):
+        return Field.from_local_data(fld.domain,
+                                     np_allreduce_sum(fld.local_data))
+    res = tuple(
+        Field.from_local_data(f.domain, np_allreduce_sum(f.local_data))
+        for f in fld.values())
+    return MultiField(fld.domain, res)
+
+
+class MetricGaussianKL_MPI(Energy):
+    def __init__(self, mean, hamiltonian, n_samples, constants=[],
+                 point_estimates=[], mirror_samples=False,
+                 _samples=None):
+        super(MetricGaussianKL_MPI, self).__init__(mean)
+
+        if not isinstance(hamiltonian, StandardHamiltonian):
+            raise TypeError
+        if hamiltonian.domain is not mean.domain:
+            raise ValueError
+        if not isinstance(n_samples, int):
+            raise TypeError
+        self._constants = list(constants)
+        self._point_estimates = list(point_estimates)
+        if not isinstance(mirror_samples, bool):
+            raise TypeError
+
+        self._hamiltonian = hamiltonian
+
+        if _samples is None:
+            lo, hi = _shareRange(n_samples, ntask, rank)
+            met = hamiltonian(Linearization.make_partial_var(
+                mean, point_estimates, True)).metric
+            _samples = []
+            for i in range(lo, hi):
+                np.random.seed(i)
+                _samples.append(met.draw_sample(from_inverse=True))
+            if mirror_samples:
+                _samples += [-s for s in _samples]
+                n_samples *= 2
+            _samples = tuple(_samples)
+        self._samples = _samples
+        self._n_samples = n_samples
+        self._lin = Linearization.make_partial_var(mean, constants)
+        v, g = None, None
+        if len(self._samples) == 0:  # hack if there are too many MPI tasks
+            tmp = self._hamiltonian(self._lin)
+            v = 0. * tmp.val.local_data
+            g = 0. * tmp.gradient
+        else:
+            for s in self._samples:
+                tmp = self._hamiltonian(self._lin+s)
+                if v is None:
+                    v = tmp.val.local_data.copy()
+                    g = tmp.gradient
+                else:
+                    v += tmp.val.local_data
+                    g = g + tmp.gradient
+        self._val = np_allreduce_sum(v)[()] / self._n_samples
+        self._grad = allreduce_sum_field(g) / self._n_samples
+        self._metric = None
+
+    def at(self, position):
+        return MetricGaussianKL_MPI(
+            position, self._hamiltonian, self._n_samples, self._constants,
+            self._point_estimates, _samples=self._samples)
+
+    @property
+    def value(self):
+        return self._val
+
+    @property
+    def gradient(self):
+        return self._grad
+
+    def _get_metric(self):
+        lin = self._lin.with_want_metric()
+        if self._metric is None:
+            if len(self._samples) == 0:  # hack if there are too many MPI tasks
+                self._metric = self._hamiltonian(lin).metric.scale(0.)
+            else:
+                mymap = map(lambda v: self._hamiltonian(lin+v).metric,
+                            self._samples)
+                self._metric = utilities.my_sum(mymap)
+                self._metric = self._metric.scale(1./self._n_samples)
+