Merge remote-tracking branch 'origin/NIFTy_5' into mpi_experiments

nifty5/__init__.py

@@ -47,8 +47,8 @@ from .operators.simple_linear_operators import (
     VdotOperator, SumReductionOperator, ConjugationOperator, Realizer,
     FieldAdapter, GeometryRemover, NullOperator)
 from .operators.energy_operators import (
-    EnergyOperator, GaussianEnergy, PoissonianEnergy, BernoulliEnergy,
-    Hamiltonian, SampledKullbachLeiblerDivergence)
+    EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood,
+    BernoulliEnergy, Hamiltonian, SampledKullbachLeiblerDivergence)
 
 from .probing import probe_with_posterior_samples, probe_diagonal, \
     StatCalculator
@@ -78,6 +78,7 @@ from .library.los_response import LOSResponse
 
 from .library.wiener_filter_curvature import WienerFilterCurvature
 from .library.correlated_fields import CorrelatedField, MfCorrelatedField
+from .library.adjust_variances import make_adjust_variances
 
 from . import extra
 

nifty5/library/adjust_variances.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 ..operators.energy_operators import Hamiltonian, InverseGammaLikelihood
+from ..operators.scaling_operator import ScalingOperator
+
+
+def make_adjust_variances(a, xi, position, samples=[], scaling=None, ic_samp=None):
+    """ Creates a Hamiltonian for constant likelihood optimizations.
+
+    Constructs a Hamiltonian to solve constant likelihood optimizations of the
+    form phi = a * xi under the constraint that phi remains constant.
+
+    Parameters
+    ----------
+    a : Operator
+        Operator which gives the amplitude when evaluated at a position
+    xi : Operator
+        Operator which gives the excitation when evaluated at a position
+    postion : Field, MultiField
+        Position of the whole problem
+    samples : Field, MultiField
+        Residual samples of the whole problem
+    scaling : Float
+        Optional rescaling of the Likelihood
+    ic_samp : Controller
+        Iteration Controller for Hamiltonian
+
+    Returns
+    -------
+    Hamiltonian
+        A Hamiltonian that can be used for further minimization
+    """
+
+    d = a*xi
+    d = (d.conjugate()*d).real
+    n = len(samples)
+    if n > 0:
+        d_eval = 0.
+        for i in range(n):
+            d_eval = d_eval + d(position + samples[i])
+        d_eval = d_eval/n
+    else:
+        d_eval = d(position)
+
+    x = (a.conjugate()*a).real
+    if scaling is not None:
+        x = ScalingOperator(scaling, x.target)(x)
+
+    return Hamiltonian(InverseGammaLikelihood(x, d_eval), ic_samp=ic_samp)

nifty5/linearization.py

@@ -93,6 +93,17 @@ class Linearization(object):
     def __rsub__(self, other):
         return (-self).__add__(other)
 
+    def __truediv__(self, other):
+        if isinstance(other, Linearization):
+            return self.__mul__(other.inverse())
+        return self.__mul__(1./other)
+
+    def __rtruediv__(self, other):
+        return self.inverse().__mul__(other)
+
+    def inverse(self):
+        return self.new(1./self._val, makeOp(-1./(self._val**2))(self._jac))
+
     def __mul__(self, other):
         from .sugar import makeOp
         if isinstance(other, Linearization):
@@ -170,11 +181,11 @@ class Linearization(object):
     @staticmethod
     def make_partial_var(field, constants, want_metric=False):
         from .operators.scaling_operator import ScalingOperator
-        from .operators.simple_linear_operators import NullOperator
+        from .operators.block_diagonal_operator import BlockDiagonalOperator
         if len(constants) == 0:
             return Linearization.make_var(field, want_metric)
         else:
             ops = [ScalingOperator(0. if key in constants else 1., dom)
                    for key, dom in field.domain.items()]
-            bdop = BlockDiagonalOperator(fielld.domain, tuple(ops))
+            bdop = BlockDiagonalOperator(field.domain, tuple(ops))
             return Linearization(field, bdop, want_metric=want_metric)

nifty5/multi_field.py

@@ -24,6 +24,7 @@ from . import utilities
 from .compat import *
 from .field import Field
 from .multi_domain import MultiDomain
+from .domain_tuple import DomainTuple
 
 
 class MultiField(object):
