Move files around

demos/newamplitudes/correlated_fields_normalized_amplitude.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.
+
+import nifty5 as ift
+
+if __name__ == '__main__':
+    import numpy as np
+    np.random.seed(42)
+    nspecs = 2
+    a = []
+    spaces = []
+    for _ in range(nspecs):
+        ndim = 2
+        sspace = ift.RGSpace(
+            np.linspace(16, 20, num=ndim).astype(np.int),
+            np.linspace(2.3, 7.99, num=ndim))
+        hspace = sspace.get_default_codomain()
+        spaces.append(sspace)
+        target = ift.PowerSpace(hspace)
+        a.append(ift.NormalizedAmplitude(target, 16, 1, 1, -3, 1, 0, 1, 0, 1))
+    tgt = ift.makeDomain(tuple(spaces))
+    op = ift.CorrelatedFieldNormAmplitude(tgt, a, 0., 1.)
+    fld = ift.from_random('normal', op.domain)
+    ift.extra.check_jacobian_consistency(op, fld)

demos/newamplitudes/normalized_amplitude.py

+import numpy as np
+
+import nifty5 as ift
+
+if __name__ == '__main__':
+    np.random.seed(42)
+    ndim = 1
+    sspace = ift.RGSpace(
+        np.linspace(16, 20, num=ndim).astype(np.int),
+        np.linspace(2.3, 7.99, num=ndim))
+    hspace = sspace.get_default_codomain()
+    target = ift.PowerSpace(hspace)
+    op = ift.NormalizedAmplitude(target, 16, 1, 1, -3, 1, 0, 1, 0, 1)
+    if not isinstance(op.target[0], ift.PowerSpace):
+        raise RuntimeError
+    fld = ift.from_random('normal', op.domain)
+    ift.extra.check_jacobian_consistency(op, fld)
+
+    pd = ift.PowerDistributor(hspace, power_space=target)
+    ht = ift.HartleyOperator(pd.target)
+    if np.testing.assert_allclose((pd @ op)(fld).integrate(), 1):
+        raise RuntimeError
+    if np.testing.assert_allclose(
+        (ht @ pd @ op)(fld).to_global_data()[ndim*(0,)], 1):
+        raise RuntimeError

demos/newamplitudes/wiener_process_integrated_amplitude.py

+import nifty5 as ift
+import numpy as np
+
+
+def mfplot(fld, A1, A2, name):
+    p = ift.Plot()
+    dom = fld.domain
+    mi, ma = np.min(fld.val), np.max(fld.val)
+    for ii in range(dom.shape[-1]):
+        extr = ift.DomainTupleFieldInserter(op.target, 1, (ii,)).adjoint
+        p.add(extr(fld), zmin=mi, zmax=ma)
+    p.add(A1)
+    p.add(A2)
+    p.add(fld)
+    p.output(name=name, xsize=14, ysize=14)
+
+
+if __name__ == '__main__':
+    np.random.seed(42)
+    sspace = ift.RGSpace((512, 512))
+    hspace = sspace.get_default_codomain()
+    target = ift.PowerSpace(hspace,
+                            ift.PowerSpace.useful_binbounds(hspace, True))
+    # A = Amplitude(target, [5, -3, 1, 0], [1, 1, 1, 0.5],
+    #               ['rest', 'smooth', 'wienersigma'])
+    # fld = ift.from_random('normal', A.domain)
+    # ift.extra.check_jacobian_consistency(A, fld)
+    # op = CorrelatedFieldNormAmplitude(sspace, A, 3, 2)
+    # plts = []
+    # p = ift.Plot()
+    # for _ in range(25):
+    #     fld = ift.from_random('normal', op.domain)
+    #     plts.append(A.force(fld))
+    #     p.add(op(fld))
+    # p.output(name='debug1.png', xsize=20, ysize=20)
+    # ift.single_plot(plts, name='debug.png')
+
+    # Multifrequency
+    sspace1 = ift.RGSpace((512, 512))
+    hspace1 = sspace1.get_default_codomain()
+    target1 = ift.PowerSpace(hspace1,
+                             ift.PowerSpace.useful_binbounds(hspace1, True))
+    A1 = ift.WPAmplitude(target1, [0, -2, 0], [1E-5, 1, 1], 1, 0.99,
+                   ['rest1', 'smooth1', 'wienersigma1'])
+    sspace2 = ift.RGSpace((20,), distances=(1e7))
+    hspace2 = sspace2.get_default_codomain()
+    target2 = ift.PowerSpace(hspace2,
+                             ift.PowerSpace.useful_binbounds(hspace2, True))
+    A2 = ift.WPAmplitude(target2, [0, -2, 0], [1E-5, 1, 1], 1, 0.99,
+                   ['rest2', 'smooth2', 'wienersigma2'])
+    op = ift.CorrelatedFieldNormAmplitude((sspace1, sspace2), (A1, A2), 3, 2)
+    for jj in range(10):
+        fld = ift.from_random('normal', op.domain)
+        mfplot(op(fld), A1.force(fld), A2.force(fld), 'debug{}.png'.format(jj))

