merge

demos/getting_started_3.py

@@ -29,10 +29,11 @@ def get_random_LOS(n_los):
 if __name__ == '__main__':
     # FIXME description of the tutorial
     np.random.seed(42)
+    np.seterr(all='raise')
     position_space = ift.RGSpace([128, 128])
 
     # Setting up an amplitude model
-    A = ift.AmplitudeModel(position_space, 16, 1, 10, -4., 1, 0., 1.)
+    A = ift.AmplitudeModel(position_space, 64, 3, 0.4, -4., 1, 1., 1.)
 
     # Building the model for a correlated signal
     harmonic_space = position_space.get_default_codomain()
@@ -40,10 +41,11 @@ if __name__ == '__main__':
     power_space = A.target[0]
     power_distributor = ift.PowerDistributor(harmonic_space, power_space)
 
-    correlated_field = ht(
-        power_distributor(A)*ift.FieldAdapter(harmonic_space, "xi"))
+    vol = harmonic_space.scalar_dvol
+    vol = ift.ScalingOperator(vol ** (-0.5),harmonic_space)
+    correlated_field = ht(vol(power_distributor(A))*ift.FieldAdapter(harmonic_space, "xi"))
     # alternatively to the block above one can do:
-    # correlated_field = ift.CorrelatedField(position_space, A)
+    #correlated_field = ift.CorrelatedField(position_space, A)
 
     # apply some nonlinearity
     signal = ift.positive_tanh(correlated_field)

nifty5/__init__.py

@@ -34,6 +34,7 @@ from .operators.inversion_enabler import InversionEnabler
 from .operators.laplace_operator import LaplaceOperator
 from .operators.linear_operator import LinearOperator
 from .operators.mask_operator import MaskOperator
+from .operators.offset_operator import OffsetOperator
 from .operators.qht_operator import QHTOperator
 from .operators.regridding_operator import RegriddingOperator
 from .operators.sampling_enabler import SamplingEnabler
@@ -73,7 +74,7 @@ from .minimization.kl_energy import KL_Energy
 from .sugar import *
 from .plot import Plot
 
-from .library.amplitude_model import AmplitudeModel, SmoothAmplitudeModel
+from .library.amplitude_model import AmplitudeModel
 from .library.inverse_gamma_model import InverseGammaModel
 from .library.los_response import LOSResponse
 

nifty5/domains/log_rg_space.py

@@ -91,22 +91,22 @@ class LogRGSpace(StructuredDomain):
         return LogRGSpace(self.shape, codomain_bindistances, self._t_0, True)
 
     def get_k_length_array(self):
-        ib = dobj.ibegin_from_shape(self._shape)
-        res = np.arange(self.local_shape[0], dtype=np.float64) + ib[0]
-        res = np.minimum(res, self.shape[0] - res)*self.bindistances[0]
-        if len(self.shape) == 1:
-            return Field.from_local_data(self, res)
-        res *= res
-        for i in range(1, len(self.shape)):
-            tmp = np.arange(self.local_shape[i], dtype=np.float64) + ib[i]
-            tmp = np.minimum(tmp, self.shape[i] - tmp)*self.bindistances[i]
-            tmp *= tmp
-            res = np.add.outer(res, tmp)
-        return Field.from_local_data(self, np.sqrt(res))
-
-    def get_expk_length_array(self):
-        # FIXME This is a hack! Only for plotting. Seems not to be the final version.
-        out = exp(self.get_k_length_array()).to_global_data_rw()
-        out[1:] = out[:-1]
-        out[0] = 0
-        return Field.from_global_data(self, out)
+        if not self.harmonic:
+            raise NotImplementedError
+        ks = self.get_k_array()
+        return Field.from_global_data(self, np.linalg.norm(ks, axis=0))
+
+    def get_k_array(self):
+        ndim = len(self.shape)
+        k_array = np.zeros((ndim,) + self.shape)
+        dist = self.bindistances
+        for i in range(ndim):
+            ks = np.zeros(self.shape[i])
+            ks[1:] = np.minimum(self.shape[i] - 1 - np.arange(self.shape[i]-1), np.arange(self.shape[i]-1)) * dist[i]
+            if self.harmonic:
+                ks[0] = np.nan
+            else:
+                ks[0] = -np.inf
+                ks[1:] += self.t_0[i]
+            k_array[i] += ks.reshape((1,)*i + (self.shape[i],) + (1,)*(ndim-i-1))
+        return k_array

nifty5/library/amplitude_model.py

@@ -23,13 +23,11 @@ import numpy as np
 from ..compat import *
 from ..domains.power_space import PowerSpace
 from ..field import Field
