improve machinery for posterior samples

demos/log_normal_wiener_filter.py

@@ -64,8 +64,7 @@ if __name__ == "__main__":
     m = ht(me[0].position)
 
     # Plotting
-    plotdict = {"xlabel": "Pixel index", "ylabel": "Pixel index",
-                "colormap": "Planck-like"}
+    plotdict = {"colormap": "Planck-like"}
     ift.plot(ht(mock_signal), name="mock_signal.png", **plotdict)
     logdata = np.log(ift.dobj.to_global_data(data.val)).reshape(signal_space.shape)
     ift.plot(ift.Field(signal_space, val=ift.dobj.from_global_data(logdata)),

demos/paper_demos/cartesian_wiener_filter.py

@@ -83,8 +83,8 @@ if __name__ == "__main__":
     noiseless_data = R(mock_signal)
     noise_amplitude = noiseless_data.val.std()/signal_to_noise
     # Setting up the noise covariance and drawing a random noise realization
-    ndiag = ift.Field.full(data_domain, noise_amplitude**2)
-    N = ift.DiagonalOperator(ndiag)
+    #ndiag = ift.Field.full(data_domain, noise_amplitude**2)
+    N = ift.ScalingOperator(noise_amplitude**2, data_domain)
     noise = ift.Field.from_random(
         domain=data_domain, random_type='normal',
         std=noise_amplitude, mean=0)
@@ -108,12 +108,5 @@ if __name__ == "__main__":
     ift.plot(ift.Field(plot_space,val=data.val), name='data.png', **plotdict)
     ift.plot(ift.Field(plot_space,val=m.val), name='map.png', **plotdict)
     # sampling the uncertainty map
-    sample_variance = ift.Field.zeros(signal_domain)
-    sample_mean = ift.Field.zeros(signal_domain)
-    n_samples = 10
-    for i in range(n_samples):
-        sample = ht(wiener_curvature.generate_posterior_sample()) + m
-        sample_variance += sample**2
-        sample_mean += sample
-    variance = sample_variance/n_samples - (sample_mean/n_samples)**2
+    mean, variance = ift.probe_with_posterior_samples(wiener_curvature, ht, 10)
     ift.plot(ift.Field(plot_space, val=ift.sqrt(variance).val), name="uncertainty.png", **plotdict)

demos/wiener_filter_via_hamiltonian.py

@@ -68,12 +68,6 @@ if __name__ == "__main__":
     sample_variance = ift.Field.zeros(s_space)
     sample_mean = ift.Field.zeros(s_space)
 
-    n_samples = 50
-    for i in range(n_samples):
-        sample = ht(D.generate_posterior_sample() + m)
-        sample_variance += sample**2
-        sample_mean += sample
-    sample_mean /= n_samples
-    sample_variance /= n_samples
-    variance = sample_variance - sample_mean**2
+    mean, variance = ift.probe_with_posterior_samples(D, ht, 50)
     ift.plot(variance, name="variance.png", colormap="Planck-like")
+    ift.plot(mean, name="mean.png", colormap="Planck-like")

nifty4/__init__.py

@@ -32,6 +32,7 @@ from .field import Field, sqrt, exp, log
 from .probing.prober import Prober
 from .probing.diagonal_prober_mixin import DiagonalProberMixin
 from .probing.trace_prober_mixin import TraceProberMixin
+from .probing.utils import probe_with_posterior_samples
 
 from .minimization.line_search import LineSearch
 from .minimization.line_search_strong_wolfe import LineSearchStrongWolfe

nifty4/field.py

@@ -217,7 +217,7 @@ class Field(object):
     def real(self):
         """ The real part of the field (data is not copied)."""
         if not np.issubdtype(self.dtype, np.complexfloating):
-            raise ValueError(".real called on a non-complex Field")
+            return self
         return Field(self._domain, self.val.real)
 
     @property

nifty4/library/wiener_filter_curvature.py

@@ -74,32 +74,8 @@ class WienerFilterCurvature(EndomorphicOperator):
             covariance.
         """
 
-        power = power_analyze(sqrt(self.S.diagonal()))
-        mock_signal = power_synthesize(power, real_signal=True)
-
-        noise = self.N.diagonal()
-
-        mock_noise = Field.from_random(random_type="normal",
-                                       domain=self.N.domain, dtype=noise.dtype)
-        mock_noise *= sqrt(noise)
-
-        mock_data = self.R(mock_signal) + mock_noise
-
-        mock_j = self.R.adjoint_times(self.N.inverse_times(mock_data))
-        mock_m = self.inverse_times(mock_j)
-        return mock_signal - mock_m
-
-    def generate_posterior_sample2(self):
-        power = self.S.diagonal()
-        mock_signal = Field.from_random(random_type="normal",
-                                        domain=self.S.domain,
-                                        dtype=power.dtype)
-        mock_signal *= sqrt(power)
-
-        noise = self.N.diagonal()
-        mock_noise = Field.from_random(random_type="normal",
-                                       domain=self.N.domain, dtype=noise.dtype)
-        mock_noise *= sqrt(noise)
+        mock_signal = self.S.generate_posterior_sample()
+        mock_noise = self.N.generate_posterior_sample()
 
         mock_data = self.R(mock_signal) + mock_noise
 

nifty4/operators/diagonal_operator.py

@@ -18,7 +18,7 @@
 
 from __future__ import division
 import numpy as np
-from ..field import Field
+from ..field import Field, sqrt
 from ..domain_tuple import DomainTuple
 from .endomorphic_operator import EndomorphicOperator
 from .. import utilities
@@ -140,3 +140,27 @@ class DiagonalOperator(EndomorphicOperator):
     def adjoint(self):
         return DiagonalOperator(self._diagonal.conjugate(), self._domain,
                                 self._spaces)
+
+    def generate_posterior_sample(self):
+        """ Generates a posterior sample from a Gaussian distribution with
+        given mean and covariance.
+
+        This method generates samples by setting up the observation and
+        reconstruction of a mock signal in order to obtain residuals of the
+        right correlation which are added to the given mean.
+
+        Returns
+        -------
+        sample : Field
+            Returns the a sample from the Gaussian of given mean and
+            covariance.
+        """
+
+        if self._spaces is not None:
+            raise ValueError("Cannot draw (yet) from this operator")
+
+        res = Field.from_random(random_type="normal",
+                                domain=self._domain,
+                                dtype=self._diagonal.dtype)
+        res *= sqrt(self._diagonal)
+        return res

nifty4/operators/scaling_operator.py

@@ -88,3 +88,23 @@ class ScalingOperator(EndomorphicOperator):
             return self.TIMES | self.ADJOINT_TIMES
         return (self.TIMES | self.ADJOINT_TIMES |
                 self.INVERSE_TIMES | self.ADJOINT_INVERSE_TIMES)
+
+    def generate_posterior_sample(self):
+        """ Generates a posterior sample from a Gaussian distribution with
+        given mean and covariance.
+
+        This method generates samples by setting up the observation and
+        reconstruction of a mock signal in order to obtain residuals of the
+        right correlation which are added to the given mean.
+
+        Returns
+        -------
+        sample : Field
+            Returns the a sample from the Gaussian of given mean and
+            covariance.
+        """
+
+        return Field.from_random(random_type="normal",
+                                 domain=self._domain,
+                                 std=np.sqrt(self._factor),
+                                 dtype=np.result_type(self._factor))

nifty4/probing/utils.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-2017 Max-Planck-Society
+#
+# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
+# and financially supported by the Studienstiftung des deutschen Volkes.
+
+from builtins import object
+
+
+class StatCalculator(object):
+    def __init__(self):
+        self._count = 0
+
+    def add(self, value):
+        self._count += 1
+        if self._count == 1:
+            self._mean = 1.*value
+            self._M2 = 0.*value
+        else:
+            delta = value - self._mean
+            self._mean += delta*(1./self._count)
+            delta2 = value - self._mean
+            self._M2 += delta*delta2
+
+    @property
+    def mean(self):
+        if self._count == 0:
+            raise RuntimeError
+        return 1.*self._mean
+
+    @property
+    def var(self):
+        if self._count < 2:
+            raise RuntimeError
+        return self._M2 * (1./(self._count-1))
+
+
+def probe_with_posterior_samples(op, post_op, nprobes):
+    sc = StatCalculator()
+    for i in range(nprobes):
+        sample = post_op(op.generate_posterior_sample())
+        sc.add(sample)
+    return sc.mean, sc.var