remove old critical filter; replace demo by new critical filter with linear nonlinearity

demos/critical_filtering.py

-import numpy as np
 import nifty4 as ift
-
+from nifty4.library.nonlinearities import Linear
+import numpy as np
 np.random.seed(42)
-# np.seterr(all="raise",under="ignore")
 
 
+def adjust_zero_mode(m0, t0):
+    mtmp = m0.to_global_data().copy()
+    zero_position = len(m0.shape)*(0,)
+    zero_mode = mtmp[zero_position]
+    mtmp[zero_position] = zero_mode / abs(zero_mode)
+    ttmp = t0.to_global_data().copy()
+    ttmp[0] += 2 * np.log(abs(zero_mode))
+    return (ift.Field.from_global_data(m0.domain, mtmp),
+            ift.Field.from_global_data(t0.domain, ttmp))
+
 if __name__ == "__main__":
+
+    noise_level = 1.
+    p_spec = (lambda k: (.5 / (k + 1) ** 3))
+
+    nonlinearity = Linear()
     # Set up position space
-    s_space = ift.RGSpace([128, 128])
-    # s_space = ift.HPSpace(32)
+    s_space = ift.RGSpace((128, 128))
+    h_space = s_space.get_default_codomain()
 
     # Define harmonic transformation and associated harmonic space
-    h_space = s_space.get_default_codomain()
-    HT = ift.HarmonicTransformOperator(h_space, s_space)
+    HT = ift.HarmonicTransformOperator(h_space, target=s_space)
 
-    # Set up power space
+    # Setting up power space
     p_space = ift.PowerSpace(h_space,
                              binbounds=ift.PowerSpace.useful_binbounds(
                                  h_space, logarithmic=True))
+    s_spec = ift.Field.full(p_space, 1.)
+    # Choosing the prior correlation structure and defining
+    # correlation operator
+    p = ift.PS_field(p_space, p_spec)
+    log_p = ift.log(p)
+    S = ift.create_power_operator(h_space, power_spectrum=s_spec)
+
+    # Drawing a sample sh from the prior distribution in harmonic space
+    sh = ift.power_synthesize(s_spec)
+
+    # Choosing the measurement instrument
+    # Instrument = SmoothingOperator(s_space, sigma=0.01)
+    mask = np.ones(s_space.shape)
+    mask[6000:8000] = 0.
+    mask = ift.Field.from_global_data(s_space, mask)
+
+    MaskOperator = ift.DiagonalOperator(mask)
+    R = ift.GeometryRemover(s_space)
+    R = R*MaskOperator
+    # R = R*HT
+    # R = R * ift.create_harmonic_smoothing_operator((harmonic_space,), 0,
+    #                                                response_sigma)
+    MeasurementOperator = R
+
+    d_space = MeasurementOperator.target
+
+    Distributor = ift.PowerDistributor(target=h_space, power_space=p_space)
+    power = Distributor(ift.exp(0.5*log_p))
+    # Creating the mock data
+    true_sky = nonlinearity(HT(power*sh))
+    noiseless_data = MeasurementOperator(true_sky)
+    noise_amplitude = noiseless_data.val.std()*noise_level
+    N = ift.ScalingOperator(noise_amplitude**2, d_space)
+    n = ift.Field.from_random(
+        domain=d_space, random_type='normal',
+        std=noise_amplitude, mean=0)
+    # Creating the mock data
+    d = noiseless_data + n
 
-    # Choose the prior correlation structure and defining correlation operator
-    p_spec = (lambda k: (.5 / (k + 1) ** 3))
-    S = ift.create_power_operator(h_space, power_spectrum=p_spec)
+    m0 = ift.power_synthesize(ift.Field.full(p_space, 1e-7))
+    t0 = ift.Field.full(p_space, -4.)
+    power0 = Distributor.times(ift.exp(0.5 * t0))
 
-    # Draw a sample sh from the prior distribution in harmonic space
-    sp = ift.PS_field(p_space, p_spec)
-    sh = ift.power_synthesize(sp, real_signal=True)
+    IC1 = ift.GradientNormController(name="IC1", iteration_limit=100,
+                                     tol_abs_gradnorm=1e-3)
+    LS = ift.LineSearchStrongWolfe(c2=0.02)
+    minimizer = ift.RelaxedNewton(IC1, line_searcher=LS)
 
-    # Choose the measurement instrument
-    Instrument = ift.ScalingOperator(1., s_space)
+    ICI = ift.GradientNormController(iteration_limit=500,
+                                     tol_abs_gradnorm=1e-3)
+    inverter = ift.ConjugateGradient(controller=ICI)
 