-from ..multi_domain import MultiDomain
-from ..operators.operator import Operator
 from ..sugar import makeOp, sqrt
 
 
 def _ceps_kernel(dof_space, k, a, k0):
-    return a**2/(1+(k/(k0*dof_space.bindistances[0]))**2)**2
+    return a**2/(1 + (k/k0)**2)**2
 
 
 def create_cepstrum_amplitude_field(domain, cepstrum):
@@ -45,22 +43,12 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
     """
 
     dim = len(domain.shape)
-    dist = domain.bindistances
     shape = domain.shape
 
-    # Prepare q_array
-    q_array = np.zeros((dim,) + shape)
-    if dim == 1:
-        ks = domain.get_k_length_array().to_global_data()
-        q_array = np.array([ks])
-    else:
-        for i in range(dim):
-            ks = np.minimum(shape[i] - np.arange(shape[i]) +
-                            1, np.arange(shape[i])) * dist[i]
-            q_array[i] += ks.reshape((1,)*i + (shape[i],) + (1,)*(dim-i-1))
+    q_array = domain.get_k_array()
 
     # Fill cepstrum field (all non-zero modes)
-    no_zero_modes = (slice(1, None),) * dim
+    no_zero_modes = (slice(1, None),)*dim
     ks = q_array[(slice(None),) + no_zero_modes]
     cepstrum_field = np.zeros(shape)
     cepstrum_field[no_zero_modes] = cepstrum(ks)
@@ -78,81 +66,55 @@ def create_cepstrum_amplitude_field(domain, cepstrum):
     return Field.from_global_data(domain, cepstrum_field)
 
 
-class AmplitudeModel(Operator):
+def CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode=True):
     '''
-    Computes a smooth power spectrum.
-    Output lives in PowerSpace.
-
     Parameters
     ----------
-
-    Npixdof : #pix in dof_space
-
     ceps_a, ceps_k0 : Smoothness parameters in ceps_kernel
                         eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
                         a = ceps_a,  k0 = ceps_k0
+    '''
+
+    from ..operators.qht_operator import QHTOperator
+    from ..operators.symmetrizing_operator import SymmetrizingOperator
+    qht = QHTOperator(target=logk_space)
+    dof_space = qht.domain[0]
+    sym = SymmetrizingOperator(logk_space)
+    kern = lambda k: _ceps_kernel(dof_space, k, ceps_a, ceps_k)
+    cepstrum = create_cepstrum_amplitude_field(dof_space, kern)
+    res = sym(qht(makeOp(sqrt(cepstrum))))
+    if not zero_mode:
+        shp = res.target.shape
+        mask = np.ones(shp)
+        mask[(0,)*len(shp)] = 0.
+        mask = makeOp(Field.from_global_data(res.target, mask))
+        res = mask(res)
+    return res
+
+
+def SlopeModel(logk_space, sm, sv, im, iv):
+    '''
+    Parameters
+    ----------
 
     sm, sv : slope_mean = expected exponent of power law (e.g. -4),
                 slope_variance (default=1)
 
-    im, iv : y-intercept_mean, y-intercept_variance  of power_slope
+    im, iv : y-intercept_mean, y-intercept_std  of power_slope
     '''
-    def __init__(self, s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
-                 keys=['tau', 'phi']):
-        from ..operators.exp_transform import ExpTransform
-        from ..operators.qht_operator import QHTOperator
-        from ..operators.slope_operator import SlopeOperator
-        from ..operators.symmetrizing_operator import SymmetrizingOperator
-
-        h_space = s_space.get_default_codomain()
-        self._exp_transform = ExpTransform(PowerSpace(h_space), Npixdof)
-        logk_space = self._exp_transform.domain[0]
-        qht = QHTOperator(target=logk_space)
-        dof_space = qht.domain[0]
-        sym = SymmetrizingOperator(logk_space)
-
-        phi_mean = np.array([sm, im])
-        phi_sig = np.array([sv, iv])
-
-        self._slope = SlopeOperator(logk_space, phi_sig)
-        self._norm_phi_mean = Field.from_global_data(self._slope.domain,
-                                                     phi_mean/phi_sig)
-
-        self._domain = MultiDomain.make({keys[0]: dof_space,
-                                         keys[1]: self._slope.domain})
-        self._target = self._exp_transform.target
-
-        kern = lambda k: _ceps_kernel(dof_space, k, ceps_a, ceps_k)
-        cepstrum = create_cepstrum_amplitude_field(dof_space, kern)
-
-        self._smooth_op = sym(qht(makeOp(sqrt(cepstrum))))
-        self._keys = tuple(keys)
-
-        self._qht = qht
-        self._ceps = makeOp(sqrt(cepstrum))
-
-    def apply(self, x):
-        self._check_input(x)
-        smooth_spec = self._smooth_op(x[self._keys[0]])
-        phi = x[self._keys[1]] + self._norm_phi_mean
-        linear_spec = self._slope(phi)
-        loglog_spec = smooth_spec + linear_spec
-        return self._exp_transform((0.5*loglog_spec).exp())
-
-    @property
-    def qht(self):
-        return self._qht
-
-    @property
-    def ceps(self):
-        return self._ceps
-
-    @property
-    def norm_phi_mean(self):
-        return self._norm_phi_mean
-
-
-class SmoothAmplitudeModel(Operator):
+
+    from ..operators.slope_operator import SlopeOperator
+    from ..operators.offset_operator import OffsetOperator
+    phi_mean = np.array([sm, im + sm*logk_space.t_0[0]])
+    phi_sig = np.array([sv, iv])
+    slope = SlopeOperator(logk_space)
+    phi_mean = Field.from_global_data(slope.domain, phi_mean)
+    phi_sig = Field.from_global_data(slope.domain, phi_sig)
+    return slope(OffsetOperator(phi_mean)(makeOp(phi_sig)))
+
+
+def AmplitudeModel(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
+                   keys=['tau', 'phi'], zero_mode=True):
     '''
     Computes a smooth power spectrum.
     Output lives in PowerSpace.
@@ -165,39 +127,26 @@ class SmoothAmplitudeModel(Operator):
     ceps_a, ceps_k0 : Smoothness parameters in ceps_kernel
                         eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
                         a = ceps_a,  k0 = ceps_k0
-    '''
-    def __init__(self, s_space, Npixdof, ceps_a, ceps_k, name='tau'):
-        from ..operators.exp_transform import ExpTransform
-        from ..operators.qht_operator import QHTOperator
 
-        h_space = s_space.get_default_codomain()
-        self._exp_transform = ExpTransform(PowerSpace(h_space), Npixdof)
-        logk_space = self._exp_transform.domain[0]
-        qht = QHTOperator(target=logk_space)
-        dof_space = qht.domain[0]
-
-        self._domain = MultiDomain.make({name: dof_space})
-        self._target = self._exp_transform.target
+    sm, sv : slope_mean = expected exponent of power law (e.g. -4),
+                slope_variance (default=1)
 
-        kern = lambda k: _ceps_kernel(dof_space, k, ceps_a, ceps_k)
-        cepstrum = create_cepstrum_amplitude_field(dof_space, kern)
+    im, iv : y-intercept_mean, y-intercept_variance  of power_slope
+    '''
 
-        self._smooth_op = qht(makeOp(sqrt(cepstrum)))
-        self._key = name
+    from ..operators.exp_transform import ExpTransform
+    from ..operators.simple_linear_operators import FieldAdapter
+    from ..operators.scaling_operator import ScalingOperator
 
-        self._qht = qht
-        self._ceps = makeOp(sqrt(cepstrum))
+    h_space = s_space.get_default_codomain()
+    et = ExpTransform(PowerSpace(h_space), Npixdof)
+    logk_space = et.domain[0]
 
-    def apply(self, x):
-        self._check_input(x)
-        smooth_spec = self._smooth_op(x[self._key])
-        loglog_spec = smooth_spec
-        return self._exp_transform((0.5*loglog_spec).exp())
+    smooth = CepstrumOperator(logk_space, ceps_a, ceps_k, zero_mode)
+    linear = SlopeModel(logk_space, sm, sv, im, iv)
 
-    @property
-    def qht(self):
-        return self._qht
+    fa_smooth = FieldAdapter(smooth.domain, keys[0])
+    fa_linear = FieldAdapter(linear.domain, keys[1])
 
-    @property
-    def ceps(self):
-        return self._ceps
+    fac = ScalingOperator(0.5, smooth.target)
+    return et((fac(smooth(fa_smooth) + linear(fa_linear))).exp())

nifty5/library/correlated_fields.py

@@ -25,6 +25,7 @@ from ..operators.contraction_operator import ContractionOperator
 from ..operators.distributors import PowerDistributor
 from ..operators.harmonic_operators import HarmonicTransformOperator
 from ..operators.simple_linear_operators import FieldAdapter
+from ..operators.scaling_operator import ScalingOperator
 
 
 def CorrelatedField(s_space, amplitude_model, name='xi'):
@@ -45,7 +46,9 @@ def CorrelatedField(s_space, amplitude_model, name='xi'):
     p_space = amplitude_model.target[0]
     power_distributor = PowerDistributor(h_space, p_space)
     A = power_distributor(amplitude_model)
-    return ht(A*FieldAdapter(h_space, name))
+    vol = h_space.scalar_dvol
+    vol = ScalingOperator(vol**(-0.5), h_space)
+    return ht(vol(A)*FieldAdapter(h_space, name))
 
 
 def MfCorrelatedField(s_space_spatial, s_space_energy, amplitude_model_spatial,

nifty5/operators/offset_operator.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-2018 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
+# and financially supported by the Studienstiftung des deutschen Volkes.
+
+from __future__ import absolute_import, division, print_function
+
+from ..compat import *
+from .operator import Operator
+
+
+class OffsetOperator(Operator):
+    def __init__(self, field):
+        self._field = field
+        self._domain = self._target = field.domain
+
+    def apply(self, x):
+        self._check_input(x)
+        return x + self._field

nifty5/operators/simple_linear_operators.py

@@ -24,7 +24,6 @@ from ..domains.unstructured_domain import UnstructuredDomain
 from ..field import Field
 from ..multi_domain import MultiDomain
 from ..multi_field import MultiField
-from ..sugar import full
 from .endomorphic_operator import EndomorphicOperator
 from .linear_operator import LinearOperator
 

nifty5/operators/slope_operator.py

@@ -46,43 +46,27 @@ class SlopeOperator(LinearOperator):
     sigmas : np.array, shape=(2,)
         The slope variance and the y-intercept variance.
     """
