diff --git a/demos/wiener_filter.py b/demos/wiener_filter.py
index b9d432b46c40bf1d5779120ac8e1b2ffa9081022..53a73e10daf5c228afb2087a7a0df50248e640fb 100644
--- a/demos/wiener_filter.py
+++ b/demos/wiener_filter.py
@@ -11,16 +11,36 @@ rank = comm.rank
if __name__ == "__main__":
distribution_strategy = 'not'
-
+
+ #Setting up physical constants
+ #total length of Interval or Volume the field lives on, e.g. in meters
+ L = 2.
+ #typical distance over which the field is correlated (in same unit as L)
+ correlation_length = 0.1
+ #variance of field in position space sqrt(<|s_x|^2>) (in unit of s)
+ field_variance = 2.
+ #smoothing length that response (in same unit as L)
+ response_sigma = 0.1
+
+ #defining resolution (pixels per dimension)
+ N_pixels = 512
+
+ #Setting up derived constants
+ k_0 = 1./correlation_length
+ #note that field_variance**2 = a*k_0/4. for this analytic form of power
+ #spectrum
+ a = field_variance**2/k_0*4.
+ pow_spec = (lambda k: a / (1 + k/k_0) ** 4)
+ pixel_width = L/N_pixels
+
# Setting up the geometry
- s_space = RGSpace([512, 512])
+ s_space = RGSpace([N_pixels, N_pixels], distances = pixel_width)
fft = FFTOperator(s_space)
h_space = fft.target[0]
p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
# Creating the mock data
- pow_spec = (lambda k: 42 / (k + 1) ** 3)
S = create_power_operator(h_space, power_spectrum=pow_spec,
distribution_strategy=distribution_strategy)
@@ -30,7 +50,7 @@ if __name__ == "__main__":
sh = sp.power_synthesize(real_signal=True)
ss = fft.inverse_times(sh)
- R = SmoothingOperator(s_space, sigma=0.1)
+ R = SmoothingOperator(s_space, sigma=response_sigma)
signal_to_noise = 1
N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
diff --git a/demos/wiener_filter_harmonic.py b/demos/wiener_filter_harmonic.py
new file mode 100644
index 0000000000000000000000000000000000000000..692f0a6479d3b26b0b21ed67cddf924c1b04bb5f
--- /dev/null
+++ b/demos/wiener_filter_harmonic.py
@@ -0,0 +1,138 @@
+from nifty import *
+from mpi4py import MPI
+import plotly.offline as py
+import plotly.graph_objs as go
+
+
+comm = MPI.COMM_WORLD
+rank = comm.rank
+
+
+def plot_maps(x, name):
+
+ trace = [None]*len(x)
+
+ keys = x.keys()
+ field = x[keys[0]]
+ domain = field.domain[0]
+ shape = len(domain.shape)
+ max_n = domain.shape[0]*domain.distances[0]
+ step = domain.distances[0]
+ x_axis = np.arange(0, max_n, step)
+
+ if shape == 1:
+ for ii in xrange(len(x)):
+ trace[ii] = go.Scatter(x= x_axis, y=x[keys[ii]].val.get_full_data(), name=keys[ii])
+ fig = go.Figure(data=trace)
+
+ py.plot(fig, filename=name)
+
+ elif shape == 2:
+ for ii in xrange(len(x)):
+ py.plot([go.Heatmap(z=x[keys[ii]].val.get_full_data())], filename=keys[ii])
+ else:
+ raise TypeError("Only 1D and 2D field plots are supported")
+
+def plot_power(x, name):
+
+ layout = go.Layout(
+ xaxis=dict(
+ type='log',
+ autorange=True
+ ),
+ yaxis=dict(
+ type='log',
+ autorange=True
+ )
+ )
+
+ trace = [None]*len(x)
+
+ keys = x.keys()
+ field = x[keys[0]]
+ domain = field.domain[0]
+ x_axis = domain.kindex
+
+ for ii in xrange(len(x)):
+ trace[ii] = go.Scatter(x= x_axis, y=x[keys[ii]].val.get_full_data(), name=keys[ii])
+
+ fig = go.Figure(data=trace, layout=layout)
+ py.plot(fig, filename=name)
+
+np.random.seed(42)
+
+
+
+if __name__ == "__main__":
+
+ distribution_strategy = 'not'
+
+ # setting spaces
+ npix = np.array([500]) # number of pixels
+ total_volume = 1. # total length
+
+ # setting signal parameters
+ lambda_s = .05 # signal correlation length
+ sigma_s = 10. # signal variance
+
+
+ #setting response operator parameters
+ length_convolution = .025
+ exposure = 1.
+
+ # calculating parameters
+ k_0 = 4. / (2 * np.pi * lambda_s)
+ a_s = sigma_s ** 2. * lambda_s * total_volume
+
+ # creation of spaces
+ # x1 = RGSpace([npix,npix], distances=total_volume / npix,
+ # zerocenter=False)
+ # k1 = RGRGTransformation.get_codomain(x1)
+
+ x1 = HPSpace(64)
+ k1 = HPLMTransformation.get_codomain(x1)
+ p1 = PowerSpace(harmonic_partner=k1, logarithmic=False)
+
+ # creating Power Operator with given spectrum
+ spec = (lambda k: a_s / (1 + (k / k_0) ** 2) ** 2)
+ p_field = Field(p1, val=spec)
+ S_op = create_power_operator(k1, spec)
+
+ # creating FFT-Operator and Response-Operator with Gaussian convolution
+ # adjust dtype_target probperly
+ Fft_op = FFTOperator(domain=x1, target=k1,
+ domain_dtype=np.float64,
+ target_dtype=np.float64)
+ R_op = ResponseOperator(x1, sigma=[length_convolution],
+ exposure=[exposure])
+
+ # drawing a random field
+ sk = p_field.power_synthesize(real_power=True, mean=0.)
+ s = Fft_op.adjoint_times(sk)
+
+ signal_to_noise = 1
+ N_op = DiagonalOperator(R_op.target, diagonal=s.var()/signal_to_noise, bare=True)
+ n = Field.from_random(domain=R_op.target,
+ random_type='normal',
+ std=s.std()/np.sqrt(signal_to_noise),
+ mean=0)
+
+ d = R_op(s) + n
+
+ # Wiener filter
+ j = Fft_op.times(R_op.adjoint_times(N_op.inverse_times(d)))
+ D = HarmonicPropagatorOperator(S=S_op, N=N_op, R=R_op)
+
+ mk = D(j)
+
+ m = Fft_op.adjoint_times(mk)
+
+ # z={}
+ # z["signal"] = s
+ # z["reconstructed_map"] = m
+ # z["data"] = d
+ # z["lambda"] = R_op(s)
+ #
+ # plot_maps(z, "Wiener_filter.html")
+
+
diff --git a/demos/wiener_filter_unit.py b/demos/wiener_filter_unit.py
new file mode 100644
index 0000000000000000000000000000000000000000..2065626903e52a9e57dd4c2a6ef5613982a50552
--- /dev/null
+++ b/demos/wiener_filter_unit.py
@@ -0,0 +1,130 @@
+from nifty import *
+from mpi4py import MPI
+import plotly.offline as py
+import plotly.graph_objs as go
+
+
+comm = MPI.COMM_WORLD
+rank = comm.rank
+
+
+def plot_maps(x, name):
+
+ trace = [None]*len(x)
+
+ keys = x.keys()
+ field = x[keys[0]]
+ domain = field.domain[0]
+ shape = len(domain.shape)
+ max_n = domain.shape[0]*domain.distances[0]
+ step = domain.distances[0]
+ x_axis = np.arange(0, max_n, step)
+
+ if shape == 1:
+ for ii in xrange(len(x)):
+ trace[ii] = go.Scatter(x= x_axis, y=x[keys[ii]].val.get_full_data(), name=keys[ii])
+ fig = go.Figure(data=trace)
+
+ py.plot(fig, filename=name)
+
+ elif shape == 2:
+ for ii in xrange(len(x)):
+ py.plot([go.Heatmap(z=x[keys[ii]].val.get_full_data())], filename=keys[ii])
+ else:
+ raise TypeError("Only 1D and 2D field plots are supported")
+
+def plot_power(x, name):
+
+ layout = go.Layout(
+ xaxis=dict(
+ type='log',
+ autorange=True
+ ),
+ yaxis=dict(
+ type='log',
+ autorange=True
+ )
+ )
+
+ trace = [None]*len(x)
+
+ keys = x.keys()
+ field = x[keys[0]]
+ domain = field.domain[0]
+ x_axis = domain.kindex
+
+ for ii in xrange(len(x)):
+ trace[ii] = go.Scatter(x= x_axis, y=x[keys[ii]].val.get_full_data(), name=keys[ii])
+
+ fig = go.Figure(data=trace, layout=layout)
+ py.plot(fig, filename=name)
+
+np.random.seed(42)
+
+if __name__ == "__main__":
+
+ distribution_strategy = 'not'
+
+ # setting spaces
+ npix = np.array([500]) # number of pixels
+ total_volume = 1. # total length
+
+ # setting signal parameters
+ lambda_s = .05 # signal correlation length
+ sigma_s = 10. # signal variance
+
+
+ #setting response operator parameters
+ length_convolution = .025
+ exposure = 1.
+
+ # calculating parameters
+ k_0 = 4. / (2 * np.pi * lambda_s)
+ a_s = sigma_s ** 2. * lambda_s * total_volume
+
+ # creation of spaces
+ x1 = RGSpace(npix, distances=total_volume / npix,
+ zerocenter=False)
+ k1 = RGRGTransformation.get_codomain(x1)
+ p1 = PowerSpace(harmonic_domain=k1, log=False)
+
+ # creating Power Operator with given spectrum
+ spec = (lambda k: a_s / (1 + (k / k_0) ** 2) ** 2)
+ p_field = Field(p1, val=spec)
+ S_op = create_power_operator(k1, spec)
+
+ # creating FFT-Operator and Response-Operator with Gaussian convolution
+ Fft_op = FFTOperator(domain=x1, target=k1,
+ domain_dtype=np.float64,
+ target_dtype=np.complex128)
+ R_op = ResponseOperator(x1, sigma=[length_convolution],
+ exposure=[exposure])
+
+ # drawing a random field
+ sk = p_field.power_synthesize(real_signal=True, mean=0.)
+ s = Fft_op.inverse_times(sk)
+
+ signal_to_noise = 1
+ N_op = DiagonalOperator(R_op.target, diagonal=s.var()/signal_to_noise, bare=True)
+ n = Field.from_random(domain=R_op.target,
+ random_type='normal',
+ std=s.std()/np.sqrt(signal_to_noise),
+ mean=0)
+
+ d = R_op(s) + n
+
+ # Wiener filter
+ j = R_op.adjoint_times(N_op.inverse_times(d))
+ D = PropagatorOperator(S=S_op, N=N_op, R=R_op)
+
+ m = D(j)
+
+ z={}
+ z["signal"] = s
+ z["reconstructed_map"] = m
+ z["data"] = d
+ z["lambda"] = R_op(s)
+
+ plot_maps(z, "Wiener_filter.html")
+
+
diff --git a/nifty/operators/__init__.py b/nifty/operators/__init__.py
index 2621c7cf0294c81d6272e28c048b0a7c81ef53da..1a8940dc1dd33b1193e00fa0341add77202e6e38 100644
--- a/nifty/operators/__init__.py
+++ b/nifty/operators/__init__.py
@@ -34,4 +34,8 @@ from projection_operator import ProjectionOperator
from propagator_operator import PropagatorOperator
+from propagator_operator import HarmonicPropagatorOperator
+
from composed_operator import ComposedOperator
+
+from response_operator import ResponseOperator
diff --git a/nifty/operators/composed_operator/composed_operator.py b/nifty/operators/composed_operator/composed_operator.py
index 9557c3c3fcb93489c82c7827f6b421355dd9909a..b2e3ecd94a8a0d42e2e04399132e5ac848b35686 100644
--- a/nifty/operators/composed_operator/composed_operator.py
+++ b/nifty/operators/composed_operator/composed_operator.py
@@ -20,6 +20,59 @@ from nifty.operators.linear_operator import LinearOperator
class ComposedOperator(LinearOperator):
+ """ NIFTY class for composed operators.
+
+ The NIFTY composed operator class combines multiple linear operators.
+
+ Parameters
+ ----------
+ operators : tuple of NIFTy Operators
+ The tuple of LinearOperators.
+
+ Attributes
+ ----------
+ domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
+ The NIFTy.space in which the operator is defined.
+ target : tuple of DomainObjects, i.e. Spaces and FieldTypes
+ The NIFTy.space in which the outcome of the operator lives
+
+ Raises
+ ------
+ TypeError
+ Raised if
+ * an element of the operator list is not an instance of the
+ LinearOperator-baseclass.
+
+ Notes
+ -----
+ Very usefull in case one has to transform a Field living over a product
+ space (see example below).
+
+ Examples
+ --------
+ Minimal example of transforming a Field living on two domains into its
+ harmonic space.
+
+ >>> x1 = RGSpace(5)
+ >>> x2 = RGSpace(10)
+ >>> k1 = RGRGTransformation.get_codomain(x1)
+ >>> k2 = RGRGTransformation.get_codomain(x2)
+ >>> FFT1 = FFTOperator(domain=x1, target=k1,
+ domain_dtype=np.float64, target_dtype=np.complex128)
+ >>> FFT2 = FFTOperator(domain=x2, target=k2,
+ domain_dtype=np.float64, target_dtype=np.complex128)
+ >>> FFT = ComposedOperator((FFT1, FFT2)
+ >>> f = Field.from_random('normal', domain=(x1,x2))
+ >>> FFT.times(f)
+
+ See Also
+ --------
+ EndomorphicOperator, ProjectionOperator,
+ DiagonalOperator, SmoothingOperator, ResponseOperator,
+ PropagatorOperator, ComposedOperator
+
+ """
+
# ---Overwritten properties and methods---
def __init__(self, operators, default_spaces=None):
super(ComposedOperator, self).__init__(default_spaces)
diff --git a/nifty/operators/diagonal_operator/diagonal_operator.py b/nifty/operators/diagonal_operator/diagonal_operator.py
index 84dd3e48e6a494db81bee7a1908ae519b95a0aa4..dec42da42dc9a788b0d64912d9a5224474f346a5 100644
--- a/nifty/operators/diagonal_operator/diagonal_operator.py
+++ b/nifty/operators/diagonal_operator/diagonal_operator.py
@@ -28,6 +28,63 @@ from nifty.operators.endomorphic_operator import EndomorphicOperator
class DiagonalOperator(EndomorphicOperator):
+ """ NIFTY class for diagonal operators.
+
+ The NIFTY DiagonalOperator class is a subclass derived from the
+ EndomorphicOperator. It multiplies an input field pixel-wise with its
+ diagonal.
+
+
+ Parameters
+ ----------
+ domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
+ The domain on which the Operator's input Field lives.
+ diagonal : {scalar, list, array, Field, d2o-object}
+ The diagonal entries of the operator.
+ bare : boolean
+ Indicates whether the input for the diagonal is bare or not
+ (default: False).
+ copy : boolean
+ Internal copy of the diagonal (default: True)
+ distribution_strategy : string
+ setting the prober distribution_strategy of the
+ diagonal (default : None). In case diagonal is d2o-object or Field,
+ their distribution_strategy is used as a fallback.
+
+ Attributes
+ ----------
+ distribution_strategy : string
+ Defines the distribution_strategy of the distributed_data_object
+ in which the diagonal entries are stored in.
+
+ Raises
+ ------
+
+ Notes
+ -----
+ The ambiguity of bare or non-bare diagonal entries is based on the choice
+ of a matrix representation of the operator in question. The naive choice
+ of absorbing the volume weights into the matrix leads to a matrix-vector
+ calculus with the non-bare entries which seems intuitive, though.
+ The choice of keeping matrix entries and volume weights separate
+ deals with the bare entries that allow for correct interpretation
+ of the matrix entries; e.g., as variance in case of an covariance operator.
+
+ Examples
+ --------
+ >>> x_space = RGSpace(5)
+ >>> D = DiagonalOperator(x_space, diagonal=[1., 3., 2., 4., 6.])
+ >>> f = Field(x_space, val=2.)
+ >>> res = D.times(f)
+ >>> res.val
+ <distributed_data_object>
+ array([ 2., 6., 4., 8., 12.])
+
+ See Also
+ --------
+ EndomorphicOperator
+
+ """
# ---Overwritten properties and methods---
@@ -64,6 +121,21 @@ class DiagonalOperator(EndomorphicOperator):
operation=lambda z: z.adjoint().__rdiv__)
def diagonal(self, bare=False, copy=True):
+ """ Returns the diagonal of the Operator.
+
+ Parameters
+ ----------
+ bare : boolean
+ Whether the returned Field values should be bare or not.
+ copy : boolean
+ Whether the returned Field should be copied or not.
+
+ Returns
+ -------
+ out : Field
+ The diagonal of the Operator.
+
+ """
if bare:
diagonal = self._diagonal.weight(power=-1)
elif copy:
@@ -73,25 +145,100 @@ class DiagonalOperator(EndomorphicOperator):
return diagonal
def inverse_diagonal(self, bare=False):
+ """ Returns the inverse-diagonal of the operator.
+
+ Parameters
+ ----------
+ bare : boolean
+ Whether the returned Field values should be bare or not.
+
+ Returns
+ -------
+ out : Field
+ The inverse of the diagonal of the Operator.
+
+ """
return 1./self.diagonal(bare=bare, copy=False)
def trace(self, bare=False):
+ """ Returns the trace the operator.
+
+ Parameters
+ ----------
+ bare : boolean
+ Whether the returned Field values should be bare or not.
+
+ Returns
+ -------
+ out : scalar
+ The trace of the Operator.
+
+ """
return self.diagonal(bare=bare, copy=False).sum()
def inverse_trace(self, bare=False):
+ """ Returns the inverse-trace of the operator.
+
+ Parameters
+ ----------
+ bare : boolean
+ Whether the returned Field values should be bare or not.
+
+ Returns
+ -------
+ out : scalar
+ The inverse of the trace of the Operator.
+
+ """
return self.inverse_diagonal(bare=bare).sum()
def trace_log(self):
+ """ Returns the trave-log of the operator.
+
+ Returns
+ -------
+ out : scalar
+ the trace of the logarithm of the Operator.
+
+ """
log_diagonal = nifty_log(self.diagonal(copy=False))
return log_diagonal.sum()
def determinant(self):
+ """ Returns the determinant of the operator.
+
+ Returns
+ -------
+ out : scalar
+ out : scalar
+ the determinant of the Operator
+
+ """
+
return self.diagonal(copy=False).val.prod()
def inverse_determinant(self):
+ """ Returns the inverse-determinant of the operator.
+
+ Returns
+ -------
+ out : scalar
+ the inverse-determinant of the Operator
+
+ """
+
return 1/self.determinant()
def log_determinant(self):
+ """ Returns the log-eterminant of the operator.
+
+ Returns
+ -------
+ out : scalar
+ the log-determinant of the Operator
+
+ """
+
return np.log(self.determinant())
# ---Mandatory properties and methods---
@@ -117,6 +264,17 @@ class DiagonalOperator(EndomorphicOperator):
@property
def distribution_strategy(self):
+ """
+ distribution_strategy : string
+ Defines the way how the diagonal operator is distributed
+ among the nodes. Available distribution_strategies are:
+ 'fftw', 'equal' and 'not'.
+
+ Notes :
+ https://arxiv.org/abs/1606.05385
+
+ """
+
return self._distribution_strategy
def _parse_distribution_strategy(self, distribution_strategy, val):
@@ -134,6 +292,20 @@ class DiagonalOperator(EndomorphicOperator):
return distribution_strategy
def set_diagonal(self, diagonal, bare=False, copy=True):
+ """ Sets the diagonal of the Operator.
+
+ Parameters
+ ----------
+ diagonal : {scalar, list, array, Field, d2o-object}
+ The diagonal entries of the operator.
+ bare : boolean
+ Indicates whether the input for the diagonal is bare or not
+ (default: False).
+ copy : boolean
+ Specifies if a copy of the input shall be made (default: True).
+
+ """
+
# use the casting functionality from Field to process `diagonal`
f = Field(domain=self.domain,
val=diagonal,
diff --git a/nifty/operators/endomorphic_operator/endomorphic_operator.py b/nifty/operators/endomorphic_operator/endomorphic_operator.py
index 0d9cc49cafc55636fdcb3776c572a3130e84a045..1340708be765601e6dfcffa274f523292b330e39 100644
--- a/nifty/operators/endomorphic_operator/endomorphic_operator.py
+++ b/nifty/operators/endomorphic_operator/endomorphic_operator.py
@@ -22,6 +22,45 @@ from nifty.operators.linear_operator import LinearOperator
class EndomorphicOperator(LinearOperator):
+ """ NIFTY class for endomorphic operators.
+
+ The NIFTY EndomorphicOperator class is a class derived from the
+ LinearOperator. By definition, domain and target are the same in
+ EndomorphicOperator.
+
+ Parameters
+ ----------
+ #TODO: Copy Parameters from LinearOperator
+
+ Attributes
+ ----------
+ domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
+ The domain on which the Operator's input Field lives.
+ target : tuple of DomainObjects, i.e. Spaces and FieldTypes
+ The domain in which the outcome of the operator lives. As the Operator
+ is endomorphic this is the same as its domain.
+ self_adjoint : boolean
+ Indicates whether the operator is self_adjoint or not.
+
+ Raises
+ ------
+ NotImplementedError
+ Raised if
+ * self_adjoint is not defined
+
+ Notes
+ -----
+
+ Examples
+ --------
+
+
+ See Also
+ --------
+ DiagonalOperator, SmoothingOperator,
+ PropagatorOperator, ProjectionOperator
+
+ """
# ---Overwritten properties and methods---
@@ -67,4 +106,7 @@ class EndomorphicOperator(LinearOperator):
@abc.abstractproperty
def self_adjoint(self):
+ """ States whether the Operator is self_adjoint or not.
+ """
+
raise NotImplementedError
diff --git a/nifty/operators/invertible_operator_mixin/invertible_operator_mixin.py b/nifty/operators/invertible_operator_mixin/invertible_operator_mixin.py
index e4e50e8d6a28696211d4206e01a6c9a0143dfcad..4786a4f9375c7e56e6b52d38005263f2b6f69057 100644
--- a/nifty/operators/invertible_operator_mixin/invertible_operator_mixin.py
+++ b/nifty/operators/invertible_operator_mixin/invertible_operator_mixin.py
@@ -22,6 +22,43 @@ from nifty.field import Field
class InvertibleOperatorMixin(object):
+ """ Mixin class to invert implicit defined operators.
+
+ To invert the application of a given implicitly defined operator on a
+ field, this class gives the necessary functionality. Inheriting
+ functionality from this class provides the derived class with the inverse
+ to the given implicitly definied application of the operator on a field.
+ (e.g. .inverse_times vs. .times and
+ .adjoint_times vs. .adjoint_inverse_times)
+
+ Parameters
+ ----------
+ inverter : Inverter
+ An instance of an Inverter class.
+ (default: ConjugateGradient)
+
+ preconditioner : LinearOperator
+ Preconditions the minimizaion problem
+
+ Attributes
+ ----------
+
+ Raises
+ ------
+
+ Notes
+ -----
+
+ Examples
+ --------
+ The PropagatorOperator inherits from InvertibleOperatorMixin.
+
+ See Also
+ --------
+ PropagatorOperator
+
+ """
+
def __init__(self, inverter=None, preconditioner=None, *args, **kwargs):
self.__preconditioner = preconditioner
if inverter is not None:
diff --git a/nifty/operators/linear_operator/linear_operator.py b/nifty/operators/linear_operator/linear_operator.py
index 24d231a48494e5f6b7f7a3fca451170b83f72060..bf0160ba964562ca2d80fceb4da95f8491642cd9 100644
--- a/nifty/operators/linear_operator/linear_operator.py
+++ b/nifty/operators/linear_operator/linear_operator.py
@@ -25,6 +25,51 @@ import nifty.nifty_utilities as utilities
class LinearOperator(Loggable, object):
+ """NIFTY base class for linear operators.
+
+ The base NIFTY operator class is an abstract class from which
+ other specific operator subclasses, including those preimplemented
+ in NIFTY (e.g. the EndomorphicOperator, ProjectionOperator,
+ DiagonalOperator, SmoothingOperator, ResponseOperator,
+ PropagatorOperator, ComposedOperator) are derived.
+
+ Parameters
+ ----------
+ default_spaces : tuple of ints *optional*
+ Defines on which space(s) of a given field the Operator acts by
+ default (default: None)
+
+ Attributes
+ ----------
+ domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
+ The domain on which the Operator's input Field lives.
+ target : tuple of DomainObjects, i.e. Spaces and FieldTypes
+ The domain in which the Operators result lives.
+ unitary : boolean
+ Indicates whether the Operator is unitary or not.
+
+ Raises
+ ------
+ NotImplementedError
+ Raised if
+ * domain is not defined
+ * target is not defined
+ * unitary is not set to (True/False)
+
+ Notes
+ -----
+ All Operators wihtin NIFTy are linear and must therefore be a subclasses of
+ the LinearOperator. A LinearOperator must have the attributes domain,
+ target and unitary to be properly defined.
+
+ See Also
+ --------
+ EndomorphicOperator, ProjectionOperator,
+ DiagonalOperator, SmoothingOperator, ResponseOperator,
+ PropagatorOperator, ComposedOperator
+
+ """
+
__metaclass__ = NiftyMeta
def __init__(self, default_spaces=None):
@@ -35,14 +80,38 @@ class LinearOperator(Loggable, object):
@abc.abstractproperty
def domain(self):
+ """
+ domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
+ The domain on which the Operator's input Field lives.
+ Every Operator which inherits from the abstract LinearOperator
+ base class must have this attribute.
+
+ """
+
raise NotImplementedError
@abc.abstractproperty
def target(self):
+ """
+ target : tuple of DomainObjects, i.e. Spaces and FieldTypes
+ The domain on which the Operator's output Field lives.
+ Every Operator which inherits from the abstract LinearOperator
+ base class must have this attribute.
+
+ """
+
raise NotImplementedError
@abc.abstractproperty
def unitary(self):
+ """
+ unitary : boolean
+ States whether the Operator is unitary or not.
+ Every Operator which inherits from the abstract LinearOperator
+ base class must have this attribute.
+
+ """
+
raise NotImplementedError
@property
@@ -56,55 +125,132 @@ class LinearOperator(Loggable, object):
def __call__(self, *args, **kwargs):
return self.times(*args, **kwargs)
- def times(self, x, spaces=None, **kwargs):
- spaces = self._check_input_compatibility(x, spaces)
+ def times(self, x, spaces=None):
+ """ Applies the Operator to a given Field.
+
+ Operator and Field have to live over the same domain.
- y = self._times(x, spaces, **kwargs)
+ Parameters
+ ----------
+ x : Field
+ The input Field.
+ spaces : tuple of ints
+ Defines on which space(s) of the given Field the Operator acts.
+
+ Returns
+ -------
+ out : Field
+ The processed Field living on the target-domain.
+
+ """
+
+ spaces = self._check_input_compatibility(x, spaces)
+ y = self._times(x, spaces)
return y
- def inverse_times(self, x, spaces=None, **kwargs):
+ def inverse_times(self, x, spaces=None):
+ """ Applies the inverse-Operator to a given Field.
+
+ Operator and Field have to live over the same domain.
+
+ Parameters
+ ----------
+ x : Field
+ The input Field.
+ spaces : tuple of ints
+ Defines on which space(s) of the given Field the Operator acts.
+
+ Returns
+ -------
+ out : Field
+ The processed Field living on the target-domain.
+
+ """
+
spaces = self._check_input_compatibility(x, spaces, inverse=True)
try:
- y = self._inverse_times(x, spaces, **kwargs)
+ y = self._inverse_times(x, spaces)
except(NotImplementedError):
if (self.unitary):
- y = self._adjoint_times(x, spaces, **kwargs)
+ y = self._adjoint_times(x, spaces)
else:
raise
return y
- def adjoint_times(self, x, spaces=None, **kwargs):
+ def adjoint_times(self, x, spaces=None):
+ """ Applies the adjoint-Operator to a given Field.
+
+ Operator and Field have to live over the same domain.
+
+ Parameters
+ ----------
+ x : Field
+ applies the Operator to the given Field
+ spaces : tuple of ints
+ defines on which space of the given Field the Operator acts
+
+ Returns
+ -------
+ out : Field
+ The processed Field living on the target-domain.
+
+ """
+
+ if self.unitary:
+ return self.inverse_times(x, spaces)
+
spaces = self._check_input_compatibility(x, spaces, inverse=True)
try:
- y = self._adjoint_times(x, spaces, **kwargs)
+ y = self._adjoint_times(x, spaces)
except(NotImplementedError):
if (self.unitary):
- y = self._inverse_times(x, spaces, **kwargs)
+ y = self._inverse_times(x, spaces)
else:
raise
return y
- # If the operator supports inverse() then the inverse adjoint is identical
- # to the adjoint inverse. We provide both names for convenience.
- def adjoint_inverse_times(self, x, spaces=None, **kwargs):
+ def adjoint_inverse_times(self, x, spaces=None):
+ """ Applies the adjoint-inverse Operator to a given Field.
+
+ Operator and Field have to live over the same domain.
+
+ Parameters
+ ----------
+ x : Field
+ applies the Operator to the given Field
+ spaces : tuple of ints
+ defines on which space of the given Field the Operator acts
+
+ Returns
+ -------
+ out : Field
+ The processed Field living on the target-domain.
+
+ Notes
+ -----
+ If the operator has an `inverse` then the inverse adjoint is identical
+ to the adjoint inverse. We provide both names for convenience.
+
+ See Also
+ --------
+
+ """
+
spaces = self._check_input_compatibility(x, spaces)
try:
- y = self._adjoint_inverse_times(x, spaces, **kwargs)
+ y = self._adjoint_inverse_times(x, spaces)
except(NotImplementedError):
- try:
- y = self._inverse_adjoint_times(x, spaces, **kwargs)
- except(NotImplementedError):
- if self.unitary:
- y = self._times(x, spaces, **kwargs)
- else:
- raise
+ if self.unitary:
+ y = self._times(x, spaces)
+ else:
+ raise
return y
- def inverse_adjoint_times(self, x, spaces=None, **kwargs):
- return adjoint_inverse_times(x, spaces, **kwargs)
+ def inverse_adjoint_times(self, x, spaces=None):
+ return self.adjoint_inverse_times(x, spaces)
def _times(self, x, spaces):
raise NotImplementedError(
@@ -129,7 +275,7 @@ class LinearOperator(Loggable, object):
def _check_input_compatibility(self, x, spaces, inverse=False):
if not isinstance(x, Field):
raise ValueError(
- "supplied object is not a `nifty.Field`.")
+ "supplied object is not a `Field`.")
if spaces is None:
spaces = self.default_spaces
diff --git a/nifty/operators/projection_operator/projection_operator.py b/nifty/operators/projection_operator/projection_operator.py
index 2217d6aa248b3b2a846a91ef8a38531ffca6849e..9a51f363634dc7bf0c1112b3ae60ee0b74b785ce 100644
--- a/nifty/operators/projection_operator/projection_operator.py
+++ b/nifty/operators/projection_operator/projection_operator.py
@@ -24,6 +24,45 @@ from nifty.operators.endomorphic_operator import EndomorphicOperator
class ProjectionOperator(EndomorphicOperator):
+ """ NIFTY class for projection operators.
+
+ The NIFTY ProjectionOperator class is a class derived from the
+ EndomorphicOperator.
+
+ Parameters
+ ----------
+ projection_field : Field
+ Field on which the operator projects
+
+ Attributes
+ ----------
+
+ Raises
+ ------
+ TypeError
+ Raised if
+ * if projection_field is not a Field
+
+ Notes
+ -----
+
+
+ Examples
+ --------
+ >>> x_space = RGSpace(5)
+ >>> f1 = Field(x_space, val=3.)
+ >>> f2 = Field(x_space, val=5.)
+ >>> P = ProjectionOperator(f1)
+ >>> res = P.times(f2)
+ >>> res.val
+ <distributed_data_object>
+ array([ 225., 225., 225., 225., 225.])
+
+ See Also
+ --------
+ EndomorphicOperator
+
+ """
# ---Overwritten properties and methods---
diff --git a/nifty/operators/propagator_operator/__init__.py b/nifty/operators/propagator_operator/__init__.py
index 48672e0631900024d2cb9ea2aeb155cd5bc7d3e0..61d962c4299059adbe8f5fab63d8f7dd8a27239b 100644
--- a/nifty/operators/propagator_operator/__init__.py
+++ b/nifty/operators/propagator_operator/__init__.py
@@ -17,3 +17,4 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
from propagator_operator import PropagatorOperator
+from harmonic_propagator_operator import HarmonicPropagatorOperator
\ No newline at end of file
diff --git a/nifty/operators/propagator_operator/harmonic_propagator_operator.py b/nifty/operators/propagator_operator/harmonic_propagator_operator.py
new file mode 100644
index 0000000000000000000000000000000000000000..3ad6501fe217cd4610ef7e1b8787aeb7f3eb4b0e
--- /dev/null
+++ b/nifty/operators/propagator_operator/harmonic_propagator_operator.py
@@ -0,0 +1,144 @@
+# NIFTy
+# Copyright (C) 2017 Theo Steininger
+#
+# Author: Theo Steininger
+#
+# 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/>.
+
+from nifty.operators import EndomorphicOperator,\
+ FFTOperator,\
+ InvertibleOperatorMixin
+
+
+class HarmonicPropagatorOperator(InvertibleOperatorMixin, EndomorphicOperator):
+ """ NIFTY Harmonic Propagator Operator D.
+
+ The propagator operator D, is known from the Wiener Filter.
+ Its inverse functional form might look like:
+ D = (S^(-1) + M)^(-1)
+ D = (S^(-1) + N^(-1))^(-1)
+ D = (S^(-1) + R^(\dagger) N^(-1) R)^(-1)
+ In contrast to the PropagatorOperator the inference is done in the
+ harmonic space.
+
+ Parameters
+ ----------
+ S : LinearOperator
+ Covariance of the signal prior.
+ M : LinearOperator
+ Likelihood contribution.
+ R : LinearOperator
+ Response operator translating signal to (noiseless) data.
+ N : LinearOperator
+ Covariance of the noise prior or the likelihood, respectively.
+ inverter : class to invert explicitly defined operators
+ (default:ConjugateGradient)
+ preconditioner : Field
+ numerical preconditioner to speed up convergence
+
+ Attributes
+ ----------
+
+ Raises
+ ------
+ ValueError
+ is raised if
+ * neither N nor M is given
+
+ Notes
+ -----
+
+ Examples
+ --------
+
+ See Also
+ --------
+ Scientific reference
+ https://arxiv.org/abs/0806.3474
+
+ """
+
+ # ---Overwritten properties and methods---
+
+ def __init__(self, S=None, M=None, R=None, N=None, inverter=None,
+ preconditioner=None):
+ """
+ Sets the standard operator properties and `codomain`, `_A1`, `_A2`,
+ and `RN` if required.
+
+ Parameters
+ ----------
+ S : operator
+ Covariance of the signal prior.
+ M : operator
+ Likelihood contribution.
+ R : operator
+ Response operator translating signal to (noiseless) data.
+ N : operator
+ Covariance of the noise prior or the likelihood, respectively.
+
+ """
+ # infer domain, and target
+ # infer domain, and target
+ if M is not None:
+ self._codomain = M.domain
+ self._likelihood = M.times
+
+ elif N is None:
+ raise ValueError("Either M or N must be given!")
+
+ elif R is not None:
+ self._codomain = R.domain
+ self._likelihood = \
+ lambda z: R.adjoint_times(N.inverse_times(R.times(z)))
+ else:
+ self._codomain = N.domain
+ self._likelihood = lambda z: N.inverse_times(z)
+
+ self._domain = S.domain
+ self._S = S
+ self._fft_S = FFTOperator(self._domain, target=self._codomain)
+
+ super(HarmonicPropagatorOperator, self).__init__(inverter=inverter,
+ preconditioner=preconditioner)
+
+ # ---Mandatory properties and methods---
+
+ @property
+ def domain(self):
+ return self._domain
+
+ @property
+ def self_adjoint(self):
+ return True
+
+ @property
+ def unitary(self):
+ return False
+
+ # ---Added properties and methods---
+ def _likelihood_times(self, x, spaces=None):
+ transformed_x = self._fft_S.times(x, spaces=spaces)
+ y = self._likelihood(transformed_x)
+ transformed_y = self._fft_S.adjoint_times(y, spaces=spaces)
+ result = x.copy_empty()
+ result.set_val(transformed_y, copy=False)
+ return result
+
+ def _inverse_times(self, x, spaces):
+ pre_result = self._S.inverse_times(x, spaces)
+ pre_result += self._likelihood_times(x)
+ result = x.copy_empty()
+ result.set_val(pre_result, copy=False)
+ return result
diff --git a/nifty/operators/propagator_operator/propagator_operator.py b/nifty/operators/propagator_operator/propagator_operator.py
index 506cce30b171102c6c954c321f87d3c310a9ff60..cb70e7b7cb9fbaa179517b8fe61abe0ac44ab1af 100644
--- a/nifty/operators/propagator_operator/propagator_operator.py
+++ b/nifty/operators/propagator_operator/propagator_operator.py
@@ -22,6 +22,59 @@ from nifty.operators import EndomorphicOperator,\
class PropagatorOperator(InvertibleOperatorMixin, EndomorphicOperator):
+ """ NIFTY Propagator Operator D.
+
+ The propagator operator D, is known from the Wiener Filter.
+ Its inverse functional form might look like:
+ D = (S^(-1) + M)^(-1)
+ D = (S^(-1) + N^(-1))^(-1)
+ D = (S^(-1) + R^(\dagger) N^(-1) R)^(-1)
+
+ Parameters
+ ----------
+ S : LinearOperator
+ Covariance of the signal prior.
+ M : LinearOperator
+ Likelihood contribution.
+ R : LinearOperator
+ Response operator translating signal to (noiseless) data.
+ N : LinearOperator
+ Covariance of the noise prior or the likelihood, respectively.
+ inverter : class to invert explicitly defined operators
+ (default:ConjugateGradient)
+ preconditioner : Field
+ numerical preconditioner to speed up convergence
+
+ Attributes
+ ----------
+
+ Raises
+ ------
+ ValueError
+ is raised if
+ * neither N nor M is given
+
+ Notes
+ -----
+
+ Examples
+ --------
+ >>> x_space = RGSpace(4)
+ >>> k_space = RGRGTransformation.get_codomain(x_space)
+ >>> f = Field(x_space, val=[2, 4, 6, 8])
+ >>> S = create_power_operator(k_space, spec=1)
+ >>> N = DiagonalOperaor(f.domain, diag=1)
+ >>> D = PropagatorOperator(S=S, N=N) # D^{-1} = S^{-1} + N^{-1}
+ >>> D(f).val
+ <distributed_data_object>
+ array([ 1., 2., 3., 4.]
+
+ See Also
+ --------
+ Scientific reference
+ https://arxiv.org/abs/0806.3474
+
+ """
# ---Overwritten properties and methods---
diff --git a/nifty/operators/response_operator/__init__.py b/nifty/operators/response_operator/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..2ae4e4775341b0c8dde48b9b0ffde30da8e14abe
--- /dev/null
+++ b/nifty/operators/response_operator/__init__.py
@@ -0,0 +1 @@
+from response_operator import ResponseOperator
diff --git a/nifty/operators/response_operator/response_operator.py b/nifty/operators/response_operator/response_operator.py
new file mode 100644
index 0000000000000000000000000000000000000000..f5ac8b05ab578ea34cf60918b2f9b675ecaa813b
--- /dev/null
+++ b/nifty/operators/response_operator/response_operator.py
@@ -0,0 +1,121 @@
+import numpy as np
+from nifty import Field,\
+ FieldArray
+from nifty.operators.linear_operator import LinearOperator
+from nifty.operators.smoothing_operator import SmoothingOperator
+from nifty.operators.composed_operator import ComposedOperator
+from nifty.operators.diagonal_operator import DiagonalOperator
+
+class ResponseOperator(LinearOperator):
+ """ NIFTy ResponseOperator (example)
+
+ This NIFTy ResponseOperator provides the user with an example how a
+ ResponseOperator can look like. It smoothes and exposes a field. The
+ outcome of the Operator is geometrically not ordered as typical data
+ set are.
+
+ Parameters
+ ----------
+ domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
+ The domains on which the operator lives. Either one space or a list
+ of spaces
+ sigma : list(np.float)
+ Defines the smoothing length of the operator for each space it lives on
+ exposure : list(np.float)
+ Defines the exposure of the operator for each space it lives on
+
+
+ Attributes
+ ----------
+
+ Raises
+ ------
+ ValueError:
+ raised if:
+ * len of sigma-list and exposure-list are not equal
+
+ Notes
+ -----
+
+ Examples
+ --------
+ >>> x1 = RGSpace(5)
+ >>> x2 = RGSpace(10)
+ >>> R = ResponseOperator(domain=(x1,x2), sigma=[.5, .25],
+ exposure=[2.,3.])
+ >>> f = Field((x1,x2), val=4.)
+ >>> R.times(f)
+ <distributed_data_object>
+ array([[ 24., 24., 24., 24., 24., 24., 24., 24., 24., 24.],
+ [ 24., 24., 24., 24., 24., 24., 24., 24., 24., 24.],
+ [ 24., 24., 24., 24., 24., 24., 24., 24., 24., 24.],
+ [ 24., 24., 24., 24., 24., 24., 24., 24., 24., 24.],
+ [ 24., 24., 24., 24., 24., 24., 24., 24., 24., 24.]])
+
+ See Also
+ --------
+
+ """
+
+ def __init__(self, domain,
+ sigma=[1.], exposure=[1.],
+ default_spaces=None):
+ super(ResponseOperator, self).__init__(default_spaces)
+
+ self._domain = self._parse_domain(domain)
+
+ shapes = len(self._domain)*[None]
+ shape_target = []
+ for ii in xrange(len(shapes)):
+ shapes[ii] = self._domain[ii].shape
+ shape_target = np.append(shape_target, self._domain[ii].shape)
+
+ self._target = self._parse_domain(FieldArray(shape_target))
+ self._sigma = sigma
+ self._exposure = exposure
+
+ self._kernel = len(self._domain)*[None]
+
+ for ii in xrange(len(self._kernel)):
+ self._kernel[ii] = SmoothingOperator(self._domain[ii],
+ sigma=self._sigma[ii])
+
+ self._composed_kernel = ComposedOperator(self._kernel)
+
+ self._exposure_op = len(self._domain)*[None]
+ if len(self._exposure_op) != len(self._kernel):
+ raise ValueError("Definition of kernel and exposure do not suit "
+ "each other")
+ else:
+ for ii in xrange(len(self._exposure_op)):
+ self._exposure_op[ii] = DiagonalOperator(
+ self._domain[ii],
+ diagonal=self._exposure[ii])
+ self._composed_exposure = ComposedOperator(self._exposure_op)
+
+ @property
+ def domain(self):
+ return self._domain
+
+ @property
+ def target(self):
+ return self._target
+
+ @property
+ def unitary(self):
+ return False
+
+ def _times(self, x, spaces):
+ res = self._composed_kernel.times(x, spaces)
+ res = self._composed_exposure.times(res, spaces)
+ # res = res.weight(power=1)
+ # removing geometric information
+ return Field(self._target, val=res.val)
+
+ def _adjoint_times(self, x, spaces):
+ # setting correct spaces
+ res = Field(self.domain, val=x.val)
+ res = self._composed_exposure.adjoint_times(res, spaces)
+ res = res.weight(power=-1)
+ res = self._composed_kernel.adjoint_times(res, spaces)
+ return res
diff --git a/nifty/operators/smoothing_operator/smoothing_operator.py b/nifty/operators/smoothing_operator/smoothing_operator.py
index fe27a16e3e14fff42142bca530c5ae63f553a5b3..abd6ea3f0bbb05a5a1b713d194a4d2c7bb049e79 100644
--- a/nifty/operators/smoothing_operator/smoothing_operator.py
+++ b/nifty/operators/smoothing_operator/smoothing_operator.py
@@ -25,6 +25,56 @@ from d2o import STRATEGIES
class SmoothingOperator(EndomorphicOperator):
+ """ NIFTY class for smoothing operators.
+
+ The NIFTy SmoothingOperator smooths Fields, with a given kernel length.
+ Fields which are not living over a PowerSpace are smoothed
+ via a gaussian convolution. Fields living over the PowerSpace are directly
+ smoothed.
+
+ Parameters
+ ----------
+ domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
+ The Space on which the operator acts
+ sigma : float
+ Sets the length of the Gaussian convolution kernel
+ log_distances : boolean
+ States whether the convolution happens on the logarithmic grid or not.
+
+ Attributes
+ ----------
+ sigma : float
+ Sets the length of the Gaussian convolution kernel
+ log_distances : boolean
+ States whether the convolution happens on the logarithmic grid or not.
+
+ Raises
+ ------
+ ValueError
+ Raised if
+ * the given domain inherits more than one space. The
+ SmoothingOperator acts only on one Space.
+
+ Notes
+ -----
+
+ Examples
+ --------
+ >>> x = RGSpace(5)
+ >>> S = SmoothingOperator(x, sigma=1.)
+ >>> f = Field(x, val=[1,2,3,4,5])
+ >>> S.times(f).val
+ <distributed_data_object>
+ array([ 3., 3., 3., 3., 3.])
+
+ See Also
+ --------
+ DiagonalOperator, SmoothingOperator,
+ PropagatorOperator, ProjectionOperator,
+ ComposedOperator
+
+ """
+
# ---Overwritten properties and methods---
def __init__(self, domain, sigma, log_distances=False,
default_spaces=None):
@@ -260,7 +310,7 @@ class SmoothingOperator(EndomorphicOperator):
transformed_x.val.set_local_data(local_transformed_x, copy=False)
- smoothed_x = Transformator.inverse_times(transformed_x, spaces=spaces)
+ smoothed_x = Transformator.adjoint_times(transformed_x, spaces=spaces)
result = x.copy_empty()
result.set_val(smoothed_x, copy=False)
diff --git a/nifty/spaces/rg_space/rg_space.py b/nifty/spaces/rg_space/rg_space.py
index cbce79be5e88828b4e4fd3421051b957615eb9dc..cb70c66452be8bc762b53ea91801ac077ee91304 100644
--- a/nifty/spaces/rg_space/rg_space.py
+++ b/nifty/spaces/rg_space/rg_space.py
@@ -53,19 +53,19 @@ class RGSpace(Space):
shape : {int, numpy.ndarray}
Number of grid points or numbers of gridpoints along each axis.
zerocenter : {bool, numpy.ndarray}, *optional*
- Whether x==0 (or k==0, respectively) is located in the center of
- the grid (or the center of each axis speparately) or not.
- (default: False).
+ Whether x==0 (or k==0, respectively) is located in the center of
+ the grid (or the center of each axis speparately) or not.
+ (default: False).
distances : {float, numpy.ndarray}, *optional*
Distance between two grid points along each axis
(default: None).
If distances==None:
if harmonic==True, all distances will be set to 1
if harmonic==False, the distance along each axis will be
- set to the inverse of the number of points along that
- axis.
+ set to the inverse of the number of points along that
+ axis.
harmonic : bool, *optional*
- Whether the space represents a grid in position or harmonic space.
+ Whether the space represents a grid in position or harmonic space.
(default: False).
Attributes
@@ -92,15 +92,6 @@ class RGSpace(Space):
def __init__(self, shape, zerocenter=False, distances=None,
harmonic=False):
- """
- Sets the attributes for an RGSpace class instance.
-
-
-
- Returns
- -------
- None
- """
self._harmonic = bool(harmonic)
super(RGSpace, self).__init__()
@@ -227,7 +218,8 @@ class RGSpace(Space):
Returns
-------
distributed_data_object
- A d2o containing the distances
+ A d2o containing the distances.
+
"""
shape = self.shape
diff --git a/nifty/sugar.py b/nifty/sugar.py
index 2eba6e7a2fa7f36e6643e61fc56a3c9cb240cfbb..5bf7dbc6f8b23f2753bebe1cde110462549dfcbb 100644
--- a/nifty/sugar.py
+++ b/nifty/sugar.py
@@ -37,8 +37,8 @@ def create_power_operator(domain, power_spectrum, dtype=None,
val=power_spectrum, dtype=dtype,
distribution_strategy=distribution_strategy)
- f = fp.power_synthesize(mean=1, std=0, real_signal=True)
+ f = fp.power_synthesize(mean=1, std=0, real_signal=False)
power_operator = DiagonalOperator(domain, diagonal=f, bare=True)
- return power_operator
+ return power_operator
\ No newline at end of file
diff --git a/test/test_init.py b/test/test_init.py
new file mode 100644
index 0000000000000000000000000000000000000000..7d23ef5c98a197262e6c75ecf90aebb078991ff8
--- /dev/null
+++ b/test/test_init.py
@@ -0,0 +1,4 @@
+from nifty import *
+
+#This tests if it is possible to import all of nifties methods. Experience shows this is not always possible.
+pass
\ No newline at end of file
diff --git a/test/test_operators/test_diagonal_operator.py b/test/test_operators/test_diagonal_operator.py
new file mode 100644
index 0000000000000000000000000000000000000000..90f1ef7c4ac2f095d400d4d60e51314b8324cc08
--- /dev/null
+++ b/test/test_operators/test_diagonal_operator.py
@@ -0,0 +1,134 @@
+import unittest
+
+import numpy as np
+from numpy.testing import assert_equal,\
+ assert_allclose,\
+ assert_approx_equal
+
+from nifty import Field,\
+ DiagonalOperator
+
+from test.common import generate_spaces
+
+from itertools import product
+from test.common import expand
+
+class DiagonalOperator_Tests(unittest.TestCase):
+ spaces = generate_spaces()
+
+ @expand(product(spaces, [True, False], [True, False]))
+ def test_property(self, space, bare, copy):
+ diag = Field.from_random('normal', domain=space)
+ D = DiagonalOperator(space, diagonal=diag)
+ if D.domain[0] != space:
+ raise TypeError
+ if D.unitary != False:
+ raise TypeError
+ if D.self_adjoint != True:
+ raise TypeError
+
+ @expand(product(spaces, [True, False], [True, False]))
+ def test_times_adjoint(self, space, bare, copy):
+ rand1 = Field.from_random('normal', domain=space)
+ rand2 = Field.from_random('normal', domain=space)
+ diag = Field.from_random('normal', domain=space)
+ D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
+ tt1 = rand1.dot(D.times(rand2))
+ tt2 = rand2.dot(D.times(rand1))
+ assert_approx_equal(tt1, tt2)
+
+ @expand(product(spaces, [True, False], [True, False]))
+ def test_times_inverse(self, space, bare, copy):
+ rand1 = Field.from_random('normal', domain=space)
+ diag = Field.from_random('normal', domain=space)
+ D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
+ tt1 = D.times(D.inverse_times(rand1))
+ assert_allclose(rand1.val.get_full_data(), tt1.val.get_full_data())
+
+ @expand(product(spaces, [True, False], [True, False]))
+ def test_times(self, space, bare, copy):
+ rand1 = Field.from_random('normal', domain=space)
+ diag = Field.from_random('normal', domain=space)
+ D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
+ tt = D.times(rand1)
+ assert_equal(tt.domain[0], space)
+
+ @expand(product(spaces, [True, False], [True, False]))
+ def test_adjoint_times(self, space, bare, copy):
+ rand1 = Field.from_random('normal', domain=space)
+ diag = Field.from_random('normal', domain=space)
+ D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
+ tt = D.adjoint_times(rand1)
+ assert_equal(tt.domain[0], space)
+
+ @expand(product(spaces, [True, False], [True, False]))
+ def test_inverse_times(self, space, bare, copy):
+ rand1 = Field.from_random('normal', domain=space)
+ diag = Field.from_random('normal', domain=space)
+ D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
+ tt = D.inverse_times(rand1)
+ assert_equal(tt.domain[0], space)
+
+ @expand(product(spaces, [True, False], [True, False]))
+ def test_adjoint_inverse_times(self, space, bare, copy):
+ rand1 = Field.from_random('normal', domain=space)
+ diag = Field.from_random('normal', domain=space)
+ D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
+ tt = D.adjoint_inverse_times(rand1)
+ assert_equal(tt.domain[0], space)
+
+ @expand(product(spaces, [True, False]))
+ def test_diagonal(self, space, copy):
+ diag = Field.from_random('normal', domain=space)
+ D = DiagonalOperator(space, diagonal=diag, copy=copy)
+ diag_op = D.diagonal()
+ assert_allclose(diag.val.get_full_data(), diag_op.val.get_full_data())
+
+ @expand(product(spaces, [True, False]))
+ def test_inverse(self, space, copy):
+ diag = Field.from_random('normal', domain=space)
+ D = DiagonalOperator(space, diagonal=diag, copy=copy)
+ diag_op = D.inverse_diagonal()
+ assert_allclose(1./diag.val.get_full_data(), diag_op.val.get_full_data())
+
+ @expand(product(spaces, [True, False]))
+ def test_trace(self, space, copy):
+ diag = Field.from_random('normal', domain=space)
+ D = DiagonalOperator(space, diagonal=diag, copy=copy)
+ trace_op = D.trace()
+ assert_allclose(trace_op, np.sum(diag.val.get_full_data()))
+
+ @expand(product(spaces, [True, False]))
+ def test_inverse_trace(self, space, copy):
+ diag = Field.from_random('normal', domain=space)
+ D = DiagonalOperator(space, diagonal=diag, copy=copy)
+ trace_op = D.inverse_trace()
+ assert_allclose(trace_op, np.sum(1./diag.val.get_full_data()))
+
+ @expand(product(spaces, [True, False]))
+ def test_trace_log(self, space, copy):
+ diag = Field.from_random('normal', domain=space)
+ D = DiagonalOperator(space, diagonal=diag, copy=copy)
+ trace_log = D.trace_log()
+ assert_allclose(trace_log, np.sum(np.log(diag.val.get_full_data())))
+
+ @expand(product(spaces, [True, False]))
+ def test_determinant(self, space, copy):
+ diag = Field.from_random('normal', domain=space)
+ D = DiagonalOperator(space, diagonal=diag, copy=copy)
+ det = D.determinant()
+ assert_allclose(det, np.prod(diag.val.get_full_data()))
+
+ @expand(product(spaces, [True, False], [True, False]))
+ def test_inverse_determinant(self, space, bare, copy):
+ diag = Field.from_random('normal', domain=space)
+ D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
+ inv_det = D.inverse_determinant()
+ assert_allclose(inv_det, 1./D.determinant())
+
+ @expand(product(spaces, [True, False], [True, False]))
+ def test_log_determinant(self, space, bare, copy):
+ diag = Field.from_random('normal', domain=space)
+ D = DiagonalOperator(space, diagonal=diag, bare=bare, copy=copy)
+ log_det = D.log_determinant()
+ assert_allclose(log_det, np.log(D.determinant()))
\ No newline at end of file