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

merge from master

parents 0c50175d 7c50e35f
import numpy as np
from nifty import RGSpace, PowerSpace, Field, FFTOperator, ComposedOperator,\
SmoothingOperator, DiagonalOperator, create_power_operator
from nifty.library import WienerFilterCurvature
from nifty import *
#import plotly.offline as pl
#import plotly.graph_objs as go
......@@ -10,36 +14,37 @@ 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
distribution_strategy = 'fftw'
# 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)
# 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)
# variance of field in position space sqrt(<|s_x|^2>) (in unit of s)
field_variance = 2.
#smoothing length of response (in same unit as L)
# smoothing length of response (in same unit as L)
response_sigma = 0.1
#defining resolution (pixels per dimension)
# defining resolution (pixels per dimension)
N_pixels = 512
#Setting up derived constants
# 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
# 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
pixel_length = L/N_pixels
# Setting up the geometry
s_space = RGSpace([N_pixels, N_pixels], distances = pixel_width)
fft = FFTOperator(s_space)
s_space = RGSpace([N_pixels, N_pixels], distances=pixel_length)
fft = FFTOperator(s_space, domain_dtype=np.float, target_dtype=np.complex)
h_space = fft.target[0]
inverse_fft = FFTOperator(h_space, target=s_space,
domain_dtype=np.complex, target_dtype=np.float)
p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
# Creating the mock data
S = create_power_operator(h_space, power_spectrum=pow_spec,
......@@ -51,6 +56,7 @@ if __name__ == "__main__":
ss = fft.inverse_times(sh)
R = SmoothingOperator(s_space, sigma=response_sigma)
R_harmonic = ComposedOperator([inverse_fft, R], default_spaces=[0, 0])
signal_to_noise = 1
N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
......@@ -63,7 +69,9 @@ if __name__ == "__main__":
# Wiener filter
j = R.adjoint_times(N.inverse_times(d))
D = PropagatorOperator(S=S, N=N, R=R)
j = R_harmonic.adjoint_times(N.inverse_times(d))
wiener_curvature = WienerFilterCurvature(S=S, N=N, R=R_harmonic)
m = wiener_curvature.inverse_times(j)
m_s = inverse_fft(m)
m = D(j)
......@@ -66,7 +66,7 @@ variable_default_field_dtype = keepers.Variable(
variable_default_distribution_strategy = keepers.Variable(
'default_distribution_strategy',
['not', 'fftw', 'equal'],
lambda z: (('pyfftw' in dependency_injector)
lambda z: (('fftw' in dependency_injector)
if z == 'fftw' else True),
genus='str')
......
......@@ -596,22 +596,15 @@ class Field(Loggable, Versionable, object):
@staticmethod
def _hermitian_decomposition(domain, val, spaces, domain_axes,
preserve_gaussian_variance=False):
# hermitianize for the first space
(h, a) = domain[spaces[0]].hermitian_decomposition(
val,
domain_axes[spaces[0]])
# hermitianize all remaining spaces using the iterative formula
for space in spaces[1:]:
(hh, ha) = domain[space].hermitian_decomposition(
h,
domain_axes[space])
(ah, aa) = domain[space].hermitian_decomposition(
a,
domain_axes[space])
c = (hh - ha - ah + aa).conjugate()
full = (hh + ha + ah + aa)
h = (full + c)/2.
a = (full - c)/2.
flipped_val = val
for space in spaces:
flipped_val = domain[space].hermitianize_inverter(
x=flipped_val,
axes=domain_axes[space])
flipped_val = flipped_val.conjugate()
h = (val + flipped_val)/2.
a = val - h
# correct variance
if preserve_gaussian_variance:
......
......@@ -117,7 +117,6 @@ class FFTOperator(LinearOperator):
super(FFTOperator, self).__init__(default_spaces)
# Initialize domain and target
self._domain = self._parse_domain(domain)
if len(self.domain) != 1:
raise ValueError("TransformationOperator accepts only exactly one "
......
......@@ -8,12 +8,21 @@ from .smoothing_operator import SmoothingOperator
class FFTSmoothingOperator(SmoothingOperator):
def _smooth(self, x, spaces, inverse):
Transformator = FFTOperator(x.domain[spaces[0]])
def __init__(self, domain, sigma, log_distances=False,
default_spaces=None):
super(FFTSmoothingOperator, self).__init__(
domain=domain,
sigma=sigma,
log_distances=log_distances,
default_spaces=default_spaces)
self._transformator_cache = {}
def _smooth(self, x, spaces, inverse):
# transform to the (global-)default codomain and perform all remaining
# steps therein
transformed_x = Transformator(x, spaces=spaces)
transformator = self._get_transformator(x.dtype)
transformed_x = transformator(x, spaces=spaces)
codomain = transformed_x.domain[spaces[0]]
coaxes = transformed_x.domain_axes[spaces[0]]
......@@ -51,9 +60,18 @@ class FFTSmoothingOperator(SmoothingOperator):
transformed_x.val.set_local_data(local_transformed_x, copy=False)
smoothed_x = Transformator.adjoint_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)
return result
def _get_transformator(self, dtype):
if dtype not in self._transformator_cache:
self._transformator_cache[dtype] = FFTOperator(
self.domain,
domain_dtype=dtype,
target_dtype=np.complex)
return self._transformator_cache[dtype]
......@@ -89,22 +89,6 @@ class LMSpace(Space):
super(LMSpace, self).__init__()
self._lmax = self._parse_lmax(lmax)
def hermitian_decomposition(self, x, axes=None):
if issubclass(x.dtype.type, np.complexfloating):
hermitian_part = x.copy_empty()
anti_hermitian_part = x.copy_empty()
hermitian_part[:] = x.real
anti_hermitian_part[:] = x.imag * 1j
else:
hermitian_part = x.copy()
anti_hermitian_part = x.copy_empty()
anti_hermitian_part[:] = 0
return (hermitian_part, anti_hermitian_part)
def hermitian_fixed_points(self):
return None
# ---Mandatory properties and methods---
def __repr__(self):
......
......@@ -100,26 +100,6 @@ class RGSpace(Space):
self._distances = self._parse_distances(distances)
self._zerocenter = self._parse_zerocenter(zerocenter)
def hermitian_decomposition(self, x, axes=None,
preserve_gaussian_variance=False):
# check axes
if axes is None:
axes = range(len(self.shape))
assert len(x.shape) >= len(self.shape), "shapes mismatch"
assert len(axes) == len(self.shape), "axes mismatch"
# compute the hermitian part
flipped_x = self._hermitianize_inverter(x, axes=axes)
flipped_x = flipped_x.conjugate()
# average x and flipped_x.
hermitian_part = x + flipped_x
hermitian_part /= 2.
# use subtraction since it is faster than flipping another time
anti_hermitian_part = (x-hermitian_part)
return (hermitian_part, anti_hermitian_part)
def hermitian_fixed_points(self):
dimensions = len(self.shape)
mid_index = np.array(self.shape)//2
......@@ -138,7 +118,7 @@ class RGSpace(Space):
fixed_points += [tuple(index * mid_index)]
return fixed_points
def _hermitianize_inverter(self, x, axes):
def hermitianize_inverter(self, x, axes):
# calculate the number of dimensions the input array has
dimensions = len(x.shape)
# prepare the slicing object which will be used for mirroring
......
......@@ -161,35 +161,34 @@ class Space(DomainObject):
raise NotImplementedError(
"There is no generic co-smoothing kernel for Space base class.")
def hermitian_decomposition(self, x, axes,
preserve_gaussian_variance=False):
""" Decomposes x into its hermitian and anti-hermitian constituents.
def hermitian_fixed_points(self):
""" Returns the array points which remain invariant under the action
of `hermitianize_inverter`
This method decomposes a field's array x into its hermitian and
antihermitian part, which corresponds to real and imaginary part
in a corresponding harmonic partner space. This is an internal function
that is mainly used for power-synthesizing and -analyzing Fields.
Returns
-------
list of index-tuples
The list contains the index-coordinates of the invariant points.
"""
return None
def hermitianize_inverter(self, x, axes):
""" Inverts/flips x in the context of Hermitian decomposition.
This method is mainly used for power-synthesizing and -analyzing
Fields.
Parameters
----------
x : distributed_data_object
The field's val object.
axes : tuple of ints
Specifies the axes of x which correspond to this space.
preserve_gaussian_variance : bool *optional*
If the hermitian decomposition is done via computing the half
sums and differences of `x` and mirrored `x`, all points except the
fixed points lose half of their variance. If `x` is complex also
they lose half of their variance since the real(/imaginary) part
gets lost.
Returns
-------
(distributed_data_object, distributed_data_object)
A tuple of two distributed_data_objects, the first being the
hermitian and the second the anti-hermitian part of x.
distributed_data_object
The Hermitian-flipped of x.
"""
raise NotImplementedError
return x
......@@ -20,7 +20,7 @@ import unittest
import numpy as np
from numpy.testing import assert_, assert_equal, assert_raises,\
assert_almost_equal
assert_almost_equal, assert_array_almost_equal
from d2o import distributed_data_object
from nifty import LMSpace
from test.common import expand
......@@ -108,15 +108,12 @@ class LMSpaceFunctionalityTests(unittest.TestCase):
for key, value in expected.iteritems():
assert_equal(getattr(l, key), value)
@expand(get_hermitian_configs())
def test_hermitian_decomposition(self, x, real, imag):
def test_hermitianize_inverter(self):
l = LMSpace(5)
assert_almost_equal(
l.hermitian_decomposition(distributed_data_object(x))[0],
real)
assert_almost_equal(
l.hermitian_decomposition(distributed_data_object(x))[1],
imag*1j)
v = distributed_data_object(global_shape=l.shape, dtype=np.complex128)
v[:] = np.random.random(l.shape) + 1j*np.random.random(l.shape)
inverted = l.hermitianize_inverter(v, axes=(0,))
assert_array_almost_equal(inverted.get_full_data(), v.get_full_data())
@expand(get_weight_configs())
def test_weight(self, x, power, axes, inplace, expected):
......
......@@ -129,7 +129,7 @@ class PowerSpaceInterfaceTest(unittest.TestCase):
class PowerSpaceConsistencyCheck(unittest.TestCase):
@expand(CONSISTENCY_CONFIGS)
def test_rhopindexConsistency(self, harmonic_partner, distribution_strategy,
binbounds, nbin,logarithmic):
binbounds, nbin, logarithmic):
if distribution_strategy == "fftw":
if not hasattr(gdi.get('fftw'), 'FFTW_MPI'):
raise SkipTest
......
......@@ -21,6 +21,8 @@ from __future__ import division
import unittest
import numpy as np
from d2o import distributed_data_object
from numpy.testing import assert_, assert_equal, assert_almost_equal
from nifty import RGSpace
from test.common import expand
......@@ -155,56 +157,24 @@ class RGSpaceFunctionalityTests(unittest.TestCase):
@expand(product([(10,), (11,), (1, 1), (4, 4), (5, 7), (8, 12), (7, 16),
(4, 6, 8), (17, 5, 3)],
[True, False]))
def test_hermitian_decomposition(self, shape, zerocenter):
def test_hermitianize_inverter(self, shape, zerocenter):
r = RGSpace(shape, harmonic=True, zerocenter=zerocenter)
v = np.empty(shape, dtype=np.complex128)
v.real = np.random.random(shape)
v.imag = np.random.random(shape)
h, a = r.hermitian_decomposition(v)
# make sure that data == h + a
# NOTE: this is only correct for preserve_gaussian_variance==False,
# but I consider this an intrinsic property of a hermitian
# decomposition.
assert_almost_equal(v, h+a)
print (h, a)
# test hermitianity of h
it = np.nditer(h, flags=['multi_index'])
while not it.finished:
i1 = it.multi_index
i2 = []
for i in range(len(i1)):
if r.zerocenter[i] and r.shape[i] % 2 != 0:
i2.append(h.shape[i]-i1[i]-1)
else:
i2.append(h.shape[i]-i1[i] if i1[i] > 0 else 0)
i2 = tuple(i2)
assert_almost_equal(h[i1], np.conj(h[i2]))
assert_almost_equal(a[i1], -np.conj(a[i2]))
it.iternext()
v = distributed_data_object(global_shape=shape, dtype=np.complex128)
v[:] = np.random.random(shape) + 1j*np.random.random(shape)
inverted = r.hermitianize_inverter(v, axes=range(len(shape)))
@expand(product([(10,), (11,), (1, 1), (4, 4), (5, 7), (8, 12), (7, 16),
(4, 6, 8), (17, 5, 3)],
[True, False]))
def test_hermitian_decomposition2(self, shape, zerocenter):
r = RGSpace(shape, harmonic=True, zerocenter=zerocenter)
v = np.random.random(shape)
h, a = r.hermitian_decomposition(v)
# make sure that data == h + a
assert_almost_equal(v, h+a)
# test hermitianity of h
it = np.nditer(h, flags=['multi_index'])
# test hermitian flipping of `inverted`
it = np.nditer(v, flags=['multi_index'])
while not it.finished:
i1 = it.multi_index
i2 = []
for i in range(len(i1)):
if r.zerocenter[i] and r.shape[i] % 2 != 0:
i2.append(h.shape[i]-i1[i]-1)
i2.append(v.shape[i]-i1[i]-1)
else:
i2.append(h.shape[i]-i1[i] if i1[i] > 0 else 0)
i2.append(v.shape[i]-i1[i] if i1[i] > 0 else 0)
i2 = tuple(i2)
assert_almost_equal(h[i1], np.conj(h[i2]))
assert_almost_equal(a[i1], -np.conj(a[i2]))
assert_almost_equal(inverted[i1], v[i2])
it.iternext()
@expand(get_distance_array_configs())
......
Supports Markdown
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