Commit 8e1151de authored by Martin Reinecke's avatar Martin Reinecke
Browse files

further adjustments

parent 85e64e48
Pipeline #17545 passed with stage
in 18 minutes and 6 seconds
...@@ -282,8 +282,8 @@ class Field(Loggable, Versionable, object): ...@@ -282,8 +282,8 @@ class Field(Loggable, Versionable, object):
# ---Powerspectral methods--- # ---Powerspectral methods---
def power_analyze(self, spaces=None, logarithmic=None, nbin=None, def power_analyze(self, spaces=None, binbounds=None,
binbounds=None, keep_phase_information=False): keep_phase_information=False):
""" Computes the square root power spectrum for a subspace of `self`. """ Computes the square root power spectrum for a subspace of `self`.
Creates a PowerSpace for the space addressed by `spaces` with the given Creates a PowerSpace for the space addressed by `spaces` with the given
...@@ -297,16 +297,8 @@ class Field(Loggable, Versionable, object): ...@@ -297,16 +297,8 @@ class Field(Loggable, Versionable, object):
spaces : int *optional* spaces : int *optional*
The subspace for which the powerspectrum shall be computed The subspace for which the powerspectrum shall be computed
(default : None). (default : None).
logarithmic : boolean *optional*
True if the output PowerSpace should use logarithmic binning.
{default : None}
nbin : int *optional*
The number of bins the resulting PowerSpace shall have
(default : None).
if nbin==None : maximum number of bins is used
binbounds : array-like *optional* binbounds : array-like *optional*
Inner bounds of the bins (default : None). Inner bounds of the bins (default : None).
Overrides nbin and logarithmic.
if binbounds==None : bins are inferred. if binbounds==None : bins are inferred.
keep_phase_information : boolean, *optional* keep_phase_information : boolean, *optional*
If False, return a real-valued result containing the power spectrum If False, return a real-valued result containing the power spectrum
...@@ -374,8 +366,6 @@ class Field(Loggable, Versionable, object): ...@@ -374,8 +366,6 @@ class Field(Loggable, Versionable, object):
parts = [self._single_power_analyze( parts = [self._single_power_analyze(
work_field=part, work_field=part,
space_index=space_index, space_index=space_index,
logarithmic=logarithmic,
nbin=nbin,
binbounds=binbounds) binbounds=binbounds)
for part in parts] for part in parts]
...@@ -387,8 +377,7 @@ class Field(Loggable, Versionable, object): ...@@ -387,8 +377,7 @@ class Field(Loggable, Versionable, object):
return result_field return result_field
@classmethod @classmethod
def _single_power_analyze(cls, work_field, space_index, logarithmic, nbin, def _single_power_analyze(cls, work_field, space_index, binbounds):
binbounds):
if not work_field.domain[space_index].harmonic: if not work_field.domain[space_index].harmonic:
raise ValueError( raise ValueError(
...@@ -407,7 +396,6 @@ class Field(Loggable, Versionable, object): ...@@ -407,7 +396,6 @@ class Field(Loggable, Versionable, object):
harmonic_domain = work_field.domain[space_index] harmonic_domain = work_field.domain[space_index]
power_domain = PowerSpace(harmonic_partner=harmonic_domain, power_domain = PowerSpace(harmonic_partner=harmonic_domain,
distribution_strategy=distribution_strategy, distribution_strategy=distribution_strategy,
logarithmic=logarithmic, nbin=nbin,
binbounds=binbounds) binbounds=binbounds)
power_spectrum = cls._calculate_power_spectrum( power_spectrum = cls._calculate_power_spectrum(
field_val=work_field.val, field_val=work_field.val,
......
...@@ -145,12 +145,6 @@ class GLSpace(Space): ...@@ -145,12 +145,6 @@ class GLSpace(Space):
return result_x return result_x
def get_distance_array(self, distribution_strategy):
raise NotImplementedError
def get_fft_smoothing_kernel_function(self, sigma):
raise NotImplementedError
# ---Added properties and methods--- # ---Added properties and methods---
@property @property
......
...@@ -119,12 +119,6 @@ class HPSpace(Space): ...@@ -119,12 +119,6 @@ class HPSpace(Space):
return result_x return result_x
def get_distance_array(self, distribution_strategy):
raise NotImplementedError
def get_fft_smoothing_kernel_function(self, sigma):
raise NotImplementedError
# ---Added properties and methods--- # ---Added properties and methods---
@property @property
......
...@@ -101,12 +101,7 @@ class PowerSpace(Space): ...@@ -101,12 +101,7 @@ class PowerSpace(Space):
""" """
nbin = int(nbin) nbin = int(nbin)
assert nbin>=3, "nbin must be at least 3" assert nbin>=3, "nbin must be at least 3"
first_bound = float(first_bound) return np.linspace(float(first_bound),float(last_bound),nbin-1)
last_bound = float(last_bound)
binbounds = np.arange(nbin-1,dtype=np.float64)/(nbin-2)
binbounds*=last_bound-first_bound
binbounds+=first_bound
return binbounds
@staticmethod @staticmethod
def logarithmic_binbounds (nbin, first_bound, last_bound): def logarithmic_binbounds (nbin, first_bound, last_bound):
...@@ -122,7 +117,11 @@ class PowerSpace(Space): ...@@ -122,7 +117,11 @@ class PowerSpace(Space):
values equidistantly spaced (in natural logarithmic scale) values equidistantly spaced (in natural logarithmic scale)
between these two. between these two.
""" """
return np.exp(linear_binbounds(np.log(first_bound),np.log(last_bound))) nbin = int(nbin)
assert nbin>=3, "nbin must be at least 3"
return np.logspace(np.log(float(first_bound)),
np.log(float(last_bound)),
nbin-1, base=np.e)
def __init__(self, harmonic_partner, distribution_strategy=None, def __init__(self, harmonic_partner, distribution_strategy=None,
binbounds=None): binbounds=None):
......
...@@ -210,6 +210,7 @@ class RGSpace(Space): ...@@ -210,6 +210,7 @@ class RGSpace(Space):
""" """
if (not self.harmonic): raise NotImplementedError
shape = self.shape shape = self.shape
# prepare the distributed_data_object # prepare the distributed_data_object
nkdict = distributed_data_object( nkdict = distributed_data_object(
...@@ -257,6 +258,7 @@ class RGSpace(Space): ...@@ -257,6 +258,7 @@ class RGSpace(Space):
return dists return dists
def get_unique_distances(self): def get_unique_distances(self):
if (not self.harmonic): raise NotImplementedError
dimensions = len(self.shape) dimensions = len(self.shape)
if dimensions == 1: # extra easy if dimensions == 1: # extra easy
maxdist = self.shape[0]//2 maxdist = self.shape[0]//2
...@@ -283,10 +285,12 @@ class RGSpace(Space): ...@@ -283,10 +285,12 @@ class RGSpace(Space):
return tmp[np.diff(np.r_[tmp, 2*tmp[-1]]) > tol] return tmp[np.diff(np.r_[tmp, 2*tmp[-1]]) > tol]
def get_natural_binbounds(self): def get_natural_binbounds(self):
if (not self.harmonic): raise NotImplementedError
tmp = self.get_unique_distances() tmp = self.get_unique_distances()
return 0.5*(tmp[:-1]+tmp[1:]) return 0.5*(tmp[:-1]+tmp[1:])
def get_fft_smoothing_kernel_function(self, sigma): def get_fft_smoothing_kernel_function(self, sigma):
if (not self.harmonic): raise NotImplementedError
return lambda x: np.exp(-2. * np.pi*np.pi * x*x * sigma*sigma) return lambda x: np.exp(-2. * np.pi*np.pi * x*x * sigma*sigma)
# ---Added properties and methods--- # ---Added properties and methods---
......
...@@ -96,10 +96,8 @@ class Space(DomainObject): ...@@ -96,10 +96,8 @@ class Space(DomainObject):
return self.__class__() return self.__class__()
def get_distance_array(self, distribution_strategy): def get_distance_array(self, distribution_strategy):
""" The distances of the pixel to zero. """ The length of the k vector for every pixel.
This method is only implemented for harmonic spaces.
This returns an array that gives for each pixel its distance to the
center of the manifolds grid.
Parameters Parameters
---------- ----------
...@@ -110,18 +108,21 @@ class Space(DomainObject): ...@@ -110,18 +108,21 @@ class Space(DomainObject):
Returns Returns
------- -------
distributed_data_object distributed_data_object
A d2o containing the distances A d2o containing the k vector lengths
""" """
raise NotImplementedError
raise NotImplementedError(
"There is no generic distance structure for Space base class.")
def get_unique_distances(self): def get_unique_distances(self):
""" Returns an array of floats containing the unique k vector lengths
for this space.
This method is only implemented for harmonic spaces.
"""
raise NotImplementedError raise NotImplementedError
def get_natural_binbounds(self): def get_natural_binbounds(self):
""" The boundaries for natural power spectrum binning. """ The boundaries for natural power spectrum binning.
This method is only implemented for harmonic spaces.
Returns Returns
------- -------
......
...@@ -18,8 +18,9 @@ ...@@ -18,8 +18,9 @@
from builtins import str from builtins import str
from parameterized import parameterized from parameterized import parameterized
from nifty import RGSpace, LMSpace, HPSpace, GLSpace, PowerSpace from nifty import Space, RGSpace, LMSpace, HPSpace, GLSpace, PowerSpace
from nifty.config import dependency_injector as gdi from nifty.config import dependency_injector as gdi
import numpy as np
def custom_name_func(testcase_func, param_num, param): def custom_name_func(testcase_func, param_num, param):
...@@ -45,3 +46,34 @@ def generate_spaces(): ...@@ -45,3 +46,34 @@ def generate_spaces():
def generate_harmonic_spaces(): def generate_harmonic_spaces():
spaces = [RGSpace(4, harmonic=True), LMSpace(5)] spaces = [RGSpace(4, harmonic=True), LMSpace(5)]
return spaces return spaces
def marco_binbounds (space, logarithmic, nbin=None):
"""Only for testing purposes. DO NOT USE IN REAL LIFE!"""
if logarithmic is None and nbin is None:
return None
if not (isinstance(space, Space) and space.harmonic):
raise ValueError("space must be a harmonic space.")
logarithmic = bool(logarithmic)
if nbin is not None:
nbin=int(nbin)
assert nbin>=3, "nbin must be at least 3"
# equidistant binning (linear or log)
# MR FIXME: this needs to improve
kindex = space.get_unique_distances()
if (logarithmic):
k = np.r_[0, np.log(kindex[1:])]
else:
k = kindex
dk = np.max(k[2:] - k[1:-1]) # minimum dk to avoid empty bins
if(nbin is None):
nbin = int((k[-1] - 0.5 * (k[2] + k[1])) /
dk - 0.5) # maximal nbin
else:
nbin = min(int(nbin), int(
(k[-1] - 0.5 * (k[2] + k[1])) / dk + 2.5))
dk = (k[-1] - 0.5 * (k[2] + k[1])) / (nbin - 2.5)
bb = np.r_[0.5 * (3 * k[1] - k[2]),
0.5 * (k[1] + k[2]) + dk * np.arange(nbin-2)]
if(logarithmic):
bb = np.exp(bb)
return bb
...@@ -21,14 +21,14 @@ import numpy as np ...@@ -21,14 +21,14 @@ import numpy as np
import nifty as ift import nifty as ift
from numpy.testing import assert_allclose from numpy.testing import assert_allclose
from itertools import product from itertools import product
from test.common import expand from test.common import expand, marco_binbounds
class LaplaceOperatorTests(unittest.TestCase): class LaplaceOperatorTests(unittest.TestCase):
@expand(product([None, False, True], [False, True], [10, 100, 1000])) @expand(product([None, False, True], [False, True], [10, 100, 1000]))
def test_Laplace(self, log1, log2, sz): def test_Laplace(self, log1, log2, sz):
s = ift.RGSpace(sz, harmonic=True) s = ift.RGSpace(sz, harmonic=True)
p = ift.PowerSpace(s, logarithmic=log1) p = ift.PowerSpace(s, binbounds=marco_binbounds(s,logarithmic=log1))
L = ift.LaplaceOperator(p, logarithmic=log2) L = ift.LaplaceOperator(p, logarithmic=log2)
arr = np.random.random(p.shape[0]) arr = np.random.random(p.shape[0])
fp = ift.Field(p, val=arr) fp = ift.Field(p, val=arr)
......
...@@ -27,7 +27,7 @@ from nifty import Field,\ ...@@ -27,7 +27,7 @@ from nifty import Field,\
SmoothingOperator SmoothingOperator
from itertools import product from itertools import product
from test.common import expand from test.common import expand, marco_binbounds
def _get_rtol(tp): def _get_rtol(tp):
...@@ -87,7 +87,7 @@ class SmoothingOperator_Tests(unittest.TestCase): ...@@ -87,7 +87,7 @@ class SmoothingOperator_Tests(unittest.TestCase):
[np.float64, np.complex128])) [np.float64, np.complex128]))
def test_smooth_regular1(self, sz, d, sigma, tp): def test_smooth_regular1(self, sz, d, sigma, tp):
tol = _get_rtol(tp) tol = _get_rtol(tp)
sp = RGSpace(sz, harmonic=True, distances=d) sp = RGSpace(sz, harmonic=False, distances=d)
smo = SmoothingOperator.make(sp, sigma=sigma) smo = SmoothingOperator.make(sp, sigma=sigma)
inp = Field.from_random(domain=sp, random_type='normal', std=1, mean=4, inp = Field.from_random(domain=sp, random_type='normal', std=1, mean=4,
dtype=tp) dtype=tp)
...@@ -99,7 +99,7 @@ class SmoothingOperator_Tests(unittest.TestCase): ...@@ -99,7 +99,7 @@ class SmoothingOperator_Tests(unittest.TestCase):
[np.float64, np.complex128])) [np.float64, np.complex128]))
def test_smooth_regular2(self, sz1, sz2, d1, d2, sigma, tp): def test_smooth_regular2(self, sz1, sz2, d1, d2, sigma, tp):
tol = _get_rtol(tp) tol = _get_rtol(tp)
sp = RGSpace([sz1, sz2], distances=[d1, d2], harmonic=True) sp = RGSpace([sz1, sz2], distances=[d1, d2], harmonic=False)
smo = SmoothingOperator.make(sp, sigma=sigma) smo = SmoothingOperator.make(sp, sigma=sigma)
inp = Field.from_random(domain=sp, random_type='normal', std=1, mean=4, inp = Field.from_random(domain=sp, random_type='normal', std=1, mean=4,
dtype=tp) dtype=tp)
...@@ -112,7 +112,8 @@ class SmoothingOperator_Tests(unittest.TestCase): ...@@ -112,7 +112,8 @@ class SmoothingOperator_Tests(unittest.TestCase):
def test_smooth_irregular1(self, sz, log, sigma, tp): def test_smooth_irregular1(self, sz, log, sigma, tp):
tol = _get_rtol(tp) tol = _get_rtol(tp)
sp = RGSpace(sz, harmonic=True) sp = RGSpace(sz, harmonic=True)
ps = PowerSpace(sp, nbin=sz, logarithmic=log) ps = PowerSpace(sp, binbounds=marco_binbounds(
sp, logarithmic=log, nbin=sz))
smo = SmoothingOperator.make(ps, sigma=sigma) smo = SmoothingOperator.make(ps, sigma=sigma)
inp = Field.from_random(domain=ps, random_type='normal', std=1, mean=4, inp = Field.from_random(domain=ps, random_type='normal', std=1, mean=4,
dtype=tp) dtype=tp)
...@@ -125,7 +126,7 @@ class SmoothingOperator_Tests(unittest.TestCase): ...@@ -125,7 +126,7 @@ class SmoothingOperator_Tests(unittest.TestCase):
def test_smooth_irregular2(self, sz1, sz2, log, sigma, tp): def test_smooth_irregular2(self, sz1, sz2, log, sigma, tp):
tol = _get_rtol(tp) tol = _get_rtol(tp)
sp = RGSpace([sz1, sz2], harmonic=True) sp = RGSpace([sz1, sz2], harmonic=True)
ps = PowerSpace(sp, logarithmic=log) ps = PowerSpace(sp, binbounds=marco_binbounds(sp, logarithmic=log))
smo = SmoothingOperator.make(ps, sigma=sigma) smo = SmoothingOperator.make(ps, sigma=sigma)
inp = Field.from_random(domain=ps, random_type='normal', std=1, mean=4, inp = Field.from_random(domain=ps, random_type='normal', std=1, mean=4,
dtype=tp) dtype=tp)
......
...@@ -25,7 +25,7 @@ from d2o import distributed_data_object ...@@ -25,7 +25,7 @@ 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,\
assert_raises assert_raises
from nifty import PowerSpace, RGSpace, Space, LMSpace from nifty import PowerSpace, RGSpace, Space, LMSpace
from test.common import expand from test.common import expand, marco_binbounds
from itertools import product, chain from itertools import product, chain
# needed to check wether fftw is available # needed to check wether fftw is available
from nifty import dependency_injector as gdi from nifty import dependency_injector as gdi
...@@ -41,7 +41,7 @@ HARMONIC_SPACES = [RGSpace((8,), harmonic=True), ...@@ -41,7 +41,7 @@ HARMONIC_SPACES = [RGSpace((8,), harmonic=True),
LMSpace(9)] LMSpace(9)]
#Try all sensible kinds of combinations of spaces, distributuion strategy and #Try all sensible kinds of combinations of spaces, distribution strategy and
#binning parameters #binning parameters
CONSISTENCY_CONFIGS_IMPLICIT = product(HARMONIC_SPACES, CONSISTENCY_CONFIGS_IMPLICIT = product(HARMONIC_SPACES,
["not", "equal", "fftw"], ["not", "equal", "fftw"],
...@@ -130,8 +130,8 @@ class PowerSpaceConsistencyCheck(unittest.TestCase): ...@@ -130,8 +130,8 @@ class PowerSpaceConsistencyCheck(unittest.TestCase):
raise SkipTest raise SkipTest
p = PowerSpace(harmonic_partner=harmonic_partner, p = PowerSpace(harmonic_partner=harmonic_partner,
distribution_strategy=distribution_strategy, distribution_strategy=distribution_strategy,
logarithmic=logarithmic, nbin=nbin, binbounds=marco_binbounds(harmonic_partner,
binbounds=binbounds) logarithmic, nbin))
assert_equal(p.pindex.flatten().bincount(), p.rho, assert_equal(p.pindex.flatten().bincount(), p.rho,
err_msg='rho is not equal to pindex degeneracy') err_msg='rho is not equal to pindex degeneracy')
...@@ -147,13 +147,13 @@ class PowerSpaceFunctionalityTest(unittest.TestCase): ...@@ -147,13 +147,13 @@ class PowerSpaceFunctionalityTest(unittest.TestCase):
with assert_raises(expected['error']): with assert_raises(expected['error']):
PowerSpace(harmonic_partner=harmonic_partner, PowerSpace(harmonic_partner=harmonic_partner,
distribution_strategy=distribution_strategy, distribution_strategy=distribution_strategy,
logarithmic=logarithmic, nbin=nbin, binbounds=marco_binbounds(harmonic_partner,
binbounds=binbounds) logarithmic, nbin))
else: else:
p = PowerSpace(harmonic_partner=harmonic_partner, p = PowerSpace(harmonic_partner=harmonic_partner,
distribution_strategy=distribution_strategy, distribution_strategy=distribution_strategy,
logarithmic=logarithmic, nbin=nbin, binbounds=marco_binbounds(harmonic_partner,
binbounds=binbounds) logarithmic, nbin))
for key, value in expected.items(): for key, value in expected.items():
if isinstance(value, np.ndarray): if isinstance(value, np.ndarray):
assert_almost_equal(getattr(p, key), value) assert_almost_equal(getattr(p, key), value)
......
...@@ -91,7 +91,8 @@ CONSTRUCTOR_CONFIGS = [ ...@@ -91,7 +91,8 @@ CONSTRUCTOR_CONFIGS = [
def get_distance_array_configs(): def get_distance_array_configs():
# for RGSpace(shape=(4, 4), distances=None, zerocenter=[False, False]) # for RGSpace(shape=(4, 4), distances=(0.25,0.25),
# zerocenter=[False, False])
cords_0 = np.ogrid[0:4, 0:4] cords_0 = np.ogrid[0:4, 0:4]
da_0 = ((cords_0[0] - 4 // 2) * 0.25)**2 da_0 = ((cords_0[0] - 4 // 2) * 0.25)**2
da_0 = np.fft.ifftshift(da_0) da_0 = np.fft.ifftshift(da_0)
...@@ -99,7 +100,7 @@ def get_distance_array_configs(): ...@@ -99,7 +100,7 @@ def get_distance_array_configs():
temp = np.fft.ifftshift(temp) temp = np.fft.ifftshift(temp)
da_0 = da_0 + temp da_0 = da_0 + temp
da_0 = np.sqrt(da_0) da_0 = np.sqrt(da_0)
# for RGSpace(shape=(4, 4), distances=None, zerocenter=[True, True]) # for RGSpace(shape=(4, 4), distances=(0.25,0.25), zerocenter=[True, True])
da_1 = ((cords_0[0] - 4 // 2) * 0.25)**2 da_1 = ((cords_0[0] - 4 // 2) * 0.25)**2
temp = ((cords_0[1] - 4 // 2) * 0.25)**2 temp = ((cords_0[1] - 4 // 2) * 0.25)**2
da_1 = da_1 + temp da_1 = da_1 + temp
...@@ -110,8 +111,8 @@ def get_distance_array_configs(): ...@@ -110,8 +111,8 @@ def get_distance_array_configs():
da_2 = da_2 + temp da_2 = da_2 + temp
da_2 = np.sqrt(da_2) da_2 = np.sqrt(da_2)
return [ return [
[(4, 4), None, [False, False], da_0], [(4, 4), (0.25, 0.25), [False, False], da_0],
[(4, 4), None, [True, True], da_1], [(4, 4), (0.25, 0.25), [True, True], da_1],
[(4, 4), (12, 12), [True, True], da_2] [(4, 4), (12, 12), [True, True], da_2]
] ]
...@@ -190,7 +191,8 @@ class RGSpaceFunctionalityTests(unittest.TestCase): ...@@ -190,7 +191,8 @@ class RGSpaceFunctionalityTests(unittest.TestCase):
@expand(get_distance_array_configs()) @expand(get_distance_array_configs())
def test_distance_array(self, shape, distances, zerocenter, expected): def test_distance_array(self, shape, distances, zerocenter, expected):
r = RGSpace(shape=shape, distances=distances, zerocenter=zerocenter) r = RGSpace(shape=shape, distances=distances, zerocenter=zerocenter,
harmonic=True)
assert_almost_equal(r.get_distance_array('not'), expected) assert_almost_equal(r.get_distance_array('not'), expected)
@expand(get_weight_configs()) @expand(get_weight_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