Commit 8b623ab3 authored by Philipp Arras's avatar Philipp Arras

Merge branch 'NIFTy_5' into amplitudes_revisited

parents 0e8af4d4 14e3add1
......@@ -34,22 +34,16 @@ if __name__ == '__main__':
# Setting up an amplitude model
A = ift.AmplitudeModel(position_space, 64, 3, 0.4, -4., 1, 1., 1.)
dummy = ift.from_random('normal', A.domain)
# Building the model for a correlated signal
harmonic_space = position_space.get_default_codomain()
ht = ift.HarmonicTransformOperator(harmonic_space, position_space)
power_space = A.target[0]
power_distributor = ift.PowerDistributor(harmonic_space, power_space)
dummy = ift.Field.from_random('normal', harmonic_space)
domain = ift.MultiDomain.union((A.domain,
ift.MultiDomain.make({
'xi': harmonic_space
})))
vol = harmonic_space.scalar_dvol
vol = ift.ScalingOperator(vol ** (-0.5),harmonic_space)
correlated_field = ht(vol(power_distributor(A))*ift.FieldAdapter(domain, "xi"))
correlated_field = ht(vol(power_distributor(A))*ift.FieldAdapter(harmonic_space, "xi"))
# alternatively to the block above one can do:
#correlated_field = ift.CorrelatedField(position_space, A)
......@@ -67,7 +61,7 @@ if __name__ == '__main__':
N = ift.ScalingOperator(noise, data_space)
# generate mock data
MOCK_POSITION = ift.from_random('normal', domain)
MOCK_POSITION = ift.from_random('normal', signal_response.domain)
data = signal_response(MOCK_POSITION) + N.draw_sample()
# set up model likelihood
......@@ -82,13 +76,13 @@ if __name__ == '__main__':
# build model Hamiltonian
H = ift.Hamiltonian(likelihood, ic_sampling)
INITIAL_POSITION = ift.full(H.domain, 0.)
INITIAL_POSITION = ift.from_random('normal', H.domain)
position = INITIAL_POSITION
plot = ift.Plot()
plot.add(signal(MOCK_POSITION), title='Ground Truth')
plot.add(R.adjoint_times(data), title='Data')
plot.add([A(MOCK_POSITION.extract(A.domain))], title='Power Spectrum')
plot.add([A.force(MOCK_POSITION)], title='Power Spectrum')
plot.output(ny=1, nx=3, xsize=24, ysize=6, name="setup.png")
# number of samples used to estimate the KL
......@@ -100,14 +94,10 @@ if __name__ == '__main__':
plot = ift.Plot()
plot.add(signal(KL.position), title="reconstruction")
plot.add(
[
A(KL.position.extract(A.domain)),
A(MOCK_POSITION.extract(A.domain))
],
title="power")
plot.add([A.force(KL.position), A.force(MOCK_POSITION)], title="power")
plot.output(ny=1, ysize=6, xsize=16, name="loop.png")
KL = ift.KL_Energy(position, H, N_samples)
plot = ift.Plot()
sc = ift.StatCalculator()
for sample in KL.samples:
......@@ -115,9 +105,8 @@ if __name__ == '__main__':
plot.add(sc.mean, title="Posterior Mean")
plot.add(ift.sqrt(sc.var), title="Posterior Standard Deviation")
powers = [A((s + KL.position).extract(A.domain)) for s in KL.samples]
powers = [A.force(s + KL.position) for s in KL.samples]
plot.add(
[A(KL.position.extract(A.domain)),
A(MOCK_POSITION.extract(A.domain))] + powers,
[A.force(KL.position), A.force(MOCK_POSITION)] + powers,
title="Sampled Posterior Power Spectrum")
plot.output(ny=1, nx=3, xsize=24, ysize=6, name="results.png")
......@@ -73,13 +73,14 @@ from .minimization.kl_energy import KL_Energy
from .sugar import *
from .plot import Plot
from .library.amplitude_model import AmplitudeModel
from .library.amplitude_model import AmplitudeModel, SmoothAmplitudeModel
from .library.inverse_gamma_model import InverseGammaModel
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 .library.adjust_variances import (make_adjust_variances,
do_adjust_variances)
from . import extra
......
......@@ -28,7 +28,21 @@ from .structured_domain import StructuredDomain
class LogRGSpace(StructuredDomain):
"""NIFTy subclass for logarithmic Cartesian grids.
Parameters
----------
shape : int or tuple of int
Number of grid points or numbers of gridpoints along each axis.
bindistances : float or tuple of float
Distance between two grid points along each axis. These are
measured on logarithmic scale and are constant therfore.
t_0 : float or tuple of float
FIXME
harmonic : bool, optional
Whether the space represents a grid in position or harmonic space.
(default: False).
"""
_needed_for_hash = ['_shape', '_bindistances', '_t_0', '_harmonic']
def __init__(self, shape, bindistances, t_0, harmonic=False):
......@@ -41,8 +55,8 @@ class LogRGSpace(StructuredDomain):
self._bindistances = tuple(bindistances)
self._t_0 = tuple(t_0)
self._dim = int(reduce(lambda x, y: x * y, self._shape))
self._dvol = float(reduce(lambda x, y: x * y, self._bindistances))
self._dim = int(reduce(lambda x, y: x*y, self._shape))
self._dvol = float(reduce(lambda x, y: x*y, self._bindistances))
@property
def harmonic(self):
......@@ -69,13 +83,12 @@ class LogRGSpace(StructuredDomain):
return np.array(self._t_0)
def __repr__(self):
return ("LogRGSpace(shape={}, harmonic={})"
.format(self.shape, self.harmonic))
return ("LogRGSpace(shape={}, harmonic={})".format(
self.shape, self.harmonic))
def get_default_codomain(self):
codomain_bindistances = 1. / (self.bindistances * self.shape)
return LogRGSpace(self.shape, codomain_bindistances,
self._t_0, True)
codomain_bindistances = 1./(self.bindistances*self.shape)
return LogRGSpace(self.shape, codomain_bindistances, self._t_0, True)
def get_k_length_array(self):
if not self.harmonic:
......@@ -97,4 +110,3 @@ class LogRGSpace(StructuredDomain):
ks[1:] += self.t_0[i]
k_array[i] += ks.reshape((1,)*i + (self.shape[i],) + (1,)*(ndim-i-1))
return k_array
......@@ -20,7 +20,7 @@ if _use_fftw:
_fft_extra_args = dict(planner_effort='FFTW_ESTIMATE', threads=nthreads)
else:
from numpy.fft import fftn, rfftn, ifftn
_fft_extra_args={}
_fft_extra_args = {}
def hartley(a, axes=None):
......
......@@ -19,11 +19,20 @@
from __future__ import absolute_import, division, print_function
from ..compat import *
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.scaling_operator import ScalingOperator
from ..operators.simple_linear_operators import FieldAdapter
def make_adjust_variances(a, xi, position, samples=[], scaling=None,
def make_adjust_variances(a,
xi,
position,
samples=[],
scaling=None,
ic_samp=None):
""" Creates a Hamiltonian for constant likelihood optimizations.
......@@ -57,13 +66,48 @@ def make_adjust_variances(a, xi, position, samples=[], scaling=None,
if n > 0:
d_eval = 0.
for i in range(n):
d_eval = d_eval + d(position + samples[i])
d_eval = d_eval + d.force(position + samples[i])
d_eval = d_eval/n
else:
d_eval = d(position)
d_eval = d.force(position)
x = (a.conjugate()*a).real
if scaling is not None:
x = ScalingOperator(scaling, x.target)(x)
return Hamiltonian(InverseGammaLikelihood(d_eval)(x), ic_samp=ic_samp)
def do_adjust_variances(position,
amplitude_model,
minimizer,
xi_key='xi',
samples=[]):
h_space = position[xi_key].domain[0]
pd = PowerDistributor(h_space, amplitude_model.target[0])
a = pd(amplitude_model)
xi = FieldAdapter(h_space, xi_key)
ham = make_adjust_variances(a, xi, position, samples=samples)
# Minimize
e = EnergyAdapter(
position.extract(a.domain), ham, constants=[], want_metric=True)
e, _ = minimizer(e)
# Update position
s_h_old = (a*xi).force(position)
position = position.to_dict()
position['xi'] = s_h_old/a(e.position)
position = MultiField.from_dict(position)
position = MultiField.union([position, e.position])
s_h_new = (a*xi).force(position)
import numpy as np
# TODO Move this into the tests
np.testing.assert_allclose(s_h_new.to_global_data(),
s_h_old.to_global_data())
return position
......@@ -143,3 +143,54 @@ class AmplitudeModel(Operator):
@property
def norm_phi_mean(self):
return self._norm_phi_mean
class SmoothAmplitudeModel(Operator):
'''
Computes a smooth power spectrum.
Output lives in PowerSpace.
Parameters
----------
Npixdof : #pix in dof_space
ceps_a, ceps_k0 : Smoothness parameters in ceps_kernel
eg. ceps_kernel(k) = (a/(1+(k/k0)**2))**2
a = ceps_a, k0 = ceps_k0
'''
def __init__(self, s_space, Npixdof, ceps_a, ceps_k, name='tau'):
from ..operators.exp_transform import ExpTransform
from ..operators.qht_operator import QHTOperator
h_space = s_space.get_default_codomain()
self._exp_transform = ExpTransform(PowerSpace(h_space), Npixdof)
logk_space = self._exp_transform.domain[0]
qht = QHTOperator(target=logk_space)
dof_space = qht.domain[0]
self._domain = MultiDomain.make({name: dof_space})
self._target = self._exp_transform.target
kern = lambda k: _ceps_kernel(dof_space, k, ceps_a, ceps_k)
cepstrum = create_cepstrum_amplitude_field(dof_space, kern)
self._smooth_op = qht(makeOp(sqrt(cepstrum)))
self._key = name
self._qht = qht
self._ceps = makeOp(sqrt(cepstrum))
def apply(self, x):
self._check_input(x)
smooth_spec = self._smooth_op(x[self._key])
loglog_spec = smooth_spec
return self._exp_transform((0.5*loglog_spec).exp())
@property
def qht(self):
return self._qht
@property
def ceps(self):
return self._ceps
......@@ -48,7 +48,7 @@ def CorrelatedField(s_space, amplitude_model, name='xi'):
A = power_distributor(amplitude_model)
vol = h_space.scalar_dvol
vol = ScalingOperator(vol ** (-0.5),h_space)
return ht(vol(A)*FieldAdapter(MultiDomain.make({name: h_space}), name))
return ht(vol(A)*FieldAdapter(h_space, name))
def MfCorrelatedField(s_space_spatial, s_space_energy, amplitude_model_spatial,
......@@ -77,4 +77,4 @@ def MfCorrelatedField(s_space_spatial, s_space_energy, amplitude_model_spatial,
a_energy = dom_distr_energy(amplitude_model_energy)
a = a_spatial*a_energy
A = pd(a)
return ht(A*FieldAdapter(MultiDomain.make({name: h_space}), name))
return ht(A*FieldAdapter(h_space, name))
......@@ -52,7 +52,7 @@ class Linearization(object):
def __getitem__(self, name):
from .operators.simple_linear_operators import FieldAdapter
return self.new(self._val[name], FieldAdapter(self.domain, name))
return self.new(self._val[name], FieldAdapter(self.domain[name], name))
def __neg__(self):
return self.new(-self._val, -self._jac,
......
......@@ -41,8 +41,8 @@ class FieldZeroPadder(LinearOperator):
if len(new_shape) != len(dom.shape):
raise ValueError("Shape mismatch")
if any([a <= b for a, b in zip(new_shape, dom.shape)]):
raise ValueError("New shape must be larger than old shape")
if any([a < b for a, b in zip(new_shape, dom.shape)]):
raise ValueError("New shape must not be smaller than old shape")
self._target = list(self._domain)
self._target[self._space] = RGSpace(new_shape, dom.distances)
self._target = DomainTuple.make(self._target)
......@@ -54,6 +54,9 @@ class FieldZeroPadder(LinearOperator):
curshp = list(self._dom(mode).shape)
tgtshp = self._tgt(mode).shape
for d in self._target.axes[self._space]:
if v.shape[d] == tgtshp[d]: # nothing to do
continue
idx = (slice(None),) * d
v, x = dobj.ensure_not_distributed(v, (d,))
......
......@@ -59,6 +59,9 @@ class Operator(NiftyMetaBase()):
def apply(self, x):
raise NotImplementedError
def force(self, x):
return self.apply(x.extract(self.domain))
def _check_input(self, x):
from ..linearization import Linearization
d = x.target if isinstance(x, Linearization) else x.domain
......
......@@ -61,4 +61,3 @@ class OuterProduct(LinearOperator):
return Field.from_global_data(
self._domain, np.tensordot(
self._field.to_global_data(), x.to_global_data(), axes))
......@@ -64,14 +64,24 @@ class Realizer(EndomorphicOperator):
class FieldAdapter(LinearOperator):
def __init__(self, dom, name):
self._target = dom[name]
"""Operator which extracts a field out of a MultiField.
Parameters
----------
target : Domain, tuple of Domain or DomainTuple:
The domain which shall be extracted out of the MultiField.
name : String
The key of the MultiField which shall be extracted.
"""
def __init__(self, target, name):
self._target = DomainTuple.make(target)
self._domain = MultiDomain.make({name: self._target})
self._capability = self.TIMES | self.ADJOINT_TIMES
def apply(self, x, mode):
self._check_input(x, mode)
if mode == self.TIMES:
return x.values()[0]
return MultiField(self._domain, (x,))
......
......@@ -62,21 +62,21 @@ class Model_Tests(unittest.TestCase):
lin2 = self.make_linearization(type2, dom2, seed)
dom = ift.MultiDomain.union((dom1, dom2))
model = ift.FieldAdapter(dom, "s1")*ift.FieldAdapter(dom, "s2")
model = ift.FieldAdapter(space, "s1")*ift.FieldAdapter(space, "s2")
pos = ift.from_random("normal", dom)
ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
model = ift.FieldAdapter(dom, "s1")+ift.FieldAdapter(dom, "s2")
model = ift.FieldAdapter(space, "s1")+ift.FieldAdapter(space, "s2")
pos = ift.from_random("normal", dom)
ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
model = ift.FieldAdapter(dom, "s1").scale(3.)
model = ift.FieldAdapter(space, "s1").scale(3.)
pos = ift.from_random("normal", dom1)
ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
model = ift.ScalingOperator(2.456, space)(
ift.FieldAdapter(dom, "s1")*ift.FieldAdapter(dom, "s2"))
ift.FieldAdapter(space, "s1")*ift.FieldAdapter(space, "s2"))
pos = ift.from_random("normal", dom)
ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
model = ift.positive_tanh(ift.ScalingOperator(2.456, space)(
ift.FieldAdapter(dom, "s1")*ift.FieldAdapter(dom, "s2")))
ift.FieldAdapter(space, "s1")*ift.FieldAdapter(space, "s2")))
pos = ift.from_random("normal", dom)
ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
pos = ift.from_random("normal", dom)
......@@ -84,7 +84,7 @@ class Model_Tests(unittest.TestCase):
ift.extra.check_value_gradient_consistency(model, pos['s2'], ntries=20)
if isinstance(space, ift.RGSpace):
model = ift.FFTOperator(space)(
ift.FieldAdapter(dom, "s1")*ift.FieldAdapter(dom, "s2"))
ift.FieldAdapter(space, "s1")*ift.FieldAdapter(space, "s2"))
pos = ift.from_random("normal", dom)
ift.extra.check_value_gradient_consistency(model, pos, ntries=20)
......
......@@ -208,7 +208,7 @@ class Consistency_Tests(unittest.TestCase):
op = ift.SymmetrizingOperator(dom, space)
ift.extra.consistency_check(op, dtype, dtype)
@expand(product([0, 2], [2, 2.7], [np.float64, np.complex128],
@expand(product([0, 2], [1, 2, 2.7], [np.float64, np.complex128],
[False, True]))
def testZeroPadder(self, space, factor, dtype, central):
dom = (ift.RGSpace(10), ift.UnstructuredDomain(13), ift.RGSpace(7, 12),
......
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