Commit ada043fa authored by Theo Steininger's avatar Theo Steininger
Browse files

Merge branch 'fix_powerspace' into 'nightly'

Enforce unique specifications for PowerSpace

See merge request !195
parents c5b2e6fb d7d9b42d
......@@ -281,8 +281,8 @@ class Field(Loggable, Versionable, object):
# ---Powerspectral methods---
def power_analyze(self, spaces=None, logarithmic=None, nbin=None,
binbounds=None, keep_phase_information=False):
def power_analyze(self, spaces=None, binbounds=None,
keep_phase_information=False):
""" Computes the square root power spectrum for a subspace of `self`.
Creates a PowerSpace for the space addressed by `spaces` with the given
......@@ -296,16 +296,8 @@ class Field(Loggable, Versionable, object):
spaces : int *optional*
The subspace for which the powerspectrum shall be computed
(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*
Inner bounds of the bins (default : None).
Overrides nbin and logarithmic.
if binbounds==None : bins are inferred.
keep_phase_information : boolean, *optional*
If False, return a real-valued result containing the power spectrum
......@@ -373,8 +365,6 @@ class Field(Loggable, Versionable, object):
parts = [self._single_power_analyze(
work_field=part,
space_index=space_index,
logarithmic=logarithmic,
nbin=nbin,
binbounds=binbounds)
for part in parts]
......@@ -386,8 +376,7 @@ class Field(Loggable, Versionable, object):
return result_field
@classmethod
def _single_power_analyze(cls, work_field, space_index, logarithmic, nbin,
binbounds):
def _single_power_analyze(cls, work_field, space_index, binbounds):
if not work_field.domain[space_index].harmonic:
raise ValueError(
......@@ -406,7 +395,6 @@ class Field(Loggable, Versionable, object):
harmonic_domain = work_field.domain[space_index]
power_domain = PowerSpace(harmonic_partner=harmonic_domain,
distribution_strategy=distribution_strategy,
logarithmic=logarithmic, nbin=nbin,
binbounds=binbounds)
power_spectrum = cls._calculate_power_spectrum(
field_val=work_field.val,
......
......@@ -17,8 +17,6 @@ def buildIdx(nr, lmax):
new_dtype = np.float32
elif nr.dtype == np.dtype('complex128'):
new_dtype = np.float64
elif nr.dtype == np.dtype('complex256'):
new_dtype = np.float128
else:
raise TypeError("dtype of nr not supported.")
......
......@@ -145,12 +145,6 @@ class GLSpace(Space):
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---
@property
......
......@@ -119,12 +119,6 @@ class HPSpace(Space):
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---
@property
......
......@@ -38,25 +38,19 @@ class PowerSpace(Space):
The distribution strategy used for the distributed_data_objects
derived from this PowerSpace, e.g. the pindex.
(default : 'not')
logarithmic : bool *optional*
True if logarithmic binning should be used (default : None).
nbin : {int, None} *optional*
The number of bins that should be used for power spectrum binning
(default : None).
if nbin == None, then nbin is set to the length of kindex.
binbounds : {list, array-like} *optional*
Boundaries between the power spectrum bins.
(If binbounds has n entries, there will be n+1 bins, the first bin
starting at -inf, the last bin ending at +inf.)
binbounds: None, or tuple/array/list of float
if None:
There will be as many bins as there are distinct k-vector lengths
in the harmonic partner space.
The "binbounds" property of the PowerSpace will also be None.
else:
the bin bounds requested for this PowerSpace. The array
must be sorted and strictly ascending. The first entry is the right
boundary of the first bin, and the last entry is the left boundary
of the last bin, i.e. thee will be len(binbounds)+1 bins in total,
with the first and last bins reaching to -+infinity, respectively.
(default : None)
if binbounds == None:
Calculates the bounds from the kindex while applying the
logarithmic and nbin keywords.
Note: if "bindounds" is not None, both "logarithmic" and "nbin" must be
None!
Note: if "binbounds", "logarithmic", and "nbin" are all None, then
"natural" binning is performed, i.e. there will be one bin for every
distinct k-vector length.
Attributes
----------
......@@ -71,13 +65,14 @@ class PowerSpace(Space):
dim : np.int
Total number of dimensionality, i.e. the number of pixels.
harmonic : bool
Specifies whether the space is a signal or harmonic space.
Always True for this space.
total_volume : np.float
The total volume of the space.
shape : tuple of np.ints
The shape of the space's data array.
binbounds : tuple or None
Boundaries between the power spectrum bins
binbounds : tuple or None
Boundaries between the power spectrum bins; None is used to indicate
natural binning
Notes
-----
......@@ -90,8 +85,74 @@ class PowerSpace(Space):
# ---Overwritten properties and methods---
@staticmethod
def linear_binbounds(nbin, first_bound, last_bound):
"""
nbin: integer
the number of bins
first_bound, last_bound: float
the k values for the right boundary of the first bin and the left
boundary of the last bin, respectively. They are given in length
units of the harmonic partner space.
This will produce a binbounds array with nbin-1 entries with
binbounds[0]=first_bound and binbounds[-1]=last_bound and the remaining
values equidistantly spaced (in linear scale) between these two.
"""
nbin = int(nbin)
assert nbin >= 3, "nbin must be at least 3"
return np.linspace(float(first_bound), float(last_bound), nbin-1)
@staticmethod
def logarithmic_binbounds(nbin, first_bound, last_bound):
"""
nbin: integer
the number of bins
first_bound, last_bound: float
the k values for the right boundary of the first bin and the left
boundary of the last bin, respectively. They are given in length
units of the harmonic partner space.
This will produce a binbounds array with nbin-1 entries with
binbounds[0]=first_bound and binbounds[-1]=last_bound and the remaining
values equidistantly spaced (in natural logarithmic scale)
between these two.
"""
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)
@staticmethod
def useful_binbounds(space, logarithmic, nbin=None):
if not (isinstance(space, Space) and space.harmonic):
raise ValueError("first argument must be a harmonic space.")
if logarithmic is None and nbin is None:
return None
nbin = None if nbin is None else int(nbin)
logarithmic = bool(logarithmic)
dists = space.get_unique_distances()
if len(dists) < 3:
raise ValueError("Space does not have enough unique k lengths")
lbound = 0.5*(dists[0]+dists[1])
rbound = 0.5*(dists[-2]+dists[-1])
dists[0] = lbound
dists[-1] = rbound
if logarithmic:
dists = np.log(dists)
binsz_min = np.max(np.diff(dists))
nbin_max = int((dists[-1]-dists[0])/binsz_min)+2
if nbin is None:
nbin = nbin_max
assert nbin >= 3, "nbin must be at least 3"
if nbin > nbin_max:
raise ValueError("nbin is too large")
if logarithmic:
return PowerSpace.logarithmic_binbounds(nbin, lbound, rbound)
else:
return PowerSpace.linear_binbounds(nbin, lbound, rbound)
def __init__(self, harmonic_partner, distribution_strategy=None,
logarithmic=None, nbin=None, binbounds=None):
binbounds=None):
super(PowerSpace, self).__init__()
self._ignore_for_hash += ['_pindex', '_kindex', '_rho']
......@@ -107,35 +168,24 @@ class PowerSpace(Space):
raise ValueError("harmonic_partner must be a harmonic space.")
self._harmonic_partner = harmonic_partner
# sanity check
if binbounds is not None and not(nbin is None and logarithmic is None):
raise ValueError(
"if binbounds is defined, nbin and logarithmic must be None")
if binbounds is not None:
binbounds = tuple(binbounds)
key = (harmonic_partner, distribution_strategy, logarithmic, nbin,
binbounds)
key = (harmonic_partner, distribution_strategy, binbounds)
if self._powerIndexCache.get(key) is None:
distance_array = \
self.harmonic_partner.get_distance_array(distribution_strategy)
temp_binbounds = self._compute_binbounds(
harmonic_partner=self.harmonic_partner,
distribution_strategy=distribution_strategy,
logarithmic=logarithmic,
nbin=nbin,
binbounds=binbounds)
temp_pindex = self._compute_pindex(
harmonic_partner=self.harmonic_partner,
distance_array=distance_array,
binbounds=temp_binbounds,
binbounds=binbounds,
distribution_strategy=distribution_strategy)
temp_rho = temp_pindex.bincount().get_full_data()
assert not np.any(temp_rho == 0), "empty bins detected"
temp_kindex = \
(temp_pindex.bincount(weights=distance_array).get_full_data() /
temp_rho)
self._powerIndexCache[key] = (temp_binbounds,
self._powerIndexCache[key] = (binbounds,
temp_pindex,
temp_kindex,
temp_rho)
......@@ -143,43 +193,6 @@ class PowerSpace(Space):
(self._binbounds, self._pindex, self._kindex, self._rho) = \
self._powerIndexCache[key]
@staticmethod
def _compute_binbounds(harmonic_partner, distribution_strategy,
logarithmic, nbin, binbounds):
if logarithmic is None and nbin is None and binbounds is None:
result = None
else:
if binbounds is not None:
bb = np.sort(np.array(binbounds))
else:
if logarithmic is not None:
logarithmic = bool(logarithmic)
if nbin is not None:
nbin = int(nbin)
# equidistant binning (linear or log)
# MR FIXME: this needs to improve
kindex = harmonic_partner.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)
result = tuple(bb)
return result
@staticmethod
def _compute_pindex(harmonic_partner, distance_array, binbounds,
distribution_strategy):
......
......@@ -125,6 +125,8 @@ class RGSpace(Space):
# return fixed_points
def hermitianize_inverter(self, x, axes):
if (not self.harmonic):
raise NotImplementedError
if nifty_configuration['harmonic_rg_base'] == 'real':
return x
else:
......@@ -210,6 +212,8 @@ class RGSpace(Space):
"""
if (not self.harmonic):
raise NotImplementedError
shape = self.shape
# prepare the distributed_data_object
nkdict = distributed_data_object(
......@@ -257,6 +261,8 @@ class RGSpace(Space):
return dists
def get_unique_distances(self):
if (not self.harmonic):
raise NotImplementedError
dimensions = len(self.shape)
if dimensions == 1: # extra easy
maxdist = self.shape[0]//2
......@@ -283,10 +289,14 @@ class RGSpace(Space):
return tmp[np.diff(np.r_[tmp, 2*tmp[-1]]) > tol]
def get_natural_binbounds(self):
if (not self.harmonic):
raise NotImplementedError
tmp = self.get_unique_distances()
return 0.5*(tmp[:-1]+tmp[1:])
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)
# ---Added properties and methods---
......
......@@ -51,21 +51,18 @@ class Space(DomainObject):
Notes
-----
`Space` is an abstract base class. In order to allow for instantiation
the methods `get_distance_array`, `total_volume` and `copy` must be
implemented as well as the abstract methods inherited from
`DomainObject`.
the methods `total_volume` and `copy` must be implemented as well as the
abstract methods inherited from `DomainObject`.
"""
def __init__(self):
super(Space, self).__init__()
@abc.abstractproperty
def harmonic(self):
""" Returns True if this space is a harmonic space.
"""
raise NotImplementedError
@abc.abstractproperty
......@@ -76,9 +73,7 @@ class Space(DomainObject):
-------
float
A real number representing the sum of all pixel volumes.
"""
raise NotImplementedError(
"There is no generic volume for the Space base class.")
......@@ -90,16 +85,12 @@ class Space(DomainObject):
-------
Space
A copy of this instance.
"""
return self.__class__()
def get_distance_array(self, distribution_strategy):
""" The distances of the pixel to zero.
This returns an array that gives for each pixel its distance to the
center of the manifolds grid.
""" The length of the k vector for every pixel.
This method is only implemented for harmonic spaces.
Parameters
----------
......@@ -110,18 +101,21 @@ class Space(DomainObject):
Returns
-------
distributed_data_object
A d2o containing the distances
A d2o containing the k vector lengths
"""
raise NotImplementedError(
"There is no generic distance structure for Space base class.")
raise NotImplementedError
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
def get_natural_binbounds(self):
""" The boundaries for natural power spectrum binning.
This method is only implemented for harmonic spaces.
Returns
-------
......@@ -164,6 +158,7 @@ class Space(DomainObject):
def hermitianize_inverter(self, x, axes):
""" Inverts/flips x in the context of Hermitian decomposition.
This method is only implemented for harmonic spaces.
This method is mainly used for power-synthesizing and -analyzing
Fields.
......
......@@ -28,7 +28,8 @@ class LaplaceOperatorTests(unittest.TestCase):
@expand(product([None, False, True], [False, True], [10, 100, 1000]))
def test_Laplace(self, log1, log2, sz):
s = ift.RGSpace(sz, harmonic=True)
p = ift.PowerSpace(s, logarithmic=log1)
p = ift.PowerSpace(s, binbounds=ift.PowerSpace.useful_binbounds(s,
logarithmic=log1))
L = ift.LaplaceOperator(p, logarithmic=log2)
arr = np.random.random(p.shape[0])
fp = ift.Field(p, val=arr)
......
......@@ -18,8 +18,7 @@
import unittest
import numpy as np
from numpy.testing import assert_equal, assert_approx_equal,\
assert_allclose
from numpy.testing import assert_approx_equal, assert_allclose
from nifty import Field,\
RGSpace,\
......@@ -36,18 +35,19 @@ def _get_rtol(tp):
else:
return 1e-5
class SmoothingOperator_Tests(unittest.TestCase):
spaces = [RGSpace(128)]
@expand(product(spaces, [0., .5, 5.], [True, False]))
def test_property(self, space, sigma, log_distances):
op = SmoothingOperator.make(space, sigma=sigma,
log_distances=log_distances)
log_distances=log_distances)
if op.domain[0] != space:
raise TypeError
if op.unitary != False:
if op.unitary:
raise ValueError
if op.self_adjoint != True:
if not op.self_adjoint:
raise ValueError
if op.sigma != sigma:
raise ValueError
......@@ -57,7 +57,7 @@ class SmoothingOperator_Tests(unittest.TestCase):
@expand(product(spaces, [0., .5, 5.], [True, False]))
def test_adjoint_times(self, space, sigma, log_distances):
op = SmoothingOperator.make(space, sigma=sigma,
log_distances=log_distances)
log_distances=log_distances)
rand1 = Field.from_random('normal', domain=space)
rand2 = Field.from_random('normal', domain=space)
tt1 = rand1.vdot(op.times(rand2))
......@@ -67,7 +67,7 @@ class SmoothingOperator_Tests(unittest.TestCase):
@expand(product(spaces, [0., .5, 5.], [False]))
def test_times(self, space, sigma, log_distances):
op = SmoothingOperator.make(space, sigma=sigma,
log_distances=log_distances)
log_distances=log_distances)
rand1 = Field(space, val=0.)
rand1.val[0] = 1.
tt1 = op.times(rand1)
......@@ -76,7 +76,7 @@ class SmoothingOperator_Tests(unittest.TestCase):
@expand(product(spaces, [0., .5, 5.], [True, False]))
def test_inverse_adjoint_times(self, space, sigma, log_distances):
op = SmoothingOperator.make(space, sigma=sigma,
log_distances=log_distances)
log_distances=log_distances)
rand1 = Field.from_random('normal', domain=space)
rand2 = Field.from_random('normal', domain=space)
tt1 = rand1.vdot(op.inverse_times(rand2))
......@@ -87,7 +87,7 @@ class SmoothingOperator_Tests(unittest.TestCase):
[np.float64, np.complex128]))
def test_smooth_regular1(self, sz, d, sigma, 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)
inp = Field.from_random(domain=sp, random_type='normal', std=1, mean=4,
dtype=tp)
......@@ -99,7 +99,7 @@ class SmoothingOperator_Tests(unittest.TestCase):
[np.float64, np.complex128]))
def test_smooth_regular2(self, sz1, sz2, d1, d2, sigma, 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)
inp = Field.from_random(domain=sp, random_type='normal', std=1, mean=4,
dtype=tp)
......@@ -112,7 +112,8 @@ class SmoothingOperator_Tests(unittest.TestCase):
def test_smooth_irregular1(self, sz, log, sigma, tp):
tol = _get_rtol(tp)
sp = RGSpace(sz, harmonic=True)
ps = PowerSpace(sp, nbin=sz, logarithmic=log)
ps = PowerSpace(sp, binbounds=PowerSpace.useful_binbounds(
sp, logarithmic=log))
smo = SmoothingOperator.make(ps, sigma=sigma)
inp = Field.from_random(domain=ps, random_type='normal', std=1, mean=4,
dtype=tp)
......@@ -125,7 +126,8 @@ class SmoothingOperator_Tests(unittest.TestCase):
def test_smooth_irregular2(self, sz1, sz2, log, sigma, tp):
tol = _get_rtol(tp)
sp = RGSpace([sz1, sz2], harmonic=True)
ps = PowerSpace(sp, logarithmic=log)
ps = PowerSpace(sp, binbounds=PowerSpace.useful_binbounds(
sp, logarithmic=log))
smo = SmoothingOperator.make(ps, sigma=sigma)
inp = Field.from_random(domain=ps, random_type='normal', std=1, mean=4,
dtype=tp)
......
......@@ -41,8 +41,8 @@ HARMONIC_SPACES = [RGSpace((8,), harmonic=True),
LMSpace(9)]
#Try all sensible kinds of combinations of spaces, distributuion strategy and
#binning parameters
# Try all sensible kinds of combinations of spaces, distribution strategy and
# binning parameters
CONSISTENCY_CONFIGS_IMPLICIT = product(HARMONIC_SPACES,
["not", "equal", "fftw"],
[None], [None, 3, 4], [True, False])
......@@ -70,14 +70,14 @@ CONSTRUCTOR_CONFIGS = [
}],
[RGSpace((8,), harmonic=True), 'not', True, None, None, {
'harmonic': True,
'shape': (2,),
'dim': 2,
'shape': (4,),
'dim': 4,
'total_volume': 8.0,
'harmonic_partner': RGSpace((8,), harmonic=True),
'binbounds': (0.70710678118654757,),
'binbounds': (0.5, 1.3228756555322954, 3.5),
'pindex': distributed_data_object([0, 1, 1, 1, 1, 1, 1, 1]),
'kindex': np.array([0., 2.28571429]),
'rho': np.array([1, 7]),
'kindex': np.array([0., 1., 2.5, 4.]),
'rho': np.array([1, 2, 4, 1]),
}],
]
......@@ -123,18 +123,20 @@ class PowerSpaceInterfaceTest(unittest.TestCase):
class PowerSpaceConsistencyCheck(unittest.TestCase):
@expand(CONSISTENCY_CONFIGS)
def test_rhopindexConsistency(self, harmonic_partner, distribution_strategy,
def test_rhopindexConsistency(self, harmonic_partner,
distribution_strategy,
binbounds, nbin, logarithmic):
if distribution_strategy == "fftw":
if not hasattr(gdi.get('fftw'), 'FFTW_MPI'):
raise SkipTest
p = PowerSpace(harmonic_partner=harmonic_partner,
distribution_strategy=distribution_strategy,
logarithmic=logarithmic, nbin=nbin,
binbounds=binbounds)
distribution_strategy=distribution_strategy,
binbounds=PowerSpace.useful_binbounds(
harmonic_partner, logarithmic, nbin))
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')
class PowerSpaceFunctionalityTest(unittest.TestCase):
@expand(CONSTRUCTOR_CONFIGS)
......@@ -147,13 +149,13 @@ class PowerSpaceFunctionalityTest(unittest.TestCase):
with assert_raises(expected['error']):
PowerSpace(harmonic_partner=harmonic_partner,
distribution_strategy=distribution_strategy,
logarithmic=logarithmic, nbin=nbin,
binbounds=binbounds)
binbounds=PowerSpace.useful_binbounds(
harmonic_partner, logarithmic, nbin))
else:
p = PowerSpace(harmonic_partner=harmonic_partner,