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

Merge branch 'master' into byebye_cython

parents d6a40b95 722948ed
Pipeline #12394 passed with stage
in 7 minutes and 22 seconds
......@@ -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)
......
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")
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")
......@@ -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
......@@ -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)
......
......@@ -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,