-    # Add a harmonic transformation to the instrument
-    R = Instrument*HT
+    for i in range(20):
+        power0 = Distributor(ift.exp(0.5*t0))
+        map0_energy = ift.library.NonlinearWienerFilterEnergy(
+            m0, d, MeasurementOperator, nonlinearity, HT, power0, N, S,
+            inverter=inverter)
 
-    noiseless_data = R(sh)
-    signal_to_noise = 1.
-    noise_amplitude = noiseless_data.val.std()/signal_to_noise
-    N = ift.ScalingOperator(noise_amplitude**2, s_space)
-    n = ift.Field.from_random("normal", s_space, std=noise_amplitude, mean=0)
+        # Minimization with chosen minimizer
+        map0_energy, convergence = minimizer(map0_energy)
+        m0 = map0_energy.position
 
-    # Create mock data
-    d = noiseless_data + n
+        # Updating parameters for correlation structure reconstruction
+        D0 = map0_energy.curvature
 
-    # The information source
-    j = R.adjoint_times(N.inverse_times(d))
-    realized_power = ift.log(ift.power_analyze(sh,
-                                               binbounds=p_space.binbounds))
-    data_power = ift.log(ift.power_analyze(HT.adjoint_times(d),
-                                           binbounds=p_space.binbounds))
-    plotdict = {"colormap": "Planck-like"}
-    ift.plot(d, name="data.png", **plotdict)
-    ift.plot(HT(sh), name="mock_signal.png", **plotdict)
+        # Initializing the power energy with updated parameters
+        power0_energy = ift.library.NonlinearPowerEnergy(
+            position=t0, d=d, N=N, xi=m0, D=D0, ht=HT,
+            Instrument=MeasurementOperator, nonlinearity=nonlinearity,
+            Distribution=Distributor, sigma=1., samples=2, inverter=inverter)
 
-    IC1 = ift.GradientNormController(name="IC1", iteration_limit=100,
-                                     tol_abs_gradnorm=1e-15)
-    minimizer = ift.RelaxedNewton(IC1)
+        power0_energy = minimizer(power0_energy)[0]
 
-    ICI = ift.GradientNormController(iteration_limit=500,
-                                     tol_abs_gradnorm=1e-10)
-    map_inverter = ift.ConjugateGradient(controller=ICI)
-
-    ICI2 = ift.GradientNormController(iteration_limit=200,
-                                      tol_abs_gradnorm=1e-10)
-    power_inverter = ift.ConjugateGradient(controller=ICI2)
-
-    # Set starting position
-    flat_power = ift.Field.full(p_space, 1e-8)
-    m0 = ift.power_synthesize(flat_power, real_signal=True)
-    t0 = ift.Field.full(p_space, -7.)
-
-    for i in range(500):
-        S0 = ift.create_power_operator(h_space, power_spectrum=ift.exp(t0))
-
-        # Initialize non-linear Wiener Filter energy
-        map_energy = ift.library.WienerFilterEnergy(
-            position=m0, d=d, R=R, N=N, S=S0, inverter=map_inverter)
-        # Solve the Wiener Filter analytically
-        D0 = map_energy.curvature
-        m0 = D0.inverse_times(j)
-        # Initialize power energy with updated parameters
-        power_energy = ift.library.CriticalPowerEnergy(
-            position=t0, m=m0, D=D0, smoothness_prior=10., samples=3,
-            inverter=power_inverter)
-
-        power_energy = minimizer(power_energy)[0]
-
-        # Set new power spectrum
-        t0 = power_energy.position
-
-        # Plot current estimate
-        ift.dobj.mprint(i)
-        if i % 50 == 0:
-            ift.plot(HT(m0), name='map.png', **plotdict)
+        # Setting new power spectrum
+        t0 = power0_energy.position
+
+        # break degeneracy between amplitude and excitation by setting
+        # excitation monopole to 1
+        m0, t0 = adjust_zero_mode(m0, t0)
+
+    ift.plot(true_sky)
+    ift.plot(nonlinearity(HT(power0*m0)),
+             title='reconstructed_sky')
+    ift.plot(MeasurementOperator.adjoint_times(d))

nifty4/library/__init__.py

-from .critical_power_energy import CriticalPowerEnergy
-from .critical_power_curvature import CriticalPowerCurvature
 from .wiener_filter_energy import WienerFilterEnergy
 from .wiener_filter_curvature import WienerFilterCurvature
 from .noise_energy import NoiseEnergy

