Commit e4cac71d authored by Martin Reinecke's avatar Martin Reinecke

fix Bernoulli

parent f54038f7
......@@ -41,8 +41,7 @@ if __name__ == '__main__':
harmonic_space = position_space.get_default_codomain()
HT = ift.HarmonicTransformOperator(harmonic_space, position_space)
domain = ift.MultiDomain.make({'xi': harmonic_space})
position = ift.from_random('normal', domain)
position = ift.from_random('normal', harmonic_space)
# Define power spectrum and amplitudes
def sqrtpspec(k):
......@@ -54,10 +53,7 @@ if __name__ == '__main__':
A = pd(a)
# Set up a sky model
xi = ift.Variable(position)['xi']
logsky_h = xi * A
logsky = HT(logsky_h)
sky = ift.PointwisePositiveTanh(logsky)
sky = lambda inp: HT(A*inp).positive_tanh()
GR = ift.GeometryRemover(position_space)
# Set up instrumental response
......@@ -65,14 +61,14 @@ if __name__ == '__main__':
# Generate mock data
d_space = R.target[0]
p = R(sky)
mock_position = ift.from_random('normal', p.position.domain)
pp = p.at(mock_position).value
p = lambda inp: R(sky(inp))
mock_position = ift.from_random('normal', harmonic_space)
pp = p(mock_position)
data = np.random.binomial(1, pp.to_global_data().astype(np.float64))
data = ift.Field.from_global_data(d_space, data)
# Compute likelihood and Hamiltonian
position = ift.from_random('normal', p.position.domain)
position = ift.from_random('normal', harmonic_space)
likelihood = ift.BernoulliEnergy(p, data)
ic_cg = ift.GradientNormController(iteration_limit=50)
ic_newton = ift.GradientNormController(name='Newton', iteration_limit=30,
......@@ -82,14 +78,15 @@ if __name__ == '__main__':
# Minimize the Hamiltonian
H = ift.Hamiltonian(likelihood, ic_sampling)
H = ift.EnergyAdapter(position, H)
H = H.make_invertible(ic_cg)
# minimizer = ift.SteepestDescent(ic_newton)
H, convergence = minimizer(H)
reconstruction = sky.at(H.position).value
reconstruction = sky(H.position)
ift.plot(reconstruction, title='reconstruction')
ift.plot(GR.adjoint_times(data), title='data')
ift.plot(sky.at(mock_position).value, title='truth')
ift.plot(sky(mock_position), title='truth')
ift.plot_finish(nx=3, xsize=16, ysize=5, title="results",
name="bernoulli.png")
......@@ -18,45 +18,21 @@
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 BernoulliEnergy(Energy):
class BernoulliEnergy(Operator):
def __init__(self, p, d):
"""
p: Model object
"""
super(BernoulliEnergy, self).__init__(p.position)
super(BernoulliEnergy, self).__init__()
self._p = p
self._d = d
p_val = self._p.value
self._value = -self._d.vdot(log(p_val)) - (1.-d).vdot(log(1.-p_val))
if isnan(self._value):
self._value = inf
metric = makeOp(1. / (p_val * (1.-p_val)))
self._gradient = self._p.jacobian.adjoint_times(metric(p_val-d))
self._metric = SandwichOperator.make(self._p.jacobian, metric)
def at(self, position):
return self.__class__(self._p.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._p(x)
v = ((-self._d)*x.log()).sum() - ((1.-self._d)*((1.-x).log())).sum()
met = makeOp(1./(x.val*(1.-x.val)))
met = SandwichOperator.make(x.jac, met)
return v.add_metric(met)
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