From 62e6b77465718505fdaa3233026e69d50e29dcc4 Mon Sep 17 00:00:00 2001
From: Philipp Arras <parras@mpa-garching.mpg.de>
Date: Wed, 19 Sep 2018 17:00:45 +0200
Subject: [PATCH] Add do_adjust_variances
---
nifty5/__init__.py | 2 +-
nifty5/library/adjust_variances.py | 44 +++++++++++++++++++++++++++++-
2 files changed, 44 insertions(+), 2 deletions(-)
diff --git a/nifty5/__init__.py b/nifty5/__init__.py
index 4d15da59f..06a553232 100644
--- a/nifty5/__init__.py
+++ b/nifty5/__init__.py
@@ -78,7 +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 .library.adjust_variances import make_adjust_variances, do_adjust_variances
from . import extra
diff --git a/nifty5/library/adjust_variances.py b/nifty5/library/adjust_variances.py
index a44da2603..2ac585907 100644
--- a/nifty5/library/adjust_variances.py
+++ b/nifty5/library/adjust_variances.py
@@ -19,11 +19,20 @@
from __future__ import absolute_import, division, print_function
from ..compat import *
+from ..minimization.energy_adapter import EnergyAdapter
+from ..multi_domain import MultiDomain
+from ..multi_field import MultiField
+from ..operators.distributors import PowerDistributor
from ..operators.energy_operators import Hamiltonian, InverseGammaLikelihood
from ..operators.scaling_operator import ScalingOperator
+from ..operators.simple_linear_operators import FieldAdapter
-def make_adjust_variances(a, xi, position, samples=[], scaling=None,
+def make_adjust_variances(a,
+ xi,
+ position,
+ samples=[],
+ scaling=None,
ic_samp=None):
""" Creates a Hamiltonian for constant likelihood optimizations.
@@ -67,3 +76,36 @@ def make_adjust_variances(a, xi, position, samples=[], scaling=None,
x = ScalingOperator(scaling, x.target)(x)
return Hamiltonian(InverseGammaLikelihood(d_eval)(x), ic_samp=ic_samp)
+
+
+def do_adjust_variances(position, xi, amplitude_model, minimizer, samples=[]):
+ h_space = xi.target[0]
+ pd = PowerDistributor(h_space, amplitude_model.target[0])
+ a = pd(amplitude_model)
+ xi = FieldAdapter(MultiDomain.make({"xi": h_space}), "xi")
+
+ axi_domain = MultiDomain.union([a.domain, xi.domain])
+ ham = make_adjust_variances(
+ a, xi, position.extract(axi_domain), samples=samples)
+ a_pos = position.extract(a.domain)
+
+ # Minimize
+ # FIXME Try also KL here
+ e = EnergyAdapter(a_pos, ham, constants=[], want_metric=True)
+ e, _ = minimizer(e)
+
+ # Update position
+ s_h_old = (a*xi)(position.extract(axi_domain))
+
+ position = position.to_dict()
+ position['xi'] = s_h_old/a(e.position)
+ position = MultiField.from_dict(position)
+ position = MultiField.union([position, e.position])
+
+ s_h_new = (a*xi)(position.extract(axi_domain))
+
+ import numpy as np
+ # TODO Move this into the tests
+ np.testing.assert_allclose(s_h_new.to_global_data(),
+ s_h_old.to_global_data())
+ return position
--
GitLab