Commit 751ba2c3 authored by Martin Reinecke's avatar Martin Reinecke
Browse files

merge fix_powerspace

parents ec41028a 3102f53d
......@@ -216,8 +216,8 @@ class Field(Loggable, 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
......@@ -231,16 +231,8 @@ class Field(Loggable, 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
......@@ -308,8 +300,6 @@ class Field(Loggable, object):
parts = [self._single_power_analyze(
work_field=part,
space_index=space_index,
logarithmic=logarithmic,
nbin=nbin,
binbounds=binbounds)
for part in parts]
......@@ -321,8 +311,7 @@ class Field(Loggable, 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(
......@@ -336,7 +325,6 @@ class Field(Loggable, object):
harmonic_domain = work_field.domain[space_index]
power_domain = PowerSpace(harmonic_partner=harmonic_domain,
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.")
......
......@@ -142,12 +142,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
......
......@@ -16,7 +16,6 @@
# NIFTy is being developed at the Max-Planck-Institut fuer Astrophysik
# and financially supported by the Studienstiftung des deutschen Volkes.
#from builtins import str
import ast
import numpy as np
......@@ -31,25 +30,19 @@ class PowerSpace(Space):
----------
harmonic_partner : Space
The harmonic Space of which this is the power space.
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
----------
......@@ -64,13 +57,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
-----
......@@ -83,8 +77,44 @@ class PowerSpace(Space):
# ---Overwritten properties and methods---
def __init__(self, harmonic_partner,
logarithmic=None, nbin=None, binbounds=None):
@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)
def __init__(self, harmonic_partner, binbounds=None):
super(PowerSpace, self).__init__()
self._ignore_for_hash += ['_pindex', '_kindex', '_rho']
......@@ -93,34 +123,22 @@ 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, logarithmic, nbin, binbounds)
key = (harmonic_partner, binbounds)
if self._powerIndexCache.get(key) is None:
distance_array = \
self.harmonic_partner.get_distance_array()
temp_binbounds = self._compute_binbounds(
harmonic_partner=self.harmonic_partner,
logarithmic=logarithmic,
nbin=nbin,
binbounds=binbounds)
temp_pindex = self._compute_pindex(
harmonic_partner=self.harmonic_partner,
distance_array=distance_array,
binbounds=temp_binbounds)
#temp_rho = temp_pindex.bincount().get_full_data()
binbounds=binbounds)
temp_rho = np.bincount(temp_pindex.flatten())
#temp_kindex = \
# (temp_pindex.bincount(weights=distance_array).get_full_data() /
# temp_rho)
temp_kindex = np.bincount(temp_pindex.flatten(),weights=distance_array.flatten())/temp_rho
self._powerIndexCache[key] = (temp_binbounds,
assert not np.any(temp_rho == 0), "empty bins detected"
temp_kindex = np.bincount(temp_pindex.flatten(),
weights=distance_array.flatten()) / temp_rho
self._powerIndexCache[key] = (binbounds,
temp_pindex,
temp_kindex,
temp_rho)
......@@ -128,45 +146,9 @@ class PowerSpace(Space):
(self._binbounds, self._pindex, self._kindex, self._rho) = \
self._powerIndexCache[key]
@staticmethod
def _compute_binbounds(harmonic_partner,
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):
pindex = np.empty(distance_array.shape,dtype=np.int)
if binbounds is None:
binbounds = harmonic_partner.get_natural_binbounds()
return np.searchsorted(binbounds, distance_array)
......
......@@ -161,6 +161,8 @@ class RGSpace(Space):
"""
if (not self.harmonic):
raise NotImplementedError
shape = self.shape
slice_of_first_dimension = slice(0, shape[0])
......@@ -192,6 +194,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
......@@ -218,10 +222,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,32 +85,31 @@ class Space(DomainObject):
-------
Space
A copy of this instance.
"""
return self.__class__()
def get_distance_array(self):
""" 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.
Returns
-------
numpy.ndarray
An array containing the distances
An array 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
-------
......@@ -158,6 +152,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.
......
......@@ -18,8 +18,8 @@
from builtins import str
from parameterized import parameterized
from nifty import RGSpace, LMSpace, HPSpace, GLSpace, PowerSpace
from nifty import Space, RGSpace, LMSpace, HPSpace, GLSpace, PowerSpace
import numpy as np
def custom_name_func(testcase_func, param_num, param):
return "%s_%s" % (
......@@ -42,3 +42,35 @@ def generate_spaces():
def generate_harmonic_spaces():
spaces = [RGSpace(4, harmonic=True), LMSpace(5)]
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
import nifty as ift
from numpy.testing import assert_allclose
from itertools import product
from test.common import expand
from test.common import expand, marco_binbounds
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=marco_binbounds(s, logarithmic=log1))
L = ift.LaplaceOperator(p, logarithmic=log2)
arr = np.random.random(p.shape[0])
fp = ift.Field(p, val=arr)
......
......@@ -27,7 +27,7 @@ from nifty import Field,\
DirectSmoothingOperator
from itertools import product
from test.common import expand
from test.common import expand, marco_binbounds
def _get_rtol(tp):
......@@ -72,7 +72,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, distances=d)
smo = FFTSmoothingOperator(sp, sigma=sigma)
inp = Field.from_random(domain=sp, random_type='normal', std=1, mean=4,
dtype=tp)
......@@ -83,7 +83,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])
smo = FFTSmoothingOperator(sp, sigma=sigma)
inp = Field.from_random(domain=sp, random_type='normal', std=1, mean=4,
dtype=tp)
......@@ -95,7 +95,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=marco_binbounds(
sp, logarithmic=log, nbin=sz))
smo = DirectSmoothingOperator(ps, sigma=sigma)
inp = Field.from_random(domain=ps, random_type='normal', std=1, mean=4,
dtype=tp)
......@@ -107,7 +108,7 @@ 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=marco_binbounds(sp, logarithmic=log))
smo = DirectSmoothingOperator(ps, sigma=sigma)
inp = Field.from_random(domain=ps, random_type='normal', std=1, mean=4,
dtype=tp)
......
......@@ -24,7 +24,7 @@ import numpy as np
from numpy.testing import assert_, assert_equal, assert_almost_equal,\
assert_raises
from nifty import PowerSpace, RGSpace, Space, LMSpace
from test.common import expand
from test.common import expand, marco_binbounds
from itertools import product, chain
HARMONIC_SPACES = [RGSpace((8,), harmonic=True),
......@@ -117,11 +117,12 @@ class PowerSpaceConsistencyCheck(unittest.TestCase):
def test_rhopindexConsistency(self, harmonic_partner,
binbounds, nbin, logarithmic):
p = PowerSpace(harmonic_partner=harmonic_partner,
logarithmic=logarithmic, nbin=nbin,
binbounds=binbounds)
binbounds=marco_binbounds(harmonic_partner,
logarithmic, nbin))
assert_equal(np.bincount(p.pindex.flatten()), 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)
......@@ -130,12 +131,12 @@ class PowerSpaceFunctionalityTest(unittest.TestCase):
if 'error' in expected:
with assert_raises(expected['error']):
PowerSpace(harmonic_partner=harmonic_partner,
logarithmic=logarithmic, nbin=nbin,
binbounds=binbounds)
binbounds=marco_binbounds(harmonic_partner,
logarithmic, nbin))
else:
p = PowerSpace(harmonic_partner=harmonic_partner,
logarithmic=logarithmic, nbin=nbin,
binbounds=binbounds)
binbounds=marco_binbounds(harmonic_partner,
logarithmic, nbin))
for key, value in expected.items():
if isinstance(value, np.ndarray):
assert_almost_equal(getattr(p, key), value)
......
......@@ -75,7 +75,7 @@ CONSTRUCTOR_CONFIGS = [
def get_distance_array_configs():
# for RGSpace(shape=(4, 4), distances=None)
# for RGSpace(shape=(4, 4), distances=(0.25,0.25))
cords_0 = np.ogrid[0:4, 0:4]
da_0 = ((cords_0[0] - 4 // 2) * 0.25)**2
da_0 = np.fft.ifftshift(da_0)
......@@ -84,7 +84,7 @@ def get_distance_array_configs():
da_0 = da_0 + temp
da_0 = np.sqrt(da_0)
return [
[(4, 4), None, da_0],
[(4, 4), (0.25, 0.25), da_0],
]
......@@ -143,7 +143,7 @@ class RGSpaceFunctionalityTests(unittest.TestCase):
@expand(get_distance_array_configs())
def test_distance_array(self, shape, distances, expected):
r = RGSpace(shape=shape, distances=distances)
r = RGSpace(shape=shape, distances=distances, harmonic=True)
assert_almost_equal(r.get_distance_array(), expected)
@expand(get_weight_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