-    def __init__(self, target, sigmas):
+
+    def __init__(self, target):
         if not isinstance(target, LogRGSpace):
             raise TypeError
+        if len(target.shape) != 1:
+            raise ValueError("Slope Operator only works for ndim == 1")
         self._domain = DomainTuple.make(UnstructuredDomain((2,)))
         self._target = DomainTuple.make(target)
         self._capability = self.TIMES | self.ADJOINT_TIMES
-
-        self._sigmas = sigmas
-        self.ndim = len(self.target[0].shape)
-        self.pos = np.zeros((self.ndim,) + self.target[0].shape)
-
-        if self.ndim == 1:
-            self.pos[0] = self.target[0].get_k_length_array().to_global_data()
-        else:
-            shape = self.target[0].shape
-            for i in range(self.ndim):
-                rng = np.arange(target.shape[i])
-                tmp = np.minimum(
-                    rng, target.shape[i]+1-rng) * target.bindistances[i]
-                self.pos[i] += tmp.reshape(
-                    (1,)*i + (shape[i],) + (1,)*(self.ndim-i-1))
+        pos = self.target[0].get_k_array() - self.target[0].t_0[0]
+        self._pos = pos[0, 1:]
 
     def apply(self, x, mode):
         self._check_input(x, mode)
-
-        # Times
+        inp = x.to_global_data()
         if mode == self.TIMES:
