Commit 34a385d9 authored by Martin Reinecke's avatar Martin Reinecke

revamp, part 1

parent 2bbd8133
......@@ -19,39 +19,6 @@
import nifty5 as ift
import numpy as np
class GaussianEnergy2(ift.Operator):
def __init__(self, mean=None, covariance=None):
super(GaussianEnergy2, self).__init__()
self._mean = mean
self._icov = None if covariance is None else covariance.inverse
def __call__(self, x):
residual = x if self._mean is None else x-self._mean
icovres = residual if self._icov is None else self._icov(residual)
res = .5 * (residual*icovres).sum()
metric = ift.SandwichOperator.make(x.jac, self._icov)
return res.add_metric(metric)
class PoissonianEnergy2(ift.Operator):
def __init__(self, op, d):
super(PoissonianEnergy2, self).__init__()
self._op = op
self._d = d
def __call__(self, x):
x = self._op(x)
res = (x - self._d*x.log()).sum()
metric = ift.SandwichOperator.make(x.jac, ift.makeOp(1./x.val))
return res.add_metric(metric)
class MyHamiltonian(ift.Operator):
def __init__(self, lh):
super(MyHamiltonian, self).__init__()
self._lh = lh
self._prior = GaussianEnergy2()
def __call__(self, x):
return self._lh(x) + self._prior(x)
class EnergyAdapter(ift.Energy):
def __init__(self, position, op):
......@@ -145,14 +112,14 @@ if __name__ == '__main__':
data = ift.Field.from_global_data(d_space, data)
# Compute likelihood and Hamiltonian
likelihood = PoissonianEnergy2(lamb, data)
likelihood = ift.PoissonianEnergy(lamb, data)
ic_cg = ift.GradientNormController(iteration_limit=50)
ic_newton = ift.GradientNormController(name='Newton', iteration_limit=50,
tol_abs_gradnorm=1e-3)
minimizer = ift.RelaxedNewton(ic_newton)
# Minimize the Hamiltonian
H = MyHamiltonian(likelihood)
H = ift.Hamiltonian(likelihood)
H = EnergyAdapter(position, H)
#ift.extra.check_value_gradient_consistency(H)
H = H.make_invertible(ic_cg)
......
......@@ -25,37 +25,6 @@ def get_random_LOS(n_los):
ends = list(np.random.uniform(0, 1, (n_los, 2)).T)
return starts, ends
class GaussianEnergy2(ift.Operator):
def __init__(self, mean=None, covariance=None):
super(GaussianEnergy2, self).__init__()
self._mean = mean
self._icov = None if covariance is None else covariance.inverse
def __call__(self, x):
residual = x if self._mean is None else x-self._mean
icovres = residual if self._icov is None else self._icov(residual)
res = .5 * (residual*icovres).sum()
metric = ift.SandwichOperator.make(x.jac, self._icov)
return res.add_metric(metric)
class MyHamiltonian(ift.Operator):
def __init__(self, lh, ic_samp=None):
super(MyHamiltonian, self).__init__()
self._lh = lh
self._prior = GaussianEnergy2()
self._ic_samp = ic_samp
def __call__(self, x):
res = self._lh(x) + self._prior(x)
if self._ic_samp is None:
return self._lh(x) + self._prior(x)
else:
lhx = self._lh(x)
prx = self._prior(x)
mtr = ift.SamplingEnabler(lhx.metric, prx.metric.inverse,
self._ic_samp, prx.metric.inverse)
return (lhx+prx).add_metric(mtr)
class EnergyAdapter(ift.Energy):
def __init__(self, position, op):
super(EnergyAdapter, self).__init__(position)
......@@ -78,15 +47,6 @@ class EnergyAdapter(ift.Energy):
def metric(self):
return self._res.metric
class FieldPicker(ift.Operator):
def __init__(self, name_dom):
self._name_dom = name_dom
def __call__(self, x):
if isinstance(x, (ift.Field, ift.MultiField)):
return x[self._name_dom]
dom = x[self._name_dom].domain
return ift.Linearization(x._val[self._name_dom], ift.FieldAdapter(dom, self._name_dom, None))
if __name__ == '__main__':
# FIXME description of the tutorial
np.random.seed(42)
......@@ -108,7 +68,7 @@ if __name__ == '__main__':
# correlated_field,_ = ift.make_correlated_field(position_space, A)
# apply some nonlinearity
signal = lambda inp: correlated_field(inp).positive_tanh()
signal = lambda inp: correlated_field(inp).positive_tanh()
# Building the Line of Sight response
LOS_starts, LOS_ends = get_random_LOS(100)
......@@ -127,7 +87,7 @@ if __name__ == '__main__':
data = signal_response(MOCK_POSITION) + N.draw_sample()
# set up model likelihood
likelihood = lambda inp: GaussianEnergy2(mean=data, covariance=N)(signal_response(inp))
likelihood = lambda inp: ift.GaussianEnergy(mean=data, covariance=N)(signal_response(inp))
# set up minimization and inversion schemes
ic_cg = ift.GradientNormController(iteration_limit=10)
......@@ -136,7 +96,7 @@ if __name__ == '__main__':
minimizer = ift.RelaxedNewton(ic_newton)
# build model Hamiltonian
H = MyHamiltonian(likelihood, ic_sampling)
H = ift.Hamiltonian(likelihood, ic_sampling)
H = EnergyAdapter(MOCK_POSITION, H)
INITIAL_POSITION = ift.from_random('normal', H.position.domain)
......@@ -160,9 +120,7 @@ if __name__ == '__main__':
position = KL.position
ift.plot(signal(position), title="reconstruction")
ift.plot([A(position), A(MOCK_POSITION)],
title="power")
ift.plot([A(position), A(MOCK_POSITION)], title="power")
ift.plot_finish(nx=2, xsize=12, ysize=6, title="loop", name="loop.png")
sc = ift.StatCalculator()
......@@ -172,7 +130,6 @@ if __name__ == '__main__':
ift.plot(ift.sqrt(sc.var), title="std deviation")
powers = [A(s+position) for s in samples]
ift.plot([A(position), A(MOCK_POSITION)]+powers,
title="power")
ift.plot([A(position), A(MOCK_POSITION)]+powers, title="power")
ift.plot_finish(nx=3, xsize=16, ysize=5, title="results",
name="results.png")
......@@ -19,48 +19,25 @@
from __future__ import absolute_import, division, print_function
from ..compat import *
from ..operator import Operator
from ..library.gaussian_energy import GaussianEnergy
from ..minimization.energy import Energy
from ..models.variable import Variable
from ..operators.sampling_enabler import SamplingEnabler
from ..utilities import memo
class Hamiltonian(Energy):
def __init__(self, lh, iteration_controller_sampling=None):
"""
lh: Likelihood (energy object)
prior:
"""
super(Hamiltonian, self).__init__(lh._position)
class Hamiltonian(Operator):
def __init__(self, lh, ic_samp=None):
super(Hamiltonian, self).__init__()
self._lh = lh
self._ic_samp = iteration_controller_sampling
self._prior = GaussianEnergy(Variable(self._position))
self._prior = GaussianEnergy()
self._ic_samp = ic_samp
def at(self, position):
return self.__class__(self._lh.at(position), self._ic_samp)
@property
@memo
def value(self):
return self._lh.value + self._prior.value
@property
@memo
def gradient(self):
return self._lh.gradient + self._prior.gradient
@property
@memo
def metric(self):
prior_mtr = self._prior.metric
def __call__(self, x):
res = self._lh(x) + self._prior(x)
if self._ic_samp is None:
return self._lh.metric + prior_mtr
return self._lh(x) + self._prior(x)
else:
return SamplingEnabler(self._lh.metric, prior_mtr.inverse,
self._ic_samp, prior_mtr.inverse)
def __str__(self):
res = 'Likelihood:\t{:.2E}\n'.format(self._lh.value)
res += 'Prior:\t\t{:.2E}'.format(self._prior.value)
return res
lhx = self._lh(x)
prx = self._prior(x)
mtr = SamplingEnabler(lhx.metric, prx.metric.inverse,
self._ic_samp, prx.metric.inverse)
return (lhx+prx).add_metric(mtr)
......@@ -19,53 +19,19 @@
from __future__ import absolute_import, division, print_function
from ..compat import *
from ..minimization.energy import Energy
from ..operator import Operator
from ..operators.sandwich_operator import SandwichOperator
from ..utilities import memo
# MR FIXME documentation incomplete
class GaussianEnergy(Energy):
def __init__(self, inp, mean=None, covariance=None):
"""
inp: Model object
value = 0.5 * s.vdot(s), i.e. a log-Gauss distribution with unit
covariance
"""
super(GaussianEnergy, self).__init__(inp._position)
self._inp = inp
class GaussianEnergy(Operator):
def __init__(self, mean=None, covariance=None):
super(GaussianEnergy, self).__init__()
self._mean = mean
self._cov = covariance
def at(self, position):
return self.__class__(self._inp.at(position), self._mean, self._cov)
@property
@memo
def _residual(self):
return self._inp.value if self._mean is None else \
self._inp.value - self._mean
@property
@memo
def _icovres(self):
return self._residual if self._cov is None else \
self._cov.inverse_times(self._residual)
@property
@memo
def value(self):
return .5 * self._residual.vdot(self._icovres).real
@property
@memo
def gradient(self):
return self._inp.jacobian.adjoint_times(self._icovres).real
@property
@memo
def metric(self):
return SandwichOperator.make(
self._inp.jacobian,
None if self._cov is None else self._cov.inverse)
self._icov = None if covariance is None else covariance.inverse
def __call__(self, x):
residual = x if self._mean is None else x-self._mean
icovres = residual if self._icov is None else self._icov(residual)
res = .5*(residual*icovres).sum()
metric = SandwichOperator.make(x.jac, self._icov)
return res.add_metric(metric)
......@@ -21,44 +21,19 @@ from __future__ import absolute_import, division, print_function
from numpy import inf, isnan
from ..compat import *
from ..minimization.energy import Energy
from ..operator import Operator
from ..operators.sandwich_operator import SandwichOperator
from ..sugar import log, makeOp
from ..sugar import makeOp
class PoissonianEnergy(Energy):
def __init__(self, lamb, d):
"""
lamb: Model object
value = 0.5 * s.vdot(s), i.e. a log-Gauss distribution with unit
covariance
"""
super(PoissonianEnergy, self).__init__(lamb.position)
self._lamb = lamb
class PoissonianEnergy(Operator):
def __init__(self, op, d):
super(PoissonianEnergy, self).__init__()
self._op = op
self._d = d
lamb_val = self._lamb.value
self._value = lamb_val.sum() - d.vdot(log(lamb_val))
if isnan(self._value):
self._value = inf
self._gradient = self._lamb.jacobian.adjoint_times(1 - d/lamb_val)
metric = makeOp(1./lamb_val)
self._metric = SandwichOperator.make(self._lamb.jacobian, metric)
def at(self, position):
return self.__class__(self._lamb.at(position), self._d)
@property
def value(self):
return self._value
@property
def gradient(self):
return self._gradient
@property
def metric(self):
return self._metric
def __call__(self, x):
x = self._op(x)
res = (x - self._d*x.log()).sum()
metric = SandwichOperator.make(x.jac, makeOp(1./x.val))
return res.add_metric(metric)
......@@ -32,7 +32,7 @@ class Energy_Tests(unittest.TestCase):
@expand(product([ift.GLSpace(15),
ift.RGSpace(64, distances=.789),
ift.RGSpace([32, 32], distances=.789)],
[ift.Tanh, ift.Exponential, ift.Linear],
["tanh", "exp", ""],
[1, 1e-2, 1e2],
[4, 78, 23]))
def testGaussianEnergy(self, space, nonlinearity, noise, seed):
......@@ -44,7 +44,7 @@ class Energy_Tests(unittest.TestCase):
pspace = ift.PowerSpace(hspace, binbounds=binbounds)
Dist = ift.PowerDistributor(target=hspace, power_space=pspace)
xi0 = ift.Field.from_random(domain=hspace, random_type='normal')
xi0_var = ift.Variable(ift.MultiField.from_dict({'xi': xi0}))['xi']
xi0_var = ift.Linearization.make_var(xi0)
def pspec(k): return 1 / (1 + k**2)**dim
pspec = ift.PS_field(pspace, pspec)
......@@ -53,16 +53,22 @@ class Energy_Tests(unittest.TestCase):
n = N.draw_sample()
s = ht(ift.makeOp(A)(xi0_var))
R = ift.ScalingOperator(10., space)
d_model = R(ift.LocalModel(s, nonlinearity()))
d = d_model.value + n
def d_model(inp):
if nonlinearity == "":
return R(ht(ift.makeOp(A)(inp)))
else:
tmp = ht(ift.makeOp(A)(inp))
nonlin = getattr(tmp, nonlinearity)
return R(nonlin())
d = d_model(xi0) + n
if noise == 1:
N = None
energy = ift.GaussianEnergy(d_model, d, N)
if isinstance(nonlinearity(), ift.Linear):
ift.extra.check_value_gradient_metric_consistency(
energy, ntries=10)
energy = lambda inp: ift.GaussianEnergy(d, N)(d_model(inp))
if nonlinearity == "":
ift.extra.check_value_gradient_metric_consistency2(
energy, xi0, ntries=10)
else:
ift.extra.check_value_gradient_consistency(
energy, ntries=10)
ift.extra.check_value_gradient_consistency2(
energy, xi0, ntries=10)
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