@@ -52,6 +53,10 @@ class MultiField(object):
     @staticmethod
     def from_dict(dict, domain=None):
         if domain is None:
+            for dd in dict.values():
+                if not isinstance(dd.domain, DomainTuple):
+                    raise TypeError('Values of dictionary need to be Fields '
+                                    'defined on DomainTuples.')
             domain = MultiDomain.make({key: v._domain
                                        for key, v in dict.items()})
         res = tuple(dict[key] if key in dict else Field(dom, 0)
@@ -61,6 +66,11 @@ class MultiField(object):
     def to_dict(self):
         return {key: val for key, val in zip(self._domain.keys(), self._val)}
 
+    def update(self, other):
+        foo = self.to_dict()
+        foo.update(other.to_dict())
+        return MultiField.from_dict(foo)
+
     def __getitem__(self, key):
         return self._val[self._domain.idx[key]]
 

nifty5/operators/energy_operators.py

@@ -113,6 +113,22 @@ class PoissonianEnergy(EnergyOperator):
         return res.add_metric(metric)
 
 
+class InverseGammaLikelihood(EnergyOperator):
+    def __init__(self, op, d):
+        self._op, self._d = op, d
+        self._domain = d.domain
+
+    def apply(self, x):
+        x = self._op(x)
+        res = 0.5*(x.log().sum() + (1./x).vdot(self._d))
+        if not isinstance(x, Linearization):
+            return Field.scalar(res)
+        if not x.want_metric:
+            return res
+        metric = SandwichOperator.make(x.jac, makeOp(0.5/(x.val**2)))
+        return res.add_metric(metric)
+
+
 class BernoulliEnergy(EnergyOperator):
     def __init__(self, p, d):
         self._p = p

test/test_energies/test_consistency.py

@@ -56,6 +56,20 @@ class Energy_Tests(unittest.TestCase):
 #         energy = ift.QuadraticEnergy(s[0], ift.makeOp(s[1]), s[2])
 #         ift.extra.check_value_gradient_consistency(energy)
 
+    @expand(
+        product([
+            ift.GLSpace(15),
+            ift.RGSpace(64, distances=.789),
+            ift.RGSpace([32, 32], distances=.789)
+        ], [4, 78, 23]))
+    def testInverseGammaLikelihood(self, space, seed):
+        model = self.make_model(space_key='s1', space=space, seed=seed)['s1']
+        d = np.random.normal(10, size=space.shape)**2
+        d = ift.Field.from_global_data(space, d)
+        energy = ift.InverseGammaLikelihood(ift.exp, d)
+        ift.extra.check_value_gradient_consistency(energy, model, tol=1e-7)
+
+
     @expand(product(
         [ift.GLSpace(15),
          ift.RGSpace(64, distances=.789),
@@ -65,7 +79,6 @@ class Energy_Tests(unittest.TestCase):
     def testPoissonian(self, space, seed):
         model = self.make_model(
             space_key='s1', space=space, seed=seed)['s1']
-        model = ift.exp(model)
         d = np.random.poisson(120, size=space.shape)
         d = ift.Field.from_global_data(space, d)
         energy = ift.PoissonianEnergy(ift.exp, d)

test/test_multi_field.py

@@ -56,3 +56,19 @@ class Test_Functionality(unittest.TestCase):
         f1 = op2(ift.full(dom, 1))
         for val in f1.values():
             assert_equal((val == 40).all(), True)
+
+    def test_update(self):
+        dom = ift.RGSpace(10)
+        f1 = ift.from_random('normal', domain=dom)
+        f2 = ift.from_random('normal', domain=dom)
+        f_new = ift.MultiField.from_dict({'dom1': f1, 'dom2': f2})
+        f3 = ift.from_random('normal', domain=dom)
+        f4 = ift.from_random('normal', domain=dom)
+        f5 = ift.from_random('normal', domain=dom)
+        f_old = ift.MultiField.from_dict({'dom1': f3, 'dom2': f4, 'dom3': f5})
+
+        updated = f_old.update(f_new).to_dict()
+        updated_true = f_old.to_dict()
+        updated_true.update(f_new.to_dict())
+        for key, val in updated.items():
+            assert_equal((val == updated_true[key]).all(), True)