Commit c6f46393 authored by Theo Steininger's avatar Theo Steininger

Merge branch 'master' into index_games2

# Conflicts:
#	nifty/spaces/power_space/power_indices.py
#	test/test_spaces/test_power_space.py
parents 0c50175d 5d5a2701
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.offline as pl
#import plotly.graph_objs as go #import plotly.graph_objs as go
...@@ -10,36 +14,37 @@ rank = comm.rank ...@@ -10,36 +14,37 @@ rank = comm.rank
if __name__ == "__main__": if __name__ == "__main__":
distribution_strategy = 'not' distribution_strategy = 'fftw'
#Setting up physical constants # Setting up physical constants
#total length of Interval or Volume the field lives on, e.g. in meters # total length of Interval or Volume the field lives on, e.g. in meters
L = 2. 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 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. 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 response_sigma = 0.1
#defining resolution (pixels per dimension) # defining resolution (pixels per dimension)
N_pixels = 512 N_pixels = 512
#Setting up derived constants # Setting up derived constants
k_0 = 1./correlation_length k_0 = 1./correlation_length
#note that field_variance**2 = a*k_0/4. for this analytic form of power # note that field_variance**2 = a*k_0/4. for this analytic form of power
#spectrum # spectrum
a = field_variance**2/k_0*4. a = field_variance**2/k_0*4.
pow_spec = (lambda k: a / (1 + k/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 # Setting up the geometry
s_space = RGSpace([N_pixels, N_pixels], distances = pixel_width) s_space = RGSpace([N_pixels, N_pixels], distances=pixel_length)
fft = FFTOperator(s_space) fft = FFTOperator(s_space, domain_dtype=np.float, target_dtype=np.complex)
h_space = fft.target[0] 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) p_space = PowerSpace(h_space, distribution_strategy=distribution_strategy)
# Creating the mock data # Creating the mock data
S = create_power_operator(h_space, power_spectrum=pow_spec, S = create_power_operator(h_space, power_spectrum=pow_spec,
...@@ -51,6 +56,7 @@ if __name__ == "__main__": ...@@ -51,6 +56,7 @@ if __name__ == "__main__":
ss = fft.inverse_times(sh) ss = fft.inverse_times(sh)
R = SmoothingOperator(s_space, sigma=response_sigma) R = SmoothingOperator(s_space, sigma=response_sigma)
R_harmonic = ComposedOperator([inverse_fft, R], default_spaces=[0, 0])
signal_to_noise = 1 signal_to_noise = 1
N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True) N = DiagonalOperator(s_space, diagonal=ss.var()/signal_to_noise, bare=True)
...@@ -63,7 +69,9 @@ if __name__ == "__main__": ...@@ -63,7 +69,9 @@ if __name__ == "__main__":
# Wiener filter # Wiener filter
j = R.adjoint_times(N.inverse_times(d)) j = R_harmonic.adjoint_times(N.inverse_times(d))
D = PropagatorOperator(S=S, N=N, R=R) 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( ...@@ -66,7 +66,7 @@ variable_default_field_dtype = keepers.Variable(
variable_default_distribution_strategy = keepers.Variable( variable_default_distribution_strategy = keepers.Variable(
'default_distribution_strategy', 'default_distribution_strategy',
['not', 'fftw', 'equal'], ['not', 'fftw', 'equal'],
lambda z: (('pyfftw' in dependency_injector) lambda z: (('fftw' in dependency_injector)
if z == 'fftw' else True), if z == 'fftw' else True),
genus='str') genus='str')
......
...@@ -596,22 +596,15 @@ class Field(Loggable, Versionable, object): ...@@ -596,22 +596,15 @@ class Field(Loggable, Versionable, object):
@staticmethod @staticmethod
def _hermitian_decomposition(domain, val, spaces, domain_axes, def _hermitian_decomposition(domain, val, spaces, domain_axes,
preserve_gaussian_variance=False): preserve_gaussian_variance=False):
# hermitianize for the first space
(h, a) = domain[spaces[0]].hermitian_decomposition( flipped_val = val
val, for space in spaces:
domain_axes[spaces[0]]) flipped_val = domain[space].hermitianize_inverter(
# hermitianize all remaining spaces using the iterative formula x=flipped_val,
for space in spaces[1:]: axes=domain_axes[space])
(hh, ha) = domain[space].hermitian_decomposition( flipped_val = flipped_val.conjugate()
h, h = (val + flipped_val)/2.
domain_axes[space]) a = val - h
(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.
# correct variance # correct variance
if preserve_gaussian_variance: if preserve_gaussian_variance:
......
...@@ -117,7 +117,6 @@ class FFTOperator(LinearOperator): ...@@ -117,7 +117,6 @@ class FFTOperator(LinearOperator):
super(FFTOperator, self).__init__(default_spaces) super(FFTOperator, self).__init__(default_spaces)
# Initialize domain and target # Initialize domain and target
self._domain = self._parse_domain(domain) self._domain = self._parse_domain(domain)
if len(self.domain) != 1: if len(self.domain) != 1:
raise ValueError("TransformationOperator accepts only exactly one " raise ValueError("TransformationOperator accepts only exactly one "
......
...@@ -8,12 +8,21 @@ from .smoothing_operator import SmoothingOperator ...@@ -8,12 +8,21 @@ from .smoothing_operator import SmoothingOperator
class FFTSmoothingOperator(SmoothingOperator): class FFTSmoothingOperator(SmoothingOperator):
def _smooth(self, x, spaces, inverse): def __init__(self, domain, sigma, log_distances=False,
Transformator = FFTOperator(x.domain[spaces[0]]) 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 # transform to the (global-)default codomain and perform all remaining
# steps therein # 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]] codomain = transformed_x.domain[spaces[0]]
coaxes = transformed_x.domain_axes[spaces[0]] coaxes = transformed_x.domain_axes[spaces[0]]
...@@ -51,9 +60,18 @@ class FFTSmoothingOperator(SmoothingOperator): ...@@ -51,9 +60,18 @@ class FFTSmoothingOperator(SmoothingOperator):
transformed_x.val.set_local_data(local_transformed_x, copy=False) 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 = x.copy_empty()
result.set_val(smoothed_x, copy=False) result.set_val(smoothed_x, copy=False)
return result 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): ...@@ -89,22 +89,6 @@ class LMSpace(Space):
super(LMSpace, self).__init__() super(LMSpace, self).__init__()
self._lmax = self._parse_lmax(lmax) 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--- # ---Mandatory properties and methods---
def __repr__(self): def __repr__(self):
......
...@@ -100,26 +100,6 @@ class RGSpace(Space): ...@@ -100,26 +100,6 @@ class RGSpace(Space):
self._distances = self._parse_distances(distances) self._distances = self._parse_distances(distances)
self._zerocenter = self._parse_zerocenter(zerocenter) 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): def hermitian_fixed_points(self):
dimensions = len(self.shape) dimensions = len(self.shape)
mid_index = np.array(self.shape)//2 mid_index = np.array(self.shape)//2
...@@ -138,7 +118,7 @@ class RGSpace(Space): ...@@ -138,7 +118,7 @@ class RGSpace(Space):
fixed_points += [tuple(index * mid_index)] fixed_points += [tuple(index * mid_index)]
return fixed_points return fixed_points
def _hermitianize_inverter(self, x, axes): def hermitianize_inverter(self, x, axes):
# calculate the number of dimensions the input array has # calculate the number of dimensions the input array has
dimensions = len(x.shape) dimensions = len(x.shape)
# prepare the slicing object which will be used for mirroring # prepare the slicing object which will be used for mirroring
......
...@@ -161,35 +161,34 @@ class Space(DomainObject): ...@@ -161,35 +161,34 @@ class Space(DomainObject):
raise NotImplementedError( raise NotImplementedError(
"There is no generic co-smoothing kernel for Space base class.") "There is no generic co-smoothing kernel for Space base class.")
def hermitian_decomposition(self, x, axes, def hermitian_fixed_points(self):
preserve_gaussian_variance=False): """ Returns the array points which remain invariant under the action
""" Decomposes x into its hermitian and anti-hermitian constituents. of `hermitianize_inverter`
This method decomposes a field's array x into its hermitian and
antihermitian part, which corresponds to real and imaginary part Returns
in a corresponding harmonic partner space. This is an internal function -------
that is mainly used for power-synthesizing and -analyzing Fields. 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 Parameters
---------- ----------
x : distributed_data_object
The field's val object.
axes : tuple of ints axes : tuple of ints
Specifies the axes of x which correspond to this space. 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 Returns
------- -------
(distributed_data_object, distributed_data_object) distributed_data_object
A tuple of two distributed_data_objects, the first being the The Hermitian-flipped of x.
hermitian and the second the anti-hermitian part of x.
""" """
raise NotImplementedError return x
...@@ -77,10 +77,13 @@ class FFTOperatorTests(unittest.TestCase): ...@@ -77,10 +77,13 @@ class FFTOperatorTests(unittest.TestCase):
b = RGSpace(dim1, zerocenter=zc2, distances=1./(dim1*d), harmonic=True) b = RGSpace(dim1, zerocenter=zc2, distances=1./(dim1*d), harmonic=True)
fft = FFTOperator(domain=a, target=b, domain_dtype=itp, fft = FFTOperator(domain=a, target=b, domain_dtype=itp,
target_dtype=_harmonic_type(itp), module=module) target_dtype=_harmonic_type(itp), module=module)
np.random.seed(16)
inp = Field.from_random(domain=a, random_type='normal', std=7, mean=3, inp = Field.from_random(domain=a, random_type='normal', std=7, mean=3,
dtype=itp) dtype=itp)
out = fft.adjoint_times(fft.times(inp)) out = fft.adjoint_times(fft.times(inp))
assert_allclose(inp.val, out.val, rtol=tol, atol=tol) assert_allclose(inp.val.get_full_data(),
out.val.get_full_data(),
rtol=tol, atol=tol)
@expand(product(["numpy", "fftw", "fftw_mpi"], @expand(product(["numpy", "fftw", "fftw_mpi"],
[10, 11], [9, 12], [False, True], [10, 11], [9, 12], [False, True],
......
...@@ -20,7 +20,7 @@ import unittest ...@@ -20,7 +20,7 @@ import unittest
import numpy as np import numpy as np
from numpy.testing import assert_, assert_equal, assert_raises,\ 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 d2o import distributed_data_object
from nifty import LMSpace from nifty import LMSpace
from test.common import expand from test.common import expand
...@@ -108,15 +108,12 @@ class LMSpaceFunctionalityTests(unittest.TestCase): ...@@ -108,15 +108,12 @@ class LMSpaceFunctionalityTests(unittest.TestCase):
for key, value in expected.iteritems(): for key, value in expected.iteritems():
assert_equal(getattr(l, key), value) assert_equal(getattr(l, key), value)
@expand(get_hermitian_configs()) def test_hermitianize_inverter(self):
def test_hermitian_decomposition(self, x, real, imag):
l = LMSpace(5) l = LMSpace(5)
assert_almost_equal( v = distributed_data_object(global_shape=l.shape, dtype=np.complex128)
l.hermitian_decomposition(distributed_data_object(x))[0], v[:] = np.random.random(l.shape) + 1j*np.random.random(l.shape)
real) inverted = l.hermitianize_inverter(v, axes=(0,))
assert_almost_equal( assert_array_almost_equal(inverted.get_full_data(), v.get_full_data())
l.hermitian_decomposition(distributed_data_object(x))[1],
imag*1j)
@expand(get_weight_configs()) @expand(get_weight_configs())
def test_weight(self, x, power, axes, inplace, expected): def test_weight(self, x, power, axes, inplace, expected):
......
...@@ -21,6 +21,8 @@ from __future__ import division ...@@ -21,6 +21,8 @@ from __future__ import division
import unittest import unittest
import numpy as np import numpy as np
from d2o import distributed_data_object
from numpy.testing import assert_, assert_equal, assert_almost_equal from numpy.testing import assert_, assert_equal, assert_almost_equal
from nifty import RGSpace from nifty import RGSpace
from test.common import expand from test.common import expand
...@@ -155,56 +157,24 @@ class RGSpaceFunctionalityTests(unittest.TestCase): ...@@ -155,56 +157,24 @@ class RGSpaceFunctionalityTests(unittest.TestCase):
@expand(product([(10,), (11,), (1, 1), (4, 4), (5, 7), (8, 12), (7, 16), @expand(product([(10,), (11,), (1, 1), (4, 4), (5, 7), (8, 12), (7, 16),
(4, 6, 8), (17, 5, 3)], (4, 6, 8), (17, 5, 3)],
[True, False])) [True, False]))
def test_hermitian_decomposition(self, shape, zerocenter): def test_hermitianize_inverter(self, shape, zerocenter):
r = RGSpace(shape, harmonic=True, zerocenter=zerocenter) r = RGSpace(shape, harmonic=True, zerocenter=zerocenter)
v = np.empty(shape, dtype=np.complex128) v = distributed_data_object(global_shape=shape, dtype=np.complex128)
v.real = np.random.random(shape) v[:] = np.random.random(shape) + 1j*np.random.random(shape)
v.imag = np.random.random(shape) inverted = r.hermitianize_inverter(v, axes=range(len(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()
@expand(product([(10,), (11,), (1, 1), (4, 4), (5, 7), (8, 12), (7, 16), # test hermitian flipping of `inverted`
(4, 6, 8), (17, 5, 3)], it = np.nditer(v, flags=['multi_index'])
[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'])
while not it.finished: while not it.finished:
i1 = it.multi_index i1 = it.multi_index
i2 = [] i2 = []
for i in range(len(i1)): for i in range(len(i1)):
if r.zerocenter[i] and r.shape[i] % 2 != 0: 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: 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) i2 = tuple(i2)
assert_almost_equal(h[i1], np.conj(h[i2])) assert_almost_equal(inverted[i1], v[i2])
assert_almost_equal(a[i1], -np.conj(a[i2]))
it.iternext() it.iternext()
@expand(get_distance_array_configs()) @expand(get_distance_array_configs())
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment