Commit 5d351e29 authored by Martin Reinecke's avatar Martin Reinecke

new convenience method for binbounds

parent da18ee57
......@@ -113,6 +113,35 @@ class PowerSpace(Space):
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, binbounds=None):
super(PowerSpace, self).__init__()
self._ignore_for_hash += ['_pindex', '_kindex', '_rho']
......@@ -179,7 +208,8 @@ class PowerSpace(Space):
binbounds=self._binbounds)
def weight(self):
return 1.
# MR FIXME: this will probably change to 1 soon
return np.asarray(self.rho, dtype=np.float64)
def get_distance_array(self):
return self.kindex.copy()
......
......@@ -44,35 +44,3 @@ 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,15 @@ import numpy as np
import nifty2go as ift
from numpy.testing import assert_allclose
from itertools import product
from test.common import expand, marco_binbounds
from test.common import expand
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, binbounds=marco_binbounds(s, logarithmic=log1))
bb = ift.PowerSpace.useful_binbounds(s, logarithmic=log1)
p = ift.PowerSpace(s, binbounds=bb)
L = ift.LaplaceOperator(p, logarithmic=log2)
arr = np.random.random(p.shape[0])
fp = ift.Field(p, val=arr)
......
......@@ -27,7 +27,7 @@ from nifty2go import Field,\
DirectSmoothingOperator
from itertools import product
from test.common import expand, marco_binbounds
from test.common import expand
def _get_rtol(tp):
......@@ -94,8 +94,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, binbounds=marco_binbounds(
sp, logarithmic=log, nbin=sz))
bb = PowerSpace.useful_binbounds(sp, logarithmic=log)
ps = PowerSpace(sp, binbounds=bb)
smo = DirectSmoothingOperator(ps, sigma=sigma)
inp = Field.from_random(domain=ps, random_type='normal', std=1, mean=4,
dtype=tp)
......@@ -107,7 +107,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, binbounds=marco_binbounds(sp, logarithmic=log))
bb = PowerSpace.useful_binbounds(sp, logarithmic=log)
ps = PowerSpace(sp, binbounds=bb)
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 nifty2go import PowerSpace, RGSpace, Space, LMSpace
from test.common import expand, marco_binbounds
from test.common import expand
from itertools import product, chain
HARMONIC_SPACES = [RGSpace((8,), harmonic=True),
......@@ -46,7 +46,7 @@ CONSISTENCY_CONFIGS = chain(CONSISTENCY_CONFIGS_IMPLICIT,
# [harmonic_partner, logarithmic, nbin, binbounds, expected]
CONSTRUCTOR_CONFIGS = [
[1, False, None, None, {'error': ValueError}],
[1, False, None, None, {'error': (ValueError, NotImplementedError)}],
[RGSpace((8,)), False, None, None, {'error': ValueError}],
[RGSpace((8,), harmonic=True), None, None, None, {
'harmonic': True,
......@@ -61,14 +61,14 @@ CONSTRUCTOR_CONFIGS = [
}],
[RGSpace((8,), harmonic=True), 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,),
'pindex': np.array([0, 1, 1, 1, 1, 1, 1, 1]),
'kindex': np.array([0., 2.28571429]),
'rho': np.array([1, 7]),
'binbounds': (0.5, 1.3228756555322954, 3.5),
'pindex': np.array([0, 1, 2, 2, 3, 2, 2, 1]),
'kindex': np.array([0., 1., 2.5, 4.]),
'rho': np.array([1, 2, 4, 1]),
}],
]
......@@ -98,9 +98,8 @@ class PowerSpaceConsistencyCheck(unittest.TestCase):
@expand(CONSISTENCY_CONFIGS)
def test_rhopindexConsistency(self, harmonic_partner,
binbounds, nbin, logarithmic):
p = PowerSpace(harmonic_partner=harmonic_partner,
binbounds=marco_binbounds(harmonic_partner,
logarithmic, nbin))
bb = PowerSpace.useful_binbounds(harmonic_partner, logarithmic, nbin)
p = PowerSpace(harmonic_partner=harmonic_partner, binbounds=bb)
assert_equal(np.bincount(p.pindex.flatten()), p.rho,
err_msg='rho is not equal to pindex degeneracy')
......@@ -112,13 +111,13 @@ class PowerSpaceFunctionalityTest(unittest.TestCase):
logarithmic, nbin, binbounds, expected):
if 'error' in expected:
with assert_raises(expected['error']):
PowerSpace(harmonic_partner=harmonic_partner,
binbounds=marco_binbounds(harmonic_partner,
logarithmic, nbin))
bb = PowerSpace.useful_binbounds(harmonic_partner,
logarithmic, nbin)
PowerSpace(harmonic_partner=harmonic_partner, binbounds=bb)
else:
p = PowerSpace(harmonic_partner=harmonic_partner,
binbounds=marco_binbounds(harmonic_partner,
logarithmic, nbin))
bb = PowerSpace.useful_binbounds(harmonic_partner,
logarithmic, nbin)
p = PowerSpace(harmonic_partner=harmonic_partner, binbounds=bb)
for key, value in expected.items():
if isinstance(value, np.ndarray):
assert_almost_equal(getattr(p, key), value)
......@@ -132,4 +131,4 @@ class PowerSpaceFunctionalityTest(unittest.TestCase):
def test_weight(self):
p = PowerSpace(harmonic_partner=RGSpace(10,harmonic=True))
assert_almost_equal(p.weight(),1.)
assert_almost_equal(p.weight(),p.rho)
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