diff --git a/nifty5/__init__.py b/nifty5/__init__.py
index 0d93e1632737f5947b607a0d73321fa06bddf487..51804da16f061a5d90d9465bcb966752df8f2e76 100644
--- a/nifty5/__init__.py
+++ b/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
diff --git a/nifty5/library/adjust_variances.py b/nifty5/library/adjust_variances.py
new file mode 100644
index 0000000000000000000000000000000000000000..c1c28ad8661ce5f9a268906234ac2cdbabc188c1
--- /dev/null
+++ b/nifty5/library/adjust_variances.py
@@ -0,0 +1,68 @@
+# 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)
diff --git a/nifty5/linearization.py b/nifty5/linearization.py
index 30a9b57bad11b5e264a0997bde6c0e0e2bbefe10..63f154e321e6c046f572a4125aad3d42f304d0ae 100644
--- a/nifty5/linearization.py
+++ b/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)
diff --git a/nifty5/multi_field.py b/nifty5/multi_field.py
index 10795f819a834a36687be3eec466f2d926c4dd43..324ee688579807437c749dfb925a0a503653303d 100644
--- a/nifty5/multi_field.py
+++ b/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]]
diff --git a/nifty5/operators/energy_operators.py b/nifty5/operators/energy_operators.py
index a68ea92527a4850bd79976cb25eab87d2e8d9d7c..dbde8b1db46cbd7da0f3ffff6a8427c5ebcc1c11 100644
--- a/nifty5/operators/energy_operators.py
+++ b/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
diff --git a/test/test_energies/test_consistency.py b/test/test_energies/test_consistency.py
index 0bf29ddc76949496e21026bf3a86ab15acee1eec..91a090c75e1941fb605d9cae0ca3260e6414f7ff 100644
--- a/test/test_energies/test_consistency.py
+++ b/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)
diff --git a/test/test_multi_field.py b/test/test_multi_field.py
index 580d9ae5591842fb0c8c5af1dba133c86dbad2a6..5a964d10cfacba292e6bcadab9754c1f8d0ed419 100644
--- a/test/test_multi_field.py
+++ b/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)