nifty5/__init__.py

@@ -77,7 +77,8 @@ from .sugar import *
 from .plot import Plot
 
 from .library.smooth_linear_amplitude import (
-    SLAmplitude, LinearSLAmplitude, CepstrumOperator)
+    SLAmplitude, LinearSLAmplitude, CepstrumOperator, WPAmplitude)
+from .library.normalized_amplitude import NormalizedAmplitude
 from .library.inverse_gamma_operator import InverseGammaOperator
 from .library.los_response import LOSResponse
 from .library.dynamic_operator import (dynamic_operator,
@@ -85,7 +86,7 @@ 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, CorrelatedFieldNormAmplitude
 from .library.adjust_variances import (make_adjust_variances_hamiltonian,
                                        do_adjust_variances)
 from .library.gridder import GridderMaker

nifty5/library/correlated_fields.py

@@ -19,10 +19,13 @@ from functools import reduce
 from operator import mul
 
 from ..domain_tuple import DomainTuple
+from ..operators.adder import Adder
 from ..operators.contraction_operator import ContractionOperator
 from ..operators.distributors import PowerDistributor
 from ..operators.harmonic_operators import HarmonicTransformOperator
-from ..operators.simple_linear_operators import ducktape
+from ..operators.operator import Operator
+from ..operators.simple_linear_operators import VdotOperator, ducktape
+from ..sugar import full
 
 
 def CorrelatedField(target, amplitude_operator, name='xi', codomain=None):
@@ -132,3 +135,81 @@ def MfCorrelatedField(target, amplitudes, name='xi'):
     # For `vol` see comment in `CorrelatedField`
     vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
     return ht(vol*A*ducktape(hsp, None, name))
+
+
+def CorrelatedFieldNormAmplitude(target,
+                                 amplitudes,
+                                 stdmean,
+                                 stdstd,
+                                 names=['xi', 'std']):
+    """Constructs an operator which turns white Gaussian excitation fields
+    into a correlated field defined on a DomainTuple with n entries and n
+    separate correlation structures.
+
+    This operator may be used as a model for multi-frequency reconstructions
+    with a correlation structure in both spatial and energy direction.
+
+    Parameters
+    ----------
+    target : Domain, DomainTuple or tuple of Domain
+        Target of the operator. Must contain exactly n spaces.
+    amplitudes: Opertor, iterable of Operator
+        Amplitude operator if n = 1 or list of n amplitude operators.
+    stdmean : float
+        Prior mean of the overall standart deviation.
+    stdstd : float
+        Prior standart deviation of the overall standart deviation.
+    names : iterable of string
+        :class:`MultiField` keys for xi-field and std-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)
+    """
+
+    amps = [
+        amplitudes,
+    ] if isinstance(amplitudes, Operator) else amplitudes
+
+    tgt = DomainTuple.make(target)
+    if len(tgt) != len(amps):
+        raise ValueError
+    stdmean, stdstd = float(stdmean), float(stdstd)
+    if stdstd <= 0:
+        raise ValueError
+
+    psp = [aa.target[0] for aa in amps]
+    hsp = DomainTuple.make([tt.get_default_codomain() for tt in tgt])
+
+    ht = HarmonicTransformOperator(hsp, target=tgt[0], space=0)
+    pd = PowerDistributor(hsp, psp[0], 0)
+
+    for i in range(1, len(amps)):
+        ht = HarmonicTransformOperator(ht.target, target=tgt[i], space=i) @ ht
+        pd = pd @ PowerDistributor(pd.domain, psp[i], space=i)
+
+    spaces = tuple(range(len(amps)))
+
+    a = ContractionOperator(pd.domain, spaces[1:]).adjoint @ amps[0]
+    for i in range(1, len(amps)):
+        a = a*(ContractionOperator(pd.domain, spaces[:i] + spaces[
+            (i + 1):]).adjoint @ amps[i])
+
+    expander = VdotOperator(full(a.target, 1.)).adjoint
+
+    Std = stdstd*ducktape(expander.domain, None, names[1])
+    Std = expander @ (Adder(full(expander.domain, stdmean)) @ Std).exp()
+
+    A = pd @ (Std*a)
+    # For `vol` see comment in `CorrelatedField`
+    vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
+    return ht(vol*A*ducktape(hsp, None, names[0]))

nifty5/library/correlated_fields_normalized_amplitude.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.
-
-import nifty5 as ift
-from functools import reduce
-from operator import mul
-
-
-def CorrelatedFieldNormAmplitude(target,
-                                 amplitudes,
-                                 stdmean,
-                                 stdstd,
-                                 names=['xi', 'std']):
-    """Constructs an operator which turns white Gaussian excitation fields
-    into a correlated field defined on a DomainTuple with n entries and n
-    separate correlation structures.
-
-    This operator may be used as a model for multi-frequency reconstructions
-    with a correlation structure in both spatial and energy direction.
-
-    Parameters
-    ----------
-    target : Domain, DomainTuple or tuple of Domain
-        Target of the operator. Must contain exactly n spaces.
-    amplitudes: Opertor, iterable of Operator
-        Amplitude operator if n = 1 or list of n amplitude operators.
-    stdmean : float
-        Prior mean of the overall standart deviation.
-    stdstd : float
-        Prior standart deviation of the overall standart deviation.
-    names : iterable of string
-        :class:`MultiField` keys for xi-field and std-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)
-    """
-
-    amps = [
-        amplitudes,
-    ] if isinstance(amplitudes, ift.Operator) else amplitudes
-
-    tgt = ift.DomainTuple.make(target)
-    if len(tgt) != len(amps):
-        raise ValueError
-    stdmean, stdstd = float(stdmean), float(stdstd)
-    if stdstd <= 0:
-        raise ValueError
-
-    psp = [aa.target[0] for aa in amps]
-    hsp = ift.DomainTuple.make([tt.get_default_codomain() for tt in tgt])
-
-    ht = ift.HarmonicTransformOperator(hsp, target=tgt[0], space=0)
-    pd = ift.PowerDistributor(hsp, psp[0], 0)
-
-    for i in range(1, len(amps)):
-        ht = ift.HarmonicTransformOperator(ht.target, target=tgt[i],
-                                           space=i) @ ht
-        pd = pd @ ift.PowerDistributor(pd.domain, psp[i], space=i)
-
-    spaces = tuple(range(len(amps)))
-
-    a = ift.ContractionOperator(pd.domain, spaces[1:]).adjoint @ amps[0]
-    for i in range(1, len(amps)):
-        a = a*(ift.ContractionOperator(pd.domain, spaces[:i] + spaces[
-            (i + 1):]).adjoint @ amps[i])
-
-    expander = ift.VdotOperator(ift.full(a.target, 1.)).adjoint
-
-    Std = stdstd*ift.ducktape(expander.domain, None, names[1])
-    Std = expander @ (
-        ift.Adder(ift.full(expander.domain, stdmean)) @ Std).exp()
-
-    A = pd @ (Std*a)
-    # For `vol` see comment in `CorrelatedField`
-    vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
-    return ht(vol*A*ift.ducktape(hsp, None, names[0]))
-
-
-if __name__ == '__main__':
-    import numpy as np
-    from normalized_amplitude import NormalizedAmplitude
-    np.random.seed(42)
-    nspecs = 2
-    a = []
-    spaces = []
-    for _ in range(nspecs):
-        ndim = 2
-        sspace = ift.RGSpace(
-            np.linspace(16, 20, num=ndim).astype(np.int),
-            np.linspace(2.3, 7.99, num=ndim))
-        hspace = sspace.get_default_codomain()
-        spaces.append(sspace)
-        target = ift.PowerSpace(hspace)
-        a.append(NormalizedAmplitude(target, 16, 1, 1, -3, 1, 0, 1, 0, 1))
-    tgt = ift.makeDomain(tuple(spaces))
-    op = CorrelatedFieldNormAmplitude(tgt, a, 0., 1.)
-    fld = ift.from_random('normal', op.domain)
-    ift.extra.check_jacobian_consistency(op, fld)

nifty5/library/normalized_amplitude.py

@@ -189,24 +189,3 @@ class SpecialSum(ift.EndomorphicOperator):
         return ift.full(self._tgt(mode), x.sum())
 
 
-if __name__ == '__main__':
-    np.random.seed(42)
-    ndim = 1
-    sspace = ift.RGSpace(
-        np.linspace(16, 20, num=ndim).astype(np.int),
-        np.linspace(2.3, 7.99, num=ndim))
-    hspace = sspace.get_default_codomain()
-    target = ift.PowerSpace(hspace)
-    op = NormalizedAmplitude(target, 16, 1, 1, -3, 1, 0, 1, 0, 1)
-    if not isinstance(op.target[0], ift.PowerSpace):
-        raise RuntimeError
-    fld = ift.from_random('normal', op.domain)
-    ift.extra.check_jacobian_consistency(op, fld)
-
-    pd = ift.PowerDistributor(hspace, power_space=target)
-    ht = ift.HartleyOperator(pd.target)
-    if np.testing.assert_allclose((pd @ op)(fld).integrate(), 1):
-        raise RuntimeError
-    if np.testing.assert_allclose(
-        (ht @ pd @ op)(fld).to_global_data()[ndim*(0,)], 1):
-        raise RuntimeError

nifty5/library/smooth_linear_amplitude.py

@@ -19,13 +19,19 @@ import numpy as np
 
 from ..domain_tuple import DomainTuple
 from ..domains.power_space import PowerSpace
+from ..domains.unstructured_domain import UnstructuredDomain
 from ..field import Field
 from ..operators.adder import Adder
+from ..operators.distributors import PowerDistributor
+from ..operators.endomorphic_operator import EndomorphicOperator
 from ..operators.exp_transform import ExpTransform
+from ..operators.linear_operator import LinearOperator
+from ..operators.operator import Operator
 from ..operators.qht_operator import QHTOperator
+from ..operators.simple_linear_operators import VdotOperator, ducktape
 from ..operators.slope_operator import SlopeOperator
 from ..operators.symmetrizing_operator import SymmetrizingOperator
-from ..sugar import makeOp
+from ..sugar import from_global_data, full, makeDomain, makeOp
 
 
 def _ceps_kernel(k, a, k0):
@@ -209,3 +215,114 @@ def LinearSLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv,
 
     # Go from loglog-space to linear-linear-space
     return et @ loglog_ampl
+
+
+class _TwoLogIntegrations(LinearOperator):
+    def __init__(self, target):
+        self._target = makeDomain(target)
+        self._domain = makeDomain(
+            UnstructuredDomain(self.target.shape[0] - 2))
+        self._capability = self.TIMES | self.ADJOINT_TIMES
+        if not isinstance(self._target[0], PowerSpace):
+            raise TypeError
+        logk_lengths = np.log(self._target[0].k_lengths[1:])
+        self._logvol = logk_lengths[1:] - logk_lengths[:-1]
+
+    def apply(self, x, mode):
+        self._check_input(x, mode)
+        if mode == self.TIMES:
+            x = x.to_global_data()
+            res = np.empty(self._target.shape)
+            res[0] = 0
+            res[1] = 0
+            res[2:] = np.cumsum(x*self._logvol)
+            res[2:] = np.cumsum(res[2:]*self._logvol)
+            return from_global_data(self._target, res)
+        else:
+            x = x.to_global_data()
+            res = np.empty(self._target.shape)
+            res[2:] = np.cumsum(x[2:][::-1])[::-1]*self._logvol
+            res[2:] = np.cumsum(res[2:][::-1])[::-1]*self._logvol
+            return from_global_data(self._domain, res[2:])
+
+
+class _Rest(LinearOperator):
+    def __init__(self, target):
+        self._target = makeDomain(target)
+        self._domain = makeDomain(UnstructuredDomain(3))
+        self._logk_lengths = np.log(self._target[0].k_lengths[1:])
+        self._logk_lengths -= self._logk_lengths[0]
+        self._capability = self.TIMES | self.ADJOINT_TIMES
+
+    def apply(self, x, mode):
+        self._check_input(x, mode)
+        x = x.to_global_data()
+        res = np.empty(self._tgt(mode).shape)
+        if mode == self.TIMES:
+            res[0] = x[0]
+            res[1:] = x[1]*self._logk_lengths + x[2]
+        else:
+            res[0] = x[0]
+            res[1] = np.vdot(self._logk_lengths, x[1:])
+            res[2] = np.sum(x[1:])
+        return from_global_data(self._tgt(mode), res)
+
+
+def LogIntegratedWienerProcess(target, means, stddevs, wienersigmastddev,
+                               wienersigmaprob, keys):
+    # means and stddevs: zm, slope, yintercept
+    # keys: rest smooth wienersigma
+    if not (len(means) == 3 and len(stddevs) == 3 and len(keys) == 3):
+        raise ValueError
+    means = np.array(means)
+    stddevs = np.array(stddevs)
+    # FIXME More checks
+    rest = _Rest(target)
+    restmeans = from_global_data(rest.domain, means)
+    reststddevs = from_global_data(rest.domain, stddevs)
+    rest = rest @ Adder(restmeans) @ makeOp(reststddevs)
+
+    expander = VdotOperator(full(target, 1.)).adjoint
+    m = means[1]
+    L = np.log(target.k_lengths[-1]) - np.log(target.k_lengths[1])
+
+    from scipy.stats import norm
+    wienermean = np.sqrt(3/L)*np.abs(m)/norm.ppf(wienersigmaprob)
+    wienermean = np.log(wienermean)
+
+    sigma = Adder(full(expander.domain, wienermean)) @ (
+        wienersigmastddev*ducktape(expander.domain, None, keys[2]))
+    sigma = expander @ sigma.exp()
+    smooth = _TwoLogIntegrations(target).ducktape(keys[1])*sigma
+    return rest.ducktape(keys[0]) + smooth
+
+
+class Normalization(Operator):
+    def __init__(self, domain):
+        self._domain = self._target = makeDomain(domain)
+        hspace = self._domain[0].harmonic_partner
+        pd = PowerDistributor(hspace, power_space=self._domain[0])
+        # TODO Does not work on sphere yet
+        self._cst = pd.adjoint(full(pd.target, hspace.scalar_dvol))
+        self._specsum = SpecialSum(self._domain)
+
+    def apply(self, x):
+        self._check_input(x)
+        return self._specsum(self._cst*x).one_over()*x
+
+
+class SpecialSum(EndomorphicOperator):
+    def __init__(self, domain):
+        self._domain = makeDomain(domain)
+        self._capability = self.TIMES | self.ADJOINT_TIMES
+
+    def apply(self, x, mode):
+        self._check_input(x, mode)
+        return full(self._tgt(mode), x.sum())
+
+
+def WPAmplitude(target, means, stddevs, wienersigmastddev, wienersigmaprob,
+              keys):
+    op = LogIntegratedWienerProcess(target, means, stddevs, wienersigmastddev,
+                                    wienersigmaprob, keys)
+    return Normalization(target) @ op.exp()

nifty5/library/wiener_process_integrated_amplitude.py

-import nifty5 as ift
-import numpy as np
-
-
-class _TwoLogIntegrations(ift.LinearOperator):
-    def __init__(self, target):
-        self._target = ift.makeDomain(target)
-        self._domain = ift.makeDomain(
-            ift.UnstructuredDomain(self.target.shape[0] - 2))
-        self._capability = self.TIMES | self.ADJOINT_TIMES
-        if not isinstance(self._target[0], ift.PowerSpace):
-            raise TypeError
-        logk_lengths = np.log(self._target[0].k_lengths[1:])
-        self._logvol = logk_lengths[1:] - logk_lengths[:-1]
-
-    def apply(self, x, mode):
-        self._check_input(x, mode)
-        if mode == self.TIMES:
-            x = x.to_global_data()
-            res = np.empty(self._target.shape)
-            res[0] = 0
-            res[1] = 0
-            res[2:] = np.cumsum(x*self._logvol)
-            res[2:] = np.cumsum(res[2:]*self._logvol)
-            return ift.from_global_data(self._target, res)
-        else:
-            x = x.to_global_data()
-            res = np.empty(self._target.shape)
-            res[2:] = np.cumsum(x[2:][::-1])[::-1]*self._logvol
-            res[2:] = np.cumsum(res[2:][::-1])[::-1]*self._logvol
-            return ift.from_global_data(self._domain, res[2:])
-
-
-class _Rest(ift.LinearOperator):
-    def __init__(self, target):
-        self._target = ift.makeDomain(target)
-        self._domain = ift.makeDomain(ift.UnstructuredDomain(3))
-        self._logk_lengths = np.log(self._target[0].k_lengths[1:])
-        self._logk_lengths -= self._logk_lengths[0]
-        self._capability = self.TIMES | self.ADJOINT_TIMES
-
-    def apply(self, x, mode):
-        self._check_input(x, mode)
-        x = x.to_global_data()
-        res = np.empty(self._tgt(mode).shape)
-        if mode == self.TIMES:
-            res[0] = x[0]
-            res[1:] = x[1]*self._logk_lengths + x[2]
-        else:
-            res[0] = x[0]
-            res[1] = np.vdot(self._logk_lengths, x[1:])
-            res[2] = np.sum(x[1:])
-        return ift.from_global_data(self._tgt(mode), res)