Commit 3996cb92 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

Merge branch 'cleanup_pa' into 'NIFTy_5'

Cleanup

See merge request ift/nifty-dev!176
parents 15e8159c 955d6828
...@@ -46,6 +46,7 @@ from .operators.outer_product_operator import OuterProduct ...@@ -46,6 +46,7 @@ from .operators.outer_product_operator import OuterProduct
from .operators.simple_linear_operators import ( from .operators.simple_linear_operators import (
VdotOperator, ConjugationOperator, Realizer, VdotOperator, ConjugationOperator, Realizer,
FieldAdapter, ducktape, GeometryRemover, NullOperator) FieldAdapter, ducktape, GeometryRemover, NullOperator)
from .operators.value_inserter import ValueInserter
from .operators.energy_operators import ( from .operators.energy_operators import (
EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood, EnergyOperator, GaussianEnergy, PoissonianEnergy, InverseGammaLikelihood,
BernoulliEnergy, Hamiltonian, AveragedEnergy) BernoulliEnergy, Hamiltonian, AveragedEnergy)
......
...@@ -34,7 +34,7 @@ def _float_or_listoffloat(inp): ...@@ -34,7 +34,7 @@ def _float_or_listoffloat(inp):
return [float(x) for x in inp] if isinstance(inp, list) else float(inp) return [float(x) for x in inp] if isinstance(inp, list) else float(inp)
def _make_dynamic_operator(domain, def _make_dynamic_operator(target,
harmonic_padding, harmonic_padding,
sm_s0, sm_s0,
sm_x0, sm_x0,
...@@ -44,8 +44,10 @@ def _make_dynamic_operator(domain, ...@@ -44,8 +44,10 @@ def _make_dynamic_operator(domain,
minimum_phase, minimum_phase,
sigc=None, sigc=None,
quant=None): quant=None):
if not isinstance(domain, RGSpace): if not isinstance(target, RGSpace):
raise TypeError("RGSpace required") raise TypeError("RGSpace required")
if not target.harmonic:
raise TypeError("Target space must be harmonic")
if not (isinstance(harmonic_padding, int) or harmonic_padding is None if not (isinstance(harmonic_padding, int) or harmonic_padding is None
or all(isinstance(ii, int) for ii in harmonic_padding)): or all(isinstance(ii, int) for ii in harmonic_padding)):
raise TypeError raise TypeError
...@@ -62,7 +64,7 @@ def _make_dynamic_operator(domain, ...@@ -62,7 +64,7 @@ def _make_dynamic_operator(domain,
if cone and (sigc is None or quant is None): if cone and (sigc is None or quant is None):
raise RuntimeError raise RuntimeError
dom = DomainTuple.make(domain) dom = DomainTuple.make(target.get_default_codomain())
ops = {} ops = {}
FFT = FFTOperator(dom) FFT = FFTOperator(dom)
Real = Realizer(dom) Real = Realizer(dom)
...@@ -134,7 +136,8 @@ def _make_dynamic_operator(domain, ...@@ -134,7 +136,8 @@ def _make_dynamic_operator(domain,
return m, ops return m, ops
def dynamic_operator(domain, def dynamic_operator(*,
target,
harmonic_padding, harmonic_padding,
sm_s0, sm_s0,
sm_x0, sm_x0,
...@@ -143,11 +146,22 @@ def dynamic_operator(domain, ...@@ -143,11 +146,22 @@ def dynamic_operator(domain,
minimum_phase=False): minimum_phase=False):
'''Constructs an operator encoding the Green's function of a linear '''Constructs an operator encoding the Green's function of a linear
homogeneous dynamic system. homogeneous dynamic system.
When evaluated, this operator returns the Green's function representation
in harmonic space. This result can be used as a convolution kernel to
construct solutions of the homogeneous stochastic differential equation
encoded in this operator. Note that if causal is True, the Green's function
is convolved with a step function in time, where the temporal axis is the
first axis of the space. In this case the resulting function only extends
up to half the length of the first axis of the space to avoid boundary
effects during convolution. If minimum_phase is true then the spectrum of
the Green's function is used to construct a corresponding minimum phase
filter.
Parameters Parameters
---------- ----------
domain : RGSpace target : RGSpace
The position space in which the Green's function shall be constructed. The harmonic space in which the Green's function shall be constructed.
harmonic_padding : None, int, list of int harmonic_padding : None, int, list of int
Amount of central padding in harmonic space in pixels. If None the Amount of central padding in harmonic space in pixels. If None the
field is not padded at all. field is not padded at all.
...@@ -159,13 +173,15 @@ def dynamic_operator(domain, ...@@ -159,13 +173,15 @@ def dynamic_operator(domain,
key for dynamics encoding parameter. key for dynamics encoding parameter.
causal : boolean causal : boolean
Whether or not the Green's function shall be causal in time. Whether or not the Green's function shall be causal in time.
Default is True.
minimum_phase: boolean minimum_phase: boolean
Whether or not the Green's function shall be a minimum phase filter. Whether or not the Green's function shall be a minimum phase filter.
Default is False.
Returns Returns
------- -------
Operator Operator
The Operator encoding the dynamic Green's function in harmonic space. The Operator encoding the dynamic Green's function in target space.
Dictionary of Operator Dictionary of Operator
A collection of sub-chains of Operator which can be used for plotting A collection of sub-chains of Operator which can be used for plotting
and evaluation. and evaluation.
...@@ -175,7 +191,7 @@ def dynamic_operator(domain, ...@@ -175,7 +191,7 @@ def dynamic_operator(domain,
The first axis of the domain is interpreted the time axis. The first axis of the domain is interpreted the time axis.
''' '''
dct = { dct = {
'domain': domain, 'target': target,
'harmonic_padding': harmonic_padding, 'harmonic_padding': harmonic_padding,
'sm_s0': sm_s0, 'sm_s0': sm_s0,
'sm_x0': sm_x0, 'sm_x0': sm_x0,
...@@ -187,7 +203,8 @@ def dynamic_operator(domain, ...@@ -187,7 +203,8 @@ def dynamic_operator(domain,
return _make_dynamic_operator(**dct) return _make_dynamic_operator(**dct)
def dynamic_lightcone_operator(domain, def dynamic_lightcone_operator(*,
target,
harmonic_padding, harmonic_padding,
sm_s0, sm_s0,
sm_x0, sm_x0,
...@@ -199,11 +216,16 @@ def dynamic_lightcone_operator(domain, ...@@ -199,11 +216,16 @@ def dynamic_lightcone_operator(domain,
minimum_phase=False): minimum_phase=False):
'''Extends the functionality of :function: dynamic_operator to a Green's '''Extends the functionality of :function: dynamic_operator to a Green's
function which is constrained to be within a light cone. function which is constrained to be within a light cone.
The resulting Green's function is constrained to be within a light cone.
This is achieved via convolution of the function with a light cone in
space-time. Thereby the first axis of the space is set to be the teporal
axis.
Parameters Parameters
---------- ----------
domain : RGSpace target : RGSpace
The position space in which the Green's function shall be constructed. The harmonic space in which the Green's function shall be constructed.
It needs to have at least two dimensions. It needs to have at least two dimensions.
harmonic_padding : None, int, list of int harmonic_padding : None, int, list of int
Amount of central padding in harmonic space in pixels. If None the Amount of central padding in harmonic space in pixels. If None the
...@@ -222,8 +244,10 @@ def dynamic_lightcone_operator(domain, ...@@ -222,8 +244,10 @@ def dynamic_lightcone_operator(domain,
Quantization of the light cone in pixels. Quantization of the light cone in pixels.
causal : boolean causal : boolean
Whether or not the Green's function shall be causal in time. Whether or not the Green's function shall be causal in time.
Default is True.
minimum_phase: boolean minimum_phase: boolean
Whether or not the Green's function shall be a minimum phase filter. Whether or not the Green's function shall be a minimum phase filter.
Default is False.
Returns Returns
------- -------
...@@ -238,10 +262,10 @@ def dynamic_lightcone_operator(domain, ...@@ -238,10 +262,10 @@ def dynamic_lightcone_operator(domain,
The first axis of the domain is interpreted the time axis. The first axis of the domain is interpreted the time axis.
''' '''
if len(domain.shape) < 2: if len(target.shape) < 2:
raise ValueError("Space must be at least 2 dimensional!") raise ValueError("Space must be at least 2 dimensional!")
dct = { dct = {
'domain': domain, 'target': target,
'harmonic_padding': harmonic_padding, 'harmonic_padding': harmonic_padding,
'sm_s0': sm_s0, 'sm_s0': sm_s0,
'sm_x0': sm_x0, 'sm_x0': sm_x0,
......
...@@ -29,35 +29,29 @@ from .linear_operator import LinearOperator ...@@ -29,35 +29,29 @@ from .linear_operator import LinearOperator
class MaskOperator(LinearOperator): class MaskOperator(LinearOperator):
"""Implementation of a mask response """Implementation of a mask response
Takes a field, applies a mask and returns the values of the field in a Takes a field, applies flags and returns the values of the field in a
UnstructuredDomain. It can be used as response operator. :class:`UnstructuredDomain`.
Parameters Parameters
---------- ----------
mask : Field flags : Field
Is converted to boolean. Where True, the input field is flagged.
""" """
def __init__(self, mask): def __init__(self, flags):
if not isinstance(mask, Field): if not isinstance(flags, Field):
raise TypeError raise TypeError
self._domain = DomainTuple.make(flags.domain)
self._domain = DomainTuple.make(mask.domain) self._flags = np.logical_not(flags.to_global_data())
self._mask = np.logical_not(mask.to_global_data()) self._target = DomainTuple.make(UnstructuredDomain(self._flags.sum()))
self._target = DomainTuple.make(UnstructuredDomain(self._mask.sum()))
self._capability = self.TIMES | self.ADJOINT_TIMES self._capability = self.TIMES | self.ADJOINT_TIMES
def data_indices(self):
if len(self.domain.shape) == 1:
return np.arange(self.domain.shape[0])[self._mask]
if len(self.domain.shape) == 2:
return np.indices(self.domain.shape).transpose((1, 2, 0))[self._mask]
def apply(self, x, mode): def apply(self, x, mode):
self._check_input(x, mode) self._check_input(x, mode)
x = x.to_global_data()
if mode == self.TIMES: if mode == self.TIMES:
res = x.to_global_data()[self._mask] res = x[self._flags]
return Field.from_global_data(self.target, res) return Field.from_global_data(self.target, res)
x = x.to_global_data()
res = np.empty(self.domain.shape, x.dtype) res = np.empty(self.domain.shape, x.dtype)
res[self._mask] = x res[self._flags] = x
res[~self._mask] = 0 res[~self._flags] = 0
return Field.from_global_data(self.domain, res) return Field.from_global_data(self.domain, res)
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
from ..domain_tuple import DomainTuple
from ..domains.unstructured_domain import UnstructuredDomain
from ..field import Field
from ..sugar import makeDomain
from .linear_operator import LinearOperator
class ValueInserter(LinearOperator):
"""Inserts one value into a field which is zero otherwise.
Parameters
----------
target : Domain, tuple of Domain or DomainTuple
index : iterable of int
The index of the target into which the value of the domain shall be
inserted.
"""
def __init__(self, target, index):
self._domain = makeDomain(UnstructuredDomain(1))
self._target = DomainTuple.make(target)
index = tuple(index)
if not all([isinstance(n, int) and n>=0 and n<self.target.shape[i] for i, n in enumerate(index)]):
raise TypeError
self._index = index
self._capability = self.TIMES | self.ADJOINT_TIMES
# Check whether index is in bounds
np.empty(self.target.shape)[self._index]
def apply(self, x, mode):
self._check_input(x, mode)
x = x.to_global_data()
if mode == self.TIMES:
res = np.zeros(self.target.shape, dtype=x.dtype)
res[self._index] = x
else:
res = np.full((1,), x[self._index], dtype=x.dtype)
return Field.from_global_data(self._tgt(mode), res)
...@@ -120,24 +120,32 @@ def testPointModel(space, seed): ...@@ -120,24 +120,32 @@ def testPointModel(space, seed):
ift.extra.check_value_gradient_consistency(model, pos, tol=1e-2, ntries=20) ift.extra.check_value_gradient_consistency(model, pos, tol=1e-2, ntries=20)
@pmp('domain', [ @pmp('target', [
ift.RGSpace(64, distances=.789), ift.RGSpace(64, distances=.789,harmonic=True),
ift.RGSpace([32, 32], distances=.789), ift.RGSpace([32, 32], distances=.789,harmonic=True),
ift.RGSpace([32, 32, 8], distances=.789) ift.RGSpace([32, 32, 8], distances=.789,harmonic=True)
]) ])
@pmp('causal', [True, False]) @pmp('causal', [True, False])
@pmp('minimum_phase', [True, False]) @pmp('minimum_phase', [True, False])
@pmp('seed', [4, 78, 23]) @pmp('seed', [4, 78, 23])
def testDynamicModel(domain, causal, minimum_phase, seed): def testDynamicModel(target, causal, minimum_phase, seed):
model, _ = ift.dynamic_operator( dct = {
domain, None, 1., 1., 'f', causal=causal, minimum_phase=minimum_phase) 'target': target,
'harmonic_padding': None,
'sm_s0': 3.,
'sm_x0': 1.,
'key': 'f',
'causal': causal,
'minimum_phase': minimum_phase
}
model, _ = ift.dynamic_operator(**dct)
S = ift.ScalingOperator(1., model.domain) S = ift.ScalingOperator(1., model.domain)
pos = S.draw_sample() pos = S.draw_sample()
# FIXME I dont know why smaller tol fails for 3D example # FIXME I dont know why smaller tol fails for 3D example
ift.extra.check_value_gradient_consistency(model, pos, tol=1e-5, ntries=20) ift.extra.check_value_gradient_consistency(model, pos, tol=1e-5, ntries=20)
if len(domain.shape) > 1: if len(target.shape) > 1:
dct = { dct = {
'domain': domain, 'target': target,
'harmonic_padding': None, 'harmonic_padding': None,
'sm_s0': 3., 'sm_s0': 3.,
'sm_x0': 1., 'sm_x0': 1.,
...@@ -148,6 +156,9 @@ def testDynamicModel(domain, causal, minimum_phase, seed): ...@@ -148,6 +156,9 @@ def testDynamicModel(domain, causal, minimum_phase, seed):
'causal': causal, 'causal': causal,
'minimum_phase': minimum_phase 'minimum_phase': minimum_phase
} }
dct['lightcone_key'] = 'c'
dct['sigc'] = 1.
dct['quant'] = 5
model, _ = ift.dynamic_lightcone_operator(**dct) model, _ = ift.dynamic_lightcone_operator(**dct)
S = ift.ScalingOperator(1., model.domain) S = ift.ScalingOperator(1., model.domain)
pos = S.draw_sample() pos = S.draw_sample()
......
...@@ -264,3 +264,17 @@ def testOuter(fdomain, domain): ...@@ -264,3 +264,17 @@ def testOuter(fdomain, domain):
f = ift.from_random('normal', fdomain) f = ift.from_random('normal', fdomain)
op = ift.OuterProduct(f, domain) op = ift.OuterProduct(f, domain)
ift.extra.consistency_check(op) ift.extra.consistency_check(op)
@pmp('sp', _h_spaces + _p_spaces + _pow_spaces)
@pmp('seed', [12, 3])
def testValueInserter(sp, seed):
np.random.seed(seed)
ind = []
for ss in sp.shape:
if ss == 1:
ind.append(0)
else:
ind.append(np.random.randint(0, ss-1))
op = ift.ValueInserter(sp, ind)
ift.extra.consistency_check(op)
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Copyright(C) 2013-2019 Max-Planck-Society
#
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
import numpy as np
import pytest
from numpy.testing import assert_
import nifty5 as ift
@pytest.mark.parametrize('sp', [
ift.RGSpace(4),
ift.PowerSpace(ift.RGSpace((4, 4), harmonic=True)),
ift.LMSpace(5),
ift.HPSpace(4),
ift.GLSpace(4)
])
@pytest.mark.parametrize('seed', [13, 2])
def test_value_inserter(sp, seed):
np.random.seed(seed)
ind = tuple([np.random.randint(0, ss - 1) for ss in sp.shape])
op = ift.ValueInserter(sp, ind)
f = ift.from_random('normal', ift.UnstructuredDomain((1,)))
inp = f.to_global_data()[0]
ret = op(f).to_global_data()
assert_(ret[ind] == inp)
assert_(np.sum(ret) == inp)
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