-            inp = x.to_global_data()
-            res = self._sigmas[-1] * inp[-1]
-            for i in range(self.ndim):
-                res += self._sigmas[i] * inp[i] * self.pos[i]
-            return Field.from_global_data(self.target, res)
-
-        # Adjoint times
-        res = np.zeros(self.domain[0].shape, dtype=x.dtype)
-        xglob = x.to_global_data()
-        res[-1] = np.sum(xglob) * self._sigmas[-1]
-        for i in range(self.ndim):
-            res[i] = np.sum(self.pos[i] * xglob) * self._sigmas[i]
-        return Field.from_global_data(self.domain, res)
+            res = np.empty(self.target.shape, dtype=x.dtype)
+            res[0] = 0
+            res[1:] = inp[1] + inp[0]*self._pos
+        else:
+            res = np.array(
+                [np.sum(self._pos*inp[1:]),
+                 np.sum(inp[1:])], dtype=x.dtype)
+        return Field.from_global_data(self._tgt(mode), res)

test/test_operators/test_adjoint.py

@@ -73,8 +73,7 @@ class Consistency_Tests(unittest.TestCase):
     def testSlopeOperator(self, args, dtype):
         tmp = ift.ExpTransform(ift.PowerSpace(args[0]), args[1], args[2])
         tgt = tmp.domain[0]
-        sig = np.array([0.3, 0.13])
-        op = ift.SlopeOperator(tgt, sig)
+        op = ift.SlopeOperator(tgt)
         ift.extra.consistency_check(op, dtype, dtype)
 
     @expand(product(_h_spaces + _p_spaces + _pow_spaces,