Commit 1eecc5ee authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'NIFTy_5' into goodbye_blue_sky

parents df3bdf4c fe28eda9
...@@ -79,7 +79,8 @@ if __name__ == '__main__': ...@@ -79,7 +79,8 @@ if __name__ == '__main__':
minimizer = ift.RelaxedNewton(ic_newton) minimizer = ift.RelaxedNewton(ic_newton)
# Minimize the Hamiltonian # Minimize the Hamiltonian
H = ift.Hamiltonian(likelihood, ic_cg) H = ift.Hamiltonian(likelihood)
H = H.makeInvertible(ic_cg)
H, convergence = minimizer(H) H, convergence = minimizer(H)
# Plot results # Plot results
......
...@@ -61,8 +61,7 @@ if __name__ == '__main__': ...@@ -61,8 +61,7 @@ if __name__ == '__main__':
minimizer = ift.RelaxedNewton(ic_newton) minimizer = ift.RelaxedNewton(ic_newton)
# build model Hamiltonian # build model Hamiltonian
H = ift.Hamiltonian(likelihood, ic_cg, H = ift.Hamiltonian(likelihood, ic_sampling)
iteration_controller_sampling=ic_sampling)
INITIAL_POSITION = ift.from_random('normal', H.position.domain) INITIAL_POSITION = ift.from_random('normal', H.position.domain)
position = INITIAL_POSITION position = INITIAL_POSITION
...@@ -78,7 +77,8 @@ if __name__ == '__main__': ...@@ -78,7 +77,8 @@ if __name__ == '__main__':
samples = [H.curvature.draw_sample(from_inverse=True) samples = [H.curvature.draw_sample(from_inverse=True)
for _ in range(N_samples)] for _ in range(N_samples)]
KL = ift.SampledKullbachLeiblerDivergence(H, samples, ic_cg) KL = ift.SampledKullbachLeiblerDivergence(H, samples)
KL = KL.makeInvertible(ic_cg)
KL, convergence = minimizer(KL) KL, convergence = minimizer(KL)
position = KL.position position = KL.position
......
...@@ -19,26 +19,23 @@ ...@@ -19,26 +19,23 @@
from ..library.gaussian_energy import GaussianEnergy from ..library.gaussian_energy import GaussianEnergy
from ..minimization.energy import Energy from ..minimization.energy import Energy
from ..models.variable import Variable from ..models.variable import Variable
from ..operators import InversionEnabler, SamplingEnabler from ..operators.sampling_enabler import SamplingEnabler
from ..utilities import memo from ..utilities import memo
class Hamiltonian(Energy): class Hamiltonian(Energy):
def __init__(self, lh, iteration_controller, def __init__(self, lh, iteration_controller_sampling=None):
iteration_controller_sampling=None):
""" """
lh: Likelihood (energy object) lh: Likelihood (energy object)
prior: prior:
""" """
super(Hamiltonian, self).__init__(lh.position) super(Hamiltonian, self).__init__(lh.position)
self._lh = lh self._lh = lh
self._ic = iteration_controller
self._ic_samp = iteration_controller_sampling self._ic_samp = iteration_controller_sampling
self._prior = GaussianEnergy(Variable(self.position)) self._prior = GaussianEnergy(Variable(self.position))
self._precond = self._prior.curvature
def at(self, position): def at(self, position):
return self.__class__(self._lh.at(position), self._ic, self._ic_samp) return self.__class__(self._lh.at(position), self._ic_samp)
@property @property
@memo @memo
...@@ -55,11 +52,10 @@ class Hamiltonian(Energy): ...@@ -55,11 +52,10 @@ class Hamiltonian(Energy):
def curvature(self): def curvature(self):
prior_curv = self._prior.curvature prior_curv = self._prior.curvature
if self._ic_samp is None: if self._ic_samp is None:
c = self._lh.curvature + prior_curv return self._lh.curvature + prior_curv
else: else:
c = SamplingEnabler(self._lh.curvature, prior_curv.inverse, return SamplingEnabler(self._lh.curvature, prior_curv.inverse,
self._ic_samp, prior_curv.inverse) self._ic_samp, prior_curv.inverse)
return InversionEnabler(c, self._ic, self._precond)
def __str__(self): def __str__(self):
res = 'Likelihood:\t{:.2E}\n'.format(self._lh.value) res = 'Likelihood:\t{:.2E}\n'.format(self._lh.value)
......
from builtins import *
from ..minimization.energy import Energy from ..minimization.energy import Energy
from ..operators.inversion_enabler import InversionEnabler from ..utilities import memo, my_sum
from ..operators.scaling_operator import ScalingOperator
from ..utilities import memo
class SampledKullbachLeiblerDivergence(Energy): class SampledKullbachLeiblerDivergence(Energy):
def __init__(self, h, res_samples, iteration_controller): def __init__(self, h, res_samples):
""" """
h: Hamiltonian h: Hamiltonian
N: Number of samples to be used N: Number of samples to be used
...@@ -13,41 +12,27 @@ class SampledKullbachLeiblerDivergence(Energy): ...@@ -13,41 +12,27 @@ class SampledKullbachLeiblerDivergence(Energy):
super(SampledKullbachLeiblerDivergence, self).__init__(h.position) super(SampledKullbachLeiblerDivergence, self).__init__(h.position)
self._h = h self._h = h
self._res_samples = res_samples self._res_samples = res_samples
self._iteration_controller = iteration_controller
self._energy_list = [] self._energy_list = tuple(h.at(self.position+ss)
for ss in res_samples: for ss in res_samples)
e = h.at(self.position+ss)
self._energy_list.append(e)
def at(self, position): def at(self, position):
return self.__class__(self._h.at(position), self._res_samples, return self.__class__(self._h.at(position), self._res_samples)
self._iteration_controller)
@property @property
@memo @memo
def value(self): def value(self):
v = self._energy_list[0].value return (my_sum(map(lambda v: v.value, self._energy_list)) /
for energy in self._energy_list[1:]: len(self._energy_list))
v += energy.value
return v / len(self._energy_list)
@property @property
@memo @memo
def gradient(self): def gradient(self):
g = self._energy_list[0].gradient return (my_sum(map(lambda v: v.gradient, self._energy_list)) /
for energy in self._energy_list[1:]: len(self._energy_list))
g += energy.gradient
return g / len(self._energy_list)
@property @property
@memo @memo
def curvature(self): def curvature(self):
# MR FIXME: This looks a bit strange... return (my_sum(map(lambda v: v.curvature, self._energy_list)) *
approx = self._energy_list[-1]._prior.curvature (1./len(self._energy_list)))
curvature_list = [e.curvature for e in self._energy_list]
op = curvature_list[0]
for curv in curvature_list[1:]:
op = op + curv
op = op * ScalingOperator(1./len(curvature_list), op.domain)
return InversionEnabler(op, self._iteration_controller, approx)
import numpy as np import numpy as np
from ..domains import PowerSpace, UnstructuredDomain from ..domains.power_space import PowerSpace
from ..domains.unstructured_domain import UnstructuredDomain
from ..field import Field from ..field import Field
from ..multi import MultiField from ..multi.multi_field import MultiField
from ..sugar import makeOp, sqrt from ..sugar import makeOp, sqrt
...@@ -26,7 +27,7 @@ def make_amplitude_model(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv, ...@@ -26,7 +27,7 @@ def make_amplitude_model(s_space, Npixdof, ceps_a, ceps_k, sm, sv, im, iv,
eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2 eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
a = ceps_a, k0 = ceps_k0 a = ceps_a, k0 = ceps_k0
sm, sv : slope_mean = expected exponent of powerlaw (e.g. -4), sm, sv : slope_mean = expected exponent of power law (e.g. -4),
slope_variance (default=1) slope_variance (default=1)
im, iv : y-intercept_mean, y-intercept_variance of power_slope im, iv : y-intercept_mean, y-intercept_variance of power_slope
......
...@@ -80,12 +80,9 @@ class DescentMinimizer(Minimizer): ...@@ -80,12 +80,9 @@ class DescentMinimizer(Minimizer):
return energy, controller.CONVERGED return energy, controller.CONVERGED
# compute a step length that reduces energy.value sufficiently # compute a step length that reduces energy.value sufficiently
try: new_energy, success = self.line_searcher.perform_line_search(
new_energy, success = self.line_searcher.perform_line_search( energy=energy, pk=self.get_descent_direction(energy),
energy=energy, pk=self.get_descent_direction(energy), f_k_minus_1=f_k_minus_1)
f_k_minus_1=f_k_minus_1)
except ValueError:
return energy, controller.ERROR
if not success: if not success:
self.reset() self.reset()
......
...@@ -129,6 +129,12 @@ class Energy(NiftyMetaBase()): ...@@ -129,6 +129,12 @@ class Energy(NiftyMetaBase()):
""" """
return None return None
def makeInvertible(self, controller, preconditioner=None):
from .iteration_controller import IterationController
if not isinstance(controller, IterationController):
raise TypeError
return CurvatureInversionEnabler(self, controller, preconditioner)
def __mul__(self, factor): def __mul__(self, factor):
from .energy_sum import EnergySum from .energy_sum import EnergySum
if isinstance(factor, (float, int)): if isinstance(factor, (float, int)):
...@@ -153,3 +159,45 @@ class Energy(NiftyMetaBase()): ...@@ -153,3 +159,45 @@ class Energy(NiftyMetaBase()):
def __neg__(self): def __neg__(self):
from .energy_sum import EnergySum from .energy_sum import EnergySum
return EnergySum.make([self], [-1.]) return EnergySum.make([self], [-1.])
class CurvatureInversionEnabler(Energy):
def __init__(self, ene, controller, preconditioner):
super(CurvatureInversionEnabler, self).__init__(ene.position)
self._energy = ene
self._controller = controller
self._preconditioner = preconditioner
def at(self, position):
if self._position.isSubsetOf(position):
return self
return CurvatureInversionEnabler(
self._energy.at(position), self._controller, self._preconditioner)
@property
def position(self):
return self._energy.position
@property
def value(self):
return self._energy.value
@property
def gradient(self):
return self._energy.gradient
@property
def curvature(self):
from ..operators.linear_operator import LinearOperator
from ..operators.inversion_enabler import InversionEnabler
curv = self._energy.curvature
if self._preconditioner is None:
precond = None
elif isinstance(self._preconditioner, LinearOperator):
precond = self._preconditioner
elif isinstance(self._preconditioner, Energy):
precond = self._preconditioner.at(self.position).curvature
return InversionEnabler(curv, self._controller, precond)
def longest_step(self, dir):
return self._energy.longest_step(dir)
...@@ -16,7 +16,8 @@ ...@@ -16,7 +16,8 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes. # and financially supported by the Studienstiftung des deutschen Volkes.
from ..utilities import memo from builtins import *
from ..utilities import memo, my_lincomb_simple, my_lincomb
from .energy import Energy from .energy import Energy
...@@ -55,26 +56,17 @@ class EnergySum(Energy): ...@@ -55,26 +56,17 @@ class EnergySum(Energy):
@property @property
@memo @memo
def value(self): def value(self):
res = self._energies[0].value * self._factors[0] return my_lincomb_simple(map(lambda v: v.value, self._energies),
for e, f in zip(self._energies[1:], self._factors[1:]): self._factors)
res += e.value * f
return res
@property @property
@memo @memo
def gradient(self): def gradient(self):
res = self._energies[0].gradient.copy() if self._factors[0] == 1. \ return my_lincomb(map(lambda v: v.gradient, self._energies),
else self._energies[0].gradient * self._factors[0] self._factors).lock()
for e, f in zip(self._energies[1:], self._factors[1:]):
res += e.gradient if f == 1. else f*e.gradient
return res.lock()
@property @property
@memo @memo
def curvature(self): def curvature(self):
res = self._energies[0].curvature if self._factors[0] == 1. \ return my_lincomb(map(lambda v: v.curvature, self._energies),
else self._energies[0].curvature * self._factors[0] self._factors)
for e, f in zip(self._energies[1:], self._factors[1:]):
res = res + (e.curvature if f == 1. else e.curvature*f)
return res
...@@ -274,9 +274,9 @@ class LinearOperator(NiftyMetaBase()): ...@@ -274,9 +274,9 @@ class LinearOperator(NiftyMetaBase()):
def _check_mode(self, mode): def _check_mode(self, mode):
if not self._validMode[mode]: if not self._validMode[mode]:
raise ValueError("invalid operator mode specified") raise NotImplementedError("invalid operator mode specified")
if mode & self.capability == 0: if mode & self.capability == 0:
raise ValueError("requested operator mode is not supported") raise NotImplementedError("requested operator mode is not supported")
def _check_input(self, x, mode): def _check_input(self, x, mode):
self._check_mode(mode) self._check_mode(mode)
......
...@@ -16,15 +16,36 @@ ...@@ -16,15 +16,36 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes. # and financially supported by the Studienstiftung des deutschen Volkes.
from builtins import next, range from builtins import *
import numpy as np import numpy as np
from itertools import product from itertools import product
import abc import abc
from future.utils import with_metaclass from future.utils import with_metaclass
from functools import reduce
__all__ = ["get_slice_list", "safe_cast", "parse_spaces", "infer_space", __all__ = ["get_slice_list", "safe_cast", "parse_spaces", "infer_space",
"memo", "NiftyMetaBase", "fft_prep", "hartley", "my_fftn_r2c", "memo", "NiftyMetaBase", "fft_prep", "hartley", "my_fftn_r2c",
"my_fftn"] "my_fftn", "my_sum", "my_lincomb_simple", "my_lincomb",
"my_product"]
def my_sum(terms):
return reduce(lambda x, y: x+y, terms)
def my_lincomb_simple(terms, factors):
terms2 = map(lambda v: v[0]*v[1], zip(terms, factors))
return my_sum(terms2)
def my_lincomb(terms, factors):
terms2 = map(lambda v: v[0] if v[1] == 1. else v[0]*v[1],
zip(terms, factors))
return my_sum(terms2)
def my_product(iterable):
return reduce(lambda x, y: x*y, iterable)
def get_slice_list(shape, axes): def get_slice_list(shape, axes):
......
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