Commit 03c31669 authored by Jakob Knollmueller's avatar Jakob Knollmueller
Browse files

renaming KL and Hamiltonian

parent 5689ba2d
......@@ -74,7 +74,7 @@ if __name__ == '__main__':
ic_sampling = ift.GradientNormController(iteration_limit=100)
# Minimize the Hamiltonian
H = ift.Hamiltonian(likelihood, ic_sampling)
H = ift.StandardHamiltonian(likelihood, ic_sampling)
H = ift.EnergyAdapter(position, H, want_metric=True)
# minimizer = ift.L_BFGS(ic_newton)
H, convergence = minimizer(H)
......@@ -99,7 +99,7 @@ if __name__ == '__main__':
minimizer = ift.NewtonCG(ic_newton)
# Compute MAP solution by minimizing the information Hamiltonian
H = ift.Hamiltonian(likelihood)
H = ift.StandardHamiltonian(likelihood)
initial_position = ift.from_random('normal', domain)
H = ift.EnergyAdapter(initial_position, H, want_metric=True)
H, convergence = minimizer(H)
......@@ -100,10 +100,10 @@ if __name__ == '__main__':
# Set up likelihood and information Hamiltonian
likelihood = ift.GaussianEnergy(mean=data, covariance=N)(signal_response)
H = ift.Hamiltonian(likelihood, ic_sampling)
H = ift.StandardHamiltonian(likelihood, ic_sampling)
initial_position = ift.MultiField.full(H.domain, 0.)
position = initial_position
initial_mean = ift.MultiField.full(H.domain, 0.)
mean = initial_mean
plot = ift.Plot()
plot.add(signal(mock_position), title='Ground Truth')
......@@ -117,9 +117,9 @@ if __name__ == '__main__':
# Draw new samples to approximate the KL five times
for i in range(5):
# Draw new samples and minimize KL
KL = ift.KL_Energy(position, H, N_samples)
KL = ift.MetricGaussianKL(mean, H, N_samples)
KL, convergence = minimizer(KL)
position = KL.position
mean = KL.position
# Plot current reconstruction
plot = ift.Plot()
......@@ -128,7 +128,7 @@ if __name__ == '__main__':
plot.output(ny=1, ysize=6, xsize=16, name="loop-{:02}.png".format(i))
# Draw posterior samples
KL = ift.KL_Energy(position, H, N_samples)
KL = ift.MetricGaussianKL(mean, H, N_samples)
sc = ift.StatCalculator()
for sample in KL.samples:
sc.add(signal(sample + KL.position))
......@@ -103,7 +103,7 @@ N = ift.DiagonalOperator(ift.from_global_data(d_space, var))
IC = ift.DeltaEnergyController(tol_rel_deltaE=1e-12, iteration_limit=200)
likelihood = ift.GaussianEnergy(d, N)(R)
Ham = ift.Hamiltonian(likelihood, IC)
Ham = ift.StandardHamiltonian(likelihood, IC)
H = ift.EnergyAdapter(params, Ham, want_metric=True)
# Minimize
......@@ -49,7 +49,7 @@ from .operators.simple_linear_operators import (
FieldAdapter, ducktape, GeometryRemover, NullOperator)
from .operators.energy_operators import (
EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood,
BernoulliEnergy, Hamiltonian, SampledKullbachLeiblerDivergence)
BernoulliEnergy, StandardHamiltonian, SampledKullbachLeiblerDivergence)
from .probing import probe_with_posterior_samples, probe_diagonal, \
......@@ -68,7 +68,7 @@ from .minimization.scipy_minimizer import (ScipyMinimizer, L_BFGS_B, ScipyCG)
from import Energy
from .minimization.quadratic_energy import QuadraticEnergy
from .minimization.energy_adapter import EnergyAdapter
from .minimization.kl_energy import KL_Energy
from .minimization.metric_gaussian_kl import MetricGaussianKL
from .sugar import *
from .plot import Plot
......@@ -19,7 +19,7 @@ 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.energy_operators import StandardHamiltonian, InverseGammaLikelihood
from ..operators.scaling_operator import ScalingOperator
from ..operators.simple_linear_operators import ducktape
......@@ -52,7 +52,7 @@ def make_adjust_variances(a,
A Hamiltonian that can be used for further minimization
......@@ -71,7 +71,7 @@ def make_adjust_variances(a,
if scaling is not None:
x = ScalingOperator(scaling,
return Hamiltonian(InverseGammaLikelihood(d_eval)(x), ic_samp=ic_samp)
return StandardHamiltonian(InverseGammaLikelihood(d_eval)(x), ic_samp=ic_samp)
def do_adjust_variances(position,
......@@ -20,31 +20,68 @@ from ..linearization import Linearization
from .. import utilities
class KL_Energy(Energy):
def __init__(self, position, h, nsamp, constants=[],
constants_samples=None, gen_mirrored_samples=False,
class MetricGaussianKL(Energy):
"""Provides the sampled Kullback-Leibler divergence between a distribution and a metric Gaussian.
The Energy object is an implementation of a scalar function including its
gradient and metric at some position.
mean : Field
The current mean of the Gaussian.
hamiltonian : Hamiltonian
The Hamiltonian of the approximated probability distribution.
n_samples : integer
The number of samples used to stochastically estimate the KL.
constants : list
A list of parameter keys that are kept constant during optimization.
point_estimates : list
A list of parameter keys for which no samples are drawn, but that are
optimized for, corresponding to point estimates of these.
mirror_samples : boolean
Whether the negative of the drawn samples are also used,
as they are equaly legitimate samples. If true, the number of used
samples doubles. Mirroring samples stabilizes the KL estimate as
extreme sample variation is counterbalanced. (default : False)
An instance of the Energy class is defined at a certain location. If one
is interested in the value, gradient or metric of the abstract energy
functional one has to 'jump' to the new position using the `at` method.
This method returns a new energy instance residing at the new position. By
this approach, intermediate results from computing e.g. the gradient can
safely be reused for e.g. the value or the metric.
Memorizing the evaluations of some quantities minimizes the computational
effort for multiple calls.
def __init__(self, mean, hamiltonian, n_sampels, constants=[],
point_estimates=None, mirror_samples=False,
super(KL_Energy, self).__init__(position)
if h.domain is not position.domain:
super(MetricGaussianKL, self).__init__(mean)
if hamiltonian.domain is not mean.domain:
raise TypeError
self._h = h
self._hamiltonian = hamiltonian
self._constants = constants
if constants_samples is None:
constants_samples = constants
self._constants_samples = constants_samples
if point_estimates is None:
point_estimates = constants
self._constants_samples = point_estimates
if _samples is None:
met = h(Linearization.make_partial_var(
position, constants_samples, True)).metric
met = hamiltonian(Linearization.make_partial_var(
mean, point_estimates, True)).metric
_samples = tuple(met.draw_sample(from_inverse=True)
for _ in range(nsamp))
if gen_mirrored_samples:
for _ in range(n_sampels))
if mirror_samples:
_samples += tuple(-s for s in _samples)
self._samples = _samples
self._lin = Linearization.make_partial_var(position, constants)
self._lin = Linearization.make_partial_var(mean, constants)
v, g = None, None
for s in self._samples:
tmp = self._h(self._lin+s)
tmp = self._hamiltonian(self._lin+s)
if v is None:
v = tmp.val.local_data[()]
g = tmp.gradient
......@@ -56,9 +93,9 @@ class KL_Energy(Energy):
self._metric = None
def at(self, position):
return KL_Energy(position, self._h, 0,
self._constants, self._constants_samples,
return MetricGaussianKL(position, self._hamiltonian, 0,
self._constants, self._constants_samples,
def value(self):
......@@ -71,7 +108,7 @@ class KL_Energy(Energy):
def _get_metric(self):
if self._metric is None:
lin = self._lin.with_want_metric()
mymap = map(lambda v: self._h(lin+v).metric, self._samples)
mymap = map(lambda v: self._hamiltonian(lin+v).metric, self._samples)
self._metric = utilities.my_sum(mymap)
self._metric = self._metric.scale(1./len(self._samples))
......@@ -147,7 +147,7 @@ class BernoulliEnergy(EnergyOperator):
return v.add_metric(met)
class Hamiltonian(EnergyOperator):
class StandardHamiltonian(EnergyOperator):
def __init__(self, lh, ic_samp=None):
self._lh = lh
self._prior = GaussianEnergy(domain=lh.domain)
......@@ -69,7 +69,7 @@ def test_hamiltonian_and_KL(field):
field = field.exp()
space = field.domain
lh = ift.GaussianEnergy(domain=space)
hamiltonian = ift.Hamiltonian(lh)
hamiltonian = ift.StandardHamiltonian(lh)
ift.extra.check_value_gradient_consistency(hamiltonian, field)
S = ift.ScalingOperator(1., space)
samps = [S.draw_sample() for i in range(3)]
Supports Markdown
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