Commit 63af9bd5 authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'master' into Tests_Operators

# Conflicts:
#	nifty/operators/diagonal_operator/diagonal_operator.py
#	nifty/operators/linear_operator/linear_operator.py
parents 336775e0 6861afd1
Pipeline #12386 passed with stage
in 6 minutes and 4 seconds
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().real)], 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(32)
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
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.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)
# z["j"] = j
#
# plot_maps(z, "Wiener_filter.html")
......@@ -85,9 +85,12 @@ if __name__ == "__main__":
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 = 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_domain=k1, log=False)
# creating Power Operator with given spectrum
......@@ -96,14 +99,15 @@ if __name__ == "__main__":
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.complex128)
target_dtype=np.float64)
R_op = ResponseOperator(x1, sigma=[length_convolution],
exposure=[exposure])
# drawing a random field
sk = p_field.power_synthesize(real_signal=True, mean=0.)
sk = p_field.power_synthesize(decompose_power=True, mean=0.)
s = Fft_op.adjoint_times(sk)
signal_to_noise = 1
......@@ -123,12 +127,12 @@ if __name__ == "__main__":
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")
# z={}
# z["signal"] = s
# z["reconstructed_map"] = m
# z["data"] = d
# z["lambda"] = R_op(s)
#
# plot_maps(z, "Wiener_filter.html")
......@@ -50,8 +50,8 @@ from spaces import *
from operators import *
from plotting import *
from probing import *
from sugar import *
import plotting
......@@ -24,6 +24,21 @@ from keepers import Loggable,\
class DomainObject(Versionable, Loggable, object):
"""The abstract class that can be used as a domain for a field.
This holds all the information and functionality a field needs to know
about its domain and how the data of the field are stored.
Attributes
----------
dim : int
Number of pixel-dimensions of the underlying data object.
shape : tuple
Shape of the array that stores the degrees of freedom for any field
on this domain.
"""
__metaclass__ = NiftyMeta
def __init__(self):
......
......@@ -20,45 +20,47 @@ from nifty.operators.linear_operator import LinearOperator
class ComposedOperator(LinearOperator):
""" NIFTY class for composed operators.
"""NIFTY class for composed operators.
The NIFTY composed operator class inherits multiple Operators of various kinds acting on
a Field living over a product space.
The NIFTY composed operator class combines multiple linear operators.
Parameters
----------
Operators : tuple(NIFTy.LinearOperators)
Contains the list of LinearOperators.
operators : tuple of NIFTy Operators
The tuple of LinearOperators.
Attributes
----------
domain : NIFTy.space
domain : tuple of DomainObjects, i.e. Spaces and FieldTypes
The NIFTy.space in which the operator is defined.
target : NIFTy.space
target : tuple of DomainObjects, i.e. Spaces and FieldTypes
The NIFTy.space in which the outcome of the operator lives
Raises
------
TypeError
Raised if
* the elements of the operator list is not an instance of the
LinearOperator-baseclass
* an element of the operator list is not an instance of the
LinearOperator-baseclass.
Notes
-----
Very usefull in case one has to transform a NIFTy.Field living over the product space. (see example below)
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.
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)
>>> 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)
......
......@@ -28,32 +28,34 @@ 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.
"""NIFTY class for diagonal operators.
The NIFTY DiagonalOperator class is a subclass derived from the
EndomorphicOperator.
Parameters
----------
domain : NIFTy.Space
The Space on which the operator acts
diagonal : {scalar, list, array, NIFTy.Field, d2o-object}
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 diagonal entries are bare or not
(default: False)
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 reused.
their distribution_strategy is used as a fallback.
Attributes
----------
distribution_strategy : string
Defines the diagonal is distributed among the nodes.
Defines the distribution_strategy of the distributed_data_object
in which the diagonal entries are stored in.
Raises
------
......@@ -71,12 +73,12 @@ class DiagonalOperator(EndomorphicOperator):
Examples
--------
>>> x_space = RGSpace(5)
>>> D = DiagonalOperator(x_space, diagonal=2.)
>>> f = Field(x_space, val=1.)
>>> 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., 2., 2., 2., 2.])
array([ 2., 6., 4., 8., 12.])
See Also
--------
......@@ -119,7 +121,7 @@ class DiagonalOperator(EndomorphicOperator):
operation=lambda z: z.adjoint().__rdiv__)
def diagonal(self, bare=False, copy=True):
""" Returns the diagonal of the operator.
""" Returns the diagonal of the Operator.
Parameters
----------
......@@ -130,11 +132,9 @@ class DiagonalOperator(EndomorphicOperator):
Returns
-------
out : NIFTy.Field
the diagonal of the Operator
out : Field
The diagonal of the Operator.
See Also
--------
"""
if bare:
diagonal = self._diagonal.weight(power=-1)
......@@ -154,13 +154,10 @@ class DiagonalOperator(EndomorphicOperator):
Returns
-------
out : NIFTy.Field
the inverse-diagonal of the Operator
out : Field
The inverse of the diagonal of the Operator.
See Also
--------
"""
return 1/self.diagonal(bare=bare, copy=False)
"""
return 1./self.diagonal(bare=bare, copy=False)
def trace(self, bare=False):
......@@ -174,10 +171,8 @@ class DiagonalOperator(EndomorphicOperator):
Returns
-------
out : scalar
the trace of the Operator
The trace of the Operator.
See Also
--------
"""
return self.diagonal(bare=bare, copy=False).sum()
......@@ -192,25 +187,19 @@ class DiagonalOperator(EndomorphicOperator):
Returns
-------
out : scalar
the inverse-trace of the Operator
The inverse of the trace of the Operator.
See Also
--------
"""
return self.inverse_diagonal(bare=bare).sum()
def trace_log(self):
""" Returns the trave-log of the operator.
Parameters
----------
Returns
-------
out : scalar
the trace-log of the Operator
the trace of the logarithm of the Operator.
See Also
--------
"""
log_diagonal = nifty_log(self.diagonal(copy=False))
return log_diagonal.sum()
......@@ -218,47 +207,38 @@ class DiagonalOperator(EndomorphicOperator):
def determinant(self):
""" Returns the determinant of the operator.
Parameters
----------
Returns
-------
out : scalar
out : scalar
the determinant of the Operator
See Also
--------
"""
return self.diagonal(copy=False).val.prod()
def inverse_determinant(self):
""" Returns the inverse-determinant of the operator.
Parameters
----------
Returns
-------
out : scalar
the inverse-determinant of the Operator
See Also
--------
"""
return 1/self.determinant()
def log_determinant(self):
""" Returns the log-eterminant of the operator.
Parameters
----------
Returns
-------
out : scalar
the log-determinant of the Operator
See Also
--------
"""
return np.log(self.determinant())
# ---Mandatory properties and methods---
......@@ -287,18 +267,14 @@ class DiagonalOperator(EndomorphicOperator):
"""
distribution_strategy : string
Defines the way how the diagonal operator is distributed
among the nodes.
Popoular ones are:
'fftw'
'all'
'local'
'slicing'
'not'
'hdf5'
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):
......@@ -320,22 +296,16 @@ class DiagonalOperator(EndomorphicOperator):
Parameters
----------
diagonal : {scalar, list, array, NIFTy.Field, d2o-object}
diagonal : {scalar, list, array, Field, d2o-object}
The diagonal entries of the operator.
bare : boolean
Indicates whether the diagonal entries are bare or not
(default: False)
Indicates whether the input for the diagonal is bare or not
(default: False).
copy : boolean
Internal copy of the diagonal (default: True)
Specifies if a copy of the input shall be made (default: True).
Returns
-------
See Also
--------
"""
# use the casting functionality from Field to process `diagonal`
f = Field(domain=self.domain,
val=diagonal,
......
......@@ -22,22 +22,23 @@ from nifty.operators.linear_operator import LinearOperator
class EndomorphicOperator(LinearOperator):
""" NIFTY class for endomorphic operators.
"""NIFTY class for endomorphic operators.
The NIFTY EndomorphicOperator class is a class derived from the
LinearOperator. Domain and target are the same in any EndomorphicOperator.
Prominent other specific operator subclasses, in NIFTy are
(e.g. DiagonalOperator, SmoothingOperator,
PropagatorOperator, ProjectionOperator).
LinearOperator. By definition, domain and target are the same in
EndomorphicOperator.
Parameters
----------
#TODO: Copy Parameters from LinearOperator
Attributes
----------
target : NIFTy.space
The NIFTy.space in which the outcome of the operator lives.
As the Operator is endomorphic this is the same as its domain.
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.
......@@ -61,33 +62,9 @@ class EndomorphicOperator(LinearOperator):
"""
# ---Overwritten properties and methods---
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 : NIFTY.Field
the input Field on which the operator acts on
spaces : integer (default: None)
defines on which space of the given Field the Operator acts
**kwargs
Additional keyword arguments get passed to the used copy_empty
routine.
Returns
-------
out : NIFTy.Field
the processed Field living on the domain space
See Also
--------
"""
if self.self_adjoint and self.unitary:
return self.times(x, spaces)
else:
......@@ -96,29 +73,6 @@ class EndomorphicOperator(LinearOperator):
spaces=spaces)
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 : NIFTY.Field
applies the Operator to the given Field
spaces : integer (default: None)
defines on which space of the given Field the Operator acts
**kwargs
Additional keyword arguments get passed to the used copy_empty
routine.
Returns
-------
out : NIFTy.Field
the processed Field living on the domain space