Commit 66934b83 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'adjust_variances_but_right' into 'NIFTy_5'

Adjust variances but right

See merge request ift/nifty-dev!100
parents 9cc4c2a6 bd7e4a8c
......@@ -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
......
# 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)
......@@ -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):
......
......@@ -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
......
......@@ -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)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment