Commit ad3efa7a authored by Martin Reinecke's avatar Martin Reinecke
Browse files

merge NIFTy_5

parents 52d9f783 d1e5a734
...@@ -19,6 +19,7 @@ from .field import Field ...@@ -19,6 +19,7 @@ from .field import Field
from .multi_field import MultiField from .multi_field import MultiField
from .operators.operator import Operator from .operators.operator import Operator
from .operators.adder import Adder
from .operators.diagonal_operator import DiagonalOperator from .operators.diagonal_operator import DiagonalOperator
from .operators.distributors import DOFDistributor, PowerDistributor from .operators.distributors import DOFDistributor, PowerDistributor
from .operators.domain_tuple_field_inserter import DomainTupleFieldInserter from .operators.domain_tuple_field_inserter import DomainTupleFieldInserter
...@@ -33,7 +34,6 @@ from .operators.field_zero_padder import FieldZeroPadder ...@@ -33,7 +34,6 @@ from .operators.field_zero_padder import FieldZeroPadder
from .operators.inversion_enabler import InversionEnabler from .operators.inversion_enabler import InversionEnabler
from .operators.linear_operator import LinearOperator from .operators.linear_operator import LinearOperator
from .operators.mask_operator import MaskOperator from .operators.mask_operator import MaskOperator
from .operators.offset_operator import OffsetOperator
from .operators.qht_operator import QHTOperator from .operators.qht_operator import QHTOperator
from .operators.regridding_operator import RegriddingOperator from .operators.regridding_operator import RegriddingOperator
from .operators.sampling_enabler import SamplingEnabler from .operators.sampling_enabler import SamplingEnabler
......
...@@ -16,7 +16,6 @@ ...@@ -16,7 +16,6 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from ..minimization.energy_adapter import EnergyAdapter from ..minimization.energy_adapter import EnergyAdapter
from ..multi_domain import MultiDomain
from ..multi_field import MultiField from ..multi_field import MultiField
from ..operators.distributors import PowerDistributor from ..operators.distributors import PowerDistributor
from ..operators.energy_operators import StandardHamiltonian, InverseGammaLikelihood from ..operators.energy_operators import StandardHamiltonian, InverseGammaLikelihood
...@@ -35,25 +34,27 @@ def make_adjust_variances(a, ...@@ -35,25 +34,27 @@ def make_adjust_variances(a,
Constructs a Hamiltonian to solve constant likelihood optimizations of the Constructs a Hamiltonian to solve constant likelihood optimizations of the
form phi = a * xi under the constraint that phi remains constant. form phi = a * xi under the constraint that phi remains constant.
FIXME xi is white.
Parameters Parameters
---------- ----------
a : Operator a : Operator
Operator which gives the amplitude when evaluated at a position Gives the amplitude when evaluated at a position.
xi : Operator xi : Operator
Operator which gives the excitation when evaluated at a position Gives the excitation when evaluated at a position.
position : Field, MultiField position : Field, MultiField
Position of the whole problem Position of the entire problem.
samples : Field, MultiField samples : Field, MultiField
Residual samples of the whole problem Residual samples of the whole problem.
scaling : Float scaling : Float
Optional rescaling of the Likelihood Optional rescaling of the Likelihood.
ic_samp : Controller ic_samp : Controller
Iteration Controller for Hamiltonian Iteration Controller for Hamiltonian.
Returns Returns
------- -------
StandardHamiltonian Energy
A Hamiltonian that can be used for further minimization A Hamiltonian that can be used for further minimization.
""" """
d = a*xi d = a*xi
...@@ -79,6 +80,9 @@ def do_adjust_variances(position, ...@@ -79,6 +80,9 @@ def do_adjust_variances(position,
minimizer, minimizer,
xi_key='xi', xi_key='xi',
samples=[]): samples=[]):
'''
FIXME
'''
h_space = position[xi_key].domain[0] h_space = position[xi_key].domain[0]
pd = PowerDistributor(h_space, amplitude_operator.target[0]) pd = PowerDistributor(h_space, amplitude_operator.target[0])
......
...@@ -146,7 +146,7 @@ def dynamic_operator(*, ...@@ -146,7 +146,7 @@ def dynamic_operator(*,
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 When evaluated, this operator returns the Green's function representation
in harmonic space. This result can be used as a convolution kernel to in harmonic space. This result can be used as a convolution kernel to
construct solutions of the homogeneous stochastic differential equation construct solutions of the homogeneous stochastic differential equation
...@@ -216,7 +216,7 @@ def dynamic_lightcone_operator(*, ...@@ -216,7 +216,7 @@ def dynamic_lightcone_operator(*,
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. 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 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 space-time. Thereby the first axis of the space is set to be the teporal
......
...@@ -20,8 +20,8 @@ import numpy as np ...@@ -20,8 +20,8 @@ import numpy as np
from ..domain_tuple import DomainTuple from ..domain_tuple import DomainTuple
from ..domains.power_space import PowerSpace from ..domains.power_space import PowerSpace
from ..field import Field from ..field import Field
from ..operators.adder import Adder
from ..operators.exp_transform import ExpTransform from ..operators.exp_transform import ExpTransform
from ..operators.offset_operator import OffsetOperator
from ..operators.qht_operator import QHTOperator from ..operators.qht_operator import QHTOperator
from ..operators.slope_operator import SlopeOperator from ..operators.slope_operator import SlopeOperator
from ..operators.symmetrizing_operator import SymmetrizingOperator from ..operators.symmetrizing_operator import SymmetrizingOperator
...@@ -29,7 +29,7 @@ from ..sugar import makeOp ...@@ -29,7 +29,7 @@ from ..sugar import makeOp
def _ceps_kernel(k, a, k0): def _ceps_kernel(k, a, k0):
return (a/(1+np.sum((k.T/k0)**2, axis=-1).T))**2 return (a/(1 + np.sum((k.T/k0)**2, axis=-1).T))**2
def CepstrumOperator(target, a, k0): def CepstrumOperator(target, a, k0):
...@@ -189,7 +189,7 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']): ...@@ -189,7 +189,7 @@ def SLAmplitude(*, target, n_pix, a, k0, sm, sv, im, iv, keys=['tau', 'phi']):
sig = np.array([sv, iv]) sig = np.array([sv, iv])
mean = Field.from_global_data(sl.domain, mean) mean = Field.from_global_data(sl.domain, mean)
sig = Field.from_global_data(sl.domain, sig) sig = Field.from_global_data(sl.domain, sig)
linear = (sl @ OffsetOperator(mean) @ makeOp(sig)).ducktape(keys[1]) linear = sl @ Adder(mean) @ makeOp(sig).ducktape(keys[1])
# Combine linear and smooth component # Combine linear and smooth component
loglog_ampl = 0.5*(smooth + linear) loglog_ampl = 0.5*(smooth + linear)
......
...@@ -15,18 +15,22 @@ ...@@ -15,18 +15,22 @@
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from ..field import Field
from ..multi_field import MultiField
from .operator import Operator from .operator import Operator
class OffsetOperator(Operator): class Adder(Operator):
"""Shifts the input by a fixed field. """Adds a fixed field.
Parameters Parameters
---------- ----------
field : Field field : Field or MultiField
The field by which the input is shifted. The field by which the input is shifted.
""" """
def __init__(self, field): def __init__(self, field):
if not isinstance(field, (Field, MultiField)):
raise TypeError
self._field = field self._field = field
self._domain = self._target = field.domain self._domain = self._target = field.domain
......
...@@ -27,8 +27,8 @@ class DiagonalOperator(EndomorphicOperator): ...@@ -27,8 +27,8 @@ class DiagonalOperator(EndomorphicOperator):
"""Represents a :class:`LinearOperator` which is diagonal. """Represents a :class:`LinearOperator` which is diagonal.
The NIFTy DiagonalOperator class is a subclass derived from the The NIFTy DiagonalOperator class is a subclass derived from the
:class:`EndomorphicOperator`. It multiplies an input field pixel-wise with its :class:`EndomorphicOperator`. It multiplies an input field pixel-wise with
diagonal. its diagonal.
Parameters Parameters
---------- ----------
......
...@@ -33,6 +33,9 @@ class ScalingOperator(EndomorphicOperator): ...@@ -33,6 +33,9 @@ class ScalingOperator(EndomorphicOperator):
Notes Notes
----- -----
:class:`Operator` supports the multiplication with a scalar. So one does
not need instantiate :class:`ScalingOperator` explicitly in most cases.
Formally, this operator always supports all operation modes (times, Formally, this operator always supports all operation modes (times,
adjoint_times, inverse_times and inverse_adjoint_times), even if `factor` adjoint_times, inverse_times and inverse_adjoint_times), even if `factor`
is 0 or infinity. It is the user's responsibility to apply the operator is 0 or infinity. It is the user's responsibility to apply the operator
......
...@@ -29,8 +29,8 @@ class SlopeOperator(LinearOperator): ...@@ -29,8 +29,8 @@ class SlopeOperator(LinearOperator):
Slope and y-intercept of this line are the two parameters which are Slope and y-intercept of this line are the two parameters which are
defined on an UnstructeredDomain (in this order) which is the domain of defined on an UnstructeredDomain (in this order) which is the domain of
the operator. Being a LogRGSpace instance each pixel has a well-defined coordinate the operator. Being a LogRGSpace instance each pixel has a well-defined
value. coordinate value.
The y-intercept is defined to be the value at t_0 of the target. The y-intercept is defined to be the value at t_0 of the target.
......
...@@ -15,9 +15,6 @@ ...@@ -15,9 +15,6 @@
# #
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik. # NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik.
from functools import reduce
from operator import mul
import numpy as np import numpy as np
from ..domain_tuple import DomainTuple from ..domain_tuple import DomainTuple
...@@ -28,7 +25,7 @@ from .linear_operator import LinearOperator ...@@ -28,7 +25,7 @@ from .linear_operator import LinearOperator
class ValueInserter(LinearOperator): class ValueInserter(LinearOperator):
"""Inserts one value into a field which is constant otherwise. """Inserts one value into a field which is zero otherwise.
Parameters Parameters
---------- ----------
...@@ -36,16 +33,11 @@ class ValueInserter(LinearOperator): ...@@ -36,16 +33,11 @@ class ValueInserter(LinearOperator):
index : iterable of int index : iterable of int
The index of the target into which the value of the domain shall be The index of the target into which the value of the domain shall be
inserted. inserted.
default_value : float
Constant value which is inserted everywhere where the input operator
is not inserted. Default is 0.
""" """
def __init__(self, target, index, default_value=0.): def __init__(self, target, index):
self._domain = makeDomain(UnstructuredDomain(1)) self._domain = makeDomain(UnstructuredDomain(1))
self._target = DomainTuple.make(target) self._target = DomainTuple.make(target)
# Type and value checks
index = tuple(index) index = tuple(index)
if not all([ if not all([
isinstance(n, int) and n >= 0 and n < self.target.shape[i] isinstance(n, int) and n >= 0 and n < self.target.shape[i]
...@@ -54,19 +46,17 @@ class ValueInserter(LinearOperator): ...@@ -54,19 +46,17 @@ class ValueInserter(LinearOperator):
raise TypeError raise TypeError
if not len(index) == len(self.target.shape): if not len(index) == len(self.target.shape):
raise ValueError raise ValueError
np.empty(self.target.shape)[index]
self._index = index self._index = index
self._dv = float(default_value)
self._dvsum = self._dv*(reduce(mul, self.target.shape) - 1)
self._capability = self.TIMES | self.ADJOINT_TIMES 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): def apply(self, x, mode):
self._check_input(x, mode) self._check_input(x, mode)
x = x.to_global_data() x = x.to_global_data()
if mode == self.TIMES: if mode == self.TIMES:
res = np.full(self.target.shape, self._dv, dtype=x.dtype) res = np.zeros(self.target.shape, dtype=x.dtype)
res[self._index] = x res[self._index] = x
else: else:
res = np.full((1,), x[self._index] + self._dvsum, dtype=x.dtype) res = np.full((1,), x[self._index], dtype=x.dtype)
return Field.from_global_data(self._tgt(mode), res) return Field.from_global_data(self._tgt(mode), res)
...@@ -17,7 +17,7 @@ ...@@ -17,7 +17,7 @@
import numpy as np import numpy as np
import pytest import pytest
from numpy.testing import assert_allclose from numpy.testing import assert_
import nifty5 as ift import nifty5 as ift
...@@ -37,17 +37,5 @@ def test_value_inserter(sp, seed): ...@@ -37,17 +37,5 @@ def test_value_inserter(sp, seed):
f = ift.from_random('normal', ift.UnstructuredDomain((1,))) f = ift.from_random('normal', ift.UnstructuredDomain((1,)))
inp = f.to_global_data()[0] inp = f.to_global_data()[0]
ret = op(f).to_global_data() ret = op(f).to_global_data()
assert_allclose(ret[ind], inp) assert_(ret[ind] == inp)
assert_allclose(np.sum(ret), inp) assert_(np.sum(ret) == inp)
def test_value_inserter_nonzero():
sp = ift.RGSpace(4)
ind = (1,)
default = 1.24
op = ift.ValueInserter(sp, ind, default)
f = ift.from_random('normal', ift.UnstructuredDomain((1,)))
inp = f.to_global_data()[0]
ret = op(f).to_global_data()
assert_allclose(ret[ind], inp)
assert_allclose(np.sum(ret), inp + 3*default)
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