nifty4/library/critical_power_curvature.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 ..operators.inversion_enabler import InversionEnabler
-from ..operators.diagonal_operator import DiagonalOperator
-
-
-def CriticalPowerCurvature(theta, T, inverter):
-    theta = DiagonalOperator(theta)
-    return InversionEnabler(T+theta, inverter, theta.inverse_times)

nifty4/library/critical_power_energy.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 ..minimization.energy import Energy
-from ..operators.smoothness_operator import SmoothnessOperator
-from ..operators.power_distributor import PowerDistributor
-from .critical_power_curvature import CriticalPowerCurvature
-from ..utilities import memo
-from .. import Field, exp
-
-
-class CriticalPowerEnergy(Energy):
-    """The Energy of the power spectrum according to the critical filter.
-
-    It describes the energy of the logarithmic amplitudes of the power spectrum
-    of a Gaussian random field under reconstruction uncertainty with smoothness
-    and inverse gamma prior. It is used to infer the correlation structure of a
-    correlated signal.
-
-    Parameters
-    ----------
-    position : Field
-        The current position of this energy. (Logarithm of power spectrum)
-    m : Field
-        The map whose power spectrum is inferred. Needs to live in harmonic
-        signal space.
-    D : EndomorphicOperator, optional
-        The curvature of the Gaussian encoding the posterior covariance.
-        If not specified, the map is assumed to be no reconstruction.
-        default : None
-    alpha : float, optional
-        The spectral prior of the inverse gamma distribution. 1.0 corresponds
-        to non-informative.
-        default : 1.0
-    q : float, optional
-        The cutoff parameter of the inverse gamma distribution. 0.0 corresponds
-        to non-informative.
-        default : 0.0
-    smoothness_prior : float, optional
-        Controls the strength of the smoothness prior
-        default : 0.0
-    logarithmic : bool, optional
-        Whether smoothness acts on linear or logarithmic scale.
-        default : True
-    samples : int, optional
-        Number of samples used for the estimation of the uncertainty
-        corrections.
-        default : 3
-    w : Field, optional
-        The contribution from the map with or without uncertainty. It is used
-        to pass on the result of the costly sampling during the minimization.
-        default : None
-    inverter : ConjugateGradient
-        The inversion strategy to invert the curvature and to generate samples.
-        default : None
-    """
-
-    def __init__(self, position, m, D=None, alpha=1.0, q=0.,
-                 smoothness_prior=0., logarithmic=True, samples=3, w=None,
-                 inverter=None):
-        super(CriticalPowerEnergy, self).__init__(position=position)
-        self.m = m
-        self.D = D
-        self.samples = samples
-        self.alpha = float(alpha)
-        self.q = float(q)
-        self._smoothness_prior = smoothness_prior
-        self._logarithmic = logarithmic
-        self.T = SmoothnessOperator(domain=self.position.domain[0],
-                                    strength=smoothness_prior,
-                                    logarithmic=logarithmic)
-        self._inverter = inverter
-
-        if w is None:
-            Dist = PowerDistributor(target=self.m.domain,
-                                    power_space=self.position.domain[0])
-            if self.D is not None:
-                # MR FIXME: we should use stuff from probing.utils for this
-                w = Field.zeros(self.position.domain, dtype=self.m.dtype)
-                for i in range(self.samples):
-                    sample = self.D.draw_sample() + self.m
-                    w += Dist.adjoint_times(abs(sample)**2)
-
-                w *= 1./self.samples
-            else:
-                w = Dist.adjoint_times(abs(self.m)**2)
-        self._w = w
-
-        self._theta = exp(-self.position) * (self.q + self._w*0.5)
-        Tt = self.T(self.position)
-
-        energy = self._theta.sum()
-        energy += self.position.sum() * (self.alpha-0.5)
-        energy += 0.5*self.position.vdot(Tt)
-        self._value = energy.real
-
-        gradient = -self._theta
-        gradient += self.alpha-0.5
-        gradient += Tt
-        self._gradient = gradient.lock()
-
-    def at(self, position):
-        return self.__class__(position, self.m, D=self.D, alpha=self.alpha,
-                              q=self.q,
-                              smoothness_prior=self._smoothness_prior,
-                              logarithmic=self._logarithmic,
-                              samples=self.samples, w=self._w,
-                              inverter=self._inverter)
-
-    @property
-    def value(self):
-        return self._value
-
-    @property
-    def gradient(self):
-        return self._gradient
-
-    @property
-    @memo
-    def curvature(self):
-        return CriticalPowerCurvature(theta=self._theta, T=self.T,
-                                      inverter=self._inverter)